Merge pull request #264 from clEsperanto/min-max-along-points

Min max along points

pyclesperanto_prototype/__init__.py

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

pyclesperanto_prototype/_tier2/__init__.py

@@ -14,7 +14,9 @@
 from ._flag_existing_labels import flag_existing_labels as flag_existing_intensities
 from ._degrees_to_radians import degrees_to_radians
 from ._gamma_correction import gamma_correction
+from ._generate_maximum_intensity_between_points_matrix import generate_maximum_intensity_between_points_matrix
 from ._generate_mean_intensity_between_points_matrix import generate_mean_intensity_between_points_matrix
+from ._generate_minimum_intensity_between_points_matrix import generate_minimum_intensity_between_points_matrix
 from ._generate_standard_deviation_intensity_between_points_matrix import generate_standard_deviation_intensity_between_points_matrix
 from ._generate_should_touch_matrix import generate_should_touch_matrix
 from ._invert import invert

pyclesperanto_prototype/_tier2/_generate_maximum_intensity_between_points_matrix.py

@@ -0,0 +1,47 @@
+from .._tier0 import execute, plugin_function, Image, create_none, create_matrix_from_pointlists, create_like
+
+
+@plugin_function(output_creator=create_none)
+def generate_maximum_intensity_between_points_matrix(intensity_image: Image, pointlist: Image, touch_matrix: Image = None,
+                                                     maximum_intensity_matrix_destination: Image = None,
+                                                     num_samples: int = 10):
+    """Determine the maximum intensity between pairs of point coordinates and
+    write them in a matrix.
+
+    Parameters
+    ----------
+    intensity_image: Image
+        image where the intensity will be measured
+    pointlist: Image
+        list of coordinates
+    touch_matrix: Image, optional
+        if only selected pairs should be measured, use this binary matrix to confige which
+    maximum_intensity_matrix_destination: Image, optional
+        matrix where the results are written ito
+    num_samples: int, optional
+        Number of samples to take along the line for averaging, default = 10
+
+    Returns
+    -------
+    average_intensity_matrix_destination
+    """
+    from .._tier1 import set
+
+    if maximum_intensity_matrix_destination is None:
+        maximum_intensity_matrix_destination = create_matrix_from_pointlists(pointlist, pointlist)
+
+    if touch_matrix is None:
+        touch_matrix = create_like(maximum_intensity_matrix_destination)
+        set(touch_matrix, 1)
+
+    parameters = {
+        "src_touch_matrix": touch_matrix,
+        "src_pointlist": pointlist,
+        "src_intensity": intensity_image,
+        "dst_maximum_intensity_matrix": maximum_intensity_matrix_destination,
+        "num_samples": int(num_samples)
+    }
+
+    execute(__file__, 'maximum_intensity_between_points_matrix_x.cl', 'maximum_intensity_between_points_matrix', touch_matrix.shape, parameters)
+
+    return maximum_intensity_matrix_destination

pyclesperanto_prototype/_tier2/_generate_mean_intensity_between_points_matrix.py

@@ -42,6 +42,6 @@ def generate_mean_intensity_between_points_matrix(intensity_image: Image, pointl
         "num_samples": int(num_samples)
     }
 
-    execute(__file__, 'mean_intensity_along_line_x.cl', 'mean_intensity_along_line', touch_matrix.shape, parameters)
+    execute(__file__, 'mean_intensity_between_points_matrix_x.cl', 'mean_intensity_between_points_matrix', touch_matrix.shape, parameters)
 
     return mean_intensity_matrix_destination

pyclesperanto_prototype/_tier2/_generate_minimum_intensity_between_points_matrix.py

