use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class PcPalmFitting method fitClusteredModel.
/**
* Fits the correlation curve with r>0 to the clustered model using the estimated density and
* precision. Parameters must be fit within a tolerance of the starting values.
*
* @param gr the correlation curve
* @param sigmaS The estimated precision
* @param proteinDensity The estimated protein density
* @param resultColour the result colour
* @return The fitted parameters [precision, density, clusterRadius, clusterDensity]
*/
@Nullable
private double[] fitClusteredModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
final ClusteredModelFunctionGradient function = new ClusteredModelFunctionGradient();
clusteredModel = function;
ImageJUtils.log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", clusteredModel.getName(), sigmaS, proteinDensity * 1e6);
clusteredModel.setLogging(true);
for (int i = offset; i < gr[0].length; i++) {
// error)
if (gr[0][i] > sigmaS * settings.fitAboveEstimatedPrecision) {
clusteredModel.addPoint(gr[0][i], gr[1][i]);
}
}
double[] parameters;
// The model is: sigma, density, range, amplitude
final double[] initialSolution = new double[] { sigmaS, proteinDensity, sigmaS * 5, 1 };
// Constrain the fitting to be close to the estimated precision (sigmaS) and protein density.
// LVM fitting does not support constrained fitting so use a bounded optimiser.
final SumOfSquaresModelFunction clusteredModelMulti = new SumOfSquaresModelFunction(clusteredModel);
final double[] x = clusteredModelMulti.x;
// Put some bounds around the initial guess. Use the fitting tolerance (in %) if provided.
final double limit = (settings.fittingTolerance > 0) ? 1 + settings.fittingTolerance / 100 : 2;
final double[] lB = new double[] { initialSolution[0] / limit, initialSolution[1] / limit, 0, 0 };
// The amplitude and range should not extend beyond the limits of the g(r) curve.
final double[] uB = new double[] { initialSolution[0] * limit, initialSolution[1] * limit, MathUtils.max(x), MathUtils.max(gr[1]) };
ImageJUtils.log("Fitting %s using a bounded search: %s < precision < %s & %s < density < %s", clusteredModel.getName(), MathUtils.rounded(lB[0], 4), MathUtils.rounded(uB[0], 4), MathUtils.rounded(lB[1] * 1e6, 4), MathUtils.rounded(uB[1] * 1e6, 4));
final PointValuePair constrainedSolution = runBoundedOptimiser(initialSolution, lB, uB, clusteredModelMulti);
if (constrainedSolution == null) {
return null;
}
parameters = constrainedSolution.getPointRef();
int evaluations = boundedEvaluations;
// Refit using a LVM
if (settings.refitWithGradients) {
ImageJUtils.log("Re-fitting %s using a gradient optimisation", clusteredModel.getName());
final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
Optimum lvmSolution;
try {
// @formatter:off
final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(parameters).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {
@Override
public double[][] value(double[] point) {
return function.jacobian(point);
}
}).build();
// @formatter:on
lvmSolution = optimizer.optimize(problem);
evaluations += lvmSolution.getEvaluations();
final double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
if (ss < constrainedSolution.getValue()) {
ImageJUtils.log("Re-fitting %s improved the SS from %s to %s (-%s%%)", clusteredModel.getName(), MathUtils.rounded(constrainedSolution.getValue(), 4), MathUtils.rounded(ss, 4), MathUtils.rounded(100 * (constrainedSolution.getValue() - ss) / constrainedSolution.getValue(), 4));
parameters = lvmSolution.getPoint().toArray();
}
} catch (final TooManyIterationsException ex) {
ImageJUtils.log("Failed to re-fit %s: Too many iterations (%s)", clusteredModel.getName(), ex.getMessage());
} catch (final ConvergenceException ex) {
ImageJUtils.log("Failed to re-fit %s: %s", clusteredModel.getName(), ex.getMessage());
}
}
clusteredModel.setLogging(false);
// Ensure the width is positive
parameters[0] = Math.abs(parameters[0]);
double ss = 0;
final double[] obs = clusteredModel.getY();
final double[] exp = clusteredModel.value(parameters);
for (int i = 0; i < obs.length; i++) {
ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
}
final double totalSumSquares = MathUtils.getTotalSumOfSquares(clusteredModel.getY());
final double adjustedR2 = MathUtils.getAdjustedCoefficientOfDetermination(ss, totalSumSquares, clusteredModel.size(), parameters.length);
final double fitSigmaS = parameters[0];
final double fitProteinDensity = parameters[1];
// The radius of the cluster domain
final double domainRadius = parameters[2];
// The density of the cluster domain
final double domainDensity = parameters[3];
// This is from the PC-PALM paper. However that paper fits the g(r)protein exponential convolved
// in 2D with the g(r)PSF. In this method we have just fit the exponential
final double nCluster = 2 * domainDensity * Math.PI * domainRadius * domainRadius * fitProteinDensity;
final double e1 = parameterDrift(sigmaS, fitSigmaS);
final double e2 = parameterDrift(proteinDensity, fitProteinDensity);
ImageJUtils.log(" %s fit: SS = %f. Adj.R^2 = %f. %d evaluations", clusteredModel.getName(), ss, adjustedR2, evaluations);
ImageJUtils.log(" %s parameters:", clusteredModel.getName());
ImageJUtils.log(" Average precision = %s nm (%s%%)", MathUtils.rounded(fitSigmaS, 4), MathUtils.rounded(e1, 4));
ImageJUtils.log(" Average protein density = %s um^-2 (%s%%)", MathUtils.rounded(fitProteinDensity * 1e6, 4), MathUtils.rounded(e2, 4));
ImageJUtils.log(" Domain radius = %s nm", MathUtils.rounded(domainRadius, 4));
ImageJUtils.log(" Domain density = %s", MathUtils.rounded(domainDensity, 4));
ImageJUtils.log(" nCluster = %s", MathUtils.rounded(nCluster, 4));
// Check the fitted parameters are within tolerance of the initial estimates
valid2 = true;
if (settings.fittingTolerance > 0 && (Math.abs(e1) > settings.fittingTolerance || Math.abs(e2) > settings.fittingTolerance)) {
ImageJUtils.log(" Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%)," + " average protein density = %g um^-2 (%s%%)", clusteredModel.getName(), MathUtils.rounded(settings.fittingTolerance, 4), fitSigmaS, MathUtils.rounded(e1, 4), fitProteinDensity * 1e6, MathUtils.rounded(e2, 4));
valid2 = false;
}
// positive
if (domainRadius < fitSigmaS) {
ImageJUtils.log(" Failed to fit %s: Domain radius is smaller than the average precision (%s < %s)", clusteredModel.getName(), MathUtils.rounded(domainRadius, 4), MathUtils.rounded(fitSigmaS, 4));
valid2 = false;
}
if (domainDensity < 0) {
ImageJUtils.log(" Failed to fit %s: Domain density is negative (%s)", clusteredModel.getName(), MathUtils.rounded(domainDensity, 4));
valid2 = false;
}
if (adjustedR2 <= randomModelAdjustedR2) {
ImageJUtils.log(" Failed to fit %s - Adjusted r^2 has decreased %s%%", clusteredModel.getName(), MathUtils.rounded((100 * (randomModelAdjustedR2 - adjustedR2) / randomModelAdjustedR2), 4));
valid2 = false;
}
addResult(clusteredModel.getName(), resultColour, valid2, fitSigmaS, fitProteinDensity, domainRadius, domainDensity, nCluster, -1, adjustedR2);
return parameters;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable 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.annotation.Nullable in project GDSC-SMLM by aherbert.
the class PcPalmMolecules method fitGaussian.
@Nullable
private static double[] fitGaussian(float[] x, float[] y) {
final WeightedObservedPoints obs = new WeightedObservedPoints();
for (int i = 0; i < x.length; i++) {
obs.add(x[i], y[i]);
}
final Collection<WeightedObservedPoint> observations = obs.toList();
final GaussianCurveFitter fitter = GaussianCurveFitter.create().withMaxIterations(2000);
final GaussianCurveFitter.ParameterGuesser guess = new GaussianCurveFitter.ParameterGuesser(observations);
double[] initialGuess = null;
try {
initialGuess = guess.guess();
return fitter.withStartPoint(initialGuess).fit(observations);
} catch (final TooManyEvaluationsException ex) {
// Use the initial estimate
return initialGuess;
} catch (final Exception ex) {
// Just in case there is another exception type, or the initial estimate failed
return null;
}
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class PcPalmAnalysis method computeAutoCorrelationCurveFft.
/**
* Compute the auto-correlation curve using FFT. Computes the correlation image and then samples
* the image at radii up to the specified length to get the average correlation at a given radius.
*
* @param im the image
* @param wp the weighted image
* @param maxRadius the max radius
* @param nmPerPixel the nm per pixel
* @param density the density
* @return { distances[], gr[], gr_se[] }
*/
@Nullable
private double[][] computeAutoCorrelationCurveFft(ImageProcessor im, ImageProcessor wp, int maxRadius, double nmPerPixel, double density) {
log("Performing FFT correlation");
final FloatProcessor corrIm = computeAutoCorrelationFft(im);
final FloatProcessor corrW = computeAutoCorrelationFft(wp);
if (corrIm == null || corrW == null) {
error("Unable to perform Fourier transform");
return null;
}
final int centre = corrIm.getHeight() / 2;
final Rectangle crop = new Rectangle(centre - maxRadius, centre - maxRadius, maxRadius * 2, maxRadius * 2);
if (settings.showCorrelationImages) {
displayCorrelation(corrIm, "Image correlation FFT", crop);
displayCorrelation(corrW, "Window correlation FFT", crop);
}
log(" Normalising correlation");
final FloatProcessor correlation = normaliseCorrelation(corrIm, corrW, density);
if (settings.showCorrelationImages) {
displayCorrelation(correlation, "Normalised correlation FFT", crop);
}
return computeRadialAverage(maxRadius, nmPerPixel, correlation);
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class TraceMatchCalculator method extractPulses.
@Nullable
private static Pulse[] extractPulses(MemoryPeakResults results) {
if (results == null) {
return null;
}
final Pulse[] pulses = new Pulse[results.size()];
final Counter i = new Counter();
results.forEach(DistanceUnit.PIXEL, (XyrResultProcedure) (x, y, result) -> pulses[i.getAndIncrement()] = new Pulse(x, y, result.getFrame(), result.getEndFrame()));
return pulses;
}
Aggregations