use of uk.ac.sussex.gdsc.smlm.model.AiryPsfModel in project GDSC-SMLM by aherbert.
the class CreateData method createPsfModel.
private PsfModel createPsfModel(List<LocalisationModelSet> localisationSets) {
// Allow reuse of the cached model
if (psfModelCache != null) {
return psfModelCache;
}
if (psfModelType == PSF_MODEL_IMAGE) {
return createImagePsf(localisationSets);
}
if (psfModelType == PSF_MODEL_ASTIGMATISM) {
astigmatismModel = AstigmatismModelManager.getModel(settings.getAstigmatismModel());
if (astigmatismModel == null) {
throw new IllegalArgumentException("Failed to load model: " + settings.getAstigmatismModel());
}
// Convert for simulation
try {
if (DoubleEquality.relativeError(astigmatismModel.getNmPerPixel(), settings.getPixelPitch()) > 1e-6) {
final String message = String.format("Astigmatism model '%s' calibration (%s nm) does not match pixel pitch (%s nm)", settings.getAstigmatismModel(), MathUtils.rounded(astigmatismModel.getNmPerPixel()), MathUtils.rounded(settings.getPixelPitch()));
// Optionally convert
final GenericDialog gd = new GenericDialog(TITLE);
gd.addMessage(TextUtils.wrap(message + ". Created data is not suitable for fitting.", 80));
gd.addMessage(TextUtils.wrap("Click OK to continue anyway (i.e. draw the spot using the " + "correct nm width on the different sized pixels).", 80));
gd.showDialog();
if (gd.wasCanceled()) {
throw new IllegalArgumentException(message);
}
// Convert to nm
astigmatismModel = AstigmatismModelManager.convert(astigmatismModel, DistanceUnit.NM, DistanceUnit.NM);
// Reset pixel pitch. This will draw the spot using the correct size on the different size
// pixels.
astigmatismModel = astigmatismModel.toBuilder().setNmPerPixel(settings.getPixelPitch()).build();
}
// Convert for simulation in pixels
astigmatismModel = AstigmatismModelManager.convert(astigmatismModel, DistanceUnit.PIXEL, DistanceUnit.PIXEL);
return new GaussianPsfModel(AstigmatismModelManager.create(astigmatismModel));
} catch (final ConversionException ex) {
// Wrap so this can be caught as the same type
throw new IllegalArgumentException(ex);
}
}
final double d = settings.getDepthOfFocus() / settings.getPixelPitch();
if (psfModelType == PSF_MODEL_GAUSSIAN) {
final double sd = getPsfSd();
final double gamma = 0;
final HoltzerAstigmatismZModel zModel = HoltzerAstigmatismZModel.create(sd, sd, gamma, d, 0, 0, 0, 0);
final GaussianPsfModel m = new GaussianPsfModel(zModel);
// m.setRange(10);
return m;
}
// Default to Airy pattern
final double width = getPsfSd() / PsfCalculator.AIRY_TO_GAUSSIAN;
final AiryPsfModel m = new AiryPsfModel(width, width, d);
m.setRing(2);
return m;
}
use of uk.ac.sussex.gdsc.smlm.model.AiryPsfModel 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()));
}
}
Aggregations