use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class DriftCalculator method calculateUsingMarkers.
/**
* Calculates drift using the feducial markers within ROI.
*
* <p>Adapted from the drift calculation method in QuickPALM.
*
* @param results the results
* @param limits the limits
* @param rois the rois
* @return the drift { dx[], dy[] }
*/
@Nullable
private double[][] calculateUsingMarkers(MemoryPeakResults results, int[] limits, Roi[] rois) {
final Spot[][] roiSpots = findSpots(results, rois, limits);
// Check we have enough data
if (roiSpots.length == 0) {
IJ.error("No peak fit results in the selected ROIs");
return null;
}
final double[] dx = new double[limits[1] + 1];
final double[] dy = new double[dx.length];
final double[] sum = new double[roiSpots.length];
final double[] weights = calculateWeights(roiSpots, dx.length, sum);
final double smoothing = updateSmoothingParameter(weights);
lastdx = null;
double change = calculateDriftUsingMarkers(roiSpots, weights, sum, dx, dy, smoothing, settings.iterations);
if (Double.isNaN(change) || tracker.isEnded()) {
return null;
}
ImageJUtils.log("Drift Calculator : Initial drift " + MathUtils.rounded(change));
for (int i = 1; i <= settings.maxIterations; i++) {
change = calculateDriftUsingMarkers(roiSpots, weights, sum, dx, dy, smoothing, settings.iterations);
if (Double.isNaN(change)) {
return null;
}
if (converged(i, change, getTotalDrift(dx, dy, weights))) {
break;
}
}
if (tracker.isEnded()) {
return null;
}
interpolate(dx, dy, weights);
plotDrift(limits, dx, dy);
saveDrift(weights, dx, dy);
return new double[][] { dx, dy };
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class MultiPathFilter method acceptAnyDoublet.
/**
* Check any new and all existing results within the multi-doublet fit results are valid. Returns
* the new results. Coordinate shift filter is disabled for the doublet results.
*
* <p>New results and validated candidates that fail the primary filter can be filtered using the
* minimal filter and sent to the store. The store can be used to determine if a fit for a
* different candidate has been performed already.
*
* @param multiPathResult the multi path result
* @param validateCandidates Set to true to validate the candidates
* @param store the store
* @param candidateId the candidate id
* @return The new results that pass the filter
*/
@Nullable
private PreprocessedPeakResult[] acceptAnyDoublet(final MultiPathFitResult multiPathResult, boolean validateCandidates, SelectedResultStore store, final int candidateId) {
final FitResult multiDoubletFitResult = multiPathResult.getMultiDoubletFitResult();
if (multiDoubletFitResult == null || multiDoubletFitResult.getResults() == null) {
return null;
}
// Doublets may drift further than single spot candidates.
// We must validate the doublet spot without shift filtering.
// Note: Only disable shift for the doublet results.
// doublets = len(multi-doublet) - len(multi) + 1
final PreprocessedPeakResult[] results = multiDoubletFitResult.getResults();
final int nDoublets = results.length - multiPathResult.getMultiFitResult().getResults().length + 1;
setupFilter(FilterValidationOption.NO_SHIFT);
validationResults = new int[results.length];
for (int i = 0; i < nDoublets; i++) {
validationResults[i] = filter.validate(results[i]);
}
restoreFilterState();
for (int i = nDoublets; i < results.length; i++) {
validationResults[i] = filter.validate(results[i]);
}
return acceptAnyInternal(candidateId, multiDoubletFitResult, validateCandidates, store);
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class MultiPathFilter method score.
/**
* Score the assignments (TP/FP) and then clear the list.
*
* @param assignments The assignments
* @param score Scores array to accumulate TP/FP scores
* @param predicted The number of predictions
* @param save Set to true to save the scored assignments
* @param actual The number of actual results in the frame
* @return the fractional assignments
*/
@Nullable
private static FractionalAssignment[] score(final ArrayList<FractionalAssignment> assignments, final double[] score, final int predicted, boolean save, int actual) {
if (assignments.isEmpty()) {
return null;
}
final FractionalAssignment[] tmp = assignments.toArray(new FractionalAssignment[0]);
final RankedScoreCalculator calc = RankedScoreCalculator.create(tmp, actual, predicted);
final double[] result = calc.score(predicted, false, save);
score[0] += result[0];
score[1] += result[1];
score[2] += result[2];
score[3] += result[3];
assignments.clear();
return calc.getScoredAssignments();
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc by aherbert.
the class FindFociBaseProcessor method findCentreGaussianFit.
/**
* Finds the centre of the image using a 2D Gaussian fit to projection along the Z-axis.
*
* @param subImage the sub image
* @param dimensions the dimensions
* @param projectionMethod (0) Average value; (1) Maximum value
* @return the centre of the image
*/
@Nullable
private int[] findCentreGaussianFit(float[] subImage, int[] dimensions, int projectionMethod) {
if (FindFoci_PlugIn.IS_GAUSSIAN_FIT_ENABLED < 1) {
return null;
}
final int blockSize = dimensions[0] * dimensions[1];
final float[] projection = new float[blockSize];
if (projectionMethod == 1) {
// Maximum value
for (int z = dimensions[2]; z-- > 0; ) {
int index = blockSize * z;
for (int i = 0; i < blockSize; i++, index++) {
if (projection[i] < subImage[index]) {
projection[i] = subImage[index];
}
}
}
} else {
// Average value
for (int z = dimensions[2]; z-- > 0; ) {
int index = blockSize * z;
for (int i = 0; i < blockSize; i++, index++) {
projection[i] += subImage[index];
}
}
for (int i = blockSize; i-- > 0; ) {
projection[i] /= dimensions[2];
}
}
final GaussianFit_PlugIn gf = new GaussianFit_PlugIn();
final double[] fitParams = gf.fit(projection, dimensions[0], dimensions[1]);
int[] centre = null;
if (fitParams != null) {
// Find the centre of mass along the z-axis
centre = convertCentre(new double[] { fitParams[2], fitParams[3] });
// Use the centre of mass along the projection axis
double com = 0;
long sum = 0;
for (int z = dimensions[2]; z-- > 0; ) {
final int index = blockSize * z;
final float value = subImage[index];
if (value > 0) {
com += z * value;
sum += value;
}
}
if (sum == 0) {
centre[2] = dimensions[2] / 2;
} else {
// Avoid clipping
centre[2] = MathUtils.clip(0, dimensions[2] - 1, round(com / sum));
}
}
return centre;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc by aherbert.
the class FindFociBaseProcessor method getSortedMaxpoints.
/**
* Find all local maxima (irrespective whether they finally qualify as maxima or not).
*
* @param pixels The image to be analyzed
* @param maxima the maxima
* @param types A byte image, same size as ip, where the maximum points are marked as MAXIMUM
* @param globalMin The image global minimum
* @param threshold The threshold below which no pixels are processed.
* @return Maxima sorted by value.
*/
@Nullable
protected Coordinate[] getSortedMaxpoints(Object pixels, int[] maxima, byte[] types, float globalMin, float threshold) {
final ArrayList<Coordinate> maxpoints = new ArrayList<>(500);
// working list for expanding local plateaus
final IntArrayList pointList = new IntArrayList();
int id = 0;
final int[] xyz = new int[3];
setPixels(pixels);
if (is2D()) {
for (int i = maxxByMaxyByMaxz; i-- > 0; ) {
if ((types[i] & (EXCLUDED | MAX_AREA | PLATEAU | NOT_MAXIMUM)) != 0) {
continue;
}
final float v = getf(i);
if (v < threshold || v == globalMin) {
continue;
}
getXy(i, xyz);
final int x = xyz[0];
final int y = xyz[1];
// Check whether we have a local maximum.
final boolean isInnerXy = (y != 0 && y != ylimit) && (x != 0 && x != xlimit);
boolean isMax = true;
boolean equalNeighbour = false;
for (int d = 8; d-- > 0; ) {
if (isInnerXy || isWithinXy(x, y, d)) {
final float vNeighbor = getf(i + offset[d]);
if (vNeighbor > v) {
isMax = false;
break;
} else if (vNeighbor == v) {
// Neighbour is equal, this is a potential plateau maximum
equalNeighbour = true;
} else {
// This is lower so cannot be a maxima
types[i + offset[d]] |= NOT_MAXIMUM;
}
}
}
if (isMax) {
id++;
if (id >= searchCapacity) {
log(() -> "The number of potential maxima exceeds the search capacity: " + searchCapacity + ". Try using a denoising/smoothing filter or increase the capacity.");
return null;
}
if (equalNeighbour) {
// Search the local area marking all equal neighbour points as maximum
if (!expandMaximum(maxima, types, globalMin, threshold, i, v, id, maxpoints, pointList)) {
// Not a true maximum, ignore this
id--;
}
} else {
types[i] |= MAXIMUM | MAX_AREA;
maxima[i] = id;
maxpoints.add(new Coordinate(getIndex(x, y), id, v));
}
}
}
} else {
for (int i = maxxByMaxyByMaxz; i-- > 0; ) {
if ((types[i] & (EXCLUDED | MAX_AREA | PLATEAU | NOT_MAXIMUM)) != 0) {
continue;
}
final float v = getf(i);
if (v < threshold || v == globalMin) {
continue;
}
getXyz(i, xyz);
final int x = xyz[0];
final int y = xyz[1];
final int z = xyz[2];
// Check whether we have a local maximum.
final boolean isInnerXy = (y != 0 && y != ylimit) && (x != 0 && x != xlimit);
final boolean isInnerXyz = (zlimit == 0) ? isInnerXy : isInnerXy && (z != 0 && z != zlimit);
boolean isMax = true;
boolean equalNeighbour = false;
for (int d = 26; d-- > 0; ) {
if (isInnerXyz || (isInnerXy && isWithinZ(z, d)) || isWithinXyz(x, y, z, d)) {
final float vNeighbor = getf(i + offset[d]);
if (vNeighbor > v) {
isMax = false;
break;
} else if (vNeighbor == v) {
// Neighbour is equal, this is a potential plateau maximum
equalNeighbour = true;
} else {
// This is lower so cannot be a maxima
types[i + offset[d]] |= NOT_MAXIMUM;
}
}
}
if (isMax) {
id++;
if (id >= searchCapacity) {
log(() -> "The number of potential maxima exceeds the search capacity: " + searchCapacity + ". Try using a denoising/smoothing filter or increase the capacity.");
return null;
}
if (equalNeighbour) {
// Search the local area marking all equal neighbour points as maximum
if (!expandMaximum(maxima, types, globalMin, threshold, i, v, id, maxpoints, pointList)) {
// Not a true maximum, ignore this
id--;
}
} else {
types[i] |= MAXIMUM | MAX_AREA;
maxima[i] = id;
maxpoints.add(new Coordinate(getIndex(x, y, z), id, v));
}
}
}
}
if (ImageJUtils.isInterrupted()) {
return null;
}
for (int i = maxxByMaxyByMaxz; i-- > 0; ) {
// reset attributes no longer needed
types[i] &= ~NOT_MAXIMUM;
}
Collections.sort(maxpoints, Coordinate::compare);
// Build a map between the original id and the new id following the sort
final int[] idMap = new int[maxpoints.size() + 1];
// Label the points
for (int i = 0; i < maxpoints.size(); i++) {
final int newId = (i + 1);
final Coordinate c = maxpoints.get(i);
final int oldId = c.id;
idMap[oldId] = newId;
c.setId(newId);
}
reassignMaxima(maxima, idMap);
return maxpoints.toArray(new Coordinate[0]);
}
Aggregations