diff --git a/scilpy/denoise/angle_aware_bilateral.cl b/scilpy/denoise/angle_aware_bilateral.cl deleted file mode 100644 index 63f8702d4..000000000 --- a/scilpy/denoise/angle_aware_bilateral.cl +++ /dev/null @@ -1,131 +0,0 @@ -/* -OpenCL kernel code for computing angle-aware bilateral filtering. -*/ - -// Compiler definitions with placeholder values -#define IN_N_COEFFS 0 -#define OUT_N_COEFFS 0 -#define N_DIRS 0 -#define SIGMA_RANGE 0 -#define IM_X_DIM 0 -#define IM_Y_DIM 0 -#define IM_Z_DIM 0 -#define H_X_DIM 0 -#define H_Y_DIM 0 -#define H_Z_DIM 0 - -int get_flat_index(const int x, const int y, - const int z, const int w, - const int xLen, - const int yLen, - const int zLen) -{ - return x + - y * xLen + - z * xLen * yLen + - w * xLen * yLen * zLen; -} - -float sf_for_direction(const int idx, const int idy, const int idz, - const int dir_id, global const float* sh_buffer, - global const float* sh_to_sf_mat) -{ - float sf_coeff = 0.0f; - for(int i = 0; i < IN_N_COEFFS; ++i) - { - const int im_index = get_flat_index(idx, idy, idz, i, - IM_X_DIM + H_X_DIM - 1, - IM_Y_DIM + H_Y_DIM - 1, - IM_Z_DIM + H_Z_DIM - 1); - const float ylmu = sh_to_sf_mat[get_flat_index(i, dir_id, 0, 0, - IN_N_COEFFS, - N_DIRS, 0)]; - sf_coeff += sh_buffer[im_index] * ylmu; - } - return sf_coeff; -} - -void sf_to_sh(const float* sf_coeffs, global const float* sf_to_sh_mat, - float* sh_coeffs) -{ - // although vector operations are supported - // in OpenCL, maximum n value is of 16. - for(int i = 0; i < OUT_N_COEFFS; ++i) - { - sh_coeffs[i] = 0.0f; - for(int u = 0; u < N_DIRS; ++u) - { - const float ylmu = sf_to_sh_mat[get_flat_index(u, i, 0, 0, - N_DIRS, - OUT_N_COEFFS, - 0)]; - sh_coeffs[i] += sf_coeffs[u] * ylmu; - } - } -} - -float range_distribution(const float iVal, const float jVal) -{ - const float x = fabs(iVal - jVal); - return exp(-pow(x, 2)/2.0f/pown((float)SIGMA_RANGE, 2)) - / SIGMA_RANGE / sqrt(2.0f * M_PI); -} - -__kernel void correlate(__global const float *sh_buffer, - __global const float *h_weights, - __global const float *sh_to_sf_mat, - __global const float *sf_to_sh_mat, - __global float *out_sh_buffer) -{ - const int idx = get_global_id(0); - const int idy = get_global_id(1); - const int idz = get_global_id(2); - - float sf_coeffs[N_DIRS]; - float norm_w; - for(int u = 0; u < N_DIRS; ++u) - { - sf_coeffs[u] = 0.0f; - const float sf_center = sf_for_direction(idx + H_X_DIM / 2, - idy + H_Y_DIM / 2, - idz + H_Z_DIM / 2, - u, sh_buffer, - sh_to_sf_mat); - norm_w = 0.0f; - for(int hx = 0; hx < H_X_DIM; ++hx) - { - for(int hy = 0; hy < H_Y_DIM; ++hy) - { - for(int hz = 0; hz < H_Z_DIM; ++hz) - { - const float sf_u = - sf_for_direction(idx + hx, idy + hy, - idz + hz, u, sh_buffer, - sh_to_sf_mat); - - const float range_w = - range_distribution(sf_center, sf_u); - - const float weight = - h_weights[get_flat_index(hx, hy, hz, u, - H_X_DIM, H_Y_DIM, - H_Z_DIM)] * range_w; - - sf_coeffs[u] += sf_u * weight; - norm_w += weight; - } - } - } - sf_coeffs[u] /= norm_w; - } - - float sh_coeffs[OUT_N_COEFFS]; - sf_to_sh(sf_coeffs, sf_to_sh_mat, sh_coeffs); - for(int i = 0; i < OUT_N_COEFFS; ++i) - { - const int out_index = get_flat_index(idx, idy, idz, i, - IM_X_DIM, IM_Y_DIM, - IM_Z_DIM); - out_sh_buffer[out_index] = sh_coeffs[i]; - } -} diff --git a/scilpy/denoise/aodf_filter.cl b/scilpy/denoise/aodf_filter.cl new file mode 100644 index 000000000..4cf707e02 --- /dev/null +++ b/scilpy/denoise/aodf_filter.cl @@ -0,0 +1,113 @@ +#define WIN_WIDTH 0 +#define SIGMA_RANGE 1.0f +#define N_DIRS 200 +#define EXCLUDE_SELF false +#define DISABLE_ANGLE false +#define DISABLE_RANGE false + +int get_flat_index(const int x, const int y, const int z, + const int w, const int x_len, + const int y_len, const int z_len) +{ + return x + y * x_len + z * x_len * y_len + + w * x_len * y_len * z_len; +} + +float range_filter(const float x) +{ + return exp(-pow(x, 2)/2.0f/pown((float)SIGMA_RANGE, 2)); +} + +// SF data is padded while out_sf isn't +__kernel void filter(__global const float* sf_data, + __global const float* nx_filter, + __global const float* uv_filter, + __global const float* uv_weights_offset, + __global const float* v_indices, + __global float* out_sf) +{ + // *output* image dimensions + const int x_len = get_global_size(0); + const int y_len = get_global_size(1); + const int z_len = get_global_size(2); + + // padded dimensions + const int x_pad_len = x_len + WIN_WIDTH - 1; + const int y_pad_len = y_len + WIN_WIDTH - 1; + const int z_pad_len = z_len + WIN_WIDTH - 1; + + // output voxel indice + const int x_ind = get_global_id(0); + const int y_ind = get_global_id(1); + const int z_ind = get_global_id(2); + + // window half width + const int win_hwidth = WIN_WIDTH / 2; + + // process each direction inside voxel + for(int ui = 0; ui < N_DIRS; ++ui) + { + // in input volume, dimensions are padded + const int xui_flat_ind = get_flat_index(x_ind + win_hwidth, + y_ind + win_hwidth, + z_ind + win_hwidth, + ui, x_pad_len, y_pad_len, + z_pad_len); + const float psi_xui = sf_data[xui_flat_ind]; + + // output value + float w_norm = 0.0f; + float tilde_psi_xui = 0.0f; +#if DISABLE_ANGLE + const int vi_in_image = ui; +#else + const int vi_offset = (int)uv_weights_offset[ui]; + const int n_ang_neigbrs = (int)uv_weights_offset[ui + 1] - vi_offset; + for(int dir_i = 0; dir_i < n_ang_neigbrs; ++dir_i) + { + const int vi_in_filter = vi_offset + dir_i; + const int vi_in_image = v_indices[vi_in_filter]; +#endif + for(int hi = 0; hi < WIN_WIDTH; ++hi) + { + for(int hj = 0; hj < WIN_WIDTH; ++hj) + { + for(int hk = 0; hk < WIN_WIDTH; ++hk) + { + const int yvi_flat_ind = + get_flat_index(x_ind + hi, y_ind + hj, z_ind + hk, + vi_in_image, x_pad_len, y_pad_len, + z_pad_len); + const float psi_yvi = sf_data[yvi_flat_ind]; +#if DISABLE_RANGE + const float r_weight = 1.0f; +#else + const float r_weight = range_filter(fabs(psi_xui - psi_yvi)); +#endif + // contains "align" weight, so direction is ui + const int y_in_nx_flat_ind = get_flat_index(hi, hj, hk, ui, + WIN_WIDTH, WIN_WIDTH, + WIN_WIDTH); + const float nx_weight = nx_filter[y_in_nx_flat_ind]; + +#if DISABLE_ANGLE + const float uv_weight = 1.0f; +#else + const float uv_weight = uv_filter[vi_in_filter]; +#endif + + const float res_weight_yvi = nx_weight * r_weight * uv_weight; + tilde_psi_xui += res_weight_yvi * psi_yvi; + w_norm += res_weight_yvi; + } + } + } +#if !DISABLE_ANGLE + } +#endif + // normalize and assign + const int output_flat_ind = get_flat_index(x_ind, y_ind, z_ind, ui, + x_len, y_len, z_len); + out_sf[output_flat_ind] = tilde_psi_xui / w_norm; + } +} \ No newline at end of file diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py index 534da1683..e63977e4d 100644 --- a/scilpy/denoise/asym_filtering.py +++ b/scilpy/denoise/asym_filtering.py @@ -1,387 +1,491 @@ # -*- coding: utf-8 -*- - import numpy as np +import logging from dipy.reconst.shm import sh_to_sf_matrix from dipy.data import get_sphere from dipy.core.sphere import Sphere from scipy.ndimage import correlate +from itertools import product as iterprod from scilpy.gpuparallel.opencl_utils import have_opencl, CLKernel, CLManager -def angle_aware_bilateral_filtering(in_sh, sh_order=8, - sh_basis='descoteaux07', - in_full_basis=False, - is_legacy=True, - sphere_str='repulsion724', - sigma_spatial=1.0, sigma_angular=1.0, - sigma_range=0.5, use_gpu=True): +def unified_filtering(sh_data, sh_order, sh_basis, is_legacy, full_basis, + sphere_str, sigma_spatial=1.0, sigma_align=0.8, + sigma_angle=None, rel_sigma_range=0.2, + win_hwidth=None, exclude_center=False, + device_type='gpu', use_opencl=True, patch_size=40): """ - Angle-aware bilateral filtering. + Unified asymmetric filtering as described in [1]. Parameters ---------- - in_sh: ndarray (x, y, z, ncoeffs) - Input SH volume. - sh_order: int, optional - Maximum SH order of input volume. - sh_basis: str, optional - Name of SH basis used. - in_full_basis: bool, optional - True if input is expressed in full SH basis. - is_legacy : bool, optional - Whether or not the SH basis is in its legacy form. - sphere_str: str, optional - Name of the DIPY sphere to use for sh to sf projection. - sigma_spatial: float, optional - Standard deviation for spatial filter. - sigma_angular: float, optional - Standard deviation for angular filter. - sigma_range: float, optional - Standard deviation for range filter. - use_gpu: bool, optional - True if GPU should be used. - - Returns - ------- - out_sh: ndarray (x, y, z, ncoeffs) - Output SH coefficient array in full SH basis. + sh_data: ndarray + SH coefficients image. + sh_order: int + Maximum order of spherical harmonics (SH) basis. + sh_basis: str + SH basis definition used for input and output SH image. + One of 'descoteaux07' or 'tournier07'. + is_legacy: bool + Whether the legacy SH basis definition should be used. + full_basis: bool + Whether the input SH basis is full or not. + sphere_str: str + Name of the DIPY sphere to use for SH to SF projection. + sigma_spatial: float or None + Standard deviation of spatial filter. Can be None to replace + by mean filter, in what case win_hwidth must be given. + sigma_align: float or None + Standard deviation of alignment filter. `None` disables + alignment filtering. + sigma_angle: float or None + Standard deviation of the angle filter. `None` disables + angle filtering. + rel_sigma_range: float or None + Standard deviation of the range filter, relative to the + range of SF amplitudes. `None` disables range filtering. + disable_spatial: bool, optional + Replace gaussian filter by a mean filter for spatial filter. + The value from `sigma_spatial` is still used for setting the + size of the filtering window. + win_hwidth: int, optional + Half-width of the filtering window. When None, the + filtering window half-width is given by (6*sigma_spatial + 1). + exclude_center: bool, optional + Assign a weight of 0 to the center voxel of the filter. + device_type: string, optional + Device on which the code should run. Choices are cpu or gpu. + use_opencl: bool, optional + Use OpenCL for software acceleration. + patch_size: int, optional + Patch size for OpenCL execution. + + References + ---------- + [1] Poirier and Descoteaux, 2024, "A Unified Filtering Method for + Estimating Asymmetric Orientation Distribution Functions", + Neuroimage, https://doi.org/10.1016/j.neuroimage.2024.120516 """ - if use_gpu and have_opencl: - return angle_aware_bilateral_filtering_gpu(in_sh, sh_order, - sh_basis, in_full_basis, - is_legacy, - sphere_str, sigma_spatial, - sigma_angular, sigma_range) - elif use_gpu and not have_opencl: - raise RuntimeError('Package pyopencl not found. Install pyopencl' - ' or set use_gpu to False.') + if sigma_spatial is None and win_hwidth is None: + raise ValueError('sigma_spatial and win_hwidth cannot both be None') + if device_type not in ['cpu', 'gpu']: + raise ValueError('Invalid device type {}. Must be cpu or gpu' + .format(device_type)) + if use_opencl and not have_opencl: + raise ValueError('pyopencl is not installed. Please install before' + ' using option use_opencl=True.') + if device_type == 'gpu' and not use_opencl: + raise ValueError('Option use_opencl must be enabled ' + 'to use device \'gpu\'.') + + sphere = get_sphere(sphere_str) + + if sigma_spatial is not None: + if sigma_spatial <= 0.0: + raise ValueError('sigma_spatial cannot be <= 0.') + # calculate half-width from sigma_spatial + half_width = int(round(3*sigma_spatial)) + if sigma_align is not None: + if sigma_align <= 0.0: + raise ValueError('sigma_align cannot be <= 0.') + if sigma_angle is not None: + if sigma_angle <= 0.0: + raise ValueError('sigma_align cannot be <= 0.') + + # overwrite half-width if win_hwidth is supplied + if win_hwidth is not None: + half_width = win_hwidth + + # filter shape computed from half_width + filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1) + + # build filters + uv_filter = _unified_filter_build_uv(sigma_angle, sphere) + nx_filter = _unified_filter_build_nx(filter_shape, sigma_spatial, + sigma_align, sphere, exclude_center) + + B = sh_to_sf_matrix(sphere, sh_order, sh_basis, full_basis, + legacy=is_legacy, return_inv=False) + _, B_inv = sh_to_sf_matrix(sphere, sh_order, sh_basis, True, + legacy=is_legacy, return_inv=True) + + # compute "real" sigma_range scaled by sf amplitudes + # if rel_sigma_range is supplied + sigma_range = None + if rel_sigma_range is not None: + if rel_sigma_range <= 0.0: + raise ValueError('sigma_rangel cannot be <= 0.') + sigma_range = rel_sigma_range * _get_sf_range(sh_data, B) + + if use_opencl: + # initialize opencl + cl_manager = _unified_filter_prepare_opencl(sigma_range, sigma_angle, + filter_shape[0], sphere, + device_type) + + return _unified_filter_call_opencl(sh_data, nx_filter, uv_filter, + cl_manager, B, B_inv, sphere, + patch_size) else: - return angle_aware_bilateral_filtering_cpu(in_sh, sh_order, - sh_basis, in_full_basis, - is_legacy, - sphere_str, sigma_spatial, - sigma_angular, sigma_range) - - -def angle_aware_bilateral_filtering_gpu(in_sh, sh_order=8, - sh_basis='descoteaux07', - in_full_basis=False, - is_legacy=True, - sphere_str='repulsion724', - sigma_spatial=1.0, - sigma_angular=1.0, - sigma_range=0.5): + return _unified_filter_call_python(sh_data, nx_filter, uv_filter, + sigma_range, B, B_inv, sphere) + + +def _unified_filter_prepare_opencl(sigma_range, sigma_angle, window_width, + sphere, device_type): """ - Angle-aware bilateral filtering using OpenCL for GPU computing. + Instantiate OpenCL context manager and compile OpenCL program. Parameters ---------- - in_sh: ndarray (x, y, z, ncoeffs) - Input SH volume. - sh_order: int, optional - Maximum SH order of input volume. - sh_basis: str, optional - Name of SH basis used. - in_full_basis: bool, optional - True if input is expressed in full SH basis. - is_legacy : bool, optional - Whether or not the SH basis is in its legacy form. - sphere_str: str, optional - Name of the DIPY sphere to use for sh to sf projection. - sigma_spatial: float, optional - Standard deviation for spatial filter. - sigma_angular: float, optional - Standard deviation for angular filter. - sigma_range: float, optional - Standard deviation for range filter. + sigma_range: float or None + Value for sigma_range. + sigma_angle: float or None + Value for sigma_angle. + window_width: int + Width of filtering window. + sphere: DIPY sphere + Sphere used for SH to SF projection. + device_type: string + Device to be used by OpenCL. Either 'cpu' or 'gpu'. Returns ------- - out_sh: ndarray (x, y, z, ncoeffs) - Output SH coefficient array in full SH basis. + cl_manager: CLManager + OpenCL manager object. """ - s_weights = _get_spatial_weights(sigma_spatial) - h_half_width = len(s_weights) // 2 + disable_range = sigma_range is None + disable_angle = sigma_angle is None + if sigma_range is None: + sigma_range = 0.0 # placeholder value for sigma_range - sphere = get_sphere(sphere_str) - a_weights = _get_angular_weights(s_weights.shape, sphere, sigma_angular) - - h_weights = s_weights[..., None] * a_weights - h_weights /= np.sum(h_weights, axis=(0, 1, 2)) - - sh_to_sf_mat = sh_to_sf_matrix(sphere, sh_order_max=sh_order, - basis_type=sh_basis, - full_basis=in_full_basis, - legacy=is_legacy, - return_inv=False) - - _, sf_to_sh_mat = sh_to_sf_matrix(sphere, sh_order_max=sh_order, - basis_type=sh_basis, - full_basis=True, - legacy=is_legacy, - return_inv=True) - - out_n_coeffs = sf_to_sh_mat.shape[1] - n_dirs = len(sphere.vertices) - volume_shape = in_sh.shape - in_sh = np.pad(in_sh, ((h_half_width, h_half_width), - (h_half_width, h_half_width), - (h_half_width, h_half_width), - (0, 0))) - - cl_kernel = CLKernel('correlate', 'denoise', 'angle_aware_bilateral.cl') - cl_kernel.set_define('IM_X_DIM', volume_shape[0]) - cl_kernel.set_define('IM_Y_DIM', volume_shape[1]) - cl_kernel.set_define('IM_Z_DIM', volume_shape[2]) - - cl_kernel.set_define('H_X_DIM', h_weights.shape[0]) - cl_kernel.set_define('H_Y_DIM', h_weights.shape[1]) - cl_kernel.set_define('H_Z_DIM', h_weights.shape[2]) - - cl_kernel.set_define('SIGMA_RANGE', float(sigma_range)) - - cl_kernel.set_define('IN_N_COEFFS', volume_shape[-1]) - cl_kernel.set_define('OUT_N_COEFFS', out_n_coeffs) - cl_kernel.set_define('N_DIRS', n_dirs) - - cl_manager = CLManager(cl_kernel, 4, 1) - cl_manager.add_input_buffer(0, in_sh) - cl_manager.add_input_buffer(1, h_weights) - cl_manager.add_input_buffer(2, sh_to_sf_mat) - cl_manager.add_input_buffer(3, sf_to_sh_mat) - - cl_manager.add_output_buffer(0, volume_shape[:3] + (out_n_coeffs,), - np.float32) - - outputs = cl_manager.run(volume_shape[:3]) - return outputs[0] - - -def angle_aware_bilateral_filtering_cpu(in_sh, sh_order=8, - sh_basis='descoteaux07', - in_full_basis=False, - is_legacy=True, - sphere_str='repulsion724', - sigma_spatial=1.0, - sigma_angular=1.0, - sigma_range=0.5): + cl_kernel = CLKernel('filter', 'denoise', 'aodf_filter.cl') + cl_kernel.set_define('WIN_WIDTH', window_width) + cl_kernel.set_define('SIGMA_RANGE', '{}f'.format(sigma_range)) + cl_kernel.set_define('N_DIRS', len(sphere.vertices)) + cl_kernel.set_define('DISABLE_ANGLE', 'true' if disable_angle else 'false') + cl_kernel.set_define('DISABLE_RANGE', 'true' if disable_range else 'false') + + return CLManager(cl_kernel, device_type) + + +def _unified_filter_build_uv(sigma_angle, sphere): """ - Angle-aware bilateral filtering on the CPU - (optionally using multiple threads). + Build the angle filter, weighted on angle between current direction u + and neighbour direction v. Parameters ---------- - in_sh: ndarray (x, y, z, ncoeffs) - Input SH volume. - sh_order: int, optional - Maximum SH order of input volume. - sh_basis: str, optional - Name of SH basis used. - in_full_basis: bool, optional - True if input is expressed in full SH basis. - is_legacy : bool, optional - Whether or not the SH basis is in its legacy form. - sphere_str: str, optional - Name of the DIPY sphere to use for sh to sf projection. - sigma_spatial: float, optional - Standard deviation for spatial filter. - sigma_angular: float, optional - Standard deviation for angular filter. - sigma_range: float, optional - Standard deviation for range filter. + sigma_angle: float + Standard deviation of filter. Values at distances greater than + sigma_angle are clipped to 0 to reduce computation time. + sphere: DIPY sphere + Sphere used for sampling the SF. Returns ------- - out_sh: ndarray (x, y, z, ncoeffs) - Output SH coefficient array in full SH basis. + weights: ndarray + Angle filter of shape (N_dirs, N_dirs). """ - # Load the sphere used for projection of SH - sphere = get_sphere(sphere_str) - - # Normalized filter for each sf direction - s_weights = _get_spatial_weights(sigma_spatial) - a_weights = _get_angular_weights(s_weights.shape, sphere, sigma_angular) - - weights = s_weights[..., None] * a_weights - weights /= np.sum(weights, axis=(0, 1, 2)) - - nb_sf = len(sphere.vertices) - B = sh_to_sf_matrix(sphere, sh_order_max=sh_order, basis_type=sh_basis, - return_inv=False, full_basis=in_full_basis, - legacy=is_legacy) - - mean_sf = np.zeros(in_sh.shape[:-1] + (nb_sf,)) - - # Apply filter to each sphere vertice - for sph_id in range(nb_sf): - w_filter = weights[..., sph_id] - - # Generate 1-channel images for directions u and -u - current_sf = np.dot(in_sh, B[:, sph_id]) - mean_sf[..., sph_id] = _correlate_spatial(current_sf, w_filter, - sigma_range) - - # Convert back to SH coefficients - _, B_inv = sh_to_sf_matrix(sphere, sh_order_max=sh_order, basis_type=sh_basis, - full_basis=True, legacy=is_legacy) - out_sh = np.array([np.dot(i, B_inv) for i in mean_sf], dtype=in_sh.dtype) - # By default, return only asymmetric SH - return out_sh + directions = sphere.vertices + if sigma_angle is not None: + dot = directions.dot(directions.T) + x = np.arccos(np.clip(dot, -1.0, 1.0)) + weights = _evaluate_gaussian_distribution(x, sigma_angle) + mask = x > (3.0*sigma_angle) + weights[mask] = 0.0 + weights /= np.sum(weights, axis=-1) + else: + weights = np.eye(len(directions)) + return weights -def _evaluate_gaussian_distribution(x, sigma): +def _unified_filter_build_nx(filter_shape, sigma_spatial, sigma_align, + sphere, exclude_center): """ - 1-dimensional 0-centered Gaussian distribution - with standard deviation sigma. + Build the combined spatial and alignment filter. Parameters ---------- - x: ndarray or float - Points where the distribution is evaluated. - sigma: float - Standard deviation. + filter_shape: tuple + Dimensions of filtering window. + sigma_spatial: float or None + Standard deviation of spatial filter. None disables Gaussian + weighting for spatial filtering. + sigma_align: float or None + Standard deviation of the alignment filter. None disables Gaussian + weighting for alignment filtering. + sphere: DIPY sphere + Sphere for SH to SF projection. + exclude_center: bool + Whether the center voxel is included in the neighbourhood. Returns ------- - out: ndarray or float - Values at x. + weights: ndarray + Combined spatial + alignment filter of shape (W, H, D, N) where + N is the number of sphere directions. """ - assert sigma > 0.0, "Sigma must be greater than 0." - cnorm = 1.0 / sigma / np.sqrt(2.0*np.pi) - return cnorm * np.exp(-x**2/2/sigma**2) + directions = sphere.vertices.astype(np.float32) + grid_directions = _get_window_directions(filter_shape).astype(np.float32) + distances = np.linalg.norm(grid_directions, axis=-1) + grid_directions[distances > 0] = grid_directions[distances > 0] /\ + distances[distances > 0][..., None] -def _get_window_directions(shape): + if sigma_spatial is None: + w_spatial = np.ones(filter_shape) + else: + w_spatial = _evaluate_gaussian_distribution(distances, sigma_spatial) + + if sigma_align is None: + w_align = np.ones(np.append(filter_shape, (len(directions),))) + else: + cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0) + theta = np.arccos(cos_theta) + theta[filter_shape[0] // 2, + filter_shape[1] // 2, + filter_shape[2] // 2] = 0.0 + w_align = _evaluate_gaussian_distribution(theta, sigma_align) + + # resulting filter + w = w_spatial[..., None] * w_align + + if exclude_center: + w[filter_shape[0] // 2, + filter_shape[1] // 2, + filter_shape[2] // 2] = 0.0 + + # normalize and return + w /= np.sum(w, axis=(0, 1, 2), keepdims=True) + return w + + +def _get_sf_range(sh_data, B_mat): """ - Get directions from center voxel to all neighbours - for a window of given shape. + Get the range of SF amplitudes for input `sh_data`. Parameters ---------- - shape: tuple - Dimensions of the window. + sh_data: ndarray + Spherical harmonics coefficients image. + B_mat: ndarray + SH to SF projection matrix. Returns ------- - grid: ndarray - Grid containing the direction from the center voxel to - the current position for all positions inside the window. + sf_range: float + Range of SF amplitudes. """ - grid = np.indices(shape) - grid = np.moveaxis(grid, 0, -1) - grid = grid - np.asarray(shape) // 2 - return grid + sf = np.array([np.dot(i, B_mat) for i in sh_data], + dtype=sh_data.dtype) + sf[sf < 0.0] = 0.0 + sf_max = np.max(sf) + sf_min = np.min(sf) + return sf_max - sf_min -def _get_spatial_weights(sigma_spatial): +def _unified_filter_call_opencl(sh_data, nx_filter, uv_filter, cl_manager, + B, B_inv, sphere, patch_size=40): """ - Compute the spatial filter, which is an isotropic Gaussian filter - of standard deviation sigma_spatial forweighting by the distance - between voxel positions. The size of the filter is given by - 6 * sigma_spatial, in order to cover the range - [-3*sigma_spatial, 3*sigma_spatial]. + Run unified filtering for asymmetric ODFs using OpenCL. Parameters ---------- - sigma_spatial: float - Standard deviation of spatial filter. + sh_data: ndarray + Input SH volume. + nx_filter: ndarray + Combined spatial and alignment filter. + uv_filter: ndarray + Angle filter. + cl_manager: CLManager + A CLManager instance. + B: ndarray + SH to SF projection matrix. + B_inv: ndarray + SF to SH projection matrix. + sphere: DIPY sphere + Sphere for SH to SF projection. + patch_size: int + Data is processed in patches of + patch_size x patch_size x patch_size. Returns ------- - spatial_weights: ndarray - Spatial filter. + out_sh: ndarray + Filtered output as SH coefficients. """ - shape = int(6 * sigma_spatial) - if shape % 2 == 0: - shape += 1 - shape = (shape, shape, shape) - - grid = _get_window_directions(shape) + uv_weights_offsets =\ + np.append([0.0], np.cumsum(np.count_nonzero(uv_filter, + axis=-1))) + v_indices = np.tile(np.arange(uv_filter.shape[1]), + (uv_filter.shape[0], 1))[uv_filter > 0.0] + flat_uv = uv_filter[uv_filter > 0.0] + + # Prepare GPU buffers + cl_manager.add_input_buffer("sf_data") # SF data not initialized + cl_manager.add_input_buffer("nx_filter", nx_filter) + cl_manager.add_input_buffer("uv_filter", flat_uv) + cl_manager.add_input_buffer("uv_weights_offsets", uv_weights_offsets) + cl_manager.add_input_buffer("v_indices", v_indices) + cl_manager.add_output_buffer("out_sf") # SF not initialized yet + + win_width = nx_filter.shape[0] + win_hwidth = win_width // 2 + volume_shape = sh_data.shape[:-1] + padded_volume_shape = tuple(np.asarray(volume_shape) + win_width - 1) + + out_sh = np.zeros(np.append(sh_data.shape[:-1], B_inv.shape[-1])) + # Pad SH data + sh_data = np.pad(sh_data, ((win_hwidth, win_hwidth), + (win_hwidth, win_hwidth), + (win_hwidth, win_hwidth), + (0, 0))) + + # process in batches + padded_patch_size = patch_size + nx_filter.shape[0] - 1 + + n_splits = np.ceil(np.asarray(volume_shape) / float(patch_size))\ + .astype(int) + splits_prod = iterprod(np.arange(n_splits[0]), + np.arange(n_splits[1]), + np.arange(n_splits[2])) + n_splits_prod = np.prod(n_splits) + for i, split_offset in enumerate(splits_prod): + logging.info('Patch {}/{}'.format(i+1, n_splits_prod)) + i, j, k = split_offset + patch_in = np.array( + [[i * patch_size, min((i*patch_size)+padded_patch_size, + padded_volume_shape[0])], + [j * patch_size, min((j*patch_size)+padded_patch_size, + padded_volume_shape[1])], + [k * patch_size, min((k*patch_size)+padded_patch_size, + padded_volume_shape[2])]]) + patch_out = np.array( + [[i * patch_size, min((i+1)*patch_size, volume_shape[0])], + [j * patch_size, min((j+1)*patch_size, volume_shape[1])], + [k * patch_size, min((k+1)*patch_size, volume_shape[2])]]) + out_shape = tuple(np.append(patch_out[:, 1] - patch_out[:, 0], + len(sphere.vertices))) + + sh_patch = sh_data[patch_in[0, 0]:patch_in[0, 1], + patch_in[1, 0]:patch_in[1, 1], + patch_in[2, 0]:patch_in[2, 1]] + + sf_patch = np.dot(sh_patch, B) + cl_manager.update_input_buffer("sf_data", sf_patch) + cl_manager.update_output_buffer("out_sf", out_shape) + out_sf = cl_manager.run(out_shape[:-1])[0] + out_sh[patch_out[0, 0]:patch_out[0, 1], + patch_out[1, 0]:patch_out[1, 1], + patch_out[2, 0]:patch_out[2, 1]] = np.dot(out_sf, B_inv) - distances = np.linalg.norm(grid, axis=-1) - spatial_weights = _evaluate_gaussian_distribution(distances, sigma_spatial) - - # normalize filter - spatial_weights /= np.sum(spatial_weights) - return spatial_weights + return out_sh -def _get_angular_weights(shape, sphere, sigma_angular): +def _unified_filter_call_python(sh_data, nx_filter, uv_filter, sigma_range, + B_mat, B_inv, sphere): """ - Compute the angular filter, weighted by the alignment between a - sphere direction and the direction to a neighbour. The parameter - sigma_angular controls the sharpness of the kernel. + Run filtering using pure python implementation. Parameters ---------- - shape: tuple - Shape of the angular filter. - sphere: dipy Sphere - Sphere on which the SH coefficeints are projected. - sigma_angular: float - Standard deviation of Gaussian distribution. + sh_data: ndarray + Input SH data. + nx_filter: ndarray + Combined spatial and alignment filter. + uv_filter: ndarray + Angle filter. + sigma_range: float or None + Standard deviation of range filter. None disables range filtering. + B_mat: ndarray + SH to SF projection matrix. + B_inv: ndarray + SF to SH projection matrix. + sphere: DIPY sphere + Sphere for SH to SF projection. Returns ------- - angular_weights: ndarray - Angular filter for each position and for each sphere direction. + out_sh: ndarray + Filtered output as SH coefficients. """ - grid_dirs = _get_window_directions(shape).astype(np.float32) - dir_norms = np.linalg.norm(grid_dirs, axis=-1) - - # normalized grid directions - grid_dirs[dir_norms > 0] /= dir_norms[dir_norms > 0][:, None] - angles = np.arccos(np.dot(grid_dirs, sphere.vertices.T)) - angles[np.logical_not(dir_norms > 0), :] = 0.0 - - angular_weights = _evaluate_gaussian_distribution(angles, sigma_angular) + nb_sf = len(sphere.vertices) + mean_sf = np.zeros(sh_data.shape[:-1] + (nb_sf,)) - # normalize filter per direction - angular_weights /= np.sum(angular_weights, axis=(0, 1, 2)) - return angular_weights + # Apply filter to each sphere vertice + for u_sph_id in range(nb_sf): + if u_sph_id % 20 == 0: + logging.info('Processing direction: {}/{}' + .format(u_sph_id, nb_sf)) + mean_sf[..., u_sph_id] = _correlate(sh_data, nx_filter, uv_filter, + sigma_range, u_sph_id, B_mat) + + out_sh = np.array([np.dot(i, B_inv) for i in mean_sf], + dtype=sh_data.dtype) + return out_sh -def _correlate_spatial(image_u, h_filter, sigma_range): +def _correlate(sh_data, nx_filter, uv_filter, sigma_range, u_index, B_mat): """ - Implementation of the correlation operation for anisotropic filtering. + Apply the filters to the SH image for the sphere direction + described by `u_index`. Parameters ---------- - image_u: ndarray (X, Y, Z) - SF image for some sphere direction u. - h_filter: ndarray (W, H, D) - 3-dimensional filter to apply. - sigma_range: float - Standard deviation of the Gaussian distribution defining the - range kernel. + sh_data: ndarray + Input SH coefficients. + nx_filter: ndarray + Combined spatial and alignment filter. + uv_filter: ndarray + Angle filter. + sigma_range: float or None + Standard deviation of range filter. None disables range filtering. + u_index: int + Index of the current sphere direction to process. + B_mat: ndarray + SH to SF projection matrix. Returns ------- - out_im: ndarray (X, Y, Z) - Filtered SF image. + out_sf: ndarray + Output SF amplitudes along the direction described by `u_index`. """ - h_w, h_h, h_d = h_filter.shape[:3] + v_indices = np.flatnonzero(uv_filter[u_index]) + nx_filter = nx_filter[..., u_index] + h_w, h_h, h_d = nx_filter.shape[:3] half_w, half_h, half_d = h_w // 2, h_h // 2, h_d // 2 - out_im = np.zeros_like(image_u) - image_u = np.pad(image_u, ((half_w, half_w), + out_sf = np.zeros(sh_data.shape[:3]) + sh_data = np.pad(sh_data, ((half_w, half_w), (half_h, half_h), - (half_d, half_d))) + (half_d, half_d), + (0, 0))) + + sf_u = np.dot(sh_data, B_mat[:, u_index]) + sf_v = np.dot(sh_data, B_mat[:, v_indices]) + uv_filter = uv_filter[u_index, v_indices] - for ii in range(out_im.shape[0]): - for jj in range(out_im.shape[1]): - for kk in range(out_im.shape[2]): - x = image_u[ii:ii+h_w, jj:jj+h_h, kk:kk+h_d]\ - - image_u[ii, jj, kk] - range_filter = _evaluate_gaussian_distribution(x, sigma_range) - res_filter = range_filter * h_filter + _get_range = _evaluate_gaussian_distribution\ + if sigma_range is not None else lambda x, _: np.ones_like(x) - out_im[ii, jj, kk] += np.sum(image_u[ii:ii+h_w, - jj:jj+h_h, - kk:kk+h_d] - * res_filter) - out_im[ii, jj, kk] /= np.sum(res_filter) + for ii in range(out_sf.shape[0]): + for jj in range(out_sf.shape[1]): + for kk in range(out_sf.shape[2]): + a = sf_v[ii:ii+h_w, jj:jj+h_h, kk:kk+h_d] + b = sf_u[ii + half_w, jj + half_h, kk + half_d] + x_range = a - b + range_filter = _get_range(x_range, sigma_range) - return out_im + # the resulting filter for the current voxel and v_index + res_filter = range_filter * nx_filter[..., None] + res_filter =\ + res_filter * np.reshape(uv_filter, + (1, 1, 1, len(uv_filter))) + out_sf[ii, jj, kk] = np.sum( + sf_v[ii:ii+h_w, jj:jj+h_h, kk:kk+h_d] * res_filter) + out_sf[ii, jj, kk] /= np.sum(res_filter) + + return out_sf def cosine_filtering(in_sh, sh_order=8, sh_basis='descoteaux07', @@ -422,7 +526,7 @@ def cosine_filtering(in_sh, sh_order=8, sh_basis='descoteaux07', sphere = get_sphere(sphere_str) # Normalized filter for each sf direction - weights = _get_weights(sphere, dot_sharpness, sigma) + weights = _get_cosine_weights(sphere, dot_sharpness, sigma) nb_sf = len(sphere.vertices) mean_sf = np.zeros(np.append(in_sh.shape[:-1], nb_sf)) @@ -459,7 +563,7 @@ def cosine_filtering(in_sh, sh_order=8, sh_basis='descoteaux07', return out_sh -def _get_weights(sphere, dot_sharpness, sigma): +def _get_cosine_weights(sphere, dot_sharpness, sigma): """ Get neighbors weight in respect to the direction to a voxel. @@ -503,3 +607,47 @@ def _get_weights(sphere, dot_sharpness, sigma): weights /= weights.reshape((-1, weights.shape[-1])).sum(axis=0) return weights + + +def _evaluate_gaussian_distribution(x, sigma): + """ + 1-dimensional 0-centered Gaussian distribution + with standard deviation sigma. + + Parameters + ---------- + x: ndarray or float + Points where the distribution is evaluated. + sigma: float + Standard deviation. + + Returns + ------- + out: ndarray or float + Values at x. + """ + assert sigma > 0.0, "Sigma must be greater than 0." + cnorm = 1.0 / sigma / np.sqrt(2.0*np.pi) + return cnorm * np.exp(-x**2/2/sigma**2) + + +def _get_window_directions(shape): + """ + Get directions from center voxel to all neighbours + for a window of given shape. + + Parameters + ---------- + shape: tuple + Dimensions of the window. + + Returns + ------- + grid: ndarray + Grid containing the direction from the center voxel to + the current position for all positions inside the window. + """ + grid = np.indices(shape) + grid = np.moveaxis(grid, 0, -1) + grid = grid - np.asarray(shape) // 2 + return grid diff --git a/scilpy/denoise/tests/test_asym_filtering.py b/scilpy/denoise/tests/test_asym_filtering.py index a574bffa5..b1405548f 100644 --- a/scilpy/denoise/tests/test_asym_filtering.py +++ b/scilpy/denoise/tests/test_asym_filtering.py @@ -1,32 +1,46 @@ import numpy as np from scilpy.denoise.asym_filtering import \ - angle_aware_bilateral_filtering_cpu, cosine_filtering + unified_filtering, cosine_filtering from scilpy.reconst.utils import get_sh_order_and_fullness from scilpy.tests.arrays import ( fodf_3x3_order8_descoteaux07, - fodf_3x3_order8_descoteaux07_filtered_bilateral, + fodf_3x3_order8_descoteaux07_filtered_unified, fodf_3x3_order8_descoteaux07_filtered_cosine) -def test_angle_aware_bilateral_filtering(): +def test_unified_asymmetric_filtering(): """ - Test angle_aware_bilateral_filtering_cpu on a simple 3x3 grid. + Test AsymmetricFilter on a simple 3x3 grid. """ in_sh = fodf_3x3_order8_descoteaux07 sh_basis = 'descoteaux07' + legacy = True sphere_str = 'repulsion100' sigma_spatial = 1.0 - sigma_angular = 1.0 - sigma_range = 1.0 + sigma_align = 0.8 + sigma_range = 0.2 + sigma_angle = 0.06 + win_hwidth = 3 + exclude_center = False + device_type = 'cpu' + use_opencl = False sh_order, full_basis = get_sh_order_and_fullness(in_sh.shape[-1]) - out = angle_aware_bilateral_filtering_cpu(in_sh, sh_order, - sh_basis, full_basis, True, - sphere_str, sigma_spatial, - sigma_angular, sigma_range) - - assert np.allclose(out, fodf_3x3_order8_descoteaux07_filtered_bilateral) + asym_sh = unified_filtering(in_sh, sh_order, sh_basis, + is_legacy=legacy, + full_basis=full_basis, + sphere_str=sphere_str, + sigma_spatial=sigma_spatial, + sigma_align=sigma_align, + sigma_angle=sigma_angle, + rel_sigma_range=sigma_range, + win_hwidth=win_hwidth, + exclude_center=exclude_center, + device_type=device_type, + use_opencl=use_opencl) + + assert np.allclose(asym_sh, fodf_3x3_order8_descoteaux07_filtered_unified) def test_cosine_filtering(): @@ -35,12 +49,13 @@ def test_cosine_filtering(): """ in_sh = fodf_3x3_order8_descoteaux07 sh_basis = 'descoteaux07' + legacy = True sphere_str = 'repulsion100' sigma_spatial = 1.0 sharpness = 1.0 sh_order, full_basis = get_sh_order_and_fullness(in_sh.shape[-1]) - out = cosine_filtering(in_sh, sh_order, sh_basis, full_basis, True, + out = cosine_filtering(in_sh, sh_order, sh_basis, full_basis, legacy, sharpness, sphere_str, sigma_spatial) assert np.allclose(out, fodf_3x3_order8_descoteaux07_filtered_cosine) diff --git a/scilpy/gpuparallel/opencl_utils.py b/scilpy/gpuparallel/opencl_utils.py index 4f093d7c6..c6f5dd716 100644 --- a/scilpy/gpuparallel/opencl_utils.py +++ b/scilpy/gpuparallel/opencl_utils.py @@ -9,23 +9,35 @@ cl, have_opencl, _ = optional_package('pyopencl') +def cl_device_type(device_type_str): + if device_type_str == 'cpu': + return cl.device_type.CPU + if device_type_str == 'gpu': + return cl.device_type.GPU + return -1 + + class CLManager(object): """ - Class for managing an OpenCL GPU program. + Class for managing an OpenCL program. Wraps a subset of pyopencl functions to simplify its - integration with python. + integration with python. The OpenCL program can be run + on the cpu or on the gpu, given the appropriate drivers + are installed. + + When multiple cpu or gpu are available, the + one that first comes up in the list of available devices + is selected. Parameters ---------- cl_kernel: CLKernel object The CLKernel containing the OpenCL program to manage. - n_inputs: int - Number of input buffers for the kernel. - n_outputs: int - Number of output buffers for the kernel. + device_type: string + The device onto which to run the program. One of 'cpu', 'gpu'. """ - def __init__(self, cl_kernel, n_inputs, n_outputs): + def __init__(self, cl_kernel, device_type='gpu'): if not have_opencl: raise RuntimeError('pyopencl is not installed. ' 'Cannot create CLManager instance.') @@ -34,8 +46,12 @@ def __init__(self, cl_kernel, n_inputs, n_outputs): logging.getLogger('pytools.persistent_dict').setLevel(logging.CRITICAL) logging.getLogger('pyopencl').setLevel(logging.CRITICAL) - self.input_buffers = [0] * n_inputs - self.output_buffers = [0] * n_outputs + self.input_buffers = [] + self.output_buffers = [] + + # maps key to index in buffers list + self.inputs_mapping = {} + self.outputs_mapping = {} # Find the best device for running GPU tasks platforms = cl.get_platforms() @@ -43,8 +59,13 @@ def __init__(self, cl_kernel, n_inputs, n_outputs): for p in platforms: devices = p.get_devices() for d in devices: - if best_device is None: - best_device = d + d_type = d.get_info(cl.device_info.TYPE) + if d_type == cl_device_type(device_type)\ + and best_device is None: + best_device = d # take the first device of right type + + if best_device is None: + raise ValueError('No device of type {} found'.format(device_type)) self.context = cl.Context(devices=[best_device]) self.queue = cl.CommandQueue(self.context) @@ -69,7 +90,7 @@ def __init__(self, buf, shape, dtype): self.shape = shape self.dtype = dtype - def add_input_buffer(self, arg_pos, arr, dtype=np.float32): + def add_input_buffer(self, key, arr=None, dtype=np.float32): """ Add an input buffer to the kernel program. Input buffers must be added in the same order as they are declared inside @@ -77,19 +98,15 @@ def add_input_buffer(self, arg_pos, arr, dtype=np.float32): Parameters ---------- - arg_pos: int - Position of the buffer in the input buffers list. + key: string + Name of the buffer in the input buffers list. Used for + referencing when updating buffers. arr: numpy ndarray Data array. dtype: dtype, optional Optional type for array data. It is recommended to use float32 whenever possible to avoid unexpected behaviours. - Returns - ------- - indice: int - Index of the input buffer in the input buffers list. - Note ---- Array is reordered as fortran array and then flattened. This is @@ -98,30 +115,95 @@ def add_input_buffer(self, arg_pos, arr, dtype=np.float32): For example, for a 3-dimensional array of shape (X, Y, Z), the flat index for position i, j, k is idx = i + j * X + z * X * Y. """ - # convert to fortran ordered, dtype array + buf = None + if arr is not None: + # convert to fortran ordered, dtype array + arr = np.asfortranarray(arr, dtype=dtype) + buf = cl.Buffer(self.context, cl.mem_flags.READ_ONLY | + cl.mem_flags.COPY_HOST_PTR, hostbuf=arr) + + if key in self.inputs_mapping.keys(): + raise ValueError('Invalid key for buffer!') + + self.inputs_mapping[key] = len(self.input_buffers) + self.input_buffers.append(buf) + + def update_input_buffer(self, key, arr, dtype=np.float32): + """ + Update an input buffer. Input buffers must first be added + to program using `add_input_buffer`. + + Parameters + ---------- + key: string + Name of the buffer in the input buffers list. + arr: numpy ndarray + Data array. + dtype: dtype, optional + Optional type for array data. It is recommended to use float32 + whenever possible to avoid unexpected behaviours. + """ + if key not in self.inputs_mapping.keys(): + raise ValueError('Invalid key for buffer!') + argpos = self.inputs_mapping[key] + arr = np.asfortranarray(arr, dtype=dtype) buf = cl.Buffer(self.context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=arr) - self.input_buffers[arg_pos] = buf + self.input_buffers[argpos] = buf - def add_output_buffer(self, arg_pos, shape, dtype=np.float32): + def add_output_buffer(self, key, shape=None, dtype=np.float32): """ - Add an output buffer. + Add an output buffer to the kernel program. Output buffers + must be added in the same order as they are declared inside + the kernel code (.cl file). Parameters ---------- - arg_pos: int - Position of the buffer in the output buffers list. + key: string + Name of the buffer in the output buffers list. Used for + referencing when updating buffers. shape: tuple Shape of the output array. dtype: dtype, optional - Data type for the output. It is recommended to keep - float32 to avoid unexpected behaviour. + Optional type for array data. It is recommended to use float32 + whenever possible to avoid unexpected behaviours. + """ + if key in self.outputs_mapping.keys(): + raise ValueError('Invalid key for buffer!') + + buf = None + if shape is not None: + buf = cl.Buffer(self.context, cl.mem_flags.WRITE_ONLY, + np.prod(shape) * np.dtype(dtype).itemsize) + + self.outputs_mapping[key] = len(self.output_buffers) + self.output_buffers.append(self.OutBuffer(buf, shape, dtype)) + + def update_output_buffer(self, key, shape, dtype=np.float32): """ + Update an output buffer. Output buffers must first be added + to program using `add_output_buffer`. + + Parameters + ---------- + key: string + Name of the buffer in the output buffers list. + shape: tuple + New shape of the output array. + dtype: dtype, optional + Optional type for array data. It is recommended to use float32 + whenever possible to avoid unexpected behaviours. + """ + if key not in self.outputs_mapping.keys(): + raise ValueError('Invalid key for buffer!') + argpos = self.outputs_mapping[key] + buf = cl.Buffer(self.context, cl.mem_flags.WRITE_ONLY, np.prod(shape) * np.dtype(dtype).itemsize) - self.output_buffers[arg_pos] = self.OutBuffer(buf, shape, dtype) + out_buf = self.OutBuffer(buf, shape, dtype) + self.output_buffers[argpos] = out_buf def run(self, global_size, local_size=None): """ @@ -174,7 +256,11 @@ class CLKernel(object): """ def __init__(self, entrypoint, module, filename): path_to_kernel = self._get_kernel_path(module, filename) - f = open(path_to_kernel, 'r') + try: + f = open(path_to_kernel, 'r') + except Exception: + raise ValueError('OpenCL file not found in {}' + .format(path_to_kernel)) self.code = f.readlines() self.entrypoint = entrypoint diff --git a/scilpy/tests/arrays.py b/scilpy/tests/arrays.py index 45021ce41..e25087d38 100644 --- a/scilpy/tests/arrays.py +++ b/scilpy/tests/arrays.py @@ -6,420 +6,420 @@ # in the FiberCup. Order 8, descoteaux07 basis. fodf_3x3_order8_descoteaux07 = np.asarray([[[ [0.26286387, -0.336828, -0.00369634, -0.25474626, - -0.01061569, 0.29750797, 0.03732642, 0.01171304, - 0.16143565, 0.00787683, 0.14300859, 0.01041067, - -0.13926704, -0.00285504, -0.28675956, 0.06106967, + -0.01061569, 0.29750797, 0.03732642, 0.01171304, + 0.16143565, 0.00787683, 0.14300859, 0.01041067, + -0.13926704, -0.00285504, -0.28675956, 0.06106967, -0.00251915, -0.0103068, -0.00626346, -0.05204236, - -0.00324461, -0.0468456, -0.00444966, 0.04529396, - 0.00119645, 0.07706978, 0.00584123, 0.08920382]], + -0.00324461, -0.0468456, -0.00444966, 0.04529396, + 0.00119645, 0.07706978, 0.00584123, 0.08920382]], [[0.27447942, -0.38676834, -0.02799141, -0.23060161, - -0.03448979, 0.09528257, 0.25290385, -0.01982788, - 0.22237155, 0.00135294, 0.13581508, 0.01650201, - -0.0275899, 0.00186949, -0.09112225, -0.08496156, - 0.00396445, -0.07162198, 0.00324492, -0.07831519, - 0.00616271, -0.05209329, -0.01479937, 0.01713125, - -0.02268516, 0.04751123, -0.03049234, 0.06638705]], - - [[0.27134237, -0.4307746, 0.0035063, -0.2592252, - -0.00330402, -0.15647785, 0.21888971, 0.00334128, - 0.22812255, 0.01228249, 0.16825177, 0.00471007, - 0.08927008, 0.00575425, 0.19901836, -0.04885881, + -0.03448979, 0.09528257, 0.25290385, -0.01982788, + 0.22237155, 0.00135294, 0.13581508, 0.01650201, + -0.0275899, 0.00186949, -0.09112225, -0.08496156, + 0.00396445, -0.07162198, 0.00324492, -0.07831519, + 0.00616271, -0.05209329, -0.01479937, 0.01713125, + -0.02268516, 0.04751123, -0.03049234, 0.06638705]], + + [[0.27134237, -0.4307746, 0.0035063, -0.2592252, + -0.00330402, -0.15647785, 0.21888971, 0.00334128, + 0.22812255, 0.01228249, 0.16825177, 0.00471007, + 0.08927008, 0.00575425, 0.19901836, -0.04885881, -0.00579403, -0.07190462, -0.00737415, -0.08707161, -0.00195517, -0.06621928, -0.00356173, -0.03289294, -0.00356576, -0.06258938, -0.00665118, -0.09943768]]], - [[[0.15886374, 0.01589506, 0.06487983, -0.02096088, - 0.02011732, 0.00589444, 0.03856898, 0.00757658, - 0.05849233, -0.11192703, 0.00713574, -0.00371081, - -0.0244258, 0.0117896, 0.03056652, -0.01305971, + [[[0.15886374, 0.01589506, 0.06487983, -0.02096088, + 0.02011732, 0.00589444, 0.03856898, 0.00757658, + 0.05849233, -0.11192703, 0.00713574, -0.00371081, + -0.0244258, 0.0117896, 0.03056652, -0.01305971, 0.00679529, -0.01269033, -0.00467887, -0.01755827, - -0.00761302, 0.02157611, 0.02873131, -0.01170345, - 0.00136947, -0.00745852, 0.02145799, 0.00045532]], + -0.00761302, 0.02157611, 0.02873131, -0.01170345, + 0.00136947, -0.00745852, 0.02145799, 0.00045532]], [[0.17820162, -0.08854709, -0.0138819, -0.08386812, - 0.04722548, -0.03020687, 0.03683474, 0.0020233, - 0.02984233, -0.07024448, 0.04935116, 0.05288185, + 0.04722548, -0.03020687, 0.03683474, 0.0020233, + 0.02984233, -0.07024448, 0.04935116, 0.05288185, -0.06069153, -0.02227145, -0.08536565, -0.00739656, - -0.00208759, -0.03538536, 0.04175738, -0.01890082, + -0.00208759, -0.03538536, 0.04175738, -0.01890082, 0.03184003, -0.04322833, -0.0197769, -0.01709824, - 0.03408453, 0.03452914, -0.00768105, 0.02879811]], + 0.03408453, 0.03452914, -0.00768105, 0.02879811]], - [[0.17783499, -0.08317029, 0.02127749, -0.05431161, - -0.0426516, 0.00668156, -0.10225131, 0.05116061, - 0.05691102, 0.06157963, 0.08037782, -0.0264888, - 0.0032487, 0.00127227, -0.00224557, 0.08135595, - -0.00129513, 0.04060476, -0.00118032, -0.00987782, + [[0.17783499, -0.08317029, 0.02127749, -0.05431161, + -0.0426516, 0.00668156, -0.10225131, 0.05116061, + 0.05691102, 0.06157963, 0.08037782, -0.0264888, + 0.0032487, 0.00127227, -0.00224557, 0.08135595, + -0.00129513, 0.04060476, -0.00118032, -0.00987782, 0.00490772, -0.03195551, -0.01908543, -0.03024822, - -0.00608355, -0.01933366, 0.00161467, -0.00429689]]], - - [[[0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0.]], - - [[0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0.]], - - [[0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0., - 0., 0., 0., 0.]]]]) - -fodf_3x3_order8_descoteaux07_filtered_bilateral = np.asarray([[[[ - 5.99904145e-02, -7.74841521e-03, -8.85743552e-04, - -2.43733904e-02, -6.71300622e-02, 1.40698629e-03, - -5.02532006e-02, -1.31764043e-03, 4.39820686e-02, - 3.04767348e-03, -2.60659751e-03, 1.19609175e-03, - 2.45643515e-03, 1.60415663e-02, -3.71781262e-04, - 1.72503732e-02, 1.87847737e-02, 1.81510421e-03, - 3.45741639e-02, -4.08627879e-03, 2.70963242e-02, - 2.48334799e-03, -2.08282714e-02, -5.04963363e-05, - -4.04026968e-02, -4.17118816e-03, 9.70829612e-05, - -2.50408771e-03, 2.67723739e-03, 1.28968880e-04, - -1.25406732e-03, -7.78307195e-03, 4.65415094e-04, - -6.77331216e-03, 5.83355784e-06, -7.46204986e-03, - 2.79709665e-03, 2.65259651e-04, -4.86151089e-03, - 9.07678667e-05, -9.86897362e-03, 8.31694950e-04, - -7.30557944e-03, -1.27246592e-04, 4.47773980e-03, - -1.67805475e-04, 9.45333969e-03, -2.06194837e-05, - 1.19741653e-02]], - - [[7.30659454e-02, 6.26644064e-04, -6.45042986e-05, - -2.44371069e-03, -8.71619020e-02, -2.97978698e-03, - -5.70290239e-02, -4.88288484e-03, 1.81396370e-02, - -1.49746749e-02, -1.82596238e-03, -3.73852432e-03, - 1.91679528e-03, 3.46388382e-03, 9.45765979e-04, - 6.35821093e-03, 4.43337556e-02, -1.14045984e-03, - 4.79787393e-02, -2.69629492e-03, 3.30809021e-02, - 3.90751032e-03, -7.69518614e-03, -1.00509094e-03, - -1.83497970e-02, 1.09524242e-02, 3.03376250e-05, - 5.84822874e-03, 1.28574628e-03, 1.31931818e-03, - -1.62226350e-03, -2.42913519e-03, -2.51710886e-03, - -1.46265067e-03, 9.64475618e-05, -1.55335427e-03, - -8.35208072e-03, 1.14634507e-04, -1.15803123e-02, - 9.78385907e-04, -1.45041153e-02, 2.11820605e-03, - -1.11267735e-02, -2.54739831e-03, 1.13908925e-03, - -1.58677391e-03, 6.85953527e-03, -4.07763058e-03, - 8.15545468e-03]], - - [[6.06366544e-02, -2.86155218e-03, -6.31169176e-04, - 2.39112903e-02, -8.10810515e-02, 6.08200360e-04, - -5.06141009e-02, -2.06508995e-03, -1.61198389e-02, - -1.83721383e-03, -7.99761006e-04, -4.98098056e-04, - 6.80259319e-04, -1.43673956e-02, 5.60409075e-04, - -1.55892906e-02, 3.52407019e-02, 1.52343110e-03, - 4.21572957e-02, 1.87781266e-03, 3.20393070e-02, - 1.71999179e-03, 8.40048303e-03, 1.82572535e-04, - 1.77166894e-02, 5.35407360e-03, -1.26213948e-03, - 7.51220125e-04, -9.24301077e-04, -9.34498808e-04, - -8.12227912e-04, 7.45251831e-03, -1.46771071e-03, - 7.26425049e-03, -2.49284089e-04, 7.85839012e-03, - -5.55945668e-03, -3.07526270e-04, -9.91274872e-03, - 2.33305790e-04, -1.34808479e-02, 7.95748091e-04, - -1.14669773e-02, -1.89505175e-03, -4.30439514e-03, - -2.42432697e-04, -4.67979922e-03, -1.23980121e-03, - -6.26678862e-03]]], - - - [[[4.68592306e-02, 3.61461743e-03, -7.63514670e-04, - -1.20844533e-02, -2.12680062e-02, 9.29568778e-03, - -2.01531071e-02, 3.47516906e-03, 9.20440534e-03, - -5.46795536e-03, 1.59709356e-04, -1.53637043e-03, - -2.87913528e-04, 6.47466282e-03, -4.50224823e-04, - 1.07221387e-02, 1.27142040e-02, 1.29577853e-03, - 2.05155038e-02, -1.90391848e-02, 1.09376471e-02, - 1.49402131e-03, -9.46249960e-03, 9.19890098e-04, - -5.18554206e-03, 4.58351245e-05, -1.76770467e-04, - 1.37229706e-03, 5.83527846e-04, 9.67807760e-04, - 1.16870730e-03, -3.93367403e-03, 1.84726513e-03, - -4.60773077e-03, -1.23348402e-04, -7.24725612e-03, - -1.76239341e-03, 1.04030103e-03, -4.23162632e-03, - 8.57213815e-04, -6.38065732e-03, 1.90766279e-05, - -5.94726402e-04, 3.32327855e-03, -1.41881029e-03, - 8.57485896e-04, 1.63106010e-03, 2.71103567e-03, - 2.36989013e-03]], - - [[5.90628671e-02, 7.12337359e-03, 6.30805066e-04, - -5.89538262e-04, -4.37513131e-02, -7.76500224e-04, - -3.43137507e-02, 5.36402775e-03, 1.64261405e-03, - -1.49954433e-02, -1.53570457e-03, -4.66620624e-03, - 1.74906011e-04, 1.35081847e-03, 1.88758144e-03, - 1.86557039e-03, 1.29785439e-02, 1.47604347e-03, - 2.23300761e-02, -1.19166709e-02, 2.09856438e-02, - 7.91100937e-03, -1.17660610e-02, -3.37538181e-03, - -1.61214712e-02, 9.81557298e-03, 9.29820184e-05, - 7.44981371e-03, 3.81985318e-04, 2.48233403e-03, - 4.68702466e-04, -2.03452776e-03, -2.59054263e-03, - 8.00222286e-04, -1.51345425e-03, 1.81863005e-03, - 2.16841328e-03, -4.22872432e-04, -6.79604466e-03, - 6.29856690e-03, -8.17551116e-03, 5.05505647e-03, - -1.05657360e-02, -2.98331814e-03, -3.59175521e-03, - 4.73267756e-03, 5.64775497e-03, -1.33188152e-03, - 5.08498298e-03]], - - [[4.97741267e-02, 3.08713766e-03, 5.63873392e-04, - 1.55436843e-02, -3.80921784e-02, 2.91715226e-03, - -2.50003788e-02, -5.73311922e-03, 4.61819744e-04, - -1.07811659e-02, -9.32497997e-04, -3.14809198e-03, - 9.03481742e-04, -7.94237810e-03, 1.94605341e-04, - -8.87043651e-03, -4.64460095e-03, 7.54078202e-03, - 2.16862677e-02, 7.15296228e-03, 2.22150527e-02, - -2.14696419e-03, -1.08946943e-03, -5.87122761e-04, - -2.62994270e-03, 9.15782450e-03, 1.61327436e-04, - 5.93658958e-03, -3.08842666e-04, 1.49215900e-03, - -1.28381218e-03, 4.03265153e-03, -2.26367799e-03, - 4.64427521e-03, -4.83266613e-04, 5.48291770e-03, - 1.11924666e-02, -2.16552214e-04, 2.52305299e-03, - 1.40712050e-03, -5.92827841e-03, 1.91162035e-03, - -8.85321885e-03, -3.83274815e-03, -5.31446078e-03, - -3.26814596e-04, -1.33576089e-03, -4.25520553e-04, - 1.38644799e-03]]], - - [[[ - 1.33088596e-02, 7.70245330e-03, 9.03457852e-04, - -4.11949105e-03, -4.39007513e-03, 2.52064292e-03, - -6.03218056e-03, 1.13561199e-03, 9.42679670e-04, - -9.34460087e-04, 2.41541137e-03, -5.82004266e-04, - -3.41451676e-03, 1.14922086e-03, 4.54078257e-04, - 2.27080762e-03, 2.60575889e-03, 9.58289863e-04, - 4.55610130e-03, -5.43882630e-03, 2.91891493e-03, - 1.00431022e-03, -2.53372884e-03, -3.86116404e-05, - -1.87465687e-03, 5.32383426e-04, 1.50538019e-04, - 9.05098395e-04, -1.81109635e-03, -1.69890726e-04, - 2.60153827e-03, -1.12953050e-03, 1.06985723e-03, - -1.62234962e-03, 2.15486359e-04, -1.12075310e-03, - -1.54880131e-04, 1.26322377e-04, -1.06276650e-03, - 4.80723343e-04, -1.65165470e-03, 1.03662292e-03, - -6.90509683e-04, 1.34792802e-04, -6.62758017e-04, - 6.82414480e-04, 5.65946760e-04, 5.43906885e-04, - 8.88275606e-04]], - - [[1.69216413e-02, 8.84666010e-03, 8.79326427e-04, - -1.15188140e-03, -8.48558842e-03, 1.44055512e-03, - -8.69017052e-03, 9.04357052e-04, -3.34932934e-04, - -5.36503671e-03, 1.06440860e-03, -2.62392070e-03, - -2.00387650e-03, -2.19755205e-04, 1.05876307e-03, - -1.70799645e-03, 1.29642088e-03, 1.17827494e-03, - 5.66355890e-03, -3.51558336e-03, 5.59002356e-03, - 8.86688817e-04, -2.01804999e-03, -3.48552502e-04, - -2.69504812e-03, 2.38875030e-03, 3.96558766e-04, - 2.72303048e-03, -1.03172620e-03, 1.79526273e-03, - 2.04984828e-03, -7.73159626e-04, 1.58370808e-04, - 2.60501679e-04, -9.53109540e-04, 5.62474486e-04, - 1.49585605e-03, -4.40141423e-05, -7.77367586e-04, - 9.91113015e-04, -1.78421857e-03, 1.38476900e-03, - -2.39602862e-03, -4.47369469e-04, -1.53124355e-03, - 9.33497953e-04, 7.39396498e-04, -2.10889604e-04, - 1.11924543e-03]], - - [[1.32409122e-02, 6.01101256e-03, 7.43026309e-04, - 4.56297819e-03, -8.77641953e-03, 7.22086091e-04, - -6.91051977e-03, -4.05665172e-04, 3.24728365e-04, - -7.13127007e-03, 6.72168797e-04, -2.20440201e-03, - 4.07174840e-04, -2.47174339e-03, -2.38993843e-04, - -2.79616112e-03, -9.28663607e-04, 1.61242764e-03, - 5.14579103e-03, 3.25095461e-04, 5.38503992e-03, - 8.93433421e-05, -1.16181579e-03, -3.91830597e-04, - -1.78770414e-03, 2.16618425e-03, 1.02836085e-03, - 3.74908286e-03, 1.19529788e-04, 1.99549249e-03, - -4.90335513e-04, 8.27447594e-04, -7.07318268e-04, - 1.03643921e-03, -3.29578680e-04, 6.04670460e-04, - 2.37059976e-03, 2.06921931e-04, 4.51941221e-04, - 7.50580440e-04, -1.43936576e-03, 5.88831780e-04, - -2.45795081e-03, -7.40358468e-04, -1.21437997e-03, - 2.45072224e-04, 2.64311350e-04, -1.81628898e-04, - 6.73057957e-04]]]]) + -0.00608355, -0.01933366, 0.00161467, -0.00429689]]], + + [[[0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0.]], + + [[0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0.]], + + [[0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0., + 0., 0., 0., 0.]]]]) + +fodf_3x3_order8_descoteaux07_filtered_unified = np.array([ + [[[2.29357947e-01, -3.40501584e-03, -6.99261017e-04, + -7.19531128e-03, -2.82538626e-01, 1.39868002e-03, + -2.13735138e-01, -6.60276345e-03, 2.40092790e-01, + -6.49857556e-03, -1.47142046e-03, -1.74323477e-03, + 1.97528288e-03, 2.74056237e-03, 1.40669684e-03, + 8.45362867e-03, 4.87161977e-02, 3.29667207e-03, + 1.56340999e-01, -4.53684466e-03, 1.36801062e-01, + 5.76189638e-03, -1.35952240e-01, 7.04539644e-03, + -2.80666044e-01, 2.01464026e-03, -3.04351606e-03, + 2.55501754e-03, 1.77060275e-03, 1.66821237e-03, + -4.83733413e-04, 1.34162234e-03, -1.88498937e-03, + -8.92315295e-04, -1.85175590e-03, -9.90733681e-03, + 7.96610117e-02, -6.93678171e-03, -1.32213544e-02, + 9.29201219e-04, -6.76385227e-02, 3.24539565e-03, + -6.21090167e-02, 1.48437903e-04, 6.24298982e-02, + -6.98251925e-03, 1.09354237e-01, -7.46142122e-03, + 1.39128340e-01]], + + [[2.23049655e-01, -6.34489485e-04, -2.46164995e-04, + 1.11559735e-03, -3.33977838e-01, -5.91911022e-03, + -2.00774725e-01, -2.42845691e-02, 7.05612574e-02, + -1.39213801e-02, -1.23520352e-03, -2.34371080e-03, + 1.91758690e-03, 8.65363592e-04, 1.02889761e-03, + 4.72826644e-03, 2.33056470e-01, -1.02386226e-03, + 1.93300973e-01, -9.51278280e-04, 1.29224787e-01, + 3.09759106e-02, -3.46205350e-02, 2.24824994e-02, + -9.29648766e-02, 1.38485242e-02, -4.48354553e-04, + 4.43312110e-03, 1.23785504e-03, 1.47203519e-04, + -1.48610994e-03, -1.48539730e-03, -1.71649289e-03, + -1.24889989e-03, 8.08452709e-04, -5.42925265e-03, + -1.06615393e-01, 5.72984074e-04, -8.65200925e-02, + 5.99094215e-04, -8.35112741e-02, 9.02314727e-04, + -5.75274392e-02, -3.27219474e-02, 1.61903376e-02, + -3.03618703e-02, 3.50629229e-02, -3.12781970e-02, + 6.09901309e-02]], + + [[2.39271129e-01, -1.29514017e-03, 2.67539048e-05, + 7.77199752e-03, -3.65455973e-01, 1.59433374e-03, + -2.30378736e-01, -1.71148930e-03, -1.31054065e-01, + -6.63373599e-03, -2.62685328e-03, -1.54850825e-03, + -1.04405158e-03, -8.82561059e-04, -8.40336616e-04, + -4.02329843e-03, 2.14542732e-01, 9.81212746e-04, + 2.25362349e-01, 1.70438428e-03, 1.68648111e-01, + 1.72316326e-03, 8.47390835e-02, -1.94316083e-03, + 1.81578935e-01, 1.24860810e-02, 8.35937367e-04, + 3.26296874e-03, 1.87537964e-03, 1.16841921e-04, + 1.56804161e-03, -3.77935733e-03, 5.18925927e-04, + -1.17750904e-03, 2.33482491e-03, 5.33718948e-03, + -6.84234291e-02, -7.23686798e-03, -9.02332211e-02, + -4.07916322e-03, -1.21786755e-01, -4.40776431e-04, + -9.72842248e-02, -2.05488288e-04, -4.92693194e-02, + 1.53163774e-03, -8.40887907e-02, 1.18279868e-03, + -1.15965859e-01]]], + + + [[[4.90583299e-02, 1.42501402e-03, -8.36557553e-04, + -1.15831440e-02, -4.22656583e-03, 1.77363767e-02, + -1.31061492e-02, 6.58267947e-03, 3.57744477e-03, + -1.19818257e-03, 7.13707057e-04, 1.00507324e-04, + -4.05178649e-05, 5.72719299e-03, -8.62838694e-04, + 6.93407399e-03, 1.28787065e-02, 3.71612264e-03, + 1.84145755e-02, -3.03712343e-02, 4.27460381e-03, + -1.05287718e-03, -6.81217935e-03, 2.91621036e-03, + 7.67712088e-03, -1.67774735e-03, -9.17083386e-04, + -7.05515351e-04, -2.57881938e-04, -5.39230247e-05, + 5.44975228e-04, -2.91389651e-03, 2.49849340e-03, + -2.18778612e-03, 4.26276808e-04, -1.88276038e-03, + -3.92622090e-03, 3.22745694e-03, -2.49579185e-03, + -5.16371150e-04, -7.14925457e-03, -3.53600795e-04, + 6.22853770e-03, 6.75767637e-03, -4.95233255e-03, + 1.00070512e-03, -5.30703269e-04, 7.75615036e-03, + 2.13095937e-04]], + + [[6.83494656e-02, 2.47590868e-03, 8.04082443e-04, + -1.07085124e-03, -4.66462691e-02, -3.14996372e-03, + -3.67278245e-02, 1.24317997e-02, 2.84904197e-04, + -3.67900721e-03, -2.45262880e-03, -1.25632506e-03, + -3.29065414e-04, 1.34005377e-03, 2.27011899e-03, + 1.12787213e-04, 1.57892279e-02, -6.42922533e-04, + 2.23863221e-02, -1.98960014e-02, 2.30945595e-02, + 1.22861979e-02, -2.09248819e-02, -7.79075565e-03, + -3.12206640e-02, 6.21542463e-05, 2.42862058e-04, + 1.12858477e-03, 1.19409512e-03, 2.01720154e-04, + 9.41874244e-04, -1.50090832e-03, -2.70566377e-03, + 2.13984489e-03, -2.74112863e-03, 4.67208319e-03, + -2.84238532e-03, 1.19886439e-03, -1.42714247e-02, + 1.35054050e-02, -1.26888266e-02, 8.93556030e-03, + -1.69535796e-02, -3.44881983e-03, -3.29416809e-03, + 1.15606070e-02, 1.51758509e-02, -2.17242328e-03, + 1.61578521e-02]], + + [[6.20684558e-02, 2.29537210e-03, 1.57782979e-04, + 1.21746858e-02, -3.75416568e-02, 3.57151832e-03, + -2.99229041e-02, -1.04691005e-02, 7.93762774e-04, + -8.24432299e-03, -8.34753632e-04, -2.42344377e-03, + 1.20899492e-03, -5.68896102e-03, 4.33085595e-04, + -4.85994374e-03, -3.41545076e-02, 1.46092528e-02, + 2.35147496e-02, 1.53551134e-02, 3.09649758e-02, + -3.32700724e-03, 7.38885638e-04, -5.33398930e-04, + -2.20765157e-03, 4.73109351e-03, 8.27559708e-04, + 4.68404939e-03, -2.84124885e-04, 1.27799454e-03, + -1.45038536e-03, 2.38853146e-03, -1.98245860e-03, + 2.58454776e-03, -6.55298954e-04, 1.39578469e-03, + 3.98648608e-02, -1.91019079e-03, 1.58616283e-02, + -2.31384838e-03, -6.71152510e-03, 8.87343339e-04, + -1.49804902e-02, -7.65802177e-03, -1.06220121e-02, + -1.99364911e-03, -4.93716710e-03, 1.93583996e-03, + 2.37014751e-03]]], + + + [[[8.56465561e-03, 7.14008949e-03, 9.31430807e-04, + -3.60053409e-03, 3.43419984e-04, 2.32896514e-03, + -3.29158998e-03, 6.33322848e-04, -1.41612004e-03, + -1.96906728e-05, 2.08462933e-03, -4.70169796e-04, + -3.18806182e-03, 7.96508703e-04, 1.69068438e-04, + 8.53310802e-04, 7.43969959e-04, 6.92280373e-04, + 1.56785435e-03, -4.36363329e-03, 8.41585833e-04, + 7.77737118e-04, -8.43136432e-04, -2.88294327e-04, + 4.52821928e-04, 2.30036728e-04, -1.17694055e-04, + 3.29426513e-04, -1.69867522e-03, -1.36557340e-04, + 2.43068996e-03, -5.02415966e-04, 1.08299347e-03, + -9.28927692e-04, -4.96423591e-05, 9.66391487e-05, + -7.75296465e-05, -7.81383154e-05, -3.84702522e-04, + 1.38600809e-04, -1.24402688e-04, 1.24713482e-03, + -4.63949171e-05, -1.84301807e-04, -5.34927571e-04, + 6.78179065e-04, -4.90089799e-04, 3.53365269e-04, + -2.82389612e-04]], + + [[1.00402920e-02, 7.58354607e-03, 8.90893563e-04, + -1.28660192e-03, -1.33702532e-03, 1.36204823e-03, + -4.63267660e-03, 6.65906190e-04, -1.62515435e-03, + -2.18256882e-03, 8.05148913e-04, -1.99928992e-03, + -1.83230347e-03, -1.09679838e-05, 9.06606032e-04, + -1.63774703e-03, -7.68655802e-04, 5.10504917e-04, + 1.39995532e-03, -2.66336265e-03, 2.41858563e-03, + 8.60507816e-05, -3.95330794e-04, -2.54730514e-04, + -7.51359487e-04, 9.56600867e-05, 2.33861137e-04, + 9.07833792e-04, -9.61201431e-04, 1.28055532e-03, + 1.78347418e-03, -4.80643827e-04, 2.17566393e-05, + 3.30969135e-04, -9.77124047e-04, 4.48248320e-04, + 7.68020370e-04, -9.77999769e-05, -4.11722548e-05, + 4.91272506e-04, 3.40984634e-04, 1.23504466e-03, + -8.25409359e-04, -9.52547943e-05, -1.01083781e-03, + 6.74211675e-04, 1.26126637e-04, -4.93598622e-04, + -4.93061797e-05]], + + [[7.32812182e-03, 4.68891712e-03, 7.34380826e-04, + 3.29422628e-03, -2.95965884e-03, 7.20124038e-04, + -3.21578559e-03, -1.79025321e-04, 1.01238881e-03, + -4.08613875e-03, 5.58987071e-04, -1.23820537e-03, + 2.40328249e-04, -1.56967169e-03, -1.39154716e-04, + -1.34541481e-03, -2.15931028e-03, 9.24919733e-04, + 1.77585704e-03, 3.15205162e-04, 2.14559276e-03, + -2.47588632e-04, -9.18796674e-04, -3.64654434e-04, + -1.39549298e-03, -1.32647212e-04, 8.84594449e-04, + 1.97716001e-03, 2.69772753e-04, 1.19616199e-03, + -3.84616915e-04, 3.64541696e-04, -7.01064946e-04, + 2.94422528e-04, -4.42401101e-04, -4.55646206e-04, + 1.19246383e-03, 4.13058489e-04, 1.06265511e-03, + 6.72362107e-04, 2.49843251e-04, 2.46282871e-04, + -6.64883205e-04, -1.99434148e-04, -3.93359154e-04, + 1.76576846e-05, 3.28426158e-04, -3.16710533e-04, + 5.22099584e-05]]]]) fodf_3x3_order8_descoteaux07_filtered_cosine = np.array([[ [[1.10403600e-01, -1.56274168e-02, -1.76220261e-03, -4.06611449e-02, - -1.30363567e-01, 3.75744696e-03, -9.97133892e-02, -3.49429274e-03, - 9.52583150e-02, 8.20238612e-03, -4.89889718e-03, 2.97695306e-03, - 4.70504991e-03, 2.74152194e-02, -1.06442787e-03, 3.09297306e-02, - 3.78406501e-02, 4.40385732e-03, 6.97472214e-02, -5.22356452e-03, - 5.54510683e-02, 4.19701656e-03, -4.61628754e-02, 1.04670478e-03, - -9.23817078e-02, -1.08910148e-02, 5.21920978e-04, -6.08726863e-03, + -1.30363567e-01, 3.75744696e-03, -9.97133892e-02, -3.49429274e-03, + 9.52583150e-02, 8.20238612e-03, -4.89889718e-03, 2.97695306e-03, + 4.70504991e-03, 2.74152194e-02, -1.06442787e-03, 3.09297306e-02, + 3.78406501e-02, 4.40385732e-03, 6.97472214e-02, -5.22356452e-03, + 5.54510683e-02, 4.19701656e-03, -4.61628754e-02, 1.04670478e-03, + -9.23817078e-02, -1.08910148e-02, 5.21920978e-04, -6.08726863e-03, 5.44869381e-03, -2.53482527e-04, -2.27573225e-03, -1.40288007e-02, - 1.35667761e-03, -1.30331232e-02, 5.35207622e-04, -1.63039610e-02, - 7.97932344e-03, 8.02639457e-04, -1.21136261e-02, -1.83430175e-03, - -2.54608366e-02, 9.91489623e-04, -1.97747572e-02, -1.75031223e-03, - 1.49128199e-02, -1.43106403e-03, 2.66094509e-02, 8.57021815e-04, + 1.35667761e-03, -1.30331232e-02, 5.35207622e-04, -1.63039610e-02, + 7.97932344e-03, 8.02639457e-04, -1.21136261e-02, -1.83430175e-03, + -2.54608366e-02, 9.91489623e-04, -1.97747572e-02, -1.75031223e-03, + 1.49128199e-02, -1.43106403e-03, 2.66094509e-02, 8.57021815e-04, 3.30743534e-02]], - [[1.40072235e-01, 3.36392275e-03, -1.12255448e-04, -5.78225308e-03, + [[1.40072235e-01, 3.36392275e-03, -1.12255448e-04, -5.78225308e-03, -1.77192188e-01, -6.19144118e-03, -1.17119585e-01, -1.16353192e-02, 3.66205780e-02, -3.46907611e-02, -3.57179047e-03, -9.03952453e-03, - 3.77892854e-03, 7.63422898e-03, 1.62939211e-03, 1.45448683e-02, - 9.98680216e-02, -2.60934868e-03, 1.02401140e-01, -3.99776516e-03, - 6.98372955e-02, 8.57505830e-03, -1.55541774e-02, -6.67052505e-04, - -3.87320367e-02, 2.70807684e-02, 4.61517902e-04, 1.48157281e-02, - 2.48281726e-03, 3.77428961e-03, -3.20783604e-03, -5.32528133e-03, - -4.66592653e-03, -4.11151389e-03, 1.53491240e-04, -6.28303818e-03, - -2.60040931e-02, 9.76421739e-04, -3.20749919e-02, 6.17639701e-04, - -3.80348770e-02, 3.50055436e-03, -2.89460088e-02, -6.21961773e-03, - 6.33647814e-03, -3.64798932e-03, 1.61175330e-02, -8.68729069e-03, + 3.77892854e-03, 7.63422898e-03, 1.62939211e-03, 1.45448683e-02, + 9.98680216e-02, -2.60934868e-03, 1.02401140e-01, -3.99776516e-03, + 6.98372955e-02, 8.57505830e-03, -1.55541774e-02, -6.67052505e-04, + -3.87320367e-02, 2.70807684e-02, 4.61517902e-04, 1.48157281e-02, + 2.48281726e-03, 3.77428961e-03, -3.20783604e-03, -5.32528133e-03, + -4.66592653e-03, -4.11151389e-03, 1.53491240e-04, -6.28303818e-03, + -2.60040931e-02, 9.76421739e-04, -3.20749919e-02, 6.17639701e-04, + -3.80348770e-02, 3.50055436e-03, -2.89460088e-02, -6.21961773e-03, + 6.33647814e-03, -3.64798932e-03, 1.61175330e-02, -8.68729069e-03, 1.87904410e-02]], - [[1.14020750e-01, -6.87045764e-03, -1.14232450e-03, 4.11661189e-02, - -1.63423278e-01, 1.42638757e-03, -1.02391114e-01, -3.10274637e-03, - -4.14102976e-02, -7.65200269e-05, -1.22479669e-03, 1.44271152e-04, - 1.02909411e-03, -2.53524591e-02, 8.46680097e-04, -2.91525407e-02, - 8.15595775e-02, 2.62140564e-03, 8.88886643e-02, 3.47353540e-03, - 6.60543440e-02, 4.57053464e-03, 2.22276194e-02, 1.59968571e-03, - 4.74753501e-02, 9.57052289e-03, -2.72297788e-03, 7.60147678e-04, - -2.24030127e-03, -2.25223735e-03, -1.58726069e-03, 1.44355725e-02, - -2.65302865e-03, 1.43794907e-02, -7.26650570e-04, 1.72575124e-02, - -2.26407370e-02, -3.96339977e-05, -2.79391904e-02, 1.96495129e-04, - -3.42417083e-02, 5.69135891e-04, -2.83188289e-02, -3.90851596e-03, + [[1.14020750e-01, -6.87045764e-03, -1.14232450e-03, 4.11661189e-02, + -1.63423278e-01, 1.42638757e-03, -1.02391114e-01, -3.10274637e-03, + -4.14102976e-02, -7.65200269e-05, -1.22479669e-03, 1.44271152e-04, + 1.02909411e-03, -2.53524591e-02, 8.46680097e-04, -2.91525407e-02, + 8.15595775e-02, 2.62140564e-03, 8.88886643e-02, 3.47353540e-03, + 6.60543440e-02, 4.57053464e-03, 2.22276194e-02, 1.59968571e-03, + 4.74753501e-02, 9.57052289e-03, -2.72297788e-03, 7.60147678e-04, + -2.24030127e-03, -2.25223735e-03, -1.58726069e-03, 1.44355725e-02, + -2.65302865e-03, 1.43794907e-02, -7.26650570e-04, 1.72575124e-02, + -2.26407370e-02, -3.96339977e-05, -2.79391904e-02, 1.96495129e-04, + -3.42417083e-02, 5.69135891e-04, -2.83188289e-02, -3.90851596e-03, -9.86179763e-03, -1.53911048e-03, -1.44494506e-02, -3.56561464e-03, -2.28537982e-02]]], - [[[7.26406236e-02, 7.02981930e-03, -1.87881430e-03, -1.61256813e-02, - -2.64719814e-02, 1.79623468e-02, -2.81063801e-02, 6.54603124e-03, - 1.45408540e-02, -1.01358115e-02, 4.89154209e-04, -2.90726497e-03, - -5.52737463e-04, 8.21197934e-03, -1.38576779e-03, 1.91516929e-02, - 1.90580115e-02, 1.39116265e-03, 3.21379718e-02, -3.41384306e-02, - 1.56080311e-02, 2.17425083e-03, -1.61882394e-02, 2.30982768e-03, - -4.68867105e-03, -6.71668742e-04, -1.23223614e-04, 2.34797057e-03, - 1.21306862e-03, 1.70700484e-03, 2.32198084e-03, -5.41906023e-03, - 4.33958103e-03, -8.72831333e-03, 4.65559442e-04, -1.77434743e-02, - -3.67747484e-03, 2.22363618e-03, -6.92884020e-03, 1.54291482e-03, - -1.17597007e-02, -7.95446846e-04, -9.29108890e-05, 6.00740874e-03, - -8.13945278e-04, 6.46410176e-04, 3.25235814e-03, 5.35093467e-03, + [[[7.26406236e-02, 7.02981930e-03, -1.87881430e-03, -1.61256813e-02, + -2.64719814e-02, 1.79623468e-02, -2.81063801e-02, 6.54603124e-03, + 1.45408540e-02, -1.01358115e-02, 4.89154209e-04, -2.90726497e-03, + -5.52737463e-04, 8.21197934e-03, -1.38576779e-03, 1.91516929e-02, + 1.90580115e-02, 1.39116265e-03, 3.21379718e-02, -3.41384306e-02, + 1.56080311e-02, 2.17425083e-03, -1.61882394e-02, 2.30982768e-03, + -4.68867105e-03, -6.71668742e-04, -1.23223614e-04, 2.34797057e-03, + 1.21306862e-03, 1.70700484e-03, 2.32198084e-03, -5.41906023e-03, + 4.33958103e-03, -8.72831333e-03, 4.65559442e-04, -1.77434743e-02, + -3.67747484e-03, 2.22363618e-03, -6.92884020e-03, 1.54291482e-03, + -1.17597007e-02, -7.95446846e-04, -9.29108890e-05, 6.00740874e-03, + -8.13945278e-04, 6.46410176e-04, 3.25235814e-03, 5.35093467e-03, 1.94505156e-03]], - [[9.81325940e-02, 1.52165204e-02, 1.24419995e-03, -1.28939516e-03, - -7.51292526e-02, -4.42181788e-03, -6.18455506e-02, 1.09682359e-02, + [[9.81325940e-02, 1.52165204e-02, 1.24419995e-03, -1.28939516e-03, + -7.51292526e-02, -4.42181788e-03, -6.18455506e-02, 1.09682359e-02, 1.47464674e-03, -3.31177366e-02, -2.92832122e-03, -1.02011317e-02, - 3.94651358e-04, 2.74508980e-03, 3.55769102e-03, 4.53198680e-03, - 1.65537934e-02, 6.31128591e-04, 3.60410320e-02, -2.02268925e-02, - 3.71610204e-02, 1.50104427e-02, -2.10404071e-02, -6.66322009e-03, - -2.82880695e-02, 2.37099397e-02, 4.67322046e-04, 1.69550184e-02, - 7.12799550e-04, 5.70549692e-03, 8.82880626e-04, -3.97702624e-03, - -4.87365758e-03, 1.15246875e-03, -2.71843590e-03, 1.70523580e-03, - 8.22004027e-03, -2.16495266e-03, -1.35291214e-02, 1.12046159e-02, - -1.70833051e-02, 9.56991767e-03, -2.20789814e-02, -5.21140434e-03, - -5.63453210e-03, 9.14118089e-03, 1.07848139e-02, -2.44076209e-03, + 3.94651358e-04, 2.74508980e-03, 3.55769102e-03, 4.53198680e-03, + 1.65537934e-02, 6.31128591e-04, 3.60410320e-02, -2.02268925e-02, + 3.71610204e-02, 1.50104427e-02, -2.10404071e-02, -6.66322009e-03, + -2.82880695e-02, 2.37099397e-02, 4.67322046e-04, 1.69550184e-02, + 7.12799550e-04, 5.70549692e-03, 8.82880626e-04, -3.97702624e-03, + -4.87365758e-03, 1.15246875e-03, -2.71843590e-03, 1.70523580e-03, + 8.22004027e-03, -2.16495266e-03, -1.35291214e-02, 1.12046159e-02, + -1.70833051e-02, 9.56991767e-03, -2.20789814e-02, -5.21140434e-03, + -5.63453210e-03, 9.14118089e-03, 1.07848139e-02, -2.44076209e-03, 7.20372163e-03]], - [[7.82262155e-02, 4.89325589e-03, 9.69035921e-04, 2.48817773e-02, - -5.77909247e-02, 5.15396685e-03, -3.73859556e-02, -1.12652597e-02, + [[7.82262155e-02, 4.89325589e-03, 9.69035921e-04, 2.48817773e-02, + -5.77909247e-02, 5.15396685e-03, -3.73859556e-02, -1.12652597e-02, 1.46133523e-03, -1.93929247e-02, -1.66098338e-03, -5.40552612e-03, 1.90785108e-03, -1.27736588e-02, -2.63140561e-04, -1.64589849e-02, - -1.67926849e-02, 1.37226384e-02, 3.40286811e-02, 1.60456988e-02, + -1.67926849e-02, 1.37226384e-02, 3.40286811e-02, 1.60456988e-02, 3.66665622e-02, -4.95865979e-03, -8.81294607e-04, -9.27511404e-04, - -5.56159406e-03, 1.85807717e-02, 4.28397297e-05, 1.14440177e-02, - -8.44202079e-04, 2.83058074e-03, -2.77968111e-03, 7.24897169e-03, - -3.92705015e-03, 8.88560041e-03, -1.14740758e-03, 1.25959122e-02, - 2.28839936e-02, -6.30383042e-04, 6.79235888e-03, 2.74534443e-03, - -1.08549693e-02, 2.89646379e-03, -1.59780244e-02, -7.49181391e-03, + -5.56159406e-03, 1.85807717e-02, 4.28397297e-05, 1.14440177e-02, + -8.44202079e-04, 2.83058074e-03, -2.77968111e-03, 7.24897169e-03, + -3.92705015e-03, 8.88560041e-03, -1.14740758e-03, 1.25959122e-02, + 2.28839936e-02, -6.30383042e-04, 6.79235888e-03, 2.74534443e-03, + -1.08549693e-02, 2.89646379e-03, -1.59780244e-02, -7.49181391e-03, -9.15354118e-03, -1.60175812e-03, -1.81015916e-03, -7.74073756e-04, 5.32445729e-03]]], - [[[1.34755663e-02, 1.28413970e-02, 1.49938841e-03, -6.47475638e-03, - 2.52029988e-03, 4.88968824e-03, -5.42289754e-03, 1.60941684e-03, - -2.21852068e-03, 2.07330837e-03, 4.86536226e-03, 1.66978774e-04, - -6.71214109e-03, 1.27351141e-03, 9.82399935e-04, 2.22677827e-03, - 2.71719965e-03, 2.06634356e-03, 2.94252961e-03, -8.67336125e-03, - 3.55625935e-04, 1.36775123e-03, -2.08127408e-03, 2.52009045e-04, - -1.49716310e-04, -3.07585778e-04, 1.29094098e-04, -1.78327891e-04, - -3.55467543e-03, -1.03702801e-03, 5.01705916e-03, -1.46068291e-03, - 2.29931595e-03, -2.32968172e-03, 7.43670432e-04, -6.92986861e-04, - -1.00498143e-03, 2.95173070e-05, -1.53450925e-03, 3.21191049e-04, - -1.58863466e-03, 2.65016963e-03, -4.11708281e-04, -7.87924032e-04, - -8.33831423e-04, 1.07967881e-03, 2.63612367e-04, 9.87115809e-04, + [[[1.34755663e-02, 1.28413970e-02, 1.49938841e-03, -6.47475638e-03, + 2.52029988e-03, 4.88968824e-03, -5.42289754e-03, 1.60941684e-03, + -2.21852068e-03, 2.07330837e-03, 4.86536226e-03, 1.66978774e-04, + -6.71214109e-03, 1.27351141e-03, 9.82399935e-04, 2.22677827e-03, + 2.71719965e-03, 2.06634356e-03, 2.94252961e-03, -8.67336125e-03, + 3.55625935e-04, 1.36775123e-03, -2.08127408e-03, 2.52009045e-04, + -1.49716310e-04, -3.07585778e-04, 1.29094098e-04, -1.78327891e-04, + -3.55467543e-03, -1.03702801e-03, 5.01705916e-03, -1.46068291e-03, + 2.29931595e-03, -2.32968172e-03, 7.43670432e-04, -6.92986861e-04, + -1.00498143e-03, 2.95173070e-05, -1.53450925e-03, 3.21191049e-04, + -1.58863466e-03, 2.65016963e-03, -4.11708281e-04, -7.87924032e-04, + -8.33831423e-04, 1.07967881e-03, 2.63612367e-04, 9.87115809e-04, 1.21380213e-03]], - [[1.61478356e-02, 1.43728777e-02, 1.93974718e-03, -3.32623522e-03, - -1.17420191e-03, 3.21700718e-03, -8.75686341e-03, 1.37117517e-03, - -4.32829335e-03, -4.42773363e-03, 2.38329302e-03, -3.33902113e-03, - -4.12497515e-03, -6.89667443e-05, 2.41897821e-03, -4.26942050e-03, - -1.74191884e-03, 1.81570155e-03, 3.17867090e-03, -5.10889033e-03, - 4.49678511e-03, -2.60849304e-04, 1.84832760e-05, 4.61673598e-04, - -1.52379129e-03, 9.91300908e-04, 9.15433033e-04, 2.68803221e-03, - -2.02188630e-03, 2.73531219e-03, 3.98791591e-03, -1.49808319e-03, - 3.31828949e-04, 6.51188325e-04, -1.63160178e-03, 2.33668240e-03, - 2.93338642e-03, -5.06191672e-05, -1.03964494e-04, 3.87231978e-04, - -8.64134554e-04, 2.56989460e-03, -2.72089308e-03, -5.28755828e-04, - -2.52279382e-03, 1.50003583e-03, 2.73936763e-04, -7.08454308e-04, + [[1.61478356e-02, 1.43728777e-02, 1.93974718e-03, -3.32623522e-03, + -1.17420191e-03, 3.21700718e-03, -8.75686341e-03, 1.37117517e-03, + -4.32829335e-03, -4.42773363e-03, 2.38329302e-03, -3.33902113e-03, + -4.12497515e-03, -6.89667443e-05, 2.41897821e-03, -4.26942050e-03, + -1.74191884e-03, 1.81570155e-03, 3.17867090e-03, -5.10889033e-03, + 4.49678511e-03, -2.60849304e-04, 1.84832760e-05, 4.61673598e-04, + -1.52379129e-03, 9.91300908e-04, 9.15433033e-04, 2.68803221e-03, + -2.02188630e-03, 2.73531219e-03, 3.98791591e-03, -1.49808319e-03, + 3.31828949e-04, 6.51188325e-04, -1.63160178e-03, 2.33668240e-03, + 2.93338642e-03, -5.06191672e-05, -1.03964494e-04, 3.87231978e-04, + -8.64134554e-04, 2.56989460e-03, -2.72089308e-03, -5.28755828e-04, + -2.52279382e-03, 1.50003583e-03, 2.73936763e-04, -7.08454308e-04, 1.54615336e-03]], - [[1.21073977e-02, 9.65647748e-03, 1.40513870e-03, 5.96978239e-03, - -5.83388001e-03, 1.18491063e-03, -7.42142705e-03, -1.61753221e-04, - 2.31758132e-03, -1.01991745e-02, 1.25833161e-03, -3.28820299e-03, + [[1.21073977e-02, 9.65647748e-03, 1.40513870e-03, 5.96978239e-03, + -5.83388001e-03, 1.18491063e-03, -7.42142705e-03, -1.61753221e-04, + 2.31758132e-03, -1.01991745e-02, 1.25833161e-03, -3.28820299e-03, 9.77997070e-04, -3.33627171e-03, -7.64652162e-04, -4.00396819e-03, - -4.28747212e-03, 2.06971554e-03, 4.22446175e-03, 1.39004952e-03, + -4.28747212e-03, 2.06971554e-03, 4.22446175e-03, 1.39004952e-03, 5.11484010e-03, -2.51029606e-04, -1.53415978e-03, -9.05476137e-04, - -3.72103599e-03, 1.62807768e-03, 2.07196281e-03, 5.50544884e-03, - 3.46282883e-04, 3.54274790e-03, -1.15406623e-03, 1.05391112e-03, - -1.00296491e-03, 9.83314237e-04, -5.52227484e-04, -2.60564565e-04, - 3.63180002e-03, 1.15833727e-03, 2.31403574e-03, 6.13974106e-04, + -3.72103599e-03, 1.62807768e-03, 2.07196281e-03, 5.50544884e-03, + 3.46282883e-04, 3.54274790e-03, -1.15406623e-03, 1.05391112e-03, + -1.00296491e-03, 9.83314237e-04, -5.52227484e-04, -2.60564565e-04, + 3.63180002e-03, 1.15833727e-03, 2.31403574e-03, 6.13974106e-04, -7.09616975e-04, -1.19250538e-04, -3.14249677e-03, -2.51516937e-04, - -7.67022089e-04, -3.86918453e-04, 7.29633770e-04, -1.56990164e-04, + -7.67022089e-04, -3.86918453e-04, 7.29633770e-04, -1.56990164e-04, 9.27285029e-04]]]]) # Bingham fit fodf_3x3_bingham = np.array([[ - [[[0.8865956, -1.01305598, 0.3570712, 4.31487453, -4.09539035, + [[[0.8865956, -1.01305598, 0.3570712, 4.31487453, -4.09539035, 1.50380415, -1.08597013], - [0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.]]], - - [[[0.85683347, -3.44257945, 0.24681098, 2.7063325, 3.0255078, - -0.39132021, 3.88427259], - [0.09353532, -2.07887656, -0.23646431, 2.54845777, 3.35888907, - -12.87836086, 1.54502972], - [0., 0., 0., 0., 0., 0., 0.]]], - - [[[0.95034369, 3.82616343, 0.66640788, 2.10370133, -2.29257285, - -0.44777004, 4.31152262], - [0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.]]]], - - [[[[0.18908378, -1.3861731, 3.03919906, -1.22601547, 1.29995979, - 2.13770894, 3.82944086], - [0.12418191, 2.74008198, -0.06444388, -1.57097804, 3.1965494, - 0.87367118, 5.53954571], - [0.12104527, 3.61536971, -2.68105852, -1.64809093, 3.14649051, + [0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.]]], + + [[[0.85683347, -3.44257945, 0.24681098, 2.7063325, 3.0255078, + -0.39132021, 3.88427259], + [0.09353532, -2.07887656, -0.23646431, 2.54845777, 3.35888907, + -12.87836086, 1.54502972], + [0., 0., 0., 0., 0., 0., 0.]]], + + [[[0.95034369, 3.82616343, 0.66640788, 2.10370133, -2.29257285, + -0.44777004, 4.31152262], + [0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.]]]], + + [[[[0.18908378, -1.3861731, 3.03919906, -1.22601547, 1.29995979, + 2.13770894, 3.82944086], + [0.12418191, 2.74008198, -0.06444388, -1.57097804, 3.1965494, + 0.87367118, 5.53954571], + [0.12104527, 3.61536971, -2.68105852, -1.64809093, 3.14649051, 4.40167953, -0.25813741]]], - [[[0.27522583, -3.48810133, 0.8774683, 1.78047152, 2.43352883, - -0.35692355, 4.94340074], - [0.20108258, 4.56868628, 3.66376182, 0.03988314, -3.4448441, - 4.25196219, 4.01761794], - [0.16118525, 1.50428971, 2.4830825, 0.17002741, 1.20848685, - -1.22853768, 7.24969027]]], + [[[0.27522583, -3.48810133, 0.8774683, 1.78047152, 2.43352883, + -0.35692355, 4.94340074], + [0.20108258, 4.56868628, 3.66376182, 0.03988314, -3.4448441, + 4.25196219, 4.01761794], + [0.16118525, 1.50428971, 2.4830825, 0.17002741, 1.20848685, + -1.22853768, 7.24969027]]], - [[[0.2612878, 4.36835975, 2.85751351, 2.79826519, -2.49064595, - -1.38304902, 5.30046931], - [0.23787563, -0.23406115, -0.25916585, 3.62612255, 5.29582601, - -3.54583996, 0.08841027], - [0.14373517, 2.05462185, 1.39279777, 0.56035226, 3.45503345, - -6.17007375, 2.66774636]]]], + [[[0.2612878, 4.36835975, 2.85751351, 2.79826519, -2.49064595, + -1.38304902, 5.30046931], + [0.23787563, -0.23406115, -0.25916585, 3.62612255, 5.29582601, + -3.54583996, 0.08841027], + [0.14373517, 2.05462185, 1.39279777, 0.56035226, 3.45503345, + -6.17007375, 2.66774636]]]], - [[[[0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.]]], + [[[[0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.]]], - [[[0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.]]], + [[[0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.]]], - [[[0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.]]]]]) + [[[0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 0.]]]]]) # Bingham array as spherical functions fodf_3x3_bingham_sf = np.array([[[[ @@ -530,46 +530,46 @@ fodf_3x3_bingham_peaks = np.array([ [[[[-0.34397547, -0.93897366, -0.00305593], - [0., 0., 0.], - [0., 0., 0.]]], + [0., 0., 0.], + [0., 0., 0.]]], - [[[0.09314371, 0.9952668, 0.02771716], - [0.73461195, 0.26645789, 0.62397555], - [0., 0., 0.]]], + [[[0.09314371, 0.9952668, 0.02771716], + [0.73461195, 0.26645789, 0.62397555], + [0., 0., 0.]]], [[[0.17614917, -0.98432624, -0.0085624], - [0., 0., 0.], - [0., 0., 0.]]]], + [0., 0., 0.], + [0., 0., 0.]]]], - [[[[0.87605833, 0.22821013, -0.42478458], - [0.04979924, -0.99058828, 0.12749469], - [0.30605549, -0.16378166, 0.9378196]]], + [[[[0.87605833, 0.22821013, -0.42478458], + [0.04979924, -0.99058828, 0.12749469], + [0.30605549, -0.16378166, 0.9378196]]], - [[[0.2244257, 0.97366253, -0.04017949], - [0.3659656, -0.46513121, 0.80605344], + [[[0.2244257, 0.97366253, -0.04017949], + [0.3659656, -0.46513121, 0.80605344], [0.84032072, -0.49375814, -0.22374984]]], - [[[0.53356359, -0.84522153, 0.03017378], - [0.5527576, 0.82792627, 0.09485319], + [[[0.53356359, -0.84522153, 0.03017378], + [0.5527576, 0.82792627, 0.09485319], [0.37296054, -0.18433035, -0.90935293]]]], - [[[[0., 0., 0.], - [0., 0., 0.], - [0., 0., 0.]]], + [[[[0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.]]], - [[[0., 0., 0.], - [0., 0., 0.], - [0., 0., 0.]]], + [[[0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.]]], - [[[0., 0., 0.], - [0., 0., 0.], - [0., 0., 0.]]]]] + [[[0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.]]]]] ) # 3D array with slices 2-7 with values from 1-6 (with borders of 0) ref_in_labels = np.zeros((10, 10, 10), dtype=np.uint16) for i in range(2, 8): - ref_in_labels[2:8, 2:8, i] = i-1 + ref_in_labels[2:8, 2:8, i] = i - 1 # Same array with all odd slices are set to 0 ref_out_labels = copy.deepcopy(ref_in_labels) diff --git a/scilpy/tracking/tracker.py b/scilpy/tracking/tracker.py index 79aed3a26..973c0c258 100644 --- a/scilpy/tracking/tracker.py +++ b/scilpy/tracking/tracker.py @@ -569,8 +569,6 @@ def _track(self): GPU streamlines generator yielding streamlines with corresponding seed positions one by one. """ - t0 = perf_counter() - # Convert theta to cos(theta) max_cos_theta = np.cos(np.deg2rad(self.theta)) @@ -594,25 +592,29 @@ def _track(self): 'true' if self.sh_interp_nn else 'false') # Create CL program - n_input_params = 8 - n_output_params = 2 - cl_manager = CLManager(cl_kernel, n_input_params, n_output_params) + cl_manager = CLManager(cl_kernel) # Input buffers # Constant input buffers - cl_manager.add_input_buffer(0, self.sh) - cl_manager.add_input_buffer(1, self.sphere.vertices) + cl_manager.add_input_buffer('sh', self.sh) + cl_manager.add_input_buffer('vertices', self.sphere.vertices) sh_order = find_order_from_nb_coeff(self.sh) B_mat = sh_to_sf_matrix(self.sphere, sh_order, self.sh_basis, return_inv=False, legacy=self.is_legacy) - cl_manager.add_input_buffer(2, B_mat) + cl_manager.add_input_buffer('b_matrix', B_mat) fodf_max = self._get_max_amplitudes(B_mat) - cl_manager.add_input_buffer(3, fodf_max) - cl_manager.add_input_buffer(4, self.mask.astype(np.float32)) + cl_manager.add_input_buffer('max_amplitudes', fodf_max) + cl_manager.add_input_buffer('mask', self.mask.astype(np.float32)) + + cl_manager.add_input_buffer('max_cos_theta', max_cos_theta) + + cl_manager.add_input_buffer('seeds') + cl_manager.add_input_buffer('randvals') - cl_manager.add_input_buffer(5, max_cos_theta) + cl_manager.add_output_buffer('out_strl') + cl_manager.add_output_buffer('out_lengths') # Generate streamlines in batches for seed_batch in self.seed_batches: @@ -624,14 +626,16 @@ def _track(self): self.max_strl_points)) # Update buffers - cl_manager.add_input_buffer(6, seed_batch) - cl_manager.add_input_buffer(7, rand_vals) + cl_manager.update_input_buffer('seeds', seed_batch) + cl_manager.update_input_buffer('randvals', rand_vals) # output streamlines buffer - cl_manager.add_output_buffer(0, (len(seed_batch), + cl_manager.update_output_buffer('out_strl', + (len(seed_batch), self.max_strl_points, 3)) # output streamlines length buffer - cl_manager.add_output_buffer(1, (len(seed_batch), 1)) + cl_manager.update_output_buffer('out_lengths', + (len(seed_batch), 1)) # Run the kernel tracks, n_points = cl_manager.run((len(seed_batch), 1, 1)) diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py index a3a31d5ca..48c4cbb53 100755 --- a/scripts/scil_sh_to_aodf.py +++ b/scripts/scil_sh_to_aodf.py @@ -4,15 +4,15 @@ Script to estimate asymmetric ODFs (aODFs) from a spherical harmonics image. Two methods are available: - * Angle-aware bilateral filtering [1] is an extension of bilateral - filtering considering the angular distance between sphere directions - for filtering 5-dimensional spatio-angular images. + * Unified filtering [1] combines four asymmetric filtering methods into + a single equation and relies on a combination of four gaussian filters. * Cosine filtering [2] is a simpler implementation using cosine distance for assigning weights to neighbours. -Angle-aware bilateral filtering can be performed on the GPU using pyopencl by -specifying --use_gpu. Make sure you have pyopencl installed to use this option. -Otherwise, the filtering runs entirely on the CPU. +Unified filtering can be accelerated using OpenCL with the option --use_opencl. +Make sure you have pyopencl installed before using this option. By default, the +OpenCL program will run on the cpu. To use a gpu instead, also specify the +option --device gpu. """ import argparse @@ -27,13 +27,13 @@ from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg, assert_inputs_exist, add_sh_basis_args, assert_outputs_exist, parse_sh_basis_arg) -from scilpy.denoise.asym_filtering import (cosine_filtering, - angle_aware_bilateral_filtering) +from scilpy.denoise.asym_filtering import (cosine_filtering, unified_filtering) EPILOG = """ -[1] Poirier et al, 2022, "Intuitive Angle-Aware Bilateral Filtering Revealing - Asymmetric Fiber ODF for Improved Tractography", ISMRM 2022 (abstract 3552) +[1] Poirier and Descoteaux, 2024, "A Unified Filtering Method for Estimating + Asymmetric Orientation Distribution Functions", Neuroimage, vol. 287, + https://doi.org/10.1016/j.neuroimage.2024.120516 [2] Poirier et al, 2021, "Investigating the Occurrence of Asymmetric Patterns in White Matter Fiber Orientation Distribution Functions", ISMRM 2021 @@ -55,39 +55,58 @@ def _build_arg_parser(): add_sh_basis_args(p) - p.add_argument('--sphere', default='repulsion724', + p.add_argument('--sphere', default='repulsion200', choices=sorted(SPHERE_FILES.keys()), help='Sphere used for the SH to SF projection. ' '[%(default)s]') - p.add_argument('--method', default='bilateral', - choices=['bilateral', 'cosine'], - help='Method for estimating asymmetric ODFs ' - '[%(default)s].\nOne of:\n' - ' \'bilateral\': Angle-aware bilateral ' - 'filtering [1].\n' - ' \'cosine\' : Cosine-based filtering [2].') + p.add_argument('--method', default='unified', + choices=['unified', 'cosine'], + help="Method for estimating asymmetric ODFs " + "[%(default)s].\nOne of:\n" + " 'unified': Unified filtering [1].\n" + " 'cosine' : Cosine-based filtering [2].") shared_group = p.add_argument_group('Shared filter arguments') shared_group.add_argument('--sigma_spatial', default=1.0, type=float, help='Standard deviation for spatial distance.' ' [%(default)s]') - trilateral_group = p.add_argument_group('Angle-aware bilateral arguments') - trilateral_group.add_argument('--sigma_angular', default=1.0, type=float, - help='Standard deviation for angular ' - 'distance. [%(default)s]') - trilateral_group.add_argument('--sigma_range', default=1.0, type=float, - help='Standard deviation for range filter.' - ' [%(default)s]') + unified_group = p.add_argument_group('Unified filter arguments') + unified_group.add_argument('--sigma_align', default=0.8, type=float, + help='Standard deviation for alignment ' + 'filter. [%(default)s]') + unified_group.add_argument('--sigma_range', default=0.2, type=float, + help='Standard deviation for range filter\n' + '*relative to SF range of image*. ' + '[%(default)s]') + unified_group.add_argument('--sigma_angle', type=float, + help='Standard deviation for angular filter\n' + '(disabled by default).') + unified_group.add_argument('--disable_spatial', action='store_true', + help='Disable spatial filtering.') + unified_group.add_argument('--disable_align', action='store_true', + help='Disable alignment filtering.') + unified_group.add_argument('--disable_range', action='store_true', + help='Disable range filtering.') + unified_group.add_argument('--include_center', action='store_true', + help='Include center voxel in neighourhood.') + unified_group.add_argument('--win_hwidth', type=int, + help='Filtering window half-width. Defaults to ' + '3*sigma_spatial.') cosine_group = p.add_argument_group('Cosine filter arguments') cosine_group.add_argument('--sharpness', default=1.0, type=float, - help='Specify sharpness factor to use for' - ' weighted average. [%(default)s]') + help='Specify sharpness factor to use for\n' + 'weighted average. [%(default)s]') - p.add_argument('--use_gpu', action='store_true', - help='Use GPU for computation.') + p.add_argument('--device', choices=['cpu', 'gpu'], default='cpu', + help='Device to use for execution. [%(default)s]') + p.add_argument('--use_opencl', action='store_true', + help='Accelerate code using OpenCL (requires pyopencl\n' + 'and a working OpenCL implementation).') + p.add_argument('--patch_size', type=int, default=40, + help='OpenCL patch size. [%(default)s]') add_verbose_arg(p) add_overwrite_arg(p) @@ -99,8 +118,15 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - if args.use_gpu and args.method == 'cosine': - parser.error('Option --use_gpu is not supported for cosine filtering.') + if args.device == 'gpu' and not args.use_opencl: + logging.warning("Device 'gpu' chosen but --use_opencl not specified. " + "Proceeding with use_opencl=True.") + # force use_opencl if option gpu is chosen + args.use_opencl = True + + if args.use_opencl and args.method == 'cosine': + parser.error('Option --use_opencl is not supported' + ' for cosine filtering.') outputs = [args.out_sh] if args.out_sym: @@ -117,17 +143,23 @@ def main(): t0 = time.perf_counter() logging.info('Filtering SH image.') - if args.method == 'bilateral': - asym_sh = angle_aware_bilateral_filtering( - data, sh_order=sh_order, - sh_basis=sh_basis, - in_full_basis=full_basis, - is_legacy=is_legacy, + if args.method == 'unified': + sigma_align = None if args.disable_align else args.sigma_align + sigma_range = None if args.disable_range else args.sigma_range + sigma_spatial = None if args.disable_spatial else args.sigma_spatial + + asym_sh = unified_filtering( + data, sh_order=sh_order, sh_basis=sh_basis, + is_legacy=is_legacy, full_basis=full_basis, sphere_str=args.sphere, - sigma_spatial=args.sigma_spatial, - sigma_angular=args.sigma_angular, - sigma_range=args.sigma_range, - use_gpu=args.use_gpu) + sigma_spatial=sigma_spatial, + sigma_align=sigma_align, + sigma_angle=args.sigma_angle, + rel_sigma_range=sigma_range, + win_hwidth=args.win_hwidth, + exclude_center=not args.include_center, + device_type=args.device, + use_opencl=args.use_opencl) else: # args.method == 'cosine' asym_sh = cosine_filtering( data, sh_order=sh_order, diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py index 8c3f683c9..af24322a3 100644 --- a/scripts/tests/test_sh_to_aodf.py +++ b/scripts/tests/test_sh_to_aodf.py @@ -8,128 +8,125 @@ import numpy as np import pytest -from scilpy import SCILPY_HOME -from scilpy.io.fetcher import fetch_data, get_testing_files_dict +from scilpy.io.dvc import pull_test_case_package +from scilpy.gpuparallel.opencl_utils import have_opencl # If they already exist, this only takes 5 seconds (check md5sum) -fetch_data(get_testing_files_dict(), keys=['fodf_filtering.zip']) -data_path = os.path.join(SCILPY_HOME, 'fodf_filtering') +test_data_root = pull_test_case_package("aodf") tmp_dir = tempfile.TemporaryDirectory() -@pytest.fixture -def mock_filtering(mocker, out_fodf): - def _mock(*args, **kwargs): - img = nib.load(out_fodf) - return img.get_fdata().astype(np.float32) - - script = 'scil_sh_to_aodf' - filtering_fn = "angle_aware_bilateral_filtering" - return mocker.patch("scripts.{}.{}".format(script, filtering_fn), - side_effect=_mock, create=True) - - def test_help_option(script_runner): ret = script_runner.run('scil_sh_to_aodf.py', '--help') assert ret.success -@pytest.mark.parametrize("in_fodf,out_fodf", - [[os.path.join(data_path, 'fodf_descoteaux07_sub.nii.gz'), - os.path.join(data_path, 'fodf_descoteaux07_sub_full.nii.gz')]]) -def test_asym_basis_output(script_runner, mock_filtering, in_fodf, out_fodf): +@pytest.mark.parametrize("in_fodf,expected_fodf", [ + [os.path.join(test_data_root, "fodf_descoteaux07_sub.nii.gz"), + os.path.join(test_data_root, + "fodf_descoteaux07_sub_unified_asym.nii.gz")]]) +def test_asym_basis_output_gpu(script_runner, in_fodf, expected_fodf): os.chdir(os.path.expanduser(tmp_dir.name)) ret = script_runner.run('scil_sh_to_aodf.py', in_fodf, 'out_fodf1.nii.gz', '--sphere', 'repulsion100', - '--sigma_angular', '1.0', + '--sigma_align', '0.8', '--sigma_spatial', '1.0', - '--sigma_range', '1.0', - '--sh_basis', 'descoteaux07', '-f', + '--sigma_range', '0.2', + '--sigma_angle', '0.06', + '--use_opencl', + '--device', 'gpu', + '--sh_basis', 'descoteaux07_legacy', '-f', + '--include_center', print_result=True, shell=True) - assert ret.success - mock_filtering.assert_called_once() - - ret_fodf = nib.load("out_fodf1.nii.gz") - test_fodf = nib.load(out_fodf) - assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata()) - - -@pytest.mark.parametrize("in_fodf,out_fodf,sym_fodf", - [[os.path.join(data_path, "fodf_descoteaux07_sub.nii.gz"), - os.path.join(data_path, "fodf_descoteaux07_sub_full.nii.gz"), - os.path.join(data_path, "fodf_descoteaux07_sub_sym.nii.gz")]]) -def test_sym_basis_output( - script_runner, mock_filtering, in_fodf, out_fodf, sym_fodf): + if have_opencl: + # if we have opencl the script should not raise an error + assert ret.success + + # output should be close to expected (but not exactly equal because + # the python implementation is float64 while gpu is float32) + ret_fodf = nib.load("out_fodf1.nii.gz") + test_fodf = nib.load(expected_fodf) + assert np.allclose(ret_fodf.get_fdata(), + test_fodf.get_fdata(), + atol=1e-6) + else: + # if we don't have opencl the script should have raised an error + assert not ret.success + + +@pytest.mark.parametrize("in_fodf,expected_fodf", [ + [os.path.join(test_data_root, "fodf_descoteaux07_sub.nii.gz"), + os.path.join(test_data_root, + "fodf_descoteaux07_sub_unified_asym.nii.gz")]]) +def test_asym_basis_output(script_runner, in_fodf, expected_fodf): os.chdir(os.path.expanduser(tmp_dir.name)) ret = script_runner.run('scil_sh_to_aodf.py', - in_fodf, - 'out_fodf2.nii.gz', - '--out_sym', 'out_sym.nii.gz', + in_fodf, 'out_fodf1.nii.gz', '--sphere', 'repulsion100', - '--sigma_angular', '1.0', + '--sigma_align', '0.8', '--sigma_spatial', '1.0', - '--sigma_range', '1.0', - '--sh_basis', 'descoteaux07', '-f', + '--sigma_range', '0.2', + '--sigma_angle', '0.06', + '--device', 'cpu', + '--sh_basis', 'descoteaux07_legacy', '-f', + '--include_center', print_result=True, shell=True) assert ret.success - mock_filtering.assert_called_once() - ret_sym_fodf = nib.load("out_sym.nii.gz") - test_sym_fodf = nib.load(sym_fodf) - assert np.allclose(ret_sym_fodf.get_fdata(), test_sym_fodf.get_fdata()) + ret_fodf = nib.load("out_fodf1.nii.gz") + test_fodf = nib.load(expected_fodf) + assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata()) -@pytest.mark.parametrize("in_fodf,out_fodf", - [[os.path.join(data_path, "fodf_descoteaux07_sub_full.nii.gz"), - os.path.join(data_path, "fodf_descoteaux07_sub_twice.nii.gz")]]) -def test_asym_input(script_runner, mock_filtering, in_fodf, out_fodf): +@pytest.mark.parametrize("in_fodf,expected_fodf", [ + [os.path.join(test_data_root, + "fodf_descoteaux07_sub_unified_asym.nii.gz"), + os.path.join(test_data_root, + "fodf_descoteaux07_sub_unified_asym_twice.nii.gz")]]) +def test_asym_input(script_runner, in_fodf, expected_fodf): os.chdir(os.path.expanduser(tmp_dir.name)) ret = script_runner.run('scil_sh_to_aodf.py', - in_fodf, - 'out_fodf3.nii.gz', + in_fodf, 'out_fodf1.nii.gz', '--sphere', 'repulsion100', - '--sigma_angular', '1.0', + '--sigma_align', '0.8', '--sigma_spatial', '1.0', - '--sigma_range', '1.0', - '--sh_basis', 'descoteaux07', '-f', + '--sigma_range', '0.2', + '--sigma_angle', '0.06', + '--device', 'cpu', + '--sh_basis', 'descoteaux07_legacy', '-f', + '--include_center', print_result=True, shell=True) assert ret.success - mock_filtering.assert_called_once() - - ret_fodf = nib.load("out_fodf3.nii.gz") - test_fodf = nib.load(out_fodf) + + ret_fodf = nib.load("out_fodf1.nii.gz") + test_fodf = nib.load(expected_fodf) assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata()) -@pytest.mark.parametrize("in_fodf,out_fodf", - [[os.path.join(data_path, 'fodf_descoteaux07_sub.nii.gz'), - os.path.join(data_path, 'fodf_descoteaux07_sub_full.nii.gz')]]) -def test_cosine_method(script_runner, mock_filtering, in_fodf, out_fodf): +@pytest.mark.parametrize("in_fodf,out_fodf", [ + [os.path.join(test_data_root, 'fodf_descoteaux07_sub.nii.gz'), + os.path.join(test_data_root, + 'fodf_descoteaux07_sub_cosine_asym.nii.gz')]]) +def test_cosine_method(script_runner, in_fodf, out_fodf): os.chdir(os.path.expanduser(tmp_dir.name)) ret = script_runner.run('scil_sh_to_aodf.py', in_fodf, 'out_fodf1.nii.gz', '--sphere', 'repulsion100', - '--method', 'cosine', - '--sh_basis', 'descoteaux07', - '-f', + '--method', 'cosine', '-f', + '--sh_basis', 'descoteaux07_legacy', print_result=True, shell=True) assert ret.success - # method cosine is fast and not mocked - mock_filtering.assert_not_called() - ret_fodf = nib.load("out_fodf1.nii.gz") test_fodf = nib.load(out_fodf) - # We expect the output to be different from the - # one obtained with angle-aware bilateral filtering - assert not np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata()) + assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata())