use of gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.
the class GaussianFit method fitMultiple.
/**
* Fits a 2D Gaussian to the given data. Fits all the specified peaks.
* <p>
* Data must be arranged in yx block order, i.e. height rows of width.
* <p>
* Note: The fit coordinates should be offset by 0.5 if the input data represents pixels
*
* @param data
* @param width
* @param height
* @param maxIndices
* Indices of the data to fit
* @param estimatedHeights
* Estimated heights for the peaks (input from smoothed data)
* @return Array containing the fitted curve data: The first value is the Background. The remaining values are
* Amplitude, PosX, PosY, StdDevX, StdDevY for each fitted peak. If elliptical fitting is performed
* the values are Amplitude, Angle, PosX, PosY, StdDevX, StdDevY for each fitted peak
* <p>
* Null if no fit is possible.
*/
private double[] fitMultiple(float[] data, int width, int height, int[] maxIndices, double[] estimatedHeights) {
if (data == null || data.length != width * height)
return null;
if (maxIndices == null || maxIndices.length == 0)
return null;
Gaussian2DFitter gf = createGaussianFitter(false);
this.fitResult = gf.fit(Utils.toDouble(data), width, height, maxIndices, estimatedHeights);
if (fitResult.getStatus() == FitStatus.OK) {
chiSquared = fitResult.getError();
double[] params = fitResult.getParameters();
convertParameters(params);
return params;
}
return null;
}
use of gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.
the class GaussianFit method fit.
/**
* Fits a single 2D Gaussian to the data. The fit is initialised at the highest
* value and then optimised.
* <p>
* Data must be arranged in yx block order, i.e. height rows of width.
* <p>
* The angle parameter is only returned if using elliptical Gaussian fitting.
* <p>
* Note: The fit coordinates should be offset by 0.5 if the input data represents pixels
*
* @return Array containing the fitted curve data: Background, Amplitude, PosX, PosY, StdDevX, StdDevY, Angle. Null
* if no fit is possible.
*/
public double[] fit(float[] data, int width, int height) {
if (data == null || data.length != width * height)
return null;
// Get the limits
float max = Float.MIN_VALUE;
int maxIndex = -1;
for (int i = data.length; i-- > 0; ) {
float f = data[i];
if (max < f) {
max = f;
maxIndex = i;
}
}
if (maxIndex < 0) {
return null;
}
Gaussian2DFitter gf = createGaussianFitter(false);
FitResult fitResult = gf.fit(Utils.toDouble(data), width, height, new int[] { maxIndex });
if (fitResult.getStatus() == FitStatus.OK) {
chiSquared = fitResult.getError();
double[] params = fitResult.getParameters();
// Check bounds
if (params[3] < 0 || params[3] >= width || params[4] < 0 || params[4] >= height)
return null;
// Re-arrange order for backwards compatibility with old code.
return new double[] { params[0], params[1], params[3], params[4], Gaussian2DFitter.fwhm2sd(params[5]), Gaussian2DFitter.fwhm2sd(params[6]), params[2] };
}
return null;
}
use of gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.
the class SpotAnalysis method updateCurrentSlice.
private void updateCurrentSlice(int slice) {
if (slice != currentSlice) {
currentSlice = slice;
double signal = getSignal(slice);
double noise = smoothSd[slice - 1];
currentLabel.setText(String.format("Frame %d: Signal = %s, SNR = %s", slice, Utils.rounded(signal, 4), Utils.rounded(signal / noise, 3)));
drawProfiles();
// Fit the PSF using a Gaussian
float[] data2 = (float[]) rawImp.getImageStack().getProcessor(slice).getPixels();
double[] data = Utils.toDouble(data2);
FitConfiguration fitConfiguration = new FitConfiguration();
fitConfiguration.setFitFunction(FitFunction.FIXED);
fitConfiguration.setBackgroundFitting(true);
fitConfiguration.setSignalStrength(0);
fitConfiguration.setCoordinateShift(rawImp.getWidth() / 4.0f);
fitConfiguration.setComputeResiduals(false);
fitConfiguration.setComputeDeviations(false);
Gaussian2DFitter gf = new Gaussian2DFitter(fitConfiguration);
double[] params = new double[7];
double psfWidth = Double.parseDouble(widthTextField.getText());
params[Gaussian2DFunction.BACKGROUND] = smoothMean[slice - 1];
params[Gaussian2DFunction.SIGNAL] = (gain * signal);
params[Gaussian2DFunction.X_POSITION] = rawImp.getWidth() / 2.0f;
params[Gaussian2DFunction.Y_POSITION] = rawImp.getHeight() / 2.0f;
params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
FitResult fitResult = gf.fit(data, rawImp.getWidth(), rawImp.getHeight(), 1, params, new boolean[1]);
if (fitResult.getStatus() == FitStatus.OK) {
params = fitResult.getParameters();
final double spotSignal = params[Gaussian2DFunction.SIGNAL] / gain;
rawFittedLabel.setText(String.format("Raw fit: Signal = %s, SNR = %s", Utils.rounded(spotSignal, 4), Utils.rounded(spotSignal / noise, 3)));
ImageROIPainter.addRoi(rawImp, slice, new PointRoi(params[Gaussian2DFunction.X_POSITION], params[Gaussian2DFunction.Y_POSITION]));
} else {
rawFittedLabel.setText("");
rawImp.setOverlay(null);
}
// Fit the PSF using a Gaussian
if (blurImp == null)
return;
data2 = (float[]) blurImp.getImageStack().getProcessor(slice).getPixels();
data = Utils.toDouble(data2);
params = new double[7];
//float psfWidth = Float.parseFloat(widthTextField.getText());
params[Gaussian2DFunction.BACKGROUND] = (float) smoothMean[slice - 1];
params[Gaussian2DFunction.SIGNAL] = (float) (gain * signal);
params[Gaussian2DFunction.X_POSITION] = rawImp.getWidth() / 2.0f;
params[Gaussian2DFunction.Y_POSITION] = rawImp.getHeight() / 2.0f;
params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
fitResult = gf.fit(data, rawImp.getWidth(), rawImp.getHeight(), 1, params, new boolean[1]);
if (fitResult.getStatus() == FitStatus.OK) {
params = fitResult.getParameters();
final double spotSignal = params[Gaussian2DFunction.SIGNAL] / gain;
blurFittedLabel.setText(String.format("Blur fit: Signal = %s, SNR = %s", Utils.rounded(spotSignal, 4), Utils.rounded(spotSignal / noise, 3)));
ImageROIPainter.addRoi(blurImp, slice, new PointRoi(params[Gaussian2DFunction.X_POSITION], params[Gaussian2DFunction.Y_POSITION]));
} else {
blurFittedLabel.setText("");
blurImp.setOverlay(null);
}
}
}
use of gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.
the class GaussianFit method runFinal.
/**
* Perform fitting using the chosen maxima. Update the overlay if successful.
*
* @param ip
* The input image
*/
private void runFinal(ImageProcessor ip) {
ip.reset();
Rectangle bounds = ip.getRoi();
// Crop to the ROI
float[] data = ImageConverter.getData(ip);
int width = bounds.width;
int height = bounds.height;
// Sort the maxima
float[] smoothData = data;
if (getSmooth() > 0) {
// Smoothing destructively modifies the data so create a copy
smoothData = Arrays.copyOf(data, width * height);
AverageFilter filter = new AverageFilter();
//filter.blockAverage(smoothData, width, height, smooth);
if (smooth <= border)
filter.stripedBlockAverageInternal(smoothData, width, height, (float) smooth);
else
filter.stripedBlockAverage(smoothData, width, height, (float) smooth);
}
Sort.sort(maxIndices, smoothData);
// Show the candidate peaks
if (maxIndices.length > 0) {
String message = String.format("Identified %d peaks", maxIndices.length);
if (isLogProgress()) {
IJ.log(message);
for (int index : maxIndices) {
IJ.log(String.format(" %.2f @ [%d,%d]", data[index], bounds.x + index % width, bounds.y + index / width));
}
}
// Check whether to run if the number of peaks is large
if (maxIndices.length > 10) {
GenericDialog gd = new GenericDialog("Warning");
gd.addMessage(message + "\nDo you want to fit?");
gd.showDialog();
if (gd.wasCanceled())
return;
}
} else {
IJ.log("No maxima identified");
return;
}
results = new IJTablePeakResults(showDeviations, imp.getTitle() + " [" + imp.getCurrentSlice() + "]");
results.begin();
// Perform the Gaussian fit
long ellapsed = 0;
if (!singleFit) {
if (isLogProgress())
IJ.log("Combined fit");
// Estimate height from smoothed data
double[] estimatedHeights = new double[maxIndices.length];
for (int i = 0; i < estimatedHeights.length; i++) estimatedHeights[i] = smoothData[maxIndices[i]];
FitConfiguration config = new FitConfiguration();
setupPeakFiltering(config);
long time = System.nanoTime();
double[] params = fitMultiple(data, width, height, maxIndices, estimatedHeights);
ellapsed = System.nanoTime() - time;
if (params != null) {
// Copy all the valid parameters into a new array
double[] validParams = new double[params.length];
int c = 0;
int validPeaks = 0;
validParams[c++] = params[0];
double[] initialParams = convertParameters(fitResult.getInitialParameters());
double[] paramsDev = convertParameters(fitResult.getParameterStdDev());
Rectangle regionBounds = new Rectangle();
int[] xpoints = new int[maxIndices.length];
int[] ypoints = new int[maxIndices.length];
int nMaxima = 0;
for (int i = 1, n = 0; i < params.length; i += 6, n++) {
int y = maxIndices[n] / width;
int x = maxIndices[n] % width;
// Check the peak is a good fit
if (filterResults && config.validatePeak(n, initialParams, params) != FitStatus.OK)
continue;
if (showFit) {
// Copy the valid parameters
validPeaks++;
for (int ii = i, j = 0; j < 6; ii++, j++) validParams[c++] = params[ii];
}
double[] peakParams = extractParams(params, i);
double[] peakParamsDev = extractParams(paramsDev, i);
addResult(bounds, regionBounds, data, peakParams, peakParamsDev, nMaxima, x, y, data[maxIndices[n]]);
// Add fit result to the overlay - Coords are updated with the region offsets in addResult
double xf = peakParams[3];
double yf = peakParams[4];
xpoints[nMaxima] = (int) (xf + 0.5);
ypoints[nMaxima] = (int) (yf + 0.5);
nMaxima++;
}
setOverlay(nMaxima, xpoints, ypoints);
// Draw the fit
if (showFit && validPeaks != 0) {
double[] pixels = new double[data.length];
EllipticalGaussian2DFunction f = new EllipticalGaussian2DFunction(validPeaks, width, height);
invertParameters(validParams);
f.initialise(validParams);
for (int x = 0; x < pixels.length; x++) pixels[x] = f.eval(x);
FloatProcessor fp = new FloatProcessor(width, height, pixels);
// Insert into a full size image
FloatProcessor fp2 = new FloatProcessor(ip.getWidth(), ip.getHeight());
fp2.insert(fp, bounds.x, bounds.y);
Utils.display(TITLE, fp2);
}
} else {
if (isLogProgress()) {
IJ.log("Failed to fit " + Utils.pleural(maxIndices.length, "peak") + getReason(fitResult));
}
imp.setOverlay(null);
}
} else {
if (isLogProgress())
IJ.log("Individual fit");
int nMaxima = 0;
int[] xpoints = new int[maxIndices.length];
int[] ypoints = new int[maxIndices.length];
// Extract each peak and fit individually
ImageExtractor ie = new ImageExtractor(data, width, height);
float[] region = null;
Gaussian2DFitter gf = createGaussianFitter(filterResults);
for (int n = 0; n < maxIndices.length; n++) {
int y = maxIndices[n] / width;
int x = maxIndices[n] % width;
long time = System.nanoTime();
Rectangle regionBounds = ie.getBoxRegionBounds(x, y, singleRegionSize);
region = ie.crop(regionBounds, region);
int newIndex = (y - regionBounds.y) * regionBounds.width + x - regionBounds.x;
if (isLogProgress()) {
IJ.log("Fitting peak " + (n + 1));
}
double[] peakParams = fitSingle(gf, region, regionBounds.width, regionBounds.height, newIndex, smoothData[maxIndices[n]]);
ellapsed += System.nanoTime() - time;
// Output fit result
if (peakParams != null) {
double[] peakParamsDev = null;
if (showDeviations) {
peakParamsDev = convertParameters(fitResult.getParameterStdDev());
}
addResult(bounds, regionBounds, data, peakParams, peakParamsDev, n, x, y, data[maxIndices[n]]);
// Add fit result to the overlay - Coords are updated with the region offsets in addResult
double xf = peakParams[3];
double yf = peakParams[4];
xpoints[nMaxima] = (int) (xf + 0.5);
ypoints[nMaxima] = (int) (yf + 0.5);
nMaxima++;
} else {
if (isLogProgress()) {
IJ.log("Failed to fit peak " + (n + 1) + getReason(fitResult));
}
}
}
// Update the overlay
if (nMaxima > 0)
setOverlay(nMaxima, xpoints, ypoints);
else
imp.setOverlay(null);
}
results.end();
if (isLogProgress())
IJ.log("Time = " + (ellapsed / 1000000.0) + "ms");
}
use of gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.
the class GaussianFit method createGaussianFitter.
private Gaussian2DFitter createGaussianFitter(boolean simpleFiltering) {
FitConfiguration config = new FitConfiguration();
config.setMaxIterations(getMaxIterations());
config.setSignificantDigits(getSignificantDigits());
config.setDelta(getDelta());
config.setInitialPeakStdDev(getInitialPeakStdDev());
config.setComputeDeviations(showDeviations);
config.setDuplicateDistance(0);
// Set-up peak filtering only for single fitting
config.setDisableSimpleFilter(!simpleFiltering);
setupPeakFiltering(config);
if (isLogProgress()) {
config.setLog(new IJLogger());
}
if (getFitCriteria() >= 0 && getFitCriteria() < FitCriteria.values().length) {
config.setFitCriteria(FitCriteria.values()[getFitCriteria()]);
} else {
config.setFitCriteria(FitCriteria.LEAST_SQUARED_ERROR);
}
if (getFitFunction() >= 0 && getFitFunction() < FitFunction.values().length) {
config.setFitFunction(FitFunction.values()[getFitFunction()]);
} else {
config.setFitFunction(FitFunction.CIRCULAR);
}
config.setBackgroundFitting(fitBackground);
return new Gaussian2DFitter(config);
}
Aggregations