Examples with StandardValueProcedure uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure used on opensource projects

Example 6 with StandardValueProcedure

use of uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionVsPsfModelTest method computesSameAsPsfModel.

private void computesSameAsPsfModel(double sum, double x0, double x1, double s0, double s1) {
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, width, height, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    a[Gaussian2DFunction.SIGNAL] = sum;
    a[Gaussian2DFunction.X_POSITION] = x0;
    a[Gaussian2DFunction.Y_POSITION] = x1;
    a[Gaussian2DFunction.X_SD] = s0;
    a[Gaussian2DFunction.Y_SD] = s1;
    final double[] o = new StandardValueProcedure().getValues(f, a);
    final GaussianPsfModel m = new GaussianPsfModel(s0, s1);
    final double[] e = new double[o.length];
    // Note that the Gaussian2DFunction has 0,0 at the centre of a pixel.
    // The model has 0.5,0.5 at the centre so add an offset.
    m.create2D(e, width, height, sum, x0 + 0.5, x1 + 0.5, null);
    for (int i = 0; i < e.length; i++) {
        // Only check where there is a reasonable amount of signal
        if (e[i] > 1e-2) {
            final double error = DoubleEquality.relativeError(e[i], o[i]);
            // uses the Apache commons implementation.
            if (error > 5e-3) {
                Assertions.fail(String.format("[%d] %s != %s  error = %f", i, Double.toString(e[i]), Double.toString(o[i]), error));
            }
        }
    }
}

Example 7 with StandardValueProcedure

use of uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure in project GDSC-SMLM by aherbert.

the class CubicSplineFunctionTest method functionComputesTargetWithAndWithoutGradientWith2Peaks.

@Test
void functionComputesTargetWithAndWithoutGradientWith2Peaks() {
    if (f2 == null) {
        return;
    }
    final StandardValueProcedure p0 = new StandardValueProcedure();
    final StandardGradient1Procedure p1 = new StandardGradient1Procedure();
    final StandardGradient2Procedure p2 = new StandardGradient2Procedure();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        // Peak 2
                        for (final double signal2 : testsignal2) {
                            for (final double cx2 : testcx2) {
                                for (final double cy2 : testcy2) {
                                    for (final double cz2 : testcz2) {
                                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1, signal2, cx2, cy2, cz2);
                                        final double[] e = p0.getValues(f1, a);
                                        final double[] o1 = p1.getValues(f1, a);
                                        final double[] o2 = p2.getValues(f1, a);
                                        Assertions.assertArrayEquals(e, o1);
                                        Assertions.assertArrayEquals(e, o2);
                                        for (int i = e.length; i-- > 0; ) {
                                            Assertions.assertArrayEquals(p1.gradients[i], p2.gradients1[i]);
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

Example 8 with StandardValueProcedure

use of uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure in project GDSC-SMLM by aherbert.

the class CubicSplineFunctionTest method functionComputesTargetWithAndWithoutGradient.

@Test
void functionComputesTargetWithAndWithoutGradient() {
    final StandardValueProcedure p0 = new StandardValueProcedure();
    final StandardGradient1Procedure p1 = new StandardGradient1Procedure();
    final StandardGradient2Procedure p2 = new StandardGradient2Procedure();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1);
                        final double[] e = p0.getValues(f1, a);
                        final double[] o1 = p1.getValues(f1, a);
                        final double[] o2 = p2.getValues(f1, a);
                        Assertions.assertArrayEquals(e, o1);
                        Assertions.assertArrayEquals(e, o2);
                        for (int i = e.length; i-- > 0; ) {
                            Assertions.assertArrayEquals(p1.gradients[i], p2.gradients1[i]);
                        }
                    }
                }
            }
        }
    }
}

Example 9 with StandardValueProcedure

use of uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure 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()));
    }
}

Example 10 with StandardValueProcedure

use of uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure in project GDSC-SMLM by aherbert.

the class DoubleDht3DTest method createData.

private static DoubleDht3D createData(double cx, double cy, double cz) {
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_ASTIGMATISM, zModel);
    final int length = size * size;
    final double[] data = new double[size * length];
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    a[Gaussian2DFunction.SIGNAL] = 1;
    a[Gaussian2DFunction.X_POSITION] = cx;
    a[Gaussian2DFunction.Y_POSITION] = cy;
    a[Gaussian2DFunction.X_SD] = 1;
    a[Gaussian2DFunction.Y_SD] = 1;
    final StandardValueProcedure p = new StandardValueProcedure();
    for (int z = 0; z < size; z++) {
        a[Gaussian2DFunction.Z_POSITION] = z - cz;
        p.getValues(f, a, data, z * length);
    }
    return new DoubleDht3D(size, size, size, data, false);
}