@@ -0,0 +1,47 @@
+from .._tier0 import execute, plugin_function, Image, create_none, create_matrix_from_pointlists, create_like
+
+
+@plugin_function(output_creator=create_none)
+def generate_minimum_intensity_between_points_matrix(intensity_image: Image, pointlist: Image, touch_matrix: Image = None,
+                                                     minimum_intensity_matrix_destination: Image = None,
+                                                     num_samples: int = 10):
+    """Determine the minimum intensity between pairs of point coordinates and
+    write them in a matrix.
+
+    Parameters
+    ----------
+    intensity_image: Image
+        image where the intensity will be measured
+    pointlist: Image
+        list of coordinates
+    touch_matrix: Image, optional
+        if only selected pairs should be measured, use this binary matrix to confige which
+    minimum_intensity_matrix_destination: Image, optional
+        matrix where the results are written ito
+    num_samples: int, optional
+        Number of samples to take along the line for averaging, default = 10
+
+    Returns
+    -------
+    average_intensity_matrix_destination
+    """
+    from .._tier1 import set
+
+    if minimum_intensity_matrix_destination is None:
+        minimum_intensity_matrix_destination = create_matrix_from_pointlists(pointlist, pointlist)
+
+    if touch_matrix is None:
+        touch_matrix = create_like(minimum_intensity_matrix_destination)
+        set(touch_matrix, 1)
+
+    parameters = {
+        "src_touch_matrix": touch_matrix,
+        "src_pointlist": pointlist,
+        "src_intensity": intensity_image,
+        "dst_minimum_intensity_matrix": minimum_intensity_matrix_destination,
+        "num_samples": int(num_samples)
+    }
+
+    execute(__file__, 'minimum_intensity_between_points_matrix_x.cl', 'minimum_intensity_between_points_matrix', touch_matrix.shape, parameters)
+
+    return minimum_intensity_matrix_destination

pyclesperanto_prototype/_tier2/_generate_standard_deviation_intensity_between_points_matrix.py

@@ -46,6 +46,6 @@ def generate_standard_deviation_intensity_between_points_matrix(intensity_image:
         "num_samples": int(num_samples)
     }
 
-    execute(__file__, 'standard_deviation_intensity_along_line_x.cl', 'standard_deviation_intensity_along_line', touch_matrix.shape, parameters)
+    execute(__file__, 'standard_deviation_intensity_between_points_matrix_x.cl', 'standard_deviation_intensity_between_points_matrix', touch_matrix.shape, parameters)
 
     return standard_deviation_intensity_matrix_destination

pyclesperanto_prototype/_tier2/maximum_intensity_between_points_matrix_x.cl

@@ -0,0 +1,59 @@
+__kernel void maximum_intensity_between_points_matrix (
+    IMAGE_src_touch_matrix_TYPE src_touch_matrix,
+    IMAGE_src_pointlist_TYPE src_pointlist,
+    IMAGE_src_intensity_TYPE src_intensity,                                            
+    IMAGE_dst_maximum_intensity_matrix_TYPE dst_maximum_intensity_matrix,
+    int num_samples
+)
+{
+  const sampler_t intsampler  = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_NONE | CLK_FILTER_NEAREST;
+
+  const int touch_x = get_global_id(0);
+  const int touch_y = get_global_id(1);
+
+  const int touching = READ_IMAGE(src_touch_matrix, intsampler, POS_src_touch_matrix_INSTANCE(touch_x, touch_y, 0, 0)).x;
+  if (touching == 0 || touch_x == 0 || touch_y == 0) {
+    WRITE_IMAGE (dst_maximum_intensity_matrix, POS_dst_maximum_intensity_matrix_INSTANCE(touch_x, touch_y,0,0), 0);
+    return;
+  }
+
+  const float x1 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_x-1, 0, 0, 0)).x;
+  const float y1 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_x-1, 1, 0, 0)).x;
+  const float z1 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_x-1, 2, 0, 0)).x;
+
+  const float x2 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_y-1, 0, 0, 0)).x;
+  const float y2 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_y-1, 1, 0, 0)).x;
+  const float z2 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_y-1, 2, 0, 0)).x;
+
+
+  float4 directionVector = (float4){x2 - x1, y2 - y1, z2 - z1, 0};
+
+  // const float len = length(directionVector);
+  directionVector.x = directionVector.x / (num_samples - 1);
+  directionVector.y = directionVector.y / (num_samples - 1);
+  directionVector.z = directionVector.z / (num_samples - 1);
+
+  if (touch_x == 1 && touch_y == 2) {
+    printf("DIR %f / %f \n", directionVector.x, directionVector.y);
+  }
+
+  int width = GET_IMAGE_WIDTH(src_intensity);
+  int height = GET_IMAGE_HEIGHT(src_intensity);
+  int depth = GET_IMAGE_DEPTH(src_intensity);
+
+  float4 position = (float4){x1, y1, z1, 0};
+  float maximum = 0;
+
+  for (int i = 0; i < num_samples; i++) {
+      POS_src_intensity_TYPE pos = POS_src_intensity_INSTANCE((int)(position.x + 0.5), (int)(position.y + 0.5), (int)(position.z + 5), 0);
+
+      float value = (float)(READ_IMAGE(src_intensity, intsampler, pos).x);
+      if (maximum < value || i == 0) {
+        maximum = value;
+      }
+
+      position = position + directionVector;
+  }
+
+  WRITE_IMAGE (dst_maximum_intensity_matrix, POS_dst_maximum_intensity_matrix_INSTANCE(touch_x, touch_y,0,0), maximum);
+}

