From 62a1f4c3b1586d802d6eeac4e2b918fbb83fb52d Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Fri, 23 Feb 2024 11:54:31 -0500
Subject: [PATCH 01/19] Copy old branch

---
 scilpy/denoise/aodf_filter.cl  | 111 +++++++++++++++++++
 scilpy/denoise/unified_asym.py | 190 +++++++++++++++++++++++++++++++++
 2 files changed, 301 insertions(+)
 create mode 100644 scilpy/denoise/aodf_filter.cl
 create mode 100644 scilpy/denoise/unified_asym.py

diff --git a/scilpy/denoise/aodf_filter.cl b/scilpy/denoise/aodf_filter.cl
new file mode 100644
index 000000000..8dd0e00be
--- /dev/null
+++ b/scilpy/denoise/aodf_filter.cl
@@ -0,0 +1,111 @@
+#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 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 = ui;
+#else
+        for(int vi = 0; vi < N_DIRS; ++vi)
+        {
+#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)
+                    {
+#if EXCLUDE_SELF
+                        if(hi == win_hwidth && hj == win_hwidth && win_hwidth == hk)
+                        {
+                            continue;
+                        }
+#endif
+                        const int yvi_flat_ind =
+                            get_flat_index(x_ind + hi, y_ind + hj,
+                                           z_ind + hk, vi, 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];
+
+                        const int uv_flat_ind = get_flat_index(ui, vi, 0, 0, N_DIRS,
+                                                               N_DIRS, 1);
+                        const float uv_weight = uv_filter[uv_flat_ind];
+
+                        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/unified_asym.py b/scilpy/denoise/unified_asym.py
new file mode 100644
index 000000000..1dd7fe689
--- /dev/null
+++ b/scilpy/denoise/unified_asym.py
@@ -0,0 +1,190 @@
+# -*- coding: utf-8 -*-
+from dipy.data import get_sphere
+from dipy.reconst.shm import sh_to_sf_matrix
+import numpy as np
+from numba import njit
+from scilpy.gpuparallel.opencl_utils import CLManager, CLKernel
+from itertools import product as iterprod
+import logging
+
+
+class AsymmetricFilter():
+    def __init__(self, sh_order, sh_basis, legacy, full_basis,
+                 sphere_str, sigma_spatial, sigma_align,
+                 sigma_angle, sigma_range, exclude_self=False,
+                 disable_spatial=False, disable_align=False,
+                 disable_angle=False, disable_range=False):
+        self.sh_order = sh_order
+        self.legacy = legacy
+        self.basis = sh_basis
+        self.full_basis = full_basis
+        self.sphere = get_sphere(sphere_str)
+        self.sigma_spatial = sigma_spatial
+        self.sigma_align = sigma_align
+        self.sigma_angle = sigma_angle
+        self.sigma_range = sigma_range
+        self.exclude_self = exclude_self
+        self.disable_range = disable_range
+        self.disable_angle = disable_angle
+
+        # won't need this if disable range
+        self.uv_filter = _build_uv_filter(self.sphere.vertices,
+                                          self.sigma_angle)
+        # sigma still controls the width of the filter
+        self.nx_filter = _build_nx_filter(self.sphere.vertices, sigma_spatial,
+                                          sigma_align, disable_spatial,
+                                          disable_align)
+        logging.info('Filter shape: {}'.format(self.nx_filter.shape[:-1]))
+
+        self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
+                                 self.basis, self.full_basis,
+                                 legacy=self.legacy, return_inv=False)
+        _, self.B_inv = sh_to_sf_matrix(self.sphere, self.sh_order,
+                                        self.basis, True, legacy=self.legacy,
+                                        return_inv=True)
+
+        # initialize gpu
+        self.cl_kernel = None
+        self.cl_manager = None
+        self._prepare_gpu()
+
+    def _prepare_gpu(self):
+        self.cl_kernel = CLKernel('filter', 'denoise', 'aodf_filter.cl')
+
+        self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
+        self.cl_kernel.set_define('SIGMA_RANGE',
+                                  '{}f'.format(self.sigma_range))
+        self.cl_kernel.set_define('N_DIRS', len(self.sphere.vertices))
+        self.cl_kernel.set_define('EXCLUDE_SELF', 'true' if self.exclude_self
+                                                  else 'false')
+        self.cl_kernel.set_define('DISABLE_ANGLE', 'true' if self.disable_angle
+                                                   else 'false')
+        self.cl_kernel.set_define('DISABLE_RANGE', 'true' if self.disable_range
+                                                   else 'false')
+        self.cl_manager = CLManager(self.cl_kernel, 3, 1)
+
+    def __call__(self, sh_data, patch_size=40):
+        # Fill const GPU buffers
+        self.cl_manager.add_input_buffer(1, self.nx_filter)
+        self.cl_manager.add_input_buffer(2, self.uv_filter)
+
+        win_width = self.nx_filter.shape[0]
+        win_hwidth = win_width // 2
+        volume_shape = sh_data.shape[:-1]
+        sf_shape = tuple(np.append(volume_shape, len(self.sphere.vertices)))
+        padded_volume_shape = tuple(np.asarray(volume_shape) + win_width - 1)
+
+        out_sh = np.zeros(np.append(sh_data.shape[:-1], self.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 + self.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(self.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]]
+            # print(patch_in, sh_patch.shape)
+            sf_patch = np.dot(sh_patch, self.B)
+            self.cl_manager.add_input_buffer(0, sf_patch)
+            self.cl_manager.add_output_buffer(0, out_shape)
+            out_sf = self.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,
+                                                             self.B_inv)
+
+        return out_sh
+
+
+@njit(cache=True)
+def _build_uv_filter(directions, sigma_angle):
+    directions = np.ascontiguousarray(directions.astype(np.float32))
+    uv_weights = np.zeros((len(directions), len(directions)), dtype=np.float32)
+
+    # 1. precompute weights on angle
+    # c'est toujours les mêmes peu importe le voxel en question
+    for u_i, u in enumerate(directions):
+        uvec = np.reshape(np.ascontiguousarray(u), (1, 3))
+        weights = np.arccos(np.clip(np.dot(uvec, directions.T), -1.0, 1.0))
+        weights = _evaluate_gaussian(sigma_angle, weights)
+        weights /= np.sum(weights)
+        uv_weights[u_i] = weights  # each line sums to 1.
+
+    return uv_weights
+
+
+@njit(cache=True)
+def _build_nx_filter(directions, sigma_spatial, sigma_align,
+                     disable_spatial, disable_align):
+    directions = np.ascontiguousarray(directions.astype(np.float32))
+
+    half_width = int(round(3 * sigma_spatial))
+    nx_weights = np.zeros((2*half_width+1, 2*half_width+1,
+                           2*half_width+1, len(directions)),
+                          dtype=np.float32)
+
+    for i in range(-half_width, half_width+1):
+        for j in range(-half_width, half_width+1):
+            for k in range(-half_width, half_width+1):
+                dxy = np.array([[i, j, k]], dtype=np.float32)
+                len_xy = np.sqrt(dxy[0, 0]**2 + dxy[0, 1]**2 + dxy[0, 2]**2)
+
+                if disable_spatial:
+                    w_spatial = 1.0
+                else:
+                    # the length controls spatial weight
+                    w_spatial = _evaluate_gaussian(sigma_spatial, len_xy)
+
+                # the direction controls the align weight
+                if i == j == k == 0 or disable_align:
+                    # hack for main direction to have maximal weight
+                    # w_align = np.ones((1, len(directions)), dtype=np.float32)
+                    w_align = np.zeros((1, len(directions)), dtype=np.float32)
+                else:
+                    dxy /= len_xy
+                    w_align = np.arccos(np.clip(np.dot(dxy, directions.T),
+                                                -1.0, 1.0))  # 1, N
+                w_align = _evaluate_gaussian(sigma_align, w_align)
+
+                nx_weights[half_width + i, half_width + j, half_width + k] =\
+                    w_align * w_spatial
+
+    # sur chaque u, le filtre doit sommer à 1
+    for ui in range(len(directions)):
+        w_sum = np.sum(nx_weights[..., ui])
+        nx_weights /= w_sum
+
+    return nx_weights
+
+
+@njit(cache=True)
+def _evaluate_gaussian(sigma, x):
+    # gaussian is not normalized
+    return np.exp(-x**2/(2*sigma**2))

From 424b274e1cacf4dcec840cbb39be247e0a490f12 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Fri, 23 Feb 2024 16:08:14 -0500
Subject: [PATCH 02/19] [WIP] Migrate code from aodf-toolkit

---
 scilpy/denoise/aodf_filter.cl      |  30 +--
 scilpy/denoise/asym_filtering.py   | 380 +++++++++++++----------------
 scilpy/denoise/unified_asym.py     | 166 +++++++------
 scilpy/gpuparallel/opencl_utils.py |  69 +++++-
 4 files changed, 340 insertions(+), 305 deletions(-)

diff --git a/scilpy/denoise/aodf_filter.cl b/scilpy/denoise/aodf_filter.cl
index 8dd0e00be..4cf707e02 100644
--- a/scilpy/denoise/aodf_filter.cl
+++ b/scilpy/denoise/aodf_filter.cl
@@ -22,6 +22,8 @@ float range_filter(const float x)
 __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
@@ -57,10 +59,14 @@ __kernel void filter(__global const float* sf_data,
         float w_norm = 0.0f;
         float tilde_psi_xui = 0.0f;
 #if DISABLE_ANGLE
-        const int vi = ui;
+        const int vi_in_image = ui;
 #else
-        for(int vi = 0; vi < N_DIRS; ++vi)
+        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)
             {
@@ -68,16 +74,10 @@ __kernel void filter(__global const float* sf_data,
                 {
                     for(int hk = 0; hk < WIN_WIDTH; ++hk)
                     {
-#if EXCLUDE_SELF
-                        if(hi == win_hwidth && hj == win_hwidth && win_hwidth == hk)
-                        {
-                            continue;
-                        }
-#endif
                         const int yvi_flat_ind =
-                            get_flat_index(x_ind + hi, y_ind + hj,
-                                           z_ind + hk, vi, x_pad_len,
-                                           y_pad_len, z_pad_len);
+                            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;
@@ -90,9 +90,11 @@ __kernel void filter(__global const float* sf_data,
                                                                     WIN_WIDTH);
                         const float nx_weight = nx_filter[y_in_nx_flat_ind];
 
-                        const int uv_flat_ind = get_flat_index(ui, vi, 0, 0, N_DIRS,
-                                                               N_DIRS, 1);
-                        const float uv_weight = uv_filter[uv_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;
diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index bb185290f..3f6ab6ee8 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -1,150 +1,186 @@
 # -*- 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,
-                                    sphere_str='repulsion724',
-                                    sigma_spatial=1.0, sigma_angular=1.0,
-                                    sigma_range=0.5, use_gpu=True):
-    """
-    Angle-aware bilateral filtering.
-
-    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.
-    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.
-    """
-    if use_gpu and have_opencl:
-        return angle_aware_bilateral_filtering_gpu(in_sh, sh_order,
-                                                   sh_basis, in_full_basis,
-                                                   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.')
-    else:
-        return angle_aware_bilateral_filtering_cpu(in_sh, sh_order,
-                                                   sh_basis, in_full_basis,
-                                                   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,
-                                        sphere_str='repulsion724',
-                                        sigma_spatial=1.0,
-                                        sigma_angular=1.0,
-                                        sigma_range=0.5):
-    """
-    Angle-aware bilateral filtering using OpenCL for GPU computing.
-
-    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.
-    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.
-
-    Returns
-    -------
-    out_sh: ndarray (x, y, z, ncoeffs)
-        Output SH coefficient array in full SH basis.
-    """
-    s_weights = _get_spatial_weights(sigma_spatial)
-    h_half_width = len(s_weights) // 2
-
-    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=sh_order,
-                                   basis_type=sh_basis,
-                                   full_basis=in_full_basis,
-                                   return_inv=False)
-
-    _, sf_to_sh_mat = sh_to_sf_matrix(sphere, sh_order=sh_order,
-                                      basis_type=sh_basis,
-                                      full_basis=True,
-                                      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]
+class AsymmetricFilter():
+    def __init__(self, sh_order, sh_basis, legacy, full_basis,
+                 sphere_str, sigma_spatial, sigma_align,
+                 sigma_angle, sigma_range, disable_spatial=False,
+                 disable_align=False, disable_angle=False,
+                 disable_range=False, device_type='gpu',
+                 j_invariance=False):
+        self.sh_order = sh_order
+        self.legacy = legacy
+        self.basis = sh_basis
+        self.full_basis = full_basis
+        self.sphere = get_sphere(sphere_str)
+        self.sigma_spatial = sigma_spatial
+        self.sigma_align = sigma_align
+        self.sigma_angle = sigma_angle
+        self.sigma_range = sigma_range
+        self.disable_range = disable_range
+        self.disable_angle = disable_angle
+        self.device_type = device_type
+
+        if device_type == 'gpu' and not have_opencl:
+            raise ValueError('device_type `gpu` requested but pyopencl'
+                             ' is not installed.\nPlease install pyopencl to'
+                             ' enable GPU acceleration.')
+
+        # won't need this if disable angle
+        self.uv_filter = self._build_uv_filter(self.sphere.vertices,
+                                               self.sigma_angle)
+
+        # sigma still controls the width of the filter
+        self.nx_filter = self._build_nx_filter(self.sphere.vertices, sigma_spatial,
+                                               sigma_align, disable_spatial,
+                                               disable_align, j_invariance)
+
+        self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
+                                 self.basis, self.full_basis,
+                                 legacy=self.legacy, return_inv=False)
+        _, self.B_inv = sh_to_sf_matrix(self.sphere, self.sh_order,
+                                        self.basis, True, legacy=self.legacy,
+                                        return_inv=True)
+
+        # initialize gpu
+        self.cl_kernel = None
+        self.cl_manager = None
+        self._prepare_gpu()
+
+    def _prepare_gpu(self):
+        self.cl_kernel = CLKernel('filter', 'filtering', 'aodf_filter.cl')
+
+        self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
+        self.cl_kernel.set_define('SIGMA_RANGE',
+                                  '{}f'.format(self.sigma_range))
+        self.cl_kernel.set_define('N_DIRS', len(self.sphere.vertices))
+        self.cl_kernel.set_define('DISABLE_ANGLE', 'true' if self.disable_angle
+                                                   else 'false')
+        self.cl_kernel.set_define('DISABLE_RANGE', 'true' if self.disable_range
+                                                   else 'false')
+        self.cl_manager = CLManager(self.cl_kernel, self.device_type)
+
+    def __call__(self, sh_data, patch_size=40):
+        uv_weights_offsets =\
+            np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter, axis=-1)))
+        v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
+                            (self.uv_filter.shape[0], 1))[self.uv_filter > 0.0]
+        flat_uv = self.uv_filter[self.uv_filter > 0.0]
+
+        # Prepare GPU buffers
+        self.cl_manager.add_input_buffer("sf_data")  # SF data not initialized
+        self.cl_manager.add_input_buffer("nx_filter", self.nx_filter)
+        self.cl_manager.add_input_buffer("uv_filter", flat_uv)
+        self.cl_manager.add_input_buffer("uv_weights_offsets", uv_weights_offsets)
+        self.cl_manager.add_input_buffer("v_indices", v_indices)
+        self.cl_manager.add_output_buffer("out_sf")  # SF not initialized yet
+
+        win_width = self.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], self.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 + self.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(self.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, self.B)
+            self.cl_manager.update_input_buffer("sf_data", sf_patch)
+            self.cl_manager.update_output_buffer("out_sf", out_shape)
+            out_sf = self.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,
+                                                             self.B_inv)
+
+        return out_sh
+
+    def _build_uv_filter(directions, sigma_angle):
+        dot = directions.dot(directions.T)
+        x = np.arccos(np.clip(dot, -1.0, 1.0))
+        weights = _evaluate_gaussian_distribution(sigma_angle, x)
+
+        mask = x > (3.0*sigma_angle)
+        weights[mask] = 0.0
+        weights /= np.sum(weights, axis=-1)
+        return weights
+
+    def _build_nx_filter(directions, sigma_spatial, sigma_align,
+                        disable_spatial=False, disable_align=False,
+                        j_invariance=False):
+        directions = directions.astype(np.float32)
+        half_width = int(round(3*sigma_spatial))
+        filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
+
+        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]
+
+        if disable_spatial:
+            w_spatial = np.ones(filter_shape)
+        else:
+            w_spatial = _evaluate_gaussian_distribution(sigma_spatial, distances)
+
+        cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
+        theta = np.arccos(cos_theta)
+        theta[half_width, half_width, half_width] = 0.0
+
+        if disable_align:
+            w_align = np.ones(np.append(filter_shape, (len(directions),)))
+        else:
+            w_align = _evaluate_gaussian_distribution(sigma_align, theta)
+
+        w = w_spatial[..., None] * w_align
+
+        if j_invariance:
+            w[half_width, half_width, half_width] = 0.0
+
+        # normalize
+        w /= np.sum(w, axis=(0,1,2), keepdims=True)
+
+        return w
 
 
 def angle_aware_bilateral_filtering_cpu(in_sh, sh_order=8,
@@ -259,74 +295,6 @@ def _get_window_directions(shape):
     return grid
 
 
-def _get_spatial_weights(sigma_spatial):
-    """
-    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].
-
-    Parameters
-    ----------
-    sigma_spatial: float
-        Standard deviation of spatial filter.
-
-    Returns
-    -------
-    spatial_weights: ndarray
-        Spatial filter.
-    """
-    shape = int(6 * sigma_spatial)
-    if shape % 2 == 0:
-        shape += 1
-    shape = (shape, shape, shape)
-
-    grid = _get_window_directions(shape)
-
-    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
-
-
-def _get_angular_weights(shape, sphere, sigma_angular):
-    """
-    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.
-
-    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.
-
-    Returns
-    -------
-    angular_weights: ndarray
-        Angular filter for each position and for each sphere direction.
-    """
-    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)
-
-    # normalize filter per direction
-    angular_weights /= np.sum(angular_weights, axis=(0, 1, 2))
-    return angular_weights
-
-
 def _correlate_spatial(image_u, h_filter, sigma_range):
     """
     Implementation of the correlation operation for anisotropic filtering.
diff --git a/scilpy/denoise/unified_asym.py b/scilpy/denoise/unified_asym.py
index 1dd7fe689..7a2f3d17a 100644
--- a/scilpy/denoise/unified_asym.py
+++ b/scilpy/denoise/unified_asym.py
@@ -2,7 +2,6 @@
 from dipy.data import get_sphere
 from dipy.reconst.shm import sh_to_sf_matrix
 import numpy as np
-from numba import njit
 from scilpy.gpuparallel.opencl_utils import CLManager, CLKernel
 from itertools import product as iterprod
 import logging
@@ -11,9 +10,10 @@
 class AsymmetricFilter():
     def __init__(self, sh_order, sh_basis, legacy, full_basis,
                  sphere_str, sigma_spatial, sigma_align,
-                 sigma_angle, sigma_range, exclude_self=False,
-                 disable_spatial=False, disable_align=False,
-                 disable_angle=False, disable_range=False):
+                 sigma_angle, sigma_range, disable_spatial=False,
+                 disable_align=False, disable_angle=False,
+                 disable_range=False, device_type='gpu',
+                 j_invariance=False):
         self.sh_order = sh_order
         self.legacy = legacy
         self.basis = sh_basis
@@ -23,17 +23,18 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis,
         self.sigma_align = sigma_align
         self.sigma_angle = sigma_angle
         self.sigma_range = sigma_range
-        self.exclude_self = exclude_self
         self.disable_range = disable_range
         self.disable_angle = disable_angle
+        self.device_type = device_type
 
-        # won't need this if disable range
+        # won't need this if disable angle
         self.uv_filter = _build_uv_filter(self.sphere.vertices,
                                           self.sigma_angle)
+
         # sigma still controls the width of the filter
         self.nx_filter = _build_nx_filter(self.sphere.vertices, sigma_spatial,
                                           sigma_align, disable_spatial,
-                                          disable_align)
+                                          disable_align, j_invariance)
         logging.info('Filter shape: {}'.format(self.nx_filter.shape[:-1]))
 
         self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
@@ -49,29 +50,36 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis,
         self._prepare_gpu()
 
     def _prepare_gpu(self):
-        self.cl_kernel = CLKernel('filter', 'denoise', 'aodf_filter.cl')
+        self.cl_kernel = CLKernel('filter', 'filtering', 'aodf_filter.cl')
 
         self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
         self.cl_kernel.set_define('SIGMA_RANGE',
                                   '{}f'.format(self.sigma_range))
         self.cl_kernel.set_define('N_DIRS', len(self.sphere.vertices))
-        self.cl_kernel.set_define('EXCLUDE_SELF', 'true' if self.exclude_self
-                                                  else 'false')
         self.cl_kernel.set_define('DISABLE_ANGLE', 'true' if self.disable_angle
                                                    else 'false')
         self.cl_kernel.set_define('DISABLE_RANGE', 'true' if self.disable_range
                                                    else 'false')
-        self.cl_manager = CLManager(self.cl_kernel, 3, 1)
+        self.cl_manager = CLManager(self.cl_kernel, self.device_type)
 
     def __call__(self, sh_data, patch_size=40):
-        # Fill const GPU buffers
-        self.cl_manager.add_input_buffer(1, self.nx_filter)
-        self.cl_manager.add_input_buffer(2, self.uv_filter)
+        uv_weights_offsets =\
+            np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter, axis=-1)))
+        v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
+                            (self.uv_filter.shape[0], 1))[self.uv_filter > 0.0]
+        flat_uv = self.uv_filter[self.uv_filter > 0.0]
+
+        # Prepare GPU buffers
+        self.cl_manager.add_input_buffer("sf_data")  # SF data not initialized
+        self.cl_manager.add_input_buffer("nx_filter", self.nx_filter)
+        self.cl_manager.add_input_buffer("uv_filter", flat_uv)
+        self.cl_manager.add_input_buffer("uv_weights_offsets", uv_weights_offsets)
+        self.cl_manager.add_input_buffer("v_indices", v_indices)
+        self.cl_manager.add_output_buffer("out_sf")  # SF not initialized yet
 
         win_width = self.nx_filter.shape[0]
         win_hwidth = win_width // 2
         volume_shape = sh_data.shape[:-1]
-        sf_shape = tuple(np.append(volume_shape, len(self.sphere.vertices)))
         padded_volume_shape = tuple(np.asarray(volume_shape) + win_width - 1)
 
         out_sh = np.zeros(np.append(sh_data.shape[:-1], self.B_inv.shape[-1]))
@@ -110,10 +118,10 @@ def __call__(self, sh_data, patch_size=40):
             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]]
-            # print(patch_in, sh_patch.shape)
+
             sf_patch = np.dot(sh_patch, self.B)
-            self.cl_manager.add_input_buffer(0, sf_patch)
-            self.cl_manager.add_output_buffer(0, out_shape)
+            self.cl_manager.update_input_buffer("sf_data", sf_patch)
+            self.cl_manager.update_output_buffer("out_sf", out_shape)
             out_sf = self.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],
@@ -123,68 +131,76 @@ def __call__(self, sh_data, patch_size=40):
         return out_sh
 
 
-@njit(cache=True)
 def _build_uv_filter(directions, sigma_angle):
-    directions = np.ascontiguousarray(directions.astype(np.float32))
-    uv_weights = np.zeros((len(directions), len(directions)), dtype=np.float32)
-
-    # 1. precompute weights on angle
-    # c'est toujours les mêmes peu importe le voxel en question
-    for u_i, u in enumerate(directions):
-        uvec = np.reshape(np.ascontiguousarray(u), (1, 3))
-        weights = np.arccos(np.clip(np.dot(uvec, directions.T), -1.0, 1.0))
-        weights = _evaluate_gaussian(sigma_angle, weights)
-        weights /= np.sum(weights)
-        uv_weights[u_i] = weights  # each line sums to 1.
+    dot = directions.dot(directions.T)
+    x = np.arccos(np.clip(dot, -1.0, 1.0))
+    weights = _evaluate_gaussian(sigma_angle, x)
 
-    return uv_weights
+    mask = x > (3.0*sigma_angle)
+    weights[mask] = 0.0
+    weights /= np.sum(weights, axis=-1)
+    return weights
 
 
-@njit(cache=True)
 def _build_nx_filter(directions, sigma_spatial, sigma_align,
-                     disable_spatial, disable_align):
-    directions = np.ascontiguousarray(directions.astype(np.float32))
-
-    half_width = int(round(3 * sigma_spatial))
-    nx_weights = np.zeros((2*half_width+1, 2*half_width+1,
-                           2*half_width+1, len(directions)),
-                          dtype=np.float32)
-
-    for i in range(-half_width, half_width+1):
-        for j in range(-half_width, half_width+1):
-            for k in range(-half_width, half_width+1):
-                dxy = np.array([[i, j, k]], dtype=np.float32)
-                len_xy = np.sqrt(dxy[0, 0]**2 + dxy[0, 1]**2 + dxy[0, 2]**2)
-
-                if disable_spatial:
-                    w_spatial = 1.0
-                else:
-                    # the length controls spatial weight
-                    w_spatial = _evaluate_gaussian(sigma_spatial, len_xy)
-
-                # the direction controls the align weight
-                if i == j == k == 0 or disable_align:
-                    # hack for main direction to have maximal weight
-                    # w_align = np.ones((1, len(directions)), dtype=np.float32)
-                    w_align = np.zeros((1, len(directions)), dtype=np.float32)
-                else:
-                    dxy /= len_xy
-                    w_align = np.arccos(np.clip(np.dot(dxy, directions.T),
-                                                -1.0, 1.0))  # 1, N
-                w_align = _evaluate_gaussian(sigma_align, w_align)
-
-                nx_weights[half_width + i, half_width + j, half_width + k] =\
-                    w_align * w_spatial
-
-    # sur chaque u, le filtre doit sommer à 1
-    for ui in range(len(directions)):
-        w_sum = np.sum(nx_weights[..., ui])
-        nx_weights /= w_sum
-
-    return nx_weights
-
-
-@njit(cache=True)
+                     disable_spatial=False, disable_align=False,
+                     j_invariance=False):
+    directions = directions.astype(np.float32)
+    half_width = int(round(3*sigma_spatial))
+    filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
+
+    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]
+
+    if disable_spatial:
+        w_spatial = np.ones(filter_shape)
+    else:
+        w_spatial = _evaluate_gaussian(sigma_spatial, distances)
+
+    cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
+    theta = np.arccos(cos_theta)
+    theta[half_width, half_width, half_width] = 0.0
+
+    if disable_align:
+        w_align = np.ones(np.append(filter_shape, (len(directions),)))
+    else:
+        w_align = _evaluate_gaussian(sigma_align, theta)
+
+    w = w_spatial[..., None] * w_align
+
+    if j_invariance:
+        w[half_width, half_width, half_width] = 0.0
+
+    # normalize
+    w /= np.sum(w, axis=(0,1,2), keepdims=True)
+
+    return w
+
+
+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
+
+
 def _evaluate_gaussian(sigma, x):
     # gaussian is not normalized
     return np.exp(-x**2/(2*sigma**2))
diff --git a/scilpy/gpuparallel/opencl_utils.py b/scilpy/gpuparallel/opencl_utils.py
index 4f093d7c6..691fd6824 100644
--- a/scilpy/gpuparallel/opencl_utils.py
+++ b/scilpy/gpuparallel/opencl_utils.py
@@ -8,6 +8,13 @@
 from dipy.utils.optpkg import optional_package
 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):
     """
