added a function to remove lighting non uniformity

src/lipase/catalog.py

@@ -1,5 +1,6 @@
 
 import json
+from imageengine import IHyperStack, IImageEngine, PixelType
 
 
 class DacMetadata(object):
@@ -74,15 +75,15 @@ class Sequence(object):
     def get_root_path(self):
         return '/'.join(self.micro_manager_metadata_file_path.split('/')[:-1])
 
-    def get_image_file_path(self, channel_index, frame_index, z_index=0):
+    def get_image_file_path(self, channel_index, frame_index, slice_index=0):
         '''
         :param int channel_index:
         :param int frame_index:
-        :param int z_index:
+        :param int slice_index:
         '''
         assert frame_index < self.num_frames
         assert channel_index < self.num_channels
-        frame_info = self.mmm['FrameKey-%d-%d-%d' % (frame_index, channel_index, z_index)]
+        frame_info = self.mmm['FrameKey-%d-%d-%d' % (frame_index, channel_index, slice_index)]
         rel_file_path = frame_info['FileName']
         return self.get_root_path() + '/' + rel_file_path
 
@@ -127,6 +128,33 @@ class Sequence(object):
         white_sequence = seqid_to_white[self.sequence_id]
         return self.catalog.sequences[white_sequence]
 
+    def as_hyperstack(self, selected_channel_ids=None, selected_frames=None, selected_slices=None):
+        """ returns a subset of the sequence as an hyperstack
+
+        :param list(str) selected_channel_ids:
+        :return IHyperStack: the resulting hyperstack
+        """
+        if selected_frames is None:
+            selected_frames = range(self.num_frames)
+
+        if selected_slices is None:
+            selected_slices = range(self.num_slices)
+
+        if selected_channel_ids is None:
+            selected_channel_indices = range(self.num_channels)
+        else:
+            selected_channel_indices = [self.get_channel_index(channel_id) for channel_id in selected_channel_ids]
+
+        pixel_type = {8: PixelType.U8, 16: PixelType.U16}[self.num_bits_per_pixels]
+
+        stack = IImageEngine.get_instance().create_hyperstack(width=self.width, height=self.height, num_slices=len(selected_slices), num_frames=len(selected_frames), num_channels=len(selected_channel_indices), pixel_type=pixel_type)
+        for frame_index in selected_frames:
+            for slice_index in selected_slices:
+                for channel_index in selected_channel_indices:
+                    image = IImageEngine.get_instance().load_image(self.get_image_file_path(channel_index=channel_index, frame_index=frame_index, slice_index=slice_index))
+                    stack.set_image(image=image, channel_index=channel_index, frame_index=frame_index, slice_index=slice_index)
+        return stack
+
 
 class ImageCatalog(object):
     def __init__(self, raw_images_root):

src/lipase/imageengine.py

@@ -7,6 +7,28 @@ import abc
 
 ABC = abc.ABCMeta('ABC', (object,), {})
 
+class Aabb(object):
+    """Axis-Aligned bounding box"""
+
+    def __init__(self, x_min, y_min, x_max, y_max):
+        self.x_min = x_min
+        self.y_min = y_min
+        self.x_max = x_max
+        self.y_max = y_max
+        
+
+class PixelType:
+    U8 = 1
+    U16 = 2
+    F32 = 3
+
+    PIXEL_TYPE_TO_NUM_BITS = {U8: 8,
+        U16: 16,
+        F32:32}
+
+    @staticmethod
+    def get_num_bits(pixel_type):
+        return PixelType.PIXEL_TYPE_TO_NUM_BITS[ pixel_type ]
 
 class IImage(ABC):
 
@@ -19,6 +41,25 @@ class IImage(ABC):
         pass
 
 
