Search in sources :

Example 1 with SeriesOpener

use of gdsc.smlm.ij.utils.SeriesOpener in project GDSC-SMLM by aherbert.

the class PeakFit method setup.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.filter.PlugInFilter#setup(java.lang.String, ij.ImagePlus)
	 */
public int setup(String arg, ImagePlus imp) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    plugin_flags = FLAGS;
    extraOptions = Utils.isExtraOptions();
    maximaIdentification = (arg != null && arg.contains("spot"));
    fitMaxima = (arg != null && arg.contains("maxima"));
    simpleFit = (arg != null && arg.contains("simple"));
    boolean runSeries = (arg != null && arg.contains("series"));
    ImageSource imageSource = null;
    if (fitMaxima) {
        imp = null;
        // The image source will be found from the peak results.
        if (!showMaximaDialog())
            return DONE;
        MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
        if (results == null || results.size() == 0) {
            IJ.error(TITLE, "No results could be loaded");
            return DONE;
        }
        // Check for single frame
        singleFrame = results.getHead().getFrame();
        for (PeakResult result : results.getResults()) {
            if (singleFrame != result.getFrame()) {
                singleFrame = 0;
                break;
            }
        }
        imageSource = results.getSource();
        plugin_flags |= NO_IMAGE_REQUIRED;
    } else if (runSeries) {
        imp = null;
        // Select input folder
        String inputDirectory;
        inputDirectory = IJ.getDirectory("Select image series ...");
        //inputDirectory = getInputDirectory("Select image series ...");
        if (inputDirectory == null)
            return DONE;
        // Load input series ...
        SeriesOpener series;
        if (extraOptions)
            series = new SeriesOpener(inputDirectory, true, numberOfThreads);
        else
            series = new SeriesOpener(inputDirectory);
        if (series.getNumberOfImages() == 0) {
            IJ.error(TITLE, "No images in the selected directory:\n" + inputDirectory);
            return DONE;
        }
        SeriesImageSource seriesImageSource = new SeriesImageSource(getName(series.getImageList()), series);
        seriesImageSource.setLogProgress(true);
        if (extraOptions) {
            numberOfThreads = Math.max(1, series.getNumberOfThreads());
            seriesImageSource.setNumberOfThreads(numberOfThreads);
        }
        imageSource = seriesImageSource;
        plugin_flags |= NO_IMAGE_REQUIRED;
    } else {
        if (imp == null) {
            IJ.noImage();
            return DONE;
        }
        // Check it is not a previous result
        if (imp.getTitle().endsWith(IJImagePeakResults.IMAGE_SUFFIX)) {
            IJImageSource tmpImageSource = null;
            // Check the image to see if it has an image source XML structure in the info property
            Object o = imp.getProperty("Info");
            Pattern pattern = Pattern.compile("Source: (<.*IJImageSource>.*<.*IJImageSource>)", Pattern.DOTALL);
            Matcher match = pattern.matcher((o == null) ? "" : o.toString());
            if (match.find()) {
                ImageSource source = ImageSource.fromXML(match.group(1));
                if (source instanceof IJImageSource) {
                    tmpImageSource = (IJImageSource) source;
                    if (!tmpImageSource.open()) {
                        tmpImageSource = null;
                    } else {
                        imp = WindowManager.getImage(tmpImageSource.getName());
                    }
                }
            }
            if (tmpImageSource == null) {
                // Look for a parent using the title
                String parentTitle = imp.getTitle().substring(0, imp.getTitle().length() - IJImagePeakResults.IMAGE_SUFFIX.length() - 1);
                ImagePlus parentImp = WindowManager.getImage(parentTitle);
                if (parentImp != null) {
                    tmpImageSource = new IJImageSource(parentImp);
                    imp = parentImp;
                }
            }
            String message = "The selected image may be a previous fit result";
            if (tmpImageSource != null) {
                // are missing
                if (!Utils.isNullOrEmpty(tmpImageSource.getName()))
                    message += " of: \n \n" + tmpImageSource.getName();
                message += " \n \nFit the parent?";
            } else
                message += " \n \nDo you want to continue?";
            YesNoCancelDialog d = new YesNoCancelDialog(null, TITLE, message);
            if (tmpImageSource == null) {
                if (!d.yesPressed())
                    return DONE;
            } else {
                if (d.yesPressed())
                    imageSource = tmpImageSource;
                if (d.cancelPressed())
                    return DONE;
            }
        }
        if (imageSource == null)
            imageSource = new IJImageSource(imp);
    }
    time = -1;
    if (!initialiseImage(imageSource, getBounds(imp), false)) {
        IJ.error(TITLE, "Failed to initialise the source image: " + imageSource.getName());
        return DONE;
    }
    int flags = showDialog(imp);
    if ((flags & DONE) == 0) {
        // Repeat so that we pass in the selected option for ignoring the bounds.
        // This should not be necessary since it is set within the readDialog method.
        //if (ignoreBoundsForNoise)
        //	initialiseImage(imageSource, bounds, ignoreBoundsForNoise);
        initialiseFitting();
    }
    return flags;
}
Also used : Pattern(java.util.regex.Pattern) Matcher(java.util.regex.Matcher) SeriesImageSource(gdsc.smlm.ij.SeriesImageSource) SeriesOpener(gdsc.smlm.ij.utils.SeriesOpener) ImagePlus(ij.ImagePlus) PeakResult(gdsc.smlm.results.PeakResult) ExtendedPeakResult(gdsc.smlm.results.ExtendedPeakResult) IJImageSource(gdsc.smlm.ij.IJImageSource) MemoryPeakResults(gdsc.smlm.results.MemoryPeakResults) YesNoCancelDialog(ij.gui.YesNoCancelDialog) AggregatedImageSource(gdsc.smlm.results.AggregatedImageSource) IJImageSource(gdsc.smlm.ij.IJImageSource) SeriesImageSource(gdsc.smlm.ij.SeriesImageSource) InterlacedImageSource(gdsc.smlm.results.InterlacedImageSource) ImageSource(gdsc.smlm.results.ImageSource)