@@ -25,7 +32,7 @@ class CLManager(object):
     n_outputs: int
         Number of output buffers for the kernel.
     """
-    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 +41,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 +54,12 @@ 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 +84,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
@@ -98,14 +113,31 @@ 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):
+        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.
 
@@ -119,9 +151,26 @@ def add_output_buffer(self, arg_pos, shape, dtype=np.float32):
             Data type for the output. It is recommended to keep
             float32 to avoid unexpected behaviour.
         """
+        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):
+        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):
         """

From 1289ab48ee0d332460d72c491b0564b33245f07d Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Mon, 26 Feb 2024 17:07:27 -0500
Subject: [PATCH 03/19] Working Python cpu, opencl gpu and opencl cpu
 implementations

---
 scilpy/denoise/asym_filtering.py   | 438 ++++++++++++++---------------
 scilpy/gpuparallel/opencl_utils.py |   6 +-
 scripts/scil_sh_to_aodf.py         |  61 ++--
 3 files changed, 267 insertions(+), 238 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index 3f6ab6ee8..8c2a73127 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -10,12 +10,48 @@
 
 
 class AsymmetricFilter():
+    """
+    Unified asymmetric filtering as described in Poirier et al, 2024.
+
+    Parameters
+    ----------
+    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'.
+    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
+        Standard deviation of spatial filter.
+    sigma_align: float
+        Standard deviation of alignment filter.
+    sigma_angle: float
+        Standard deviation of the angle filter.
+    rel_sigma_range: float
+        Standard deviation of the range filter, relative to the
+        range of SF amplitudes.
+    disable_spatial: bool, optional
+        Disable spatial filtering. Replaces the
+        gaussian weights by a constant value.
+    disable_align: bool, optional
+        Disable alignment filtering.
+    disable_angle: bool, optional
+        Disable angle filtering, which considers neighbour *sphere
+        directions* in the estimation of the a-ODF.
+    disable_range: bool, optional
+        Disable range filtering.
+    """
     def __init__(self, sh_order, sh_basis, legacy, full_basis,
                  sphere_str, sigma_spatial, sigma_align,
-                 sigma_angle, sigma_range, disable_spatial=False,
+                 sigma_angle, rel_sigma_range, disable_spatial=False,
                  disable_align=False, disable_angle=False,
-                 disable_range=False, device_type='gpu',
-                 j_invariance=False):
+                 disable_range=False, j_invariance=False,
+                 device_type='gpu', use_opencl=True):
         self.sh_order = sh_order
         self.legacy = legacy
         self.basis = sh_basis
@@ -24,24 +60,28 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis,
         self.sigma_spatial = sigma_spatial
         self.sigma_align = sigma_align
         self.sigma_angle = sigma_angle
-        self.sigma_range = sigma_range
+        self.rel_sigma_range = rel_sigma_range
         self.disable_range = disable_range
         self.disable_angle = disable_angle
+        self.disable_spatial = disable_spatial
+        self.disable_align = disable_align
         self.device_type = device_type
-
-        if device_type == 'gpu' and not have_opencl:
-            raise ValueError('device_type `gpu` requested but pyopencl'
-                             ' is not installed.\nPlease install pyopencl to'
-                             ' enable GPU acceleration.')
-
-        # won't need this if disable angle
-        self.uv_filter = self._build_uv_filter(self.sphere.vertices,
-                                               self.sigma_angle)
-
-        # sigma still controls the width of the filter
-        self.nx_filter = self._build_nx_filter(self.sphere.vertices, sigma_spatial,
-                                               sigma_align, disable_spatial,
-                                               disable_align, j_invariance)
+        self.use_opencl = use_opencl
+        self.j_invariance = j_invariance
+
+        if device_type not in ['cpu', 'gpu']:
+            raise ValueError('Invalid device type {}. '
+                             'Must be one of \{cpu, gpu\}')
+        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\'.')
+
+        # build filters
+        self.uv_filter = self._build_uv_filter()
+        self.nx_filter = self._build_nx_filter()
 
         self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
                                  self.basis, self.full_basis,
@@ -50,17 +90,17 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis,
                                         self.basis, True, legacy=self.legacy,
                                         return_inv=True)
 
-        # initialize gpu
-        self.cl_kernel = None
-        self.cl_manager = None
-        self._prepare_gpu()
+        if self.use_opencl:
+            # initialize opencl
+            self.cl_kernel = None
+            self.cl_manager = None
 
-    def _prepare_gpu(self):
-        self.cl_kernel = CLKernel('filter', 'filtering', 'aodf_filter.cl')
+    def _prepare_opencl(self, sigma_range):
+        self.cl_kernel = CLKernel('filter', 'denoise', 'aodf_filter.cl')
 
         self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
         self.cl_kernel.set_define('SIGMA_RANGE',
-                                  '{}f'.format(self.sigma_range))
+                                  '{}f'.format(sigma_range))
         self.cl_kernel.set_define('N_DIRS', len(self.sphere.vertices))
         self.cl_kernel.set_define('DISABLE_ANGLE', 'true' if self.disable_angle
                                                    else 'false')
@@ -68,7 +108,65 @@ def _prepare_gpu(self):
                                                    else 'false')
         self.cl_manager = CLManager(self.cl_kernel, self.device_type)
 
-    def __call__(self, sh_data, patch_size=40):
+    def _build_uv_filter(self):
+        directions = self.sphere.vertices
+        if not self.disable_angle:
+            dot = directions.dot(directions.T)
+            x = np.arccos(np.clip(dot, -1.0, 1.0))
+            weights = _evaluate_gaussian_distribution(x, self.sigma_angle)
+            mask = x > (3.0*self.sigma_angle)
+            weights[mask] = 0.0
+            weights /= np.sum(weights, axis=-1)
+        else:
+            weights = np.eye(len(directions))
+        return weights
+
+    def _build_nx_filter(self):
+        directions = self.sphere.vertices.astype(np.float32)
+        half_width = int(round(3*self.sigma_spatial))
+        filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
+
+        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]
+
+        if self.disable_spatial:
+            w_spatial = np.ones(filter_shape)
+        else:
+            w_spatial = _evaluate_gaussian_distribution(distances,
+                                                        self.sigma_spatial)
+
+        cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
+        theta = np.arccos(cos_theta)
+        theta[half_width, half_width, half_width] = 0.0
+
+        if self.disable_align:
+            w_align = np.ones(np.append(filter_shape, (len(directions),)))
+        else:
+            w_align = _evaluate_gaussian_distribution(theta, self.sigma_align)
+
+        w = w_spatial[..., None] * w_align
+
+        if self.j_invariance:
+            w[half_width, half_width, half_width] = 0.0
+
+        # normalize
+        w /= np.sum(w, axis=(0,1,2), keepdims=True)
+
+        return w
+
+    def _get_sf_range(self, sh_data):
+        sf_max = 0.0
+        sf_min = np.inf
+        for sph_id in range(len(self.sphere.vertices)):
+            sf = np.dot(sh_data, self.B[:, sph_id])
+            sf[sf < 0.0] = 0.0
+            sf_max = max(sf_max, np.max(sf))
+            sf_min = min(sf_min, np.min(sf))
+        return sf_max - sf_min
+
+    def _call_opencl(self, sh_data, patch_size):
         uv_weights_offsets =\
             np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter, axis=-1)))
         v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
@@ -136,206 +234,64 @@ def __call__(self, sh_data, patch_size=40):
 
         return out_sh
 
-    def _build_uv_filter(directions, sigma_angle):
-        dot = directions.dot(directions.T)
-        x = np.arccos(np.clip(dot, -1.0, 1.0))
-        weights = _evaluate_gaussian_distribution(sigma_angle, x)
-
-        mask = x > (3.0*sigma_angle)
-        weights[mask] = 0.0
-        weights /= np.sum(weights, axis=-1)
-        return weights
-
-    def _build_nx_filter(directions, sigma_spatial, sigma_align,
-                        disable_spatial=False, disable_align=False,
-                        j_invariance=False):
-        directions = directions.astype(np.float32)
-        half_width = int(round(3*sigma_spatial))
-        filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
-
-        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]
-
-        if disable_spatial:
-            w_spatial = np.ones(filter_shape)
-        else:
-            w_spatial = _evaluate_gaussian_distribution(sigma_spatial, distances)
-
-        cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
-        theta = np.arccos(cos_theta)
-        theta[half_width, half_width, half_width] = 0.0
-
-        if disable_align:
-            w_align = np.ones(np.append(filter_shape, (len(directions),)))
-        else:
-            w_align = _evaluate_gaussian_distribution(sigma_align, theta)
-
-        w = w_spatial[..., None] * w_align
-
-        if j_invariance:
-            w[half_width, half_width, half_width] = 0.0
-
-        # normalize
-        w /= np.sum(w, axis=(0,1,2), keepdims=True)
-
-        return w
-
-
-def angle_aware_bilateral_filtering_cpu(in_sh, sh_order=8,
-                                        sh_basis='descoteaux07',
-                                        in_full_basis=False,
-                                        sphere_str='repulsion724',
-                                        sigma_spatial=1.0,
-                                        sigma_angular=1.0,
-                                        sigma_range=0.5):
-    """
-    Angle-aware bilateral filtering on the CPU
-    (optionally using multiple threads).
-
-    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.
-    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.
-
-    Returns
-    -------
-    out_sh: ndarray (x, y, z, ncoeffs)
-        Output SH coefficient array in full SH basis.
-    """
-    # 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=sh_order, basis_type=sh_basis,
-                        return_inv=False, full_basis=in_full_basis)
-
-    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=sh_order, basis_type=sh_basis,
-                               full_basis=True)
-    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
-
-
-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.
+    def _call_purepython(self, sh_data, sigma_range):
+        nb_sf = len(self.sphere.vertices)
+        mean_sf = np.zeros(sh_data.shape[:-1] + (nb_sf,))
 
-    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)
+        # Apply filter to each sphere vertice
+        for u_sph_id in range(nb_sf):
+            logging.info('Processing sphere direction: {}'
+                         .format(u_sph_id))
+            mean_sf[..., u_sph_id] = self._correlate(sh_data, sigma_range,
+                                                     u_sph_id)
 
+        out_sh = np.array([np.dot(i, self.B_inv) for i in mean_sf],
+                          dtype=sh_data.dtype)
+        return out_sh
 
-def _get_window_directions(shape):
-    """
-    Get directions from center voxel to all neighbours
-    for a window of given shape.
+    def _correlate(self, sh_data, sigma_range, u_index):
+        v_indices = np.flatnonzero(self.uv_filter[u_index])
+        nx_filter = self.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(sh_data.shape[:3])
+        sh_data = np.pad(sh_data, ((half_w, half_w),
+                                   (half_h, half_h),
+                                   (half_d, half_d),
+                                   (0, 0)))
 