pyclesperanto_prototype/_tier2/mean_intensity_along_line_x.cl → pyclesperanto_prototype/_tier2/mean_intensity_between_points_matrix_x.cl

@@ -1,4 +1,4 @@
-__kernel void mean_intensity_along_line (
+__kernel void mean_intensity_between_points_matrix (
     IMAGE_src_touch_matrix_TYPE src_touch_matrix,
     IMAGE_src_pointlist_TYPE src_pointlist,
     IMAGE_src_intensity_TYPE src_intensity,                                            
@@ -50,10 +50,6 @@ __kernel void mean_intensity_along_line (
       float value = (float)(READ_IMAGE(src_intensity, intsampler, pos).x);
       sum = sum + value;
 
-      if (touch_x == 1 && touch_y == 2) {
-        printf("position %f / %f \n", position.x, position.y);     
-        printf("pos %d / %d \n", pos.x, pos.y);     
-      }
       position = position + directionVector;
   }
 

pyclesperanto_prototype/_tier2/minimum_intensity_between_points_matrix_x.cl

@@ -0,0 +1,59 @@
+__kernel void minimum_intensity_between_points_matrix (
+    IMAGE_src_touch_matrix_TYPE src_touch_matrix,
+    IMAGE_src_pointlist_TYPE src_pointlist,
+    IMAGE_src_intensity_TYPE src_intensity,                                            
+    IMAGE_dst_minimum_intensity_matrix_TYPE dst_minimum_intensity_matrix,
+    int num_samples
+)
+{
+  const sampler_t intsampler  = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_NONE | CLK_FILTER_NEAREST;
+
+  const int touch_x = get_global_id(0);
+  const int touch_y = get_global_id(1);
+
+  const int touching = READ_IMAGE(src_touch_matrix, intsampler, POS_src_touch_matrix_INSTANCE(touch_x, touch_y, 0, 0)).x;
+  if (touching == 0 || touch_x == 0 || touch_y == 0) {
+    WRITE_IMAGE (dst_minimum_intensity_matrix, POS_dst_minimum_intensity_matrix_INSTANCE(touch_x, touch_y,0,0), 0);
+    return;
+  }
+
+  const float x1 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_x-1, 0, 0, 0)).x;
+  const float y1 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_x-1, 1, 0, 0)).x;
+  const float z1 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_x-1, 2, 0, 0)).x;
+
+  const float x2 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_y-1, 0, 0, 0)).x;
+  const float y2 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_y-1, 1, 0, 0)).x;
+  const float z2 = READ_IMAGE(src_pointlist, intsampler, POS_src_pointlist_INSTANCE(touch_y-1, 2, 0, 0)).x;
+
+
+  float4 directionVector = (float4){x2 - x1, y2 - y1, z2 - z1, 0};
+
+  // const float len = length(directionVector);
+  directionVector.x = directionVector.x / (num_samples - 1);
+  directionVector.y = directionVector.y / (num_samples - 1);
+  directionVector.z = directionVector.z / (num_samples - 1);
+
+  if (touch_x == 1 && touch_y == 2) {
+    printf("DIR %f / %f \n", directionVector.x, directionVector.y);
+  }
+
+  int width = GET_IMAGE_WIDTH(src_intensity);
+  int height = GET_IMAGE_HEIGHT(src_intensity);
+  int depth = GET_IMAGE_DEPTH(src_intensity);
+
+  float4 position = (float4){x1, y1, z1, 0};
+  float minimum = 0;
+
+  for (int i = 0; i < num_samples; i++) {
+      POS_src_intensity_TYPE pos = POS_src_intensity_INSTANCE((int)(position.x + 0.5), (int)(position.y + 0.5), (int)(position.z + 5), 0);
+
+      float value = (float)(READ_IMAGE(src_intensity, intsampler, pos).x);
+      if (minimum > value || i == 0) {
+        minimum = value;
+      }
+
+      position = position + directionVector;
+  }
+
+  WRITE_IMAGE (dst_minimum_intensity_matrix, POS_dst_minimum_intensity_matrix_INSTANCE(touch_x, touch_y,0,0), minimum);
+}

