use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.
the class BenchmarkFilterAnalysis method depthAnalysis.
/**
* Depth analysis.
*
* @param allAssignments The assignments generated from running the filter (or null)
* @param filter the filter
* @return the assignments
*/
@Nullable
private List<FractionalAssignment[]> depthAnalysis(List<FractionalAssignment[]> allAssignments, DirectFilter filter) {
if (!settings.depthRecallAnalysis || simulationParameters.fixedDepth) {
return null;
}
// Build a histogram of the number of spots at different depths
final double[] depths = fitResultData.depthStats.getValues();
double[] limits = MathUtils.limits(depths);
final int bins = HistogramPlot.getBinsSqrtRule(depths.length);
final double[][] h1 = HistogramPlot.calcHistogram(depths, limits[0], limits[1], bins);
final double[][] h2 = HistogramPlot.calcHistogram(fitResultData.depthFitStats.getValues(), limits[0], limits[1], bins);
// manually to get the results that pass.
if (allAssignments == null) {
allAssignments = getAssignments(filter);
}
double[] depths2 = new double[results.size()];
int count = 0;
for (final FractionalAssignment[] assignments : allAssignments) {
if (assignments == null) {
continue;
}
for (final FractionalAssignment assignment : assignments) {
final CustomFractionalAssignment c = (CustomFractionalAssignment) assignment;
depths2[count++] = c.peak.getZPosition();
}
}
depths2 = Arrays.copyOf(depths2, count);
// Build a histogram using the same limits
final double[][] h3 = HistogramPlot.calcHistogram(depths2, limits[0], limits[1], bins);
// Convert pixel depth to nm
for (int i = 0; i < h1[0].length; i++) {
h1[0][i] *= simulationParameters.pixelPitch;
}
limits[0] *= simulationParameters.pixelPitch;
limits[1] *= simulationParameters.pixelPitch;
// Produce a histogram of the number of spots at each depth
final String title1 = TITLE + " Depth Histogram";
final Plot plot1 = new Plot(title1, "Depth (nm)", "Frequency");
plot1.setLimits(limits[0], limits[1], 0, MathUtils.max(h1[1]));
plot1.setColor(Color.black);
plot1.addPoints(h1[0], h1[1], Plot.BAR);
plot1.addLabel(0, 0, "Black = Spots; Blue = Fitted; Red = Filtered");
plot1.setColor(Color.blue);
plot1.addPoints(h1[0], h2[1], Plot.BAR);
plot1.setColor(Color.red);
plot1.addPoints(h1[0], h3[1], Plot.BAR);
plot1.setColor(Color.magenta);
ImageJUtils.display(title1, plot1, wo);
// Interpolate
final double halfBinWidth = (h1[0][1] - h1[0][0]) * 0.5;
// Remove final value of the histogram as this is at the upper limit of the range (i.e. count
// zero)
h1[0] = Arrays.copyOf(h1[0], h1[0].length - 1);
h1[1] = Arrays.copyOf(h1[1], h1[0].length);
h2[1] = Arrays.copyOf(h2[1], h1[0].length);
h3[1] = Arrays.copyOf(h3[1], h1[0].length);
// TODO : Fix the smoothing since LOESS sometimes does not work.
// Perhaps allow configuration of the number of histogram bins and the smoothing bandwidth.
// Use minimum of 3 points for smoothing
// Ensure we use at least x% of data
final double bandwidth = Math.max(3.0 / h1[0].length, 0.15);
final LoessInterpolator loess = new LoessInterpolator(bandwidth, 1);
final PolynomialSplineFunction spline1 = loess.interpolate(h1[0], h1[1]);
final PolynomialSplineFunction spline2 = loess.interpolate(h1[0], h2[1]);
final PolynomialSplineFunction spline3 = loess.interpolate(h1[0], h3[1]);
// Use a second interpolator in case the LOESS fails
final LinearInterpolator lin = new LinearInterpolator();
final PolynomialSplineFunction spline1b = lin.interpolate(h1[0], h1[1]);
final PolynomialSplineFunction spline2b = lin.interpolate(h1[0], h2[1]);
final PolynomialSplineFunction spline3b = lin.interpolate(h1[0], h3[1]);
// Increase the number of points to show a smooth curve
final double[] points = new double[bins * 5];
limits = MathUtils.limits(h1[0]);
final double interval = (limits[1] - limits[0]) / (points.length - 1);
final double[] v = new double[points.length];
final double[] v2 = new double[points.length];
final double[] v3 = new double[points.length];
for (int i = 0; i < points.length - 1; i++) {
points[i] = limits[0] + i * interval;
v[i] = getSplineValue(spline1, spline1b, points[i]);
v2[i] = getSplineValue(spline2, spline2b, points[i]);
v3[i] = getSplineValue(spline3, spline3b, points[i]);
points[i] += halfBinWidth;
}
// Final point on the limit of the spline range
final int ii = points.length - 1;
v[ii] = getSplineValue(spline1, spline1b, limits[1]);
v2[ii] = getSplineValue(spline2, spline2b, limits[1]);
v3[ii] = getSplineValue(spline3, spline3b, limits[1]);
points[ii] = limits[1] + halfBinWidth;
// Calculate recall
for (int i = 0; i < v.length; i++) {
v2[i] = v2[i] / v[i];
v3[i] = v3[i] / v[i];
}
final double halfSummaryDepth = settings.summaryDepth * 0.5;
final String title2 = TITLE + " Depth Histogram (normalised)";
final Plot plot2 = new Plot(title2, "Depth (nm)", "Recall");
plot2.setLimits(limits[0] + halfBinWidth, limits[1] + halfBinWidth, 0, MathUtils.min(1, MathUtils.max(v2)));
plot2.setColor(Color.black);
plot2.addLabel(0, 0, "Blue = Fitted; Red = Filtered");
plot2.setColor(Color.blue);
plot2.addPoints(points, v2, Plot.LINE);
plot2.setColor(Color.red);
plot2.addPoints(points, v3, Plot.LINE);
plot2.setColor(Color.magenta);
if (-halfSummaryDepth - halfBinWidth >= limits[0]) {
plot2.drawLine(-halfSummaryDepth, 0, -halfSummaryDepth, getSplineValue(spline3, spline3b, -halfSummaryDepth - halfBinWidth) / getSplineValue(spline1, spline1b, -halfSummaryDepth - halfBinWidth));
}
if (halfSummaryDepth - halfBinWidth <= limits[1]) {
plot2.drawLine(halfSummaryDepth, 0, halfSummaryDepth, getSplineValue(spline3, spline3b, halfSummaryDepth - halfBinWidth) / getSplineValue(spline1, spline1b, halfSummaryDepth - halfBinWidth));
}
ImageJUtils.display(title2, plot2, wo);
return allAssignments;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class CameraModelFisherInformationAnalysis method createExponents.
@Nullable
private double[] createExponents() {
final int n = 1 + Math.max(0, settings.getSubDivisions());
final double h = 1.0 / n;
final double minExp = settings.getMinExponent();
final double maxExp = settings.getMaxExponent();
final double size = (maxExp - minExp) * n + 1;
if (size > 100) {
final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Number of exponents is " + Math.ceil(size) + ". OK to continue?");
gd.showDialog();
if (gd.wasCanceled()) {
return null;
}
}
final TDoubleArrayList list = new TDoubleArrayList();
for (int i = 0; ; i++) {
final double e = minExp + i * h;
list.add(e);
if (e >= settings.getMaxExponent()) {
break;
}
}
return list.toArray();
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class BlinkEstimator method fit.
/**
* Fit the dark time to counts of molecules curve. Only use the first n fitted points.
*
* <p>Calculates:<br> N = The number of photoblinking molecules in the sample<br> nBlink = The
* average number of blinks per flourophore<br> tOff = The off-time
*
* @param td The dark time
* @param ntd The counts of molecules
* @param numberOfFittedPointsSetting the number of fitted points
* @param log Write the fitting results to the ImageJ log window
* @return The fitted parameters [N, nBlink, tOff], or null if no fit was possible
*/
@Nullable
public double[] fit(double[] td, double[] ntd, int numberOfFittedPointsSetting, boolean log) {
blinkingModel = new BlinkingFunction();
blinkingModel.setLogging(true);
for (int i = 0; i < numberOfFittedPointsSetting; i++) {
blinkingModel.addPoint(td[i], ntd[i]);
}
final LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(INITIAL_STEP_BOUND_FACTOR, COST_RELATIVE_TOLERANCE, PAR_RELATIVE_TOLERANCE, ORTHO_TOLERANCE, THRESHOLD);
try {
final double[] obs = blinkingModel.getY();
// @formatter:off
final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(1000).start(new double[] { ntd[0], 0.1, td[1] }).target(obs).weight(new DiagonalMatrix(blinkingModel.getWeights())).model(blinkingModel, blinkingModel::jacobian).build();
// @formatter:on
blinkingModel.setLogging(false);
final Optimum optimum = optimiser.optimize(problem);
final double[] parameters = optimum.getPoint().toArray();
double mean = 0;
for (final double d : obs) {
mean += d;
}
mean /= obs.length;
double ssResiduals = 0;
double ssTotal = 0;
for (int i = 0; i < obs.length; i++) {
ssTotal += (obs[i] - mean) * (obs[i] - mean);
}
// This is true if the weights are 1
ssResiduals = optimum.getResiduals().dotProduct(optimum.getResiduals());
r2 = 1 - ssResiduals / ssTotal;
adjustedR2 = getAdjustedCoefficientOfDetermination(ssResiduals, ssTotal, obs.length, parameters.length);
if (log) {
ImageJUtils.log(" Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, MathUtils.rounded(r2, 4), MathUtils.rounded(adjustedR2, 4));
ImageJUtils.log(" N=%s, nBlink=%s, tOff=%s (%s frames)", MathUtils.rounded(parameters[0], 4), MathUtils.rounded(parameters[1], 4), MathUtils.rounded(parameters[2], 4), MathUtils.rounded(parameters[2] / msPerFrame, 4));
}
return parameters;
} catch (final TooManyIterationsException ex) {
if (log) {
ImageJUtils.log(" Failed to fit %d points: Too many iterations: (%s)", blinkingModel.size(), ex.getMessage());
}
} catch (final ConvergenceException ex) {
if (log) {
ImageJUtils.log(" Failed to fit %d points", blinkingModel.size());
}
}
return null;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class BenchmarkFilterAnalysis method scoreAnalysis.
/**
* Score analysis.
*
* @param allAssignments The assignments generated from running the filter (or null)
* @param filter the filter
* @return the assignments
*/
@Nullable
private ArrayList<FractionalAssignment[]> scoreAnalysis(ArrayList<FractionalAssignment[]> allAssignments, DirectFilter filter) {
if (!settings.scoreAnalysis) {
return null;
}
// Build a histogram of the fitted spots that were available to be scored
final double[] signal = fitResultData.signalFactorStats.getValues();
final double[] distance = fitResultData.distanceStats.getValues();
double[] limits1;
if (fitSignalFactor > 0 && settings.upperSignalFactor > 0) {
final double range = fitSignalFactor * settings.upperSignalFactor / 100.0;
limits1 = new double[] { -range, range };
} else {
limits1 = MathUtils.limits(signal);
// Prevent the auto-range being too big
final double bound = 3;
if (limits1[0] < -bound) {
limits1[0] = -bound;
}
if (limits1[1] > bound) {
limits1[1] = bound;
}
}
double[] limits2;
if (spotFitResults.distanceInPixels > 0 && settings.upperMatchDistance > 0) {
final double range = simulationParameters.pixelPitch * spotFitResults.distanceInPixels * settings.upperMatchDistance / 100.0;
limits2 = new double[] { 0, range };
} else {
limits2 = MathUtils.limits(distance);
}
final int bins = HistogramPlot.getBinsSqrtRule(signal.length);
final double[][] h1 = HistogramPlot.calcHistogram(signal, limits1[0], limits1[1], bins);
final double[][] h2 = HistogramPlot.calcHistogram(distance, limits2[0], limits2[1], bins);
// Run the filter manually to get the results that pass.
if (allAssignments == null) {
allAssignments = getAssignments(filter);
}
double[] signal2 = new double[results.size()];
double[] distance2 = new double[results.size()];
int count = 0;
double sumSignal = 0;
double sumDistance = 0;
for (final FractionalAssignment[] assignments : allAssignments) {
if (assignments == null) {
continue;
}
for (int i = 0; i < assignments.length; i++) {
final CustomFractionalAssignment c = (CustomFractionalAssignment) assignments[i];
sumDistance += distance2[count] = c.distToTarget;
sumSignal += signal2[count] = c.getSignalFactor();
count++;
}
}
signal2 = Arrays.copyOf(signal2, count);
distance2 = Arrays.copyOf(distance2, count);
// Build a histogram using the same limits
final double[][] h1b = HistogramPlot.calcHistogram(signal2, limits1[0], limits1[1], bins);
final double[][] h2b = HistogramPlot.calcHistogram(distance2, limits2[0], limits2[1], bins);
// Since the distance and signal factor are computed for all fits (single, multi, doublet)
// there will be far more of them so we normalise and just plot the histogram profile.
double s1 = 0;
double s2 = 0;
double s1b = 0;
double s2b = 0;
for (int i = 0; i < h1b[0].length; i++) {
s1 += h1[1][i];
s2 += h2[1][i];
s1b += h1b[1][i];
s2b += h2b[1][i];
}
for (int i = 0; i < h1b[0].length; i++) {
h1[1][i] /= s1;
h2[1][i] /= s2;
h1b[1][i] /= s1b;
h2b[1][i] /= s2b;
}
// Draw distance histogram first
final String title2 = TITLE + " Distance Histogram";
final Plot plot2 = new Plot(title2, "Distance (nm)", "Frequency");
plot2.setLimits(limits2[0], limits2[1], 0, MathUtils.maxDefault(MathUtils.max(h2[1]), h2b[1]));
plot2.setColor(Color.black);
plot2.addLabel(0, 0, String.format("Blue = Fitted (%s); Red = Filtered (%s)", MathUtils.rounded(fitResultData.distanceStats.getMean()), MathUtils.rounded(sumDistance / count)));
plot2.setColor(Color.blue);
plot2.addPoints(h2[0], h2[1], Plot.BAR);
plot2.setColor(Color.red);
plot2.addPoints(h2b[0], h2b[1], Plot.BAR);
ImageJUtils.display(title2, plot2, wo);
// Draw signal factor histogram
final String title1 = TITLE + " Signal Factor Histogram";
final Plot plot1 = new Plot(title1, "Signal Factor", "Frequency");
plot1.setLimits(limits1[0], limits1[1], 0, MathUtils.maxDefault(MathUtils.max(h1[1]), h1b[1]));
plot1.setColor(Color.black);
plot1.addLabel(0, 0, String.format("Blue = Fitted (%s); Red = Filtered (%s)", MathUtils.rounded(fitResultData.signalFactorStats.getMean()), MathUtils.rounded(sumSignal / count)));
plot1.setColor(Color.blue);
plot1.addPoints(h1[0], h1[1], Plot.BAR);
plot1.setColor(Color.red);
plot1.addPoints(h1b[0], h1b[1], Plot.BAR);
ImageJUtils.display(title1, plot1, wo);
return allAssignments;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class BenchmarkFilterAnalysis method scoreFilters.
@Nullable
private FilterScoreResult[] scoreFilters(FilterSet filterSet, boolean createTextResult) {
if (filterSet.size() == 0) {
return null;
}
initialiseScoring(filterSet);
FilterScoreResult[] scoreResults = new FilterScoreResult[filterSet.size()];
if (scoreResults.length == 1) {
// No need to multi-thread this
scoreResults[0] = scoreFilter((DirectFilter) filterSet.getFilters().get(0), defaultMinimalFilter, createTextResult, coordinateStore);
} else {
// Multi-thread score all the result
final int nThreads = getThreads(scoreResults.length);
final BlockingQueue<ScoreJob> jobs = new ArrayBlockingQueue<>(nThreads * 2);
final List<Thread> threads = new LinkedList<>();
final Ticker ticker = ImageJUtils.createTicker(scoreResults.length, nThreads, "Scoring Filters");
for (int i = 0; i < nThreads; i++) {
final ScoreWorker worker = new ScoreWorker(jobs, scoreResults, createTextResult, (coordinateStore == null) ? null : coordinateStore.newInstance(), ticker);
final Thread t = new Thread(worker);
threads.add(t);
t.start();
}
int index = 0;
for (final Filter filter : filterSet.getFilters()) {
if (IJ.escapePressed()) {
break;
}
put(jobs, new ScoreJob((DirectFilter) filter, index++));
}
// Finish all the worker threads by passing in a null job
for (int i = 0; i < threads.size(); i++) {
put(jobs, new ScoreJob(null, -1));
}
// Wait for all to finish
for (int i = 0; i < threads.size(); i++) {
try {
threads.get(i).join();
} catch (final InterruptedException ex) {
Logger.getLogger(BenchmarkFilterAnalysis.class.getName()).log(Level.WARNING, "Interrupted!", ex);
Thread.currentThread().interrupt();
throw new ConcurrentRuntimeException("Unexpected interruption", ex);
}
}
threads.clear();
ImageJUtils.finished();
// In case the threads were interrupted
if (ImageJUtils.isInterrupted()) {
scoreResults = null;
}
}
finishScoring();
return scoreResults;
}
Aggregations