morphological_chan_vese

pyclesperanto_prototype/_tier1/__init__.py

@@ -68,6 +68,7 @@
 from ._greater_or_equal import greater_or_equal
 from ._greater_or_equal_constant import greater_or_equal_constant
 from ._hessian_eigenvalues import hessian_eigenvalues
+from ._inferior_superior import inferior_superior
 from ._laplace_box import laplace_box
 from ._laplace_diamond import laplace_diamond
 from ._logarithm import logarithm
@@ -161,6 +162,7 @@
 from ._sum_x_projection import sum_x_projection
 from ._sum_y_projection import sum_y_projection
 from ._sum_z_projection import sum_z_projection
+from ._superior_inferior import superior_inferior
 from ._transpose_xy import transpose_xy
 from ._transpose_xz import transpose_xz
 from ._transpose_yz import transpose_yz

pyclesperanto_prototype/_tier1/_inferior_superior.py

@@ -0,0 +1,40 @@
+from .._tier0 import execute, create_binary_like
+from .._tier0 import plugin_function
+from .._tier0 import Image
+
+@plugin_function(categories=['binary processing'], output_creator=create_binary_like)
+def inferior_superior(source : Image, destination : Image = None) -> Image:
+    """Dilates the image respectively with a number of kernels and takes the minimum
+    value across all results for each pixel.
+
+    Parameters
+    ----------
+    source : Image
+    destination : Image, optional
+
+    Returns
+    -------
+    destination
+
+    Examples
+    --------
+    >>> import pyclesperanto_prototype as cle
+    >>> cle.inferior_superior(source, destination)
+
+    References
+    ----------
+    Implemented in inf_sup function in scikit morphological snakes:
+    .. [1] https://github.com/scikit-image/scikit-image/blob/00177e14097237ef20ed3141ed454bc81b308f82/skimage/segmentation/morphsnakes.py
+    """
+    import numpy as np
+    if source.dtype != np.uint8:
+        # the read function in the kernel below is custom and only supports images of type uint8
+        source = source.astype(np.uint8)
+
+    parameters = {
+        "src":source,
+        "dst":destination
+    }
+
+    execute(__file__, './inferior_superior_' + str(len(destination.shape)) + 'd_x.cl', 'inferior_superior_' + str(len(destination.shape)) + 'd', destination.shape, parameters)
+    return destination

pyclesperanto_prototype/_tier1/_superior_inferior.py

@@ -0,0 +1,40 @@
+from .._tier0 import execute, create_binary_like
+from .._tier0 import plugin_function
+from .._tier0 import Image
+
+@plugin_function(categories=['binary processing'], output_creator=create_binary_like)
+def superior_inferior(source : Image, destination : Image = None) -> Image:
+    """Erodes the image respectively with a number of kernels and takes the maximum
+    value across all results for each pixel.
+
+    Parameters
+    ----------
+    source : Image
+    destination : Image, optional
+
+    Returns
+    -------
+    destination
+
+    Examples
+    --------
+    >>> import pyclesperanto_prototype as cle
+    >>> cle.inferior_superior(source, destination)
+
+    References
+    ----------
+    Implemented in sup_inf function in scikit morphological snakes:
+    .. [1] https://github.com/scikit-image/scikit-image/blob/00177e14097237ef20ed3141ed454bc81b308f82/skimage/segmentation/morphsnakes.py
+    """
+    import numpy as np
+    if source.dtype != np.uint8:
+        # the read function in the kernel below is custom and only supports images of type uint8
+        source = source.astype(np.uint8)
+
+    parameters = {
+        "src":source,
+        "dst":destination
+    }
+
+    execute(__file__, './superior_inferior_' + str(len(destination.shape)) + 'd_x.cl', 'superior_inferior_' + str(len(destination.shape)) + 'd', destination.shape, parameters)
+    return destination

pyclesperanto_prototype/_tier1/inferior_superior_2d_x.cl

