Search in sources :

Example 6 with CameraModel

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

the class SpotFinderPreview method run.

@Override
public void run(ImageProcessor ip) {
    if (refreshing) {
        return;
    }
    final Rectangle bounds = ip.getRoi();
    // Only do this if the settings changed
    final Calibration calibration = fitConfig.getCalibration();
    final FitEngineSettings fitEngineSettings = config.getFitEngineSettings();
    final PSF psf = fitConfig.getPsf();
    boolean newCameraModel = filter == null;
    if (!calibration.equals(lastCalibration)) {
        newCameraModel = true;
        // Set a camera model.
        // We have to set the camera type too to avoid configuration errors.
        CameraModel cameraModel = CameraModelManager.load(fitConfig.getCameraModelName());
        if (cameraModel == null) {
            cameraModel = new FakePerPixelCameraModel(0, 1, 1);
            fitConfig.setCameraType(CameraType.EMCCD);
        } else {
            fitConfig.setCameraType(CameraType.SCMOS);
            // Support cropped origin selection.
            final Rectangle sourceBounds = IJImageSource.getBounds(imp);
            cameraModel = PeakFit.cropCameraModel(cameraModel, sourceBounds, null, true);
            if (cameraModel == null) {
                gd.getPreviewCheckbox().setState(false);
                return;
            }
        }
        fitConfig.setCameraModel(cameraModel);
    }
    if (newCameraModel || !fitEngineSettings.equals(lastFitEngineSettings) || !psf.equals(lastPsf)) {
        // Configure a jury filter
        if (config.getDataFilterType() == DataFilterType.JURY && !PeakFit.configureDataFilter(config, PeakFit.FLAG_NO_SAVE)) {
            gd.getPreviewCheckbox().setState(false);
            return;
        }
        try {
            filter = config.createSpotFilter();
        } catch (final Exception ex) {
            filter = null;
            this.imp.setOverlay(overlay);
            // Required for ImageJ to disable the preview
            throw new IllegalStateException("Unable to create spot filter", ex);
        }
        ImageJUtils.log(filter.getDescription());
    }
    lastCalibration = calibration;
    lastFitEngineSettings = fitEngineSettings;
    lastPsf = psf;
    // This code can probably be removed since the crop is done above.
    if (fitConfig.getCameraTypeValue() == CameraType.SCMOS_VALUE) {
        // Instead just warn if the roi cannot be extracted from the selected model
        // or there is a mismatch
        final Rectangle modelBounds = fitConfig.getCameraModel().getBounds();
        if (modelBounds != null) {
            if (!modelBounds.contains(bounds)) {
                // @formatter:off
                ImageJUtils.log("WARNING: Camera model bounds [x=%d,y=%d,width=%d,height=%d]" + " does not contain image target bounds [x=%d,y=%d,width=%d,height=%d]", modelBounds.x, modelBounds.y, modelBounds.width, modelBounds.height, bounds.x, bounds.y, bounds.width, bounds.height);
            // @formatter:on
            // Warn if the model bounds are mismatched than the image as this may be an incorrect
            // selection for the camera model
            } else if (modelBounds.x != 0 || modelBounds.y != 0 || modelBounds.width > ip.getWidth() || modelBounds.height > ip.getHeight()) {
                // @formatter:off
                ImageJUtils.log("WARNING: Probably an incorrect camera model!\n" + "Model bounds [x=%d,y=%d,width=%d,height=%d]\n" + "do not match the image target bounds [width=%d,height=%d].", modelBounds.x, modelBounds.y, modelBounds.width, modelBounds.height, ip.getWidth(), ip.getHeight());
            // @formatter:on
            }
        }
    }
    run(ip, filter);
}
Also used : PSF(uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF) FakePerPixelCameraModel(uk.ac.sussex.gdsc.smlm.model.camera.FakePerPixelCameraModel) CameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CameraModel) Rectangle(java.awt.Rectangle) Calibration(uk.ac.sussex.gdsc.smlm.data.config.CalibrationProtos.Calibration) FitEngineSettings(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.FitEngineSettings) FakePerPixelCameraModel(uk.ac.sussex.gdsc.smlm.model.camera.FakePerPixelCameraModel)

Example 7 with CameraModel

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

the class BenchmarkFit method runFit.

