From 608fefce0a21309294e4df89fdded44967e13556 Mon Sep 17 00:00:00 2001 From: Guillaume Raffy Date: Tue, 7 Apr 2020 14:15:54 +0200 Subject: [PATCH] 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 | 5 +- src/lipase/circsymdetector.py | 124 +++++++++++++++----- src/lipase/imageengine.py | 16 ++- src/lipase/imagej/ijimageengine.py | 13 +- 4 files changed, 116 insertions(+), 42 deletions(-) diff --git a/src/ij-plugins/Ipr/Lipase/Radial_Profile.py b/src/ij-plugins/Ipr/Lipase/Radial_Profile.py index 20fb184..d7a3361 100644 --- a/src/ij-plugins/Ipr/Lipase/Radial_Profile.py +++ b/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 diff --git a/src/lipase/circsymdetector.py b/src/lipase/circsymdetector.py index 89e1a1a..bac51a5 100644 --- a/src/lipase/circsymdetector.py +++ b/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): - return self.num_projectors + #def get_num_projectors(self): + # return self.num_projectors - def get_projector_kernel(self, projector_index): - assert projector_index < self.get_num_projectors() + def get_projector_kernel(self, radius_index, angular_sector_index): + """ 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" - print(type(mean_value)) + 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" 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) - 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) - for projector_index in range(projector_base.get_num_projectors()): - projector, center_of_filter = projector_base.get_projector_kernel(projector_index) - ie.debugger.on_image(projector, 'projector_%d' % (projector_index)) - print(type(center_of_filter)) - projection = ie.filter2D(src_image, dst_type=PixelType.F32, kernel=projector, anchor=(center_of_filter['x'], center_of_filter['y'])) - ie.debugger.on_image(projection, 'projection_%d' % (projector_index)) - radial_profile_image.set_image(projection, frame_index=0, slice_index=0, channel_index=projector_index) + 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.num_radial_sectors, num_slices=1, num_frames=1, pixel_type=PixelType.F32) + for radius_index in range(projector_base.num_radial_sectors): + 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) + for angular_sector_index in range(projector_base.num_angular_sectors): + arc_kernel, center_of_filter = projector_base.get_projector_kernel(radius_index, angular_sector_index) + ie.debugger.on_image(arc_kernel, 'projector_%d_%d' % (radius_index, angular_sector_index)) + projection = ie.filter2D(src_image, dst_type=PixelType.F32, kernel=arc_kernel, anchor=(center_of_filter['x'], center_of_filter['y'])) + # 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 diff --git a/src/lipase/imageengine.py b/src/lipase/imageengine.py index 5d05e62..34be352 100644 --- a/src/lipase/imageengine.py +++ b/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(self.get_num_images()) + # print(frame_index) + # 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 diff --git a/src/lipase/imagej/ijimageengine.py b/src/lipase/imagej/ijimageengine.py index 4ca73af..7dfc52e 100644 --- a/src/lipase/imagej/ijimageengine.py +++ b/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))