lipase/preprocessing.py

"""preprocessing of synchrotron images based on telemosToolbox."""

from ij import IJ
from ij.plugin import ImageCalculator
from catalog import ImageCatalog, Sequence


def imagej_run_image_command(image, command, options):
    """performs the given imagej command on the given image

    :param ImagePlus image:
    :param str command: imagej command (eg "Gray Morphology")
    :param str options: imagej options (eg "radius=1 type=square operator=open")

    wrapper around IJ.run (https://imagej.nih.gov/ij/developer/api/ij/IJ.html#run-ij.ImagePlus-java.lang.String-java.lang.String-) which raises an exception on error
    """
    IJ.run(image, command, options)
    error_message = IJ.getErrorMessage()
    if error_message is not None:
        raise Exception('The command "%s" with options "%s" failed because of the following error : %s' % (command, options, error_message))


def compute_max(images_file_path):
    """Computes for each pixel position the maximum at this position in all input images.

    :param list(str) images_file_path:
    :rtype ImagePlus:
    """
    assert len(images_file_path) > 1
    max_image = IJ.openImage(images_file_path[0])
    print('max_image', max_image)

    for image_file_path in images_file_path[2:-1]:
        other_image = IJ.openImage(image_file_path)
        print('other_image', other_image)
        ic = ImageCalculator()
        ic.run("max", max_image, other_image)
    print('max_image', max_image)
    return max_image


def replace_border(image, band_width):
    """Overwrites the outer band of the image by duplicating the value of the pixels that touch this outer band

    :param ImagePlus image:
    :param int band_width: width of the outer band, in pixels
    :rtype ImagePlus:

    adaptation for the following code from matlab telemosToolbx's estimatedwhiteFluoImageTelemos function
        outer = white_estimate;
        white_estimate = white_estimate(4:(end-3),4:(end-3));
        outer(1:3,4:(end-3))=repmat(white_estimate(1,:),3,1);  # top
        outer((end-2):end,4:(end-3))=repmat(white_estimate(end,:),3,1); # bottom
        outer(:,1:3)=repmat(outer(:,4),1,3); # left
        outer(:,(end-2):end)=repmat(outer(:,end-3),1,3); # right
        white_estimate = outer;
        clear outer;
    """
    raise NotImplementedError()
    return image


def perform_gray_morphology(image, operator, structuring_element_shape, structuring_element_radius):
    """
    :param ImagePlus image:
    :param str operator: eg 'open'
    :param str structuring_element_shape: eg 'square'
    :param int structuring_element_radius:
    """
    assert operator not in ['fast open', 'fast erode'], "as of 13/09/2019, fast operators such as 'fast erode' seem broken in fiji (the resulting image is not at all similar to their slow equivalent)"
    
    processor = image.getProcessor()
    convert_to_byte = False
    if processor.getBitDepth() != 8:
        convert_to_byte = True
        min_value = processor.getMin()
        max_value = processor.getMax()
        print("warning: downgrading image to byte as imagej's Gray Morphology processing doesn't handle 16bit images (range=[%d, %d])" % (min_value, max_value))
        do_scaling = True
        image.setProcessor(processor.convertToByte(do_scaling))
    print("before gray morphology")
    assert structuring_element_radius < 11, "the radius of the structuring element is too big (%d); using it with Fiji's 'Gray Morphology' tool would result in very long computations." % structuring_element_radius
    imagej_run_image_command(image, "Gray Morphology", "radius=%d type=%s operator=%s" % (structuring_element_radius, structuring_element_shape, operator))
    print("after gray morphology")


