more work towards graffy/lipase#3 : replaced circle kernels with circle arc kernels, as circle arcs are needed to compute the variance of the image along the circles. Note that at the moment the code doesn't yet compute the variance.

This commit is not intended to change the radial profile image.

src/ij-plugins/Ipr/Lipase/Radial_Profile.py

@@ -1,5 +1,6 @@
 #@ ImagePlus (label="the input image") INPUT_IMAGE
 #@ Float (label="maximal radius", value=10.0, min=0.0, max=100.0, style="slider") MAX_RADIUS
+#@ Int (label="number of angular sectors", value=4, min=0, max=100, style="slider") NUM_ANGULAR_SECTORS
 
 #@output ImagePlus RADIAL_PROFILE
 """This script is supposed to be launched from fiji's jython interpreter
@@ -25,13 +26,13 @@ def run_script():
     global INPUT_IMAGE  # pylint:disable=global-variable-not-assigned
     global MAX_RADIUS  # pylint:disable=global-variable-not-assigned
     global RADIAL_PROFILE  # pylint:disable=global-variable-not-assigned
-
+    global NUM_ANGULAR_SECTORS  # pylint:disable=global-variable-not-assigned
     IImageEngine.set_instance(IJImageEngine())
 
     src_image = IImageEngine.get_instance().create_image(width=1, height=1, pixel_type=PixelType.U8)
     src_image.ij_image = INPUT_IMAGE  # pylint: disable=undefined-variable
 
-    detector = CircularSymmetryDetector(max_radius=MAX_RADIUS)  # pylint: disable=undefined-variable
+    detector = CircularSymmetryDetector(max_radius=MAX_RADIUS, num_angular_sectors=NUM_ANGULAR_SECTORS)  # pylint: disable=undefined-variable
     radial_profiles = detector.compute_radial_profiles(src_image)
 
     RADIAL_PROFILE = radial_profiles.hyperstack

src/lipase/circsymdetector.py

@@ -2,7 +2,8 @@
     an image processing technique to produce a 1D signal for each pixel : this 1D signal is obtained by projecting the neighborhood of this pixel on a set of bases (each base is used as a convilution kernel)
     https://subversion.ipr.univ-rennes1.fr/repos/main/projects/antipode/src/python/antipode/texori.py
 """
-from imageengine import IImageEngine, PixelType
+from imageengine import IImageEngine, PixelType, StackImageFeeder
+import math
 
 class IProjectorBase(object):
     """
@@ -48,64 +49,120 @@ def create_circle_image(image_size, circle_radius, circle_pos, circle_thickness,
             image.set_pixel(x, y, pixel_value)
     return image
 
-class CircularSymmetryProjectorBase(IProjectorBase):
+def create_arc_image(image_size, circle_pos, radius_range, angle_range, background_value=0.0, circle_value=1.0):
     """
-        generates a base of circles 
+    :param dict image_size:
+    """
+    assert radius_range[0] < radius_range[1]
+    assert angle_range[0] < angle_range[1]
+    assert 0.0 <= angle_range[0] <= math.pi * 2.0
+    assert 0.0 <= angle_range[1] <= math.pi * 2.0
+    ie = IImageEngine.get_instance()
+    width = image_size['width']
+    height = image_size['height']
+    assert isinstance(width, int)
+    assert isinstance(height, int)
+    circle_pos_x = float(circle_pos['x'])
+    circle_pos_y = float(circle_pos['y'])
+    r_min_square = radius_range[0] * radius_range[0]
+    r_max_square = radius_range[1] * radius_range[1]
+    image = ie.create_image(width=width, height=height, pixel_type=PixelType.F32)
+    for y in range(height):
+        for x in range(width):
+            dx = x - circle_pos_x
+            dy = y - circle_pos_y
+            r_square = (dx * dx + dy * dy)
+            pixel_value = background_value
+            if r_min_square - 0.001 <= r_square < r_max_square:
+                angle = math.atan2(dy, dx)
+                if angle < 0.0:
+                    angle += math.pi * 2.0
+                assert 0.0 <= angle <= math.pi * 2.0, "unexpected value for angle : %f" % angle
+                if angle_range[0] <= angle < angle_range[1]:
+                    pixel_value = circle_value
+                    # if 0.1 < angle < 0.2 and r_square > 2.0:
+                    #    print("x", x, "y", y, "dx", dx, "dy", dy, "angle", angle)
+                    #    assert False
+            image.set_pixel(x, y, pixel_value)
+    return image
+
+class CircularSymmetryProjectorBase(object):
+    """
+        generates a base of circles arc kernels.
     """
 
-    def __init__(self, max_radius, oversampling_scale=2):
+    def __init__(self, max_radius, num_angular_sectors, num_radial_sectors, oversampling_scale=2):
         """
         :param int max_radius: the biggest circle radius in the set of projectors