+class IHyperStack(ABC):
+
+    @abc.abstractmethod
+    def width(self):
+        pass
+
+    @abc.abstractmethod
+    def height(self):
+        pass
+
+    @abc.abstractmethod
+    def set_image(self, image, frame_index=0, slice_index=0, channel_index=0):
+        """
+        :param IImage image:
+        """
+        pass
+
+
+
 class IImageEngine(ABC):
     """
     The active image processor.
@@ -46,6 +87,18 @@ class IImageEngine(ABC):
         assert IImageEngine.__instance is not None
         return IImageEngine.__instance
 
+    @abc.abstractmethod
+    def create_hyperstack(self, width, height, num_channels, num_slices, num_frames, pixel_type):
+        """
+        """
+
+    @abc.abstractmethod
+    def load_image(self, src_image_file_path):
+        """
+        :param str src_image_file_path: eg './white_estimate.tiff'
+        :return IImage:
+        """
+
     @abc.abstractmethod
     def save_as_tiff(self, image, out_file_path):
         """

src/lipase/imagej/ijimageengine.py

@@ -1,7 +1,7 @@
 """Image processing using imagej.
 """
 
-from lipase.imageengine import IImage, IImageEngine
+from ..imageengine import IImage, IHyperStack, IImageEngine, PixelType
 from ij import IJ  # pylint: disable=import-error
 from ij.plugin import ImageCalculator  # pylint: disable=import-error
 from ijopencv.ij import ImagePlusMatConverter  # pylint: disable=import-error
@@ -22,10 +22,29 @@ class IJImage(IImage):
         self.ij_image = ij_image
 
     def width(self):
-        raise NotImplementedError
+        return self.ij_image.width
 
     def height(self):
-        raise NotImplementedError
+        return self.ij_image.height
+
+
+class IJHyperStack(IHyperStack):
+
+    def __init__(self, width, height, num_channels, num_slices, num_frames, pixel_type):
+        """
+        """
+        stack_name = ''
+        self.hyperstack = IJ.createHyperStack(stack_name, width, height, num_channels, num_slices, num_frames, PixelType.get_num_bits(pixel_type))
+
+    def set_image(self, image, frame_index=0, slice_index=0, channel_index=0):
+        self.hyperstack.setPositionWithoutUpdate(channel_index + 1, slice_index + 1, frame_index + 1)
+        self.hyperstack.setProcessor(image.ij_image.getProcessor())
+
+    def width(self):
+        return self.hyperstack.width
+
+    def height(self):
+        return self.hyperstack.height
         
 
 def imagej_run_image_command(image, command, options):
@@ -123,9 +142,21 @@ class IJImageEngine(IImageEngine):
     def missing(self):
         pass
 
+    def create_hyperstack(self, width, height, num_channels, num_slices, num_frames, pixel_type):
+        return IJHyperStack(width, height, num_channels, num_slices, num_frames, pixel_type)
+
+    def load_image(self, src_image_file_path):
+        image_plus = IJ.openImage(src_image_file_path)
+        return IJImage(self, image_plus)
+
     def save_as_tiff(self, image, out_file_path):
         IJ.saveAsTiff(image.ij_image, out_file_path)
 
