[Parallel Histogram Matching for Stacks] #Groovy #Fiji #ImageJ

Histogram_Matching_Parallel.groovy

//@File (label="Directory with file(s)",  style="directory" ) theDir
//@Integer (Label="Channel To Correct") ch_correct
//@Boolean (label="Use Normalization Parameters Below") is_normalize
//@Integer (Label="Minimum") the_min
//@Integer (Label="Maximum") the_max

/**
 * Image stack equalizer for image folder
 * 
 * The purpose of this script is to use the histogram of the first slice of the first image
 * in a folder as a reference and use Histogram Matching on all other slices in all other images
 * 
 * This was originally used to fix a flickering brightfield acquisition to make all images the same
 * brightness and contrast.
 * 
 * It normalizes the first image between the Minimum and Maximum value provided by the user before starting
 * if the checkbox is set.
 * Author: Olivier Burri, BioImage Analyst, BiImaging And Optics Platform (BIOP) 
 * 
 * Licensed under CC-BY-SA https://creativecommons.org/licenses/by-sa/4.0/
 */
 
import histogram2.HistogramMatcher
import ij.plugin.ChannelSplitter
import ij.plugin.RGBStackMerge
import ij.ImagePlus
import ij.IJ

import groovy.io.FileType
import groovyx.gpars.GParsExecutorsPool


import loci.plugins.BF;


// List all the files and keep only the ND ones
def allFilesList = []
// Get all the files in all the subfolders
theDir.eachFile(FileType.FILES) { file ->
	if(file.name.endsWith(".tif"))
		allFilesList << file
}


def savedir = new File(theDir.getAbsolutePath()+'/Equalized/')
savedir.mkdir()


def ref_hist

allFilesList.eachWithIndex{ file, i ->
	def imp = IJ.openImage(file.getAbsolutePath())
	
		print("\nProcessing Image "+imp.getTitle())
		String imageTitle = file.getName().substring(0,file.getName().size()-4)
		def image = imp;
		def all_images = []
		if(imp.getNChannels() > 1) {
			all_images = ChannelSplitter.split(imp)
			image = all_images[ch_correct-1] 
		}
			
		if (i == 0) {
			def ref_ip = image.getStack().getProcessor(1).convertToFloat().duplicate()
			if (is_normalize) { 
				def stats = ref_ip.getStatistics()
				ref_ip.subtract(stats.min)
				ref_ip.multiply(1/(stats.max-stats.min))
				ref_ip.multiply(the_max-the_min)
				ref_ip.add(the_min)
				
			}
			ref_ip = ref_ip.convertToShort(false)
			ref_hist = ref_ip.getHistogram()
		}
		def matchedImage= matchSliceHistograms(ref_hist, image)
		def newImage = matchedImage
		if(imp.getNChannels() > 1) {
			all_images[ch_correct-1] = matchedImage	
			newImage = RGBStackMerge.mergeChannels(all_images, false)
			newImage.setDisplayMode(imp.getDisplayMode())
			newImage.setLuts(imp.getLuts())
		} 
		

		// Make sure it is still the right hyperstack dimensions
		
		IJ.saveAs(newImage, "tif", savedir.getAbsolutePath()+"/"+imageTitle+".tif") 
		newImage.close()
}

return;

ImagePlus matchSliceHistograms(int[] refhist, ImagePlus imp) {
	def stack = imp.getStack()
	def final_imp = IJ.createHyperStack("Nobody cares",imp.getWidth(), imp.getHeight(), imp.getNChannels(), imp.getNSlices(), imp.getNFrames(), imp.getBitDepth())
	def matcher = new HistogramMatcher()
	def cores = Runtime.getRuntime().availableProcessors()
	GParsExecutorsPool.withPool(cores) {
		(1..stack.getSize()).eachParallel{
			print("\n\tProcessing Slice "+it)
			def hist_old = stack.getProcessor(it).getHistogram()
			def newHist = matcher.matchHistograms(hist_old, refhist)
			def ip = stack.getProcessor(it).duplicate()
			ip.applyTable(newHist)
			stack.setProcessor(ip, it)
		}
	}
	final_imp.setStack(stack)
	return final_imp
}