Richardson Lucy TV deconvolution

Hi,

The issue I raised is because I read the article on ImageJ list about the deconvolution. The Richardson Lucy TV deconvolution algorithm was introduced by Brian Northan (@bnorthan) on the last ImageJ conference. It seems to be a promising method for image deconvolution from the open source. Is this method able to do deconvolution for images acquired from confocal microscopes?

The codes are on the Github now and will continue to be maintained, though, I could not find it on Fiji itself or its update site. As SteveWest said on ImageJ list “Is there any implementation of it in Fiji/ImageJ for us to try out?”.

All the best,
Chin

It is contained in Fiji, but not accessible via a menu command. But you can find it using Plugins > Utilities > Find Ops…:

It is only visible in developer view, though.

To use it in a script, you can do (adapted from the presentation slide):

# @OpService ops
# @Dataset image
# @Dataset psf
# @OUTPUT Dataset deconvolved

deconvolved=ops.deconvolve().richardsonLucyTV(image, psf, 10, 1.0);
3 Likes

Hi Jan Eglinger,

Thanks for the answer.
I have no idea which language is supposed to be selected when I try to put them in a “script”. Please correct me if I used wrong method. I put these lines into a new macro window, then I chose Python as its “language” and run.

It is the key point if I did use a right script to deconvolve the image at the beginning. Please let me know if you have any suggestion. I appreciate your information.

Best,
chin

Yes, I forgot to state that the script adapted from the presentation is written in Python. But apart from the comment syntax, that doesn’t make a difference, since it is calling the Java API of ImageJ-Ops.

In Groovy, the script would look like this:

// @OpService ops
// @Dataset image
// @Dataset psf
// @OUTPUT Dataset deconvolved

deconvolved=ops.deconvolve().richardsonLucyTV(image, psf, 10, 1.0)
1 Like

Hi @chin

If you want to use the method with a confocal microscope you need to have a PSF of your microscope. You could image and crop beads (measured PSF). Or you could generate a theoretical PSF.

The ‘Ops’ project is “foundational”, it is meant to provide a common infrastructure for higher level mage processing apps. It is also meant to be usable in mutliple platforms (like ImageJ, KNIME, etc). However (currently) it is better suited for programmers and script writers.

I recommend also trying Deconvolution Lab. It looks like they are releasing a new version soon. I will be trying this out myself when it is available.

It would be great if you could share images and Point Spread Function(s). I think we should perhaps start a new Wiki with deconvolution examples using ops, Deconvolution Lab, and other plugins. We can then get a good idea of the usability and performance of the various implementations.

I am happy to take ownership of the wiki page and work on example scripts (using Ops), and document how to process images with other plugins. I just need images that people are able to share.

Brian

4 Likes

Hi Brian,

Thank you for your suggestion.

I would like to try DeconvolutionLab or its new version, though, it needs some parameters setting and use a PSF generator. I have been using AutoQuant for deconvolution for a long time, and my colleagues want to try some open source method. But we did not have experience to use the DeconvolutionLab.

I uploaded some PSF images and regular images for your test.
https://file.io/hcHOPg
https://file.io/KVT4xR
Please let me know if you need further material. I appreciate your any experience.

Best,
chin

@bnorthan I added ops RL-TV to my convolution/deconvolution demo javascript.
Interesting difference in results vs. old iterative 3D deconvolution plugin
RL-TV has better noise squashing but +ve image intensity offset at high spatial frequencies…

Took me a while to figure out the @parameters stuff need to go at the very top,
that I need to input Img not ImagePlus,
and also the syntax of the ops command, and how many and which parameters to use.

How do I use the acceleration and boundary conditions and other fancy parameters?

3 Likes

Thanks for sharing this great illustrative script, @chalkie666!


From the Script parameters documentation:

Any parameters after the first non-parameter line will not be recognized.

If that isn’t clear enough, please feel free to make it more visible, e.g. by different formatting or an info box.