class WhiteEstimator(object):

    def __init__(self, open_size, close_size, average_size):
        self.open_size = open_size
        self.close_size = close_size
        self.average_size = average_size

    def _remove_particles(self, white_estimate):
        perform_gray_morphology(white_estimate, operator='open', structuring_element_shape='square', structuring_element_radius=(self.open_size + 1) / 2)
        perform_gray_morphology(white_estimate, operator='close', structuring_element_shape='square', structuring_element_radius=(self.close_size + 1) / 2)

    def estimate_white(self, sequences, channel_ids, dark=None):
        """Estimation of the white fluorescence image shape of synchrotron light from experimental images of Telemos microscope.

        :param list(Sequence) sequences: the sequences to consider
        :param list(str) channel_ids: the channels to consider, eg ['DM300_327-353_fluo']
        :param WhiteEstimatorSettings white_estimator_settings: 

        :param ImagePlus or None dark:
        :rtype ImagePlus:
        Code adapted from matlab telemosToolbx's estimatedwhiteFluoImageTelemos function

        % this function is specific of Telemos images (DISCO line SOLEIL)
            % acquired using micromanager software linked to imageJ

        %% input

        %  nothing : interactive function

        %% output

        %  nothing : the final reference fluorescence image is saved on disk

        %% principe
        %  NB: MICROMANAGER save images in a structured file folder architecture :
        %
        % root folder : name given by the user
        %
        %   subfolders roi(n)_tile1 : n folders for each selected roi
        %           or pos(n) : n folders for each selected position
        %
        %   display_and_comments.txt : file describing the channels acquired
        %
        % in each subfolder roi(n)_tile1 or Pos(n) :
        %       images files with name img_00000000(n)_DM300_327-353_00(p).tif
        %               n = number of time
        %               p = z focal plane
        %               DM300_327-353 = name of the channel
        %       metadata.txt : files describing all metadata associated with the
        %       acquisition : x, y, z position, camera settings .... etc.
        %
        % (Nb: depending on the date of acquisition, an extension fluo is found or
        % not in the name of fluorescence image. This extension is given in
        % contrast to the visible image.)
        %
        % IMAGES : fluorescence images are considered
        %
        %
        % IMPORTANT:
        %       it is expected that in the folders at least several images cover the whole field of view
        %       illuminated by the synchrtron light
        %       IF IT IS NOT THE CASE, IT WON'T WORK
        %
        % BASIC: the estimated white image is the MAX of the selected images
        %       when interactive mode is selected, only first time and first z are
        %       proposed
        %       for automatic computing, all images are taken into account
        %
        % WHAT IS DONE:
        %        0 - the user chosse the channels that are considered  in the estimation of the
        % white fluorescence image
        %
        %       1 - spurious pixels (three lines and columns around the image) are not taken
        % in filtering and replaced by line and column number 4 and end-3
        %
        %       2 - Filtering
        % the estimated dark image is filtered :
        %        opening to remove  white regions brighter than the low frequency signal that should correspond to  synchrotron light shape (should be small)
        %        closing to remove black regions that remains (should be small to avoid synchrotron light shape deformation)
        %        average filtering

        % the image computed should be used  to
        %   - to find the best white z plane from the reference white stack
        %   acquired for ex. using the so called "Matthieu lame"%
        %   - correct images for intensities

        %% use
        %       estimatedwhiteFluoImageTelemos

        %% Comments
        %   adapted from proposals 20161050 and 20171187


        %% Author
        % MF Devaux
        % INRA BIA
        % PVPP

        %% date
        % 5 octobre 2017:
        % 23 mars 2018
        % 3 septembre 2018 : comments and general case
        % 16 avril 2019: name of function and default filtering values
        % 29 avril 2019 : new comments and spurious pixels
        % 27 mai 2019 : offset dark
        %     4 juin 2019 : replace exist by isfolder or isfile
        % 14 juin : close figure at the end
        """
        images_file_path = []
        for sequence in sequences:
            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):
                        images_file_path.append(sequence.get_image_file_path(channel_index, frame_index, z_index))

        white_estimate = compute_max(images_file_path)
        
        # modify spurious pixels on the side of the images
        try:
            replace_border(white_estimate, 3)
        except NotImplementedError as error:
            print('warning: replace_outer_frame is not implemented yet')

        self._remove_particles(white_estimate)
        
        return white_estimate


