use of uk.ac.sussex.gdsc.smlm.results.count.Counter in project GDSC-SMLM by aherbert.
the class SpotInspector method mouseClicked.
private void mouseClicked(MouseEvent event) {
if (id != currentId.get()) {
return;
}
// Show the result that was double clicked in the result table
if (event.getClickCount() > 1) {
final int rank = textPanel.getSelectionStart() + 1;
// Show the spot that was double clicked
final ImagePlus imp = WindowManager.getImage(TITLE);
if (imp != null && rank > 0 && rank <= imp.getStackSize()) {
imp.setSlice(rank);
if (imp.getWindow() != null) {
imp.getWindow().toFront();
}
final PeakResult r = rankedResults.get(rank - 1).peakResult;
final TypeConverter<DistanceUnit> dc = results.getDistanceConverter(DistanceUnit.PIXEL);
final float rx = dc.convert(r.getXPosition());
final float ry = dc.convert(r.getYPosition());
final int x = (int) rx;
final int y = (int) ry;
// Find bounds
final int minX = x - settings.radius;
final int minY = y - settings.radius;
// Require the Shift key to add all spots
if (!event.isShiftDown()) {
// Add the single clicked spot
imp.setRoi(new OffsetPointRoi(rx - minX, ry - minY));
return;
}
// Add all the spots
final int maxX = x + settings.radius + 1;
final int maxY = y + settings.radius + 1;
// Create ROIs
final HashSet<Point2D.Float> spots = new HashSet<>();
results.forEach(DistanceUnit.PIXEL, (XyResultProcedure) (xp, yp) -> {
if (xp > minX && xp < maxX && yp > minY && yp < maxY) {
// Use only unique points
spots.add(new Point2D.Float(xp - minX, yp - minY));
}
});
final int points = spots.size();
final float[] ox = new float[points];
final float[] oy = new float[points];
final Counter c = new Counter();
spots.forEach(p -> {
ox[c.getCount()] = p.x;
oy[c.getAndIncrement()] = p.y;
});
imp.setRoi(new OffsetPointRoi(ox, oy, points));
}
}
}
use of uk.ac.sussex.gdsc.smlm.results.count.Counter in project GDSC-SMLM by aherbert.
the class PsfCreator method runUsingFitting.
private void runUsingFitting() {
if (!showFittingDialog()) {
return;
}
if (!loadConfiguration()) {
return;
}
final BasePoint[] spots = getSpots(0, true);
if (spots.length == 0) {
IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
return;
}
final ImageStack stack = getImageStack();
final int width = imp.getWidth();
final int height = imp.getHeight();
final int currentSlice = imp.getSlice();
// Adjust settings for a single maxima
config.setIncludeNeighbours(false);
final ArrayList<double[]> centres = new ArrayList<>(spots.length);
final int iterations = 1;
final LoessInterpolator loess = new LoessInterpolator(settings.getSmoothing(), iterations);
// TODO - The fitting routine may not produce many points. In this instance the LOESS
// interpolator
// fails to smooth the data very well. A higher bandwidth helps this but perhaps
// try a different smoothing method.
// For each spot
ImageJUtils.log(TITLE + ": " + imp.getTitle());
ImageJUtils.log("Finding spot locations...");
ImageJUtils.log(" %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
final StoredDataStatistics averageSd = new StoredDataStatistics();
final StoredDataStatistics averageA = new StoredDataStatistics();
final Statistics averageRange = new Statistics();
final MemoryPeakResults allResults = new MemoryPeakResults();
allResults.setCalibration(fitConfig.getCalibration());
allResults.setPsf(fitConfig.getPsf());
allResults.setName(TITLE);
allResults.setBounds(new Rectangle(0, 0, width, height));
MemoryPeakResults.addResults(allResults);
for (int n = 1; n <= spots.length; n++) {
final BasePoint spot = spots[n - 1];
final int x = (int) spot.getX();
final int y = (int) spot.getY();
final MemoryPeakResults results = fitSpot(stack, width, height, x, y);
allResults.add(results);
if (results.size() < 5) {
ImageJUtils.log(" Spot %d: Not enough fit results %d", n, results.size());
continue;
}
// Get the results for the spot centre and width
final double[] z = new double[results.size()];
final double[] xCoord = new double[z.length];
final double[] yCoord = new double[z.length];
final double[] sd;
final double[] a;
final Counter counter = new Counter();
// We have fit the results so they will be in the preferred units
results.forEach(new PeakResultProcedure() {
@Override
public void execute(PeakResult peak) {
final int i = counter.getAndIncrement();
z[i] = peak.getFrame();
xCoord[i] = peak.getXPosition() - x;
yCoord[i] = peak.getYPosition() - y;
}
});
final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
wp.getW();
sd = SimpleArrayUtils.toDouble(wp.wx);
final HeightResultProcedure hp = new HeightResultProcedure(results, IntensityUnit.COUNT);
hp.getH();
a = SimpleArrayUtils.toDouble(hp.heights);
// Smooth the amplitude plot
final double[] smoothA = loess.smooth(z, a);
// Find the maximum amplitude
int maximumIndex = findMaximumIndex(smoothA);
// Find the range at a fraction of the max. This is smoothed to find the X/Y centre
int start = 0;
int stop = smoothA.length - 1;
final double limit = smoothA[maximumIndex] * settings.getAmplitudeFraction();
for (int j = 0; j < smoothA.length; j++) {
if (smoothA[j] > limit) {
start = j;
break;
}
}
for (int j = smoothA.length; j-- > 0; ) {
if (smoothA[j] > limit) {
stop = j;
break;
}
}
averageRange.add(stop - start + 1);
// Extract xy centre coords and smooth
double[] smoothX = new double[stop - start + 1];
double[] smoothY = new double[smoothX.length];
double[] smoothSd = new double[smoothX.length];
final double[] newZ = new double[smoothX.length];
for (int j = start, k = 0; j <= stop; j++, k++) {
smoothX[k] = xCoord[j];
smoothY[k] = yCoord[j];
smoothSd[k] = sd[j];
newZ[k] = z[j];
}
smoothX = loess.smooth(newZ, smoothX);
smoothY = loess.smooth(newZ, smoothY);
smoothSd = loess.smooth(newZ, smoothSd);
// Since the amplitude is not very consistent move from this peak to the
// lowest width which is the in-focus spot.
maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
// Find the centre at the amplitude peak
final double cx = smoothX[maximumIndex] + x;
final double cy = smoothY[maximumIndex] + y;
int cz = (int) newZ[maximumIndex];
double csd = smoothSd[maximumIndex];
double ca = smoothA[maximumIndex + start];
// The average should weight the SD using the signal for each spot
averageSd.add(smoothSd[maximumIndex]);
averageA.add(ca);
if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) {
ImageJUtils.log(" Spot %d was ignored", n);
continue;
}
// Store result - it may have been moved interactively
maximumIndex += this.slice - cz;
cz = (int) newZ[maximumIndex];
csd = smoothSd[maximumIndex];
ca = smoothA[maximumIndex + start];
ImageJUtils.log(" Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f", n, cx, cy, cz, csd, ca);
centres.add(new double[] { cx, cy, cz, csd, n });
}
if (settings.getInteractiveMode()) {
imp.setSlice(currentSlice);
imp.setOverlay(null);
// Hide the amplitude and spot plots
ImageJUtils.hide(TITLE_AMPLITUDE);
ImageJUtils.hide(TITLE_PSF_PARAMETERS);
}
if (centres.isEmpty()) {
final String msg = "No suitable spots could be identified";
ImageJUtils.log(msg);
IJ.error(TITLE, msg);
return;
}
// Find the limits of the z-centre
int minz = (int) centres.get(0)[2];
int maxz = minz;
for (final double[] centre : centres) {
if (minz > centre[2]) {
minz = (int) centre[2];
} else if (maxz < centre[2]) {
maxz = (int) centre[2];
}
}
IJ.showStatus("Creating PSF image");
// Create a stack that can hold all the data.
final ImageStack psf = createStack(stack, minz, maxz, settings.getMagnification());
// For each spot
final Statistics stats = new Statistics();
boolean ok = true;
for (int i = 0; ok && i < centres.size(); i++) {
final double increment = 1.0 / (stack.getSize() * centres.size());
setProgress((double) i / centres.size());
final double[] centre = centres.get(i);
// Extract the spot
final float[][] spot = new float[stack.getSize()][];
Rectangle regionBounds = null;
for (int slice = 1; slice <= stack.getSize(); slice++) {
final ImageExtractor ie = ImageExtractor.wrap((float[]) stack.getPixels(slice), width, height);
if (regionBounds == null) {
regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius);
}
spot[slice - 1] = ie.crop(regionBounds);
}
if (regionBounds == null) {
// Empty stack
continue;
}
final int n = (int) centre[4];
final float b = getBackground(n, spot);
if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
ImageJUtils.log(" Spot %d was ignored", n);
continue;
}
stats.add(b);
// Adjust the centre using the crop
centre[0] -= regionBounds.x;
centre[1] -= regionBounds.y;
// This takes a long time so this should track progress
ok = addToPsf(maxz, settings.getMagnification(), psf, centre, spot, regionBounds, increment, settings.getCentreEachSlice());
}
if (settings.getInteractiveMode()) {
ImageJUtils.hide(TITLE_INTENSITY);
}
IJ.showProgress(1);
if (!ok || stats.getN() == 0) {
return;
}
final double avSd = getAverage(averageSd, averageA, 2);
ImageJUtils.log(" Average background = %.2f, Av. SD = %s px", stats.getMean(), MathUtils.rounded(avSd, 4));
normalise(psf, maxz, avSd * settings.getMagnification(), false);
IJ.showProgress(1);
psfImp = ImageJUtils.display(TITLE_PSF, psf);
psfImp.setSlice(maxz);
psfImp.resetDisplayRange();
psfImp.updateAndDraw();
final double[][] fitCom = new double[2][psf.getSize()];
Arrays.fill(fitCom[0], Double.NaN);
Arrays.fill(fitCom[1], Double.NaN);
final double fittedSd = fitPsf(psf, loess, maxz, averageRange.getMean(), fitCom);
// Compute the drift in the PSF:
// - Use fitted centre if available; otherwise find CoM for each frame
// - express relative to the average centre
final double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / settings.getMagnification());
final double[] slice = SimpleArrayUtils.newArray(psf.getSize(), 1, 1.0);
final String title = TITLE + " CoM Drift";
final Plot plot = new Plot(title, "Slice", "Drift (nm)");
plot.addLabel(0, 0, "Red = X; Blue = Y");
// double[] limitsX = Maths.limits(com[0]);
// double[] limitsY = Maths.limits(com[1]);
final double[] limitsX = getLimits(com[0]);
final double[] limitsY = getLimits(com[1]);
plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1]));
plot.setColor(Color.red);
plot.addPoints(slice, com[0], Plot.DOT);
plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(slice, com[1], Plot.DOT);
plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE);
ImageJUtils.display(title, plot);
// TODO - Redraw the PSF with drift correction applied.
// This means that the final image should have no drift.
// This is relevant when combining PSF images. It doesn't matter too much for simulations
// unless the drift is large.
// Add Image properties containing the PSF details
final double fwhm = getFwhm(psf, maxz);
psfImp.setProperty("Info", ImagePsfHelper.toString(ImagePsfHelper.create(maxz, nmPerPixel / settings.getMagnification(), settings.getNmPerSlice(), stats.getN(), fwhm, createNote())));
ImageJUtils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, " + "PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, MathUtils.rounded(nmPerPixel / settings.getMagnification(), 3), MathUtils.rounded(settings.getNmPerSlice(), 3), stats.getN(), MathUtils.rounded(fittedSd * nmPerPixel, 4), MathUtils.rounded(fwhm));
createInteractivePlots(psf, maxz, nmPerPixel / settings.getMagnification(), fittedSd * nmPerPixel);
IJ.showStatus("");
}
use of uk.ac.sussex.gdsc.smlm.results.count.Counter in project GDSC-SMLM by aherbert.
the class PeakFit method getSingleFrame.
/**
* Gets the single frame containing all the results (if they are all in a single frame), else 0.
*
* @param results the results (must not be empty)
* @return the single frame (or zero)
*/
private static int getSingleFrame(MemoryPeakResults results) {
final FrameCounter counter = new FrameCounter(results.getFirstFrame());
// The counter will return true (stop execution) if a new frame
results.forEach((PeakResultProcedureX) peakResult -> counter.advance(peakResult.getFrame()));
if (counter.currentFrame() != counter.previousFrame()) {
return 0;
}
return counter.currentFrame();
}
use of uk.ac.sussex.gdsc.smlm.results.count.Counter in project GDSC-SMLM by aherbert.
the class Fire method createImages.
/**
* Creates the images to use for the FIRE calculation. This must be called after
* {@link #initialise(MemoryPeakResults, MemoryPeakResults)}.
*
* @param fourierImageScale the fourier image scale (set to zero to auto compute)
* @param imageSize the image size
* @param useSignal Use the localisation signal to weight the intensity. The default uses a value
* of 1 per localisation.
* @return the fire images
*/
public FireImages createImages(double fourierImageScale, int imageSize, boolean useSignal) {
if (results == null) {
return null;
}
final SignalProvider signalProvider = (useSignal && (results.hasIntensity())) ? new PeakSignalProvider() : new FixedSignalProvider();
// Draw images using the existing IJ routines.
final Rectangle bounds = new Rectangle((int) Math.ceil(dataBounds.getWidth()), (int) Math.ceil(dataBounds.getHeight()));
final ResultsImageSettings.Builder builder = ResultsImageSettings.newBuilder().setImageType(ResultsImageType.DRAW_NONE).setWeighted(true).setEqualised(false).setImageMode(ResultsImageMode.IMAGE_ADD);
if (fourierImageScale > 0) {
builder.setImageSizeMode(ResultsImageSizeMode.SCALED);
builder.setScale(fourierImageScale);
} else {
builder.setImageSizeMode(ResultsImageSizeMode.IMAGE_SIZE);
builder.setImageSize(imageSize);
}
ImageJImagePeakResults image1 = createPeakResultsImage(bounds, builder, "IP1");
ImageJImagePeakResults image2 = createPeakResultsImage(bounds, builder, "IP2");
final float minx = (float) dataBounds.getX();
final float miny = (float) dataBounds.getY();
if (this.results2 != null) {
// Two image comparison
final ImageJImagePeakResults i1 = image1;
results.forEach((PeakResultProcedure) result -> {
final float x = result.getXPosition() - minx;
final float y = result.getYPosition() - miny;
i1.add(x, y, signalProvider.getSignal(result));
});
final ImageJImagePeakResults i2 = image2;
results2.forEach((PeakResultProcedure) result -> {
final float x = result.getXPosition() - minx;
final float y = result.getYPosition() - miny;
i2.add(x, y, signalProvider.getSignal(result));
});
} else {
// Block sampling.
// Ensure we have at least 2 even sized blocks.
int blockSize = Math.min(results.size() / 2, Math.max(1, settings.blockSize));
int nblocks = (int) Math.ceil((double) results.size() / blockSize);
while (nblocks <= 1 && blockSize > 1) {
blockSize /= 2;
nblocks = (int) Math.ceil((double) results.size() / blockSize);
}
if (nblocks <= 1) {
// This should not happen since the results should contain at least 2 localisations
return null;
}
if (blockSize != settings.blockSize) {
IJ.log(pluginTitle + " Warning: Changed block size to " + blockSize);
}
final Counter i = new Counter();
final Counter block = new Counter();
final int finalBlockSize = blockSize;
final PeakResult[][] blocks = new PeakResult[nblocks][blockSize];
results.forEach((PeakResultProcedure) result -> {
if (i.getCount() == finalBlockSize) {
block.increment();
i.reset();
}
blocks[block.getCount()][i.getAndIncrement()] = result;
});
// Truncate last block
blocks[block.getCount()] = Arrays.copyOf(blocks[block.getCount()], i.getCount());
final int[] indices = SimpleArrayUtils.natural(block.getCount() + 1);
if (settings.randomSplit) {
MathArrays.shuffle(indices);
}
for (final int index : indices) {
// Split alternating so just rotate
final ImageJImagePeakResults image = image1;
image1 = image2;
image2 = image;
for (final PeakResult p : blocks[index]) {
final float x = p.getXPosition() - minx;
final float y = p.getYPosition() - miny;
image.add(x, y, signalProvider.getSignal(p));
}
}
}
image1.end();
final ImageProcessor ip1 = image1.getImagePlus().getProcessor();
image2.end();
final ImageProcessor ip2 = image2.getImagePlus().getProcessor();
if (settings.maxPerBin > 0 && signalProvider instanceof FixedSignalProvider) {
// We can eliminate over-sampled pixels
for (int i = ip1.getPixelCount(); i-- > 0; ) {
if (ip1.getf(i) > settings.maxPerBin) {
ip1.setf(i, settings.maxPerBin);
}
if (ip2.getf(i) > settings.maxPerBin) {
ip2.setf(i, settings.maxPerBin);
}
}
}
return new FireImages(ip1, ip2, nmPerUnit / image1.getScale());
}
use of uk.ac.sussex.gdsc.smlm.results.count.Counter in project GDSC-SMLM by aherbert.
the class BenchmarkFilterAnalysis method getCoordinates.
private static TIntObjectHashMap<UniqueIdPeakResult[]> getCoordinates(MemoryPeakResults results) {
final TIntObjectHashMap<UniqueIdPeakResult[]> coords = new TIntObjectHashMap<>();
if (results.size() > 0) {
// Do not use HashMap directly to build the coords object since there
// will be many calls to getEntry(). Instead sort the results and use
// a new list for each time point
results.sort();
final Counter uniqueId = new Counter();
final FrameCounter counter = new FrameCounter();
final LocalList<PeakResult> tmp = new LocalList<>();
// Add the results to the lists
results.forEach((PeakResultProcedure) result -> {
if (counter.advanceAndReset(result.getFrame()) && !tmp.isEmpty()) {
coords.put(counter.previousFrame(), tmp.toArray(new UniqueIdPeakResult[0]));
tmp.clear();
}
tmp.add(new UniqueIdPeakResult(tmp.size(), uniqueId.getAndIncrement(), result));
});
if (!tmp.isEmpty()) {
coords.put(counter.currentFrame(), tmp.toArray(new UniqueIdPeakResult[0]));
}
}
return coords;
}
Aggregations