add das

src/beamformer/das.py

@@ -0,0 +1,139 @@
+from threading import Thread
+from typing import Callable
+
+import cupy as cp
+import numpy as np
+
+from ds.Config import DeviceConfig
+from imaging.dist import refraction_dist, direct_dist
+from imaging.kernels import dist_mat_to_yids, dist_mat_to_yids_pwi
+
+
+def repeat_range_in_axis(shape: tuple, axis: int,p=cp):
+    idx = [None for _ in shape]
+    idx[axis] = slice(None)
+    idx = tuple(idx)
+    return p.zeros(shape, dtype=p.uint8) + p.arange(shape[axis], dtype=p.uint8)[idx]
+
+
+def gen_pwi(dist_mat: cp.ndarray, dev_cfg=DeviceConfig()):
+    col_idx = repeat_range_in_axis((dev_cfg.rows, dev_cfg.cols, dev_cfg.cols), 2)
+    row_idx = dist_mat_to_yids_pwi(dist_mat, dev_cfg)()
+    last_available_index = cp.where(row_idx[:, 128, 128] == (dev_cfg.rows))[0][0]
+
+    def pwi_(mat_in):
+        return mat_in[col_idx, row_idx].sum(axis=2).T[:, :last_available_index]
+
+    return pwi_, row_idx
+
+
+class TFM:
+    def __init__(self, device_cfg=DeviceConfig()):
+        self.device_cfg = device_cfg
+        self.yks: list[cp.ndarray] = []
+        self.xks: list[cp.ndarray] = []
+        self.eks: list[cp.ndarray] = []
+        self.canvas: list[cp.ndarray] = []
+        self.inputs: list[cp.ndarray] = []
+        self.s2 = 8
+        for gpu_id in range(4):
+            with cp.cuda.Device(gpu_id):
+                self.canvas.append(cp.zeros((self.device_cfg.rows, 64), dtype=cp.int32))
+                self.eks.append(repeat_range_in_axis((device_cfg.rows, self.s2, device_cfg.cols, device_cfg.cols), 2))
+                self.xks.append(repeat_range_in_axis((device_cfg.rows, self.s2, device_cfg.cols, device_cfg.cols), 3))
+
+    def load_yids(self, gkx: Callable[[int], cp.ndarray]):
+        self.yks.clear()
+        for gpu_id in range(4):
+            with cp.cuda.Device(gpu_id):
+                self.yks.append(gkx(gpu_id * 64))
+
+    def load_bscan(self, bscan_mat_cpu):
+        self.inputs.clear()
+        for gpu_id in range(4):
+            with cp.cuda.Device(gpu_id):
+                self.inputs.append(cp.asarray(bscan_mat_cpu))
+
+    def get_idx_mat(self):
+        t = np.zeros((
+            self.device_cfg.rows,
+            self.device_cfg.cols,
+            self.device_cfg.cols,
+            self.device_cfg.cols
+        ), dtype=np.uint16)
+
+        def send_thread(target, arr, idx):
+            target[:, idx * 64:(idx + 1) * 64, :, :] = arr[i].get()
+
+        ts = []
+        for i in range(4):
+            th = Thread(target=send_thread, args=(t, self.yks, i))
+            th.start()
+            ts.append(th)
+        for th in ts:
+            th.join()
+        return t
+
+    def __call__(self, bscan_mat_cpu=None):
+        if bscan_mat_cpu is not None:
+            self.load_bscan(bscan_mat_cpu)
+        assert self.inputs.__len__() > 0
+        ts = []
+        canvas_cpu = np.zeros((self.device_cfg.rows, self.device_cfg.cols))
+        for gpu_id in range(4):
+            with cp.cuda.Device(gpu_id):
+                for i in range(int(64 / self.s2)):
+                    i_start = i * self.s2
+                    i_end = (i + 1) * self.s2
+                    a1 = self.eks[gpu_id]
+                    a2 = self.xks[gpu_id]
+                    a3 = self.yks[gpu_id][:, i_start:i_end, :, :]
+                    if self.inputs[gpu_id].shape.__len__() == 2:
+                        res = self.inputs[gpu_id][a2, a3]
+                    else:
+                        res = self.inputs[gpu_id][a1, a2, a3]
+                    res = res.sum(axis=2)
+                    res = res.sum(axis=2)
+                    self.canvas[gpu_id][:, i_start:i_end] = res
+
+            def send_thread(canvas_, gpu_id_):
+                canvas_cpu[:, gpu_id_ * 64:(gpu_id_ + 1) * 64] = canvas_[gpu_id_].get()
+
+            t = Thread(target=send_thread, args=(self.canvas, gpu_id))
+            ts.append(t)
+            t.start()
+        for t in ts:
+            t.join()
+        return canvas_cpu.T
+
+
+def test1():
+    device_cfg = DeviceConfig()
+    das = TFM()
+    das.load_yids()
+    das(np.ones((device_cfg.cols, device_cfg.cols, device_cfg.rows), dtype=np.int16))
+
+
+def test2():
+    device_cfg = DeviceConfig()
+    tfm = TFM()
+    tfm.load_yids(dist_mat_to_yids(direct_dist()))
+    r = tfm(np.ones((device_cfg.cols, device_cfg.cols, device_cfg.rows), dtype=np.int16))
+    print(r.shape)
+
+
+def test3():
+    device_cfg = DeviceConfig()
+    pwi, pwi_row_idx = gen_pwi(refraction_dist())
+    pwi(cp.ones((device_cfg.cols, device_cfg.rows), dtype=np.int16))
+
+def test4():
+    device_cfg = DeviceConfig()
+    pwi, pwi_row_idx = gen_pwi(refraction_dist())
+
+
+if __name__ == '__main__':
+    # test1()
+    # test2()
+    # test3()
+    test4()