-    Parameters
-    ----------
-    shape: tuple
-        Dimensions of the window.
+        sf_u = np.dot(sh_data, self.B[:, u_index])
+        sf_v  = np.dot(sh_data, self.B[:, v_indices])
+        uv_filter = self.uv_filter[u_index, v_indices]
 
-    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
+        _get_range = _evaluate_gaussian_distribution\
+            if not self.disable_range else lambda x, _ : np.ones_like(x)
 
+        for ii in range(out_im.shape[0]):
+            for jj in range(out_im.shape[1]):
+                for kk in range(out_im.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)
 
-def _correlate_spatial(image_u, h_filter, sigma_range):
-    """
-    Implementation of the correlation operation for anisotropic filtering.
+                    # 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)))
+                    # print(res_filter.shape)
+                    out_im[ii, jj, kk] = np.sum(
+                        sf_v[ii:ii+h_w, jj:jj+h_h, kk:kk+h_d] * res_filter)
+                    out_im[ii, jj, kk] /= np.sum(res_filter)
 
-    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.
+        return out_im
 
-    Returns
-    -------
-    out_im: ndarray (X, Y, Z)
-        Filtered SF image.
-    """
-    h_w, h_h, h_d = h_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),
-                               (half_h, half_h),
-                               (half_d, half_d)))
-
-    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
-
-                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)
-
-    return out_im
+    def __call__(self, sh_data, patch_size=40):
+        sigma_range = self.rel_sigma_range * self._get_sf_range(sh_data)
+        if self.use_opencl:
+            self._prepare_opencl(sigma_range)
+            return self._call_opencl(sh_data, patch_size)
+        else:
+            return self._call_purepython(sh_data, sigma_range)
 
 
 def cosine_filtering(in_sh, sh_order=8, sh_basis='descoteaux07',
@@ -374,7 +330,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))
@@ -409,7 +365,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.
 
@@ -453,3 +409,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
\ No newline at end of file
diff --git a/scilpy/gpuparallel/opencl_utils.py b/scilpy/gpuparallel/opencl_utils.py
index 691fd6824..75e149fc7 100644
--- a/scilpy/gpuparallel/opencl_utils.py
+++ b/scilpy/gpuparallel/opencl_utils.py
@@ -223,7 +223,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/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index ffda00ddf..6d64a86c2 100644
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -27,8 +27,7 @@
 from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg,
                              assert_inputs_exist, add_sh_basis_args,
                              assert_outputs_exist)