@@ -0,0 +1,83 @@
+__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+
+#define READ_IMAGE_ZERO_OUTSIDE(a,b,c) read_buffer2duc_zero_outside(GET_IMAGE_WIDTH(a),GET_IMAGE_HEIGHT(a),GET_IMAGE_DEPTH(a),a,b,c)
+
+inline uchar2 read_buffer2duc_zero_outside(int read_buffer_width, int read_buffer_height, int read_buffer_depth, __global uchar * buffer_var, sampler_t sampler, int2 position )
+{
+    int2 pos = (int2){position.x, position.y};
+    int pos_in_buffer = pos.x + pos.y * read_buffer_width;
+    if (pos.x < 0 || pos.x >= read_buffer_width || pos.y < 0 || pos.y >= read_buffer_height) {
+        return (uchar2){0, 0};
+    }
+    return (uchar2){buffer_var[pos_in_buffer],0};
+}
+
+__kernel void inferior_superior_2d (
+    IMAGE_src_TYPE  src,
+    IMAGE_dst_TYPE  dst
+)
+{
+  const int x = get_global_id(0);
+  const int y = get_global_id(1);
+
+  const int2 pos = (int2){x,y};
+
+  float value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, pos).x;
+
+  // if value is already 1, dilate will return 1
+  if (value == 1) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1));
+    return;
+  }
+
+  /* Dilate with kernel [[1, 0, 0], 
+                         [0, 1, 0], 
+                         [0, 0, 1]] */
+  value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, 1})).x;
+  if (value == 0) {
+    value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, -1})).x;
+    if (value == 0) {
+      WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+      return;
+    }
+  }
+
+  /* Dilate with kernel [[0, 1, 0], 
+                         [0, 1, 0], 
+                         [0, 1, 0]] */
+  value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){0, 1})).x;
+    if (value == 0) {
+      value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){0, -1})).x;
+      if (value == 0) {
+        WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+        return;
+      }
+    }
+
+  /* Dilate with kernel [[0, 0, 1], 
+                         [0, 1, 0], 
+                         [1, 0, 0]] */
+  value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, 1})).x;
+    if (value == 0) {
+      value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, -1})).x;
+      if (value == 0) {
+        WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+        return;
+      }
+    }
+
+  /* Dilate with kernel [[0, 0, 0], 
+                         [1, 1, 1], 
+                         [0, 0, 0]] */
+  value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, 0})).x;
+    if (value == 0) {
+      value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, 0})).x;
+      if (value == 0) {
+        WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+        return;
+      }
+    }
+
+  // If all dilates are 1 then return 1
+  WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1));
+}

pyclesperanto_prototype/_tier1/inferior_superior_3d_x.cl

@@ -0,0 +1,194 @@
+__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+
+
+#define READ_IMAGE_ZERO_OUTSIDE(a,b,c) read_buffer3duc_zero_outside(GET_IMAGE_WIDTH(a),GET_IMAGE_HEIGHT(a),GET_IMAGE_DEPTH(a),a,b,c)
+
+inline uchar2 read_buffer3duc_zero_outside(int read_buffer_width, int read_buffer_height, int read_buffer_depth, __global uchar * buffer_var, sampler_t sampler, int4 position )
+{
+    int4 pos = (int4){position.x, position.y, position.z, 0};
+    int pos_in_buffer = pos.x + pos.y * read_buffer_width + pos.z * read_buffer_width * read_buffer_height;
+    if (pos.x < 0 || pos.x >= read_buffer_width || pos.y < 0 || pos.y >= read_buffer_height || pos.z < 0 || pos.z >= read_buffer_depth) {
+        return (uchar2){0, 0};
+    }
+    return (uchar2){buffer_var[pos_in_buffer],0};
+}
+
+__kernel void inferior_superior_3d (
+    IMAGE_src_TYPE  src,
+    IMAGE_dst_TYPE  dst
+)
+{
+  const int x = get_global_id(0);
+  const int y = get_global_id(1);
+  const int z = get_global_id(2);
+
+  const int4 pos = (int4){x, y, z, 0};
+
+  float value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, pos).x;
+
+  // if value is already 0, erode will return 0
+  if (value != 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1));
+    return;
+  }
+
+  // printf("pixel at coord x%i y%i z%i\n", x, y, z);
+
+
+  // P0
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, 0, 0})).x;
+        // printf("value %i\n", value);
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P1
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, 0, j, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P2
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){0, i, j, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P3
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, j, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P4
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){j, i, -i, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P5
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, i, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P6
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, -i, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P7
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, i, j, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // P8
+  for (int i = -1; i <= 1; i++) {
+      for (int j = -1; j <= 1; j++) {
+        value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, -i, j, 0})).x;
+        if (value != 0) {
+          break;
+        }
+      }
+      if (value != 0) {
+        break;
+      }
+    }
+  if (value == 0) {
+    WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0));
+    return;
+  }
+
+  // If all erodes are 0 then return 0
+  WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1));
+}