src/beamformer/dist.py

@@ -0,0 +1,69 @@
+"""
+dist function return mat f(y,x),
+where f(y,x) is the pixel distance in echo pulse RF signal between point(0,0) and point(x,y).
+"""
+import numpy as np
+from scipy.optimize import fsolve
+
+from ds.Config import DeviceConfig
+from utils.filename import DS
+
+
+def refraction_dist(dev_cfg=DeviceConfig()):
+    v1 = dev_cfg.v1
+    v1x = dev_cfg.v1x
+    v2 = dev_cfg.v2
+    y1 = dev_cfg.y1
+    SECONDS_PER_Y_PIX = dev_cfg.second_per_y_pix
+    METER_PER_X_PIX = dev_cfg.meter_per_x_pix
+
+    save_path = DS / f'{y1}_{v1x}_{v2}_fslove.npy'
+    if save_path.exists():
+        return np.load(str(save_path))
+
+    def gf(x, y2):
+        def f(arg):
+            θ1 = arg[0]
+            return [
+                y1 * np.tan(θ1) + v2 * SECONDS_PER_Y_PIX * y2 * np.tan(
+                    np.arcsin(v2 / v1x * np.sin(θ1))) - METER_PER_X_PIX * x]
+
+        return f
+
+    iter_start = np.arcsin(v1x / v2 * np.sin(np.deg2rad(90)))
+    iter_start = iter_start // 1e-7 / 1e+7
+    m = np.zeros((dev_cfg.rows, dev_cfg.cols - 1))
+    for y2_ in range(m.shape[0]):
+        for x_ in range(m.shape[1]):
+            m[y2_, x_] = fsolve(gf(x_ + 1, y2_ + 1), np.array([iter_start]))[0]
+
+    dmat = np.zeros((dev_cfg.rows, dev_cfg.cols))
+    t1 = (np.arange(dev_cfg.rows)) * SECONDS_PER_Y_PIX
+    y2_mat = (t1 * v2)[:, np.newaxis]
+    x1_mat = y1 * np.tan(m)
+    x2_mat = y2_mat * np.tan(np.arcsin((v2 / v1x) * np.sin(m)))
+    d1 = np.hypot(x1_mat, y1)
+    d2 = np.hypot(x2_mat, y2_mat)
+    dmat[:, 1:] = (d1 / v1 + d2 / v2)
+    dmat[:, 0] = y1 / v1 + t1
+    dmat /= SECONDS_PER_Y_PIX
+    np.save(save_path.__str__(), dmat)
+    return dmat
+
+
+def direct_dist(dev_cfg=DeviceConfig(), p=np):
+    y, x = p.ogrid[:dev_cfg.rows, :dev_cfg.cols]
+    y1_meter = dev_cfg.y1
+    y2_meter = y * dev_cfg.second_per_y_pix * dev_cfg.v2
+    y_meter = y1_meter + y2_meter
+    x_meter = x * dev_cfg.meter_per_x_pix
+    xy_meter = p.hypot(y_meter, x_meter)
+    xy1_meter = xy_meter * (y1_meter / y_meter)
+    xy2_meter = xy_meter * (y2_meter / y_meter)
+    m = xy1_meter / dev_cfg.v1 + xy2_meter / dev_cfg.v2
+    return m / dev_cfg.second_per_y_pix
+
+
+if __name__ == '__main__':
+    print(refraction_dist())
+    print(direct_dist())