-from scilpy.denoise.asym_filtering import (cosine_filtering,
-                                           angle_aware_bilateral_filtering)
+from scilpy.denoise.asym_filtering import cosine_filtering, AsymmetricFilter
 
 
 EPILOG = """
@@ -74,20 +73,37 @@ def _build_arg_parser():
                                    ' [%(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]')
+    trilateral_group.add_argument('--sigma_align', default=0.8, type=float,
+                                  help='Standard deviation for alignment '
+                                       'filter. [%(default)s]')
+    trilateral_group.add_argument('--sigma_range', default=0.2, type=float,
+                                  help='Standard deviation for range filter '
+                                       '*relative to SF range of image*. '
+                                       '[%(default)s]')
+    trilateral_group.add_argument('--sigma_angle', default=0.1, type=float,
+                                  help='Standard deviation for angular filter.' 
+                                       '[%(default)s]')
+    trilateral_group.add_argument('--disable_spatial', action='store_true',
+                                  help='Disable spatial filtering.')
+    trilateral_group.add_argument('--disable_align', action='store_true',
+                                  help='Disable alignment filtering.')
+    trilateral_group.add_argument('--disable_angle', action='store_true',
+                                  help='Disable angle filtering, i.e. do not '
+                                       'include\nneighbour sphere directions in'
+                                       ' filtering.')
+    trilateral_group.add_argument('--disable_range', action='store_true',
+                                  help='Disable range filtering.')
 
     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]')
 
-    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\n(requires pyopencl'
+                        ' and a working OpenCL implementation).')
 
     add_verbose_arg(p)
     add_overwrite_arg(p)
@@ -101,7 +117,10 @@ def main():
     if args.verbose:
         logging.getLogger().setLevel(logging.INFO)
 
-    if args.use_gpu and args.method == 'cosine':
+    if args.device == 'gpu' and not args.use_opencl:
+        parser.error('Option --use_opencl is required for device \'gpu\'.')
+
+    if args.use_opencl and args.method == 'cosine':
         parser.error('Option --use_gpu is not supported for cosine filtering.')
 
     outputs = [args.out_sh]
@@ -119,15 +138,21 @@ 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=args.sh_basis,
-            in_full_basis=full_basis,
+        asym_filter = AsymmetricFilter(
+            sh_order=sh_order, sh_basis=args.sh_basis,
+            legacy=True, 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_align=args.sigma_align,
+            sigma_angle=args.sigma_angle,
+            rel_sigma_range=args.sigma_range,
+            disable_spatial=args.disable_spatial,
+            disable_align=args.disable_align,
+            disable_range=args.disable_range,
+            disable_angle=args.disable_angle,
+            device_type=args.device,
+            use_opencl=args.use_opencl)
+        asym_sh = asym_filter(data)
     else:  # args.method == 'cosine'
         asym_sh = cosine_filtering(
             data, sh_order=sh_order,

From 0563a14880f41dc1b4d147dbfb67510c480e8ce2 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 10:17:24 -0500
Subject: [PATCH 04/19] Fix tracking code

---
 scilpy/gpuparallel/opencl_utils.py | 67 ++++++++++++++++++++++--------
 scilpy/tracking/tracker.py         | 34 ++++++++-------
 2 files changed, 68 insertions(+), 33 deletions(-)

diff --git a/scilpy/gpuparallel/opencl_utils.py b/scilpy/gpuparallel/opencl_utils.py
index 75e149fc7..a014bb9cd 100644
--- a/scilpy/gpuparallel/opencl_utils.py
+++ b/scilpy/gpuparallel/opencl_utils.py
@@ -18,19 +18,23 @@ def cl_device_type(device_type_str):
 
 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, device_type='gpu'):
         if not have_opencl:
@@ -92,19 +96,15 @@ def add_input_buffer(self, key, arr=None, 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
@@ -127,6 +127,20 @@ def add_input_buffer(self, key, arr=None, dtype=np.float32):
         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]
@@ -139,17 +153,20 @@ def update_input_buffer(self, key, arr, 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!')
@@ -163,6 +180,20 @@ def add_output_buffer(self, key, shape=None, dtype=np.float32):
         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]
diff --git a/scilpy/tracking/tracker.py b/scilpy/tracking/tracker.py
index d81a7c902..afbac2c06 100644
--- a/scilpy/tracking/tracker.py
+++ b/scilpy/tracking/tracker.py
@@ -566,8 +566,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))
 
@@ -591,25 +589,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)
-        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:
@@ -621,14 +623,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))

From b4a821d5af44d81ae0b30f1566d6336c4606d1bf Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 11:05:43 -0500
Subject: [PATCH 05/19] Better argparsing

---
 scilpy/denoise/asym_filtering.py | 86 ++++++++++++++++----------------
 scripts/scil_sh_to_aodf.py       | 73 ++++++++++++---------------
 2 files changed, 75 insertions(+), 84 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index 8c2a73127..db8ab1f02 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -27,44 +27,43 @@ class AsymmetricFilter():
     sphere_str: str
         Name of the DIPY sphere to use for SH to SF projection.
     sigma_spatial: float
-        Standard deviation of spatial filter.
-    sigma_align: float
-        Standard deviation of alignment filter.
-    sigma_angle: float
-        Standard deviation of the angle filter.
-    rel_sigma_range: float
+        Standard deviation of spatial filter. Also controls the
+        filtering window size (6*sigma_spatial + 1).
+    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.
+        range of SF amplitudes. `None` disables range filtering.
     disable_spatial: bool, optional
-        Disable spatial filtering. Replaces the
-        gaussian weights by a constant value.
-    disable_align: bool, optional
-        Disable alignment filtering.
-    disable_angle: bool, optional
-        Disable angle filtering, which considers neighbour *sphere
-        directions* in the estimation of the a-ODF.
-    disable_range: bool, optional
-        Disable range filtering.
+        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.
+    j_invariance: 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.
     """
-    def __init__(self, sh_order, sh_basis, legacy, full_basis,
-                 sphere_str, sigma_spatial, sigma_align,
-                 sigma_angle, rel_sigma_range, disable_spatial=False,
-                 disable_align=False, disable_angle=False,
-                 disable_range=False, j_invariance=False,
-                 device_type='gpu', use_opencl=True):
+    def __init__(self, sh_order, sh_basis, legacy, full_basis, sphere_str,
+                 sigma_spatial=1.0, sigma_align=0.8, sigma_angle=None,
+                 rel_sigma_range=0.2, disable_spatial=False,
+                 j_invariance=False, device_type='gpu',
+                 use_opencl=True):
         self.sh_order = sh_order
         self.legacy = legacy
         self.basis = sh_basis
         self.full_basis = full_basis
         self.sphere = get_sphere(sphere_str)
-        self.sigma_spatial = sigma_spatial
-        self.sigma_align = sigma_align
-        self.sigma_angle = sigma_angle
         self.rel_sigma_range = rel_sigma_range
-        self.disable_range = disable_range
-        self.disable_angle = disable_angle
+        self.disable_range = rel_sigma_range is None
+        self.disable_angle = sigma_angle is None
         self.disable_spatial = disable_spatial
-        self.disable_align = disable_align
+        self.disable_align = sigma_align is None
         self.device_type = device_type
         self.use_opencl = use_opencl
         self.j_invariance = j_invariance
@@ -80,8 +79,8 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis,
                              'to use device \'gpu\'.')
 
         # build filters
-        self.uv_filter = self._build_uv_filter()
-        self.nx_filter = self._build_nx_filter()
+        self.uv_filter = self._build_uv_filter(sigma_angle)
+        self.nx_filter = self._build_nx_filter(sigma_spatial, sigma_align)
 
         self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
                                  self.basis, self.full_basis,
@@ -108,22 +107,22 @@ def _prepare_opencl(self, sigma_range):
                                                    else 'false')
         self.cl_manager = CLManager(self.cl_kernel, self.device_type)
 
-    def _build_uv_filter(self):
+    def _build_uv_filter(self, sigma_angle):
         directions = self.sphere.vertices
         if not self.disable_angle:
             dot = directions.dot(directions.T)
             x = np.arccos(np.clip(dot, -1.0, 1.0))
-            weights = _evaluate_gaussian_distribution(x, self.sigma_angle)
-            mask = x > (3.0*self.sigma_angle)
+            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 _build_nx_filter(self):
+    def _build_nx_filter(self, sigma_spatial, sigma_align):
         directions = self.sphere.vertices.astype(np.float32)
-        half_width = int(round(3*self.sigma_spatial))
+        half_width = int(round(3*sigma_spatial))
         filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
 
         grid_directions = _get_window_directions(filter_shape).astype(np.float32)
@@ -134,17 +133,15 @@ def _build_nx_filter(self):
         if self.disable_spatial:
             w_spatial = np.ones(filter_shape)
         else:
-            w_spatial = _evaluate_gaussian_distribution(distances,
-                                                        self.sigma_spatial)
-
-        cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
-        theta = np.arccos(cos_theta)
-        theta[half_width, half_width, half_width] = 0.0
+            w_spatial = _evaluate_gaussian_distribution(distances, sigma_spatial)
 
         if self.disable_align:
             w_align = np.ones(np.append(filter_shape, (len(directions),)))
         else:
-            w_align = _evaluate_gaussian_distribution(theta, self.sigma_align)
+            cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
+            theta = np.arccos(cos_theta)
+            theta[half_width, half_width, half_width] = 0.0
+            w_align = _evaluate_gaussian_distribution(theta, sigma_align)
 
         w = w_spatial[..., None] * w_align
 
@@ -240,8 +237,9 @@ def _call_purepython(self, sh_data, sigma_range):
 
         # Apply filter to each sphere vertice
         for u_sph_id in range(nb_sf):
-            logging.info('Processing sphere direction: {}'
-                         .format(u_sph_id))
+            if u_sph_id % 20 == 0:
+                logging.info('Processing direction: {}/{}'
+                             .format(u_sph_id, nb_sf))
             mean_sf[..., u_sph_id] = self._correlate(sh_data, sigma_range,
                                                      u_sph_id)
 
diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 6d64a86c2..6536a2cfb 100644
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -4,15 +4,14 @@
 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. 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, specify the option --device gpu.
 """
 
 import argparse
@@ -31,8 +30,8 @@
 
 
 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
 
 [2] Poirier et al, 2021, "Investigating the Occurrence of Asymmetric Patterns
     in White Matter Fiber Orientation Distribution Functions", ISMRM 2021
@@ -59,40 +58,35 @@ def _build_arg_parser():
                    help='Sphere used for the SH to SF projection. '
                         '[%(default)s]')
 
-    p.add_argument('--method', default='bilateral',
-                   choices=['bilateral', 'cosine'],
+    p.add_argument('--method', default='unified',
+                   choices=['unified', 'cosine'],
                    help='Method for estimating asymmetric ODFs '
                         '[%(default)s].\nOne of:\n'
-                        '    \'bilateral\': Angle-aware bilateral '
-                        'filtering [1].\n'
-                        '    \'cosine\'  : Cosine-based filtering [2].')
+                        '    \'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_align', default=0.8, type=float,
-                                  help='Standard deviation for alignment '
-                                       'filter. [%(default)s]')
-    trilateral_group.add_argument('--sigma_range', default=0.2, type=float,
-                                  help='Standard deviation for range filter '
-                                       '*relative to SF range of image*. '
-                                       '[%(default)s]')
-    trilateral_group.add_argument('--sigma_angle', default=0.1, type=float,
-                                  help='Standard deviation for angular filter.' 
-                                       '[%(default)s]')
-    trilateral_group.add_argument('--disable_spatial', action='store_true',
-                                  help='Disable spatial filtering.')
-    trilateral_group.add_argument('--disable_align', action='store_true',
-                                  help='Disable alignment filtering.')
-    trilateral_group.add_argument('--disable_angle', action='store_true',
-                                  help='Disable angle filtering, i.e. do not '
-                                       'include\nneighbour sphere directions in'
-                                       ' filtering.')
-    trilateral_group.add_argument('--disable_range', action='store_true',
-                                  help='Disable range filtering.')
+    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 '
+                                    '*relative to SF range of image*. '
+                                    '[%(default)s]')
+    unified_group.add_argument('--sigma_angle', type=float,
+                               help='Standard deviation for angular filter. '
+                                    '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.')
 
     cosine_group = p.add_argument_group('Cosine filter arguments')
     cosine_group.add_argument('--sharpness', default=1.0, type=float,
@@ -137,19 +131,18 @@ def main():
 
     t0 = time.perf_counter()
     logging.info('Filtering SH image.')
-    if args.method == 'bilateral':
+    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
         asym_filter = AsymmetricFilter(
             sh_order=sh_order, sh_basis=args.sh_basis,
             legacy=True, full_basis=full_basis,
             sphere_str=args.sphere,
             sigma_spatial=args.sigma_spatial,
-            sigma_align=args.sigma_align,
+            sigma_align=sigma_align,
             sigma_angle=args.sigma_angle,
-            rel_sigma_range=args.sigma_range,
+            rel_sigma_range=sigma_range,
             disable_spatial=args.disable_spatial,
-            disable_align=args.disable_align,
-            disable_range=args.disable_range,
-            disable_angle=args.disable_angle,
             device_type=args.device,
             use_opencl=args.use_opencl)
         asym_sh = asym_filter(data)

From a550d8e46cbbbbbdb26bb7901855f915176f0c56 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 11:42:19 -0500
Subject: [PATCH 06/19] Documentation everywhere

---
 scilpy/denoise/asym_filtering.py | 149 ++++++++++++++++++++++++++++---
 scripts/scil_sh_to_aodf.py       |   3 +-
 2 files changed, 141 insertions(+), 11 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index db8ab1f02..71623e4f6 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -11,7 +11,7 @@
 
 class AsymmetricFilter():
     """
-    Unified asymmetric filtering as described in Poirier et al, 2024.
+    Unified asymmetric filtering as described in [1].
 
     Parameters
     ----------
