use of uk.ac.sussex.gdsc.smlm.model.ImagePsfModel in project GDSC-SMLM by aherbert.
the class CreateData method reportAndSaveFittingLimits.
/**
* Output the theoretical limits for fitting a Gaussian and store the benchmark settings.
*
* @param dist The distribution
*/
private void reportAndSaveFittingLimits(SpatialDistribution dist) {
ImageJUtils.log(TITLE + " Benchmark");
final double a = settings.getPixelPitch();
final double[] xyz = dist.next().clone();
final int size = settings.getSize();
final double offset = size * 0.5;
for (int i = 0; i < 2; i++) {
xyz[i] += offset;
}
// Get the width for the z-depth by using the PSF Model
final PsfModel psf = createPsfModel(xyz);
psfModelCache = psf;
double sd0;
double sd1;
if (psf instanceof GaussianPsfModel) {
sd0 = ((GaussianPsfModel) psf).getS0(xyz[2]);
sd1 = ((GaussianPsfModel) psf).getS1(xyz[2]);
} else if (psf instanceof AiryPsfModel) {
psf.create3D((double[]) null, size, size, 1, xyz[0], xyz[1], xyz[2], null);
sd0 = ((AiryPsfModel) psf).getW0() * AiryPattern.FACTOR;
sd1 = ((AiryPsfModel) psf).getW1() * AiryPattern.FACTOR;
} else if (psf instanceof ImagePsfModel) {
psf.create3D((double[]) null, size, size, 1, xyz[0], xyz[1], xyz[2], null);
sd0 = ((ImagePsfModel) psf).getHwhm0() / Gaussian2DFunction.SD_TO_HWHM_FACTOR;
sd1 = ((ImagePsfModel) psf).getHwhm1() / Gaussian2DFunction.SD_TO_HWHM_FACTOR;
} else {
throw new IllegalStateException("Unknown PSF: " + psf.getClass().getSimpleName());
}
final double sd = Gaussian2DPeakResultHelper.getStandardDeviation(sd0, sd1) * a;
ImageJUtils.log("X = %s nm : %s px", MathUtils.rounded(xyz[0] * a), MathUtils.rounded(xyz[0], 6));
ImageJUtils.log("Y = %s nm : %s px", MathUtils.rounded(xyz[1] * a), MathUtils.rounded(xyz[1], 6));
ImageJUtils.log("Z = %s nm : %s px", MathUtils.rounded(xyz[2] * a), MathUtils.rounded(xyz[2], 6));
ImageJUtils.log("Width (s) = %s nm : %s px", MathUtils.rounded(sd), MathUtils.rounded(sd / a));
final double sa = PsfCalculator.squarePixelAdjustment(sd, a);
ImageJUtils.log("Adjusted Width (sa) = %s nm : %s px", MathUtils.rounded(sa), MathUtils.rounded(sa / a));
ImageJUtils.log("Signal (N) = %s - %s photons", MathUtils.rounded(settings.getPhotonsPerSecond()), MathUtils.rounded(settings.getPhotonsPerSecondMaximum()));
boolean emCcd;
double totalGain;
final double qe = getQuantumEfficiency();
final double noise = getNoiseEstimateInPhotoelectrons(qe);
double readNoise;
if (CalibrationProtosHelper.isCcdCameraType(settings.getCameraType())) {
final CreateDataSettingsHelper helper = new CreateDataSettingsHelper(settings);
emCcd = helper.isEmCcd;
totalGain = helper.getTotalGainSafe();
// Store read noise in ADUs
readNoise = settings.getReadNoise() * ((settings.getCameraGain() > 0) ? settings.getCameraGain() : 1);
} else if (settings.getCameraType() == CameraType.SCMOS) {
// Assume sCMOS amplification is like a CCD for the precision computation.
emCcd = false;
// Not required for sCMOS
totalGain = 0;
readNoise = 0;
} else {
throw new IllegalStateException("Unknown camera type: " + settings.getCameraType());
}
// The precision calculation is dependent on the model. The classic Mortensen formula
// is for a Gaussian Mask Estimator. Use other equation for MLE. The formula provided
// for WLSE requires an offset to the background used to stabilise the fitting. This is
// not implemented (i.e. we used an offset of zero) and in this case the WLSE precision
// is the same as MLE with the caveat of numerical instability.
// Apply QE directly to simulated photons to allow computation of bounds
// with the captured photons...
final double min = settings.getPhotonsPerSecond() * qe;
final double max = settings.getPhotonsPerSecondMaximum() * qe;
final double lowerP = Gaussian2DPeakResultHelper.getPrecision(a, sd, max, noise, emCcd);
final double upperP = Gaussian2DPeakResultHelper.getPrecision(a, sd, min, noise, emCcd);
final double lowerMlP = Gaussian2DPeakResultHelper.getMLPrecision(a, sd, max, noise, emCcd);
final double upperMlP = Gaussian2DPeakResultHelper.getMLPrecision(a, sd, min, noise, emCcd);
final double lowerN = getPrecisionN(a, sd, min, MathUtils.pow2(noise), emCcd);
final double upperN = getPrecisionN(a, sd, max, MathUtils.pow2(noise), emCcd);
if (settings.getCameraType() == CameraType.SCMOS) {
ImageJUtils.log("sCMOS camera background estimate uses an average read noise");
}
ImageJUtils.log("Effective background noise = %s photo-electron " + "[includes read noise and background photons]", MathUtils.rounded(noise));
ImageJUtils.log("Localisation precision (LSE): %s - %s nm : %s - %s px", MathUtils.rounded(lowerP), MathUtils.rounded(upperP), MathUtils.rounded(lowerP / a), MathUtils.rounded(upperP / a));
ImageJUtils.log("Localisation precision (MLE): %s - %s nm : %s - %s px", MathUtils.rounded(lowerMlP), MathUtils.rounded(upperMlP), MathUtils.rounded(lowerMlP / a), MathUtils.rounded(upperMlP / a));
ImageJUtils.log("Signal precision: %s - %s photo-electrons", MathUtils.rounded(lowerN), MathUtils.rounded(upperN));
// Wrap to a function
final PsfModelGradient1Function f = new PsfModelGradient1Function(psf, size, size);
// Set parameters
final double[] params = new double[5];
// No background when computing the SNR
// params[0] = settings.getBackground() * qe;
params[1] = min;
System.arraycopy(xyz, 0, params, 2, 3);
// Compute SNR using mean signal at 50%. Assume the region covers the entire PSF.
final double[] v = new StandardValueProcedure().getValues(f, params);
final double u = FunctionHelper.getMeanValue(v, 0.5);
final double u0 = MathUtils.max(v);
// Store the benchmark settings when not using variable photons
if (Double.compare(min, max) == 0) {
ImageJUtils.log("50%% PSF SNR : %s : Peak SNR : %s", MathUtils.rounded(u / noise), MathUtils.rounded(u0 / noise));
// Compute the true CRLB using the fisher information
createLikelihoodFunction();
// Compute Fisher information
final UnivariateLikelihoodFisherInformationCalculator c = new UnivariateLikelihoodFisherInformationCalculator(f, fiFunction);
// Get limits: Include background in the params
params[0] = settings.getBackground() * qe;
final FisherInformationMatrix m = c.compute(params);
// Report and store the limits
final double[] crlb = m.crlbSqrt();
if (crlb != null) {
ImageJUtils.log("Localisation precision (CRLB): B=%s, I=%s photons", MathUtils.rounded(crlb[0]), MathUtils.rounded(crlb[1]));
ImageJUtils.log("Localisation precision (CRLB): X=%s, Y=%s, Z=%s nm : %s,%s,%s px", MathUtils.rounded(crlb[2] * a), MathUtils.rounded(crlb[3] * a), MathUtils.rounded(crlb[4] * a), MathUtils.rounded(crlb[2]), MathUtils.rounded(crlb[3]), MathUtils.rounded(crlb[4]));
}
benchmarkParameters = new BenchmarkParameters(settings.getParticles(), sd, a, settings.getPhotonsPerSecond(), xyz[0], xyz[1], xyz[2], settings.getBias(), totalGain, qe, readNoise, settings.getCameraType(), settings.getCameraModelName(), cameraModel.getBounds(), settings.getBackground(), noise, lowerN, lowerP, lowerMlP, createPsf(sd / a), crlb);
} else {
// SNR will just scale
final double scale = max / min;
ImageJUtils.log("50%% PSF SNR : %s - %s : Peak SNR : %s - %s", MathUtils.rounded(u / noise), MathUtils.rounded(scale * u / noise), MathUtils.rounded(u0 / noise), MathUtils.rounded(scale * u0 / noise));
ImageJUtils.log("Warning: Benchmark settings are only stored in memory when the number of photons is " + "fixed. Min %s != Max %s", MathUtils.rounded(settings.getPhotonsPerSecond()), MathUtils.rounded(settings.getPhotonsPerSecondMaximum()));
}
}
use of uk.ac.sussex.gdsc.smlm.model.ImagePsfModel in project GDSC-SMLM by aherbert.
the class CreateData method createImagePsf.
/**
* Create a PSF model from the image that contains all the z-slices needed to draw the given
* localisations.
*
* @param localisationSets the localisation sets
* @return the image PSF model
*/
private ImagePsfModel createImagePsf(List<LocalisationModelSet> localisationSets) {
final ImagePlus imp = WindowManager.getImage(settings.getPsfImageName());
if (imp == null) {
IJ.error(TITLE, "Unable to create the PSF model from image: " + settings.getPsfImageName());
return null;
}
try {
final ImagePSF psfSettings = ImagePsfHelper.fromString(imp.getProperty("Info").toString());
if (psfSettings == null) {
throw new IllegalStateException("Unknown PSF settings for image: " + imp.getTitle());
}
// Check all the settings have values
if (psfSettings.getPixelSize() <= 0) {
throw new IllegalStateException("Missing nmPerPixel calibration settings for image: " + imp.getTitle());
}
if (psfSettings.getPixelDepth() <= 0) {
throw new IllegalStateException("Missing nmPerSlice calibration settings for image: " + imp.getTitle());
}
if (psfSettings.getCentreImage() <= 0) {
throw new IllegalStateException("Missing zCentre calibration settings for image: " + imp.getTitle());
}
if (psfSettings.getFwhm() <= 0) {
throw new IllegalStateException("Missing FWHM calibration settings for image: " + imp.getTitle());
}
// To save memory construct the Image PSF using only the slices that are within
// the depth of field of the simulation
double minZ = Double.POSITIVE_INFINITY;
double maxZ = Double.NEGATIVE_INFINITY;
for (final LocalisationModelSet l : localisationSets) {
for (final LocalisationModel m : l.getLocalisations()) {
final double z = m.getZ();
if (minZ > z) {
minZ = z;
}
if (maxZ < z) {
maxZ = z;
}
}
}
final int nSlices = imp.getStackSize();
// z-centre should be an index and not the ImageJ slice number so subtract 1
final int zCentre = psfSettings.getCentreImage() - 1;
// Calculate the start/end slices to cover the depth of field
// This logic must match the ImagePSFModel.
final double unitsPerSlice = psfSettings.getPixelDepth() / settings.getPixelPitch();
// We assume the PSF was imaged axially with increasing z-stage position (moving the stage
// closer to the objective). Thus higher z-coordinate are for higher slice numbers.
int lower = (int) Math.round(minZ / unitsPerSlice) + zCentre;
int upper = (int) Math.round(maxZ / unitsPerSlice) + zCentre;
// Add extra to the range so that gradients can be computed.
lower--;
upper++;
upper = (upper < 0) ? 0 : (upper >= nSlices) ? nSlices - 1 : upper;
lower = (lower < 0) ? 0 : (lower >= nSlices) ? nSlices - 1 : lower;
// Image PSF requires the z-centre for normalisation
if (!(lower <= zCentre && upper >= zCentre)) {
// Ensure we include the zCentre
lower = Math.min(lower, zCentre);
upper = Math.max(upper, zCentre);
}
final double noiseFraction = 1e-3;
final float[][] image = extractImageStack(imp, lower, upper);
final ImagePsfModel model = new ImagePsfModel(image, zCentre - lower, psfSettings.getPixelSize() / settings.getPixelPitch(), unitsPerSlice, noiseFraction);
// Add the calibrated centres. The map will not be null
final Map<Integer, Offset> map = psfSettings.getOffsetsMap();
if (!map.isEmpty()) {
final int sliceOffset = lower + 1;
for (final Entry<Integer, Offset> entry : map.entrySet()) {
model.setRelativeCentre(entry.getKey() - sliceOffset, entry.getValue().getCx(), entry.getValue().getCy());
}
} else {
// Use the CoM if present
final double cx = psfSettings.getXCentre();
final double cy = psfSettings.getYCentre();
if (cx != 0 || cy != 0) {
for (int slice = 0; slice < image.length; slice++) {
model.setCentre(slice, cx, cy);
}
}
}
// Initialise the HWHM table so that it can be cloned
model.initialiseHwhm();
return model;
} catch (final Exception ex) {
IJ.error(TITLE, "Unable to create the image PSF model:\n" + ex.getMessage());
return null;
}
}
use of uk.ac.sussex.gdsc.smlm.model.ImagePsfModel in project GDSC-SMLM by aherbert.
the class PsfDrift method createImagePsf.
private ImagePsfModel createImagePsf(int lower, int upper, double scale) {
final int zCentre = psfSettings.getCentreImage();
final double unitsPerPixel = 1.0 / scale;
// So we can move from -depth to depth
final double unitsPerSlice = 1;
// Extract data uses index not slice number as arguments so subtract 1
final double noiseFraction = 1e-3;
final float[][] image = CreateData.extractImageStack(imp, lower - 1, upper - 1);
final ImagePsfModel model = new ImagePsfModel(image, zCentre - lower, unitsPerPixel, unitsPerSlice, noiseFraction);
// Add the calibrated centres
final Map<Integer, Offset> oldOffset = psfSettings.getOffsetsMap();
if (settings.useOffset && !oldOffset.isEmpty()) {
final int sliceOffset = lower;
for (final Entry<Integer, Offset> entry : oldOffset.entrySet()) {
model.setRelativeCentre(entry.getKey() - sliceOffset, entry.getValue().getCx(), entry.getValue().getCy());
}
} else {
// Use the CoM if present
final double cx = psfSettings.getXCentre();
final double cy = psfSettings.getYCentre();
if (cx != 0 || cy != 0) {
for (int slice = 0; slice < image.length; slice++) {
model.setCentre(slice, cx, cy);
}
}
}
return model;
}
use of uk.ac.sussex.gdsc.smlm.model.ImagePsfModel in project GDSC-SMLM by aherbert.
the class PsfDrift method computeDrift.
private void computeDrift() {
// Create a grid of XY offset positions between 0-1 for PSF insert
final double[] grid = new double[settings.gridSize];
for (int i = 0; i < grid.length; i++) {
grid[i] = (double) i / settings.gridSize;
}
// Configure fitting region
final int w = 2 * settings.regionSize + 1;
centrePixel = w / 2;
// Check region size using the image PSF
final double newPsfWidth = imp.getWidth() / settings.scale;
if (Math.ceil(newPsfWidth) > w) {
ImageJUtils.log(TITLE + ": Fitted region size (%d) is smaller than the scaled PSF (%.1f)", w, newPsfWidth);
}
// Create robust PSF fitting settings
final double a = psfSettings.getPixelSize() * settings.scale;
final double sa = PsfCalculator.squarePixelAdjustment(psfSettings.getPixelSize() * (psfSettings.getFwhm() / Gaussian2DFunction.SD_TO_FWHM_FACTOR), a);
fitConfig.setInitialPeakStdDev(sa / a);
fitConfig.setBackgroundFitting(settings.backgroundFitting);
fitConfig.setNotSignalFitting(false);
fitConfig.setComputeDeviations(false);
fitConfig.setDisableSimpleFilter(true);
// Create the PSF over the desired z-depth
final int depth = (int) Math.round(settings.zDepth / psfSettings.getPixelDepth());
int startSlice = psfSettings.getCentreImage() - depth;
int endSlice = psfSettings.getCentreImage() + depth;
final int nSlices = imp.getStackSize();
startSlice = MathUtils.clip(1, nSlices, startSlice);
endSlice = MathUtils.clip(1, nSlices, endSlice);
final ImagePsfModel psf = createImagePsf(startSlice, endSlice, settings.scale);
final int minz = startSlice - psfSettings.getCentreImage();
final int maxz = endSlice - psfSettings.getCentreImage();
final int nZ = maxz - minz + 1;
final int gridSize2 = grid.length * grid.length;
total = nZ * gridSize2;
// Store all the fitting results
final int nStartPoints = getNumberOfStartPoints();
results = new double[total * nStartPoints][];
// TODO - Add ability to iterate this, adjusting the current offset in the PSF
// each iteration
// Create a pool of workers
final int threadCount = Prefs.getThreads();
final Ticker ticker = ImageJUtils.createTicker(total, threadCount, "Fitting...");
final BlockingQueue<Job> jobs = new ArrayBlockingQueue<>(threadCount * 2);
final List<Thread> threads = new LinkedList<>();
for (int i = 0; i < threadCount; i++) {
final Worker worker = new Worker(jobs, psf, w, fitConfig, ticker);
final Thread t = new Thread(worker);
threads.add(t);
t.start();
}
// Fit
outer: for (int z = minz, i = 0; z <= maxz; z++) {
for (int x = 0; x < grid.length; x++) {
for (int y = 0; y < grid.length; y++, i++) {
if (IJ.escapePressed()) {
break outer;
}
put(jobs, new Job(z, grid[x], grid[y], i));
}
}
}
// If escaped pressed then do not need to stop the workers, just return
if (ImageJUtils.isInterrupted()) {
ImageJUtils.finished();
return;
}
// Finish all the worker threads by passing in a null job
for (int i = 0; i < threads.size(); i++) {
put(jobs, new Job());
}
// Wait for all to finish
for (int i = 0; i < threads.size(); i++) {
try {
threads.get(i).join();
} catch (final InterruptedException ex) {
Thread.currentThread().interrupt();
throw new ConcurrentRuntimeException("Unexpected interrupt", ex);
}
}
threads.clear();
ImageJUtils.finished();
// Plot the average and SE for the drift curve
// Plot the recall
final double[] zPosition = new double[nZ];
final double[] avX = new double[nZ];
final double[] seX = new double[nZ];
final double[] avY = new double[nZ];
final double[] seY = new double[nZ];
final double[] recall = new double[nZ];
for (int z = minz, i = 0; z <= maxz; z++, i++) {
final Statistics statsX = new Statistics();
final Statistics statsY = new Statistics();
for (int s = 0; s < nStartPoints; s++) {
int resultPosition = i * gridSize2 + s * total;
final int endResultPosition = resultPosition + gridSize2;
while (resultPosition < endResultPosition) {
if (results[resultPosition] != null) {
statsX.add(results[resultPosition][0]);
statsY.add(results[resultPosition][1]);
}
resultPosition++;
}
}
zPosition[i] = z * psfSettings.getPixelDepth();
avX[i] = statsX.getMean();
seX[i] = statsX.getStandardError();
avY[i] = statsY.getMean();
seY[i] = statsY.getStandardError();
recall[i] = (double) statsX.getN() / (nStartPoints * gridSize2);
}
// Find the range from the z-centre above the recall limit
int centre = 0;
for (int slice = startSlice, i = 0; slice <= endSlice; slice++, i++) {
if (slice == psfSettings.getCentreImage()) {
centre = i;
break;
}
}
if (recall[centre] < settings.recallLimit) {
return;
}
int start = centre;
int end = centre;
for (int i = centre; i-- > 0; ) {
if (recall[i] < settings.recallLimit) {
break;
}
start = i;
}
for (int i = centre; ++i < recall.length; ) {
if (recall[i] < settings.recallLimit) {
break;
}
end = i;
}
final int iterations = 1;
LoessInterpolator loess = null;
if (settings.smoothing > 0) {
loess = new LoessInterpolator(settings.smoothing, iterations);
}
final double[][] smoothx = displayPlot("Drift X", "X (nm)", zPosition, avX, seX, loess, start, end);
final double[][] smoothy = displayPlot("Drift Y", "Y (nm)", zPosition, avY, seY, loess, start, end);
displayPlot("Recall", "Recall", zPosition, recall, null, null, start, end);
windowOrganiser.tile();
// Ask the user if they would like to store them in the image
final GenericDialog gd = new GenericDialog(TITLE);
gd.enableYesNoCancel();
gd.hideCancelButton();
startSlice = psfSettings.getCentreImage() - (centre - start);
endSlice = psfSettings.getCentreImage() + (end - centre);
ImageJUtils.addMessage(gd, "Save the drift to the PSF?\n \nSlices %d (%s nm) - %d (%s nm) above recall limit", startSlice, MathUtils.rounded(zPosition[start]), endSlice, MathUtils.rounded(zPosition[end]));
gd.addMessage("Optionally average the end points to set drift outside the limits.\n" + "(Select zero to ignore)");
gd.addSlider("Number_of_points", 0, 10, settings.positionsToAverage);
gd.showDialog();
if (gd.wasOKed()) {
settings.positionsToAverage = Math.abs((int) gd.getNextNumber());
final Map<Integer, Offset> oldOffset = psfSettings.getOffsetsMap();
final boolean useOldOffset = settings.useOffset && !oldOffset.isEmpty();
final LocalList<double[]> offset = new LocalList<>();
final double pitch = psfSettings.getPixelSize();
int index = 0;
for (int i = start, slice = startSlice; i <= end; slice++, i++) {
index = findCentre(zPosition[i], smoothx, index);
if (index == -1) {
ImageJUtils.log("Failed to find the offset for depth %.2f", zPosition[i]);
continue;
}
// The offset should store the difference to the centre in pixels so divide by the pixel
// pitch
double cx = smoothx[1][index] / pitch;
double cy = smoothy[1][index] / pitch;
if (useOldOffset) {
final Offset o = oldOffset.get(slice);
if (o != null) {
cx += o.getCx();
cy += o.getCy();
}
}
offset.add(new double[] { slice, cx, cy });
}
addMissingOffsets(startSlice, endSlice, nSlices, offset);
final Offset.Builder offsetBuilder = Offset.newBuilder();
final ImagePSF.Builder imagePsfBuilder = psfSettings.toBuilder();
for (final double[] o : offset) {
final int slice = (int) o[0];
offsetBuilder.setCx(o[1]);
offsetBuilder.setCy(o[2]);
imagePsfBuilder.putOffsets(slice, offsetBuilder.build());
}
imagePsfBuilder.putNotes(TITLE, String.format("Solver=%s, Region=%d", PeakFit.getSolverName(fitConfig), settings.regionSize));
imp.setProperty("Info", ImagePsfHelper.toString(imagePsfBuilder));
}
}
use of uk.ac.sussex.gdsc.smlm.model.ImagePsfModel in project GDSC-SMLM by aherbert.
the class PsfDrift method showHwhm.
private void showHwhm() {
// Build a list of suitable images
final List<String> titles = createImageList(false);
if (titles.isEmpty()) {
IJ.error(TITLE, "No suitable PSF images");
return;
}
final GenericDialog gd = new GenericDialog(TITLE);
gd.addMessage("Approximate the volume of the PSF as a Gaussian and\n" + "compute the equivalent Gaussian width.");
settings = Settings.load();
gd.addChoice("PSF", titles.toArray(new String[0]), settings.title);
gd.addCheckbox("Use_offset", settings.useOffset);
gd.addSlider("Smoothing", 0, 0.5, settings.smoothing);
gd.addHelp(HelpUrls.getUrl("psf-hwhm"));
gd.showDialog();
if (gd.wasCanceled()) {
return;
}
settings.title = gd.getNextChoice();
settings.useOffset = gd.getNextBoolean();
settings.smoothing = gd.getNextNumber();
settings.save();
imp = WindowManager.getImage(settings.title);
if (imp == null) {
IJ.error(TITLE, "No PSF image for image: " + settings.title);
return;
}
psfSettings = getPsfSettings(imp);
if (psfSettings == null) {
IJ.error(TITLE, "No PSF settings for image: " + settings.title);
return;
}
final int size = imp.getStackSize();
final ImagePsfModel psf = createImagePsf(1, size, 1);
final double[] w0 = psf.getAllHwhm0();
final double[] w1 = psf.getAllHwhm1();
// Get current centre
final int centre = psfSettings.getCentreImage();
// Extract valid values (some can be NaN)
double[] sw0 = new double[w0.length];
double[] sw1 = new double[w1.length];
final TDoubleArrayList s0 = new TDoubleArrayList(w0.length);
final TDoubleArrayList s1 = new TDoubleArrayList(w0.length);
int c0 = 0;
int c1 = 0;
for (int i = 0; i < w0.length; i++) {
if (Double.isFinite(w0[i])) {
s0.add(i + 1);
sw0[c0++] = w0[i];
}
if (Double.isFinite(w1[i])) {
s1.add(i + 1);
sw1[c1++] = w1[i];
}
}
if (c0 == 0 && c1 == 0) {
IJ.error(TITLE, "No computed HWHM for image: " + settings.title);
return;
}
double[] slice0 = s0.toArray();
sw0 = Arrays.copyOf(sw0, c0);
double[] slice1 = s1.toArray();
sw1 = Arrays.copyOf(sw1, c1);
// Smooth
if (settings.smoothing > 0) {
final LoessInterpolator loess = new LoessInterpolator(settings.smoothing, 1);
sw0 = loess.smooth(slice0, sw0);
sw1 = loess.smooth(slice1, sw1);
}
final TDoubleArrayList minWx = new TDoubleArrayList();
final TDoubleArrayList minWy = new TDoubleArrayList();
for (int i = 0; i < w0.length; i++) {
double weight = 0;
if (Double.isFinite(w0[i])) {
if (Double.isFinite(w1[i])) {
weight = w0[i] * w1[i];
} else {
weight = w0[i] * w0[i];
}
} else if (Double.isFinite(w1[i])) {
weight = w1[i] * w1[i];
}
if (weight != 0) {
minWx.add(i + 1);
minWy.add(Math.sqrt(weight));
}
}
// Smooth the combined line
final double[] cx = minWx.toArray();
double[] cy = minWy.toArray();
if (settings.smoothing > 0) {
final LoessInterpolator loess = new LoessInterpolator(settings.smoothing, 1);
cy = loess.smooth(cx, cy);
}
final int newCentre = SimpleArrayUtils.findMinIndex(cy);
// Convert to FWHM
final double fwhm = psfSettings.getFwhm();
// Widths are in pixels
final String title = TITLE + " HWHM";
final Plot plot = new Plot(title, "Slice", "HWHM (px)");
double[] limits = MathUtils.limits(sw0);
limits = MathUtils.limits(limits, sw1);
final double maxY = limits[1] * 1.05;
plot.setLimits(1, size, 0, maxY);
plot.setColor(Color.red);
plot.addPoints(slice0, sw0, Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(slice1, sw1, Plot.LINE);
plot.setColor(Color.magenta);
plot.addPoints(cx, cy, Plot.LINE);
plot.setColor(Color.black);
plot.addLabel(0, 0, "X=red; Y=blue, Combined=Magenta");
final PlotWindow pw = ImageJUtils.display(title, plot);
// Show a non-blocking dialog to allow the centre to be updated ...
// Add a label and dynamically update when the centre is moved.
final NonBlockingExtendedGenericDialog gd2 = new NonBlockingExtendedGenericDialog(TITLE);
final double scale = psfSettings.getPixelSize();
// @formatter:off
ImageJUtils.addMessage(gd2, "Update the PSF information?\n \n" + "Current z-centre = %d, FHWM = %s px (%s nm)\n", centre, MathUtils.rounded(fwhm), MathUtils.rounded(fwhm * scale));
// @formatter:on
gd2.addSlider("z-centre", cx[0], cx[cx.length - 1], newCentre);
final TextField tf = gd2.getLastTextField();
gd2.addMessage("");
gd2.addAndGetButton("Reset", event -> tf.setText(Integer.toString(newCentre)));
final Label label = gd2.getLastLabel();
gd2.addCheckbox("Update_centre", settings.updateCentre);
gd2.addCheckbox("Update_HWHM", settings.updateHwhm);
gd2.enableYesNoCancel();
gd2.hideCancelButton();
final UpdateDialogListener dl = new UpdateDialogListener(cx, cy, maxY, newCentre, scale, pw, label);
gd2.addDialogListener(dl);
gd.addHelp(HelpUrls.getUrl("psf-hwhm"));
gd2.showDialog();
if (gd2.wasOKed() && (settings.updateCentre || settings.updateHwhm)) {
final ImagePSF.Builder b = psfSettings.toBuilder();
if (settings.updateCentre) {
b.setCentreImage(dl.centre);
}
if (settings.updateHwhm) {
b.setFwhm(dl.getFwhm());
}
imp.setProperty("Info", ImagePsfHelper.toString(b));
}
}
Aggregations