private void runFit() {
    // Initialise the answer.
    answer[Gaussian2DFunction.BACKGROUND] = benchmarkParameters.getBackground();
    answer[Gaussian2DFunction.SIGNAL] = benchmarkParameters.getSignal();
    answer[Gaussian2DFunction.X_POSITION] = benchmarkParameters.x;
    answer[Gaussian2DFunction.Y_POSITION] = benchmarkParameters.y;
    answer[Gaussian2DFunction.Z_POSITION] = benchmarkParameters.z;
    answer[Gaussian2DFunction.X_SD] = benchmarkParameters.sd / benchmarkParameters.pixelPitch;
    answer[Gaussian2DFunction.Y_SD] = benchmarkParameters.sd / benchmarkParameters.pixelPitch;
    // Set up the fit region. Always round down since 0.5 is the centre of the pixel.
    final int x = (int) benchmarkParameters.x;
    final int y = (int) benchmarkParameters.y;
    region = new Rectangle(x - regionSize, y - regionSize, 2 * regionSize + 1, 2 * regionSize + 1);
    if (!new Rectangle(0, 0, imp.getWidth(), imp.getHeight()).contains(region)) {
        // Check if it is incorrect by only 1 pixel
        if (region.width <= imp.getWidth() + 1 && region.height <= imp.getHeight() + 1) {
            ImageJUtils.log("Adjusting region %s to fit within image bounds (%dx%d)", region.toString(), imp.getWidth(), imp.getHeight());
            region = new Rectangle(0, 0, imp.getWidth(), imp.getHeight());
        } else {
            IJ.error(TITLE, "Fit region does not fit within the image");
            return;
        }
    }
    // Adjust the centre & account for 0.5 pixel offset during fitting
    answer[Gaussian2DFunction.X_POSITION] -= (region.x + 0.5);
    answer[Gaussian2DFunction.Y_POSITION] -= (region.y + 0.5);
    // Configure for fitting
    fitConfig.setBackgroundFitting(backgroundFitting);
    fitConfig.setNotSignalFitting(!signalFitting);
    fitConfig.setComputeDeviations(false);
    // Create the camera model
    CameraModel cameraModel = fitConfig.getCameraModel();
    // Crop for speed. Reset origin first so the region is within the model
    cameraModel.setOrigin(0, 0);
    cameraModel = cameraModel.crop(region, false);
    final ImageStack stack = imp.getImageStack();
    final int totalFrames = benchmarkParameters.frames;
    // Create a pool of workers
    final int nThreads = Prefs.getThreads();
    final BlockingQueue<Integer> jobs = new ArrayBlockingQueue<>(nThreads * 2);
    final List<Worker> workers = new LinkedList<>();
    final List<Thread> threads = new LinkedList<>();
    final Ticker ticker = ImageJUtils.createTicker(totalFrames, nThreads, "Fitting frames ...");
    for (int i = 0; i < nThreads; i++) {
        final Worker worker = new Worker(jobs, stack, region, fitConfig, cameraModel, ticker);
        final Thread t = new Thread(worker);
        workers.add(worker);
        threads.add(t);
        t.start();
    }
    // Store all the fitting results
    results = new double[totalFrames * startPoints.length][];
    resultsTime = new long[results.length];
    // Fit the frames
    for (int i = 0; i < totalFrames; i++) {
        // Only fit if there were simulated photons
        if (benchmarkParameters.framePhotons[i] > 0) {
            put(jobs, i);
        }
    }
    // Finish all the worker threads by passing in a null job
    for (int i = 0; i < threads.size(); i++) {
        put(jobs, -1);
    }
    // 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(ex);
        }
    }
    threads.clear();
    if (hasOffsetXy()) {
        ImageJUtils.log(TITLE + ": CoM within start offset = %d / %d (%s%%)", comValid.intValue(), totalFrames, MathUtils.rounded((100.0 * comValid.intValue()) / totalFrames));
    }
    ImageJUtils.finished("Collecting results ...");
    // Collect the results
    Statistics[] stats = null;
    for (int i = 0; i < workers.size(); i++) {
        final Statistics[] next = workers.get(i).stats;
        if (stats == null) {
            stats = next;
            continue;
        }
        for (int j = 0; j < next.length; j++) {
            stats[j].add(next[j]);
        }
    }
    workers.clear();
    Objects.requireNonNull(stats, "No statistics were computed");
    // Show a table of the results
    summariseResults(stats, cameraModel);
    // Optionally show histograms
    if (showHistograms) {
        IJ.showStatus("Calculating histograms ...");
        final WindowOrganiser windowOrganiser = new WindowOrganiser();
        final double[] convert = getConversionFactors();
        final HistogramPlotBuilder builder = new HistogramPlotBuilder(TITLE).setNumberOfBins(histogramBins);
        for (int i = 0; i < NAMES.length; i++) {
            if (displayHistograms[i] && convert[i] != 0) {
                // We will have to convert the values...
                final double[] tmp = ((StoredDataStatistics) stats[i]).getValues();
                for (int j = 0; j < tmp.length; j++) {
                    tmp[j] *= convert[i];
                }
                final StoredDataStatistics tmpStats = StoredDataStatistics.create(tmp);
                builder.setData(tmpStats).setName(NAMES[i]).setPlotLabel(String.format("%s +/- %s", MathUtils.rounded(tmpStats.getMean()), MathUtils.rounded(tmpStats.getStandardDeviation()))).show(windowOrganiser);
            }
        }
        windowOrganiser.tile();
    }
    if (saveRawData) {
        final String dir = ImageJUtils.getDirectory("Data_directory", rawDataDirectory);
        if (dir != null) {
            saveData(stats, dir);
        }
    }
    IJ.showStatus("");
}
Also used : CameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CameraModel) ImageStack(ij.ImageStack) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) Rectangle(java.awt.Rectangle) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) WindowOrganiser(uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) LinkedList(java.util.LinkedList) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue)