@bnorthan @chin @imagejan Linked here is a minimal-ish python script to test ops RL-TV deconv on some ImageJ sample data, using Bob D’s Diffraction PSF 3D plugin to make a theoretical (optimistic) point spread function image.

@bnorthan I was surprised that I have to do the following object and type conversion stuff?
Is it as expected?
I thought using ops that kind of thing happens magically in the background?
Am I doing it wrong? I get the feeling I am doing it a hard way…
I was expecting to simply pass the name of the open image windows to the ops…?

# get ImagePlus from IJ1 image window for channel 1 of the sample image. 
IJ.selectWindow("PSF")
psfimp = IJ.getImage()
IJ.selectWindow("C1-confocal-series.tif")
imageimp = IJ.getImage()

# wrap IJ1 ImagePlus into IJ2 Img for use in ops
psfimg = ImageJFunctions.wrap(psfimp)
imageimg = ImageJFunctions.wrap(imageimp)
3 Likes

Well, ideally, you’d just use @Parameters to get these images in the correct form (@ImgPlus or @ImagePlus as needed) and you’d avoid any call to ij.IJ. And to display an image, you’d declare it as an @OUTPUT and it will be displayed when the script finishes.


Of course, in your case, the interactive workflow requires that you open and display images during script execution, so harvesting the parameters in the script header won’t work here.

However, you should be able to use a Service for that, e.g. for opening an image:

// @IOService io
// @OUTPUT ImgPlus img

img = io.open("/path/to/your/image");

(This however didn’t work for URLs in my tests, see this issue. There’s a discussion about remote locations here, maybe @gab1one or @dietzc can comment.)

or for displaying the image within the script:

// @IOService io
// @UIService ui

img = io.open("/path/to/your/image");

ui.show(img);

// do some more things

Your demo script with the interactive explanations would be an ideal candidate for a tutorial in the style of an IPython or Beaker notebook, like the ones in imagej-tutorials @etadobson is currently working on.

2 Likes

More help and ideas on the tutorials would certainly be awesome. Right now, @panovr is the one taking lead on them, with help from @etarena as her time allows (although she is quite busy doing ImageJ workshops lately!).

Hi @chin

Sorry I took so long to reply. I got busy with a couple of other things. I tried to upload your files and I get a ‘404 not found’ message. Is it possible the link expired because I waited to long?? Let me know.

Thanks

Brian

@chalkie666

Here is a groovy script that shows how to use deconvolution ops with acceleration parameter. The call is really annoying because the signature is too long, you have to pass in a bunch of nulls (for boundary conditions and such) before getting to the acceleration parameter. I need to figure out how to make better signatures.

I have been distracted with other projects the last few months, but I am going to try and make an effort to get back to this, and put a little bit in each day, even if it is only 20 minutes. Next week I will work out an example with boundary conditions.

If anybody has a chance to run the script let me know if there are issues, or if there are ways it can be improved.

Thanks

// @OpService ops
// @UIService ui
// @ConvertService convert
// @DatasetService data

import ij.IJ;
import net.imglib2.FinalDimensions
import net.imglib2.type.numeric.real.FloatType
import net.imagej.ImgPlus
import net.imglib2.img.Img
import net.imglib2.img.display.imagej.ImageJFunctions

def plotProfile() {
	IJ.resetMinAndMax();
	IJ.makeLine(0, 32, 511, 32);
	IJ.run("Plot Profile");
}

blank=ops.create().img(new FinalDimensions(512,512),new FloatType())
noisy=ops.create().img(new FinalDimensions(512,512),new FloatType())

formula = "50 * (Math.sin(2*Math.PI*0.1*Math.pow(3,p[0]/149.8)*p[0]/149.8 )+1)+1"
exponentialChirp = ops.image().equation(blank, formula)
ui.show(exponentialChirp)

plotProfile();

// generate diffraction based psf
IJ.run("Diffraction PSF 3D", "index=1.520 numerical=1.42 wavelength=510 "
+ "longitudinal=0 image=10 slice=200 width,=512 height,=512 depth,=1 "
+ "normalization=[Sum of pixel values = 1] title=PSF");