@@ -48,6 +48,12 @@ class AsymmetricFilter():
         Device on which the code should run. Choices are cpu or gpu.
     use_opencl: bool, optional
         Use OpenCL for software acceleration.
+
+    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
     """
     def __init__(self, sh_order, sh_basis, legacy, full_basis, sphere_str,
                  sigma_spatial=1.0, sigma_align=0.8, sigma_angle=None,
@@ -95,6 +101,15 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis, sphere_str,
             self.cl_manager = None
 
     def _prepare_opencl(self, sigma_range):
+        """
+        Instantiate OpenCL context manager and compile OpenCL program.
+
+        Parameters
+        ----------
+        sigma_range: float
+            Value for sigma_range. Will be cast to float. Must be provided even
+            when range filtering is disabled for the OpenCL program to compile.
+        """
         self.cl_kernel = CLKernel('filter', 'denoise', 'aodf_filter.cl')
 
         self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
@@ -108,6 +123,21 @@ def _prepare_opencl(self, sigma_range):
         self.cl_manager = CLManager(self.cl_kernel, self.device_type)
 
     def _build_uv_filter(self, sigma_angle):
+        """
+        Build the angle filter.
+
+        Parameters
+        ----------
+        sigma_angle: float
+            Standard deviation of filter. Values at distances greater than
+            sigma_angle are clipped to 0 to reduce computation time.
+            Ignored when self.disable_angle is True.
+
+        Returns
+        -------
+        weights: ndarray
+            Angle filter of shape (N_dirs, N_dirs).
+        """
         directions = self.sphere.vertices
         if not self.disable_angle:
             dot = directions.dot(directions.T)
@@ -121,6 +151,24 @@ def _build_uv_filter(self, sigma_angle):
         return weights
 
     def _build_nx_filter(self, sigma_spatial, sigma_align):
+        """
+        Build the combined spatial and alignment filter.
+
+        Parameters
+        ----------
+        sigma_spatial: float
+            Standard deviation of spatial filter. Also controls the
+            filtering window size (6*sigma_spatial + 1).
+        sigma_align: float
+            Standard deviation of the alignment filter. Ignored when
+            self.disable_align is True.
+
+        Returns
+        -------
+        weights: ndarray
+            Combined spatial + alignment filter of shape (W, H, D, N) where
+            N is the number of sphere directions.
+        """
         directions = self.sphere.vertices.astype(np.float32)
         half_width = int(round(3*sigma_spatial))
         filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
@@ -154,6 +202,19 @@ def _build_nx_filter(self, sigma_spatial, sigma_align):
         return w
 
     def _get_sf_range(self, sh_data):
+        """
+        Get the range of SF amplitudes for input `sh_data`.
+
+        Parameters
+        ----------
+        sh_data: ndarray
+            Spherical harmonics coefficients image.
+
+        Returns
+        -------
+        sf_range: float
+            Range of SF amplitudes.
+        """
         sf_max = 0.0
         sf_min = np.inf
         for sph_id in range(len(self.sphere.vertices)):
@@ -164,6 +225,22 @@ def _get_sf_range(self, sh_data):
         return sf_max - sf_min
 
     def _call_opencl(self, sh_data, patch_size):
+        """
+        Run filtering using opencl.
+
+        Parameters
+        ----------
+        sh_data: ndarray
+            Input SH volume.
+        patch_size: int
+            Data is processed in patches of
+            patch_size x patch_size x patch_size.
+
+        Returns
+        -------
+        out_sh: ndarray
+            Filtered output as SH coefficients.
+        """
         uv_weights_offsets =\
             np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter, axis=-1)))
         v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
@@ -232,6 +309,22 @@ def _call_opencl(self, sh_data, patch_size):
         return out_sh
 
     def _call_purepython(self, sh_data, sigma_range):
+        """
+        Run filtering using pure python implementation.
+
+        Parameters
+        ----------
+        sh_data: ndarray
+            Input SH data.
+        sigma_range: float
+            Standard deviation of range filter. Ignored when
+            self.disable_range is True.
+
+        Returns
+        -------
+        out_sh: ndarray
+            Filtered output as SH coefficients.
+        """
         nb_sf = len(self.sphere.vertices)
         mean_sf = np.zeros(sh_data.shape[:-1] + (nb_sf,))
 
@@ -248,11 +341,30 @@ def _call_purepython(self, sh_data, sigma_range):
         return out_sh
 
     def _correlate(self, sh_data, sigma_range, u_index):
+        """
+        Apply the filters to the SH image for the sphere direction
+        described by `u_index`.
+
+        Parameters
+        ----------
+        sh_data: ndarray
+            Input SH coefficients.
+        sigma_range: float
+            Standard deviation of range filter. Ignored when
+            self.disable_range is True.
+        u_index: int
+            Index of the current sphere direction to process.
+
+        Returns
+        -------
+        out_sf: ndarray
+            Output SF amplitudes along the direction described by `u_index`.
+        """
         v_indices = np.flatnonzero(self.uv_filter[u_index])
         nx_filter = self.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(sh_data.shape[:3])
+        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),
@@ -265,9 +377,9 @@ def _correlate(self, sh_data, sigma_range, u_index):
         _get_range = _evaluate_gaussian_distribution\
             if not self.disable_range else lambda x, _ : np.ones_like(x)
 
-        for ii in range(out_im.shape[0]):
-            for jj in range(out_im.shape[1]):
-                for kk in range(out_im.shape[2]):
+        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
@@ -276,15 +388,32 @@ def _correlate(self, sh_data, sigma_range, u_index):
                     # 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)))
-                    # print(res_filter.shape)
-                    out_im[ii, jj, kk] = np.sum(
+                    out_sf[ii, jj, kk] = np.sum(
                         sf_v[ii:ii+h_w, jj:jj+h_h, kk:kk+h_d] * res_filter)
-                    out_im[ii, jj, kk] /= np.sum(res_filter)
+                    out_sf[ii, jj, kk] /= np.sum(res_filter)
 
-        return out_im
+        return out_sf
 
     def __call__(self, sh_data, patch_size=40):
-        sigma_range = self.rel_sigma_range * self._get_sf_range(sh_data)
+        """
+        Run filtering.
+
+        Parameters
+        ----------
+        sh_data: ndarray
+            Input SH coefficients.
+        patch_size: int, optional
+            Data is processed in patches of
+            patch_size x patch_size x patch_size.
+
+        Returns
+        -------
+        out_sh: ndarray
+            Filtered output as SH coefficients.
+        """
+        sigma_range = 0.0
+        if not self.disable_range:
+            sigma_range = self.rel_sigma_range * self._get_sf_range(sh_data)
         if self.use_opencl:
             self._prepare_opencl(sigma_range)
             return self._call_opencl(sh_data, patch_size)
diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 6536a2cfb..57576e71f 100644
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -31,7 +31,8 @@
 
 EPILOG = """
 [1] Poirier and Descoteaux, 2024, "A Unified Filtering Method for Estimating
-    Asymmetric Orientation Distribution Functions", Neuroimage, vol. 287
+    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

From 103c71b70a3f203baa713f460f084c5f552019d9 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 11:47:35 -0500
Subject: [PATCH 07/19] Remove angle_aware_bilateral.cl and module duplicate

---
 scilpy/denoise/angle_aware_bilateral.cl | 131 ---------------
 scilpy/denoise/unified_asym.py          | 206 ------------------------
 2 files changed, 337 deletions(-)
 delete mode 100644 scilpy/denoise/angle_aware_bilateral.cl
 delete mode 100644 scilpy/denoise/unified_asym.py

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/unified_asym.py b/scilpy/denoise/unified_asym.py
deleted file mode 100644
index 7a2f3d17a..000000000
--- a/scilpy/denoise/unified_asym.py
+++ /dev/null
@@ -1,206 +0,0 @@
-# -*- coding: utf-8 -*-
-from dipy.data import get_sphere
-from dipy.reconst.shm import sh_to_sf_matrix
-import numpy as np
-from scilpy.gpuparallel.opencl_utils import CLManager, CLKernel
-from itertools import product as iterprod
-import logging
-
-
-class AsymmetricFilter():
-    def __init__(self, sh_order, sh_basis, legacy, full_basis,
-                 sphere_str, sigma_spatial, sigma_align,
-                 sigma_angle, sigma_range, disable_spatial=False,
-                 disable_align=False, disable_angle=False,
-                 disable_range=False, device_type='gpu',
-                 j_invariance=False):
-        self.sh_order = sh_order
-        self.legacy = legacy
-        self.basis = sh_basis
-        self.full_basis = full_basis
-        self.sphere = get_sphere(sphere_str)
-        self.sigma_spatial = sigma_spatial
-        self.sigma_align = sigma_align
-        self.sigma_angle = sigma_angle
-        self.sigma_range = sigma_range
-        self.disable_range = disable_range
-        self.disable_angle = disable_angle
-        self.device_type = device_type
-
-        # won't need this if disable angle
-        self.uv_filter = _build_uv_filter(self.sphere.vertices,
-                                          self.sigma_angle)
-
-        # sigma still controls the width of the filter
-        self.nx_filter = _build_nx_filter(self.sphere.vertices, sigma_spatial,
-                                          sigma_align, disable_spatial,
-                                          disable_align, j_invariance)
-        logging.info('Filter shape: {}'.format(self.nx_filter.shape[:-1]))
-
-        self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
-                                 self.basis, self.full_basis,
-                                 legacy=self.legacy, return_inv=False)
-        _, self.B_inv = sh_to_sf_matrix(self.sphere, self.sh_order,
-                                        self.basis, True, legacy=self.legacy,
-                                        return_inv=True)
-
-        # initialize gpu
-        self.cl_kernel = None
-        self.cl_manager = None
-        self._prepare_gpu()
-
-    def _prepare_gpu(self):
-        self.cl_kernel = CLKernel('filter', 'filtering', 'aodf_filter.cl')
-
-        self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
-        self.cl_kernel.set_define('SIGMA_RANGE',
-                                  '{}f'.format(self.sigma_range))
-        self.cl_kernel.set_define('N_DIRS', len(self.sphere.vertices))
-        self.cl_kernel.set_define('DISABLE_ANGLE', 'true' if self.disable_angle
-                                                   else 'false')
-        self.cl_kernel.set_define('DISABLE_RANGE', 'true' if self.disable_range
-                                                   else 'false')
-        self.cl_manager = CLManager(self.cl_kernel, self.device_type)
-
-    def __call__(self, sh_data, patch_size=40):
-        uv_weights_offsets =\
-            np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter, axis=-1)))
-        v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
-                            (self.uv_filter.shape[0], 1))[self.uv_filter > 0.0]
-        flat_uv = self.uv_filter[self.uv_filter > 0.0]
-
-        # Prepare GPU buffers
-        self.cl_manager.add_input_buffer("sf_data")  # SF data not initialized
-        self.cl_manager.add_input_buffer("nx_filter", self.nx_filter)
-        self.cl_manager.add_input_buffer("uv_filter", flat_uv)
-        self.cl_manager.add_input_buffer("uv_weights_offsets", uv_weights_offsets)
-        self.cl_manager.add_input_buffer("v_indices", v_indices)
-        self.cl_manager.add_output_buffer("out_sf")  # SF not initialized yet
-
-        win_width = self.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], self.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 + self.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(self.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, self.B)
-            self.cl_manager.update_input_buffer("sf_data", sf_patch)
-            self.cl_manager.update_output_buffer("out_sf", out_shape)
-            out_sf = self.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,
-                                                             self.B_inv)
-
-        return out_sh
-
-
-def _build_uv_filter(directions, sigma_angle):
-    dot = directions.dot(directions.T)
-    x = np.arccos(np.clip(dot, -1.0, 1.0))
-    weights = _evaluate_gaussian(sigma_angle, x)
-
-    mask = x > (3.0*sigma_angle)
-    weights[mask] = 0.0
-    weights /= np.sum(weights, axis=-1)
-    return weights
-
-
-def _build_nx_filter(directions, sigma_spatial, sigma_align,
-                     disable_spatial=False, disable_align=False,
-                     j_invariance=False):
-    directions = directions.astype(np.float32)
-    half_width = int(round(3*sigma_spatial))
-    filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
-
-    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]
-
-    if disable_spatial:
-        w_spatial = np.ones(filter_shape)
-    else:
-        w_spatial = _evaluate_gaussian(sigma_spatial, distances)
-
-    cos_theta = np.clip(grid_directions.dot(directions.T), -1.0, 1.0)
-    theta = np.arccos(cos_theta)
-    theta[half_width, half_width, half_width] = 0.0
-
-    if disable_align:
-        w_align = np.ones(np.append(filter_shape, (len(directions),)))
-    else:
-        w_align = _evaluate_gaussian(sigma_align, theta)
-
-    w = w_spatial[..., None] * w_align
-
-    if j_invariance:
-        w[half_width, half_width, half_width] = 0.0
-
-    # normalize
-    w /= np.sum(w, axis=(0,1,2), keepdims=True)
-
-    return w
-
-
-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
-
-
-def _evaluate_gaussian(sigma, x):
-    # gaussian is not normalized
-    return np.exp(-x**2/(2*sigma**2))

From 161d759233c7c1bbb10bbfe85d5bffb7e31be07a Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 12:45:18 -0500
Subject: [PATCH 08/19] Fix PEP8 for file asym_filtering.py