Example 2 with SeriesOpener

use of gdsc.smlm.ij.utils.SeriesOpener in project GDSC-SMLM by aherbert.

the class MeanVarianceTest method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (Utils.isExtraOptions()) {
        ImagePlus imp = WindowManager.getCurrentImage();
        if (imp.getStackSize() > 1) {
            GenericDialog gd = new GenericDialog(TITLE);
            gd.addMessage("Perform single image analysis on the current image?");
            gd.addNumericField("Bias", _bias, 0);
            gd.showDialog();
            if (gd.wasCanceled())
                return;
            singleImage = true;
            _bias = Math.abs(gd.getNextNumber());
        } else {
            IJ.error(TITLE, "Single-image mode requires a stack");
            return;
        }
    }
    List<ImageSample> images;
    String inputDirectory = "";
    if (singleImage) {
        IJ.showStatus("Loading images...");
        images = getImages();
        if (images.size() == 0) {
            IJ.error(TITLE, "Not enough images for analysis");
            return;
        }
    } else {
        inputDirectory = IJ.getDirectory("Select image series ...");
        if (inputDirectory == null)
            return;
        SeriesOpener series = new SeriesOpener(inputDirectory, false, 0);
        series.setVariableSize(true);
        if (series.getNumberOfImages() < 3) {
            IJ.error(TITLE, "Not enough images in the selected directory");
            return;
        }
        if (!IJ.showMessageWithCancel(TITLE, String.format("Analyse %d images, first image:\n%s", series.getNumberOfImages(), series.getImageList()[0]))) {
            return;
        }
        IJ.showStatus("Loading images");
        images = getImages(series);
        if (images.size() < 3) {
            IJ.error(TITLE, "Not enough images for analysis");
            return;
        }
        if (images.get(0).exposure != 0) {
            IJ.error(TITLE, "First image in series must have exposure 0 (Bias image)");
            return;
        }
    }
    boolean emMode = (arg != null && arg.contains("em"));
    GenericDialog gd = new GenericDialog(TITLE);
    gd.addMessage("Set the output options:");
    gd.addCheckbox("Show_table", showTable);
    gd.addCheckbox("Show_charts", showCharts);
    if (emMode) {
        // Ask the user for the camera gain ...
        gd.addMessage("Estimating the EM-gain requires the camera gain without EM readout enabled");
        gd.addNumericField("Camera_gain (ADU/e-)", cameraGain, 4);
    }
    gd.showDialog();
    if (gd.wasCanceled())
        return;
    showTable = gd.getNextBoolean();
    showCharts = gd.getNextBoolean();
    if (emMode) {
        cameraGain = gd.getNextNumber();
    }
    IJ.showStatus("Computing mean & variance");
    final double nImages = images.size();
    for (int i = 0; i < images.size(); i++) {
        IJ.showStatus(String.format("Computing mean & variance %d/%d", i + 1, images.size()));
        images.get(i).compute(singleImage, i / nImages, (i + 1) / nImages);
    }
    IJ.showProgress(1);
    IJ.showStatus("Computing results");
    // Allow user to input multiple bias images
    int start = 0;
    Statistics biasStats = new Statistics();
    Statistics noiseStats = new Statistics();
    final double bias;
    if (singleImage) {
        bias = _bias;
    } else {
        while (start < images.size()) {
            ImageSample sample = images.get(start);
            if (sample.exposure == 0) {
                biasStats.add(sample.means);
                for (PairSample pair : sample.samples) {
                    noiseStats.add(pair.variance);
                }
                start++;
            } else
                break;
        }
        bias = biasStats.getMean();
    }
    // Get the mean-variance data
    int total = 0;
    for (int i = start; i < images.size(); i++) total += images.get(i).samples.size();
    if (showTable && total > 2000) {
        gd = new GenericDialog(TITLE);
        gd.addMessage("Table output requires " + total + " entries.\n \nYou may want to disable the table.");
        gd.addCheckbox("Show_table", showTable);
        gd.showDialog();
        if (gd.wasCanceled())
            return;
        showTable = gd.getNextBoolean();
    }
    TextWindow results = (showTable) ? createResultsWindow() : null;
    double[] mean = new double[total];
    double[] variance = new double[mean.length];
    Statistics gainStats = (singleImage) ? new StoredDataStatistics(total) : new Statistics();
    final WeightedObservedPoints obs = new WeightedObservedPoints();
    for (int i = (singleImage) ? 0 : start, j = 0; i < images.size(); i++) {
        StringBuilder sb = (showTable) ? new StringBuilder() : null;
        ImageSample sample = images.get(i);
        for (PairSample pair : sample.samples) {
            if (j % 16 == 0)
                IJ.showProgress(j, total);
            mean[j] = pair.getMean();
            variance[j] = pair.variance;
            // Gain is in ADU / e
            double gain = variance[j] / (mean[j] - bias);
            gainStats.add(gain);
            obs.add(mean[j], variance[j]);
            if (emMode) {
                gain /= (2 * cameraGain);
            }
            if (showTable) {
                sb.append(sample.title).append("\t");
                sb.append(sample.exposure).append("\t");
                sb.append(pair.slice1).append("\t");
                sb.append(pair.slice2).append("\t");
                sb.append(IJ.d2s(pair.mean1, 2)).append("\t");
                sb.append(IJ.d2s(pair.mean2, 2)).append("\t");
                sb.append(IJ.d2s(mean[j], 2)).append("\t");
                sb.append(IJ.d2s(variance[j], 2)).append("\t");
                sb.append(Utils.rounded(gain, 4)).append("\n");
            }
            j++;
        }
        if (showTable)
            results.append(sb.toString());
    }
    IJ.showProgress(1);
    if (singleImage) {
        StoredDataStatistics stats = (StoredDataStatistics) gainStats;
        Utils.log(TITLE);
        if (emMode) {
            double[] values = stats.getValues();
            MathArrays.scaleInPlace(0.5, values);
            stats = new StoredDataStatistics(values);
        }
        if (showCharts) {
            // Plot the gain over time
            String title = TITLE + " Gain vs Frame";
            Plot2 plot = new Plot2(title, "Slice", "Gain", Utils.newArray(gainStats.getN(), 1, 1.0), stats.getValues());
            PlotWindow pw = Utils.display(title, plot);
            // Show a histogram
            String label = String.format("Mean = %s, Median = %s", Utils.rounded(stats.getMean()), Utils.rounded(stats.getMedian()));
            int id = Utils.showHistogram(TITLE, stats, "Gain", 0, 1, 100, true, label);
            if (Utils.isNewWindow()) {
                Point point = pw.getLocation();
                point.x = pw.getLocation().x;
                point.y += pw.getHeight();
                WindowManager.getImage(id).getWindow().setLocation(point);
            }
        }
        Utils.log("Single-image mode: %s camera", (emMode) ? "EM-CCD" : "Standard");
        final double gain = stats.getMedian();
        if (emMode) {
            final double totalGain = gain;
            final double emGain = totalGain / cameraGain;
            Utils.log("  Gain = 1 / %s (ADU/e-)", Utils.rounded(cameraGain, 4));
            Utils.log("  EM-Gain = %s", Utils.rounded(emGain, 4));
            Utils.log("  Total Gain = %s (ADU/e-)", Utils.rounded(totalGain, 4));
        } else {
            cameraGain = gain;
            Utils.log("  Gain = 1 / %s (ADU/e-)", Utils.rounded(cameraGain, 4));
        }
    } else {
        IJ.showStatus("Computing fit");
        // Sort
        int[] indices = rank(mean);
        mean = reorder(mean, indices);
        variance = reorder(variance, indices);
        // Compute optimal coefficients.
        // a - b x
        final double[] init = { 0, 1 / gainStats.getMean() };
        final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(2).withStartPoint(init);
        final double[] best = fitter.fit(obs.toList());
        // Construct the polynomial that best fits the data.
        final PolynomialFunction fitted = new PolynomialFunction(best);
        if (showCharts) {
            // Plot mean verses variance. Gradient is gain in ADU/e.
            String title = TITLE + " results";
            Plot2 plot = new Plot2(title, "Mean", "Variance");
            double[] xlimits = Maths.limits(mean);
            double[] ylimits = Maths.limits(variance);
            double xrange = (xlimits[1] - xlimits[0]) * 0.05;
            if (xrange == 0)
                xrange = 0.05;
            double yrange = (ylimits[1] - ylimits[0]) * 0.05;
            if (yrange == 0)
                yrange = 0.05;
            plot.setLimits(xlimits[0] - xrange, xlimits[1] + xrange, ylimits[0] - yrange, ylimits[1] + yrange);
            plot.setColor(Color.blue);
            plot.addPoints(mean, variance, Plot2.CROSS);
            plot.setColor(Color.red);
            plot.addPoints(new double[] { mean[0], mean[mean.length - 1] }, new double[] { fitted.value(mean[0]), fitted.value(mean[mean.length - 1]) }, Plot2.LINE);
            Utils.display(title, plot);
        }
        final double avBiasNoise = Math.sqrt(noiseStats.getMean());
        Utils.log(TITLE);
        Utils.log("  Directory = %s", inputDirectory);
        Utils.log("  Bias = %s +/- %s (ADU)", Utils.rounded(bias, 4), Utils.rounded(avBiasNoise, 4));
        Utils.log("  Variance = %s + %s * mean", Utils.rounded(best[0], 4), Utils.rounded(best[1], 4));
        if (emMode) {
            final double emGain = best[1] / (2 * cameraGain);
            // Noise is standard deviation of the bias image divided by the total gain (in ADU/e-)
            final double totalGain = emGain * cameraGain;
            Utils.log("  Read Noise = %s (e-) [%s (ADU)]", Utils.rounded(avBiasNoise / totalGain, 4), Utils.rounded(avBiasNoise, 4));
            Utils.log("  Gain = 1 / %s (ADU/e-)", Utils.rounded(1 / cameraGain, 4));
            Utils.log("  EM-Gain = %s", Utils.rounded(emGain, 4));
            Utils.log("  Total Gain = %s (ADU/e-)", Utils.rounded(totalGain, 4));
        } else {
            // Noise is standard deviation of the bias image divided by the gain (in ADU/e-)
            cameraGain = best[1];
            final double readNoise = avBiasNoise / cameraGain;
            Utils.log("  Read Noise = %s (e-) [%s (ADU)]", Utils.rounded(readNoise, 4), Utils.rounded(readNoise * cameraGain, 4));
            Utils.log("  Gain = 1 / %s (ADU/e-)", Utils.rounded(1 / cameraGain, 4));
        }
    }
    IJ.showStatus("");
}
Also used : StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) PlotWindow(ij.gui.PlotWindow) PolynomialFunction(org.apache.commons.math3.analysis.polynomials.PolynomialFunction) SeriesOpener(gdsc.smlm.ij.utils.SeriesOpener) Plot2(ij.gui.Plot2) Point(java.awt.Point) ImagePlus(ij.ImagePlus) StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) Statistics(gdsc.core.utils.Statistics) Point(java.awt.Point) PolynomialCurveFitter(org.apache.commons.math3.fitting.PolynomialCurveFitter) WeightedObservedPoints(org.apache.commons.math3.fitting.WeightedObservedPoints) TextWindow(ij.text.TextWindow) GenericDialog(ij.gui.GenericDialog)

Aggregations

SeriesOpener (gdsc.smlm.ij.utils.SeriesOpener)2 ImagePlus (ij.ImagePlus)2 Statistics (gdsc.core.utils.Statistics)1 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)1 IJImageSource (gdsc.smlm.ij.IJImageSource)1 SeriesImageSource (gdsc.smlm.ij.SeriesImageSource)1 AggregatedImageSource (gdsc.smlm.results.AggregatedImageSource)1 ExtendedPeakResult (gdsc.smlm.results.ExtendedPeakResult)1 ImageSource (gdsc.smlm.results.ImageSource)1 InterlacedImageSource (gdsc.smlm.results.InterlacedImageSource)1 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)1 PeakResult (gdsc.smlm.results.PeakResult)1 GenericDialog (ij.gui.GenericDialog)1 Plot2 (ij.gui.Plot2)1 PlotWindow (ij.gui.PlotWindow)1 YesNoCancelDialog (ij.gui.YesNoCancelDialog)1 TextWindow (ij.text.TextWindow)1 Point (java.awt.Point)1 Matcher (java.util.regex.Matcher)1 Pattern (java.util.regex.Pattern)1