Example 8 with CameraModel

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

the class PsfCreator method runUsingAlignment.

private void runUsingAlignment() {
    if (!showAlignmentDialog()) {
        return;
    }
    boxRadius = (int) Math.ceil(settings.getRadius());
    final CalibrationReader calibration = new CalibrationReader(settings.getCalibration());
    // Limit this
    final int halfBoxRadius = boxRadius / 2;
    settings.setAnalysisWindow(Math.min(settings.getAnalysisWindow(), halfBoxRadius));
    // Find the selected PSF spots x,y,z centre
    // We offset the centre to the middle of pixel.
    BasePoint[] centres = getSpots(0.5f, false);
    if (centres.length == 0) {
        IJ.error(TITLE, "No PSFs");
        return;
    }
    CameraModel cameraModel = null;
    if (calibration.isScmos()) {
        cameraModel = CameraModelManager.load(calibration.getCameraModelName());
        if (cameraModel == null) {
            IJ.error(TITLE, "No camera model");
            return;
        }
        cameraModel = PeakFit.cropCameraModel(cameraModel, IJImageSource.getBounds(imp), null, true);
    } else {
        cameraModel = new CcdCameraModel(calibration.getBias(), 1);
    }
    // Extract the image data for processing as float
    final float[][] image = CreateData.extractImageStack(imp, 0, imp.getStackSize() - 1);
    for (final float[] data : image) {
        cameraModel.removeBiasAndGain(data);
    }
    zSelector = new PsfCentreSelector();
    // Relocate the initial centres
    ImageJUtils.showStatus("Relocating initial centres");
    centres = relocateCentres(image, centres);
    if (centres == null) {
        return;
    }
    zRadius = (int) Math.ceil(settings.getAlignmentZRadius());
    // Check the region overlap in 3D and exclude overlapping PSFs
    boolean[] bad = findSpotOverlap(centres, null);
    centres = getNonBadSpots(centres, bad);
    if (centres.length == 0) {
        IJ.error(TITLE, "No PSFs without neighbours within the box region");
        return;
    }
    // Multi-thread for speed
    if (threadPool == null) {
        threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
    }
    // Extract each PSF into a scaled PSF
    ImageJUtils.showStatus(String.format("[%d] Extracting PSFs", 0));
    ExtractedPsf[] psfs = extractPsfs(image, centres);
    Point location = null;
    // Iterate until centres have converged
    boolean converged = false;
    for (int iter = 0; !converged && iter < settings.getMaxIterations(); iter++) {
        if (ImageJUtils.isInterrupted()) {
            return;
        }
        // Combine all PSFs
        ImageJUtils.showStatus(String.format("[%d] Aligning PSFs", iter + 1));
        final ExtractedPsf combined = combine(psfs);
        combined.createProjections();
        // Get the current combined z-centre.
        // This is used to get the centre of mass for repositioning.
        // It also effects the alignment so do it for the first iteration.
        zSelector.setPsf(combined);
        if (iter == 0) {
            // TODO - check if the z-centre should be guessed here.
            // We assume that the combined PSF may be easier to guess if the initial
            // guess for each individual PSF was OK. It may not be necessary since all
            // the PSFs are combined around their z-centres. Once alignment has
            // started we skip this step.
            zSelector.analyse();
            zSelector.guessZCentre();
        }
        if (settings.getInteractiveMode()) {
            if (iter != 0) {
                zSelector.analyse();
            }
            // zSelector.guessZCentre();
            final double dz = zSelector.run("Update combined PSF z-centre", true, false, false, null);
            if (dz < 0) {
                return;
            }
        }
        // Align each to the combined PSF
        final float[][] translation = align(combined, psfs);
        if (ImageJUtils.isInterrupted()) {
            return;
        }
        // Find the new centre using the old centre plus the alignment shift
        for (int j = 0; j < psfs.length; j++) {
            centres[j] = psfs[j].updateCentre(translation[j]);
            // Update to get the correct scale
            translation[j][0] = centres[j].getX() - psfs[j].centre.getX();
            translation[j][1] = centres[j].getY() - psfs[j].centre.getY();
            translation[j][2] = centres[j].getZ() - psfs[j].centre.getZ();
            ImageJUtils.log("[%d] Centre %d : Shift X = %s : Shift Y = %s : Shift Z = %s", iter, j + 1, rounder.toString(translation[j][0]), rounder.toString(translation[j][1]), rounder.toString(translation[j][2]));
        }
        final boolean[] excluded = new boolean[psfs.length];
        if (checkAlignments) {
            combined.show(TITLE_PSF);
            // Ask about each centre in turn.
            // Update Point ROI using float coordinates and set image slice to
            // correct z-centre.
            // imp.saveRoi();
            imp.killRoi();
            final ImageCanvas ic = imp.getCanvas();
            // ic.setMagnification(16);
            int reject = 0;
            final float box = boxRadius + 0.5f;
            final int n = imp.getStackSize();
            for (int j = 0; j < centres.length; j++) {
                psfs[j].show(TITLE_SPOT_PSF);
                final Overlay o = new Overlay();
                o.add(createRoi(psfs[j].centre.getX(), psfs[j].centre.getY(), Color.RED));
                final float cx = centres[j].getX();
                final float cy = centres[j].getY();
                o.add(createRoi(cx, cy, Color.GREEN));
                final Roi roi = new Roi(cx - box, cy - box, 2 * box, 2 * box);
                o.add(roi);
                // The centre is absolute within the original stack
                imp.setSlice(MathUtils.clip(1, n, Math.round(centres[j].getZ())));
                final Rectangle r = ic.getSrcRect();
                final int x = centres[j].getXint();
                final int y = centres[j].getYint();
                if (!r.contains(x, y)) {
                    r.x = x - r.width / 2;
                    r.y = y - r.height / 2;
                    ic.setSourceRect(r);
                }
                imp.setOverlay(o);
                imp.updateAndDraw();
                final NonBlockingExtendedGenericDialog gd = new NonBlockingExtendedGenericDialog(TITLE);
                ImageJUtils.addMessage(gd, "Shift X = %s\nShift Y = %s\nShift Z = %s", rounder.toString(translation[j][0]), rounder.toString(translation[j][1]), rounder.toString(translation[j][2]));
                final int spotIndex = j;
                gd.addAndGetButton("Exclude spot", event -> {
                    if (excluded[spotIndex]) {
                        ImageJUtils.log("Included spot %d", spotIndex + 1);
                        excluded[spotIndex] = false;
                    } else {
                        ImageJUtils.log("Excluded spot %d", spotIndex + 1);
                        excluded[spotIndex] = true;
                    }
                });
                gd.enableYesNoCancel("Accept", "Reject");
                if (location != null) {
                    gd.setLocation(location.x, location.y);
                }
                gd.showDialog();
                if (gd.wasCanceled()) {
                    resetImp();
                    return;
                }
                final boolean failed = excluded[spotIndex] || !gd.wasOKed();
                if (failed) {
                    reject++;
                    centres[j] = psfs[j].centre;
                    // For RMSD computation
                    Arrays.fill(translation[j], 0f);
                }
                location = gd.getLocation();
            }
            resetImp();
            if (reject == psfs.length) {
                IJ.error(TITLE, "No PSF translations were accepted");
                return;
            }
        }
        bad = findSpotOverlap(centres, excluded);
        final int badCount = count(bad);
        final int excludedCount = count(excluded);
        int ok = bad.length - badCount - excludedCount;
        if (ok < bad.length) {
            if (badCount != 0 && settings.getInteractiveMode()) {
                final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
                gd.addMessage("Warning: Regions now overlap!");
                gd.addMessage("OK = " + TextUtils.pleural(ok, "PSF"));
                gd.addMessage("Overlapping = " + TextUtils.pleural(badCount, "PSF"));
                // gd.addMessage("Excluded = " + TextUtils.pleural(excludedCount, "PSF"));
                gd.enableYesNoCancel("Exclude", "Include");
                if (location != null) {
                    gd.setLocation(location.x, location.y);
                }
                gd.showDialog();
                if (gd.wasCanceled()) {
                    resetImp();
                    return;
                }
                if (!gd.wasOKed()) {
                    // allow bad spots
                    Arrays.fill(bad, false);
                    ok = bad.length;
                }
                location = gd.getLocation();
            }
            if (ok == 0) {
                IJ.error(TITLE, "No PSFs remaining");
                resetImp();
                return;
            }
        }
        // Merge bad and excluded to get new centres
        for (int i = 0; i < bad.length; i++) {
            if (excluded[i]) {
                bad[i] = true;
            }
        }
        ok = bad.length - count(bad);
        final BasePoint[] newCentres = getNonBadSpots(centres, bad);
        // Find the change in centres
        final double[] rmsd = new double[2];
        for (int j = 0; j < psfs.length; j++) {
            if (bad[j]) {
                continue;
            }
            rmsd[0] += MathUtils.pow2(translation[j][0]) + MathUtils.pow2(translation[j][1]);
            rmsd[1] += MathUtils.pow2(translation[j][2]);
        }
        for (int j = 0; j < 2; j++) {
            rmsd[j] = Math.sqrt(rmsd[j] / ok);
        }
        ImageJUtils.showStatus(String.format("[%d] Checking combined PSF", iter + 1));
        // Compute CoM shift using the current z-centre and z-window
        final double[] shift = combined.getCentreOfMassXyShift(zSelector.getCentreSlice());
        final double shiftd = Math.sqrt(shift[0] * shift[0] + shift[1] * shift[1]);
        ImageJUtils.log("[%d] RMSD XY = %s : RMSD Z = %s : Combined CoM shift = %s,%s (%s)", iter, rounder.toString(rmsd[0]), rounder.toString(rmsd[1]), rounder.toString(shift[0]), rounder.toString(shift[1]), rounder.toString(shiftd));
        if (settings.getInteractiveMode()) {
            // Ask if OK to continue?
            final GenericDialog gd = new GenericDialog(TITLE);
            ImageJUtils.addMessage(gd, "RMSD XY = %s\nRMSD Z = %s\nCombined CoM shift = %s,%s (%s)", rounder.toString(rmsd[0]), rounder.toString(rmsd[1]), rounder.toString(shift[0]), rounder.toString(shift[1]), rounder.toString(shiftd));
            // Check if we can do more iterations
            if (iter + 1 < settings.getMaxIterations()) {
                gd.enableYesNoCancel("Continue", "Converged");
            } else {
                gd.setOKLabel("Converged");
            }
            gd.showDialog();
            if (gd.wasCanceled()) {
                return;
            }
            converged = !gd.wasOKed();
        } else {
            // Check convergence thresholds
            converged = rmsd[0] < settings.getRmsdXyThreshold() && rmsd[1] < settings.getRmsdZThreshold() && shiftd < settings.getComShiftThreshold();
        }
        // For the next round we move to the non-overlapping spots
        centres = newCentres;
        // Update the centres using the centre-of-mass of the combined PSF
        centres = updateUsingCentreOfMassXyShift(shift, shiftd, combined, centres);
        // Extract each PSF into a scaled PSF
        ImageJUtils.showStatus(String.format("[%d] Extracting PSFs", iter + 1));
        psfs = extractPsfs(image, centres);
    }
    // Combine all
    ExtractedPsf combined = combine(psfs);
    // Show an interactive dialog for cropping the PSF and choosing the
    // final output
    final PsfOutputSelector cropSelector = new PsfOutputSelector(combined);
    combined = cropSelector.run();
    if (combined == null) {
        return;
    }
    if (settings.getUpdateRoi()) {
        final float[] ox = new float[centres.length];
        final float[] oy = new float[centres.length];
        for (int i = 0; i < centres.length; i++) {
            ox[i] = centres[i].getX();
            oy[i] = centres[i].getY();
        }
        imp.setRoi(new OffsetPointRoi(ox, oy));
    }
    // For an image PSF we can just enlarge the PSF and window.
    // For a CSpline then we already have the 3D cubic spline function.
    // However we want to post-process the function to allow windowing and
    // normalisation. So we enlarge by 3 in each dimension.
    // The CSpline can be created by solving the coefficients for the
    // 4x4x4 (64) sampled points on each node.
    int magnification;
    if (settings.getOutputType() == OUTPUT_TYPE_IMAGE_PSF) {
        magnification = settings.getPsfMagnification();
    } else {
        magnification = 3;
    }
    // Enlarge the combined PSF for final processing
    ExtractedPsf finalPsf = combined.enlarge(magnification, threadPool);
    // Show a dialog to collect final z-centre interactively
    ImageJUtils.showStatus("Analysing PSF");
    zSelector.setPsf(finalPsf);
    zSelector.analyse();
    // zSelector.guessZCentre(); // No need to guess the centre
    final double dz = zSelector.run("Finalise PSF", true, true, true, null);
    if (dz < 0) {
        return;
    }
    zCentre = zSelector.getCentreSlice();
    if (settings.getCropToZCentre()) {
        finalPsf = finalPsf.cropToZCentre(zCentre);
        // Back to 1-based index
        zCentre = finalPsf.stackZCentre + 1;
    }
    // When click ok the background is subtracted from the PSF
    // All pixels below the background are set to zero
    // Apply a Tukey window to roll-off to zero at the outer pixels
    ImageJUtils.showStatus("Windowing PSF");
    final double[] wx = ImageWindow.tukeyEdge(finalPsf.maxx, settings.getWindow());
    final double[] wz = ImageWindow.tukeyEdge(finalPsf.psf.length, settings.getWindow());
    // Normalisation so the max intensity frame is one
    final float[][] psf = finalPsf.psf;
    final int maxz = psf.length;
    final double[] sum = new double[maxz];
    for (int z = 0; z < maxz; z++) {
        sum[z] = applyWindow(psf[z], z, wx, wz, zSelector.background);
    }
    // Smooth the intensity
    ImageJUtils.showStatus("Normalising PSF");
    final Smoother smoother = zSelector.ssmoother;
    final double[] ssum = smoother.smooth(sum).getDSmooth();
    // Compute normalisation and apply.
    SimpleArrayUtils.multiply(ssum, 1.0 / MathUtils.max(ssum));
    for (int z = 0; z < psf.length; z++) {
        if (sum[z] != 0) {
            SimpleArrayUtils.multiply(psf[z], ssum[z] / sum[z]);
        }
        sum[z] = MathUtils.sum(psf[z]);
    }
    // Show the final intensity profile
    final double[] slice = SimpleArrayUtils.newArray(maxz, 1, 1.0);
    final Plot plot = new Plot(TITLE_SIGNAL, "Slice", "Signal");
    final double[] range = MathUtils.limits(sum);
    plot.setLimits(1, maxz, range[0], range[1]);
    plot.setColor(Color.black);
    plot.addPoints(slice, sum, Plot.LINE);
    ImageJUtils.display(TITLE_SIGNAL, plot);
    // Create a new extracted PSF and show
    ImageJUtils.showStatus("Displaying PSF");
    magnification = finalPsf.magnification;
    finalPsf = new ExtractedPsf(psf, finalPsf.maxx, finalPsf.centre, magnification);
    finalPsf.createProjections();
    psfOut = finalPsf.show(TITLE_PSF, zCentre);
    psfImp = psfOut[0];
    // Add image info
    final int imageCount = centres.length;
    final ImagePSF.Builder imagePsf = ImagePsfHelper.create(zCentre, nmPerPixel / magnification, settings.getNmPerSlice() / magnification, imageCount, 0, createNote()).toBuilder();
    // Add the CoM
    // Find the XY centre around the z centre
    final double[] com = getCentreOfMassXy(finalPsf.psf, finalPsf.maxx, finalPsf.maxy, zCentre - 1, settings.getComWindow(), getComXyBorder(finalPsf.maxx, finalPsf.maxy));
    imagePsf.setXCentre(com[0]);
    imagePsf.setYCentre(com[1]);
    imagePsf.setZCentre(zCentre - 1);
    psfImp.setProperty("Info", ImagePsfHelper.toString(imagePsf));
    psfImp.setRoi(new OffsetPointRoi(com[0], com[1]));
    psfImp.setSlice(zCentre);
    psfImp.resetDisplayRange();
    psfImp.updateAndDraw();
    ImageJUtils.log("Final Centre-of-mass = %s,%s\n", rounder.toString(com[0]), rounder.toString(com[1]));
    ImageJUtils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images\n", psfImp.getTitle(), zCentre, MathUtils.rounded(imagePsf.getPixelSize(), 3), MathUtils.rounded(imagePsf.getPixelDepth(), 3), imageCount);
    if (settings.getOutputType() == OUTPUT_TYPE_CSPLINE) {
        // Ask this again as it is important
        // if (TextUtils.isNullOrEmpty(settings.getSplineFilename()))
        // {
        final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
        gd.addFilenameField("Spline_filename", settings.getSplineFilename());
        gd.showDialog(true);
        if (gd.wasCanceled()) {
            return;
        }
        settings.setSplineFilename(gd.getNextString());
        // }
        if (!TextUtils.isNullOrEmpty(settings.getSplineFilename())) {
            // Save the result ...
            IJ.showStatus("Creating cubic spline");
            final CubicSplinePsf cubicSplinePsf = CubicSplineManager.createCubicSpline(imagePsf, psfImp.getImageStack(), settings.getSinglePrecision());
            IJ.showStatus("Saving cubic spline");
            CubicSplineManager.save(cubicSplinePsf, settings.getSplineFilename());
            final String msg = "Spline saved to " + settings.getSplineFilename();
            IJ.showStatus(msg);
            IJ.log(msg);
            // To leave the status message
            return;
        }
    }
    IJ.showStatus("");
}
Also used : CameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CameraModel) CcdCameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CcdCameraModel) CcdCameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CcdCameraModel) BasePoint(uk.ac.sussex.gdsc.core.match.BasePoint) Rectangle(java.awt.Rectangle) CubicSplinePsf(uk.ac.sussex.gdsc.smlm.ij.plugins.CubicSplineManager.CubicSplinePsf) ImageCanvas(ij.gui.ImageCanvas) ImagePSF(uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.ImagePSF) GenericDialog(ij.gui.GenericDialog) NonBlockingExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) Overlay(ij.gui.Overlay) OffsetPointRoi(uk.ac.sussex.gdsc.core.ij.gui.OffsetPointRoi) Plot(ij.gui.Plot) NonBlockingExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog) Point(java.awt.Point) BasePoint(uk.ac.sussex.gdsc.core.match.BasePoint) NonBlockingExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) CalibrationReader(uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader) OffsetPointRoi(uk.ac.sussex.gdsc.core.ij.gui.OffsetPointRoi) Roi(ij.gui.Roi) Point(java.awt.Point) BasePoint(uk.ac.sussex.gdsc.core.match.BasePoint)