src/beamformer/kernels.py

@@ -0,0 +1,101 @@
+"""
+yfi: y fire index
+yfs: y fire shape
+xfi: x fire index
+xfs: x fire shape
+xri: x receive index
+xrs: x receive shape
+"""
+import cupy as cp
+
+from ds.Config import DeviceConfig
+from imaging.dist import direct_dist
+
+
+def dist_mat_to_yids(dist_mat: cp.ndarray, dev_cfg=DeviceConfig()):
+    def yids(offset: int) -> cp.ndarray:
+        yfs = dev_cfg.rows
+        xfs = 64
+        xss = dev_cfg.cols
+        xrs = dev_cfg.cols
+        kernel = cp.RawKernel(f"""
+                extern "C" __global__ void k(unsigned short int *result, float * dist_mat, int offset)
+                {{
+                    const long long cols = {dev_cfg.cols};
+                    
+                    const long long yfi = blockIdx.x;
+                    const long long yfs = gridDim.x;
+                    
+                    const long long xfi = blockIdx.y;
+                    const long long xfs = gridDim.y;
+                    
+                    const long long xsi = blockIdx.z;
+                    const long long xss = gridDim.z;
+                    
+                    const long long xri = threadIdx.x;
+                    const long long xrs = blockDim.x;
+
+                    long long idx = xri + xsi*xrs + xfi*xrs*xss + yfi*xrs*xss*xfs;
+                    float r = dist_mat[yfi*cols+abs(xsi-(xfi + offset))] + dist_mat[yfi*cols+abs(xri-(xfi + offset))];
+                    result[idx] = min(__float2int_rd(r), (int)yfs);
+                }}
+        """, 'k')
+        result = cp.zeros(yfs * xfs * xss * xrs, dtype=cp.uint16)
+        kernel((yfs, xfs, xss), (xrs, 1, 1), (
+            result,
+            cp.asarray(dist_mat.astype(cp.float32)),
+            cp.int32(offset)
+        ))
+        result = result.reshape((yfs, xfs, xss, xrs))
+        return result
+
+    return yids
+
+
+def dist_mat_to_yids_pwi(dist_mat: cp.ndarray, dev_cfg=DeviceConfig()):
+    """
+    min in cuda is for __/\\__
+    :param dist_mat: cp.ndarray from dist.py
+    :param dev_cfg: DeviceConfig
+    :return: pixel distance (yf,xf,xr)
+    """
+
+    def yids() -> cp.ndarray:
+        yfs = dev_cfg.rows
+        xfs = dev_cfg.cols
+        xrs = dev_cfg.cols
+        kernel = cp.RawKernel(f"""
+                extern "C" __global__ void k(unsigned short int *result, float * dist_mat)
+                {{
+                    const long long cols = {dev_cfg.cols};
+                    
+                    const long long yfi = blockIdx.x;
+                    const long long yfs = gridDim.x;
+                    
+                    const long long xfi = blockIdx.y;
+                    const long long xfs = gridDim.y;
+                    
+                    const long long xri = threadIdx.x;
+                    const long long xrs = blockDim.x;
+
+                    long long idx = xri + xfi*xrs + yfi*xrs*xfs; //(yf,xf,xr)
+                    float r = dist_mat[yfi*cols] + dist_mat[yfi*cols+abs(xri-xfi)];
+                    result[idx] = min(__float2int_rd(r), (int)yfs); 
+                }}
+        """, 'k')
+        result = cp.zeros(yfs * xfs * xrs, dtype=cp.uint16)
+        kernel((yfs, xfs), (xrs, 1, 1), (result, cp.asarray(dist_mat.astype(cp.float32))))
+        result = result.reshape((yfs, xfs, xrs))
+        return result
+
+    return yids
+
+
+def test1():
+    f = dist_mat_to_yids(direct_dist())
+    r = f(0)
+    print(r.flags['C_CONTIGUOUS'], r.flags['F_CONTIGUOUS'])
+
+
+if __name__ == '__main__':
+    test1()

