Merge pull request #440 from clEsperanto/hessian_gaussian_derivative

Ridge detector and some other functions

Author: StRigaud

CMakeLists.txt

@@ -1,8 +1,8 @@
 cmake_minimum_required(VERSION 3.20)
 
-project(CLIc VERSION 0.17.1)
+project(CLIc VERSION 0.18.0)
 
-set(kernel_version_tag "3.2.0"  CACHE STRING "clEsperanto kernel version tag") # "3.1.1"
+set(kernel_version_tag "3.3.0"  CACHE STRING "clEsperanto kernel version tag")
 set(eigen_lib_version_tag "3.4.0" CACHE STRING "Eigen library version tag")
 
 # if not set, set the default build type to Release

clic/include/execution.hpp

@@ -51,6 +51,29 @@ execute_separable(const Device::Pointer &      device,
                   const std::array<float, 3> & sigma,
                   const std::array<int, 3> &   radius) -> void;
 
+/**
+ * @brief Execution function for separable kernel operation (with with extra parameters)
+ *
+ * This is currently only used for Gaussian derivative kernels.
+ * to be fused / generalized with the previous function in future
+ *
+ * @param device Device pointer
+ * @param kernel Kernel function name and code
+ * @param src Source array
+ * @param dst Destination array
+ * @param sigma Sigma value for the kernel
+ * @param radius Radius value for the kernel
+ * @param orders Orders for the kernel operation
+ */
+auto
+execute_separable(const Device::Pointer &      device,
+                  const KernelInfo &           kernel,
+                  const Array::Pointer &       src,
+                  const Array::Pointer &       dst,
+                  const std::array<float, 3> & sigma,
+                  const std::array<int, 3> &   radius,
+                  const std::array<int, 3> &   orders) -> void;
+
 /**
  * @brief Execute a kernel using native framework
  * @param device Device pointer

clic/include/tier1.hpp

@@ -745,6 +745,36 @@ gaussian_blur_func(const Device::Pointer & device,
                    float                   sigma_y,
                    float                   sigma_z) -> Array::Pointer;
 
+/**
+ * @name gaussian_derivative
+ * @brief Convolve the image with a gaussian derivate. The filter kernel can have nonisotropic sigma and order.
+ * The implementation is done separable. In case a sigma equals zero, the direction is not filtered.
+ * If all orders are zero, the filtering is equivalent to a Gaussian blur.
+ *
+ * @param device Device to perform the operation on. [const Device::Pointer &]
+ * @param src Input image to process. [const Array::Pointer &]
+ * @param dst Output result image. [Array::Pointer ( = None )]
+ * @param sigma_x Sigma value along the x axis. [float ( = 0 )]
+ * @param sigma_y Sigma value along the y axis. [float ( = 0 )]
+ * @param sigma_z Sigma value along the z axis. [float ( = 0 )]
+ * @param order_x Order of derivation along the x axis. [int ( = 0 )]
+ * @param order_y Order of derivation along the y axis. [int ( = 0 )]
+ * @param order_z Order of derivation along the z axis. [int ( = 0 )]
+ * @return Array::Pointer
+ *
+ * @note 'filter', 'in assistant'
+ */
+auto
+gaussian_derivative_func(const Device::Pointer & device,
+                         const Array::Pointer &  src,
+                         Array::Pointer          dst,
+                         float                   sigma_x,
+                         float                   sigma_y,
+                         float                   sigma_z,
+                         int                     order_x,
+                         int                     order_y,
+                         int                     order_z) -> Array::Pointer;
+
 
 /**
  * @name generate_distance_matrix
@@ -906,10 +936,14 @@ greater_or_equal_constant_func(const Device::Pointer & device,
 /**
  * @name hessian_eigenvalues
  * @brief Computes the eigenvalues of the hessian matrix of a 2d or 3d image. Hessian matrix or 2D images: [Ixx, Ixy]
- * [Ixy, Iyy] Hessian matrix for 3D images: [Ixx, Ixy, Ixz] [Ixy, Iyy, Iyz] [Ixz, Iyz, Izz] Ixx denotes the second
+ * [Ixy, Iyy] Hessian matrix for 3D images: [Ixx, Ixy, Ixz] [Ixy, Iyy, Iyz] [Ixz, Iyz, Izz], Ixx denotes the second
  * derivative in x. Ixx and Iyy are calculated by convolving the image with the 1d kernel [1 2 1]. Ixy is calculated by
- * a convolution with the 2d kernel: [ 0.25 0 0.25] [ 0 0 0] [0.25 0 0.25] Note: This is the only clesperanto function
- * that returns multiple images. This API might be subject to change in the future. Consider using
+ * a convolution with the 2d kernel: [ 0.25 0 0.25] [ 0 0 0] [0.25 0 0.25]
+ *
+ * The function return the list of eigenvalues as images, by decreasing order. The first image is the largest
+ * eigenvalue,
+ *
+ * Note: This function returns multiple images. This API might be subject to change in the future. Consider using
  * small_hessian_eigenvalue() and/or large_hessian_eigenvalue() instead which return only one image.
  *
  * @param device Device to perform the operation on. [const Device::Pointer &]

clic/include/tier2.hpp

@@ -1105,7 +1105,9 @@ top_hat_func(const Device::Pointer & device,
 
 /**
  * @name extended_depth_of_focus_variance_projection
- * @brief Extended depth of focus projection maximizing local pixel intensity variance.
+ * @brief Depth projection using the local variance maxima to determine the best focus plane.
+ * The radius parameter control the local variance calculation, and the sigma apply a
+ * gaussian blur for smoothness of the projection.
  *
  * @param device Device to perform the operation on. [const Device::Pointer &]
  * @param src Input image to process. [const Array::Pointer &]
@@ -1125,6 +1127,51 @@ extended_depth_of_focus_variance_projection_func(const Device::Pointer & device,
                                                  float                   radius_y,
                                                  float                   sigma) -> Array::Pointer;
 
+/**
+ * @name extended_depth_of_focus_sobel_projection
+ * @brief Depth projection using the local sobel gradient magnitude maxima to determine the best focus plane.
+ * The sigma apply a gaussian blur for smoothness of the projection.
+ *
+ * @param device Device to perform the operation on. [const Device::Pointer &]
+ * @param src Input image to process. [const Array::Pointer &]
+ * @param dst Output result image. [Array::Pointer ( = None )]
+ * @param sigma Sigma for Gaussian blur. [float ( = 5 )]
+ * @return Array::Pointer
+ *
+ * @note 'projection'
+ */
+auto
+extended_depth_of_focus_sobel_projection_func(const Device::Pointer & device,
+                                              const Array::Pointer &  src,
+                                              Array::Pointer          dst,
+                                              float                   sigma) -> Array::Pointer;
+
+/**
+ * @name hessian_gaussian_eigenvalues
+ * @brief Determines the Hessian matrix eigenvalues using the gaussian derivative method and returns the small,
+ * middle and large eigenvalue images.
+ *
+ * The function return the list of eigenvalues as images, by decreasing order. The first image is the largest
+ * eigenvalue,
+ *
+ * @param device Device to perform the operation on. [const Device::Pointer &]
+ * @param src Input image to process. [const Array::Pointer &]
+ * @param small_eigenvalue Output result image for the small eigenvalue. [Array::Pointer ( = None )]
+ * @param middle_eigenvalue Output result image for the middle eigenvalue. [Array::Pointer ( = None )]
+ * @param large_eigenvalue Output result image for the large eigenvalue. [Array::Pointer ( = None )]
+ * @param sigma Sigma of the Gaussian kernel. [float ( = 1 )]
+ * @return std::vector<Array::Pointer>
+ *
+ */
+auto
+hessian_gaussian_eigenvalues_func(const Device::Pointer & device,
+                                  const Array::Pointer &  src,
+                                  Array::Pointer          small_eigenvalue,
+                                  Array::Pointer          middle_eigenvalue,
+                                  Array::Pointer          large_eigenvalue,
+                                  float                   sigma) -> std::vector<Array::Pointer>;
+
+
 } // namespace cle::tier2
 
 #endif // __INCLUDE_TIER2_HPP