class InteractiveWhiteEstimator(WhiteEstimator):

    def __init__(self, open_size, close_size, average_size):
        WhiteEstimator.__init__(self, open_size, close_size, average_size)

    def _remove_particles(self, white):
        """shows the image to a user so that he can visually check that the particles have been removed correctly in the given image
        """
        raise NotImplementedError()
        # part_rem_settings_are_good = False
        # while part_rem_settings_are_good == False:
        #     # perform opening to remove white particles
        #     IJ.run( white_estimate, "Gray Morphology", "radius=1 type=square operator=open" )
        #     # run("Gray Morphology", "radius=1 type=square operator=open");
        #     # IJ.run( input_image_plus_copy, "Skeletonize (2D/3D)", "" )
        #     # perform 
        #     if white_estimator_settings.particles_are_removed():
        #         part_rem_settings_are_good = true
        #     else
        #         white_estimator_settings.ask


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.

    :param Image white_estimate:
    :param Sequence white_z_stack:
    % this function is specific of Telemos images (DISCO line SOLEIL)
        % acquired using micromanager software linked to imageJ

    %% input
    %  nothing : interactive function

    %% output
    %    zmax : focal plane

    %% principe
    % images were acquired for a given roi which position is found in
    % metadata.txt file
    % white and dark  images were recorded for the full field of view of the DISCO
    % TELEMOS camera
    %
    % correlation between a estimated  white fluorescence image estalished from
    % actual acquisitions of a given sample and all z plane acquired for the reference TELEMOS white
    % image (Matthieu slide or any fluorescent homogeneous image) :

    % The  estimated  white fluorescence image  is generallly obtained bt using function whiteFluoImageTelemosestimation
    % This is not compulsary as any homogenous sample image hat can roughly show the shape of illumination can be used to find
    % the white reference image
    %
    %
    % the z plane for which the maximum correlation is observed between estimated white and reference white images is retained.
    % the white image is then offsetted (its corresponding Dark is subtracted) and copied in the subfolder <WhiteReference> of the sample
    % rootfolder to show that it has been specifically selected for the considered experiment
    %The matlab corrcoeff function is used
    %
    % correlation coefficients are saved in a file called
    % 'corr.txt' in the subfolder 'WhiteReference'
    %
    %
    % expected input folder hierarchy:
    %
    %      >sampleFolder
    %           > PosFolders or RoiFolders
    %           > WhiteEstimate
    %       >darkFolder
    %       >darkFolder.smooth
    %       >whiteFolder
    %           > darkFolderforWhite
    %           > darkFolderforWhite.smooth
    %       >whiteFolder.smooth
    %
    %
    % expected output folder hierarchy:
    %
    %      >sampleFolder
    %           > PosFolders or RoiFolders
    %           > WhiteEstimate
    %           > WhiteReference
    %                   > white after offset
    %       >darkFolder
    %       >darkFolder.smooth
    %       >whiteFolder
    %           > darkFolderforWhite
    %           > darkFolderforWhite.smooth
    %       >whiteFolder.smooth





    %% use
    % [zmax]=findWhiteReferenceImageTelemos

    %% Comments
    %   written for  proposal 20161050
    %   adapted for proposal 20170043


    %% Author
    % MF Devaux
    % INRA BIA
    % PVPP

    %% date
    % 5 octobre 2017:
    % 15 decembre 2017 : adapted to take into account the roi known from
    % metadata
    % 27 fevrier 2018 : comments details
    % 4 septembre 2018: comments and check and naming of folders
    % 12 mars 2019 : save white reference with the same size as white estimate
    % 16 avril 2019 : include diagonals to check the relevance of white
    % reference
    % 20 mai 2019 : track folder
    % 27 mai 2019 : offset dark
    %     4 juin 2019 : replace exist by isfolder or isfile
    """
    raise NotImplementedError()


def test_preprocessing():
    """Test preprocessing."""
    catalog = ImageCatalog('/Users/graffy/ownCloud/ipr/lipase/raw-images')
    sequence = 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)
    white_estimate = white_estimator.estimate_white([sequence], ['DM300_327-353_fluo'])
    find_white_reference_image(white_estimate, sequence.get_white())


if __name__ == '__main__':
    test_preprocessing()