use of org.apache.commons.math3.geometry.partitioning.Region in project GDSC-SMLM by aherbert.
the class FIRE method showInputDialog.
private boolean showInputDialog() {
ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Compute the resolution using Fourier Ring Correlation");
gd.addHelp(About.HELP_URL);
// Build a list of all images with a region ROI
List<String> titles = new LinkedList<String>();
if (WindowManager.getWindowCount() > 0) {
for (int imageID : WindowManager.getIDList()) {
ImagePlus imp = WindowManager.getImage(imageID);
if (imp != null && imp.getRoi() != null && imp.getRoi().isArea())
titles.add(imp.getTitle());
}
}
ResultsManager.addInput(gd, inputOption, InputSource.MEMORY);
ResultsManager.addInput(gd, "Input2", inputOption2, InputSource.NONE, InputSource.MEMORY);
if (!titles.isEmpty())
gd.addCheckbox((titles.size() == 1) ? "Use_ROI" : "Choose_ROI", chooseRoi);
gd.showDialog();
if (gd.wasCanceled())
return false;
inputOption = ResultsManager.getInputSource(gd);
inputOption2 = ResultsManager.getInputSource(gd);
if (!titles.isEmpty())
chooseRoi = gd.getNextBoolean();
if (!titles.isEmpty() && chooseRoi) {
if (titles.size() == 1) {
roiImage = titles.get(0);
Recorder.recordOption("Image", roiImage);
} else {
String[] items = titles.toArray(new String[titles.size()]);
gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Select the source image for the ROI");
gd.addChoice("Image", items, roiImage);
gd.showDialog();
if (gd.wasCanceled())
return false;
roiImage = gd.getNextChoice();
}
ImagePlus imp = WindowManager.getImage(roiImage);
roiBounds = imp.getRoi().getBounds();
roiImageWidth = imp.getWidth();
roiImageHeight = imp.getHeight();
} else {
roiBounds = null;
}
return true;
}
use of org.apache.commons.math3.geometry.partitioning.Region in project GDSC-SMLM by aherbert.
the class DensityImage method logDensityResults.
/**
* Output a log message of the results including the average density for localisations and the expected average.
*
* @param results
* @param density
* @param radius
* @param filtered
* @return
*/
private SummaryStatistics logDensityResults(MemoryPeakResults results, int[] density, float radius, int filtered) {
float region = (float) (radius * radius * ((useSquareApproximation) ? 4 : Math.PI));
Rectangle bounds = results.getBounds();
float area = bounds.width * bounds.height;
float expected = results.size() * region / area;
SummaryStatistics summary = new SummaryStatistics();
for (int i = 0; i < results.size(); i++) {
summary.addValue(density[i]);
}
DensityManager dm = createDensityManager(results);
// Compute this using the input density scores since the radius is the same.
final double l = (useSquareApproximation) ? dm.ripleysLFunction(radius) : dm.ripleysLFunction(density, radius);
String msg = String.format("Density %s : N=%d, %.0fpx : Radius=%s : L(r) - r = %s : E = %s, Obs = %s (%sx)", results.getName(), summary.getN(), area, rounded(radius), rounded(l - radius), rounded(expected), rounded(summary.getMean()), rounded(summary.getMean() / expected));
if (filterLocalisations)
msg += String.format(" : Filtered=%d (%s%%)", filtered, rounded(filtered * 100.0 / density.length));
IJ.log(msg);
return summary;
}
use of org.apache.commons.math3.geometry.partitioning.Region in project GDSC-SMLM by aherbert.
the class BenchmarkSpotFilter method summariseResults.
private BenchmarkFilterResult summariseResults(TIntObjectHashMap<FilterResult> filterResults, FitEngineConfiguration config, MaximaSpotFilter spotFilter, boolean relativeDistances, boolean batchSummary) {
BenchmarkFilterResult filterResult = new BenchmarkFilterResult(filterResults, config, spotFilter);
// Note:
// Although we can compute the TP/FP score as each additional spot is added
// using the RankedScoreCalculator this is not applicable to the PeakFit method.
// The method relies on all spot candidates being present in order to make a
// decision to fit the candidate as a multiple. So scoring the filter candidates using
// for example the top 10 may get a better score than if all candidates were scored
// and the scores accumulated for the top 10, it is not how the algorithm will use the
// candidate set. I.e. It does not use the top 10, then top 20 to refine the fit, etc.
// (the method is not iterative) .
// We require an assessment of how a subset of the scored candidates
// in ranked order contributes to the overall score, i.e. are the candidates ranked
// in the correct order, those most contributing to the match to the underlying data
// should be higher up and those least contributing will be at the end.
// TODO We could add some smart filtering of candidates before ranking. This would
// allow assessment of the candidate set handed to PeakFit. E.g. Threshold the image
// and only use candidates that are in the foreground region.
double[][] cumul = histogramFailures(filterResult);
// Create the overall match score
final double[] total = new double[3];
final ArrayList<ScoredSpot> allSpots = new ArrayList<BenchmarkSpotFilter.ScoredSpot>();
filterResults.forEachValue(new TObjectProcedure<FilterResult>() {
public boolean execute(FilterResult result) {
total[0] += result.result.getTP();
total[1] += result.result.getFP();
total[2] += result.result.getFN();
allSpots.addAll(Arrays.asList(result.spots));
return true;
}
});
double tp = total[0], fp = total[1], fn = total[2];
FractionClassificationResult allResult = new FractionClassificationResult(tp, fp, 0, fn);
// The number of actual results
final double n = (tp + fn);
StringBuilder sb = new StringBuilder();
double signal = (simulationParameters.minSignal + simulationParameters.maxSignal) * 0.5;
// Create the benchmark settings and the fitting settings
sb.append(imp.getStackSize()).append("\t");
final int w = lastAnalysisBorder.width;
final int h = lastAnalysisBorder.height;
sb.append(w).append("\t");
sb.append(h).append("\t");
sb.append(Utils.rounded(n)).append("\t");
double density = (n / imp.getStackSize()) / (w * h) / (simulationParameters.a * simulationParameters.a / 1e6);
sb.append(Utils.rounded(density)).append("\t");
sb.append(Utils.rounded(signal)).append("\t");
sb.append(Utils.rounded(simulationParameters.s)).append("\t");
sb.append(Utils.rounded(simulationParameters.a)).append("\t");
sb.append(Utils.rounded(simulationParameters.depth)).append("\t");
sb.append(simulationParameters.fixedDepth).append("\t");
sb.append(Utils.rounded(simulationParameters.gain)).append("\t");
sb.append(Utils.rounded(simulationParameters.readNoise)).append("\t");
sb.append(Utils.rounded(simulationParameters.b)).append("\t");
sb.append(Utils.rounded(simulationParameters.b2)).append("\t");
// Compute the noise
double noise = simulationParameters.b2;
if (simulationParameters.emCCD) {
// The b2 parameter was computed without application of the EM-CCD noise factor of 2.
//final double b2 = backgroundVariance + readVariance
// = simulationParameters.b + readVariance
// This should be applied only to the background variance.
final double readVariance = noise - simulationParameters.b;
noise = simulationParameters.b * 2 + readVariance;
}
sb.append(Utils.rounded(signal / Math.sqrt(noise))).append("\t");
sb.append(Utils.rounded(simulationParameters.s / simulationParameters.a)).append("\t");
sb.append(config.getDataFilterType()).append("\t");
//sb.append(spotFilter.getName()).append("\t");
sb.append(spotFilter.getSearch()).append("\t");
sb.append(spotFilter.getBorder()).append("\t");
sb.append(Utils.rounded(spotFilter.getSpread())).append("\t");
sb.append(config.getDataFilter(0)).append("\t");
final double param = config.getSmooth(0);
final double hwhmMin = config.getHWHMMin();
if (relativeDistances) {
sb.append(Utils.rounded(param * hwhmMin)).append("\t");
sb.append(Utils.rounded(param)).append("\t");
} else {
sb.append(Utils.rounded(param)).append("\t");
sb.append(Utils.rounded(param / hwhmMin)).append("\t");
}
sb.append(spotFilter.getDescription()).append("\t");
sb.append(lastAnalysisBorder.x).append("\t");
sb.append(MATCHING_METHOD[matchingMethod]).append("\t");
sb.append(Utils.rounded(lowerMatchDistance)).append("\t");
sb.append(Utils.rounded(matchDistance)).append("\t");
sb.append(Utils.rounded(lowerSignalFactor)).append("\t");
sb.append(Utils.rounded(upperSignalFactor));
resultPrefix = sb.toString();
// Add the results
sb.append("\t");
// Rank the scored spots by intensity
Collections.sort(allSpots);
// Produce Recall, Precision, Jaccard for each cut of the spot candidates
double[] r = new double[allSpots.size() + 1];
double[] p = new double[r.length];
double[] j = new double[r.length];
double[] c = new double[r.length];
double[] truePositives = new double[r.length];
double[] falsePositives = new double[r.length];
double[] intensity = new double[r.length];
// Note: fn = n - tp
tp = fp = 0;
int i = 1;
p[0] = 1;
FastCorrelator corr = new FastCorrelator();
double lastC = 0;
double[] i1 = new double[r.length];
double[] i2 = new double[r.length];
int ci = 0;
SimpleRegression regression = new SimpleRegression(false);
for (ScoredSpot s : allSpots) {
if (s.match) {
// Score partial matches as part true-positive and part false-positive.
// TP can be above 1 if we are allowing multiple matches.
tp += s.getScore();
fp += s.antiScore();
// Just use a rounded intensity for now
final double spotIntensity = s.getIntensity();
final long v1 = (long) Math.round(spotIntensity);
final long v2 = (long) Math.round(s.intensity);
regression.addData(spotIntensity, s.intensity);
i1[ci] = spotIntensity;
i2[ci] = s.intensity;
ci++;
corr.add(v1, v2);
lastC = corr.getCorrelation();
} else
fp++;
r[i] = (double) tp / n;
p[i] = (double) tp / (tp + fp);
// (tp+fp+fn) == (fp+n) since tp+fn=n;
j[i] = (double) tp / (fp + n);
c[i] = lastC;
truePositives[i] = tp;
falsePositives[i] = fp;
intensity[i] = s.getIntensity();
i++;
}
i1 = Arrays.copyOf(i1, ci);
i2 = Arrays.copyOf(i2, ci);
final double slope = regression.getSlope();
sb.append(Utils.rounded(slope)).append("\t");
addResult(sb, allResult, c[c.length - 1]);
// Output the match results when the recall achieves the fraction of the maximum.
double target = r[r.length - 1];
if (recallFraction < 100)
target *= recallFraction / 100.0;
int fractionIndex = 0;
while (fractionIndex < r.length && r[fractionIndex] < target) {
fractionIndex++;
}
if (fractionIndex == r.length)
fractionIndex--;
addResult(sb, new FractionClassificationResult(truePositives[fractionIndex], falsePositives[fractionIndex], 0, n - truePositives[fractionIndex]), c[fractionIndex]);
// Output the match results at the maximum jaccard score
int maxIndex = 0;
for (int ii = 1; ii < r.length; ii++) {
if (j[maxIndex] < j[ii])
maxIndex = ii;
}
addResult(sb, new FractionClassificationResult(truePositives[maxIndex], falsePositives[maxIndex], 0, n - truePositives[maxIndex]), c[maxIndex]);
sb.append(Utils.rounded(time / 1e6));
// Calculate AUC (Average precision == Area Under Precision-Recall curve)
final double auc = AUCCalculator.auc(p, r);
// Compute the AUC using the adjusted precision curve
// which uses the maximum precision for recall >= r
final double[] maxp = new double[p.length];
double max = 0;
for (int k = maxp.length; k-- > 0; ) {
if (max < p[k])
max = p[k];
maxp[k] = max;
}
final double auc2 = AUCCalculator.auc(maxp, r);
sb.append("\t").append(Utils.rounded(auc));
sb.append("\t").append(Utils.rounded(auc2));
// Output the number of fit failures that must be processed to capture fractions of the true positives
if (cumul[0].length != 0) {
sb.append("\t").append(Utils.rounded(getFailures(cumul, 0.80)));
sb.append("\t").append(Utils.rounded(getFailures(cumul, 0.90)));
sb.append("\t").append(Utils.rounded(getFailures(cumul, 0.95)));
sb.append("\t").append(Utils.rounded(getFailures(cumul, 0.99)));
sb.append("\t").append(Utils.rounded(cumul[0][cumul[0].length - 1]));
} else
sb.append("\t\t\t\t\t");
BufferedTextWindow resultsTable = getTable(batchSummary);
resultsTable.append(sb.toString());
// Store results
filterResult.auc = auc;
filterResult.auc2 = auc2;
filterResult.r = r;
filterResult.p = p;
filterResult.j = j;
filterResult.c = c;
filterResult.maxIndex = maxIndex;
filterResult.fractionIndex = fractionIndex;
filterResult.cumul = cumul;
filterResult.slope = slope;
filterResult.i1 = i1;
filterResult.i2 = i2;
filterResult.intensity = intensity;
filterResult.relativeDistances = relativeDistances;
filterResult.time = time;
return filterResult;
}
use of org.apache.commons.math3.geometry.partitioning.Region in project GDSC-SMLM by aherbert.
the class PSFDrift method computeDrift.
private void computeDrift() {
// Create a grid of XY offset positions between 0-1 for PSF insert
final double[] grid = new double[gridSize];
for (int i = 0; i < grid.length; i++) grid[i] = (double) i / gridSize;
// Configure fitting region
final int w = 2 * regionSize + 1;
centrePixel = w / 2;
// Check region size using the image PSF
double newPsfWidth = (double) imp.getWidth() / scale;
if (Math.ceil(newPsfWidth) > w)
Utils.log(TITLE + ": Fitted region size (%d) is smaller than the scaled PSF (%.1f)", w, newPsfWidth);
// Create robust PSF fitting settings
final double a = psfSettings.nmPerPixel * scale;
final double sa = PSFCalculator.squarePixelAdjustment(psfSettings.nmPerPixel * (psfSettings.fwhm / Gaussian2DFunction.SD_TO_FWHM_FACTOR), a);
fitConfig.setInitialPeakStdDev(sa / a);
fitConfig.setBackgroundFitting(backgroundFitting);
fitConfig.setNotSignalFitting(false);
fitConfig.setComputeDeviations(false);
fitConfig.setDisableSimpleFilter(true);
// Create the PSF over the desired z-depth
int depth = (int) Math.round(zDepth / psfSettings.nmPerSlice);
int startSlice = psfSettings.zCentre - depth;
int endSlice = psfSettings.zCentre + depth;
int nSlices = imp.getStackSize();
startSlice = (startSlice < 1) ? 1 : (startSlice > nSlices) ? nSlices : startSlice;
endSlice = (endSlice < 1) ? 1 : (endSlice > nSlices) ? nSlices : endSlice;
ImagePSFModel psf = createImagePSF(startSlice, endSlice);
int minz = startSlice - psfSettings.zCentre;
int maxz = endSlice - psfSettings.zCentre;
final int nZ = maxz - minz + 1;
final int gridSize2 = grid.length * grid.length;
total = nZ * gridSize2;
// Store all the fitting results
int nStartPoints = getNumberOfStartPoints();
results = new double[total * nStartPoints][];
// TODO - Add ability to iterate this, adjusting the current offset in the PSF
// each iteration
// Create a pool of workers
int nThreads = Prefs.getThreads();
BlockingQueue<Job> jobs = new ArrayBlockingQueue<Job>(nThreads * 2);
List<Worker> workers = new LinkedList<Worker>();
List<Thread> threads = new LinkedList<Thread>();
for (int i = 0; i < nThreads; i++) {
Worker worker = new Worker(jobs, psf, w, fitConfig);
Thread t = new Thread(worker);
workers.add(worker);
threads.add(t);
t.start();
}
// Fit
Utils.showStatus("Fitting ...");
final int step = Utils.getProgressInterval(total);
outer: for (int z = minz, i = 0; z <= maxz; z++) {
for (int x = 0; x < grid.length; x++) for (int y = 0; y < grid.length; y++, i++) {
if (IJ.escapePressed()) {
break outer;
}
put(jobs, new Job(z, grid[x], grid[y], i));
if (i % step == 0) {
IJ.showProgress(i, total);
}
}
}
// If escaped pressed then do not need to stop the workers, just return
if (Utils.isInterrupted()) {
IJ.showProgress(1);
return;
}
// Finish all the worker threads by passing in a null job
for (int i = 0; i < threads.size(); i++) {
put(jobs, new Job());
}
// Wait for all to finish
for (int i = 0; i < threads.size(); i++) {
try {
threads.get(i).join();
} catch (InterruptedException e) {
e.printStackTrace();
}
}
threads.clear();
IJ.showProgress(1);
IJ.showStatus("");
// Plot the average and SE for the drift curve
// Plot the recall
double[] zPosition = new double[nZ];
double[] avX = new double[nZ];
double[] seX = new double[nZ];
double[] avY = new double[nZ];
double[] seY = new double[nZ];
double[] recall = new double[nZ];
for (int z = minz, i = 0; z <= maxz; z++, i++) {
Statistics statsX = new Statistics();
Statistics statsY = new Statistics();
for (int s = 0; s < nStartPoints; s++) {
int resultPosition = i * gridSize2 + s * total;
final int endResultPosition = resultPosition + gridSize2;
while (resultPosition < endResultPosition) {
if (results[resultPosition] != null) {
statsX.add(results[resultPosition][0]);
statsY.add(results[resultPosition][1]);
}
resultPosition++;
}
}
zPosition[i] = z * psfSettings.nmPerSlice;
avX[i] = statsX.getMean();
seX[i] = statsX.getStandardError();
avY[i] = statsY.getMean();
seY[i] = statsY.getStandardError();
recall[i] = (double) statsX.getN() / (nStartPoints * gridSize2);
}
// Find the range from the z-centre above the recall limit
int centre = 0;
for (int slice = startSlice, i = 0; slice <= endSlice; slice++, i++) {
if (slice == psfSettings.zCentre) {
centre = i;
break;
}
}
if (recall[centre] < recallLimit)
return;
int start = centre, end = centre;
for (int i = centre; i-- > 0; ) {
if (recall[i] < recallLimit)
break;
start = i;
}
for (int i = centre; ++i < recall.length; ) {
if (recall[i] < recallLimit)
break;
end = i;
}
int iterations = 1;
LoessInterpolator loess = null;
if (smoothing > 0)
loess = new LoessInterpolator(smoothing, iterations);
double[][] smoothx = displayPlot("Drift X", "X (nm)", zPosition, avX, seX, loess, start, end);
double[][] smoothy = displayPlot("Drift Y", "Y (nm)", zPosition, avY, seY, loess, start, end);
displayPlot("Recall", "Recall", zPosition, recall, null, null, start, end);
WindowOrganiser wo = new WindowOrganiser();
wo.tileWindows(idList);
// Ask the user if they would like to store them in the image
GenericDialog gd = new GenericDialog(TITLE);
gd.enableYesNoCancel();
gd.hideCancelButton();
startSlice = psfSettings.zCentre - (centre - start);
endSlice = psfSettings.zCentre + (end - centre);
gd.addMessage(String.format("Save the drift to the PSF?\n \nSlices %d (%s nm) - %d (%s nm) above recall limit", startSlice, Utils.rounded(zPosition[start]), endSlice, Utils.rounded(zPosition[end])));
gd.addMessage("Optionally average the end points to set drift outside the limits.\n(Select zero to ignore)");
gd.addSlider("Number_of_points", 0, 10, positionsToAverage);
gd.showDialog();
if (gd.wasOKed()) {
positionsToAverage = Math.abs((int) gd.getNextNumber());
ArrayList<PSFOffset> offset = new ArrayList<PSFOffset>();
final double pitch = psfSettings.nmPerPixel;
int j = 0, jj = 0;
for (int i = start, slice = startSlice; i <= end; slice++, i++) {
j = findCentre(zPosition[i], smoothx, j);
if (j == -1) {
Utils.log("Failed to find the offset for depth %.2f", zPosition[i]);
continue;
}
// The offset should store the difference to the centre in pixels so divide by the pixel pitch
double cx = smoothx[1][j] / pitch;
double cy = smoothy[1][j] / pitch;
jj = findOffset(slice, jj);
if (jj != -1) {
cx += psfSettings.offset[jj].cx;
cy += psfSettings.offset[jj].cy;
}
offset.add(new PSFOffset(slice, cx, cy));
}
addMissingOffsets(startSlice, endSlice, nSlices, offset);
psfSettings.offset = offset.toArray(new PSFOffset[offset.size()]);
psfSettings.addNote(TITLE, String.format("Solver=%s, Region=%d", PeakFit.getSolverName(fitConfig), regionSize));
imp.setProperty("Info", XmlUtils.toXML(psfSettings));
}
}
use of org.apache.commons.math3.geometry.partitioning.Region in project GDSC-SMLM by aherbert.
the class SpotInspector method run.
/*
* (non-Javadoc)
*
* @see ij.plugin.PlugIn#run(java.lang.String)
*/
public void run(String arg) {
SMLMUsageTracker.recordPlugin(this.getClass(), arg);
if (MemoryPeakResults.isMemoryEmpty()) {
IJ.error(TITLE, "No localisations in memory");
return;
}
if (!showDialog())
return;
// Load the results
results = ResultsManager.loadInputResults(inputOption, false);
if (results == null || results.size() == 0) {
IJ.error(TITLE, "No results could be loaded");
IJ.showStatus("");
return;
}
// Check if the original image is open
ImageSource source = results.getSource();
if (source == null) {
IJ.error(TITLE, "Unknown original source image");
return;
}
source = source.getOriginal();
if (!source.open()) {
IJ.error(TITLE, "Cannot open original source image: " + source.toString());
return;
}
final float stdDevMax = getStandardDeviation(results);
if (stdDevMax < 0) {
// TODO - Add dialog to get the initial peak width
IJ.error(TITLE, "Fitting configuration (for initial peak width) is not available");
return;
}
// Rank spots
rankedResults = new ArrayList<PeakResultRank>(results.size());
final double a = results.getNmPerPixel();
final double gain = results.getGain();
final boolean emCCD = results.isEMCCD();
for (PeakResult r : results.getResults()) {
float[] score = getScore(r, a, gain, emCCD, stdDevMax);
rankedResults.add(new PeakResultRank(r, score[0], score[1]));
}
Collections.sort(rankedResults);
// Prepare results table. Get bias if necessary
if (showCalibratedValues) {
// Get a bias if required
Calibration calibration = results.getCalibration();
if (calibration.getBias() == 0) {
ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Calibrated results requires a camera bias");
gd.addNumericField("Camera_bias (ADUs)", calibration.getBias(), 2);
gd.showDialog();
if (!gd.wasCanceled()) {
calibration.setBias(Math.abs(gd.getNextNumber()));
}
}
}
IJTablePeakResults table = new IJTablePeakResults(false, results.getName(), true);
table.copySettings(results);
table.setTableTitle(TITLE);
table.setAddCounter(true);
table.setShowCalibratedValues(showCalibratedValues);
table.begin();
// Add a mouse listener to jump to the frame for the clicked line
textPanel = table.getResultsWindow().getTextPanel();
// We must ignore old instances of this class from the mouse listeners
id = ++currentId;
textPanel.addMouseListener(this);
// Add results to the table
int n = 0;
for (PeakResultRank rank : rankedResults) {
rank.rank = n++;
PeakResult r = rank.peakResult;
table.add(r.getFrame(), r.origX, r.origY, r.origValue, r.error, r.noise, r.params, r.paramsStdDev);
}
table.end();
if (plotScore || plotHistogram) {
// Get values for the plots
float[] xValues = null, yValues = null;
double yMin, yMax;
int spotNumber = 0;
xValues = new float[rankedResults.size()];
yValues = new float[xValues.length];
for (PeakResultRank rank : rankedResults) {
xValues[spotNumber] = spotNumber + 1;
yValues[spotNumber++] = recoverScore(rank.score);
}
// Set the min and max y-values using 1.5 x IQR
DescriptiveStatistics stats = new DescriptiveStatistics();
for (float v : yValues) stats.addValue(v);
if (removeOutliers) {
double lower = stats.getPercentile(25);
double upper = stats.getPercentile(75);
double iqr = upper - lower;
yMin = FastMath.max(lower - iqr, stats.getMin());
yMax = FastMath.min(upper + iqr, stats.getMax());
IJ.log(String.format("Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax));
} else {
yMin = stats.getMin();
yMax = stats.getMax();
IJ.log(String.format("Data range: %f - %f", yMin, yMax));
}
plotScore(xValues, yValues, yMin, yMax);
plotHistogram(yValues, yMin, yMax);
}
// Extract spots into a stack
final int w = source.getWidth();
final int h = source.getHeight();
final int size = 2 * radius + 1;
ImageStack spots = new ImageStack(size, size, rankedResults.size());
// To assist the extraction of data from the image source, process them in time order to allow
// frame caching. Then set the appropriate slice in the result stack
Collections.sort(rankedResults, new Comparator<PeakResultRank>() {
public int compare(PeakResultRank o1, PeakResultRank o2) {
if (o1.peakResult.getFrame() < o2.peakResult.getFrame())
return -1;
if (o1.peakResult.getFrame() > o2.peakResult.getFrame())
return 1;
return 0;
}
});
for (PeakResultRank rank : rankedResults) {
PeakResult r = rank.peakResult;
// Extract image
// Note that the coordinates are relative to the middle of the pixel (0.5 offset)
// so do not round but simply convert to int
final int x = (int) (r.params[Gaussian2DFunction.X_POSITION]);
final int y = (int) (r.params[Gaussian2DFunction.Y_POSITION]);
// Extract a region but crop to the image bounds
int minX = x - radius;
int minY = y - radius;
int maxX = FastMath.min(x + radius + 1, w);
int maxY = FastMath.min(y + radius + 1, h);
int padX = 0, padY = 0;
if (minX < 0) {
padX = -minX;
minX = 0;
}
if (minY < 0) {
padY = -minY;
minY = 0;
}
int sizeX = maxX - minX;
int sizeY = maxY - minY;
float[] data = source.get(r.getFrame(), new Rectangle(minX, minY, sizeX, sizeY));
// Prevent errors with missing data
if (data == null)
data = new float[sizeX * sizeY];
ImageProcessor spotIp = new FloatProcessor(sizeX, sizeY, data, null);
// Pad if necessary, i.e. the crop is too small for the stack
if (padX > 0 || padY > 0 || sizeX < size || sizeY < size) {
ImageProcessor spotIp2 = spotIp.createProcessor(size, size);
spotIp2.insert(spotIp, padX, padY);
spotIp = spotIp2;
}
int slice = rank.rank + 1;
spots.setPixels(spotIp.getPixels(), slice);
spots.setSliceLabel(Utils.rounded(rank.originalScore), slice);
}
source.close();
ImagePlus imp = Utils.display(TITLE, spots);
imp.setRoi((PointRoi) null);
// Make bigger
for (int i = 10; i-- > 0; ) imp.getWindow().getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
}
Aggregations