clic/include/tier3.hpp

@@ -419,6 +419,30 @@ statistics_of_background_and_labelled_pixels_func(const Device::Pointer & device
                                                   Array::Pointer          intensity,
                                                   Array::Pointer          label) -> StatisticsMap;
 
+/**
+ * @name sato_filter
+ * @brief Applies the multi-scale ridge detection Sato filter.
+ * This filter was first introduced by Sato et al. in 1998 (https://doi.org/10.1016/S1361-8415(98)80009-1)
+ *
+ * @param device Device to perform the operation on. [const Device::Pointer &]
+ * @param src Input image to process. [const Array::Pointer &]
+ * @param dst Output result image. [Array::Pointer ( = None )]
+ * @param sigma_minimum Minimum sigma value for the filter. [float ( = 1 )]
+ * @param sigma_maximum Maximum sigma value for the filter. [float ( = 3 )]
+ * @param sigma_step Step size for the sigma values. [float ( = 1 )]
+ * @return Array::Pointer
+ *
+ * @note 'filter', 'in assistant'
+ * @see https://doi.org/10.1016/S1361-8415(98)80009-1
+ */
+auto
+sato_filter_func(const Device::Pointer & device,
+                 const Array::Pointer &  src,
+                 Array::Pointer          dst,
+                 float                   sigma_minimum,
+                 float                   sigma_maximum,
+                 float                   sigma_step) -> Array::Pointer;
+
 } // namespace cle::tier3
 
 #endif // __INCLUDE_TIER3_HPP

