use of gdsc.smlm.engine.filter.DistanceResultFilter in project GDSC-SMLM by aherbert.
the class FitWorker method run.
/**
* Locate all the peaks in the image specified by the fit job
* <p>
* WARNING: The FitWorker fits a sub-region of the data for each maxima. It then updates the FitResult parameters
* with an offset reflecting the position. The initialParameters are not updated with this offset unless configured.
*
* @param job
* The fit job
*/
public void run(FitJob job) {
final long start = System.nanoTime();
job.start();
this.job = job;
benchmarking = false;
this.slice = job.slice;
// Used for debugging
//if (logger == null) logger = new gdsc.fitting.logging.ConsoleLogger();
// Crop to the ROI
cc = new CoordinateConverter(job.bounds);
final int width = cc.dataBounds.width;
final int height = cc.dataBounds.height;
borderLimitX = width - border;
borderLimitY = height - border;
data = job.data;
FitParameters params = job.getFitParameters();
this.endT = (params != null) ? params.endT : -1;
candidates = indentifySpots(job, width, height, params);
if (candidates.size == 0) {
finish(job, start);
return;
}
fittedBackground = 0;
// Always get the noise and store it with the results.
if (params != null && !Float.isNaN(params.noise)) {
noise = params.noise;
fitConfig.setNoise(noise);
} else if (calculateNoise) {
noise = estimateNoise(width, height);
fitConfig.setNoise(noise);
}
//System.out.printf("Slice %d : Noise = %g\n", slice, noise);
if (logger != null)
logger.info("Slice %d: Noise = %f", slice, noise);
final ImageExtractor ie = new ImageExtractor(data, width, height);
double[] region = null;
if (params != null && params.fitTask == FitTask.MAXIMA_IDENITIFICATION) {
final float sd0 = (float) fitConfig.getInitialPeakStdDev0();
final float sd1 = (float) fitConfig.getInitialPeakStdDev1();
for (int n = 0; n < candidates.getSize(); n++) {
// Find the background using the perimeter of the data.
// TODO - Perhaps the Gaussian Fitter should be used to produce the initial estimates but no actual fit done.
// This would produce coords using the centre-of-mass.
final Candidate candidate = candidates.get(n);
final int x = candidate.x;
final int y = candidate.y;
final Rectangle regionBounds = ie.getBoxRegionBounds(x, y, fitting);
region = ie.crop(regionBounds, region);
final float b = (float) Gaussian2DFitter.getBackground(region, regionBounds.width, regionBounds.height, 1);
// Offset the coords to the centre of the pixel. Note the bounds will be added later.
// Subtract the background to get the amplitude estimate then convert to signal.
final float amplitude = candidate.intensity - ((relativeIntensity) ? 0 : b);
final float signal = (float) (amplitude * 2.0 * Math.PI * sd0 * sd1);
final float[] peakParams = new float[] { b, signal, 0, x + 0.5f, y + 0.5f, sd0, sd1 };
final int index = y * width + x;
sliceResults.add(createResult(cc.fromDataToGlobalX(x), cc.fromDataToGlobalY(y), data[index], 0, noise, peakParams, null, n));
}
} else {
initialiseFitting();
// Smooth the data to provide initial background estimates
final BlockAverageDataProcessor processor = new BlockAverageDataProcessor(1, 1);
final float[] smoothedData = processor.process(data, width, height);
final ImageExtractor ie2 = new ImageExtractor(smoothedData, width, height);
// Allow the results to be filtered for certain peaks
if (params != null && params.filter != null) {
resultFilter = new DistanceResultFilter(params.filter, params.distanceThreshold, candidates.getlength());
//filter = new OptimumDistanceResultFilter(params.filter, params.distanceThreshold, maxIndices.length);
}
// The SpotFitter is used to create a dynamic MultiPathFitResult object.
// This is then passed to a multi-path filter. Thus the same fitting decision process
// is used when benchmarking and when running on actual data.
// Note: The SpotFitter labels each PreprocessedFitResult using the offset in the FitResult object.
// The initial params and deviations can then be extracted for the results that pass the filter.
MultiPathFilter filter;
IMultiPathFitResults multiPathResults = this;
SelectedResultStore store = this;
coordinateStore = coordinateStore.resize(width, height);
if (params != null && params.fitTask == FitTask.BENCHMARKING) {
// Run filtering as normal. However in the event that a candidate is missed or some
// results are not generated we must generate them. This is done in the complete(int)
// method if we set the benchmarking flag.
benchmarking = true;
// Filter using the benchmark filter
filter = params.benchmarkFilter;
if (filter == null) {
// Create a default filter using the standard FitConfiguration to ensure sensible fits
// are stored as the current slice results.
// Note the current fit configuration for benchmarking may have minimal filtering settings
// so we do not use that object.
final FitConfiguration tmp = new FitConfiguration();
final double residualsThreshold = 0.4;
filter = new MultiPathFilter(tmp, createMinimalFilter(), residualsThreshold);
}
} else {
// Filter using the configuration
filter = new MultiPathFilter(fitConfig, createMinimalFilter(), config.getResidualsThreshold());
}
// If we are benchmarking then do not generate results dynamically since we will store all
// results in the fit job
dynamicMultiPathFitResult = new DynamicMultiPathFitResult(ie, ie2, !benchmarking);
// Debug where the fit config may be different between benchmarking and fitting
if (slice == -1) {
SettingsManager.saveFitEngineConfiguration(config, String.format("/tmp/config.%b.xml", benchmarking));
Utils.write(String.format("/tmp/filter.%b.xml", benchmarking), filter.toXML());
//filter.setDebugFile(String.format("/tmp/fitWorker.%b.txt", benchmarking));
StringBuilder sb = new StringBuilder();
sb.append((benchmarking) ? ((gdsc.smlm.results.filter.Filter) filter.getFilter()).toXML() : fitConfig.getSmartFilter().toXML()).append("\n");
sb.append(((gdsc.smlm.results.filter.Filter) filter.getMinimalFilter()).toXML()).append("\n");
sb.append(filter.residualsThreshold).append("\n");
sb.append(config.getFailuresLimit()).append("\n");
sb.append(fitConfig.getDuplicateDistance()).append("\n");
if (spotFilter != null)
sb.append(spotFilter.getDescription()).append("\n");
sb.append("MaxCandidate = ").append(candidates.getSize()).append("\n");
for (int i = 0; i < candidates.list.length; i++) {
sb.append(String.format("Fit %d [%d,%d = %.1f]\n", i, candidates.get(i).x, candidates.get(i).y, candidates.get(i).intensity));
}
Utils.write(String.format("/tmp/candidates.%b.xml", benchmarking), sb.toString());
}
filter.select(multiPathResults, config.getFailuresLimit(), true, store, coordinateStore);
if (logger != null)
logger.info("Slice %d: %d / %d", slice, success, candidates.getSize());
// Result filter post-processing
if (resultFilter != null) {
resultFilter.finalise();
job.setResults(resultFilter.getResults());
job.setIndices(resultFilter.getMaxIndices());
for (int i = 0; i < resultFilter.getFilteredCount(); i++) {
job.setFitResult(i, resultFilter.getFitResults()[i]);
}
sliceResults.clear();
sliceResults.addAll(resultFilter.getResults());
}
}
// Add the ROI bounds to the fitted peaks
final float offsetx = cc.dataBounds.x;
final float offsety = cc.dataBounds.y;
for (int i = 0; i < sliceResults.size(); i++) {
final PeakResult result = sliceResults.get(i);
result.params[Gaussian2DFunction.X_POSITION] += offsetx;
result.params[Gaussian2DFunction.Y_POSITION] += offsety;
}
this.results.addAll(sliceResults);
finish(job, start);
}
Aggregations