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))