[Run Fiji Local Thickness in Parallel] this script can crop blob-like objects in an image using `Find Maxima` and runs `Local Thickness` on each cropped image. #fiji #imagej #imagesc #localthickness #BIOP

Local_Thickness_Parallel.groovy

/**
 * This script computes local thickness on binary images which (optionally) contain an ROI to be analyzed
 * The goal is to run local thickness in parallel to detect blob like structures that are quite far apart from each other. 
 * To achieve this, we scale down the original image, blur it, detect maxima
 * Each maxima is then used as the center of a box to crop the original image
 * Each cropped image is processed in parallel for local thickness
 * Reply to an Image.sc forum post: https://forum.image.sc/t/local-thickness-run-time/42143/
 * Licensed under BSD 3 license
 * Copyright Olivier Burri 2020
 * Copyright 2020 Olivier Burri
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#@ImagePlus original (label="Image to process")
#@Double scale (label="Resizing scale", value=0.25)
#@Integer tolerance (label="Max finder tolerance", value=150)
#@Integer crop_box (label="Crop Box Size", value=15)

// Get image information for scaling image later
def w = original.getWidth()
def h = original.getHeight()

// Capture ROI to limit maximum detection inside it
def roi = original.getRoi()

// Get rid of the ROI before cropping to avoid the scaled image being cropped.
original.killRoi()
// Scale the ROI if there is one using RoiScaler
def scaled_roi = roi != null ?  scaled_roi = RoiScaler.scale( roi, scale, scale, false ) : null

// Scale the image using the resize function. Need final width and height for this function
def scaled = original.resize( Math.round( w * scale ) as int, Math.round( h * scale ) as int, "bilinear" )

// Set the ROI on the rescaled image before using find maxima
if (scaled_roi != null) scaled.setRoi( scaled_roi )

// Not sure it is needed, but ensure that the image does not have an inverting LUT
// Also make sure that the spots are white on black
if ( scaled.isInvertedLut() ) scaled.getProcessor().invertLut()

// Blur to avoid counting broken spots twice
IJ.run(scaled, "Gaussian Blur...", "sigma=2")

// Maximum finder which returns a java polygon
def final maxima = new MaximumFinder().getMaxima( scaled.getProcessor(), tolerance, true )

// After locating the points, we can create the crop boxes for each point on the original image
def rois = (0..maxima.npoints-1).collect { i ->
	x = ( (double) maxima.xpoints[i] + 0.5 ) / scale - crop_box/2 
	y = ( (double) maxima.ypoints[i] + 0.5 ) / scale - crop_box/2
	print(".")
	return new Roi( x, y, crop_box, crop_box )
}

// We can crop the images by giving a list of Rois to the crop method. 
def images = original.crop( rois as Roi[] )

// Make sure that the ROI is reset on the original image
original.setRoi(roi)

// Parallel processing: Run Local Thickness
def locthk = null
GParsPool.withPool(10) {
locthk = images.collectParallel{ image ->
	println( "Computing local thickness for image: " + image )
	return LocalThickness.getLocalThickness( image )
}
}

// Stack all images together and get the StackStatistics histogram
locthk = ImagesToStack.run( locthk as ImagePlus[] )

//locthk.show()
// Get the histogram as a list so we can remove elements
def histogram = new StackStatistics( locthk, 256, 1.0, 256.0 ).histogram()

// Make a montage to see what is happening
def stack_size = locthk.getStackSize()
def cols = Math.round( Math.sqrt(stack_size ) )
def rows = cols*2 < stack_size ? cols+1 : cols

def montage = new MontageMaker().makeMontage2( locthk, cols as int, cols as int, 1.0, 1, stack_size, 1, 0, false)

montage.show()
IJ.run(montage, "Fire", "")
IJ.resetMinAndMax(montage)

def plot = new Plot( "Scack Histogram of " + original.getTitle(), "Bins", "Counts" )
def bins = (1..histogram.size()-1).collect{ it } as double[]
plot.add( "bar", bins, histogram )
plot.show()

// Convenience class to compute local thicknesses via scripting without showing the image
class LocalThickness {
	public static ImagePlus cleanUpLocalThickness(ImagePlus imp) {
		Clean_Up_Local_Thickness cult = new Clean_Up_Local_Thickness();
		cult.setup("", imp);
		cult.runSilent = true;
		cult.run(imp.getProcessor());
		return cult.getResultImage();
	}

	public static ImagePlus distanceRidge(ImagePlus imp) {
		Distance_Ridge dr = new Distance_Ridge();
		dr.setup("", imp);
		dr.runSilent = true;
		dr.run(imp.getProcessor());
		return dr.getResultImage();
	}

	public static ImagePlus distanceMap(ImagePlus imp) {
		EDT_S1D edt = new EDT_S1D();
		edt.setup("", imp);
		edt.runSilent = true;
		edt.showOptions = false;
		edt.run(imp.getProcessor());
		return edt.getResultImage();
	}

	public static ImagePlus localThickness( ImagePlus imp ) {
		Local_Thickness_Parallel lt = new Local_Thickness_Parallel();
		lt.setup("", imp);
		lt.runSilent = true;
		lt.run(imp.getProcessor());
		return lt.getResultImage();
	}

	public static ImagePlus getLocalThickness( ImagePlus imp ) {
		return localThickness( distanceRidge( distanceMap( imp ) ) );
	}
}

import ij.ImagePlus
import ij.IJ
import ij.gui.Plot
import ij.gui.Roi
import ij.plugin.filter.MaximumFinder
import ij.plugin.RoiScaler
import ij.plugin.ImagesToStack
import ij.plugin.MontageMaker
import ij.process.StackStatistics
import sc.fiji.localThickness.*
import groovyx.gpars.*