use of uk.ac.sussex.gdsc.smlm.results.procedures.WidthResultProcedure in project GDSC-SMLM by aherbert.
the class CreateData method showSummary.
private double showSummary(List<? extends FluorophoreSequenceModel> fluorophores, List<LocalisationModel> localisations) {
IJ.showStatus("Calculating statistics ...");
final Statistics[] stats = new Statistics[NAMES.length];
for (int i = 0; i < stats.length; i++) {
stats[i] = (settings.getShowHistograms() || alwaysRemoveOutliers[i]) ? new StoredDataStatistics() : new Statistics();
}
// Find the largest timepoint
final ImagePlus outputImp = WindowManager.getImage(benchmarkImageId);
int frameCount;
if (outputImp == null) {
sortLocalisationsByTime(localisations);
frameCount = localisations.get(localisations.size() - 1).getTime();
} else {
frameCount = outputImp.getStackSize();
}
final int[] countHistogram = new int[frameCount + 1];
// Use the localisations that were drawn to create the sampled on/off times
rebuildNeighbours(localisations);
// Assume that there is at least one localisation
final LocalisationModel first = localisations.get(0);
// The current localisation
int currentId = first.getId();
// The last time this localisation was on
int lastT = first.getTime();
// Number of blinks
int blinks = 0;
// On-time of current pulse
int currentT = 0;
double signal = 0;
final double centreOffset = settings.getSize() * 0.5;
// Used to convert the sampled times in frames into seconds
final double framesPerSecond = 1000.0 / settings.getExposureTime();
// final double gain = new CreateDataSettingsHelper(settings).getTotalGainSafe();
for (final LocalisationModel l : localisations) {
final double[] data = l.getData();
if (data == null) {
throw new IllegalStateException("No localisation data. This should not happen!");
}
final double noise = data[1];
final double sx = data[2];
final double sy = data[3];
final double intensityInPhotons = data[4];
// Q. What if the noise is zero, i.e. no background photon / read noise?
// Just ignore it at current. This is only an approximation to the SNR estimate
// if this is not a Gaussian spot.
final double snr = Gaussian2DPeakResultHelper.getMeanSignalUsingP05(intensityInPhotons, sx, sy) / noise;
stats[SIGNAL].add(intensityInPhotons);
stats[NOISE].add(noise);
if (noise != 0) {
stats[SNR].add(snr);
}
// if (l.isContinuous())
if (l.getNext() != null && l.getPrevious() != null) {
stats[SIGNAL_CONTINUOUS].add(intensityInPhotons);
if (noise != 0) {
stats[SNR_CONTINUOUS].add(snr);
}
}
final int id = l.getId();
// Check if this a new fluorophore
if (currentId != id) {
// Add previous fluorophore
stats[SAMPLED_BLINKS].add(blinks);
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
stats[TOTAL_SIGNAL].add(signal);
// Reset
blinks = 0;
currentT = 1;
currentId = id;
signal = intensityInPhotons;
} else {
signal += intensityInPhotons;
// Check if the current fluorophore pulse is broken (i.e. a blink)
if (l.getTime() - 1 > lastT) {
blinks++;
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
currentT = 1;
stats[SAMPLED_T_OFF].add(((l.getTime() - 1) - lastT) / framesPerSecond);
} else {
// Continuous on-time
currentT++;
}
}
lastT = l.getTime();
countHistogram[lastT]++;
stats[X].add((l.getX() - centreOffset) * settings.getPixelPitch());
stats[Y].add((l.getY() - centreOffset) * settings.getPixelPitch());
stats[Z].add(l.getZ() * settings.getPixelPitch());
}
// Final fluorophore
stats[SAMPLED_BLINKS].add(blinks);
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
stats[TOTAL_SIGNAL].add(signal);
// Samples per frame
for (int t = 1; t < countHistogram.length; t++) {
stats[SAMPLES].add(countHistogram[t]);
}
if (fluorophores != null) {
for (final FluorophoreSequenceModel f : fluorophores) {
stats[BLINKS].add(f.getNumberOfBlinks());
// On-time
for (final double t : f.getOnTimes()) {
stats[T_ON].add(t);
}
// Off-time
for (final double t : f.getOffTimes()) {
stats[T_OFF].add(t);
}
}
} else {
// show no blinks
stats[BLINKS].add(0);
stats[T_ON].add(1);
}
if (results != null) {
// Convert depth-of-field to pixels
final double depth = settings.getDepthOfField() / settings.getPixelPitch();
try {
// Get widths
final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
wp.getW();
stats[WIDTH].add(wp.wx);
} catch (final DataException ex) {
ImageJUtils.log("Unable to compute width: " + ex.getMessage());
}
try {
// Get z depth
final StandardResultProcedure sp = new StandardResultProcedure(results, DistanceUnit.PIXEL);
sp.getXyz();
// Get precision
final PrecisionResultProcedure pp = new PrecisionResultProcedure(results);
pp.getPrecision();
stats[PRECISION].add(pp.precisions);
for (int i = 0; i < pp.size(); i++) {
if (Math.abs(sp.z[i]) < depth) {
stats[PRECISION_IN_FOCUS].add(pp.precisions[i]);
}
}
} catch (final DataException ex) {
ImageJUtils.log("Unable to compute LSE precision: " + ex.getMessage());
}
// Compute density per frame. Multi-thread for speed
if (settings.getDensityRadius() > 0) {
final int threadCount = Prefs.getThreads();
final Ticker ticker = ImageJUtils.createTicker(results.getLastFrame(), threadCount, "Calculating density ...");
final ExecutorService threadPool = Executors.newFixedThreadPool(threadCount);
final List<Future<?>> futures = new LinkedList<>();
final TFloatArrayList coordsX = new TFloatArrayList();
final TFloatArrayList coordsY = new TFloatArrayList();
final Statistics densityStats = stats[DENSITY];
final float radius = (float) (settings.getDensityRadius() * getHwhm());
final Rectangle bounds = results.getBounds();
final double area = (double) bounds.width * bounds.height;
// Store the density for each result.
final int[] allDensity = new int[results.size()];
final FrameCounter counter = results.newFrameCounter();
results.forEach((PeakResultProcedure) result -> {
if (counter.advance(result.getFrame())) {
counter.increment(runDensityCalculation(threadPool, futures, coordsX, coordsY, densityStats, radius, area, allDensity, counter.getCount(), ticker));
}
coordsX.add(result.getXPosition());
coordsY.add(result.getYPosition());
});
runDensityCalculation(threadPool, futures, coordsX, coordsY, densityStats, radius, area, allDensity, counter.getCount(), ticker);
ConcurrencyUtils.waitForCompletionUnchecked(futures);
threadPool.shutdown();
ImageJUtils.finished();
// Split results into singles (density = 0) and clustered (density > 0)
final MemoryPeakResults singles = copyMemoryPeakResults("No Density");
final MemoryPeakResults clustered = copyMemoryPeakResults("Density");
counter.reset();
results.forEach((PeakResultProcedure) result -> {
final int density = allDensity[counter.getAndIncrement()];
result.setOrigValue(density);
if (density == 0) {
singles.add(result);
} else {
clustered.add(result);
}
});
}
}
final StringBuilder sb = new StringBuilder();
sb.append(datasetNumber).append('\t');
if (settings.getCameraType() == CameraType.SCMOS) {
sb.append("sCMOS (").append(settings.getCameraModelName()).append(") ");
final Rectangle bounds = cameraModel.getBounds();
sb.append(" ").append(bounds.x).append(",").append(bounds.y);
final int size = settings.getSize();
sb.append(" ").append(size).append("x").append(size);
} else if (CalibrationProtosHelper.isCcdCameraType(settings.getCameraType())) {
sb.append(CalibrationProtosHelper.getName(settings.getCameraType()));
final int size = settings.getSize();
sb.append(" ").append(size).append("x").append(size);
if (settings.getCameraType() == CameraType.EMCCD) {
sb.append(" EM=").append(settings.getEmGain());
}
sb.append(" CG=").append(settings.getCameraGain());
sb.append(" RN=").append(settings.getReadNoise());
sb.append(" B=").append(settings.getBias());
} else {
throw new IllegalStateException();
}
sb.append(" QE=").append(settings.getQuantumEfficiency()).append('\t');
sb.append(settings.getPsfModel());
if (psfModelType == PSF_MODEL_IMAGE) {
sb.append(" Image").append(settings.getPsfImageName());
} else if (psfModelType == PSF_MODEL_ASTIGMATISM) {
sb.append(" model=").append(settings.getAstigmatismModel());
} else {
sb.append(" DoF=").append(MathUtils.rounded(settings.getDepthOfFocus()));
if (settings.getEnterWidth()) {
sb.append(" SD=").append(MathUtils.rounded(settings.getPsfSd()));
} else {
sb.append(" λ=").append(MathUtils.rounded(settings.getWavelength()));
sb.append(" NA=").append(MathUtils.rounded(settings.getNumericalAperture()));
}
}
sb.append('\t');
sb.append((fluorophores == null) ? localisations.size() : fluorophores.size()).append('\t');
sb.append(stats[SAMPLED_BLINKS].getN() + (int) stats[SAMPLED_BLINKS].getSum()).append('\t');
sb.append(localisations.size()).append('\t');
sb.append(frameCount).append('\t');
sb.append(MathUtils.rounded(areaInUm)).append('\t');
sb.append(MathUtils.rounded(localisations.size() / (areaInUm * frameCount), 4)).append('\t');
sb.append(MathUtils.rounded(getHwhm(), 4)).append('\t');
double sd = getPsfSd();
sb.append(MathUtils.rounded(sd, 4)).append('\t');
sd *= settings.getPixelPitch();
final double sa = PsfCalculator.squarePixelAdjustment(sd, settings.getPixelPitch()) / settings.getPixelPitch();
sb.append(MathUtils.rounded(sa, 4)).append('\t');
// Width not valid for the Image PSF.
// Q. Is this true? We can approximate the FHWM for a spot-like image PSF.
final int nStats = (psfModelType == PSF_MODEL_IMAGE) ? stats.length - 1 : stats.length;
for (int i = 0; i < nStats; i++) {
final double centre = (alwaysRemoveOutliers[i]) ? ((StoredDataStatistics) stats[i]).getStatistics().getPercentile(50) : stats[i].getMean();
sb.append(MathUtils.rounded(centre, 4)).append('\t');
}
createSummaryTable().accept(sb.toString());
// Show histograms
if (settings.getShowHistograms() && !java.awt.GraphicsEnvironment.isHeadless()) {
IJ.showStatus("Calculating histograms ...");
final boolean[] chosenHistograms = getChoosenHistograms();
final WindowOrganiser wo = new WindowOrganiser();
final HistogramPlotBuilder builder = new HistogramPlotBuilder(TITLE);
for (int i = 0; i < NAMES.length; i++) {
if (chosenHistograms[i]) {
builder.setData((StoredDataStatistics) stats[i]).setName(NAMES[i]).setIntegerBins(integerDisplay[i]).setRemoveOutliersOption((settings.getRemoveOutliers() || alwaysRemoveOutliers[i]) ? 2 : 0).setNumberOfBins(settings.getHistogramBins()).show(wo);
}
}
wo.tile();
}
IJ.showStatus("");
return stats[SIGNAL].getMean();
}
use of uk.ac.sussex.gdsc.smlm.results.procedures.WidthResultProcedure in project GDSC-SMLM by aherbert.
the class PsfCreator method fitPsf.
/**
* Fit the new PSF image and show a graph of the amplitude/width.
*
* @param psfStack the psf stack
* @param loess the loess
* @param cz the cz
* @param averageRange the average range
* @param fitCom the fit com
* @return The width of the PSF in the z-centre
*/
private double fitPsf(ImageStack psfStack, LoessInterpolator loess, int cz, double averageRange, final double[][] fitCom) {
IJ.showStatus("Fitting final PSF");
// is not appropriate for a normalised PSF.
if (fitConfig.getFitSolver() != FitSolver.LVM_LSE) {
ImageJUtils.log(" " + FitProtosHelper.getName(fitConfig.getFitSolver()) + " is not appropriate for final PSF fitting.");
ImageJUtils.log(" Switching to Least Square Estimation");
fitConfig.setFitSolver(FitSolver.LVM_LSE);
if (settings.getInteractiveMode()) {
// This assumes the LVM does not need the calibration
PeakFit.configureFitSolver(config, null, null, 0);
}
}
// Update the box radius since this is used in the fitSpot method.
boxRadius = psfStack.getWidth() / 2;
final int x = boxRadius;
final int y = boxRadius;
final double shift = fitConfig.getCoordinateShiftFactor();
// Scale the PSF
final PSF.Builder localPsf = fitConfig.getPsf().toBuilder();
for (int i = 0; i < localPsf.getParametersCount(); i++) {
final PSFParameter param = localPsf.getParameters(i);
if (param.getUnit() == PSFParameterUnit.DISTANCE) {
final PSFParameter.Builder b = param.toBuilder();
b.setValue(b.getValue() * settings.getMagnification());
localPsf.setParameters(i, b);
}
}
fitConfig.setPsf(localPsf.build());
// Need to be updated after the widths have been set
fitConfig.setCoordinateShiftFactor(shift);
fitConfig.setBackgroundFitting(false);
// Since the PSF will be normalised remove the camera calibration
fitConfig.setCameraType(CameraType.CAMERA_TYPE_NA);
fitConfig.setMinPhotons(0);
fitConfig.setBias(0);
fitConfig.setGain(1);
// No complex filtering so we get a fit. It should be easy to fit anyway.
fitConfig.setPrecisionThreshold(0);
fitConfig.setDirectFilter(null);
// fitConfig.setDisableSimpleFilter(true);
// Use this for debugging the fit
// fitConfig.setLog(uk.ac.sussex.gdsc.core.ij.ImageJPluginLoggerHelper.getDefaultLogger());
final MemoryPeakResults results = fitSpot(psfStack, psfStack.getWidth(), psfStack.getHeight(), x, y);
if (results.size() < 5) {
ImageJUtils.log(" Final PSF: Not enough fit results %d", results.size());
return 0;
}
// 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 = new double[z.length];
final double[] a = new double[z.length];
// Set limits for the fit
final float maxWidth = (float) (Math.max(fitConfig.getInitialXSd(), fitConfig.getInitialYSd()) * settings.getMagnification() * 4);
// PSF is normalised to 1
final float maxSignal = 2;
final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
wp.getWxWy();
final HeightResultProcedure hp = new HeightResultProcedure(results, IntensityUnit.COUNT);
hp.getH();
final Counter counter = new Counter();
final Counter counterOk = new Counter();
// We have fit the results so they will be in the preferred units
results.forEach((PeakResultProcedure) peak -> {
int index = counter.getAndIncrement();
final float w = Math.max(wp.wx[index], wp.wy[index]);
if (peak.getIntensity() > maxSignal || w > maxWidth) {
return;
}
index = counterOk.getAndIncrement();
z[index] = peak.getFrame();
fitCom[0][peak.getFrame() - 1] = xCoord[index] = peak.getXPosition() - x;
fitCom[1][peak.getFrame() - 1] = yCoord[index] = peak.getYPosition() - y;
sd[index] = w;
a[index] = hp.heights[index];
});
// Truncate
final double[] z2 = Arrays.copyOf(z, counter.getCount());
final double[] xCoord2 = Arrays.copyOf(xCoord, z2.length);
final double[] yCoord2 = Arrays.copyOf(yCoord, z2.length);
final double[] sd2 = Arrays.copyOf(sd, z2.length);
final double[] a2 = Arrays.copyOf(a, z2.length);
// Extract the average smoothed range from the individual fits
final int r = (int) Math.ceil(averageRange / 2);
int start = 0;
int stop = z2.length - 1;
for (int j = 0; j < z2.length; j++) {
if (z2[j] > cz - r) {
start = j;
break;
}
}
for (int j = z2.length; j-- > 0; ) {
if (z2[j] < cz + r) {
stop = j;
break;
}
}
// 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];
double[] smoothA = new double[smoothX.length];
final double[] newZ = new double[smoothX.length];
int smoothCzIndex = 0;
for (int j = start, k = 0; j <= stop; j++, k++) {
smoothX[k] = xCoord2[j];
smoothY[k] = yCoord2[j];
smoothSd[k] = sd2[j];
smoothA[k] = a2[j];
newZ[k] = z2[j];
if (newZ[k] == cz) {
smoothCzIndex = k;
}
}
smoothX = loess.smooth(newZ, smoothX);
smoothY = loess.smooth(newZ, smoothY);
smoothSd = loess.smooth(newZ, smoothSd);
smoothA = loess.smooth(newZ, smoothA);
// Update the widths and positions using the magnification
final double scale = 1.0 / settings.getMagnification();
for (int j = 0; j < xCoord2.length; j++) {
xCoord2[j] *= scale;
yCoord2[j] *= scale;
sd2[j] *= scale;
}
for (int j = 0; j < smoothX.length; j++) {
smoothX[j] *= scale;
smoothY[j] *= scale;
smoothSd[j] *= scale;
}
showPlots(z2, a2, newZ, smoothA, xCoord2, yCoord2, sd2, newZ, smoothX, smoothY, smoothSd, cz);
// Store the data for replotting
this.z = z2;
this.amplitude = a2;
this.smoothAz = newZ;
this.smoothA = smoothA;
this.xCoord = xCoord2;
this.yCoord = yCoord2;
this.sd = sd2;
this.newZ = newZ;
this.smoothX = smoothX;
this.smoothY = smoothY;
this.smoothSd = smoothSd;
// maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
return smoothSd[smoothCzIndex];
}
use of uk.ac.sussex.gdsc.smlm.results.procedures.WidthResultProcedure 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.procedures.WidthResultProcedure in project GDSC-SMLM by aherbert.
the class SpotInspector method run.
@Override
public void run(String arg) {
SmlmUsageTracker.recordPlugin(this.getClass(), arg);
if (MemoryPeakResults.isMemoryEmpty()) {
IJ.error(TITLE, "No localisations in memory");
return;
}
if (!showDialog()) {
return;
}
// Load the results
results = ResultsManager.loadInputResults(settings.inputOption, false, DistanceUnit.PIXEL, null);
if (MemoryPeakResults.isEmpty(results)) {
IJ.error(TITLE, "No results could be loaded");
IJ.showStatus("");
return;
}
// Check if the original image is open
ImageSource source = results.getSource();
if (source == null) {
IJ.error(TITLE, "Unknown original source image");
return;
}
source = source.getOriginal();
source.setReadHint(ReadHint.NONSEQUENTIAL);
if (!source.open()) {
IJ.error(TITLE, "Cannot open original source image: " + source.toString());
return;
}
final float stdDevMax = getStandardDeviation(results);
if (stdDevMax < 0) {
// TODO - Add dialog to get the initial peak width
IJ.error(TITLE, "Fitting configuration (for initial peak width) is not available");
return;
}
// Rank spots
rankedResults = new ArrayList<>(results.size());
// Data for the sorting
final PrecisionResultProcedure pp;
if (settings.sortOrderIndex == 1) {
pp = new PrecisionResultProcedure(results);
pp.getPrecision();
} else {
pp = null;
}
// Build procedures to get:
// Shift = position in pixels - originXY
final StandardResultProcedure sp;
if (settings.sortOrderIndex == 9) {
sp = new StandardResultProcedure(results, DistanceUnit.PIXEL);
sp.getXyr();
} else {
sp = null;
}
// SD = gaussian widths only for Gaussian PSFs
final WidthResultProcedure wp;
if (settings.sortOrderIndex >= 6 && settings.sortOrderIndex <= 8) {
wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
wp.getWxWy();
} else {
wp = null;
}
// Amplitude for Gaussian PSFs
final HeightResultProcedure hp;
if (settings.sortOrderIndex == 2) {
hp = new HeightResultProcedure(results, IntensityUnit.PHOTON);
hp.getH();
} else {
hp = null;
}
final Counter c = new Counter();
results.forEach((PeakResultProcedure) result -> {
final float[] score = getScore(result, c.getAndIncrement(), pp, sp, wp, hp, stdDevMax);
rankedResults.add(new PeakResultRank(result, score[0], score[1]));
});
Collections.sort(rankedResults, PeakResultRank::compare);
// Prepare results table
final ImageJTablePeakResults table = new ImageJTablePeakResults(false, results.getName(), true);
table.copySettings(results);
table.setTableTitle(TITLE);
table.setAddCounter(true);
table.setShowZ(results.is3D());
// TODO - Add to settings
table.setShowFittingData(true);
table.setShowNoiseData(true);
if (settings.showCalibratedValues) {
table.setDistanceUnit(DistanceUnit.NM);
table.setIntensityUnit(IntensityUnit.PHOTON);
} else {
table.setDistanceUnit(DistanceUnit.PIXEL);
table.setIntensityUnit(IntensityUnit.COUNT);
}
table.begin();
// Add a mouse listener to jump to the frame for the clicked line
textPanel = table.getResultsWindow().getTextPanel();
// We must ignore old instances of this class from the mouse listeners
id = currentId.incrementAndGet();
textPanel.addMouseListener(new MouseAdapter() {
@Override
public void mouseClicked(MouseEvent event) {
SpotInspector.this.mouseClicked(event);
}
});
// Add results to the table
int count = 0;
for (final PeakResultRank rank : rankedResults) {
rank.rank = count++;
table.add(rank.peakResult);
}
table.end();
if (settings.plotScore || settings.plotHistogram) {
// Get values for the plots
float[] xValues = null;
float[] yValues = null;
double yMin;
double yMax;
int spotNumber = 0;
xValues = new float[rankedResults.size()];
yValues = new float[xValues.length];
for (final PeakResultRank rank : rankedResults) {
xValues[spotNumber] = spotNumber + 1;
yValues[spotNumber++] = recoverScore(rank.score);
}
// Set the min and max y-values using 1.5 x IQR
final DescriptiveStatistics stats = new DescriptiveStatistics();
for (final float v : yValues) {
stats.addValue(v);
}
if (settings.removeOutliers) {
final double lower = stats.getPercentile(25);
final double upper = stats.getPercentile(75);
final double iqr = upper - lower;
yMin = Math.max(lower - iqr, stats.getMin());
yMax = Math.min(upper + iqr, stats.getMax());
IJ.log(String.format("Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax));
} else {
yMin = stats.getMin();
yMax = stats.getMax();
IJ.log(String.format("Data range: %f - %f", yMin, yMax));
}
plotScore(xValues, yValues, yMin, yMax);
plotHistogram(yValues);
}
// Extract spots into a stack
final int w = source.getWidth();
final int h = source.getHeight();
final int size = 2 * settings.radius + 1;
final ImageStack spots = new ImageStack(size, size, rankedResults.size());
// To assist the extraction of data from the image source, process them in time order to allow
// frame caching. Then set the appropriate slice in the result stack
Collections.sort(rankedResults, (o1, o2) -> Integer.compare(o1.peakResult.getFrame(), o2.peakResult.getFrame()));
for (final PeakResultRank rank : rankedResults) {
final PeakResult r = rank.peakResult;
// Extract image
// Note that the coordinates are relative to the middle of the pixel (0.5 offset)
// so do not round but simply convert to int
final int x = (int) (r.getXPosition());
final int y = (int) (r.getYPosition());
// Extract a region but crop to the image bounds
int minX = x - settings.radius;
int minY = y - settings.radius;
final int maxX = Math.min(x + settings.radius + 1, w);
final int maxY = Math.min(y + settings.radius + 1, h);
int padX = 0;
int padY = 0;
if (minX < 0) {
padX = -minX;
minX = 0;
}
if (minY < 0) {
padY = -minY;
minY = 0;
}
final int sizeX = maxX - minX;
final int sizeY = maxY - minY;
float[] data = source.get(r.getFrame(), new Rectangle(minX, minY, sizeX, sizeY));
// Prevent errors with missing data
if (data == null) {
data = new float[sizeX * sizeY];
}
ImageProcessor spotIp = new FloatProcessor(sizeX, sizeY, data, null);
// Pad if necessary, i.e. the crop is too small for the stack
if (padX > 0 || padY > 0 || sizeX < size || sizeY < size) {
final ImageProcessor spotIp2 = spotIp.createProcessor(size, size);
spotIp2.insert(spotIp, padX, padY);
spotIp = spotIp2;
}
final int slice = rank.rank + 1;
spots.setPixels(spotIp.getPixels(), slice);
spots.setSliceLabel(MathUtils.rounded(rank.originalScore), slice);
}
source.close();
// Reset to the rank order
Collections.sort(rankedResults, PeakResultRank::compare);
final ImagePlus imp = ImageJUtils.display(TITLE, spots);
imp.setRoi((PointRoi) null);
// Make bigger
for (int i = 10; i-- > 0; ) {
imp.getWindow().getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
}
}
use of uk.ac.sussex.gdsc.smlm.results.procedures.WidthResultProcedure in project GDSC-SMLM by aherbert.
the class AstigmatismModelManager method plotData.
private boolean plotData() {
if (results.size() <= imp.getStackSize() / 2) {
IJ.error(TITLE, "Not enough fit results " + results.size());
return false;
}
final double umPerSlice = pluginSettings.getNmPerSlice() / 1000.0;
// final double nmPerPixel = results.getNmPerPixel();
z = new double[results.size()];
x = new double[z.length];
y = new double[z.length];
intensity = new double[z.length];
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() * umPerSlice;
x[i] = (peak.getXPosition() - cx);
y[i] = (peak.getYPosition() - cy);
intensity[i] = peak.getIntensity();
}
});
final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
wp.getWxWy();
sx = SimpleArrayUtils.toDouble(wp.wx);
sy = SimpleArrayUtils.toDouble(wp.wy);
final WindowOrganiser wo = new WindowOrganiser();
plot(wo, z, "Intensity (photon)", intensity, "Intensity", null, null);
xyPlot = plot(wo, z, "Position (px)", x, "X", y, "Y");
widthPlot = plot(wo, z, "Width (px)", sx, "Sx", sy, "Sy");
wo.tile();
return true;
}
Aggregations