clic/include/tier4.hpp

@@ -409,6 +409,20 @@ standard_deviation_intensity_map_func(const Device::Pointer & device,
                                       Array::Pointer          dst) -> Array::Pointer;
 
 
+/**
+ * @name percentile
+ * @brief Computes the percentile value of an image.
+ *
+ * @param device Device to perform the operation on. [const Device::Pointer &]
+ * @param src Input image. [const Array::Pointer &]
+ * @param percentile Percentile to compute. [float ( = 50.0 )]
+ * @return float
+ *
+ */
+auto
+percentile_func(const Device::Pointer & device, const Array::Pointer & src, const float percentile) -> float;
+
+
 } // namespace cle::tier4
 
 #endif // __INCLUDE_TIER4_HPP

clic/src/execution.cpp

@@ -277,6 +277,42 @@ execute(const Device::Pointer & device,
 }
 
 
+auto
+execute_separable(const Device::Pointer &      device,
+                  const KernelInfo &           kernel,
+                  const Array::Pointer &       src,
+                  const Array::Pointer &       dst,
+                  const std::array<float, 3> & sigma,
+                  const std::array<int, 3> &   radius,
+                  const std::array<int, 3> &   orders) -> void
+{
+  Array::check_ptr(src, "Error: 'src' is null. Please ensure the provided parameters before calling this function.");
+  Array::check_ptr(dst, "Error: 'dst' is null. Please ensure the provided parameters before calling this function.");
+
+  const RangeArray global_range = { dst->width(), dst->height(), dst->depth() };
+
+  auto tmp1 = Array::create(dst);
+  auto tmp2 = Array::create(dst);
+
+  auto execute_if_needed = [&](int dim, int idx, auto & input, auto & output) {
+    if (dim > 1 && sigma[idx] > 0)
+    {
+      const ParameterList parameters = { { "src", input },     { "dst", output },   { "dim", idx },
+                                         { "N", radius[idx] }, { "s", sigma[idx] }, { "order", orders[idx] } };
+      execute(device, kernel, parameters, global_range);
+    }
+    else
+    {
+      input->copyTo(output);
+    }
+  };
+
+  execute_if_needed(dst->width(), 0, src, tmp1);
+  execute_if_needed(dst->height(), 1, tmp1, tmp2);
+  execute_if_needed(dst->depth(), 2, tmp2, dst);
+}
+
+
 auto
 execute_separable(const Device::Pointer &      device,
                   const KernelInfo &           kernel,

clic/src/tier1/gaussian_blur.cpp

@@ -4,6 +4,7 @@
 #include "utils.hpp"
 
 #include "cle_gaussian_blur_separable.h"
+#include "cle_gaussian_derivative_separable.h"
 
 namespace cle::tier1
 {
@@ -29,11 +30,42 @@ gaussian_blur_func(const Device::Pointer & device,
                     temp,
                     dst,
                     { sigma_x, sigma_y, sigma_z },
-                    { sigma2kernelsize(sigma_x), sigma2kernelsize(sigma_y), sigma2kernelsize(sigma_z) });
+                    { sigma2kernelsize(sigma_x), sigma2kernelsize(sigma_y), sigma2kernelsize(sigma_z) },
+                    { 0, 0, 0 });
   return dst;
 }
 