src/beamformer/process.py

@@ -0,0 +1,44 @@
+from ds.Config import ImagingConfig
+from imaging.ScanData import ScanData
+import cupy as cp
+
+from imaging.das import TFM
+
+tfm_cache = [None]
+
+
+def pwi_process(s: ScanData, icfg: ImagingConfig, pwi):
+    return (s
+            .filter_max_persent(icfg.bscan_max/1000, bid=True)
+            .sum(cond=s.d == 3)
+            .dct(icfg.dct_start, icfg.dct_end)
+            .call(lambda m: m.astype(cp.int16))
+            .call(pwi)
+            .call(cp.asarray, order='C')
+            .argrelextrema()
+            .conv_guass(b=icfg.beta * 0.01)
+            .clip(icfg.focus_start, icfg.focus_end)
+            .filter_max_persent(icfg.focus_max/100, mmax=icfg.focus_mmax if icfg.uafm else None)
+            .time_gain_compensation_linear_max(icfg.tgcl, mmax=icfg.focus_mmax if icfg.uafm else None)
+            .cpu()
+            .get()
+            )
+
+
+def tfm_process(s: ScanData, icfg: ImagingConfig, disable_cache: bool, tfm: TFM):
+    # print(icfg.changed_field, icfg.changed_field in ['dct_start', 'dct_end', 'bscan_max'], disable_cache)
+    if icfg.changed_field in ['dct_start', 'dct_end', 'bscan_max'] or tfm_cache[0] is None or disable_cache:
+        tfm_cache[0] = (s.dct(icfg.dct_start, icfg.dct_end)
+                        .filter_max_persent(icfg.bscan_max/1000, bid=True)
+                        .call(lambda m: m.astype(cp.int16))
+                        .call(tfm))
+    return (tfm_cache[0]
+            .call(cp.asarray, order='C')
+            .argrelextrema(axis=1)
+            .conv_guass(b=icfg.beta * 0.01, axis=1)
+            .clip(icfg.focus_start, icfg.focus_end)
+            .filter_max_persent(icfg.focus_max, mmax=icfg.focus_mmax if icfg.uafm else None)
+            .time_gain_compensation_linear_max(icfg.tgcl, mmax=icfg.focus_mmax if icfg.uafm else None)
+            .cpu()
+            .get()
+            )

src/utils/Config.py

@@ -0,0 +1,54 @@
+import dataclasses
+import json
+
+from utils.filename import LAST_CONFIG, CONFIG_FOLDER
+
+
+@dataclasses.dataclass
+class ImagingConfig:
+    focus_start: int = 0
+    focus_end: int = 1500
+    focus_max: int = 100
+    focus_min: int = 0
+    focus_mmax: int = 8192
+    uafm: bool = False
+
+    bscan_start: int = 0
+    bscan_end: int = 1500
+    bscan_max: int = 100
+    bscan_min: int = 0
+    bscan_mmax: int = 8192
+
+    dct_start: int = 300
+    dct_end: int = 550
+
+    eid: int = 0
+    beta: int = 1
+
+    tgcl: int = 0
+
+    changed_field: str = None
+    changed_field_orig_value: str | int = None
+
+    @staticmethod
+    def from_file(input_format, focus_mode, folder=CONFIG_FOLDER):
+        try:
+            return ImagingConfig(**json.load((folder / f'icfg_{input_format}_{focus_mode}.json').open()))
+        except Exception as e:
+            return ImagingConfig()
+
+    def save(self, input_format: str, focus_mode: str, folder=CONFIG_FOLDER):
+        json.dump(self.__dict__, (folder / f'icfg_{input_format}_{focus_mode}.json').open('w'))
+
+
+@dataclasses.dataclass
+class DeviceConfig:
+    second_per_y_pix: float = 2e-8
+    meter_per_x_pix: float = 2e-4
+
+    rows: int = 1500
+    cols: int = 256
+    v1: int = 1915
+    v1x: int = 1915
+    v2: int = 5900
+    y1: float = 1e-3

