From 0d8ee30b4ca2685c78a14374ba47c6dcd6d73e30 Mon Sep 17 00:00:00 2001 From: FMarmoreo Date: Tue, 7 May 2024 15:33:33 +0200 Subject: [PATCH] Vectorization of convolution functions in utils/convolve.py --- cobrawap/pipeline/utils/convolve.py | 202 +++++++++++++--------------- 1 file changed, 91 insertions(+), 111 deletions(-) diff --git a/cobrawap/pipeline/utils/convolve.py b/cobrawap/pipeline/utils/convolve.py index 2b5d850a..3a1c9fe3 100644 --- a/cobrawap/pipeline/utils/convolve.py +++ b/cobrawap/pipeline/utils/convolve.py @@ -1,74 +1,76 @@ -import numpy as np import itertools -import warnings +import numpy as np from scipy.ndimage import convolve as conv +from scipy.signal import convolve2d from types import SimpleNamespace +import warnings simple_3x3 = SimpleNamespace( - x=np.array([[-0, 0, 0], - [-1, 0, 1], - [-0, 0, 0]], dtype=float), - y=np.array([[-0, -1, -0], - [ 0, 0, 0], - [ 0, 1, 0]], dtype=float)) + x=np.array([[-0, 0, 0], + [-1, 0, 1], + [-0, 0, 0]], dtype=float), + y=np.array([[-0, -1, -0], + [0, 0, 0], + [0, 1, 0]], dtype=float)) prewitt_3x3 = SimpleNamespace( - x=np.array([[-1, 0, 1], - [-1, 0, 1], - [-1, 0, 1]], dtype=float), - y=np.array([[-1, -1, -1], - [ 0, 0, 0], - [ 1, 1, 1]], dtype=float)) + x=np.array([[-1, 0, 1], + [-1, 0, 1], + [-1, 0, 1]], dtype=float), + y=np.array([[-1, -1, -1], + [0, 0, 0], + [1, 1, 1]], dtype=float)) scharr_3x3 = SimpleNamespace( - x=np.array([[ -3, 0, 3], - [-10, 0, 10], - [ -3, 0, 3]], dtype=float), - y=np.array([[-3, -10, -3], - [ 0, 0, 0], - [ 3, 10, 3]], dtype=float)) + x=np.array([[-3, 0, 3], + [-10, 0, 10], + [-3, 0, 3]], dtype=float), + y=np.array([[-3, -10, -3], + [0, 0, 0], + [3, 10, 3]], dtype=float)) sobel_3x3 = SimpleNamespace( - x=np.array([[-1, 0, 1], - [-2, 0, 2], - [-1, 0, 1]], dtype=float), - y=np.array([[-1, -2, -1], - [ 0, 0, 0], - [ 1, 2, 1]], dtype=float)) + x=np.array([[-1, 0, 1], + [-2, 0, 2], + [-1, 0, 1]], dtype=float), + y=np.array([[-1, -2, -1], + [0, 0, 0], + [1, 2, 1]], dtype=float)) sobel_5x5 = SimpleNamespace( - x=np.array([[ -5, -4, 0, 4, 5], - [ -8,-10, 0, 10, 8], - [-10,-20, 0, 20, 10], - [ -8,-10, 0, 10, 8], - [ -5, -4, 0, 4, 5]], dtype=float), - y=np.array([[ -5, -8, -10, -8, -5], - [ -4,-10, -20,-10, 4], - [ 0, 0, 0, 0, 0], - [ 4, 10, 20, 10, 4], - [ 5, 8, 10, 8, 5]], dtype=float)) + x=np.array([[-5, -4, 0, 4, 5], + [-8, -10, 0, 10, 8], + [-10, -20, 0, 20, 10], + [-8, -10, 0, 10, 8], + [-5, -4, 0, 4, 5]], dtype=float), + y=np.array([[-5, -8, -10, -8, -5], + [-4, -10, -20, -10, 4], + [0, 0, 0, 0, 0], + [4, 10, 20, 10, 4], + [5, 8, 10, 8, 5]], dtype=float)) sobel_7x7 = SimpleNamespace( - x=np.array([[-3/18, -2/13, -1/10, 0, 1/10, 2/13, 3/18], - [-3/13, -2/8 , -1/5 , 0, 1/5 , 2/8 , 3/13], - [-3/10, -2/5 , -1/2 , 0, 1/2 , 2/5 , 3/10], - [-3/9 , -2/4 , -1/1 , 0, 1/1 , 2/4 , 3/9 ], - [-3/10, -2/5 , -1/2 , 0, 1/2 , 2/5 , 3/10], - [-3/13, -2/8 , -1/5 , 0, 1/5 , 2/8 , 3/13], - [-3/18, -2/13, -1/10, 0, 1/10, 2/13, 3/18]], dtype=float), - y=np.array([[-3/18, -3/13, -3/10, -3/9, -3/10, -3/13, -3/18], - [-2/13, -2/8 , -2/5 , -2/4, -2/5 , -2/8 , -2/13], - [-1/10, -1/5 , -1/2 , -1/1, -1/2 , -1/5 , -1/10], - [ 0, 0, 0 , 0, 0 , 0 , 0], - [ 1/10, 1/5 , 1/2 , 1/1, 1/2 , 1/5 , 1/10], - [ 2/13, 2/8 , 2/5 , 2/4, 2/5 , 2/8 , 2/13], - [ 3/18, 3/13, 3/10, 3/9, 3/10, 3/13, 3/18]], dtype=float)) + x=np.array([[-3 / 18, -2 / 13, -1 / 10, 0, 1 / 10, 2 / 13, 3 / 18], + [-3 / 13, -2 / 8, -1 / 5, 0, 1 / 5, 2 / 8, 3 / 13], + [-3 / 10, -2 / 5, -1 / 2, 0, 1 / 2, 2 / 5, 3 / 10], + [-3 / 9, -2 / 4, -1 / 1, 0, 1 / 1, 2 / 4, 3 / 9], + [-3 / 10, -2 / 5, -1 / 2, 0, 1 / 2, 2 / 5, 3 / 10], + [-3 / 13, -2 / 8, -1 / 5, 0, 1 / 5, 2 / 8, 3 / 13], + [-3 / 18, -2 / 13, -1 / 10, 0, 1 / 10, 2 / 13, 3 / 18]], dtype=float), + y=np.array([[-3 / 18, -3 / 13, -3 / 10, -3 / 9, -3 / 10, -3 / 13, -3 / 18], + [-2 / 13, -2 / 8, -2 / 5, -2 / 4, -2 / 5, -2 / 8, -2 / 13], + [-1 / 10, -1 / 5, -1 / 2, -1 / 1, -1 / 2, -1 / 5, -1 / 10], + [0, 0, 0, 0, 0, 0, 0], + [1 / 10, 1 / 5, 1 / 2, 1 / 1, 1 / 2, 1 / 5, 1 / 10], + [2 / 13, 2 / 8, 2 / 5, 2 / 4, 2 / 5, 2 / 8, 2 / 13], + [3 / 18, 3 / 13, 3 / 10, 3 / 9, 3 / 10, 3 / 13, 3 / 18]], dtype=float)) for kernel in [simple_3x3, prewitt_3x3, scharr_3x3, sobel_3x3, sobel_5x5, sobel_7x7]: kernel.x /= np.sum(np.abs(kernel.x)) kernel.y /= np.sum(np.abs(kernel.y)) + def get_kernel(kernel_name): if kernel_name.lower() in ['simple_3x3', 'simple']: return simple_3x3 @@ -83,72 +85,50 @@ def get_kernel(kernel_name): elif kernel_name.lower() in ['sobel_7x7']: return sobel_7x7 else: - warnings.warn(f'Deriviative name {kernel_name} is not implemented, ' - + 'using sobel filter instead.') + warnings.warn(f'Derivative name {kernel_name} is not implemented, ' + + 'using sobel filter instead.') return sobel_3x3 -def nan_conv2d(frame, kernel, kernel_center=None): - dx, dy = kernel.shape - dimx, dimy = frame.shape - dframe = np.empty((dimx, dimy))*np.nan - - if kernel_center is None: - kernel_center = [int((dim-1)/2) for dim in kernel.shape] - - # inverse kernel to mimic behavior or regular convolution algorithm - k = kernel[::-1, ::-1] - ci = dx - 1 - kernel_center[0] - cj = dy - 1 - kernel_center[1] +def nan_conv2d(image, kernel): + mask = 1 - np.isnan(image) + image_ = np.nan_to_num(image) + conve = convolve2d(image_, kernel, mode='same') + sum_ker = convolve2d(mask, kernel, mode='same') + pro = np.multiply(sum_ker, image) + quot = convolve2d(mask, abs(kernel), mode='same') + conve = np.divide((conve - pro), quot) + conve[np.isnan(image)] = np.nan + return conve - # loop over each frame site - for i,j in zip(*np.where(np.isfinite(frame))): - site = frame[i,j] - # loop over kernel window for frame site - window = np.zeros((dx,dy), dtype=float)*np.nan - for di,dj in itertools.product(range(dx), range(dy)): +def norm_angle(p): + # maps angle to [-pi, pi] + return -np.mod(p + np.pi, 2 * np.pi) + np.pi - # kernelsite != 0, framesite within borders and != nan - if k[di,dj] and 0 <= i+di-ci < dimx and 0 <= j+dj-cj < dimy \ - and np.isfinite(frame[i+di-ci,j+dj-cj]): - sign = -1*np.sign(k[di,dj]) - window[di,dj] = sign * (site - frame[i+di-ci,j+dj-cj]) - xi, yi = np.where(np.logical_not(np.isnan(window))) - if np.sum(np.logical_not(np.isnan(window))) > dx*dy/10: - dframe[i,j] = np.average(window[xi,yi], weights=abs(k[xi,yi])) - return dframe - -norm_angle = lambda p: -np.mod(p + np.pi, 2*np.pi) + np.pi - -def phase_conv2d(frame, kernel, kernel_center=None): +def phase_conv2d(image, kernel): + mask = 1 - np.isnan(image) + image = np.nan_to_num(image) + filters = [] dx, dy = kernel.shape - dimx, dimy = frame.shape - dframe = np.zeros_like(frame) - - if kernel_center is None: - kernel_center = [int((dim-1)/2) for dim in kernel.shape] - - # inverse kernel to mimic behavior or regular convolution algorithm - k = kernel[::-1, ::-1] - ci = dx - 1 - kernel_center[0] - cj = dy - 1 - kernel_center[1] - - # loop over kernel window for each frame site - for i,j in zip(*np.where(np.isfinite(frame))): - phase = frame[i,j] - dphase = np.zeros((dx,dy), dtype=float) - - for di,dj in itertools.product(range(dx), range(dy)): - - # kernelsite != 0, framesite within borders and != nan - if k[di,dj] and i+di-ci < dimx and j+dj-cj < dimy \ - and np.isfinite(frame[i+di-ci,j+dj-cj]): - sign = -1*np.sign(k[di,dj]) - # pos = clockwise from phase to frame[..] - dphase[di,dj] = sign*norm_angle(phase-frame[i+di-ci,j+dj-cj]) - - if dphase.any(): - dframe[i,j] = np.average(dphase, weights=abs(k)) / np.pi - return dframe + ci, cj = [int((dim - 1) / 2) for dim in kernel.shape] + for di, dj in itertools.product(range(dx), range(dy)): + if (di, dj) != (ci, cj): + a = np.zeros((dx, dy)) + a[di, dj] = 1 + filters.append(a) + n_channels = len(filters) + dimx, dimy = image.shape + conv = np.zeros((n_channels, dimx, dimy)) + for i in range(n_channels): + channel = convolve2d(image, filters[i], mode='same') + mean_ = convolve2d(mask, filters[i], mode='same') * image + conv[i] = channel - mean_ + conv = norm_angle(conv) + for i in range(n_channels): + conv[i] = conv[i] * np.sum(filters[i] * kernel) + conv = np.sum(conv, axis=0) + quot = np.sum(abs(kernel)) * np.pi + out = conv / quot + return out * mask