use of uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder in project GDSC-SMLM by aherbert.
the class CreateData method showSummary.
private double showSummary(List<? extends FluorophoreSequenceModel> fluorophores, List<LocalisationModel> localisations) {
IJ.showStatus("Calculating statistics ...");
final Statistics[] stats = new Statistics[NAMES.length];
for (int i = 0; i < stats.length; i++) {
stats[i] = (settings.getShowHistograms() || alwaysRemoveOutliers[i]) ? new StoredDataStatistics() : new Statistics();
}
// Find the largest timepoint
final ImagePlus outputImp = WindowManager.getImage(benchmarkImageId);
int frameCount;
if (outputImp == null) {
sortLocalisationsByTime(localisations);
frameCount = localisations.get(localisations.size() - 1).getTime();
} else {
frameCount = outputImp.getStackSize();
}
final int[] countHistogram = new int[frameCount + 1];
// Use the localisations that were drawn to create the sampled on/off times
rebuildNeighbours(localisations);
// Assume that there is at least one localisation
final LocalisationModel first = localisations.get(0);
// The current localisation
int currentId = first.getId();
// The last time this localisation was on
int lastT = first.getTime();
// Number of blinks
int blinks = 0;
// On-time of current pulse
int currentT = 0;
double signal = 0;
final double centreOffset = settings.getSize() * 0.5;
// Used to convert the sampled times in frames into seconds
final double framesPerSecond = 1000.0 / settings.getExposureTime();
// final double gain = new CreateDataSettingsHelper(settings).getTotalGainSafe();
for (final LocalisationModel l : localisations) {
final double[] data = l.getData();
if (data == null) {
throw new IllegalStateException("No localisation data. This should not happen!");
}
final double noise = data[1];
final double sx = data[2];
final double sy = data[3];
final double intensityInPhotons = data[4];
// Q. What if the noise is zero, i.e. no background photon / read noise?
// Just ignore it at current. This is only an approximation to the SNR estimate
// if this is not a Gaussian spot.
final double snr = Gaussian2DPeakResultHelper.getMeanSignalUsingP05(intensityInPhotons, sx, sy) / noise;
stats[SIGNAL].add(intensityInPhotons);
stats[NOISE].add(noise);
if (noise != 0) {
stats[SNR].add(snr);
}
// if (l.isContinuous())
if (l.getNext() != null && l.getPrevious() != null) {
stats[SIGNAL_CONTINUOUS].add(intensityInPhotons);
if (noise != 0) {
stats[SNR_CONTINUOUS].add(snr);
}
}
final int id = l.getId();
// Check if this a new fluorophore
if (currentId != id) {
// Add previous fluorophore
stats[SAMPLED_BLINKS].add(blinks);
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
stats[TOTAL_SIGNAL].add(signal);
// Reset
blinks = 0;
currentT = 1;
currentId = id;
signal = intensityInPhotons;
} else {
signal += intensityInPhotons;
// Check if the current fluorophore pulse is broken (i.e. a blink)
if (l.getTime() - 1 > lastT) {
blinks++;
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
currentT = 1;
stats[SAMPLED_T_OFF].add(((l.getTime() - 1) - lastT) / framesPerSecond);
} else {
// Continuous on-time
currentT++;
}
}
lastT = l.getTime();
countHistogram[lastT]++;
stats[X].add((l.getX() - centreOffset) * settings.getPixelPitch());
stats[Y].add((l.getY() - centreOffset) * settings.getPixelPitch());
stats[Z].add(l.getZ() * settings.getPixelPitch());
}
// Final fluorophore
stats[SAMPLED_BLINKS].add(blinks);
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
stats[TOTAL_SIGNAL].add(signal);
// Samples per frame
for (int t = 1; t < countHistogram.length; t++) {
stats[SAMPLES].add(countHistogram[t]);
}
if (fluorophores != null) {
for (final FluorophoreSequenceModel f : fluorophores) {
stats[BLINKS].add(f.getNumberOfBlinks());
// On-time
for (final double t : f.getOnTimes()) {
stats[T_ON].add(t);
}
// Off-time
for (final double t : f.getOffTimes()) {
stats[T_OFF].add(t);
}
}
} else {
// show no blinks
stats[BLINKS].add(0);
stats[T_ON].add(1);
}
if (results != null) {
// Convert depth-of-field to pixels
final double depth = settings.getDepthOfField() / settings.getPixelPitch();
try {
// Get widths
final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
wp.getW();
stats[WIDTH].add(wp.wx);
} catch (final DataException ex) {
ImageJUtils.log("Unable to compute width: " + ex.getMessage());
}
try {
// Get z depth
final StandardResultProcedure sp = new StandardResultProcedure(results, DistanceUnit.PIXEL);
sp.getXyz();
// Get precision
final PrecisionResultProcedure pp = new PrecisionResultProcedure(results);
pp.getPrecision();
stats[PRECISION].add(pp.precisions);
for (int i = 0; i < pp.size(); i++) {
if (Math.abs(sp.z[i]) < depth) {
stats[PRECISION_IN_FOCUS].add(pp.precisions[i]);
}
}
} catch (final DataException ex) {
ImageJUtils.log("Unable to compute LSE precision: " + ex.getMessage());
}
// Compute density per frame. Multi-thread for speed
if (settings.getDensityRadius() > 0) {
final int threadCount = Prefs.getThreads();
final Ticker ticker = ImageJUtils.createTicker(results.getLastFrame(), threadCount, "Calculating density ...");
final ExecutorService threadPool = Executors.newFixedThreadPool(threadCount);
final List<Future<?>> futures = new LinkedList<>();
final TFloatArrayList coordsX = new TFloatArrayList();
final TFloatArrayList coordsY = new TFloatArrayList();
final Statistics densityStats = stats[DENSITY];
final float radius = (float) (settings.getDensityRadius() * getHwhm());
final Rectangle bounds = results.getBounds();
final double area = (double) bounds.width * bounds.height;
// Store the density for each result.
final int[] allDensity = new int[results.size()];
final FrameCounter counter = results.newFrameCounter();
results.forEach((PeakResultProcedure) result -> {
if (counter.advance(result.getFrame())) {
counter.increment(runDensityCalculation(threadPool, futures, coordsX, coordsY, densityStats, radius, area, allDensity, counter.getCount(), ticker));
}
coordsX.add(result.getXPosition());
coordsY.add(result.getYPosition());
});
runDensityCalculation(threadPool, futures, coordsX, coordsY, densityStats, radius, area, allDensity, counter.getCount(), ticker);
ConcurrencyUtils.waitForCompletionUnchecked(futures);
threadPool.shutdown();
ImageJUtils.finished();
// Split results into singles (density = 0) and clustered (density > 0)
final MemoryPeakResults singles = copyMemoryPeakResults("No Density");
final MemoryPeakResults clustered = copyMemoryPeakResults("Density");
counter.reset();
results.forEach((PeakResultProcedure) result -> {
final int density = allDensity[counter.getAndIncrement()];
result.setOrigValue(density);
if (density == 0) {
singles.add(result);
} else {
clustered.add(result);
}
});
}
}
final StringBuilder sb = new StringBuilder();
sb.append(datasetNumber).append('\t');
if (settings.getCameraType() == CameraType.SCMOS) {
sb.append("sCMOS (").append(settings.getCameraModelName()).append(") ");
final Rectangle bounds = cameraModel.getBounds();
sb.append(" ").append(bounds.x).append(",").append(bounds.y);
final int size = settings.getSize();
sb.append(" ").append(size).append("x").append(size);
} else if (CalibrationProtosHelper.isCcdCameraType(settings.getCameraType())) {
sb.append(CalibrationProtosHelper.getName(settings.getCameraType()));
final int size = settings.getSize();
sb.append(" ").append(size).append("x").append(size);
if (settings.getCameraType() == CameraType.EMCCD) {
sb.append(" EM=").append(settings.getEmGain());
}
sb.append(" CG=").append(settings.getCameraGain());
sb.append(" RN=").append(settings.getReadNoise());
sb.append(" B=").append(settings.getBias());
} else {
throw new IllegalStateException();
}
sb.append(" QE=").append(settings.getQuantumEfficiency()).append('\t');
sb.append(settings.getPsfModel());
if (psfModelType == PSF_MODEL_IMAGE) {
sb.append(" Image").append(settings.getPsfImageName());
} else if (psfModelType == PSF_MODEL_ASTIGMATISM) {
sb.append(" model=").append(settings.getAstigmatismModel());
} else {
sb.append(" DoF=").append(MathUtils.rounded(settings.getDepthOfFocus()));
if (settings.getEnterWidth()) {
sb.append(" SD=").append(MathUtils.rounded(settings.getPsfSd()));
} else {
sb.append(" λ=").append(MathUtils.rounded(settings.getWavelength()));
sb.append(" NA=").append(MathUtils.rounded(settings.getNumericalAperture()));
}
}
sb.append('\t');
sb.append((fluorophores == null) ? localisations.size() : fluorophores.size()).append('\t');
sb.append(stats[SAMPLED_BLINKS].getN() + (int) stats[SAMPLED_BLINKS].getSum()).append('\t');
sb.append(localisations.size()).append('\t');
sb.append(frameCount).append('\t');
sb.append(MathUtils.rounded(areaInUm)).append('\t');
sb.append(MathUtils.rounded(localisations.size() / (areaInUm * frameCount), 4)).append('\t');
sb.append(MathUtils.rounded(getHwhm(), 4)).append('\t');
double sd = getPsfSd();
sb.append(MathUtils.rounded(sd, 4)).append('\t');
sd *= settings.getPixelPitch();
final double sa = PsfCalculator.squarePixelAdjustment(sd, settings.getPixelPitch()) / settings.getPixelPitch();
sb.append(MathUtils.rounded(sa, 4)).append('\t');
// Width not valid for the Image PSF.
// Q. Is this true? We can approximate the FHWM for a spot-like image PSF.
final int nStats = (psfModelType == PSF_MODEL_IMAGE) ? stats.length - 1 : stats.length;
for (int i = 0; i < nStats; i++) {
final double centre = (alwaysRemoveOutliers[i]) ? ((StoredDataStatistics) stats[i]).getStatistics().getPercentile(50) : stats[i].getMean();
sb.append(MathUtils.rounded(centre, 4)).append('\t');
}
createSummaryTable().accept(sb.toString());
// Show histograms
if (settings.getShowHistograms() && !java.awt.GraphicsEnvironment.isHeadless()) {
IJ.showStatus("Calculating histograms ...");
final boolean[] chosenHistograms = getChoosenHistograms();
final WindowOrganiser wo = new WindowOrganiser();
final HistogramPlotBuilder builder = new HistogramPlotBuilder(TITLE);
for (int i = 0; i < NAMES.length; i++) {
if (chosenHistograms[i]) {
builder.setData((StoredDataStatistics) stats[i]).setName(NAMES[i]).setIntegerBins(integerDisplay[i]).setRemoveOutliersOption((settings.getRemoveOutliers() || alwaysRemoveOutliers[i]) ? 2 : 0).setNumberOfBins(settings.getHistogramBins()).show(wo);
}
}
wo.tile();
}
IJ.showStatus("");
return stats[SIGNAL].getMean();
}
use of uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder in project GDSC-SMLM by aherbert.
the class BenchmarkSpotFit method showDoubleHistogram.
private double[] showDoubleHistogram(StoredDataStatistics[][] stats, final int index, WindowOrganiser wo, double[][] matchScores) {
final String xLabel = filterCriteria[index].name;
LowerLimit lower = filterCriteria[index].lower;
UpperLimit upper = filterCriteria[index].upper;
double[] jaccard = null;
double[] metric = null;
double maxJaccard = 0;
if (index <= FILTER_PRECISION && (settings.showFilterScoreHistograms || upper.requiresJaccard || lower.requiresJaccard)) {
// Jaccard score verses the range of the metric
for (final double[] d : matchScores) {
if (!Double.isFinite(d[index])) {
System.out.printf("Error in fit data [%d]: %s%n", index, d[index]);
}
}
// Do not use Double.compare(double, double) so we get exceptions in the sort for inf/nan
Arrays.sort(matchScores, (o1, o2) -> {
if (o1[index] < o2[index]) {
return -1;
}
if (o1[index] > o2[index]) {
return 1;
}
return 0;
});
final int scoreIndex = FILTER_PRECISION + 1;
final int n = results.size();
double tp = 0;
double fp = 0;
jaccard = new double[matchScores.length + 1];
metric = new double[jaccard.length];
for (int k = 0; k < matchScores.length; k++) {
final double score = matchScores[k][scoreIndex];
tp += score;
fp += (1 - score);
jaccard[k + 1] = tp / (fp + n);
metric[k + 1] = matchScores[k][index];
}
metric[0] = metric[1];
maxJaccard = MathUtils.max(jaccard);
if (settings.showFilterScoreHistograms) {
final String title = TITLE + " Jaccard " + xLabel;
final Plot plot = new Plot(title, xLabel, "Jaccard");
plot.addPoints(metric, jaccard, Plot.LINE);
// Remove outliers
final double[] limitsx = MathUtils.limits(metric);
final Percentile p = new Percentile();
final double l = p.evaluate(metric, 25);
final double u = p.evaluate(metric, 75);
final double iqr = 1.5 * (u - l);
limitsx[1] = Math.min(limitsx[1], u + iqr);
plot.setLimits(limitsx[0], limitsx[1], 0, MathUtils.max(jaccard));
ImageJUtils.display(title, plot, wo);
}
}
// [0] is all
// [1] is matches
// [2] is no match
final StoredDataStatistics s1 = stats[0][index];
final StoredDataStatistics s2 = stats[1][index];
final StoredDataStatistics s3 = stats[2][index];
if (s1.getN() == 0) {
return new double[4];
}
final DescriptiveStatistics d = s1.getStatistics();
double median = 0;
Plot plot = null;
String title = null;
if (settings.showFilterScoreHistograms) {
median = d.getPercentile(50);
final String label = String.format("n = %d. Median = %s nm", s1.getN(), MathUtils.rounded(median));
final HistogramPlot histogramPlot = new HistogramPlotBuilder(TITLE, s1, xLabel).setMinBinWidth(filterCriteria[index].minBinWidth).setRemoveOutliersOption((filterCriteria[index].restrictRange) ? 1 : 0).setPlotLabel(label).build();
final PlotWindow plotWindow = histogramPlot.show(wo);
if (plotWindow == null) {
IJ.log("Failed to show the histogram: " + xLabel);
return new double[4];
}
title = plotWindow.getTitle();
// Reverse engineer the histogram settings
plot = histogramPlot.getPlot();
final double[] xvalues = histogramPlot.getPlotXValues();
final int bins = xvalues.length;
final double yMin = xvalues[0];
final double binSize = xvalues[1] - xvalues[0];
final double yMax = xvalues[0] + (bins - 1) * binSize;
if (s2.getN() > 0) {
final double[] values = s2.getValues();
final double[][] hist = HistogramPlot.calcHistogram(values, yMin, yMax, bins);
if (hist[0].length > 0) {
plot.setColor(Color.red);
plot.addPoints(hist[0], hist[1], Plot.BAR);
ImageJUtils.display(title, plot);
}
}
if (s3.getN() > 0) {
final double[] values = s3.getValues();
final double[][] hist = HistogramPlot.calcHistogram(values, yMin, yMax, bins);
if (hist[0].length > 0) {
plot.setColor(Color.blue);
plot.addPoints(hist[0], hist[1], Plot.BAR);
ImageJUtils.display(title, plot);
}
}
}
// Do cumulative histogram
final double[][] h1 = MathUtils.cumulativeHistogram(s1.getValues(), true);
final double[][] h2 = MathUtils.cumulativeHistogram(s2.getValues(), true);
final double[][] h3 = MathUtils.cumulativeHistogram(s3.getValues(), true);
if (settings.showFilterScoreHistograms) {
title = TITLE + " Cumul " + xLabel;
plot = new Plot(title, xLabel, "Frequency");
// Find limits
double[] xlimit = MathUtils.limits(h1[0]);
xlimit = MathUtils.limits(xlimit, h2[0]);
xlimit = MathUtils.limits(xlimit, h3[0]);
// Restrict using the inter-quartile range
if (filterCriteria[index].restrictRange) {
final double q1 = d.getPercentile(25);
final double q2 = d.getPercentile(75);
final double iqr = (q2 - q1) * 2.5;
xlimit[0] = MathUtils.max(xlimit[0], median - iqr);
xlimit[1] = MathUtils.min(xlimit[1], median + iqr);
}
plot.setLimits(xlimit[0], xlimit[1], 0, 1.05);
plot.addPoints(h1[0], h1[1], Plot.LINE);
plot.setColor(Color.red);
plot.addPoints(h2[0], h2[1], Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(h3[0], h3[1], Plot.LINE);
}
// Determine the maximum difference between the TP and FP
double maxx1 = 0;
double maxx2 = 0;
double max1 = 0;
double max2 = 0;
// We cannot compute the delta histogram, or use percentiles
if (s2.getN() == 0) {
upper = UpperLimit.ZERO;
lower = LowerLimit.ZERO;
}
final boolean requireLabel = (settings.showFilterScoreHistograms && filterCriteria[index].requireLabel);
if (requireLabel || upper.requiresDeltaHistogram() || lower.requiresDeltaHistogram()) {
if (s2.getN() != 0 && s3.getN() != 0) {
final LinearInterpolator li = new LinearInterpolator();
final PolynomialSplineFunction f1 = li.interpolate(h2[0], h2[1]);
final PolynomialSplineFunction f2 = li.interpolate(h3[0], h3[1]);
for (final double x : h1[0]) {
if (x < h2[0][0] || x < h3[0][0]) {
continue;
}
try {
final double v1 = f1.value(x);
final double v2 = f2.value(x);
final double diff = v2 - v1;
if (diff > 0) {
if (max1 < diff) {
max1 = diff;
maxx1 = x;
}
} else if (max2 > diff) {
max2 = diff;
maxx2 = x;
}
} catch (final OutOfRangeException ex) {
// Because we reached the end
break;
}
}
}
}
if (plot != null) {
// We use bins=1 on charts where we do not need a label
if (requireLabel) {
final String label = String.format("Max+ %s @ %s, Max- %s @ %s", MathUtils.rounded(max1), MathUtils.rounded(maxx1), MathUtils.rounded(max2), MathUtils.rounded(maxx2));
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
}
ImageJUtils.display(title, plot, wo);
}
// Now compute the bounds using the desired limit
double lowerBound;
double upperBound;
switch(lower) {
case MAX_NEGATIVE_CUMUL_DELTA:
// Switch to percentiles if we have no delta histogram
if (maxx2 < 0) {
lowerBound = maxx2;
break;
}
// fall-through
case ONE_PERCENT:
lowerBound = getPercentile(h2, 0.01);
break;
case MIN:
lowerBound = getPercentile(h2, 0.0);
break;
case ZERO:
lowerBound = 0;
break;
case HALF_MAX_JACCARD_VALUE:
lowerBound = getXValue(metric, jaccard, maxJaccard * 0.5);
break;
default:
throw new IllegalStateException("Missing lower limit method");
}
switch(upper) {
case MAX_POSITIVE_CUMUL_DELTA:
// Switch to percentiles if we have no delta histogram
if (maxx1 > 0) {
upperBound = maxx1;
break;
}
// fall-through
case NINETY_NINE_PERCENT:
upperBound = getPercentile(h2, 0.99);
break;
case NINETY_NINE_NINE_PERCENT:
upperBound = getPercentile(h2, 0.999);
break;
case ZERO:
upperBound = 0;
break;
case MAX_JACCARD2:
upperBound = getXValue(metric, jaccard, maxJaccard) * 2;
// System.out.printf("MaxJ = %.4f @ %.3f\n", maxJ, u / 2);
break;
default:
throw new IllegalStateException("Missing upper limit method");
}
final double min = getPercentile(h1, 0);
final double max = getPercentile(h1, 1);
return new double[] { lowerBound, upperBound, min, max };
}
use of uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder in project GDSC-SMLM by aherbert.
the class PcPalmMolecules method plot.
@Nullable
private static double[][] plot(DoubleData stats, String label, boolean integerBins) {
if (integerBins) {
// The histogram is not need for the return statement
new HistogramPlotBuilder(TITLE, stats, label).setMinBinWidth(1).show();
return null;
}
// Show a cumulative histogram so that the bin size is not relevant
final double[][] hist = MathUtils.cumulativeHistogram(stats.values(), false);
// Create the axes
final double[] xValues = hist[0];
final double[] yValues = hist[1];
// Plot
final String title = TITLE + " " + label;
final Plot plot = new Plot(title, label, "Frequency");
plot.addPoints(xValues, yValues, Plot.LINE);
ImageJUtils.display(title, plot);
return hist;
}
use of uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder in project GDSC-SMLM by aherbert.
the class SpotFinderPreview method run.
private void run(ImageProcessor ip, MaximaSpotFilter filter) {
if (refreshing) {
return;
}
currentSlice = imp.getCurrentSlice();
final Rectangle bounds = ip.getRoi();
// Crop to the ROI
FloatProcessor fp = ip.crop().toFloat(0, null);
float[] data = (float[]) fp.getPixels();
final int width = fp.getWidth();
final int height = fp.getHeight();
// Store the mean bias and gain of the region data.
// This is used to correctly overlay the filtered data on the original image.
double bias = 0;
double gain = 1;
boolean adjust = false;
// Set weights
final CameraModel cameraModel = fitConfig.getCameraModel();
if (!(cameraModel instanceof FakePerPixelCameraModel)) {
// This should be done on the normalised data
final float[] w = cameraModel.getNormalisedWeights(bounds);
filter.setWeights(w, width, height);
data = data.clone();
if (data.length < ip.getPixelCount()) {
adjust = true;
bias = MathUtils.sum(cameraModel.getBias(bounds)) / data.length;
gain = MathUtils.sum(cameraModel.getGain(bounds)) / data.length;
}
cameraModel.removeBiasAndGain(bounds, data);
}
final Spot[] spots = filter.rank(data, width, height);
data = filter.getPreprocessedData();
final int size = spots.length;
if (topNScrollBar != null) {
topNScrollBar.setMaximum(size);
selectScrollBar.setMaximum(size);
}
fp = new FloatProcessor(width, height, data);
final FloatProcessor out = new FloatProcessor(ip.getWidth(), ip.getHeight());
out.copyBits(ip, 0, 0, Blitter.COPY);
if (adjust) {
fp.multiply(gain);
fp.add(bias);
}
out.insert(fp, bounds.x, bounds.y);
final double min = fp.getMin();
final double max = fp.getMax();
out.setMinAndMax(min, max);
final Overlay o = new Overlay();
o.add(new ImageRoi(0, 0, out));
if (label != null) {
// Get results for frame
final Coordinate[] actual = ResultsMatchCalculator.getCoordinates(actualCoordinates, imp.getCurrentSlice());
final Coordinate[] predicted = new Coordinate[size];
for (int i = 0; i < size; i++) {
predicted[i] = new BasePoint(spots[i].x + bounds.x, spots[i].y + bounds.y);
}
// Compute assignments
final LocalList<FractionalAssignment> fractionalAssignments = new LocalList<>(3 * predicted.length);
final double matchDistance = settings.distance * fitConfig.getInitialPeakStdDev();
final RampedScore score = RampedScore.of(matchDistance, matchDistance * settings.lowerDistance / 100, false);
final double dmin = matchDistance * matchDistance;
final int nActual = actual.length;
final int nPredicted = predicted.length;
for (int j = 0; j < nPredicted; j++) {
// Centre in the middle of the pixel
final float x = predicted[j].getX() + 0.5f;
final float y = predicted[j].getY() + 0.5f;
// Any spots that match
for (int i = 0; i < nActual; i++) {
final double dx = (x - actual[i].getX());
final double dy = (y - actual[i].getY());
final double d2 = dx * dx + dy * dy;
if (d2 <= dmin) {
final double d = Math.sqrt(d2);
final double s = score.score(d);
if (s == 0) {
continue;
}
double distance = 1 - s;
if (distance == 0) {
// In the case of a match below the distance thresholds
// the distance will be 0. To distinguish between candidates all below
// the thresholds just take the closest.
// We know d2 is below dmin so we subtract the delta.
distance -= (dmin - d2);
}
// Store the match
fractionalAssignments.add(new ImmutableFractionalAssignment(i, j, distance, s));
}
}
}
final FractionalAssignment[] assignments = fractionalAssignments.toArray(new FractionalAssignment[0]);
// Compute matches
final RankedScoreCalculator calc = RankedScoreCalculator.create(assignments, nActual - 1, nPredicted - 1);
final boolean save = settings.showTP || settings.showFP;
final double[] calcScore = calc.score(nPredicted, settings.multipleMatches, save);
final ClassificationResult result = RankedScoreCalculator.toClassificationResult(calcScore, nActual);
// Compute AUC and max jaccard (and plot)
final double[][] curve = RankedScoreCalculator.getPrecisionRecallCurve(assignments, nActual, nPredicted);
final double[] precision = curve[0];
final double[] recall = curve[1];
final double[] jaccard = curve[2];
final double auc = AucCalculator.auc(precision, recall);
// Show scores
final String scoreLabel = String.format("Slice=%d, AUC=%s, R=%s, Max J=%s", imp.getCurrentSlice(), MathUtils.rounded(auc), MathUtils.rounded(result.getRecall()), MathUtils.rounded(MathUtils.maxDefault(0, jaccard)));
setLabel(scoreLabel);
// Plot
String title = TITLE + " Performance";
Plot plot = new Plot(title, "Spot Rank", "");
final double[] rank = SimpleArrayUtils.newArray(precision.length, 0, 1.0);
plot.setLimits(0, nPredicted, 0, 1.05);
plot.setColor(Color.blue);
plot.addPoints(rank, precision, Plot.LINE);
plot.setColor(Color.red);
plot.addPoints(rank, recall, Plot.LINE);
plot.setColor(Color.black);
plot.addPoints(rank, jaccard, Plot.LINE);
plot.setColor(Color.black);
plot.addLabel(0, 0, scoreLabel);
final WindowOrganiser windowOrganiser = new WindowOrganiser();
ImageJUtils.display(title, plot, 0, windowOrganiser);
title = TITLE + " Precision-Recall";
plot = new Plot(title, "Recall", "Precision");
plot.setLimits(0, 1, 0, 1.05);
plot.setColor(Color.red);
plot.addPoints(recall, precision, Plot.LINE);
plot.drawLine(recall[recall.length - 1], precision[recall.length - 1], recall[recall.length - 1], 0);
plot.setColor(Color.black);
plot.addLabel(0, 0, scoreLabel);
ImageJUtils.display(title, plot, 0, windowOrganiser);
windowOrganiser.tile();
// Create Rois for TP and FP
if (save) {
final double[] matchScore = RankedScoreCalculator.getMatchScore(calc.getScoredAssignments(), nPredicted);
int matches = 0;
for (int i = 0; i < matchScore.length; i++) {
if (matchScore[i] != 0) {
matches++;
}
}
if (settings.showTP) {
final float[] x = new float[matches];
final float[] y = new float[x.length];
int count = 0;
for (int i = 0; i < matchScore.length; i++) {
if (matchScore[i] != 0) {
final BasePoint p = (BasePoint) predicted[i];
x[count] = p.getX() + 0.5f;
y[count] = p.getY() + 0.5f;
count++;
}
}
addRoi(0, o, x, y, count, Color.green);
}
if (settings.showFP) {
final float[] x = new float[nPredicted - matches];
final float[] y = new float[x.length];
int count = 0;
for (int i = 0; i < matchScore.length; i++) {
if (matchScore[i] == 0) {
final BasePoint p = (BasePoint) predicted[i];
x[count] = p.getX() + 0.5f;
y[count] = p.getY() + 0.5f;
count++;
}
}
addRoi(0, o, x, y, count, Color.red);
}
}
} else {
final WindowOrganiser wo = new WindowOrganiser();
// Option to show the number of neighbours within a set pixel box radius
final int[] count = spotFilterHelper.countNeighbours(spots, width, height, settings.neighbourRadius);
// Show as histogram the totals...
new HistogramPlotBuilder(TITLE, StoredData.create(count), "Neighbours").setIntegerBins(true).setPlotLabel("Radius = " + settings.neighbourRadius).show(wo);
// TODO - Draw n=0, n=1 on the image overlay
final LUT lut = LutHelper.createLut(LutColour.FIRE_LIGHT);
// These are copied by the ROI
final float[] x = new float[1];
final float[] y = new float[1];
// Plot the intensity
final double[] intensity = new double[size];
final double[] rank = SimpleArrayUtils.newArray(size, 1, 1.0);
final int top = (settings.topN > 0) ? settings.topN : size;
final int size_1 = size - 1;
for (int i = 0; i < size; i++) {
intensity[i] = spots[i].intensity;
if (i < top) {
x[0] = spots[i].x + bounds.x + 0.5f;
y[0] = spots[i].y + bounds.y + 0.5f;
final Color c = LutHelper.getColour(lut, size_1 - i, size);
addRoi(0, o, x, y, 1, c, 2, 1);
}
}
final String title = TITLE + " Intensity";
final Plot plot = new Plot(title, "Rank", "Intensity");
plot.setColor(Color.blue);
plot.addPoints(rank, intensity, Plot.LINE);
if (settings.topN > 0 && settings.topN < size) {
plot.setColor(Color.magenta);
plot.drawLine(settings.topN, 0, settings.topN, intensity[settings.topN - 1]);
}
if (settings.select > 0 && settings.select < size) {
plot.setColor(Color.yellow);
final int index = settings.select - 1;
final double in = intensity[index];
plot.drawLine(settings.select, 0, settings.select, in);
x[0] = spots[index].x + bounds.x + 0.5f;
y[0] = spots[index].y + bounds.y + 0.5f;
final Color c = LutHelper.getColour(lut, size_1 - settings.select, size);
addRoi(0, o, x, y, 1, c, 3, 3);
plot.setColor(Color.black);
plot.addLabel(0, 0, "Selected spot intensity = " + MathUtils.rounded(in));
}
ImageJUtils.display(title, plot, 0, wo);
wo.tile();
}
imp.setOverlay(o);
}
use of uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder in project GDSC-SMLM by aherbert.
the class TcPalmAnalysis method plotHistogram.
/**
* Plot a histogram of the extracted statistic.
*
* @param show set to true to show the histogram
* @param wo the window organiser
* @param clusters the clusters
* @param name the name
* @param function the function to extract the plotted statistic
* @param predicate the predicate to filter the stream of data
* @param action the action to use on the histogram plot (to set non-standard options)
*/
private static void plotHistogram(boolean show, WindowOrganiser wo, LocalList<ClusterData> clusters, String name, ToDoubleFunction<ClusterData> function, DoublePredicate predicate, Consumer<HistogramPlotBuilder> action) {
if (!show) {
return;
}
final StoredData data = new StoredData(clusters.size());
DoubleStream stream = clusters.stream().mapToDouble(function);
if (predicate != null) {
stream = stream.filter(predicate);
}
stream.forEach(data::add);
final HistogramPlotBuilder builder = new HistogramPlotBuilder(TITLE, data, name);
action.accept(builder);
builder.show(wo);
}
Aggregations