pyclesperanto_prototype/_tier2/standard_deviation_intensity_along_line_x.cl → pyclesperanto_prototype/_tier2/standard_deviation_intensity_between_points_matrix_x.cl

@@ -1,4 +1,4 @@
-__kernel void standard_deviation_intensity_along_line (
+__kernel void standard_deviation_intensity_between_points_matrix (
     IMAGE_src_touch_matrix_TYPE src_touch_matrix,
     IMAGE_src_pointlist_TYPE src_pointlist,
     IMAGE_src_intensity_TYPE src_intensity,                                            
@@ -54,10 +54,6 @@ __kernel void standard_deviation_intensity_along_line (
       sum = sum + value;
       sq_sum = sq_sum + value * value;
 
-      if (touch_x == 1 && touch_y == 2) {
-        printf("position %f / %f \n", position.x, position.y);     
-        printf("pos %d / %d \n", pos.x, pos.y);     
-      }
       position = position + directionVector;
   }
   const float mean = sum / num_samples;

setup.cfg

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

tests/test_generate_maximum_intensity_between_points_matrix.py

@@ -0,0 +1,27 @@
+def test_generate_maximum_intensity_between_points_matrix():
+    import pyclesperanto_prototype as cle
+    import numpy as np
+
+    image = np.zeros((10, 10))
+    image[5:, :5] = 5
+    image[5:, 5:] = 10
+
+    coords = np.asarray([
+        [1, 1],  # first point
+        [8, 1],  # second ...
+        [1, 8],
+        [8, 8]
+    ]).T
+
+    touch_matrix = cle.generate_distance_matrix(coords, coords) > 0
+
+    average_intensity_matrix = cle.generate_maximum_intensity_between_points_matrix(image, coords, touch_matrix)
+
+    reference = np.asarray([
+        [0.,  0.,  0.,  0., 0.],
+        [0.,  0.,  0.,  5., 10.],
+        [0.,  0.,  0.,  5., 10.],
+        [0.,  5.,  5.,  0., 10],
+        [0.,  10., 10., 10, 0.]])
+
+    assert cle.array_equal(reference, average_intensity_matrix)

tests/test_generate_mean_intensity_between_points_matrix.py

@@ -1,4 +1,4 @@
-def test_generate_average_intensity_between_points_matrix():
+def test_generate_mean_intensity_between_points_matrix():
     import pyclesperanto_prototype as cle
     import numpy as np
 

tests/test_generate_minimum_intensity_between_points_matrix.py

@@ -0,0 +1,27 @@
+def test_generate_minimum_intensity_between_points_matrix():
+    import pyclesperanto_prototype as cle
+    import numpy as np
+
+    image = np.zeros((10, 10))
+    image[5:, :5] = 5
+    image[5:, 5:] = 10
+
+    coords = np.asarray([
+        [1, 1],  # first point
+        [8, 1],  # second ...
+        [1, 8],
+        [8, 8]
+    ]).T
+
+    touch_matrix = cle.generate_distance_matrix(coords, coords) > 0
+
+    average_intensity_matrix = cle.generate_minimum_intensity_between_points_matrix(image, coords, touch_matrix)
+
+    reference = np.asarray([
+        [0.,  0.,  0.,  0., 0.],
+        [0.,  0.,  0.,  0,  0.],
+        [0.,  0.,  0.,  0,  0.],
+        [0.,  0.,  0.,  0., 5.],
+        [0.,  0.,  0.,  5., 0.]])
+
+    assert cle.array_equal(reference, average_intensity_matrix)