+        :param int num_angular_sectors: defines how many sectors we want in each circle
         :param int oversampling_scale: oversampling is used to generate antialased circles. The higher this value, the better the quality of antialiasing is
         """
         super(CircularSymmetryProjectorBase, self).__init__()
         assert max_radius > 0
-        assert oversampling_scale > 0
         self.max_radius = max_radius
+        self.num_radial_sectors = num_radial_sectors
+        self.radial_step = self.max_radius / self.num_radial_sectors
+ 
+        self.num_angular_sectors = num_angular_sectors
+        self.angular_step = math.pi * 2.0 / self.num_angular_sectors
+        assert oversampling_scale > 0
         self.oversampling_scale = oversampling_scale
-        self.num_projectors = int(max_radius) + 1
 
-    def get_num_projectors(self):
+    #def get_num_projectors(self):
-        return self.num_projectors
+    #    return self.num_projectors
         
-    def get_projector_kernel(self, projector_index):
+    def get_projector_kernel(self, radius_index, angular_sector_index):
-        assert projector_index < self.get_num_projectors()
+        """ Computes a kernel containing a patch in the shape of the arc of a circle.
+
+        All patches together form a disc with a radius equal to max_radius.
+        The arc corresponds to the the patch defined by radius_index and angular_sector_index
+
+        :param int radius_index: index of the patch in the radial direction
+
+        """
+        assert radius_index < self.num_radial_sectors
+        assert angular_sector_index < self.num_angular_sectors
+        start_angle = angular_sector_index * self.angular_step
+        end_angle = start_angle + self.angular_step
         oversampling_is_handled = False  # TODO: handle oversampling to increase quality
         if not oversampling_is_handled and self.oversampling_scale != 1  :
             raise NotImplementedError("at the moment, oversampling is not yet implemented")  
-        radius = projector_index
+        radius_index
         image_size = int(self.max_radius) * 2 + 1
         circle_pos = {'x': self.max_radius, 'y': self.max_radius}
-        oversampled_circle = create_circle_image(
+        oversampled_arc = create_arc_image(
             image_size={'width': image_size * self.oversampling_scale, 'height': image_size*self.oversampling_scale},
-            circle_radius=radius * self.oversampling_scale,
             circle_pos={'x': self.max_radius* self.oversampling_scale, 'y': self.max_radius* self.oversampling_scale},
-            circle_thickness=1 * self.oversampling_scale)
+            radius_range=(float(radius_index)* self.radial_step * self.oversampling_scale, float(radius_index + 1) * self.radial_step * self.oversampling_scale),
+            angle_range=(start_angle, end_angle),)
         if self.oversampling_scale == 1:
-            circle_image = oversampled_circle
+            circle_image = oversampled_arc
         else:
             if oversampling_is_handled:
-                circle_image = oversampled_circle.resample(width=image_size, height=image_size)
+                circle_image = oversampled_arc.resample(width=image_size, height=image_size)
         mean_value = circle_image.get_mean_value()
-        assert mean_value > 0.0, "unexpected mean value of this circle image : this circle might be empty, or worse, have negative values, which doesn't make sense"
+        assert mean_value > 0.0, "unexpected mean value of this arc image : this circle might be empty, or worse, have negative values, which doesn't make sense"
-        print(type(mean_value))
         num_pixels = image_size * image_size
         sum_of_pixel_values = mean_value * num_pixels
         # we want each circle to have a total weight of 1.0, regardless their radius
         circle_image.scale_values(1.0/sum_of_pixel_values)
         
-
-        
         anchor_point = {'x': int(circle_pos['x']), 'y': int(circle_pos['y'])}
         return circle_image, anchor_point
 
 
 class CircularSymmetryDetector:
-    def __init__(self, max_radius):
+    def __init__(self, max_radius, num_angular_sectors, num_radial_sectors=None):
         """