+    def divide(self, numerator_image, denominator_image):
+        ic = ImageCalculator()
+        result_image = ic.run("Divide create float", numerator_image.ij_image, denominator_image.ij_image)
+        return IJImage(self, result_image)
+
     def compute_max(self, images_file_path):
         """Computes for each pixel position the maximum at this position in all input images.
 

src/lipase/preprocessing.py

@@ -1,7 +1,7 @@
 """preprocessing of synchrotron images based on telemosToolbox."""
 
 from catalog import ImageCatalog, Sequence
-from imageengine import IImageEngine
+from imageengine import IImageEngine, PixelType
 
 class WhiteEstimator(object):
 
@@ -122,8 +122,8 @@ class WhiteEstimator(object):
             for channel_id in channel_ids:
                 channel_index = sequence.get_channel_index(channel_id)
                 for frame_index in range(sequence.num_frames):
-                    for z_index in range(sequence.num_slices):
+                    for slice_index in range(sequence.num_slices):
-                        images_file_path.append(sequence.get_image_file_path(channel_index, frame_index, z_index))
+                        images_file_path.append(sequence.get_image_file_path(channel_index, frame_index, slice_index))
 
         white_estimate = IImageEngine.get_instance().compute_max(images_file_path)
         
@@ -160,6 +160,64 @@ class InteractiveWhiteEstimator(WhiteEstimator):
         #         white_estimator_settings.ask
 
 
+def correct_non_uniform_lighting(non_uniform_sequence, channel_id, white_estimator=WhiteEstimator(open_size=75, close_size=75, average_size=75)):
+    """ correction of non uniform lighting
+    :param Sequence non_uniform_sequence: the input sequence, which is expected to have a non uniform lighting
+    :param str channel_id:
+    :param WhiteEstimator white_estimator: the method used to estimate the light image
+    :return IHyperStack: the hyperstack of images after correction
+    """
+    white_estimate = white_estimator.estimate_white([non_uniform_sequence], [channel_id])
+    uniform_stack = IImageEngine.get_instance().create_hyperstack(width=white_estimate.width(), height=white_estimate.height(), num_slices=1, num_frames=non_uniform_sequence.num_frames, num_channels=1, pixel_type=PixelType.F32)
+
+    for frame_index in range(non_uniform_sequence.num_frames):
+        non_uniform_frame = IImageEngine.get_instance().load_image(non_uniform_sequence.get_image_file_path(channel_index=0, frame_index=frame_index, slice_index=0))
+        uniform_frame = IImageEngine.get_instance().divide(non_uniform_frame, white_estimate)
+        uniform_stack.set_image(frame_index=frame_index, image=uniform_frame)
+    return uniform_stack
+
+# def get_sequence_as_hyperstack(sequence, selected_channel_ids=None, selected_frames=None, selected_slices=None):
+#     """ returns a subset of the sequence as an hyperstack
+#     :param IHyperStack non_uniform_sequence: the input sequence
+#     :param list(str) selected_channel_ids:
+#     :return IHyperStack: the resulting hyperstack
+#     """
+#     if selected_frames is None:
+#         selected_frames = range(sequence.num_frames)
+
+#     if selected_slices is None:
+#         selected_slices = range(sequence.num_slices)
+
+#     if selected_channel_ids is None:
+#         selected_channel_indices = range(sequence.num_channels)
+#     else:
+#         selected_channel_indices = [sequence.get_channel_index(channel_id) for channel_id in selected_channel_ids]
+
+#     pixel_type = {8: PixelType.U8, 16: PixelType.U16}[sequence.num_bits_per_pixels]
+
+#     stack = IImageEngine.get_instance().create_hyperstack(width=sequence.width, height=sequence.height, num_slices=len(selected_slices), num_frames=len(selected_frames), num_channels=len(selected_channel_indices), pixel_type=pixel_type)
+#     for frame_index in selected_frames:
+#         for slice_index in selected_slices:
+#             for channel_index in selected_channel_indices:
+#                 image = IImageEngine.get_instance().load_image(sequence.get_image_file_path(channel_index=channel_index, frame_index=frame_index, slice_index=slice_index))
+#                 stack.set_image(image=image, channel_index=channel_index, frame_index=frame_index, slice_index=slice_index)
+#     return stack
+
+# def correct_non_uniform_lighting(non_uniform_sequence, channel_id, white_estimator=WhiteEstimator(open_size=75, close_size=75, average_size=75)):
+#     """ correction of non uniform lighting
+#     :param IHyperStack non_uniform_sequence: the input sequence, which is expected to have a non uniform lighting
+#     :param WhiteEstimator white_estimator: the method used to estimate the light image
+#     :return IHyperStack: the hyperstack of images after correction
+#     """
+#     white_estimate = white_estimator.estimate_white([non_uniform_sequence], [channel_id])
+#     uniform_stack = IImageEngine.get_instance().create_hyperstack(width=white_estimate.width(), height=white_estimate.height(), num_slices=1, num_frames=non_uniform_sequence.num_frames, num_channels=1, pixel_type=PixelType.F32)
+
+#     for frame_index in range(non_uniform_sequence.num_frames):
+#         non_uniform_frame = IImageEngine.get_instance().load_image(non_uniform_sequence.get_image_file_path(channel_index=0, frame_index=frame_index, slice_index=0))
+#         uniform_frame = IImageEngine.get_instance().divide(non_uniform_frame, white_estimate)
+#         uniform_stack.set_image(frame_index=frame_index, image=uniform_frame)
+#     return uniform_stack
+
 def find_white_reference_image(white_estimate, white_z_stack):
     """Find the white reference Image among a z stack of reference image the most correlated to the white image estimated from a serie of acquisition.
 

src/lipase/traps_remover.py

@@ -0,0 +1,32 @@
+
+from catalog import ImageCatalog, Sequence
+# from imageengine import Toto
+from imageengine import IImageEngine
+from imageengine import PixelType
+from preprocessing import WhiteEstimator
+
+def remove_traps(sequence, channel_id, trap_aabb):
+    """Remove the traps in the input sequence
+
+    :param Sequence sequence:
+    :param str channel_id:
+    :param Aabb trap_aabb:
+    """
+    white_estimator = WhiteEstimator(open_size=75, close_size=75, average_size=75)
+    white_estimate = white_estimator.estimate_white([sequence], [channel_id])
+    print(white_estimate)
+
+    uniform_stack = IImageEngine.get_instance().create_hyperstack(width=white_estimate.width(), height=white_estimate.height(), num_slices=1, num_frames=sequence.num_frames, num_channels=1, pixel_type=PixelType.F32)
+
+    for frame_index in range(sequence.num_frames):
+        non_uniform_frame = IImageEngine.get_instance().load_image(sequence.get_image_file_path(channel_index=0, frame_index=frame_index, slice_index=0))
+        uniform_frame = IImageEngine.get_instance().divide(non_uniform_frame, white_estimate)
+        uniform_stack.set_image(frame_index=frame_index, image=uniform_frame)
+
+    #num_expected_traps = 13  # 13 traps are completely visible in the 
+    #assert num_traps = num_expected_traps
+    
+    #non_uniform_stack = sequence.as_stack()
+    #uniform_stack = IImageEngine.get_instance().divide(non_uniform_stack, white_estimate)
+    # IImageEngine.get_instance().save_as_tiff(white_estimate, './white_estimate.tiff')
+    

tests/test0001.py

@@ -7,10 +7,11 @@
 
 import unittest  # unittest2 doesn't exist in fiji
 import sys
-from lipase.catalog import ImageCatalog, Sequence
+from lipase.imageengine import IImageEngine, PixelType, Aabb
-from lipase.imageengine import IImageEngine
 from lipase.imagej.ijimageengine import IJImageEngine, IJImage
-from lipase.preprocessing import WhiteEstimator
+from lipase.preprocessing import WhiteEstimator, correct_non_uniform_lighting
+from lipase.traps_remover import remove_traps
+from lipase.catalog import ImageCatalog, Sequence
 
 class TestLipase(unittest.TestCase):
 
@@ -28,8 +29,6 @@ class TestLipase(unittest.TestCase):
         print("uninitializing TestLipase instance")
         self.catalog = None
 
-    
-
     def test_estimate_white(self):
         sequence = self.catalog.sequences['res_soleil2018/GGH/GGH_2018_cin2_phiG_I_327_vis_-40_1/Pos2']
         white_estimator = WhiteEstimator(open_size=75, close_size=75, average_size=75)
@@ -40,6 +39,18 @@ class TestLipase(unittest.TestCase):
         # assert False, "hellooooo"
         print('end of test_estimate_white')
 
+    def test_uniform_lighting_correction(self):
+        non_uniform_sequence = self.catalog.sequences['res_soleil2018/GGH/GGH_2018_cin2_phiG_I_327_vis_-40_1/Pos0']
+        uniform_sequence = correct_non_uniform_lighting(non_uniform_sequence, 'DM300_nofilter_vis', white_estimator=WhiteEstimator(open_size=75, close_size=75, average_size=75))
+
+    # def test_traps_removal(self):
+    #     sequence = self.catalog.sequences['res_soleil2018/GGH/GGH_2018_cin2_phiG_I_327_vis_-40_1/Pos0']
+    #     x_min = 423
+    #     x_max = 553
+    #     y_min = 419
+    #     y_max = 533
+    #     trap_aabb = Aabb(x_min, y_min, x_max, y_max)
+    #     remove_traps(sequence, 'DM300_nofilter_vis', trap_aabb)
 
 def run_script():