Example 9 with CameraModel

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

the class PeakFit method run.

/**
 * Locate the peaks in the configured image source. Results are saved to the configured output.
 *
 * <p>This must be called after initialisation with an image source. Note that each call to this
 * method must be preceded with initialisation to prepare the image and output options.
 */
@SuppressWarnings("null")
public void run() {
    if (source == null) {
        return;
    }
    final int totalFrames = source.getFrames();
    final ImageStack stack = (extraSettings.showProcessedFrames) ? new ImageStack(bounds.width, bounds.height) : null;
    // Do not crop the region from the source if the bounds match the source dimensions
    final Rectangle cropBounds = (bounds.x == 0 && bounds.y == 0 && bounds.width == source.getWidth() && bounds.height == source.getHeight()) ? null : bounds;
    // Use the FitEngine to allow multi-threading.
    final FitEngine engine = createFitEngine(getNumberOfThreads(totalFrames));
    if (engine == null) {
        return;
    }
    final int step = ImageJUtils.getProgressInterval(totalFrames);
    // To pre-process data for noise estimation
    boolean isFitCameraCounts = false;
    CameraModel cameraModel = null;
    if (ignoreBoundsForNoise) {
        isFitCameraCounts = fitConfig.isFitCameraCounts();
        cameraModel = fitConfig.getCameraModel();
    }
    runTime = System.nanoTime();
    boolean shutdown = false;
    int slice = 0;
    final String format = String.format("Slice: %%d / %d (Results=%%d)", totalFrames);
    while (!shutdown) {
        // Noise can optionally be estimated from the entire frame
        float[] data = (ignoreBoundsForNoise) ? source.next() : source.next(cropBounds);
        if (data == null) {
            break;
        }
        if (++slice % step == 0) {
            final int frames = slice;
            if (ImageJUtils.showStatus(() -> String.format(format, frames, results.size()))) {
                IJ.showProgress(slice, totalFrames);
            }
        }
        float noise = Float.NaN;
        if (ignoreBoundsForNoise) {
            // We must pre-process the data before noise estimation
            final float[] data2 = data.clone();
            if (isFitCameraCounts) {
                cameraModel.removeBias(data2);
            } else {
                cameraModel.removeBiasAndGain(data2);
            }
            noise = FitWorker.estimateNoise(data2, source.getWidth(), source.getHeight(), config.getNoiseMethod());
            // Crop the data to the region
            data = ImageJImageConverter.getData(data, source.getWidth(), source.getHeight(), bounds, null);
        }
        if (stack != null) {
            stack.addSlice(String.format("Frame %d - %d", source.getStartFrameNumber(), source.getEndFrameNumber()), data);
        }
        // Get the frame number from the source to allow for interlaced and aggregated data
        engine.run(createJob(source.getStartFrameNumber(), source.getEndFrameNumber(), data, bounds, noise));
        shutdown = escapePressed();
    }
    engine.end(shutdown);
    time = engine.getTime();
    runTime = System.nanoTime() - runTime;
    if (stack != null) {
        ImageJUtils.display("Processed frames", stack);
    }
    showResults();
    source.close();
}
Also used : FitEngine(uk.ac.sussex.gdsc.smlm.engine.FitEngine) CameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CameraModel) PerPixelCameraModel(uk.ac.sussex.gdsc.smlm.model.camera.PerPixelCameraModel) ImageStack(ij.ImageStack) Rectangle(java.awt.Rectangle)

Aggregations

CameraModel (uk.ac.sussex.gdsc.smlm.model.camera.CameraModel)9 Rectangle (java.awt.Rectangle)7 ImageStack (ij.ImageStack)4 PerPixelCameraModel (uk.ac.sussex.gdsc.smlm.model.camera.PerPixelCameraModel)4 GenericDialog (ij.gui.GenericDialog)3 Overlay (ij.gui.Overlay)3 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)3 CalibrationReader (uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader)3 ImagePlus (ij.ImagePlus)2 Plot (ij.gui.Plot)2 Roi (ij.gui.Roi)2 LUT (ij.process.LUT)2 Color (java.awt.Color)2 IJImageSource (uk.ac.sussex.gdsc.smlm.ij.IJImageSource)2 ImageSource (uk.ac.sussex.gdsc.smlm.results.ImageSource)2 IJ (ij.IJ)1 Prefs (ij.Prefs)1 WindowManager (ij.WindowManager)1 ImageCanvas (ij.gui.ImageCanvas)1 ImageRoi (ij.gui.ImageRoi)1