Examples with PsfModel uk.ac.sussex.gdsc.smlm.model.PsfModel used on opensource projects

Example 1 with PsfModel

use of uk.ac.sussex.gdsc.smlm.model.PsfModel in project GDSC-SMLM by aherbert.

the class CreateData method drawImage.

// StoredDataStatistics rawPhotons = new StoredDataStatistics();
// StoredDataStatistics drawPhotons = new StoredDataStatistics();
// private synchronized void addRaw(double d)
// {
// //rawPhotons.add(d);
// }
// 
// private synchronized void addDraw(double d)
// {
// //drawPhotons.add(d);
// }
/**
 * Create an image from the localisations using the configured PSF width. Draws a new stack image.
 *
 * <p>Note that the localisations are filtered using the signal. The input list of localisations
 * will be updated.
 *
 * @param localisationSets the localisation sets
 * @return The localisations
 */
private List<LocalisationModel> drawImage(final List<LocalisationModelSet> localisationSets) {
    if (localisationSets.isEmpty()) {
        return null;
    }
    // Create a new list for all localisation that are drawn (i.e. pass the signal filters)
    List<LocalisationModelSet> newLocalisations = Collections.synchronizedList(new ArrayList<LocalisationModelSet>(localisationSets.size()));
    photonsRemoved = new AtomicInteger();
    removedT1 = new AtomicInteger();
    removedTn = new AtomicInteger();
    photonStats = new SummaryStatistics();
    // Add drawn spots to memory
    results = new MemoryPeakResults();
    final CalibrationWriter c = new CalibrationWriter();
    c.setDistanceUnit(DistanceUnit.PIXEL);
    c.setIntensityUnit(IntensityUnit.PHOTON);
    c.setNmPerPixel(settings.getPixelPitch());
    c.setExposureTime(settings.getExposureTime());
    c.setCameraType(settings.getCameraType());
    if (settings.getCameraType() == CameraType.SCMOS) {
        c.setCameraModelName(settings.getCameraModelName());
    } else {
        final CreateDataSettingsHelper helper = new CreateDataSettingsHelper(settings);
        c.setCountPerPhoton(helper.getTotalGainSafe());
        c.setBias(settings.getBias());
        c.setReadNoise(helper.getReadNoiseInCounts());
        c.setQuantumEfficiency(helper.getQuantumEfficiency());
    }
    results.setCalibration(c.getCalibration());
    results.setSortAfterEnd(true);
    results.begin();
    maxT = localisationSets.get(localisationSets.size() - 1).getTime();
    // Display image
    final ImageStack stack = new ImageStack(settings.getSize(), settings.getSize(), maxT);
    final double psfSd = getPsfSd();
    if (psfSd <= 0) {
        return null;
    }
    final PsfModel psfModel = createPsfModel(localisationSets);
    if (psfModel == null) {
        return null;
    }
    // Create the camera noise model
    createPerPixelCameraModelData(cameraModel);
    createBackgroundPixels();
    final int threadCount = Prefs.getThreads();
    IJ.showStatus("Drawing image ...");
    // Multi-thread for speed
    final PeakResults syncResults = SynchronizedPeakResults.create(results, threadCount);
    final ExecutorService threadPool = Executors.newFixedThreadPool(threadCount);
    final List<Future<?>> futures = new LinkedList<>();
    // Count all the frames to process
    final Ticker ticker = ImageJUtils.createTicker(maxT, threadCount);
    // Collect statistics on the number of photons actually simulated
    // Process all frames
    int index = 0;
    int lastT = -1;
    for (final LocalisationModelSet l : localisationSets) {
        if (ImageJUtils.isInterrupted()) {
            break;
        }
        if (l.getTime() != lastT) {
            lastT = l.getTime();
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, index, lastT, psfModel.copy(), syncResults, stack, poissonNoise, createRandomGenerator(), ticker)));
        }
        index++;
    }
    // Finish processing data
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    futures.clear();
    if (ImageJUtils.isInterrupted()) {
        IJ.showProgress(1);
        return null;
    }
    // Do all the frames that had no localisations
    float[] limits = null;
    for (int t = 1; t <= maxT; t++) {
        if (ImageJUtils.isInterrupted()) {
            break;
        }
        final Object pixels = stack.getPixels(t);
        if (pixels == null) {
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, maxT, t, null, syncResults, stack, poissonNoise, createRandomGenerator(), ticker)));
        } else if (limits == null) {
            limits = MathUtils.limits((float[]) pixels);
        }
    }
    // Finish
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    futures.clear();
    threadPool.shutdown();
    IJ.showProgress(1);
    if (ImageJUtils.isInterrupted() || limits == null) {
        return null;
    }
    results.end();
    if (photonsRemoved.get() > 0) {
        ImageJUtils.log("Removed %d localisations with less than %.1f rendered photons", photonsRemoved.get(), settings.getMinPhotons());
    }
    if (removedT1.get() > 0) {
        ImageJUtils.log("Removed %d localisations with no neighbours @ SNR %.2f", removedT1.get(), settings.getMinSnrT1());
    }
    if (removedTn.get() > 0) {
        ImageJUtils.log("Removed %d localisations with valid neighbours @ SNR %.2f", removedTn.get(), settings.getMinSnrTN());
    }
    if (photonStats.getN() > 0) {
        ImageJUtils.log("Average photons rendered = %s +/- %s", MathUtils.rounded(photonStats.getMean()), MathUtils.rounded(photonStats.getStandardDeviation()));
    }
    // System.out.printf("rawPhotons = %f\n", rawPhotons.getMean());
    // System.out.printf("drawPhotons = %f\n", drawPhotons.getMean());
    // new HistogramPlotBuilder("draw photons", drawPhotons, "photons", true, 0, 1000);
    // Update with all those localisation that have been drawn
    localisationSets.clear();
    localisationSets.addAll(newLocalisations);
    newLocalisations = null;
    IJ.showStatus("Displaying image ...");
    ImageStack newStack = stack;
    if (!settings.getRawImage()) {
        // Get the global limits and ensure all values can be represented
        final Object[] imageArray = stack.getImageArray();
        limits = MathUtils.limits((float[]) imageArray[0]);
        for (int j = 1; j < imageArray.length; j++) {
            limits = MathUtils.limits(limits, (float[]) imageArray[j]);
        }
        // float limits0 = limits[0];
        // Leave bias in place
        final float limits0 = 0;
        // Check if the image will fit in a 16-bit range
        if ((limits[1] - limits0) < 65535) {
            // Convert to 16-bit
            newStack = new ImageStack(stack.getWidth(), stack.getHeight(), stack.getSize());
            // Account for rounding
            final float min = (float) (limits0 - 0.5);
            for (int j = 0; j < imageArray.length; j++) {
                final float[] image = (float[]) imageArray[j];
                final short[] pixels = new short[image.length];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = (short) (image[k] - min);
                }
                newStack.setPixels(pixels, j + 1);
                // Free memory
                imageArray[j] = null;
                // Attempt to stay within memory (check vs 32MB)
                if (MemoryUtils.getFreeMemory() < 33554432L) {
                    MemoryUtils.runGarbageCollectorOnce();
                }
            }
            for (int k = 2; k-- > 0; ) {
                limits[k] = (float) Math.floor(limits[k] - min);
            }
        } else {
            // Keep as 32-bit but round to whole numbers
            for (int j = 0; j < imageArray.length; j++) {
                final float[] pixels = (float[]) imageArray[j];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = Math.round(pixels[k]);
                }
            }
            for (int k = 2; k-- > 0; ) {
                limits[k] = Math.round(limits[k]);
            }
        }
    }
    // Show image
    final ImagePlus imp = ImageJUtils.display(CREATE_DATA_IMAGE_TITLE, newStack);
    final ij.measure.Calibration cal = new ij.measure.Calibration();
    String unit = "nm";
    double unitPerPixel = settings.getPixelPitch();
    if (unitPerPixel > 100) {
        unit = "um";
        unitPerPixel /= 1000.0;
    }
    cal.setUnit(unit);
    cal.pixelHeight = cal.pixelWidth = unitPerPixel;
    final Rectangle bounds = cameraModel.getBounds();
    if (bounds != null) {
        cal.xOrigin = -bounds.x;
        cal.yOrigin = -bounds.y;
    }
    imp.setCalibration(cal);
    imp.setDimensions(1, 1, newStack.getSize());
    imp.setDisplayRange(limits[0], limits[1]);
    imp.updateAndDraw();
    saveImage(imp);
    final IJImageSource imageSource = new IJImageSource(imp);
    // Shift simulation image source to correct location
    results.setSource(imageSource);
    results.setName(CREATE_DATA_IMAGE_TITLE + " (" + TITLE + ")");
    // Bounds are relative to the image source
    results.setBounds(new Rectangle(settings.getSize(), settings.getSize()));
    results.setPsf(createPsf(psfSd));
    MemoryPeakResults.addResults(results);
    setBenchmarkResults(imp, results);
    if (benchmarkMode && benchmarkParameters != null) {
        benchmarkParameters.setPhotons(results);
    }
    final List<LocalisationModel> localisations = toLocalisations(localisationSets);
    savePulses(localisations, results);
    // Saved the fixed and moving localisations into different datasets
    saveFixedAndMoving(results);
    saveCompoundMolecules(results);
    return localisations;
}

Example 2 with PsfModel

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