---
 scilpy/denoise/asym_filtering.py | 33 +++++++++++++++++++-------------
 1 file changed, 20 insertions(+), 13 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index 71623e4f6..ea0d8b086 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -75,8 +75,8 @@ def __init__(self, sh_order, sh_basis, legacy, full_basis, sphere_str,
         self.j_invariance = j_invariance
 
         if device_type not in ['cpu', 'gpu']:
-            raise ValueError('Invalid device type {}. '
-                             'Must be one of \{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.')
@@ -173,15 +173,18 @@ def _build_nx_filter(self, sigma_spatial, sigma_align):
         half_width = int(round(3*sigma_spatial))
         filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
 
-        grid_directions = _get_window_directions(filter_shape).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]
+        grid_directions[distances > 0] =\
+            grid_directions[distances > 0] /\
+            distances[distances > 0][..., None]
 
         if self.disable_spatial:
             w_spatial = np.ones(filter_shape)
         else:
-            w_spatial = _evaluate_gaussian_distribution(distances, sigma_spatial)
+            w_spatial = _evaluate_gaussian_distribution(distances,
+                                                        sigma_spatial)
 
         if self.disable_align:
             w_align = np.ones(np.append(filter_shape, (len(directions),)))
@@ -197,7 +200,7 @@ def _build_nx_filter(self, sigma_spatial, sigma_align):
             w[half_width, half_width, half_width] = 0.0
 
         # normalize
-        w /= np.sum(w, axis=(0,1,2), keepdims=True)
+        w /= np.sum(w, axis=(0, 1, 2), keepdims=True)
 
         return w
 
@@ -242,7 +245,8 @@ def _call_opencl(self, sh_data, patch_size):
             Filtered output as SH coefficients.
         """
         uv_weights_offsets =\
-            np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter, axis=-1)))
+            np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter,
+                                                        axis=-1)))
         v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
                             (self.uv_filter.shape[0], 1))[self.uv_filter > 0.0]
         flat_uv = self.uv_filter[self.uv_filter > 0.0]
@@ -251,7 +255,8 @@ def _call_opencl(self, sh_data, patch_size):
         self.cl_manager.add_input_buffer("sf_data")  # SF data not initialized
         self.cl_manager.add_input_buffer("nx_filter", self.nx_filter)
         self.cl_manager.add_input_buffer("uv_filter", flat_uv)
-        self.cl_manager.add_input_buffer("uv_weights_offsets", uv_weights_offsets)
+        self.cl_manager.add_input_buffer("uv_weights_offsets",
+                                         uv_weights_offsets)
         self.cl_manager.add_input_buffer("v_indices", v_indices)
         self.cl_manager.add_output_buffer("out_sf")  # SF not initialized yet
 
@@ -371,11 +376,11 @@ def _correlate(self, sh_data, sigma_range, u_index):
                                    (0, 0)))
 
         sf_u = np.dot(sh_data, self.B[:, u_index])
-        sf_v  = np.dot(sh_data, self.B[:, v_indices])
+        sf_v = np.dot(sh_data, self.B[:, v_indices])
         uv_filter = self.uv_filter[u_index, v_indices]
 
         _get_range = _evaluate_gaussian_distribution\
-            if not self.disable_range else lambda x, _ : np.ones_like(x)
+            if not self.disable_range else lambda x, _: np.ones_like(x)
 
         for ii in range(out_sf.shape[0]):
             for jj in range(out_sf.shape[1]):
@@ -387,7 +392,9 @@ def _correlate(self, sh_data, sigma_range, u_index):
 
                     # 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)))
+                    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)
@@ -579,4 +586,4 @@ def _get_window_directions(shape):
     grid = np.indices(shape)
     grid = np.moveaxis(grid, 0, -1)
     grid = grid - np.asarray(shape) // 2
-    return grid
\ No newline at end of file
+    return grid

From c1bae95841341e478b730177e95a18778fec492d Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 12:47:25 -0500
Subject: [PATCH 09/19] PEP8 opencl_utils.py

---
 scilpy/gpuparallel/opencl_utils.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/scilpy/gpuparallel/opencl_utils.py b/scilpy/gpuparallel/opencl_utils.py
index a014bb9cd..c6f5dd716 100644
--- a/scilpy/gpuparallel/opencl_utils.py
+++ b/scilpy/gpuparallel/opencl_utils.py
@@ -8,6 +8,7 @@
 from dipy.utils.optpkg import optional_package
 cl, have_opencl, _ = optional_package('pyopencl')
 
+
 def cl_device_type(device_type_str):
     if device_type_str == 'cpu':
         return cl.device_type.CPU
@@ -24,7 +25,7 @@ class CLManager(object):
     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.
@@ -59,7 +60,8 @@ def __init__(self, cl_kernel, device_type='gpu'):
             devices = p.get_devices()
             for d in devices:
                 d_type = d.get_info(cl.device_info.TYPE)
-                if d_type == cl_device_type(device_type) and best_device is None:
+                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:

From 6d2f0b02f770dae56451509c4e83be328f3725fc Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 27 Feb 2024 15:44:59 -0500
Subject: [PATCH 10/19] add unit tests

---
 scilpy/denoise/tests/test_asym_filtering.py |  37 +-
 scilpy/tests/arrays.py                      | 762 ++++++++++----------
 2 files changed, 406 insertions(+), 393 deletions(-)

diff --git a/scilpy/denoise/tests/test_asym_filtering.py b/scilpy/denoise/tests/test_asym_filtering.py
index b626a0ed8..6decfbe5b 100644
--- a/scilpy/denoise/tests/test_asym_filtering.py
+++ b/scilpy/denoise/tests/test_asym_filtering.py
@@ -1,32 +1,45 @@
 import numpy as np
 
 from scilpy.denoise.asym_filtering import \
-        angle_aware_bilateral_filtering_cpu, cosine_filtering
+        AsymmetricFilter, 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'
     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
+    disable_spatial = False
+    j_invariance = 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,
-                                              sphere_str, sigma_spatial,
-                                              sigma_angular, sigma_range)
-
-    assert np.allclose(out, fodf_3x3_order8_descoteaux07_filtered_bilateral)
+    asym_filter = AsymmetricFilter(sh_order, sh_basis, legacy=True,
+                                   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,
+                                   disable_spatial=disable_spatial,
+                                   j_invariance=j_invariance,
+                                   device_type=device_type,
+                                   use_opencl=use_opencl)
+    out = asym_filter(in_sh)
+
+    assert np.allclose(out, fodf_3x3_order8_descoteaux07_filtered_unified)
 
 
 def test_cosine_filtering():
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)

From 5bc4302fe3ae732178940288c8c5243dc1063343 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Wed, 28 Feb 2024 17:02:59 -0500
Subject: [PATCH 11/19] Beautify argparser --help output

---
 scripts/scil_sh_to_aodf.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 72e648f34..c6cc24e74 100755
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -76,11 +76,11 @@ def _build_arg_parser():
                                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 '
+                               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. '
+                               help='Standard deviation for angular filter.\n'
                                     'Disabled by default.')
     unified_group.add_argument('--disable_spatial', action='store_true',
                                help='Disable spatial filtering.')
@@ -91,14 +91,14 @@ def _build_arg_parser():
 
     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('--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\n(requires pyopencl'
-                        ' and a working OpenCL implementation).')
+                   help='Accelerate code using OpenCL (requires pyopencl\n'
+                        'and a working OpenCL implementation).')
 
     add_verbose_arg(p)
     add_overwrite_arg(p)
@@ -133,6 +133,7 @@ def main():
     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
+        # instantiate asymmetric filter
         asym_filter = AsymmetricFilter(
             sh_order=sh_order, sh_basis=args.sh_basis,
             legacy=True, full_basis=full_basis,
@@ -144,7 +145,9 @@ def main():
             disable_spatial=args.disable_spatial,
             device_type=args.device,
             use_opencl=args.use_opencl)
+        # filter the input image
         asym_sh = asym_filter(data)
+
     else:  # args.method == 'cosine'
         asym_sh = cosine_filtering(
             data, sh_order=sh_order,

From 6c831174ccf216d0ec8c828c05a6f5a86d6f3478 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Wed, 28 Feb 2024 17:03:25 -0500
Subject: [PATCH 12/19] Fixed tests (waiting on updated test data)

---
 scripts/tests/test_sh_to_aodf.py | 95 +++++++++-----------------------
 1 file changed, 27 insertions(+), 68 deletions(-)

diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py
index a5d9eb469..9c13016dc 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -16,119 +16,78 @@
 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(data_path, "fodf_descoteaux07_sub.nii.gz"),
+      os.path.join(data_path,
+                   "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_fodf1.nii.gz',
                             '--sphere', 'repulsion100',
-                            '--sigma_angular', '1.0',
+                            '--sigma_align', '0.8',
                             '--sigma_spatial', '1.0',
-                            '--sigma_range', '1.0',
+                            '--sigma_range', '0.2',
+                            '--sigma_angle', '0.06',
+                            '--device', 'cpu',
                             '--sh_basis', 'descoteaux07', '-f',
                             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)
+    test_fodf = nib.load(expected_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):
+@pytest.mark.parametrize("in_fodf,expected_fodf",
+    [[os.path.join(data_path,
+                   "fodf_descoteaux07_sub_unified_asym.nii.gz"),
+      os.path.join(data_path,
+                   "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_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',
+                            '--sigma_range', '0.2',
+                            '--sigma_angle', '0.06',
+                            '--device', 'cpu',
                             '--sh_basis', 'descoteaux07', '-f',
                             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())
 
-
-@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):
-    os.chdir(os.path.expanduser(tmp_dir.name))
-
-    ret = script_runner.run('scil_sh_to_aodf.py',
-                            in_fodf,
-                            'out_fodf3.nii.gz',
-                            '--sphere', 'repulsion100',
-                            '--sigma_angular', '1.0',
-                            '--sigma_spatial', '1.0',
-                            '--sigma_range', '1.0',
-                            '--sh_basis', 'descoteaux07', '-f',
-                            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):
+      os.path.join(data_path, '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',
+                            '--method', 'cosine', '-f',
                             '--sh_basis', 'descoteaux07',
-                            '-f',
                             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())

From e37c4252a101afa018a830de450c5d1e49bcc59b Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 5 Mar 2024 15:08:36 -0500
Subject: [PATCH 13/19] Convert class to functions + set defaults to paper
 defaults

---
 scilpy/denoise/asym_filtering.py            | 799 +++++++++++---------
 scilpy/denoise/tests/test_asym_filtering.py |  33 +-
 scripts/scil_sh_to_aodf.py                  |  28 +-
 scripts/tests/test_sh_to_aodf.py            |   2 +
 4 files changed, 459 insertions(+), 403 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index ea0d8b086..3ebcc748d 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -9,12 +9,18 @@
 from scilpy.gpuparallel.opencl_utils import have_opencl, CLKernel, CLManager
 
 
-class AsymmetricFilter():
+def unified_filtering(sh_data, sh_order, sh_basis, 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):
     """
     Unified asymmetric filtering as described in [1].
 
     Parameters
     ----------
+    sh_data: ndarray
+        SH coefficients image.
     sh_order: int
         Maximum order of spherical harmonics (SH) basis.
     sh_basis: str
@@ -26,9 +32,9 @@ class AsymmetricFilter():
         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
-        Standard deviation of spatial filter. Also controls the
-        filtering window size (6*sigma_spatial + 1).
+    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.
@@ -42,12 +48,17 @@ class AsymmetricFilter():
         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.
-    j_invariance: bool, optional
+    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
     ----------
@@ -55,377 +66,415 @@ class AsymmetricFilter():
         Estimating Asymmetric Orientation Distribution Functions",
         Neuroimage, https://doi.org/10.1016/j.neuroimage.2024.120516
     """
-    def __init__(self, sh_order, sh_basis, legacy, full_basis, sphere_str,
-                 sigma_spatial=1.0, sigma_align=0.8, sigma_angle=None,
-                 rel_sigma_range=0.2, disable_spatial=False,
-                 j_invariance=False, device_type='gpu',
-                 use_opencl=True):
-        self.sh_order = sh_order
-        self.legacy = legacy
-        self.basis = sh_basis
-        self.full_basis = full_basis
-        self.sphere = get_sphere(sphere_str)
-        self.rel_sigma_range = rel_sigma_range
-        self.disable_range = rel_sigma_range is None
-        self.disable_angle = sigma_angle is None
-        self.disable_spatial = disable_spatial
-        self.disable_align = sigma_align is None
-        self.device_type = device_type
-        self.use_opencl = use_opencl
-        self.j_invariance = j_invariance
-
-        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\'.')
-
-        # build filters
-        self.uv_filter = self._build_uv_filter(sigma_angle)
-        self.nx_filter = self._build_nx_filter(sigma_spatial, sigma_align)
-
-        self.B = sh_to_sf_matrix(self.sphere, self.sh_order,
-                                 self.basis, self.full_basis,
-                                 legacy=self.legacy, return_inv=False)
-        _, self.B_inv = sh_to_sf_matrix(self.sphere, self.sh_order,
-                                        self.basis, True, legacy=self.legacy,
-                                        return_inv=True)
-
-        if self.use_opencl:
-            # initialize opencl
-            self.cl_kernel = None
-            self.cl_manager = None
-
-    def _prepare_opencl(self, sigma_range):
-        """
-        Instantiate OpenCL context manager and compile OpenCL program.
-
-        Parameters
-        ----------
-        sigma_range: float
-            Value for sigma_range. Will be cast to float. Must be provided even
-            when range filtering is disabled for the OpenCL program to compile.
-        """
-        self.cl_kernel = CLKernel('filter', 'denoise', 'aodf_filter.cl')
-
-        self.cl_kernel.set_define('WIN_WIDTH', self.nx_filter.shape[0])
-        self.cl_kernel.set_define('SIGMA_RANGE',
-                                  '{}f'.format(sigma_range))
-        self.cl_kernel.set_define('N_DIRS', len(self.sphere.vertices))
-        self.cl_kernel.set_define('DISABLE_ANGLE', 'true' if self.disable_angle
-                                                   else 'false')
-        self.cl_kernel.set_define('DISABLE_RANGE', 'true' if self.disable_range
-                                                   else 'false')
-        self.cl_manager = CLManager(self.cl_kernel, self.device_type)
-
-    def _build_uv_filter(self, sigma_angle):
-        """
-        Build the angle filter.
-
-        Parameters
-        ----------
-        sigma_angle: float
-            Standard deviation of filter. Values at distances greater than
-            sigma_angle are clipped to 0 to reduce computation time.
-            Ignored when self.disable_angle is True.
-
-        Returns
-        -------
-        weights: ndarray
-            Angle filter of shape (N_dirs, N_dirs).
-        """
-        directions = self.sphere.vertices
-        if not self.disable_angle:
-            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 _build_nx_filter(self, sigma_spatial, sigma_align):
-        """
-        Build the combined spatial and alignment filter.
-
-        Parameters
-        ----------
-        sigma_spatial: float
-            Standard deviation of spatial filter. Also controls the
-            filtering window size (6*sigma_spatial + 1).
-        sigma_align: float
-            Standard deviation of the alignment filter. Ignored when
-            self.disable_align is True.
-
-        Returns
-        -------
-        weights: ndarray
-            Combined spatial + alignment filter of shape (W, H, D, N) where
-            N is the number of sphere directions.
-        """
-        directions = self.sphere.vertices.astype(np.float32)
+    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)
+
+    # calculate half-width from sigma_spatial if supplied
+    if sigma_spatial is not None:
         half_width = int(round(3*sigma_spatial))
-        filter_shape = (half_width*2+1, half_width*2+1, half_width*2+1)
-
-        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]
-
-        if self.disable_spatial:
-            w_spatial = np.ones(filter_shape)
-        else:
-            w_spatial = _evaluate_gaussian_distribution(distances,
-                                                        sigma_spatial)
-
-        if self.disable_align:
-            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[half_width, half_width, half_width] = 0.0
-            w_align = _evaluate_gaussian_distribution(theta, sigma_align)
-
-        w = w_spatial[..., None] * w_align
-
-        if self.j_invariance:
-            w[half_width, half_width, half_width] = 0.0
-
-        # normalize
-        w /= np.sum(w, axis=(0, 1, 2), keepdims=True)
-
-        return w
-
-    def _get_sf_range(self, sh_data):
-        """
-        Get the range of SF amplitudes for input `sh_data`.
-
-        Parameters
-        ----------
-        sh_data: ndarray
-            Spherical harmonics coefficients image.
-
-        Returns
-        -------
-        sf_range: float
-            Range of SF amplitudes.
-        """
-        sf_max = 0.0
-        sf_min = np.inf
-        for sph_id in range(len(self.sphere.vertices)):
-            sf = np.dot(sh_data, self.B[:, sph_id])
-            sf[sf < 0.0] = 0.0
-            sf_max = max(sf_max, np.max(sf))
-            sf_min = min(sf_min, np.min(sf))
-        return sf_max - sf_min
-
-    def _call_opencl(self, sh_data, patch_size):
-        """
-        Run filtering using opencl.
-
-        Parameters
-        ----------
-        sh_data: ndarray
-            Input SH volume.
-        patch_size: int
-            Data is processed in patches of
-            patch_size x patch_size x patch_size.
-
-        Returns
-        -------
-        out_sh: ndarray
-            Filtered output as SH coefficients.
-        """
-        uv_weights_offsets =\
-            np.append([0.0], np.cumsum(np.count_nonzero(self.uv_filter,
-                                                        axis=-1)))
-        v_indices = np.tile(np.arange(self.uv_filter.shape[1]),
-                            (self.uv_filter.shape[0], 1))[self.uv_filter > 0.0]
-        flat_uv = self.uv_filter[self.uv_filter > 0.0]
-
-        # Prepare GPU buffers
-        self.cl_manager.add_input_buffer("sf_data")  # SF data not initialized
-        self.cl_manager.add_input_buffer("nx_filter", self.nx_filter)
-        self.cl_manager.add_input_buffer("uv_filter", flat_uv)
-        self.cl_manager.add_input_buffer("uv_weights_offsets",
-                                         uv_weights_offsets)
-        self.cl_manager.add_input_buffer("v_indices", v_indices)
-        self.cl_manager.add_output_buffer("out_sf")  # SF not initialized yet
-
-        win_width = self.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], self.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 + self.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(self.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, self.B)
-            self.cl_manager.update_input_buffer("sf_data", sf_patch)
-            self.cl_manager.update_output_buffer("out_sf", out_shape)
-            out_sf = self.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,
-                                                             self.B_inv)
-
-        return out_sh
-
-    def _call_purepython(self, sh_data, sigma_range):
-        """
-        Run filtering using pure python implementation.
-
-        Parameters
-        ----------
-        sh_data: ndarray
-            Input SH data.
-        sigma_range: float
-            Standard deviation of range filter. Ignored when
-            self.disable_range is True.
-
-        Returns
-        -------
-        out_sh: ndarray
-            Filtered output as SH coefficients.
-        """
-        nb_sf = len(self.sphere.vertices)
-        mean_sf = np.zeros(sh_data.shape[:-1] + (nb_sf,))
-
-        # 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] = self._correlate(sh_data, sigma_range,
-                                                     u_sph_id)
-
-        out_sh = np.array([np.dot(i, self.B_inv) for i in mean_sf],
-                          dtype=sh_data.dtype)
-        return out_sh
-
-    def _correlate(self, sh_data, sigma_range, u_index):
-        """
-        Apply the filters to the SH image for the sphere direction
-        described by `u_index`.
-
-        Parameters
-        ----------
-        sh_data: ndarray
-            Input SH coefficients.
-        sigma_range: float
-            Standard deviation of range filter. Ignored when
-            self.disable_range is True.
-        u_index: int
-            Index of the current sphere direction to process.
-
-        Returns
-        -------
-        out_sf: ndarray
-            Output SF amplitudes along the direction described by `u_index`.
-        """
-        v_indices = np.flatnonzero(self.uv_filter[u_index])
-        nx_filter = self.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_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),
-                                   (0, 0)))
-
-        sf_u = np.dot(sh_data, self.B[:, u_index])
-        sf_v = np.dot(sh_data, self.B[:, v_indices])
-        uv_filter = self.uv_filter[u_index, v_indices]
-
-        _get_range = _evaluate_gaussian_distribution\
-            if not self.disable_range else lambda x, _: np.ones_like(x)
-
-        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)
-
-                    # 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 __call__(self, sh_data, patch_size=40):
-        """
-        Run filtering.
-
-        Parameters
-        ----------
-        sh_data: ndarray
-            Input SH coefficients.
-        patch_size: int, optional
-            Data is processed in patches of
-            patch_size x patch_size x patch_size.
-
-        Returns
-        -------
-        out_sh: ndarray
-            Filtered output as SH coefficients.
-        """
-        sigma_range = 0.0
-        if not self.disable_range:
-            sigma_range = self.rel_sigma_range * self._get_sf_range(sh_data)
-        if self.use_opencl:
-            self._prepare_opencl(sigma_range)
-            return self._call_opencl(sh_data, patch_size)
-        else:
-            return self._call_purepython(sh_data, sigma_range)
+    # 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=legacy, return_inv=False)
+    _, B_inv = sh_to_sf_matrix(sphere, sh_order, sh_basis, True,
+                               legacy=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:
+        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 _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):
+    """
+    Instantiate OpenCL context manager and compile OpenCL program.
+
+    Parameters
+    ----------
+    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
+    -------
+    cl_manager: CLManager
+        OpenCL manager object.
+    """
+    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
+
+    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):
+    """
+    Build the angle filter, weighted on angle between current direction u
+    and neighbour direction v.
+
+    Parameters
+    ----------
+    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
+    -------
+    weights: ndarray
+        Angle filter of shape (N_dirs, N_dirs).
+    """
+    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 _unified_filter_build_nx(filter_shape, sigma_spatial, sigma_align,
+                             sphere, exclude_center):
+    """
+    Build the combined spatial and alignment filter.
+
+    Parameters
+    ----------
+    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
+    -------
+    weights: ndarray
+        Combined spatial + alignment filter of shape (W, H, D, N) where
+        N is the number of sphere directions.
+    """
+    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]
+
+    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 the range of SF amplitudes for input `sh_data`.
+
+    Parameters
+    ----------
+    sh_data: ndarray
+        Spherical harmonics coefficients image.
+    B_mat: ndarray
+        SH to SF projection matrix.
+
+    Returns
+    -------
+    sf_range: float
+        Range of SF amplitudes.
+    """
+    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 _unified_filter_call_opencl(sh_data, nx_filter, uv_filter, cl_manager,
+                                B, B_inv, sphere, patch_size=40):
+    """
+    Run unified filtering for asymmetric ODFs using OpenCL.
+
+    Parameters
+    ----------
+    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
+    -------
+    out_sh: ndarray
+        Filtered output as SH coefficients.
+    """
+    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)
+
+    return out_sh
+
+
+def _unified_filter_call_python(sh_data, nx_filter, uv_filter, sigma_range,
+                                B_mat, B_inv, sphere):
+    """
+    Run filtering using pure python implementation.
+
+    Parameters
+    ----------
+    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
+    -------
+    out_sh: ndarray
+        Filtered output as SH coefficients.
+    """
+    nb_sf = len(sphere.vertices)
+    mean_sf = np.zeros(sh_data.shape[:-1] + (nb_sf,))
+
+    # 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(sh_data, nx_filter, uv_filter, sigma_range, u_index, B_mat):
+    """
+    Apply the filters to the SH image for the sphere direction
+    described by `u_index`.
+
+    Parameters
+    ----------
+    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_sf: ndarray
+        Output SF amplitudes along the direction described by `u_index`.
+    """
+    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_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),
+                               (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]
+
+    _get_range = _evaluate_gaussian_distribution\
+        if sigma_range is not None else lambda x, _: np.ones_like(x)
+
+    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)
+
+                # 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',
diff --git a/scilpy/denoise/tests/test_asym_filtering.py b/scilpy/denoise/tests/test_asym_filtering.py
index 6decfbe5b..8b3e0dd91 100644
--- a/scilpy/denoise/tests/test_asym_filtering.py
+++ b/scilpy/denoise/tests/test_asym_filtering.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from scilpy.denoise.asym_filtering import \
-        AsymmetricFilter, 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,
@@ -20,26 +20,25 @@ def test_unified_asymmetric_filtering():
     sigma_align = 0.8
     sigma_range = 0.2
     sigma_angle = 0.06
-    disable_spatial = False
-    j_invariance = False
+    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])
-    asym_filter = AsymmetricFilter(sh_order, sh_basis, legacy=True,
-                                   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,
-                                   disable_spatial=disable_spatial,
-                                   j_invariance=j_invariance,
-                                   device_type=device_type,
-                                   use_opencl=use_opencl)
-    out = asym_filter(in_sh)
-
-    assert np.allclose(out, fodf_3x3_order8_descoteaux07_filtered_unified)
+    asym_sh = unified_filtering(in_sh, sh_order, sh_basis, legacy=True,
+                                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():
diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index c6cc24e74..9e3e1eb07 100755
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -26,7 +26,7 @@
 from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg,
                              assert_inputs_exist, add_sh_basis_args,
                              assert_outputs_exist)
-from scilpy.denoise.asym_filtering import cosine_filtering, AsymmetricFilter
+from scilpy.denoise.asym_filtering import cosine_filtering, unified_filtering
 
 
 EPILOG = """
@@ -54,7 +54,7 @@ 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]')
@@ -80,14 +80,19 @@ def _build_arg_parser():
                                     '*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.')
+                               help='Standard deviation for angular filter '
+                                    '(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,
@@ -99,6 +104,8 @@ def _build_arg_parser():
     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)
@@ -133,21 +140,20 @@ def main():
     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
-        # instantiate asymmetric filter
-        asym_filter = AsymmetricFilter(
+        sigma_spatial = None if args.disable_spatial else args.sigma_spatial
+
+        asym_sh = unified_filtering(data,
             sh_order=sh_order, sh_basis=args.sh_basis,
             legacy=True, full_basis=full_basis,
             sphere_str=args.sphere,
-            sigma_spatial=args.sigma_spatial,
+            sigma_spatial=sigma_spatial,
             sigma_align=sigma_align,
             sigma_angle=args.sigma_angle,
             rel_sigma_range=sigma_range,
-            disable_spatial=args.disable_spatial,
+            win_hwidth=args.win_hwidth,
+            exclude_center=not(args.include_center),
             device_type=args.device,
             use_opencl=args.use_opencl)
-        # filter the input image
-        asym_sh = asym_filter(data)
-
     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 9c13016dc..f5344a076 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -37,6 +37,7 @@ def test_asym_basis_output(script_runner, in_fodf, expected_fodf):
                             '--sigma_angle', '0.06',
                             '--device', 'cpu',
                             '--sh_basis', 'descoteaux07', '-f',
+                            '--include_center',
                             print_result=True, shell=True)
 
     assert ret.success
@@ -63,6 +64,7 @@ def test_asym_input(script_runner, in_fodf, expected_fodf):
                             '--sigma_angle', '0.06',
                             '--device', 'cpu',
                             '--sh_basis', 'descoteaux07', '-f',
+                            '--include_center',
                             print_result=True, shell=True)
 
     assert ret.success

From d2a92cac3cb1fd93e60a260b1802f6117ab30c1b Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Tue, 5 Mar 2024 15:23:44 -0500
Subject: [PATCH 14/19] Repair tests

---
 scilpy/denoise/tests/test_asym_filtering.py | 6 ++++--
 scripts/scil_sh_to_aodf.py                  | 2 +-
 scripts/tests/test_sh_to_aodf.py            | 6 +++---
 3 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/scilpy/denoise/tests/test_asym_filtering.py b/scilpy/denoise/tests/test_asym_filtering.py
index 368fd9004..80b67d486 100644
--- a/scilpy/denoise/tests/test_asym_filtering.py
+++ b/scilpy/denoise/tests/test_asym_filtering.py
@@ -15,6 +15,7 @@ def test_unified_asymmetric_filtering():
     """
     in_sh = fodf_3x3_order8_descoteaux07
     sh_basis = 'descoteaux07'