// get the psf
psf=IJ.getImage()

// convert to imglib img
psf=ImageJFunctions.wrapFloat(psf)

// convolve chirp and psf
convolved=ops.filter().convolve(exponentialChirp, psf)
ui.show("convolved", convolved)

plotProfile()

// add noise
noisy=ops.filter().addPoissonNoise(noisy,convolved)
ui.show("noisy", noisy)

plotProfile()

// deconvolve
deconvolved=ops.deconvolve().richardsonLucy(noisy, psf,10)
ui.show("deconvolved", deconvolved)
plotProfile()

deconvolved_accelerated=ops.deconvolve().richardsonLucy(noisy, psf, null, null, null, null, null, 10, false, true);
ui.show("deconvolved accelerated", deconvolved_accelerated)
plotProfile()

deconvolved_rltv=ops.deconvolve().richardsonLucyTV(noisy, psf, 10, 0.01);
ui.show("deconvolved rltv", deconvolved_rltv)
plotProfile()

deconvolved_rltv_accelerated=ops.deconvolve().richardsonLucyTV(noisy, psf, null, null, null, null, null, 10, false, true, 0.01);
ui.show("deconvolved rltv accelerated", deconvolved_rltv_accelerated)
plotProfile()
4 Likes

Hi Brain,

Yes, the link expired but I am out of office right now. I will upload those images to another site next Monday. Sorry for the missing files.

Best,
chin

1 Like

Hi Brian,

You could download images include for beads and cells now.
I uploaded two sets of images acquired by using confocal microscopes.

http://www.mediafire.com/file/5kh9uc2hwzhytmt
http://www.mediafire.com/file/iob4baxahvjg8tp

Best,
chin

@chin

Thanks. I got them. Are these confocal images??

Edit: Just re-read your original post. They are confocal images.

Hi @chin

Here is a script that should work on your images. It creates a Gaussian to use as the PSF. You may need to experiment with different Gaussian sizes.

Is anyone aware of a better way to calculate a confocal PSF in ImageJ?? I believe PSF generator is only widefield.

Keep in mind ops is still in beta. I recommend also trying other plugins and verifying the results match. I will probably do that myself this week by testing against ‘Deconvolution Lab’.

# @OpService ops
# @UIService ui
# @ImgPlus img
# @Integer sxy
# @Integer sz
# @Integer numIterations
# @OUTPUT ImgPlus deconvolved

# convert to float (TODO: make sure deconvolution op works on other types)
img_float=ops.convert().float32(img)

# create and show the gaussian kernel
psf=ops.create().kernelGauss([sxy, sxy, sz])
ui.show(psf)

# deconvolve
deconvolved=ops.deconvolve().richardsonLucyTV(img_float, psf, numIterations, 0.01);
3 Likes

@bnorthan I will have a go at testing this groovy script as soon as I can. Looks nice! Cheers!

@chin @bnorthan in a very rough sense the true closed pinhole confocal PSF (not 1Airy unit… That’s not truly confocal but often used for good optical sectioning), is the square of the WF psf. That’s what I used in my demo script to show you still need to Deconvolve confocal images. 1 airy unit pinhole confocal PSF might be somewhere inbetween WF and true confocal. It’s possible to calculate theoretically. But I dont know how. Anyway its better to use a real measured and noise cleaned Psf.

1 Like

@ctrueden @bnorthan @hinerm
Just so you know, I noticed that running RLTV from ops finder, then selecting the PSF image in the drop down, causes IJ to flog the processors and the hard disk, then eventually crash IJ.

http://imagej.net/Op_Finder says
"The play button essentially automates the process of turning an Op into a script: optional parameters are discarded and required parameters are annotated. Because of this, Ops with arcane or unusual parameters may fail to run because the framework does not know how to provide them."
So I guess that’s not a surprise :wink:
Let me know if i can be any help testing fixes for that.