-        :para float max_radius:
+        :param float max_radius: in pixels
+        :param int num_angular_sectors: the angle range [0;2 pi] is sampled into num_angular_sectors
+        :param int or None num_radial_sectors: the radius range [0;max_radius] is sampled into num_radial_sectors. If not defined, the number of radial sectors is set so that a radial sector has the width of 1 pixel
         """
         self.max_radius = max_radius
+        self.num_angular_sectors = num_angular_sectors
+        if num_radial_sectors is None:
+            num_radial_sectors = int(max_radius)
+        self.num_radial_sectors = num_radial_sectors
 
     def compute_radial_profiles(self, src_image):
         """ Computes for each pixel the radial profile (with this pixel as center)
@@ -115,13 +172,18 @@ class CircularSymmetryDetector:
         """
         ie = IImageEngine.get_instance()
         ie.debugger.on_image(src_image, 'src_image')
-        projector_base = CircularSymmetryProjectorBase(self.max_radius, oversampling_scale=1)
+        projector_base = CircularSymmetryProjectorBase(self.max_radius, num_angular_sectors=self.num_angular_sectors, num_radial_sectors=self.num_radial_sectors, oversampling_scale=1)
-        radial_profile_image = ie.create_hyperstack(width=src_image.get_width(), height=src_image.get_height(), num_channels=projector_base.get_num_projectors(), num_slices=1, num_frames=1, pixel_type=PixelType.F32)
+        radial_profile_image = ie.create_hyperstack(width=src_image.get_width(), height=src_image.get_height(), num_channels=projector_base.num_radial_sectors, num_slices=1, num_frames=1, pixel_type=PixelType.F32)
-        for projector_index in range(projector_base.get_num_projectors()):
+        for radius_index in range(projector_base.num_radial_sectors):
-            projector, center_of_filter = projector_base.get_projector_kernel(projector_index)
+            circle_stack = ie.create_hyperstack(width=src_image.get_width(), height=src_image.get_height(), num_channels=projector_base.num_angular_sectors, num_slices=1, num_frames=1, pixel_type=PixelType.F32)
-            ie.debugger.on_image(projector, 'projector_%d' % (projector_index))
+            for angular_sector_index in range(projector_base.num_angular_sectors):
-            print(type(center_of_filter))
+                arc_kernel, center_of_filter = projector_base.get_projector_kernel(radius_index, angular_sector_index)
-            projection = ie.filter2D(src_image, dst_type=PixelType.F32, kernel=projector, anchor=(center_of_filter['x'], center_of_filter['y']))
+                ie.debugger.on_image(arc_kernel, 'projector_%d_%d' % (radius_index, angular_sector_index))
-            ie.debugger.on_image(projection, 'projection_%d' % (projector_index))
+                projection = ie.filter2D(src_image, dst_type=PixelType.F32, kernel=arc_kernel, anchor=(center_of_filter['x'], center_of_filter['y']))
-            radial_profile_image.set_image(projection, frame_index=0, slice_index=0, channel_index=projector_index)
+                # ie.debugger.on_image(projection, 'projection_%d' % (angular_sector_index))
+                circle_stack.set_image(projection, frame_index=0, slice_index=0, channel_index=angular_sector_index)
+            # compute the image in which each pixel contains the mean value of the source image in the given circular region around the pixel
+            circle_mean_image = ie.compute_mean(StackImageFeeder(circle_stack))
+
+            radial_profile_image.set_image(circle_mean_image, frame_index=0, slice_index=0, channel_index=radius_index)
         return radial_profile_image

src/lipase/imageengine.py

@@ -180,8 +180,8 @@ class IImageFeeder(ABC):
         for image in it:
             hyperstack.set_image(image, frame_index=frame_index)
             frame_index += 1
-        print(frame_index)
+        # print(frame_index)
-        print(self.get_num_images())
+        # print(self.get_num_images())
         
         assert frame_index == self.get_num_images()
         return hyperstack
@@ -264,7 +264,7 @@ class StackImageFeeder(IImageFeeder):
         return image
 
     def get_num_images(self):
-        return self.hyperstack.num_frames()
+        return self.hyperstack.num_frames() * self.hyperstack.num_slices() * self.hyperstack.num_channels()
 
 
 class IImageProcessingDebugger(ABC):
@@ -307,7 +307,7 @@ class FileBasedDebugger(IImageProcessingDebugger):
         mkdir_p(self.debug_images_root_path)
 
     def on_image(self, image, image_id):
-        print('FileBasedDebugger.on_image : image_id = %s' % image_id)
+        # print('FileBasedDebugger.on_image : image_id = %s' % image_id)
         ie = IImageEngine.get_instance()
         ie.save_as_tiff(image, '%s/%s.tiff' % (self.debug_images_root_path, image_id))
 
@@ -426,6 +426,14 @@ class IImageEngine(ABC):
         :rtype IImage:
         """
 
