How to use LeastSquaresBuilder in jython imagej script

Dear All,

I would like to use LeastSquaresBuilder from org.apache.commons.maths3 for fitting a spot on a image with a 2D gaussian. I know how to fit a 2D gaussian with scipy but for a Imagej script I cannot use the scipy library and need to use the LeastSquaresBuilder. If you have a example with the use of LeastSquaresBuilder in jython script, it will be wonderful…
Please let me know

I have no time to give you a full example. But this link might help you showing examples in Java.
The conversion to Jython should’nt be hard:

Also in this forum you find a lot of useful Jython examples and link:

You would first have to put the jars for the math package from apache and its dependencies in your fiji jars folder.
Might not be so straightforward :sweat_smile:

Scijava which is shipped with Fiji has a couple of classes for function fitting on images, it might be a bit of work but at least you can test it right away.

An import of one of those classes would look like that in jython
from net.imglib2.algorithm.localization import FitFunction

I havent used it myself, but some people like @tinevez or @imagejan might have, they might be rather busy at the moment though…

For general (jython) scripting documentation in Fiji I would advise

Hello @philgi (thank you @LThomas for the heads-up).
I made this a long time ago for Mohamed El Beheiry.
I think it is a good idea to all this directly in Jython rather than implementing LSQ loop in Jython.
Tell us if you need help.

Thank you @LThomas and @tinevez to point me out these classes from Imglib2. I still do not know well enough this new library and am still using the ImageJ1 library.
Here is my code using these classes

def fit2DGaussian(imp_) :
center = array([imp_.width / 2, imp_.height / 2], 'l')
cpoint = Point(center)
data = LocalizationUtils.gatherObservationData(ImagePlusAdapter.wrapShort(imp_), cpoint, center)
params = MLGaussianEstimator(2.0, 2).initializeFit(cpoint, data)

LevenbergMarquardtSolver.solve(data.X, params, data.I, Gaussian(), 1e-3, 1e-1, 300)
return params

with these classes:

from net.imglib2 import Point
from  net.imglib2.algorithm.localization import Gaussian, LevenbergMarquardtSolver, LocalizationUtils, MLGaussianEstimator, Observation
from net.imglib2.img import ImagePlusAdapter

but it does not work. Do you have any idea why ?
thanks by advance


Let’s make it work like we did before.
Can you send me your script and a test image, and I will play with it on my side?
It worked for us before.

Hello @philgi

Ok, here is a version of the script that works. I put comments and explanations in the script itself:

from math import sqrt

from ij import  IJ, ImagePlus,WindowManager 

from net.imglib2 import Point as Pt
from net.imglib2.util import Util
from net.imglib2.algorithm.localization import PeakFitter, Gaussian, LevenbergMarquardtSolver, LocalizationUtils, MLGaussianEstimator, Observation
from net.imglib2.img import ImagePlusAdapter

from jarray import zeros,array

def fit2DGaussian( img, peaks ) :

	# The core class is PeakFitter.
	# It needs to know what solver to use (here LM), what function to fit (here,
	# a plain symmetric Gaussian), and an estimator that can generate
	# some estimates of the fit parameters as initial values for the fit (here, the 
	# MLGaussianEstimator, configured with a typical sigma of 2 and on 2D).
	fitter = PeakFitter( img, peaks, LevenbergMarquardtSolver(), Gaussian(), MLGaussianEstimator( 2., 2 ) )
	print( fitter.toString() )
	if not fitter.checkInput() or not fitter.process():
		print( 'Error with the peak fitter: %s.' % fitter.getErrorMessage() )

	return fitter.getResult()


imp = WindowManager.getCurrentImage()

# Wrap into an ImgLib2 image.
img = ImagePlusAdapter.wrap( imp )

# You need to provide the initial guess of the position of peaks.
# It is a list of 'Localizable', for instance 'Point's.
center = array([imp.width / 2, imp.height / 2], 'l')
cpoint = Pt(center)
peaks = [ cpoint ] # Here we have just one point.
# This list can be provided automatically, e.g. by a LoG detector.

# Fit by a symmetric 2D Gaussian.
results = fit2DGaussian( img, peaks )

# The results are returned as a dictionary of the peak objects as 
# keys, and with the fit output as values.
# In our case we fit A × exp( -b × ∑ (xᵢ - x₀ᵢ)² ) and the fit output
# is a 4-double array made of x, y, A and b.

for peak in results:
	print( '\nFound a peak around %s. Fit results:' % Util.printCoordinates( peak ) )
	params = results[ peak ]
	print( ' - X position: %9.3f' % params[0] )
	print( ' - Y position: %9.3f' % params[1] )
	print( ' - Amplitude:  %9.3f' % params[2] )
	print( ' - Sigma:      %9.3f' % ( 1. / sqrt(2. * params[3]) ) )

When I run it on the image you sent:
Gauss.tif (9.6 KB)

I get the following output:

Started at Thu Mar 26 23:27:53 CET 2020
PeakFitter configured to:
 - fit a Gaussian function A × exp( -b × ∑ (xᵢ - x₀ᵢ)² )
 - on 1 peaks
 - in image FloatImagePlus [49x49]
 - using Maximum-likelihood estimator for symetric gaussian peaks
 - and Levenberg-Marquardt least-square curve fitting algorithm
 - allocating 8 threads.

Found a peak around (24.0, 24.0). Fit results:
 - X position:    24.992
 - Y position:    25.004
 - Amplitude:      0.195
 - Sigma:          2.804

Tell me if this is what you need.


Hi @tinevez

as usual, you did a good job. It is working as I want.
Just the equation used for the fit: A × exp( -b × ∑ (xᵢ - x₀ᵢ)² ) does not propose an offset because the background is not necessary equal to 0. So I suppose that it is possible to create its own function like this one
f(x,y) = Aexp(-(x-x_m)^2/(2s_x^2))exp(-(y-y_m)^2/(2s_y^2))+offset
with the interface FunctionFitter. That’s right.

Let me know…
And thank you for the code it’s perfect…

True, exactly.
But I imposed such functions return the gradient and the hessian to be used with many solvers, and that can be a pain to implement.
What about background removal instead?

This is the interface you want to implement:

ok thank you. You are right, it is better to remove the background.

It just wanted to know how to proceed if i need to implement a new function in the future…

Your code is working as it is and I will use it …


If someone has some time, I think it would be a great single-particle tracking scripting use case to put in the documentation.
Just not sure where this shoud go, the wiki or the imageJ tuto repository ? Or does it exist already ?
Maybe @imagejan, @ctrueden could tell :wink: