gdsc.smlm.ij.plugins.FilterResults.java Source code

Java tutorial

Introduction

Here is the source code for gdsc.smlm.ij.plugins.FilterResults.java

Source

package gdsc.smlm.ij.plugins;

/*----------------------------------------------------------------------------- 
 * GDSC SMLM Software
 * 
 * Copyright (C) 2013 Alex Herbert
 * Genome Damage and Stability Centre
 * University of Sussex, UK
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *---------------------------------------------------------------------------*/

import gdsc.smlm.function.gaussian.Gaussian2DFunction;
import gdsc.smlm.ij.plugins.ResultsManager.InputSource;
import gdsc.smlm.ij.settings.FilterSettings;
import gdsc.smlm.ij.settings.GlobalSettings;
import gdsc.smlm.ij.settings.SettingsManager;
import gdsc.smlm.model.MaskDistribution;
import gdsc.smlm.results.MemoryPeakResults;
import gdsc.smlm.results.PeakResult;
import ij.IJ;
import ij.ImagePlus;
import ij.WindowManager;
import ij.gui.GenericDialog;
import ij.plugin.PlugIn;
import ij.process.ByteProcessor;

import java.util.ArrayList;

import org.apache.commons.math3.util.FastMath;

/**
 * Filters PeakFit results that are stored in memory using various fit criteria.
 */
public class FilterResults implements PlugIn {
    private static final String TITLE = "Filter Results";
    private static String inputOption = "";

    private MemoryPeakResults results;

    private FilterSettings filterSettings = new FilterSettings();

    private double nmPerPixel, gain;
    private boolean emCCD;

    // Used to pass data from analyseResults() to checkLimits()
    private float minDrift = Float.MAX_VALUE;
    private float maxDrift = 0;
    private float minSignal = Float.MAX_VALUE;
    private float maxSignal = 0;
    private float minSNR = Float.MAX_VALUE;
    private float maxSNR = 0;
    private double minPrecision = Float.MAX_VALUE;
    private double maxPrecision = 0;
    private double averageWidth = 0;
    private float minWidth = Float.MAX_VALUE;
    private float maxWidth = 0;

    /*
     * (non-Javadoc)
     * 
     * @see ij.plugin.PlugIn#run(java.lang.String)
     */
    public void run(String arg) {
        if (MemoryPeakResults.countMemorySize() == 0) {
            IJ.error(TITLE, "There are no fitting results in memory");
            return;
        }

        // Show a dialog allowing the results set to be filtered
        GenericDialog gd = new GenericDialog(TITLE);
        gd.addMessage("Select a dataset to filter");
        ResultsManager.addInput(gd, inputOption, InputSource.MEMORY);
        gd.showDialog();
        if (gd.wasCanceled())
            return;
        inputOption = ResultsManager.getInputSource(gd);
        results = ResultsManager.loadInputResults(inputOption, false);
        if (results == null || results.size() == 0) {
            IJ.error(TITLE, "No results could be loaded");
            IJ.showStatus("");
            return;
        }

        analyseResults();

        if (!showDialog())
            return;

        filterResults();
    }

    /**
     * Analyse the results and determine the range for each filter
     */
    private void analyseResults() {
        IJ.showStatus("Analysing results ...");

        nmPerPixel = results.getNmPerPixel();
        gain = results.getGain();
        emCCD = results.isEMCCD();

        double maxVariance = maxPrecision * maxPrecision;
        double minVariance = minPrecision * minPrecision;

        int size = results.size();
        int i = 0;
        for (PeakResult result : results.getResults()) {
            if (i % 64 == 0)
                IJ.showProgress(i, size);

            final float drift = getDrift(result);
            if (maxDrift < drift)
                maxDrift = drift;
            if (minDrift > drift)
                minDrift = drift;

            final float signal = result.getSignal();
            if (maxSignal < signal)
                maxSignal = signal;
            if (minSignal > signal)
                minSignal = signal;

            final float snr = getSNR(result);
            if (maxSNR < snr)
                maxSNR = snr;
            if (minSNR > snr)
                minSNR = snr;

            // Use variance to avoid sqrt()
            final double variance = getVariance(result);
            if (maxVariance < variance)
                maxVariance = variance;
            if (minVariance > variance)
                minVariance = variance;

            final float width = getWidth(result);
            averageWidth += width;
            if (maxWidth < width)
                maxWidth = width;
            if (minWidth > width)
                minWidth = width;
        }
        averageWidth /= results.size();

        maxPrecision = Math.sqrt(maxVariance);
        minPrecision = Math.sqrt(minVariance);

        IJ.showProgress(1);
        IJ.showStatus("");
    }

    private float getDrift(PeakResult result) {
        float drift = FastMath.max(Math.abs(result.origX - result.params[Gaussian2DFunction.X_POSITION]),
                Math.abs(result.origY - result.params[Gaussian2DFunction.Y_POSITION]));
        return drift;
    }

    private float getSNR(PeakResult result) {
        if (result.noise <= 0)
            return 0;
        return result.getSignal() / result.noise;
    }

    private double getVariance(PeakResult result) {
        return PeakResult.getVariance(nmPerPixel, result.getSD() * nmPerPixel, result.getSignal() / gain,
                result.noise / gain, emCCD);
    }

    private float getWidth(PeakResult result) {
        // The X-width should be the largest (major axis)
        // Q. Should a filter be used for the Y-width too?
        return result.params[Gaussian2DFunction.X_SD];
    }

    /**
     * Check that none of the filter values are outside the limits
     */
    private void checkLimits() {
        if (filterSettings.maxDrift > maxDrift || filterSettings.maxDrift < minDrift)
            filterSettings.maxDrift = maxDrift;

        if (filterSettings.minSignal > maxSignal || filterSettings.minSignal < minSignal)
            filterSettings.minSignal = minSignal;

        if (filterSettings.minSNR > maxSNR || filterSettings.minSNR < minSNR)
            filterSettings.minSNR = minSNR;

        if (filterSettings.maxPrecision > maxPrecision || filterSettings.maxPrecision < minPrecision)
            filterSettings.maxPrecision = maxPrecision;

        if (filterSettings.minWidth > maxWidth || filterSettings.minWidth < minWidth)
            filterSettings.minWidth = minWidth;

        if (filterSettings.maxWidth > maxWidth || filterSettings.maxWidth < minWidth)
            filterSettings.maxWidth = maxWidth;

        if (filterSettings.minWidth > filterSettings.maxWidth) {
            float tmp = filterSettings.maxWidth;
            filterSettings.maxWidth = filterSettings.minWidth;
            filterSettings.minWidth = tmp;
        }
    }

    /**
     * Apply the filters to the data
     */
    private void filterResults() {
        checkLimits();

        MemoryPeakResults newResults = new MemoryPeakResults();
        newResults.copySettings(results);
        newResults.setName(results.getName() + " Filtered");

        // Initialise the mask
        ByteProcessor mask = getMask(filterSettings.maskTitle);
        MaskDistribution maskFilter = null;
        final float centreX = results.getBounds().width / 2.0f;
        final float centreY = results.getBounds().height / 2.0f;
        if (mask != null) {
            double scaleX = (double) results.getBounds().width / mask.getWidth();
            double scaleY = (double) results.getBounds().height / mask.getHeight();
            maskFilter = new MaskDistribution((byte[]) mask.getPixels(), mask.getWidth(), mask.getHeight(), 0,
                    scaleX, scaleY);
        }

        int i = 0;
        final int size = results.size();
        final double maxVariance = filterSettings.maxPrecision * filterSettings.maxPrecision;
        for (PeakResult result : results.getResults()) {
            if (i % 64 == 0)
                IJ.showProgress(i, size);

            if (getDrift(result) > filterSettings.maxDrift)
                continue;

            if (result.getSignal() < filterSettings.minSignal)
                continue;

            if (getSNR(result) < filterSettings.minSNR)
                continue;

            if (getVariance(result) > maxVariance)
                continue;

            final float width = getWidth(result);
            if (width < filterSettings.minWidth || width > filterSettings.maxWidth)
                continue;

            if (maskFilter != null) {
                // Check the coordinates are inside the mask
                double[] xy = new double[] { result.getXPosition() - centreX, result.getYPosition() - centreY };
                if (!maskFilter.isWithinXY(xy))
                    continue;
            }

            // Passed all filters. Add to the results
            newResults.add(result);
        }

        IJ.showProgress(1);
        IJ.showStatus(newResults.size() + " Filtered localisations");
        MemoryPeakResults.addResults(newResults);
    }

    private ByteProcessor getMask(String maskTitle) {
        ImagePlus imp = WindowManager.getImage(maskTitle);
        if (imp != null) {
            return (ByteProcessor) imp.getProcessor().convertToByte(false);
        }
        return null;
    }

    private boolean showDialog() {
        GenericDialog gd = new GenericDialog(TITLE);
        gd.addHelp(About.HELP_URL);

        GlobalSettings gs = SettingsManager.loadSettings();
        filterSettings = gs.getFilterSettings();

        checkLimits();

        gd.addSlider("Max_drift", minDrift, maxDrift, filterSettings.maxDrift);
        gd.addSlider("Min_Signal", minSignal, maxSignal, filterSettings.minSignal);
        gd.addSlider("Min_SNR", minSNR, maxSNR, filterSettings.minSNR);
        gd.addSlider("Min_Precision", minPrecision, maxPrecision, filterSettings.maxPrecision);

        // TODO - If calibrated present the widths in nm
        gd.addMessage("Average Width = " + IJ.d2s(averageWidth, 3));

        gd.addSlider("Min_Width", minWidth, maxWidth, filterSettings.minWidth);
        gd.addSlider("Max_Width", minWidth, maxWidth, filterSettings.maxWidth);

        // Get a list of potential mask images
        String[] items = getImageList();
        gd.addChoice("Mask", items, filterSettings.maskTitle);

        gd.showDialog();

        if (gd.wasCanceled())
            return false;

        filterSettings.maxDrift = (float) gd.getNextNumber();
        filterSettings.minSignal = (float) gd.getNextNumber();
        filterSettings.minSNR = (float) gd.getNextNumber();
        filterSettings.maxPrecision = (float) gd.getNextNumber();
        filterSettings.minWidth = (float) gd.getNextNumber();
        filterSettings.maxWidth = (float) gd.getNextNumber();
        filterSettings.maskTitle = gd.getNextChoice();

        return SettingsManager.saveSettings(gs);
    }

    /**
     * Build a list of all the image names.
     * 
     * @return The list of images
     */
    public static String[] getImageList() {
        ArrayList<String> newImageList = new ArrayList<String>();
        newImageList.add("[None]");

        for (int id : getIDList()) {
            ImagePlus imp = WindowManager.getImage(id);
            if (imp == null)
                continue;
            if (!imp.getProcessor().isBinary())
                continue;
            newImageList.add(imp.getTitle());
        }

        return newImageList.toArray(new String[0]);
    }

    public static int[] getIDList() {
        int[] list = WindowManager.getIDList();
        return (list != null) ? list : new int[0];
    }
}