Merge pull request #292 from clEsperanto/mode_filters

Mode filters

pyclesperanto_prototype/__init__.py

@@ -10,5 +10,5 @@
 from ._tier10 import *
 from ._tier11 import *
 
-__version__ = "0.23.4"
+__version__ = "0.23.5"
 __common_alias__ = "cle"

pyclesperanto_prototype/_tier1/__init__.py

@@ -99,6 +99,8 @@
 from ._minimum_x_projection import minimum_x_projection
 from ._minimum_y_projection import minimum_y_projection
 from ._minimum_z_projection import minimum_z_projection
+from ._mode_box import mode_box
+from ._mode_sphere import mode_sphere
 from ._modulo_images import modulo_images
 from ._modulo_images import modulo_images as mod
 from ._modulo_images import modulo_images as remainder

pyclesperanto_prototype/_tier1/_mode_box.py

@@ -0,0 +1,57 @@
+import warnings
+
+from .._tier0 import radius_to_kernel_size, execute, create_none, create
+from .._tier0 import plugin_function
+from .._tier0 import Image
+
+
+@plugin_function(categories=['label processing', 'in assistant'], output_creator=create_none)
+def mode_box(source : Image, destination : Image = None, radius_x : int = 1, radius_y : int = 1, radius_z : int = 1) -> Image:
+    """Computes the local mode of a pixels box shaped neighborhood.
+    This can be used to post-process and locally correct semantic segmentation results.
+
+    The box is specified by its half-width and half-height (radius).
+    For technical reasons, the intensities must lie within a range from 0 to 255.
+    In case multiple values have maximum frequency, the smallest one is returned.
+    
+    Parameters
+    ----------
+    source : Image
+    destination : Image, optional
+    radius_x : Number, optional
+    radius_y : Number, optional
+    radius_z : Number, optional
+    
+    Returns
+    -------
+    destination
+    
+    Examples
+    --------
+    >>> import pyclesperanto_prototype as cle
+    >>> cle.mode_box(source, destination, radius_x, radius_y, radius_z)
+    
+    """
+    import numpy as np
+    if source.dtype != np.uint8:
+        warnings.warn("mode_box supports values between 0 and 255 only. Use pixel type uint8 for the source image to mute this warning.")
+
+    if destination is None:
+        destination = create(source, dtype=np.uint8)
+
+    kernel_size_x = radius_to_kernel_size(radius_x)
+    kernel_size_y = radius_to_kernel_size(radius_y)
+    kernel_size_z = radius_to_kernel_size(radius_z)
+
+    parameters = {
+        "dst":destination,
+        "src":source,
+        "Nx":int(kernel_size_x),
+        "Ny":int(kernel_size_y)
+    }
+
+    if (len(destination.shape) == 3):
+        parameters.update({"Nz":int(kernel_size_z)})
+    execute(__file__, 'mode_box_' + str(len(destination.shape)) + 'd_x.cl', 'mode_box_' + str(len(destination.shape)) + 'd', destination.shape, parameters)
+
+    return destination

pyclesperanto_prototype/_tier1/_mode_sphere.py

@@ -0,0 +1,57 @@
+import warnings
+
+from .._tier0 import radius_to_kernel_size, execute, create_none, create
+from .._tier0 import plugin_function
+from .._tier0 import Image
+
+
+@plugin_function(categories=['label processing', 'in assistant'], output_creator=create_none)
+def mode_sphere(source : Image, destination : Image = None, radius_x : int = 1, radius_y : int = 1, radius_z : int = 1) -> Image:
+    """Computes the local mode of a pixels sphere shaped neighborhood.
+    This can be used to post-process and locally correct semantic segmentation results.
+
+    The sphere is specified by its half-width and half-height (radius).
+    For technical reasons, the intensities must lie within a range from 0 to 255.
+    In case multiple values have maximum frequency, the smallest one is returned.
+    
+    Parameters
+    ----------
+    source : Image
+    destination : Image, optional
+    radius_x : Number, optional
+    radius_y : Number, optional
+    radius_z : Number, optional
+    
+    Returns
+    -------
+    destination
+    
+    Examples
+    --------
+    >>> import pyclesperanto_prototype as cle
+    >>> cle.mode_sphere(source, destination, radius_x, radius_y, radius_z)
+    
+    """
+    import numpy as np
+    if source.dtype != np.uint8:
+        warnings.warn("mode_sphere supports values between 0 and 255 only. Use pixel type uint8 for the source image to mute this warning.")
+
+    if destination is None:
+        destination = create(source, dtype=np.uint8)
+
+    kernel_size_x = radius_to_kernel_size(radius_x)
+    kernel_size_y = radius_to_kernel_size(radius_y)
+    kernel_size_z = radius_to_kernel_size(radius_z)
+
+    parameters = {
+        "dst":destination,
+        "src":source,
+        "Nx":int(kernel_size_x),
+        "Ny":int(kernel_size_y)
+    }
+
+    if (len(destination.shape) == 3):
+        parameters.update({"Nz":int(kernel_size_z)})
+    execute(__file__, 'mode_sphere_' + str(len(destination.shape)) + 'd_x.cl', 'mode_sphere_' + str(len(destination.shape)) + 'd', destination.shape, parameters)
+
+    return destination

pyclesperanto_prototype/_tier1/mode_box_2d_x.cl

@@ -0,0 +1,49 @@
+__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+
+__kernel void mode_box_2d
+(
+  IMAGE_dst_TYPE dst,
+  IMAGE_src_TYPE src,
+  const int Nx,
+  const int Ny
+)
+{
+  const int i = get_global_id(0);
+  const int j = get_global_id(1);
+  const int2 coord = (int2){i,j};
+
+  long histogram[256];
+  for (int h = 0; h < 256; h++){
+    histogram[h]=0;
+  }
+
+  const int4   e = (int4)  { (Nx-1)/2, (Ny-1)/2, 0, 0 };
+
+  for (int x = -e.x; x <= e.x; x++) {
+    for (int y = -e.y; y <= e.y; y++) {
+      int x1 = coord.x + x;
+      int x2 = coord.y + y;
+
+      if (x1 < 0 || x2 < 0|| x1 >= GET_IMAGE_WIDTH(src) || x2 >= GET_IMAGE_HEIGHT(src)) {
+        continue;
+      }
+
+      const int2 pos = (int2){x1,x2};
+
+      histogram[(int)READ_src_IMAGE(src,sampler,pos).x]++;
+    }
+  }
+
+  long max_value = 0;
+  int max_pos = 0;
+  for (int i = 0; i < 256; i++){
+    if (max_value < histogram[i]){
+      max_value = histogram[i];
+      max_pos = i;
+    }
+  }
+
+  WRITE_dst_IMAGE(dst, coord, CONVERT_dst_PIXEL_TYPE(max_pos));
+}
+
+

pyclesperanto_prototype/_tier1/mode_box_3d_x.cl

@@ -0,0 +1,55 @@
+__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+
+__kernel void mode_box_3d
+(
+  IMAGE_dst_TYPE dst,
+  IMAGE_src_TYPE src,
+  const int Nx,
+  const int Ny,
+  const int Nz
+)
+{
+  const int i = get_global_id(0);
+  const int j = get_global_id(1);
+  const int k = get_global_id(2);
+  const int4 coord = (int4){i,j,k,0};
+
+  long histogram[256];
+  for (int h = 0; h < 256; h++){
+    histogram[h]=0;
+  }
+
+
+  const int4   e = (int4)  { (Nx-1)/2, (Ny-1)/2, (Nz-1)/2, 0 };
+
+  for (int x = -e.x; x <= e.x; x++) {
+    for (int y = -e.y; y <= e.y; y++) {
+      for (int z = -e.z; z <= e.z; z++) {
+        int x1 = coord.x + x;
+        int x2 = coord.y + y;
+        int x3 = coord.z + z;
+
+        if (x1 < 0 || x2 < 0 || x3 < 0 || x1 >= GET_IMAGE_WIDTH(src) || x2 >= GET_IMAGE_HEIGHT(src) || x3 >= GET_IMAGE_DEPTH(src)) {
+          continue;
+        }
+
+        const int4 pos = (int4){x1,x2,x3,0};
+        int value_res = (int)READ_src_IMAGE(src,sampler,pos).x;
+        histogram[value_res]++;
+      }
+    }
+  }
+
+
+  long max_value = 0;
+  int max_pos = 0;
+  for (int h = 0; h < 256; h++){
+    if (max_value < histogram[h]){
+      max_value = histogram[h];
+      max_pos = h;
+    }
+  }
+
+  WRITE_dst_IMAGE(dst, coord, CONVERT_dst_PIXEL_TYPE(max_pos));
+}
+

pyclesperanto_prototype/_tier1/mode_sphere_2d_x.cl

