use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class DriftCalculator method calculateUsingFrames.
/**
* Calculates drift using images from N consecutive frames aligned to the overall image.
*
* @param results the results
* @param limits the limits
* @param reconstructionSize the reconstruction size
* @return the drift { dx[], dy[] }
*/
@Nullable
private double[][] calculateUsingFrames(MemoryPeakResults results, int[] limits, int reconstructionSize) {
// Extract the localisations into blocks of N consecutive frames
final BlockPeakResultProcedure p = new BlockPeakResultProcedure(settings);
results.sort();
results.forEach(p);
final List<List<Localisation>> blocks = p.blocks;
if (blocks.size() <= 1) {
tracker.log("ERROR : Require at least 2 images for drift calculation");
return null;
}
// Check the final block has enough localisations
final List<Localisation> nextBlock = p.nextBlock;
if (nextBlock.size() < settings.minimimLocalisations) {
blocks.remove(blocks.size() - 1);
if (blocks.size() <= 1) {
tracker.log("ERROR : Require at least 2 images for drift calculation");
return null;
}
final List<Localisation> combinedBlock = blocks.get(blocks.size() - 1);
combinedBlock.addAll(nextBlock);
}
// Find the average time point for each block
final int[] blockT = new int[blocks.size()];
int time = 0;
for (final List<Localisation> block : blocks) {
long sum = 0;
for (final Localisation r : block) {
sum += r.time;
}
blockT[time++] = (int) (sum / block.size());
}
// Calculate a scale to use when constructing the images for alignment
final Rectangle bounds = results.getBounds(true);
final ResultsImageSettings.Builder builder = ResultsImageSettings.newBuilder().setImageSizeMode(ResultsImageSizeMode.IMAGE_SIZE).setImageSize(reconstructionSize);
final float scale = ImagePeakResultsFactory.getScale(builder, bounds, 1);
executor = Executors.newFixedThreadPool(Prefs.getThreads());
final double[] dx = new double[limits[1] + 1];
final double[] dy = new double[dx.length];
final double[] originalDriftTimePoints = getOriginalDriftTimePoints(dx, blockT);
lastdx = null;
final double smoothing = updateSmoothingParameter(originalDriftTimePoints);
double change = calculateDriftUsingFrames(blocks, blockT, bounds, scale, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
if (Double.isNaN(change) || tracker.isEnded()) {
return null;
}
plotDrift(limits, dx, dy);
ImageJUtils.log("Drift Calculator : Initial drift " + MathUtils.rounded(change));
for (int i = 1; i <= settings.maxIterations; i++) {
change = calculateDriftUsingFrames(blocks, blockT, bounds, scale, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
if (Double.isNaN(change)) {
return null;
}
plotDrift(limits, dx, dy);
if (converged(i, change, getTotalDrift(dx, dy, originalDriftTimePoints))) {
break;
}
}
if (tracker.isEnded()) {
return null;
}
plotDrift(limits, dx, dy);
return new double[][] { dx, dy };
}
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 ssTotal = 0;
for (final double o : obs) {
ssTotal += (o - mean) * (o - mean);
}
// This is true if the weights are 1
final double 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;
}
Aggregations