use of gdsc.smlm.engine.FitEngine in project GDSC-SMLM by aherbert.
the class PeakFit method createFitEngine.
/**
* Creates a fitting engine using the current configuration.
*
* @param numberOfThreads
* @param queue
* @param queueSize
* @return The fiting engine
*/
public FitEngine createFitEngine(int numberOfThreads, FitQueue queue, int queueSize) {
PeakResults r = results;
if (results.numberOfOutputs() == 1)
// Reduce to single object for speed
r = results.toArray()[0];
// Update the configuration
updateFitConfiguration(config);
FitEngine engine = new FitEngine(config, r, numberOfThreads, queue, queueSize);
// Write settings out to the IJ log
if (resultsSettings.logProgress) {
IJ.log("-=-=-=-");
IJ.log("Peak Fit");
IJ.log("-=-=-=-");
Utils.log("Initial Peak SD = %s,%s", Utils.rounded(fitConfig.getInitialPeakStdDev0()), Utils.rounded(fitConfig.getInitialPeakStdDev1()));
SpotFilter spotFilter = engine.getSpotFilter();
IJ.log("Spot Filter = " + spotFilter.getDescription());
int w = 2 * engine.getFitting() + 1;
Utils.log("Fit window = %d x %d", w, w);
if (!fitConfig.isDisableSimpleFilter()) {
IJ.log("Coordinate shift = " + Utils.rounded(config.getFitConfiguration().getCoordinateShift()));
IJ.log("Signal strength = " + Utils.rounded(fitConfig.getSignalStrength()));
}
if (fitConfig.isDirectFilter())
IJ.log("Smart filter = " + fitConfig.getSmartFilter().getDescription());
if (extraOptions)
IJ.log("Noise = " + Utils.rounded(fitConfig.getNoise()));
IJ.log("Width factor = " + Utils.rounded(fitConfig.getWidthFactor()));
IJ.log("-=-=-=-");
}
return engine;
}
use of gdsc.smlm.engine.FitEngine in project GDSC-SMLM by aherbert.
the class PeakFit method run.
/**
* Locate the peaks in the configured image source. Results are saved to the configured output.
* <p>
* This must be called after initialisation with an image source. Note that each call to this method must be
* preceded with initialisation to prepare the image and output options.
*/
public void run() {
if (source == null)
return;
int totalFrames = source.getFrames();
ImageStack stack = null;
if (showProcessedFrames)
stack = new ImageStack(bounds.width, bounds.height);
// Do not crop the region from the source if the bounds match the source dimensions
Rectangle cropBounds = (bounds.x == 0 && bounds.y == 0 && bounds.width == source.getWidth() && bounds.height == source.getHeight()) ? null : bounds;
// Use the FitEngine to allow multi-threading.
FitEngine engine = createFitEngine(getNumberOfThreads(totalFrames));
final int step = Utils.getProgressInterval(totalFrames);
runTime = System.nanoTime();
boolean shutdown = false;
int slice = 0;
while (!shutdown) {
// Noise can optionally be estimated from the entire frame
float[] data = (ignoreBoundsForNoise) ? source.next() : source.next(cropBounds);
if (data == null)
break;
if (++slice % step == 0) {
if (Utils.showStatus("Slice: " + slice + " / " + totalFrames))
IJ.showProgress(slice, totalFrames);
}
float noise = Float.NaN;
if (ignoreBoundsForNoise) {
noise = FitWorker.estimateNoise(data, source.getWidth(), source.getHeight(), config.getNoiseMethod());
// Crop the data to the region
data = ImageConverter.getData(data, source.getWidth(), source.getHeight(), bounds, null);
}
if (showProcessedFrames) {
stack.addSlice(String.format("Frame %d - %d", source.getStartFrameNumber(), source.getEndFrameNumber()), data);
}
// Get the frame number from the source to allow for interlaced and aggregated data
engine.run(createJob(source.getStartFrameNumber(), source.getEndFrameNumber(), data, bounds, noise));
if (escapePressed())
shutdown = true;
}
engine.end(shutdown);
time = engine.getTime();
runTime = System.nanoTime() - runTime;
if (showProcessedFrames)
Utils.display("Processed frames", stack);
showResults();
source.close();
}
use of gdsc.smlm.engine.FitEngine in project GDSC-SMLM by aherbert.
the class PSFEstimator method calculateStatistics.
private boolean calculateStatistics(PeakFit fitter, double[] params, double[] params_dev) {
debug(" Fitting PSF");
swapStatistics();
// Create the fit engine using the PeakFit plugin
FitConfiguration fitConfig = config.getFitConfiguration();
fitConfig.setInitialAngle((float) params[0]);
fitConfig.setInitialPeakStdDev0((float) params[1]);
fitConfig.setInitialPeakStdDev1((float) params[2]);
ImageStack stack = imp.getImageStack();
Rectangle roi = stack.getProcessor(1).getRoi();
ImageSource source = new IJImageSource(imp);
// Allow interlaced data by wrapping the image source
if (interlacedData) {
source = new InterlacedImageSource(source, dataStart, dataBlock, dataSkip);
}
// Allow frame aggregation by wrapping the image source
if (integrateFrames > 1) {
source = new AggregatedImageSource(source, integrateFrames);
}
fitter.initialiseImage(source, roi, true);
fitter.addPeakResults(this);
fitter.initialiseFitting();
FitEngine engine = fitter.createFitEngine();
// Use random slices
int[] slices = new int[stack.getSize()];
for (int i = 0; i < slices.length; i++) slices[i] = i + 1;
Random rand = new Random();
rand.shuffle(slices);
IJ.showStatus("Fitting ...");
// Use multi-threaded code for speed
int i;
for (i = 0; i < slices.length; i++) {
int slice = slices[i];
//debug(" Processing slice = %d\n", slice);
IJ.showProgress(size(), settings.numberOfPeaks);
ImageProcessor ip = stack.getProcessor(slice);
// stack processor does not set the bounds required by ImageConverter
ip.setRoi(roi);
FitJob job = new FitJob(slice, ImageConverter.getData(ip), roi);
engine.run(job);
if (sampleSizeReached() || Utils.isInterrupted()) {
break;
}
}
if (Utils.isInterrupted()) {
IJ.showProgress(1);
engine.end(true);
return false;
}
// Wait until we have enough results
while (!sampleSizeReached() && !engine.isQueueEmpty()) {
IJ.showProgress(size(), settings.numberOfPeaks);
try {
Thread.sleep(50);
} catch (InterruptedException e) {
break;
}
}
// End now if we have enough samples
engine.end(sampleSizeReached());
IJ.showStatus("");
IJ.showProgress(1);
// This count will be an over-estimate given that the provider is ahead of the consumer
// in this multi-threaded system
debug(" Processed %d/%d slices (%d peaks)", i, slices.length, size());
setParams(ANGLE, params, params_dev, sampleNew[ANGLE]);
setParams(X, params, params_dev, sampleNew[X]);
setParams(Y, params, params_dev, sampleNew[Y]);
if (settings.showHistograms) {
int[] idList = new int[NAMES.length];
int count = 0;
boolean requireRetile = false;
for (int ii = 0; ii < 3; ii++) {
if (sampleNew[ii].getN() == 0)
continue;
StoredDataStatistics stats = new StoredDataStatistics(sampleNew[ii].getValues());
idList[count++] = Utils.showHistogram(TITLE, stats, NAMES[ii], 0, 0, settings.histogramBins, "Mean = " + Utils.rounded(stats.getMean()) + ". Median = " + Utils.rounded(sampleNew[ii].getPercentile(50)));
requireRetile = requireRetile || Utils.isNewWindow();
}
if (requireRetile && count > 0) {
new WindowOrganiser().tileWindows(Arrays.copyOf(idList, count));
}
}
if (size() < 2) {
log("ERROR: Insufficient number of fitted peaks, terminating ...");
return false;
}
return true;
}
use of gdsc.smlm.engine.FitEngine in project GDSC-SMLM by aherbert.
the class TraceMolecules method fitTraces.
private void fitTraces(MemoryPeakResults results, Trace[] traces) {
// Check if the original image is open and the fit configuration can be extracted
ImageSource source = results.getSource();
if (source == null)
return;
if (!source.open())
return;
FitEngineConfiguration config = (FitEngineConfiguration) XmlUtils.fromXML(results.getConfiguration());
if (config == null)
return;
// Show a dialog asking if the traces should be refit
ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Do you want to fit the traces as a single peak using a combined image?");
gd.addCheckbox("Fit_closest_to_centroid", !fitOnlyCentroid);
gd.addSlider("Distance_threshold", 0.01, 3, distanceThreshold);
gd.addSlider("Expansion_factor", 1, 4.5, expansionFactor);
// Allow fitting settings to be adjusted
FitConfiguration fitConfig = config.getFitConfiguration();
gd.addMessage("--- Gaussian fitting ---");
String[] filterTypes = SettingsManager.getNames((Object[]) DataFilterType.values());
gd.addChoice("Spot_filter_type", filterTypes, filterTypes[config.getDataFilterType().ordinal()]);
String[] filterNames = SettingsManager.getNames((Object[]) DataFilter.values());
gd.addChoice("Spot_filter", filterNames, filterNames[config.getDataFilter(0).ordinal()]);
gd.addSlider("Smoothing", 0, 2.5, config.getSmooth(0));
gd.addSlider("Search_width", 0.5, 2.5, config.getSearch());
gd.addSlider("Border", 0.5, 2.5, config.getBorder());
gd.addSlider("Fitting_width", 2, 4.5, config.getFitting());
String[] solverNames = SettingsManager.getNames((Object[]) FitSolver.values());
gd.addChoice("Fit_solver", solverNames, solverNames[fitConfig.getFitSolver().ordinal()]);
String[] functionNames = SettingsManager.getNames((Object[]) FitFunction.values());
gd.addChoice("Fit_function", functionNames, functionNames[fitConfig.getFitFunction().ordinal()]);
String[] criteriaNames = SettingsManager.getNames((Object[]) FitCriteria.values());
gd.addChoice("Fit_criteria", criteriaNames, criteriaNames[fitConfig.getFitCriteria().ordinal()]);
gd.addNumericField("Significant_digits", fitConfig.getSignificantDigits(), 0);
gd.addNumericField("Coord_delta", fitConfig.getDelta(), 4);
gd.addNumericField("Lambda", fitConfig.getLambda(), 4);
gd.addNumericField("Max_iterations", fitConfig.getMaxIterations(), 0);
gd.addNumericField("Fail_limit", config.getFailuresLimit(), 0);
gd.addCheckbox("Include_neighbours", config.isIncludeNeighbours());
gd.addSlider("Neighbour_height", 0.01, 1, config.getNeighbourHeightThreshold());
gd.addSlider("Residuals_threshold", 0.01, 1, config.getResidualsThreshold());
//gd.addSlider("Duplicate_distance", 0, 1.5, fitConfig.getDuplicateDistance());
gd.addMessage("--- Peak filtering ---\nDiscard fits that shift; are too low; or expand/contract");
gd.addCheckbox("Smart_filter", fitConfig.isSmartFilter());
gd.addCheckbox("Disable_simple_filter", fitConfig.isDisableSimpleFilter());
gd.addSlider("Shift_factor", 0.01, 2, fitConfig.getCoordinateShiftFactor());
gd.addNumericField("Signal_strength", fitConfig.getSignalStrength(), 2);
gd.addNumericField("Min_photons", fitConfig.getMinPhotons(), 0);
gd.addSlider("Min_width_factor", 0, 0.99, fitConfig.getMinWidthFactor());
gd.addSlider("Width_factor", 1.01, 5, fitConfig.getWidthFactor());
gd.addNumericField("Precision", fitConfig.getPrecisionThreshold(), 2);
gd.addCheckbox("Debug_failures", debugFailures);
gd.showDialog();
if (!gd.wasOKed()) {
source.close();
return;
}
// Get parameters for the fit
fitOnlyCentroid = !gd.getNextBoolean();
distanceThreshold = (float) gd.getNextNumber();
expansionFactor = (float) gd.getNextNumber();
config.setDataFilterType(gd.getNextChoiceIndex());
config.setDataFilter(gd.getNextChoiceIndex(), Math.abs(gd.getNextNumber()), 0);
config.setSearch(gd.getNextNumber());
config.setBorder(gd.getNextNumber());
config.setFitting(gd.getNextNumber());
fitConfig.setFitSolver(gd.getNextChoiceIndex());
fitConfig.setFitFunction(gd.getNextChoiceIndex());
fitConfig.setFitCriteria(gd.getNextChoiceIndex());
fitConfig.setSignificantDigits((int) gd.getNextNumber());
fitConfig.setDelta(gd.getNextNumber());
fitConfig.setLambda(gd.getNextNumber());
fitConfig.setMaxIterations((int) gd.getNextNumber());
config.setFailuresLimit((int) gd.getNextNumber());
config.setIncludeNeighbours(gd.getNextBoolean());
config.setNeighbourHeightThreshold(gd.getNextNumber());
config.setResidualsThreshold(gd.getNextNumber());
fitConfig.setSmartFilter(gd.getNextBoolean());
fitConfig.setDisableSimpleFilter(gd.getNextBoolean());
fitConfig.setCoordinateShiftFactor(gd.getNextNumber());
fitConfig.setSignalStrength(gd.getNextNumber());
fitConfig.setMinPhotons(gd.getNextNumber());
fitConfig.setMinWidthFactor(gd.getNextNumber());
fitConfig.setWidthFactor(gd.getNextNumber());
fitConfig.setPrecisionThreshold(gd.getNextNumber());
// Check arguments
try {
Parameters.isAboveZero("Distance threshold", distanceThreshold);
Parameters.isAbove("Expansion factor", expansionFactor, 1);
Parameters.isAboveZero("Search_width", config.getSearch());
Parameters.isAboveZero("Fitting_width", config.getFitting());
Parameters.isAboveZero("Significant digits", fitConfig.getSignificantDigits());
Parameters.isAboveZero("Delta", fitConfig.getDelta());
Parameters.isAboveZero("Lambda", fitConfig.getLambda());
Parameters.isAboveZero("Max iterations", fitConfig.getMaxIterations());
Parameters.isPositive("Failures limit", config.getFailuresLimit());
Parameters.isPositive("Neighbour height threshold", config.getNeighbourHeightThreshold());
Parameters.isPositive("Residuals threshold", config.getResidualsThreshold());
Parameters.isPositive("Coordinate Shift factor", fitConfig.getCoordinateShiftFactor());
Parameters.isPositive("Signal strength", fitConfig.getSignalStrength());
Parameters.isPositive("Min photons", fitConfig.getMinPhotons());
Parameters.isPositive("Min width factor", fitConfig.getMinWidthFactor());
Parameters.isPositive("Width factor", fitConfig.getWidthFactor());
Parameters.isPositive("Precision threshold", fitConfig.getPrecisionThreshold());
} catch (IllegalArgumentException e) {
IJ.error(TITLE, e.getMessage());
source.close();
return;
}
debugFailures = gd.getNextBoolean();
if (!PeakFit.configureSmartFilter(globalSettings, filename))
return;
if (!PeakFit.configureDataFilter(globalSettings, filename, false))
return;
if (!PeakFit.configureFitSolver(globalSettings, filename, false))
return;
// Adjust settings for a single maxima
config.setIncludeNeighbours(false);
fitConfig.setDuplicateDistance(0);
// Create a fit engine
MemoryPeakResults refitResults = new MemoryPeakResults();
refitResults.copySettings(results);
refitResults.setName(results.getName() + " Trace Fit");
refitResults.setSortAfterEnd(true);
refitResults.begin();
// No border since we know where the peaks are and we must not miss them due to truncated searching
FitEngine engine = new FitEngine(config, refitResults, Prefs.getThreads(), FitQueue.BLOCKING);
// Either : Only fit the centroid
// or : Extract a bigger region, allowing all fits to run as normal and then
// find the correct spot using Euclidian distance.
// Set up the limits
final double stdDev = FastMath.max(fitConfig.getInitialPeakStdDev0(), fitConfig.getInitialPeakStdDev1());
float fitWidth = (float) (stdDev * config.getFitting() * ((fitOnlyCentroid) ? 1 : expansionFactor));
IJ.showStatus("Refitting traces ...");
List<JobItem> jobItems = new ArrayList<JobItem>(traces.length);
int singles = 0;
int fitted = 0;
for (int n = 0; n < traces.length; n++) {
Trace trace = traces[n];
if (n % 32 == 0)
IJ.showProgress(n, traces.length);
// Skip traces with one peak
if (trace.size() == 1) {
singles++;
// Use the synchronized method to avoid thread clashes with the FitEngine
refitResults.addSync(trace.getHead());
continue;
}
Rectangle bounds = new Rectangle();
double[] combinedNoise = new double[1];
float[] data = buildCombinedImage(source, trace, fitWidth, bounds, combinedNoise, false);
if (data == null)
continue;
// Fit the combined image
FitParameters params = new FitParameters();
params.noise = (float) combinedNoise[0];
float[] centre = trace.getCentroid();
if (fitOnlyCentroid) {
int newX = (int) Math.round(centre[0]) - bounds.x;
int newY = (int) Math.round(centre[1]) - bounds.y;
params.maxIndices = new int[] { newY * bounds.width + newX };
} else {
params.filter = new ArrayList<float[]>();
params.filter.add(new float[] { centre[0] - bounds.x, centre[1] - bounds.y });
params.distanceThreshold = distanceThreshold;
}
// This is not needed since the bounds are passed using the FitJob
//params.setOffset(new float[] { bounds.x, bounds.y });
int startT = trace.getHead().getFrame();
params.endT = trace.getTail().getFrame();
ParameterisedFitJob job = new ParameterisedFitJob(n, params, startT, data, bounds);
jobItems.add(new JobItem(job, trace, centre));
engine.run(job);
fitted++;
}
engine.end(false);
IJ.showStatus("");
IJ.showProgress(1);
// Check the success ...
FitStatus[] values = FitStatus.values();
int[] statusCount = new int[values.length + 1];
ArrayList<String> names = new ArrayList<String>(Arrays.asList(SettingsManager.getNames((Object[]) values)));
names.add(String.format("No maxima within %.2f of centroid", distanceThreshold));
int separated = 0;
int success = 0;
final int debugLimit = 3;
for (JobItem jobItem : jobItems) {
int id = jobItem.getId();
ParameterisedFitJob job = jobItem.job;
Trace trace = jobItem.trace;
int[] indices = job.getIndices();
FitResult fitResult = null;
int status;
if (indices.length < 1) {
status = values.length;
} else if (indices.length > 1) {
// Choose the first OK result. This is all that matters for the success reporting
for (int n = 0; n < indices.length; n++) {
if (job.getFitResult(n).getStatus() == FitStatus.OK) {
fitResult = job.getFitResult(n);
break;
}
}
// Otherwise use the closest failure.
if (fitResult == null) {
final float[] centre = traces[id].getCentroid();
double minD = Double.POSITIVE_INFINITY;
for (int n = 0; n < indices.length; n++) {
// Since the fit has failed we use the initial parameters
final double[] params = job.getFitResult(n).getInitialParameters();
final double dx = params[Gaussian2DFunction.X_POSITION] - centre[0];
final double dy = params[Gaussian2DFunction.Y_POSITION] - centre[1];
final double d = dx * dx + dy * dy;
if (minD > d) {
minD = d;
fitResult = job.getFitResult(n);
}
}
}
status = fitResult.getStatus().ordinal();
} else {
fitResult = job.getFitResult(0);
status = fitResult.getStatus().ordinal();
}
// All jobs have only one peak
statusCount[status]++;
// Debug why any fits failed
if (fitResult == null || fitResult.getStatus() != FitStatus.OK) {
refitResults.addAll(trace.getPoints());
separated += trace.size();
if (debugFailures) {
FitStatus s = (fitResult == null) ? FitStatus.UNKNOWN : fitResult.getStatus();
// Only display the first n per category to limit the number of images
double[] noise = new double[1];
if (statusCount[status] <= debugLimit) {
Rectangle bounds = new Rectangle();
buildCombinedImage(source, trace, fitWidth, bounds, noise, true);
float[] centre = trace.getCentroid();
Utils.display(String.format("Trace %d (n=%d) : x=%f,y=%f", id, trace.size(), centre[0], centre[1]), slices);
switch(s) {
case INSUFFICIENT_PRECISION:
float precision = (Float) fitResult.getStatusData();
IJ.log(String.format("Trace %d (n=%d) : %s = %f", id, trace.size(), names.get(status), precision));
break;
case INSUFFICIENT_SIGNAL:
if (noise[0] == 0)
noise[0] = getCombinedNoise(trace);
float snr = (Float) fitResult.getStatusData();
IJ.log(String.format("Trace %d (n=%d) : %s = %f (noise=%.2f)", id, trace.size(), names.get(status), snr, noise[0]));
break;
case COORDINATES_MOVED:
case OUTSIDE_FIT_REGION:
case WIDTH_DIVERGED:
float[] shift = (float[]) fitResult.getStatusData();
IJ.log(String.format("Trace %d (n=%d) : %s = %.3f,%.3f", id, trace.size(), names.get(status), shift[0], shift[1]));
break;
default:
IJ.log(String.format("Trace %d (n=%d) : %s", id, trace.size(), names.get(status)));
break;
}
}
}
} else {
success++;
if (debugFailures) {
// Only display the first n per category to limit the number of images
double[] noise = new double[1];
if (statusCount[status] <= debugLimit) {
Rectangle bounds = new Rectangle();
buildCombinedImage(source, trace, fitWidth, bounds, noise, true);
float[] centre = trace.getCentroid();
Utils.display(String.format("Trace %d (n=%d) : x=%f,y=%f", id, trace.size(), centre[0], centre[1]), slices);
}
}
}
}
IJ.log(String.format("Trace fitting : %d singles : %d / %d fitted : %d separated", singles, success, fitted, separated));
if (separated > 0) {
IJ.log("Reasons for fit failure :");
// Start at i=1 to skip FitStatus.OK
for (int i = 1; i < statusCount.length; i++) {
if (statusCount[i] != 0)
IJ.log(" " + names.get(i) + " = " + statusCount[i]);
}
}
refitResults.end();
MemoryPeakResults.addResults(refitResults);
source.close();
}
use of gdsc.smlm.engine.FitEngine in project GDSC-SMLM by aherbert.
the class PSFCreator method fitSpot.
private MemoryPeakResults fitSpot(ImageStack stack, final int width, final int height, final int x, final int y) {
Rectangle regionBounds = null;
// Create a fit engine
MemoryPeakResults results = new MemoryPeakResults();
results.setSortAfterEnd(true);
results.begin();
FitEngine engine = new FitEngine(config, results, Prefs.getThreads(), FitQueue.BLOCKING);
List<ParameterisedFitJob> jobItems = new ArrayList<ParameterisedFitJob>(stack.getSize());
for (int slice = 1; slice <= stack.getSize(); slice++) {
// Extract the region from each frame
ImageExtractor ie = new ImageExtractor((float[]) stack.getPixels(slice), width, height);
if (regionBounds == null)
regionBounds = ie.getBoxRegionBounds(x, y, boxRadius);
float[] region = ie.crop(regionBounds);
// Fit only a spot in the centre
FitParameters params = new FitParameters();
params.maxIndices = new int[] { boxRadius * regionBounds.width + boxRadius };
ParameterisedFitJob job = new ParameterisedFitJob(slice, params, slice, region, regionBounds);
jobItems.add(job);
engine.run(job);
}
engine.end(false);
results.end();
return results;
}
Aggregations