From 4337b10a87bc2496fec2031d95ef234c77421194 Mon Sep 17 00:00:00 2001 From: flandre Date: Mon, 13 Jan 2025 14:20:54 +0800 Subject: [PATCH] add das --- src/beamformer/das.py | 139 +++++++++++++++ src/beamformer/dist.py | 69 ++++++++ src/beamformer/kernels.py | 101 +++++++++++ src/beamformer/process.py | 44 +++++ src/utils/Config.py | 54 ++++++ src/{ => utils}/Msg.py | 0 src/{ => utils}/RfFile.py | 0 src/utils/ScanData.py | 348 ++++++++++++++++++++++++++++++++++++++ test/main.ipynb | 37 ++++ test/test.py | 0 10 files changed, 792 insertions(+) create mode 100644 src/beamformer/das.py create mode 100644 src/beamformer/dist.py create mode 100644 src/beamformer/kernels.py create mode 100644 src/beamformer/process.py create mode 100644 src/utils/Config.py rename src/{ => utils}/Msg.py (100%) rename src/{ => utils}/RfFile.py (100%) create mode 100644 src/utils/ScanData.py create mode 100644 test/main.ipynb create mode 100644 test/test.py diff --git a/src/beamformer/das.py b/src/beamformer/das.py new file mode 100644 index 0000000..63e4ccc --- /dev/null +++ b/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() diff --git a/src/beamformer/dist.py b/src/beamformer/dist.py new file mode 100644 index 0000000..56435cb --- /dev/null +++ b/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()) diff --git a/src/beamformer/kernels.py b/src/beamformer/kernels.py new file mode 100644 index 0000000..4203ae1 --- /dev/null +++ b/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() diff --git a/src/beamformer/process.py b/src/beamformer/process.py new file mode 100644 index 0000000..5885127 --- /dev/null +++ b/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() + ) diff --git a/src/utils/Config.py b/src/utils/Config.py new file mode 100644 index 0000000..9f18ed9 --- /dev/null +++ b/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 diff --git a/src/Msg.py b/src/utils/Msg.py similarity index 100% rename from src/Msg.py rename to src/utils/Msg.py diff --git a/src/RfFile.py b/src/utils/RfFile.py similarity index 100% rename from src/RfFile.py rename to src/utils/RfFile.py diff --git a/src/utils/ScanData.py b/src/utils/ScanData.py new file mode 100644 index 0000000..10b5662 --- /dev/null +++ b/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() diff --git a/test/main.ipynb b/test/main.ipynb new file mode 100644 index 0000000..54f657b --- /dev/null +++ b/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 +} diff --git a/test/test.py b/test/test.py new file mode 100644 index 0000000..e69de29