+    legacy=True
     sphere_str = 'repulsion100'
     sigma_spatial = 1.0
     sigma_align = 0.8
@@ -26,7 +27,7 @@ def test_unified_asymmetric_filtering():
     use_opencl = False
 
     sh_order, full_basis = get_sh_order_and_fullness(in_sh.shape[-1])
-    asym_sh = unified_filtering(in_sh, sh_order, sh_basis, legacy=True,
+    asym_sh = unified_filtering(in_sh, sh_order, sh_basis, legacy=legacy,
                                 full_basis=full_basis,
                                 sphere_str=sphere_str,
                                 sigma_spatial=sigma_spatial,
@@ -47,12 +48,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/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 23fbd10ac..41ca10546 100755
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -144,7 +144,7 @@ def main():
         sigma_spatial = None if args.disable_spatial else args.sigma_spatial
 
         asym_sh = unified_filtering(data,
-            sh_order=sh_order, sh_basis=args.sh_basis,
+            sh_order=sh_order, sh_basis=sh_basis,
             legacy=is_legacy, full_basis=full_basis,
             sphere_str=args.sphere,
             sigma_spatial=sigma_spatial,
diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py
index f5344a076..4ae719dc1 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -36,7 +36,7 @@ def test_asym_basis_output(script_runner, in_fodf, expected_fodf):
                             '--sigma_range', '0.2',
                             '--sigma_angle', '0.06',
                             '--device', 'cpu',
-                            '--sh_basis', 'descoteaux07', '-f',
+                            '--sh_basis', 'descoteaux07_legacy', '-f',
                             '--include_center',
                             print_result=True, shell=True)
 
@@ -63,7 +63,7 @@ def test_asym_input(script_runner, in_fodf, expected_fodf):
                             '--sigma_range', '0.2',
                             '--sigma_angle', '0.06',
                             '--device', 'cpu',
-                            '--sh_basis', 'descoteaux07', '-f',
+                            '--sh_basis', 'descoteaux07_legacy', '-f',
                             '--include_center',
                             print_result=True, shell=True)
 
@@ -84,7 +84,7 @@ def test_cosine_method(script_runner, in_fodf, out_fodf):
                             in_fodf, 'out_fodf1.nii.gz',
                             '--sphere', 'repulsion100',
                             '--method', 'cosine', '-f',
-                            '--sh_basis', 'descoteaux07',
+                            '--sh_basis', 'descoteaux07_legacy',
                             print_result=True, shell=True)
 
     assert ret.success

From efa1d6f88344839ef764f4281fc5045adef06de4 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Wed, 6 Mar 2024 13:31:35 -0500
Subject: [PATCH 15/19] Use DVC for testing

---
 scripts/tests/test_sh_to_aodf.py | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py
index 73b2473a8..3e69927b5 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -8,12 +8,10 @@
 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
 
 # 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()
 
 
@@ -23,8 +21,8 @@ def test_help_option(script_runner):
 
 
 @pytest.mark.parametrize("in_fodf,expected_fodf",
-    [[os.path.join(data_path, "fodf_descoteaux07_sub.nii.gz"),
-      os.path.join(data_path,
+    [[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))
@@ -49,9 +47,9 @@ def test_asym_basis_output(script_runner, in_fodf, expected_fodf):
 
 
 @pytest.mark.parametrize("in_fodf,expected_fodf",
-    [[os.path.join(data_path,
+    [[os.path.join(test_data_root,
                    "fodf_descoteaux07_sub_unified_asym.nii.gz"),
-      os.path.join(data_path,
+      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))
@@ -76,8 +74,9 @@ def test_asym_input(script_runner, in_fodf, expected_fodf):
 
 
 @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_cosine_asym.nii.gz')]])
+    [[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))
 

From ce84cf35f031b625caea864cca3f5b64da6a45c8 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Wed, 6 Mar 2024 13:43:40 -0500
Subject: [PATCH 16/19] Make PEP8 happy

---
 scilpy/denoise/tests/test_asym_filtering.py |  4 ++--
 scripts/scil_sh_to_aodf.py                  |  6 ++---
 scripts/tests/test_sh_to_aodf.py            | 26 ++++++++++-----------
 3 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/scilpy/denoise/tests/test_asym_filtering.py b/scilpy/denoise/tests/test_asym_filtering.py
index 80b67d486..543f6fed9 100644
--- a/scilpy/denoise/tests/test_asym_filtering.py
+++ b/scilpy/denoise/tests/test_asym_filtering.py
@@ -15,13 +15,13 @@ def test_unified_asymmetric_filtering():
     """
     in_sh = fodf_3x3_order8_descoteaux07
     sh_basis = 'descoteaux07'
-    legacy=True
+    legacy = True
     sphere_str = 'repulsion100'
     sigma_spatial = 1.0
     sigma_align = 0.8
     sigma_range = 0.2
     sigma_angle = 0.06
-    win_hwidth=3
+    win_hwidth = 3
     exclude_center = False
     device_type = 'cpu'
     use_opencl = False
diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 41ca10546..10655f447 100755
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -143,8 +143,8 @@ def main():
         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,
+        asym_sh = unified_filtering(
+            data, sh_order=sh_order, sh_basis=sh_basis,
             legacy=is_legacy, full_basis=full_basis,
             sphere_str=args.sphere,
             sigma_spatial=sigma_spatial,
@@ -152,7 +152,7 @@ def main():
             sigma_angle=args.sigma_angle,
             rel_sigma_range=sigma_range,
             win_hwidth=args.win_hwidth,
-            exclude_center=not(args.include_center),
+            exclude_center=not args.include_center,
             device_type=args.device,
             use_opencl=args.use_opencl)
     else:  # args.method == 'cosine'
diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py
index 3e69927b5..447cbcd1b 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -20,10 +20,10 @@ def test_help_option(script_runner):
     assert 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")]])
+@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))
 
@@ -46,11 +46,11 @@ def test_asym_basis_output(script_runner, in_fodf, expected_fodf):
     assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata())
 
 
-@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")]])
+@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))
 
@@ -73,10 +73,10 @@ def test_asym_input(script_runner, in_fodf, expected_fodf):
     assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata())
 
 
-@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')]])
+@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))
 

From e999a6e7d449d713b71e84770f646f50c9497345 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Mon, 11 Mar 2024 16:39:34 -0400
Subject: [PATCH 17/19] Answer review comments

---
 scilpy/denoise/asym_filtering.py            | 8 ++++----
 scilpy/denoise/tests/test_asym_filtering.py | 3 ++-
 scripts/scil_sh_to_aodf.py                  | 4 ++--
 3 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index cd396c158..89a6f05ad 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -9,7 +9,7 @@
 from scilpy.gpuparallel.opencl_utils import have_opencl, CLKernel, CLManager
 
 
-def unified_filtering(sh_data, sh_order, sh_basis, legacy, full_basis,
+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,
@@ -26,7 +26,7 @@ def unified_filtering(sh_data, sh_order, sh_basis, legacy, full_basis,
     sh_basis: str
         SH basis definition used for input and output SH image.
         One of 'descoteaux07' or 'tournier07'.
-    legacy: bool
+    is_legacy: bool
         Whether the legacy SH basis definition should be used.
     full_basis: bool
         Whether the input SH basis is full or not.
@@ -96,9 +96,9 @@ def unified_filtering(sh_data, sh_order, sh_basis, legacy, full_basis,
                                          sigma_align, sphere, exclude_center)
 
     B = sh_to_sf_matrix(sphere, sh_order, sh_basis, full_basis,
-                        legacy=legacy, return_inv=False)
+                        legacy=is_legacy, return_inv=False)
     _, B_inv = sh_to_sf_matrix(sphere, sh_order, sh_basis, True,
-                               legacy=legacy, return_inv=True)
+                               legacy=is_legacy, return_inv=True)
 
     # compute "real" sigma_range scaled by sf amplitudes
     # if rel_sigma_range is supplied
diff --git a/scilpy/denoise/tests/test_asym_filtering.py b/scilpy/denoise/tests/test_asym_filtering.py
index 543f6fed9..b1405548f 100644
--- a/scilpy/denoise/tests/test_asym_filtering.py
+++ b/scilpy/denoise/tests/test_asym_filtering.py
@@ -27,7 +27,8 @@ def test_unified_asymmetric_filtering():
     use_opencl = False
 
     sh_order, full_basis = get_sh_order_and_fullness(in_sh.shape[-1])
-    asym_sh = unified_filtering(in_sh, sh_order, sh_basis, legacy=legacy,
+    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,
diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 10655f447..31d4991a6 100755
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -121,7 +121,7 @@ def main():
         parser.error('Option --use_opencl is required for device \'gpu\'.')
 
     if args.use_opencl and args.method == 'cosine':
-        parser.error('Option --use_gpu is not supported for cosine filtering.')
+        parser.error('Option --use_opencl is not supported for cosine filtering.')
 
     outputs = [args.out_sh]
     if args.out_sym:
@@ -145,7 +145,7 @@ def main():
 
         asym_sh = unified_filtering(
             data, sh_order=sh_order, sh_basis=sh_basis,
-            legacy=is_legacy, full_basis=full_basis,
+            is_legacy=is_legacy, full_basis=full_basis,
             sphere_str=args.sphere,
             sigma_spatial=sigma_spatial,
             sigma_align=sigma_align,

From 05a3767238f56f0d7da3ade5a60d0ccaf0304d2e Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Wed, 13 Mar 2024 13:49:20 -0400
Subject: [PATCH 18/19] Answers comments from Emmanuelle

---
 scilpy/denoise/asym_filtering.py | 13 +++++++++++-
 scripts/scil_sh_to_aodf.py       | 25 ++++++++++++++---------
 scripts/tests/test_sh_to_aodf.py | 35 ++++++++++++++++++++++++++++++++
 3 files changed, 62 insertions(+), 11 deletions(-)

diff --git a/scilpy/denoise/asym_filtering.py b/scilpy/denoise/asym_filtering.py
index 89a6f05ad..c0a63b3e4 100644
--- a/scilpy/denoise/asym_filtering.py
+++ b/scilpy/denoise/asym_filtering.py
@@ -80,9 +80,18 @@ def unified_filtering(sh_data, sh_order, sh_basis, is_legacy, full_basis,
 
     sphere = get_sphere(sphere_str)
 
-    # calculate half-width from sigma_spatial if supplied
     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
@@ -104,6 +113,8 @@ def unified_filtering(sh_data, sh_order, sh_basis, is_legacy, full_basis,
     # 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:
diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py
index 31d4991a6..48c4cbb53 100755
--- a/scripts/scil_sh_to_aodf.py
+++ b/scripts/scil_sh_to_aodf.py
@@ -9,9 +9,10 @@
     * Cosine filtering [2] is a simpler implementation using cosine distance
       for assigning weights to neighbours.
 
-Unified filtering can be accelerated using 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, specify the option --device gpu.
+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
@@ -61,10 +62,10 @@ def _build_arg_parser():
 
     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].')
+                   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,
@@ -80,7 +81,7 @@ def _build_arg_parser():
                                     '*relative to SF range of image*. '
                                     '[%(default)s]')
     unified_group.add_argument('--sigma_angle', type=float,
-                               help='Standard deviation for angular filter '
+                               help='Standard deviation for angular filter\n'
                                     '(disabled by default).')
     unified_group.add_argument('--disable_spatial', action='store_true',
                                help='Disable spatial filtering.')
@@ -118,10 +119,14 @@ def main():
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
     if args.device == 'gpu' and not args.use_opencl:
-        parser.error('Option --use_opencl is required for device \'gpu\'.')
+        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.')
+        parser.error('Option --use_opencl is not supported'
+                     ' for cosine filtering.')
 
     outputs = [args.out_sh]
     if args.out_sym:
diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py
index 447cbcd1b..4c0ff26fd 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -9,6 +9,7 @@
 import pytest
 
 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)
 test_data_root = pull_test_case_package("aodf")
@@ -20,6 +21,40 @@ def test_help_option(script_runner):
     assert 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_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_align', '0.8',
+                            '--sigma_spatial', '1.0',
+                            '--sigma_range', '0.2',
+                            '--sigma_angle', '0.06',
+                            '--use_opencl',
+                            '--device', 'gpu',
+                            '--sh_basis', 'descoteaux07_legacy', '-f',
+                            '--include_center',
+                            print_result=True, shell=True)
+
+    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,

From 7555a077cc785650de8939d2412bda0d813005b5 Mon Sep 17 00:00:00 2001
From: CHrlS98 <charles.poirier@usherbrooke.ca>
Date: Wed, 13 Mar 2024 13:50:25 -0400
Subject: [PATCH 19/19] Fix PEP8

---
 scripts/tests/test_sh_to_aodf.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/scripts/tests/test_sh_to_aodf.py b/scripts/tests/test_sh_to_aodf.py
index 4c0ff26fd..af24322a3 100644
--- a/scripts/tests/test_sh_to_aodf.py
+++ b/scripts/tests/test_sh_to_aodf.py
@@ -49,7 +49,9 @@ def test_asym_basis_output_gpu(script_runner, in_fodf, expected_fodf):
         # 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)
+        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