+    @abc.abstractmethod
+    def compute_mean(self, image_feeder):
+        """Compute for each pixel position the mean value at this position in all input images.
+
+        :param IImageFeeder image_feeder:
+        :rtype IImage:
+        """
+
     @abc.abstractmethod
     def mean_filter(self, image, radius):
         """Each pixel becomes an average of its neighbours within the given radius

src/lipase/imagej/ijimageengine.py

@@ -314,10 +314,8 @@ class IJImageEngine(IImageEngine):
         # for image_file_path in images_file_path[2:-1]:
         for other_image in it:
             # other_image = IJ.openImage(image_file_path)
-            print('other_image', other_image)
             ic = ImageCalculator()
             ic.run("max", max_image, other_image.ij_image)
-        print('max_image', max_image)
         return IJImage(self, max_image)
 
     def compute_median(self, image_feeder):
@@ -335,12 +333,19 @@ class IJImageEngine(IImageEngine):
         # median_image = IJ.getImage()  # get the currently selected image
         return IJImage(self, median_image)
 
+    def compute_mean(self, image_feeder):
+        hyperstack = image_feeder.create_hyperstack()
+        # https://imagej.nih.gov/ij/developer/api/ij/plugin/ZProjector.html
+        projector = ZProjector()
+        mean_image = projector.run(hyperstack.hyperstack, 'average')
+        # imagej_run_image_command(image=hyperstack.hyperstack, command="Z Project...", options="projection=Mean")
+        return IJImage(self, mean_image)
+
     def mean_filter(self, image, radius):
         """Implement interface method."""
         IJ.run(image.ij_image, "Mean...", "radius=%d" % radius)
 
     def filter2D(self, src_image, dst_type, kernel, anchor=(-1, -1)):
-        print(type(anchor), anchor)
         if anchor == (-1, -1):
             assert kernel.get_width() % 2 == 1 and kernel.get_height() % 2 == 1, "kernel sizes are expected to be odd if you want the anchor to be at the center of the kernel"
 
@@ -358,7 +363,6 @@ class IJImageEngine(IImageEngine):
         # warning ! trying to guess the arguments from opencv's documentation is time consuming as the error messages are misleading (in the following call, javacpp complains that the 1st argument is not a org.bytedeco.javacpp.opencv_core$Mat, while it is ! The problem comes from the order of the arguments). So, instead of guessing, the found accepted signatures of filter2D are in https://github.com/bytedeco/javacpp-presets/blob/master/opencv/src/gen/java/org/bytedeco/opencv/global/opencv_imgproc.java
         # opencv_imgproc.morphologyEx(cv_src_image, opencv_imgproc.MORPH_OPEN, struct_element, dst_image)
         # TypeError: morphologyEx(): 1st arg can't be coerced to org.bytedeco.javacpp.opencv_core$GpuMat, org.bytedeco.javacpp.opencv_core$Mat, org.bytedeco.javacpp.opencv_core$UMat
-        print('before opencv_imgproc.filter2D')
         delta = 0
         border_type = opencv_core.BORDER_DEFAULT
         if dst_type is None:
@@ -366,7 +370,6 @@ class IJImageEngine(IImageEngine):
         else:
             ddepth = PIXEL_TYPE_TO_CV_PIXEL_TYPE[dst_type]
         opencv_imgproc.filter2D(cv_src_image, cv_dst_image, ddepth, cv_kernel, cv_anchor, delta, border_type)
-        print('after opencv_imgproc.filter2D')
 
         mat2imp = MatImagePlusConverter()
         dst_image.ij_image.setProcessor(mat2imp.toImageProcessor(cv_dst_image))