-// generate_angle_matrix_func
-// generate_binary_overlap_matrix_func
+auto
+gaussian_derivative_func(const Device::Pointer & device,
+                         const Array::Pointer &  src,
+                         Array::Pointer          dst,
+                         float                   sigma_x,
+                         float                   sigma_y,
+                         float                   sigma_z,
+                         int                     order_x,
+                         int                     order_y,
+                         int                     order_z) -> Array::Pointer
+{
+  tier0::create_like(src, dst, dType::FLOAT);
+  Array::Pointer temp = src;
+  if (temp->dtype() != dType::FLOAT)
+  {
+    temp = Array::create(dst);
+    tier1::copy_func(device, src, temp);
+  }
+
+  constexpr int        truncate = 8;
+  std::array<float, 3> sigmas = { std::max({ sigma_x, 0.0f }),
+                                  std::max({ sigma_y, 0.0f }),
+                                  std::max({ sigma_z, 0.0f }) };
+  std::array<int, 3>   radii = { static_cast<int>(truncate * sigmas[0] + 0.5f),
+                                 static_cast<int>(truncate * sigmas[1] + 0.5f),
+                                 static_cast<int>(truncate * sigmas[2] + 0.5f) };
+  std::array<int, 3>   orders = { std::min(order_x, 2), std::min(order_y, 2), std::min(order_z, 2) };
+
+  const KernelInfo kernel = { "gaussian_derivative_separable", kernel::gaussian_derivative_separable };
+  execute_separable(device, kernel, temp, dst, sigmas, radii, orders);
+  return dst;
+}
 
 } // namespace cle::tier1

clic/src/tier1/hessian_eigenvalues.cpp

@@ -36,9 +36,9 @@ hessian_eigenvalues_func(const Device::Pointer & device,
   execute(device, kernel, params, range);
   if (src->depth() == 1)
   {
-    return { small_eigenvalue, large_eigenvalue };
+    return { large_eigenvalue, small_eigenvalue };
   }
-  return { small_eigenvalue, middle_eigenvalue, large_eigenvalue };
+  return { large_eigenvalue, middle_eigenvalue, small_eigenvalue };
 }
 
 // inferior_superior(const Device::Pointer & device, const Array::Pointer & src, Array::Pointer dst) -> Array::Pointer

clic/src/tier1/maximum_filter.cpp

@@ -39,7 +39,8 @@ maximum_filter_func(const Device::Pointer & device,
                       src,
                       dst,
                       { static_cast<float>(radius_x), static_cast<float>(radius_y), static_cast<float>(radius_z) },
-                      { r_x, r_y, r_z });
+                      { r_x, r_y, r_z },
+                      { 0, 0, 0 });
   }
   return dst;
 }