#BIOP #PSF #Deconvolution

ComputePSF.groovy

/**
 * This script takes a 3D input a stack with beads and a stack of the same size
 * with single points placed at the center of detected beads / spots
 * this is the second step to compute a clean PSF with distillation.
 * 
 * Example for this script:
 * - Image = https://zenodo.org/records/14203207/files/30x1000_FS-025_FW-60_AC-170.tif
 * - Bead Image = the result of the DetectBeads script
 * 
 * Link to useful threads:
 * https://forum.image.sc/t/clij-deconvolution-regularization-fiji-crash/55155/20
 * Intro to deconvolution by Brian Northan:
 * https://forum.image.sc/t/getting-started-with-deconvolution-in-imagej/32916/1
 * 
 * This script is a modified version of the scripts created by Brian in the forum post above
 * */

#@OpService ops

#@ConvertService cs
#@output ImgPlus thresholded
#@output ImgPlus psf

#@ImagePlus beads_imp
#@ImagePlus points_imp

// To work on files, you can use the parameters below instead
//#@File imageFile
//def beads_imp = IJ.openImage(imageFile.getAbsolutePath());
//def points_imp = IJ.openImage(imageFile.getAbsolutePath()+"_points.tif");

// reference
// (note. there used to be good references on the Huygens web site, but that seems only accessible by password now)
// http://www.leica-microsystems.com/science-lab/measuring-the-3d-sted-psf-with-a-new-type-of-fluorescent-beads/

// assumptions
// 1. assuming a 3D (x,yz) dataset
// 2. assuming the image contains sub-resolution beads
//    note.  It is expected some beads may "clump" into larger objects,
//    the script may still work under these conditions, but the cleaner the
//    the bead image the better the PSF.
// 3. We use integer locations for bead centroids.  In reality beads are centered
//    at sub-resolution locations.  However as long as we have enough beads the error
//    in position should "even out".


def points = ImagePlusAdapter.wrapImgPlus(points_imp) //cs.convert(points_imp, ImgPlus.class)
def beads = ImagePlusAdapter.wrapImgPlus(beads_imp) //cs.convert(beads_imp, ImgPlus.class)


// get dimensions of bead image
def xSize = beads.dimension(0)
def ySize = beads.dimension(1)
def zSize = beads.dimension(2)

// call connected components to label each connected region
def labeling = ops.labeling().cca(points, StructuringElement.FOUR_CONNECTED)

// get the index image (each object will have a unique gray level)
def labelingIndex = labeling.getIndexImg()

// get the collection of regions and loop through them
def regions = new LabelRegions(labeling)

def randomAccess = points.randomAccess()
for (region in regions) {

	// get the center of mass of the region
	center = region.getCenterOfMass()

	// place a point at the bead centroid
	randomAccess.setPosition(new long[]{(long)center.getFloatPosition(0), (long)center.getFloatPosition(1), (long)center.getFloatPosition(2)})
	randomAccess.get().setReal(255.0)
}

// convert images to float
def beads32 = ops.convert().float32(beads)
def points32 = ops.convert().float32(points)

// use "PSF Distilling" (reverse deconvolution) to solve for PSF
def image = beads32
def kernel = points32
def borderSize = null
def obfInput = null
def obfKernel = null
def outType = new FloatType()
def fftType = null
def maxIterations = 10
def nonCirculant = false
def accelerate = true // The default is False.
def psf = ops.create().img(image, outType)
ops.deconvolve().richardsonLucy(psf, beads32, points32, borderSize, obfInput, obfKernel, outType, fftType, maxIterations, nonCirculant, accelerate)

import net.imglib2.img.display.imagej.ImageJFunctions

def output_psf = ImageJFunctions.wrap(psf, "title") // wraps the ImgPlus as an ImagePlus
output_psf.show()

IJ.run(output_psf, "Properties...", "channels=1 slices="+beads_imp.getNSlices()+" frames=1 pixel_width=1.0000 pixel_height=1.0000 voxel_depth=1.0000");
IJ.run(output_psf, "Enhance Contrast", "saturated=0.35");
output_psf.setRoi(116,121,31,25);
output_psf.setSlice(65);
IJ.run(output_psf, "Enhance Contrast", "saturated=0.35");

// To save the output file automatically, uncomment the lines below
//IJ.saveAs(output_psf, "Tiff", imageFile.getAbsolutePath()+"_psf.tif");
//IJ.run("Close All", "");


import ij.IJ
import ij.ImageJ
import net.imagej.ImgPlus

import net.imglib2.algorithm.labeling.ConnectedComponents.StructuringElement
import net.imglib2.roi.labeling.LabelRegions
import net.imglib2.type.numeric.real.FloatType
import net.imglib2.img.ImagePlusAdapter // https://gist.github.com/GenevieveBuckley/460d0abc7c1b13eee983187b955330ba#file-convert_imageplus_to_imgplus-py