src/utils/ScanData.py

@@ -0,0 +1,348 @@
+from pathlib import Path
+
+import numpy as np
+
+import cupy as cp
+import cupyx
+import cv2
+import matplotlib
+import scipy
+import cupyx.scipy.signal
+import scipy.signal
+import cupyx.scipy.fft
+import cupyx.scipy.ndimage
+from cupyx.scipy.fft import dctn, idctn
+from scipy.stats import norm as norms
+
+
+def skip(f):
+    def wrapper(self, **kwargs):
+        if 'cond' not in kwargs:
+            return f(self, **kwargs)
+        if kwargs['cond']:
+            del kwargs['cond']
+            return f(self, **kwargs)
+        else:
+            return self
+
+    return wrapper
+
+
+class ScanData:
+    @staticmethod
+    def from_buffer(buf: bytes | memoryview, shape, p=np):
+        return ScanData(p.frombuffer(buf, dtype=p.int16, count=np.prod(shape)).reshape(shape))
+
+    @staticmethod
+    def from_file(filename: Path | str, shape, p=np):
+        filename = Path(filename)
+        return ScanData.from_buffer(filename.read_bytes(), shape, p)
+
+    @staticmethod
+    def from_npy(filename: Path, p=np):
+        return ScanData(p.load(filename))
+
+    def __init__(self, mat: cp.ndarray):
+        if isinstance(mat, np.ndarray):
+            self.p = np
+        elif isinstance(mat, cp.ndarray):
+            self.p = cp
+        else:
+            raise Exception('wrong type mat!!')
+        self.m = mat
+        self.d = mat.shape.__len__()
+
+    def save(self, path: Path):
+        self.p.save(path.__str__(), self.m)
+
+    # @property
+    # def p(self) -> np:
+    #     if self._p == 'np':
+    #         return np
+    #     if self._p == 'cp':
+    #         return cp
+
+    def get(self):
+        return self.m
+
+    def apply(self, f):
+        f(self.m)
+        return self
+
+    def call(self, f, *args, **kwargs):
+        return ScanData(f(self.m, *args, **kwargs))
+
+    def copy(self):
+        return ScanData(self.m.copy())
+
+    def cpu(self):
+        if self.p == cp:
+            return ScanData(self.m.get())
+        if self.p == np:
+            return self.copy()
+
+    def gpu(self):
+        if self.p == np:
+            return ScanData(cp.asarray(self.m))
+        if self.p == cp:
+            return self.copy()
+
+    def norm(self):
+        m = self.m.astype(self.p.float32)
+        m -= m.min()
+        mmax = m.max()
+        if mmax == 0:
+            return ScanData(self.p.zeros_like(m))
+        m /= mmax
+        return ScanData(m)
+
+    @skip
+    def sum(self, axis=0):
+        return ScanData(self.m.sum(axis=axis))
+
+    def rotate90(self):
+        return ScanData(cv2.rotate(self.m, cv2.ROTATE_90_CLOCKWISE))
+
+    def filter_min(self, t, v=None, bid=False):
+        m = self.m
+        if v is None:
+            v = t
+        if bid:
+            m[(0 <= m) & (m <= t)] = v
+            m[(-t <= m) & (m <= 0)] = -v
+        else:
+            m[(m <= t)] = v
+        return self
+
+    def filter_max(self, t, v=None, bid=False):
+        m = self.m.copy()
+        if v is None:
+            v = t
+        if bid:
+            m[m >= t] = v
+            m[-m >= t] = -v
+        else:
+            m[m >= t] = v
+        return ScanData(m)
+
+    def filter_max_persent(self, t, v=None, bid=False, mmax=None):
+        if mmax is None:
+            mmax = self.m.max()
+        else:
+            self.m[0, 0] = mmax
+        return self.filter_max(mmax * t, v, bid)
+
+    def astype(self, t):
+        return ScanData(self.m.astype(t))
+
+    def dctn_cp(self):
+        return ScanData(dctn(self.m))
+
+    def dctn_cp_(self):
+        self.m = dctn(self.m)
+        return self
+
+    def idctn_cp_(self):
+        self.m = idctn(self.m)
+        return self
+
+    def dct(self, mmin, mmax):
+        dct_ = scipy.fft.dct
+        idct = scipy.fft.idct
+        if self.p == cp:
+            dct_ = cupyx.scipy.fft.dct
+            idct = cupyx.scipy.fft.idct
+        m_dct = dct_(self.m)
+        if self.d == 2:
+            m_dct[:, mmax:] = 0
+            m_dct[:, :mmin] = 0
+        elif self.d == 3:
+            m_dct[:, :, mmax:] = 0
+            m_dct[:, :, :mmin] = 0
+        else:
+            raise NotImplementedError()
+        return ScanData(idct(m_dct))
+
+    def argrelextrema(self, axis=1):
+        arg = scipy.signal.argrelextrema
+        m = self.m
+        p = self.p
+        if p == cp:
+            arg = cupyx.scipy.signal.argrelextrema
+        rm = p.zeros_like(m)
+        indies1 = arg(m, p.greater, axis=axis)
+        bl = p.zeros_like(m, dtype=p.bool_)
+        indies2 = arg(m, p.less, axis=axis)
+        bl2 = p.zeros_like(m, dtype=p.bool_)
+        bl[indies1] = True
+        bl2[indies2] = True
+        i1 = bl & (m > 0)
+        i2 = bl2 & (m < 0)
+        idx = i1 | i2
+        rm[idx] = m[idx].__abs__()
+        return ScanData(rm)
+
+    def conv_guass(self, b=0.01, axis=1):
+        cv = scipy.ndimage.convolve1d
+        m = self.m
+        p = self.p
+        if p == cp:
+            cv = cupyx.scipy.ndimage.convolve1d
+        rv = norms(loc=0, scale=b)
+        x2 = np.arange(-1, 1.1, 0.1)
+        w = rv.pdf(x2)
+        if p == cp:
+            w = cp.asarray(w)
+        rm = cv(m, w, axis=axis)
+        return ScanData(rm)
+
+    def clip(self, start, end):
+        return ScanData(self.m[:, start:end])
+
+    def print_contiguous(self):
+        print(self.m.flags['C_CONTIGUOUS'], self.m.flags['F_CONTIGUOUS'])
+        return self
+
+    def rgb(self):
+        return ScanData((self.norm().m * 255).astype(np.uint8)[:, :, None].repeat(3, axis=2))
+
+    def pseudo_color(self):
+        m = self.m * 0.7
+        p = self.p
+        h = m
+        s = p.zeros_like(h) + 1
+        v = p.zeros_like(h) + 1
+        hsv = p.stack((h, s, v), axis=2)
+        if p == cp:
+            rgb = matplotlib.colors.hsv_to_rgb(hsv.get())
+        else:
+            rgb = matplotlib.colors.hsv_to_rgb(hsv)
+        return ScanData((rgb * 255).astype(np.uint8))
+
+    def pseudo_color_rgb(self):
+        m = (1 - self.m) * 0.7
+        p = self.p
+        h = m
+        s = p.zeros_like(h) + 1
+        v = p.zeros_like(h) + 1
+        hsv = p.stack((h, s, v), axis=2)
+        if p == cp:
+            rgb = matplotlib.colors.hsv_to_rgb(hsv.get())
+        else:
+            rgb = matplotlib.colors.hsv_to_rgb(hsv)
+        return ScanData((rgb * 255).astype(np.uint8))
+
+    def grey_color(self):
+        m = self.norm().m
+        return ScanData((m * 255).astype(np.uint8))
+
+    def time_gain_compensation_linear(self, scale: float, start: int = 0):
+        h = self.m.shape[-1]
+        addend = self.p.zeros((1, h), dtype=np.int64)
+        addend[:, start:] = self.p.arange(h - start) * scale
+        self.m += addend
+        return self
+
+    def time_gain_compensation_linear_float(self, scale: float, start: int = 0):
+        h = self.m.shape[-1]
+        addend = self.p.zeros((1, h), dtype=self.p.float32)
+        addend[:, start:] = (self.p.arange(h - start) * scale) + 1
+        print(addend)
+        return ScanData(self.m * addend)
+
+    def time_gain_compensation_linear_max(self, scale: float, mmax: int | None = None, start: int = 0):
+        if scale == 0:
+            return self
+        if mmax is None:
+            mmax = self.m.max()
+        h = self.m.shape[-1]
+        self.m = self.m.astype(np.float64)
+        mmax_arr = self.p.zeros(h) + mmax
+        mmax_arr[start:] -= self.p.arange(h - start) * scale
+        for i in range(h):
+            # a[1, a[1, :] > 99] = 99
+            self.m[self.m[:, i] > mmax_arr[i], i] = mmax_arr[i]
+            self.m[:, i] *= (mmax / mmax_arr[i])
+        self.m[self.m > mmax] = mmax
+        self.m = self.m.astype(np.int64)
+        return self
+
+    def hib(self):
+        if self.p == cp:
+            hilbert = cupyx.scipy.signal.hilbert
+        else:
+            hilbert = scipy.signal.hilbert
+        return ScanData(hilbert(self.m))
+
+    def abs(self):
+        return ScanData(self.p.abs(self.m))
+
+    def log(self):
+        return ScanData(self.p.log(self.m))
+
+    def show(self, cmap='grey', w=40, h=20, no_axis=False, hsv=False):
+        from matplotlib import pyplot as plt
+        # fig = plt.figure(figsize=(w, h))
+        fig, ax = plt.subplots(figsize=(w, h))
+        if no_axis:
+            ax.get_xaxis().set_ticks([])
+            ax.get_yaxis().set_ticks([])
+        image = self.m if self.p == np else self.m.get()
+        if hsv:
+            image = self.norm().pseudo_color_rgb().m
+        plt.imshow(image, cmap=cmap)
+        return self
+
+    def plot_center_line(self):
+        from matplotlib import pyplot as plt
+        plt.figure(figsize=(40, 5))
+        m = self.m if self.p == np else self.m.get()
+        plt.plot(m[128, :])
+        return self
+
+    def windows(self, fn, center_f, width_f, sps_f):
+        fft = scipy.fft.fft if self.p == np else cupyx.scipy.fft.fft
+        ifft = scipy.fft.ifft if self.p == np else cupyx.scipy.fft.ifft
+
+        seq_len = self.m.shape[-1]
+        per_f = sps_f / seq_len
+
+        center = int(center_f / per_f)
+        width = int(width_f / 2 / per_f)
+
+        start = center - width
+        end = center + width
+
+        fft_seq = fft(self.m)
+        window = fn(end - start)
+        filter = np.zeros(seq_len) + window[0]
+        filter[start:end] = window
+        filter[seq_len - end:seq_len - start] = window
+        fft_seq *= filter
+        return ScanData(self.p.real(ifft(fft_seq)))
+
+    def windows_plot(self):
+        pass
+
+    def dctn(self, xmin, xmax, ymin, ymax):
+        dctn = scipy.fft.dctn if self.p == np else cupyx.scipy.fft.dctn
+        idctn = scipy.fft.idctn if self.p == np else cupyx.scipy.fft.idctn
+
+        dctm = dctn(self.m)
+        dctm[xmin:xmax, ymin:ymax] = 0
+        return ScanData(idctn(dctm))
+
+    def resize_XY(self, xdyd):
+        print(self.m.shape[0])
+        if self.p != np:
+            raise NotImplementedError()
+        return ScanData(cv2.resize(self.m, (int(self.m.shape[1] * xdyd), self.m.shape[0])))
+
+
+def test1():
+    s = ScanData(np.array([1, 1, 2])).norm().call(np.add, 1)
+    print(s.m)
+
+
+if __name__ == '__main__':
+    test1()

test/main.ipynb

@@ -0,0 +1,37 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "initial_id",
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    ""
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}