@@ -0,0 +1,64 @@
+__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+
+__kernel void mode_sphere_2d
+(
+  IMAGE_dst_TYPE dst,
+  IMAGE_src_TYPE src,
+  const int Nx,
+  const int Ny
+)
+{
+  const int i = get_global_id(0);
+  const int j = get_global_id(1);
+  const int2 coord = (int2){i,j};
+
+  int array_size = Nx * Ny;
+  IMAGE_dst_PIXEL_TYPE array[MAX_ARRAY_SIZE];
+
+
+  const int4   e = (int4)  { (Nx-1)/2, (Ny-1)/2, 0, 0 };
+
+  float aSquared = e.x * e.x;
+  float bSquared = e.y * e.y;
+    if (aSquared == 0) {
+        aSquared = FLT_MIN;
+    }
+    if (bSquared == 0) {
+        bSquared = FLT_MIN;
+    }
+
+  long histogram[256];
+  for (int h = 0; h < 256; h++){
+    histogram[h]=0;
+  }
+
+  for (int x = -e.x; x <= e.x; x++) {
+    float xSquared = x * x;
+    for (int y = -e.y; y <= e.y; y++) {
+      float ySquared = y * y;
+      if (xSquared / aSquared + ySquared / bSquared <= 1.0) {
+        int x1 = coord.x + x;
+        int x2 = coord.y + y;
+
+        if (x1 < 0 || x2 < 0|| x1 >= GET_IMAGE_WIDTH(src) || x2 >= GET_IMAGE_HEIGHT(src)) {
+          continue;
+        }
+
+        const int2 pos = (int2){x1,x2};
+
+        histogram[(int)READ_src_IMAGE(src,sampler,pos).x]++;
+      }
+    }
+  }
+
+  long max_value = 0;
+  int max_pos = 0;
+  for (int h = 0; h < 256; h++){
+    if (max_value < histogram[h]){
+      max_value = histogram[h];
+      max_pos = h;
+    }
+  }
+
+  WRITE_dst_IMAGE(dst, coord, CONVERT_dst_PIXEL_TYPE(max_pos));
+}

pyclesperanto_prototype/_tier1/mode_sphere_3d_x.cl

@@ -0,0 +1,76 @@
+__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+
+
+__kernel void mode_sphere_3d
+(
+  IMAGE_dst_TYPE dst,
+  IMAGE_src_TYPE src,
+  const int Nx,
+  const int Ny,
+  const int Nz
+)
+{
+  const int i = get_global_id(0);
+  const int j = get_global_id(1);
+  const int k = get_global_id(2);
+  const int4 coord = (int4){i,j,k,0};
+
+  int array_size = Nx * Ny * Nz;
+  IMAGE_dst_PIXEL_TYPE array[MAX_ARRAY_SIZE];
+
+  // centers
+  const int4   e = (int4)  {(Nx-1)/2, (Ny-1)/2, (Nz-1)/2, 0 };
+
+  long histogram[256];
+  for (int h = 0; h < 256; h++){
+    histogram[h]=0;
+  }
+
+  float aSquared = e.x * e.x;
+  float bSquared = e.y * e.y;
+  float cSquared = e.z * e.z;
+    if (aSquared == 0) {
+        aSquared = FLT_MIN;
+    }
+    if (bSquared == 0) {
+        bSquared = FLT_MIN;
+    }
+    if (cSquared == 0) {
+        cSquared = FLT_MIN;
+    }
+
+  for (int x = -e.x; x <= e.x; x++) {
+    float xSquared = x * x;
+    for (int y = -e.y; y <= e.y; y++) {
+      float ySquared = y * y;
+      for (int z = -e.z; z <= e.z; z++) {
+        float zSquared = z * z;
+        if (xSquared / aSquared + ySquared / bSquared + zSquared / cSquared <= 1.0) {
+          int x1 = coord.x + x;
+          int x2 = coord.y + y;
+          int x3 = coord.z + z;
+
+          if (x1 < 0 || x2 < 0 || x3 < 0 || x1 >= GET_IMAGE_WIDTH(src) || x2 >= GET_IMAGE_HEIGHT(src) || x3 >= GET_IMAGE_DEPTH(src)) {
+            continue;
+          }
+
+          const int4 pos = (int4){x1,x2,x3,0};
+          int value_res = (int)READ_src_IMAGE(src,sampler,pos).x;
+          histogram[value_res]++;
+        }
+      }
+    }
+  }
+
+  long max_value = 0;
+  int max_pos = 0;
+  for (int h = 0; h < 256; h++){
+    if (max_value < histogram[h]){
+      max_value = histogram[h];
+      max_pos = h;
+    }
+  }
+
+  WRITE_dst_IMAGE(dst, coord, CONVERT_dst_PIXEL_TYPE(max_pos));
+}
+

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pyclesperanto_prototype
-version = 0.23.4
+version = 0.23.5
 author = Robert Haase
 author_email = [email protected]
 url = https://github.com/clEsperanto/pyclesperanto_prototype