Examples with Frequency org.apache.commons.math3.stat.Frequency used on opensource projects

Example 16 with Frequency

use of org.apache.commons.math3.stat.Frequency in project ma-modules-public by infiniteautomation.

the class PointValueFftCalculator method streamData.

/* (non-Javadoc)
	 * @see com.serotonin.m2m2.web.mvc.rest.v1.model.pointValue.PointValueTimeStream#streamData(java.io.Writer)
	 */
@Override
public void streamData(JsonGenerator jgen) {
    this.setupDates();
    DateTime startTime = new DateTime(from);
    DateTime endTime = new DateTime(to);
    FftGenerator generator = this.calculate(startTime, endTime);
    double[] fftData = generator.getValues();
    double sampleRateHz = 1000d / generator.getAverageSamplePeriodMs();
    double dataLength = (double) fftData.length;
    // Output The Real Steady State Mangitude
    try {
        jgen.writeStartObject();
        // Amplitude
        jgen.writeNumberField("value", fftData[0]);
        double frequency = 0;
        jgen.writeNumberField("frequency", frequency);
        jgen.writeEndObject();
    } catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
    double realComponent, imaginaryComponent, frequency;
    if (fftData.length % 2 == 0) {
        for (int i = 2; i < fftData.length / 2; i++) {
            try {
                realComponent = fftData[i * 2];
                imaginaryComponent = fftData[2 * i + 1];
                Complex c = new Complex(realComponent, imaginaryComponent);
                jgen.writeStartObject();
                // Amplitude
                jgen.writeNumberField("value", c.abs());
                if (this.returnFrequency)
                    frequency = (double) i * sampleRateHz / dataLength;
                else
                    frequency = 1d / ((double) i * sampleRateHz / dataLength);
                jgen.writeNumberField("frequency", frequency);
                jgen.writeEndObject();
            } catch (Exception e) {
                // TODO Auto-generated catch block
                e.printStackTrace();
            }
        }
    } else {
        for (int i = 2; i < (fftData.length - 1) / 2; i++) {
            try {
                realComponent = fftData[i * 2];
                imaginaryComponent = fftData[2 * i + 1];
                Complex c = new Complex(realComponent, imaginaryComponent);
                jgen.writeStartObject();
                // Amplitude
                jgen.writeNumberField("value", c.abs());
                if (this.returnFrequency)
                    frequency = (double) i * sampleRateHz / dataLength;
                else
                    frequency = 1d / ((double) i * sampleRateHz / dataLength);
                jgen.writeNumberField("frequency", frequency);
                jgen.writeEndObject();
            } catch (Exception e) {
                // TODO Auto-generated catch block
                e.printStackTrace();
            }
        }
        // Write the last value out as it isn't in order in the array
        try {
            realComponent = fftData[fftData.length / 2];
            imaginaryComponent = fftData[1];
            Complex c = new Complex(realComponent, imaginaryComponent);
            jgen.writeStartObject();
            // Amplitude
            jgen.writeNumberField("value", c.abs());
            if (this.returnFrequency)
                frequency = (double) (((fftData.length - 1) / 2) - 1) * sampleRateHz / dataLength;
            else
                frequency = 1d / ((double) (((fftData.length - 1) / 2) - 1) * sampleRateHz / dataLength);
            jgen.writeNumberField("frequency", frequency);
            jgen.writeEndObject();
        } catch (Exception e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
    }
}

Example 17 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    settings = Settings.load();
    // Saved by reference so just save now
    settings.save();
    // Read in multiple traced datasets
    // All datasets must have the same pixel pitch and exposure time
    // Get parameters
    // Convert datasets to tracks
    // For each track compute the 4 local track features using the configured window
    // 
    // Optional:
    // Fit a multi-variate Gaussian mixture model to the data
    // (using the configured number of components/populations)
    // Assign each point in the track using the model.
    // Smooth the assignments.
    // 
    // The alternative is to use the localisation category to assign populations.
    // 
    // Plot histograms of each track parameter, coloured by component
    final List<MemoryPeakResults> combinedResults = new LocalList<>();
    if (!showInputDialog(combinedResults)) {
        return;
    }
    final boolean hasCategory = showHasCategoryDialog(combinedResults);
    if (!showDialog(hasCategory)) {
        return;
    }
    ImageJUtils.log(TITLE + "...");
    final List<Trace> tracks = getTracks(combinedResults, settings.window, settings.minTrackLength);
    if (tracks.isEmpty()) {
        IJ.error(TITLE, "No tracks. Please check the input data and min track length setting.");
        return;
    }
    final Calibration cal = combinedResults.get(0).getCalibration();
    final CalibrationReader cr = new CalibrationReader(cal);
    // Use micrometer / second
    final TypeConverter<DistanceUnit> distanceConverter = cr.getDistanceConverter(DistanceUnit.UM);
    final double exposureTime = cr.getExposureTime() / 1000.0;
    final Pair<int[], double[][]> trackData = extractTrackData(tracks, distanceConverter, exposureTime, hasCategory);
    final double[][] data = trackData.getValue();
    // Histogram the raw data.
    final Array2DRowRealMatrix raw = new Array2DRowRealMatrix(data, false);
    final WindowOrganiser wo = new WindowOrganiser();
    // Store the histogram data for plotting the components
    final double[][] columns = new double[FEATURE_NAMES.length][];
    final double[][] limits = new double[FEATURE_NAMES.length][];
    // Get column data
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        columns[i] = raw.getColumn(i);
        if (i == FEATURE_D) {
            // Plot using a logarithmic scale
            SimpleArrayUtils.apply(columns[i], Math::log10);
        }
        limits[i] = MathUtils.limits(columns[i]);
    }
    // Compute histogram bins
    final int[] bins = new int[FEATURE_NAMES.length];
    if (settings.histogramBins > 0) {
        Arrays.fill(bins, settings.histogramBins);
    } else {
        for (int i = 0; i < FEATURE_NAMES.length; i++) {
            bins[i] = HistogramPlot.getBins(StoredData.create(columns[i]), BinMethod.FD);
        }
        // Use the maximum so all histograms look the same
        Arrays.fill(bins, MathUtils.max(bins));
    }
    // Compute plots
    final Plot[] plots = new Plot[FEATURE_NAMES.length];
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        final double[][] hist = HistogramPlot.calcHistogram(columns[i], limits[i][0], limits[i][1], bins[i]);
        plots[i] = new Plot(TITLE + " " + FEATURE_NAMES[i], getFeatureLabel(i, i == FEATURE_D), "Frequency");
        plots[i].addPoints(hist[0], hist[1], Plot.BAR);
        ImageJUtils.display(plots[i].getTitle(), plots[i], ImageJUtils.NO_TO_FRONT, wo);
    }
    wo.tile();
    // The component for each data point
    int[] component;
    // The number of components
    int numComponents;
    // Data used to fit the Gaussian mixture model
    double[][] fitData;
    // The fitted model
    MixtureMultivariateGaussianDistribution model;
    if (hasCategory) {
        // Use the category as the component.
        // No fit data and no output model
        fitData = null;
        model = null;
        // The component is stored at the end of the raw track data.
        final int end = data[0].length - 1;
        component = Arrays.stream(data).mapToInt(d -> (int) d[end]).toArray();
        numComponents = MathUtils.max(component) + 1;
        // In the EM algorithm the probability of each data point is computed and normalised to
        // sum to 1. The normalised probabilities are averaged to create the weights.
        // Note the probability of each data point uses the previous weight and the algorithm
        // iterates.
        // This is not a fitted model but the input model so use
        // zero weights to indicate no fitting was performed.
        final double[] weights = new double[numComponents];
        // Remove the trailing component to show the 'model' in a table.
        createModelTable(Arrays.stream(data).map(d -> Arrays.copyOf(d, end)).toArray(double[][]::new), weights, component);
    } else {
        // Multivariate Gaussian mixture EM
        // Provide option to not use the anomalous exponent in the population mix.
        int sortDimension = SORT_DIMENSION;
        if (settings.ignoreAlpha) {
            // Remove index 0. This shifts the sort dimension.
            sortDimension--;
            fitData = Arrays.stream(data).map(d -> Arrays.copyOfRange(d, 1, d.length)).toArray(double[][]::new);
        } else {
            fitData = SimpleArrayUtils.deepCopy(data);
        }
        final MultivariateGaussianMixtureExpectationMaximization mixed = fitGaussianMixture(fitData, sortDimension);
        if (mixed == null) {
            IJ.error(TITLE, "Failed to fit a mixture model");
            return;
        }
        model = sortComponents(mixed.getFittedModel(), sortDimension);
        // For the best model, assign to the most likely population.
        component = assignData(fitData, model);
        // Table of the final model using the original data (i.e. not normalised)
        final double[] weights = model.getWeights();
        numComponents = weights.length;
        createModelTable(data, weights, component);
    }
    // Output coloured histograms of the populations.
    final LUT lut = LutHelper.createLut(settings.lutIndex);
    IntFunction<Color> colourMap;
    if (LutHelper.getColour(lut, 0).equals(Color.BLACK)) {
        colourMap = i -> LutHelper.getNonZeroColour(lut, i, 0, numComponents - 1);
    } else {
        colourMap = i -> LutHelper.getColour(lut, i, 0, numComponents - 1);
    }
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        // Extract the data for each component
        final double[] col = columns[i];
        final Plot plot = plots[i];
        for (int n = 0; n < numComponents; n++) {
            final StoredData feature = new StoredData();
            for (int j = 0; j < component.length; j++) {
                if (component[j] == n) {
                    feature.add(col[j]);
                }
            }
            if (feature.size() == 0) {
                continue;
            }
            final double[][] hist = HistogramPlot.calcHistogram(feature.values(), limits[i][0], limits[i][1], bins[i]);
            // Colour the points
            plot.setColor(colourMap.apply(n));
            plot.addPoints(hist[0], hist[1], Plot.BAR);
        }
        plot.updateImage();
    }
    createTrackDataTable(tracks, trackData, fitData, model, component, cal, colourMap);
// Analysis.
// Assign the original localisations to their track component.
// Q. What about the start/end not covered by the window?
// Save tracks as a dataset labelled with the sub-track ID?
// Output for the bound component and free components track parameters.
// Compute dwell times.
// Other ...
// Track analysis plugin:
// Extract all continuous segments of the same component.
// Produce MSD plot with error bars.
// Fit using FBM model.
}

Example 18 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method depthAnalysis.

/**
 * Depth analysis.
 *
 * @param allAssignments The assignments generated from running the filter (or null)
 * @param filter the filter
 * @return the assignments
 */
@Nullable
private ArrayList<FractionalAssignment[]> depthAnalysis(ArrayList<FractionalAssignment[]> allAssignments, DirectFilter filter) {
    if (!settings.depthRecallAnalysis || simulationParameters.fixedDepth) {
        return null;
    }
    // Build a histogram of the number of spots at different depths
    final double[] depths = fitResultData.depthStats.getValues();
    double[] limits = MathUtils.limits(depths);
    final int bins = HistogramPlot.getBinsSqrtRule(depths.length);
    final double[][] h1 = HistogramPlot.calcHistogram(depths, limits[0], limits[1], bins);
    final double[][] h2 = HistogramPlot.calcHistogram(fitResultData.depthFitStats.getValues(), limits[0], limits[1], bins);
    // manually to get the results that pass.
    if (allAssignments == null) {
        allAssignments = getAssignments(filter);
    }
    double[] depths2 = new double[results.size()];
    int count = 0;
    for (final FractionalAssignment[] assignments : allAssignments) {
        if (assignments == null) {
            continue;
        }
        for (int i = 0; i < assignments.length; i++) {
            final CustomFractionalAssignment c = (CustomFractionalAssignment) assignments[i];
            depths2[count++] = c.peak.getZPosition();
        }
    }
    depths2 = Arrays.copyOf(depths2, count);
    // Build a histogram using the same limits
    final double[][] h3 = HistogramPlot.calcHistogram(depths2, limits[0], limits[1], bins);
    // Convert pixel depth to nm
    for (int i = 0; i < h1[0].length; i++) {
        h1[0][i] *= simulationParameters.pixelPitch;
    }
    limits[0] *= simulationParameters.pixelPitch;
    limits[1] *= simulationParameters.pixelPitch;
    // Produce a histogram of the number of spots at each depth
    final String title1 = TITLE + " Depth Histogram";
    final Plot plot1 = new Plot(title1, "Depth (nm)", "Frequency");
    plot1.setLimits(limits[0], limits[1], 0, MathUtils.max(h1[1]));
    plot1.setColor(Color.black);
    plot1.addPoints(h1[0], h1[1], Plot.BAR);
    plot1.addLabel(0, 0, "Black = Spots; Blue = Fitted; Red = Filtered");
    plot1.setColor(Color.blue);
    plot1.addPoints(h1[0], h2[1], Plot.BAR);
    plot1.setColor(Color.red);
    plot1.addPoints(h1[0], h3[1], Plot.BAR);
    plot1.setColor(Color.magenta);
    ImageJUtils.display(title1, plot1, wo);
    // Interpolate
    final double halfBinWidth = (h1[0][1] - h1[0][0]) * 0.5;
    // Remove final value of the histogram as this is at the upper limit of the range (i.e. count
    // zero)
    h1[0] = Arrays.copyOf(h1[0], h1[0].length - 1);
    h1[1] = Arrays.copyOf(h1[1], h1[0].length);
    h2[1] = Arrays.copyOf(h2[1], h1[0].length);
    h3[1] = Arrays.copyOf(h3[1], h1[0].length);
    // TODO : Fix the smoothing since LOESS sometimes does not work.
    // Perhaps allow configuration of the number of histogram bins and the smoothing bandwidth.
    // Use minimum of 3 points for smoothing
    // Ensure we use at least x% of data
    final double bandwidth = Math.max(3.0 / h1[0].length, 0.15);
    final LoessInterpolator loess = new LoessInterpolator(bandwidth, 1);
    final PolynomialSplineFunction spline1 = loess.interpolate(h1[0], h1[1]);
    final PolynomialSplineFunction spline2 = loess.interpolate(h1[0], h2[1]);
    final PolynomialSplineFunction spline3 = loess.interpolate(h1[0], h3[1]);
    // Use a second interpolator in case the LOESS fails
    final LinearInterpolator lin = new LinearInterpolator();
    final PolynomialSplineFunction spline1b = lin.interpolate(h1[0], h1[1]);
    final PolynomialSplineFunction spline2b = lin.interpolate(h1[0], h2[1]);
    final PolynomialSplineFunction spline3b = lin.interpolate(h1[0], h3[1]);
    // Increase the number of points to show a smooth curve
    final double[] points = new double[bins * 5];
    limits = MathUtils.limits(h1[0]);
    final double interval = (limits[1] - limits[0]) / (points.length - 1);
    final double[] v = new double[points.length];
    final double[] v2 = new double[points.length];
    final double[] v3 = new double[points.length];
    for (int i = 0; i < points.length - 1; i++) {
        points[i] = limits[0] + i * interval;
        v[i] = getSplineValue(spline1, spline1b, points[i]);
        v2[i] = getSplineValue(spline2, spline2b, points[i]);
        v3[i] = getSplineValue(spline3, spline3b, points[i]);
        points[i] += halfBinWidth;
    }
    // Final point on the limit of the spline range
    final int ii = points.length - 1;
    v[ii] = getSplineValue(spline1, spline1b, limits[1]);
    v2[ii] = getSplineValue(spline2, spline2b, limits[1]);
    v3[ii] = getSplineValue(spline3, spline3b, limits[1]);
    points[ii] = limits[1] + halfBinWidth;
    // Calculate recall
    for (int i = 0; i < v.length; i++) {
        v2[i] = v2[i] / v[i];
        v3[i] = v3[i] / v[i];
    }
    final double halfSummaryDepth = settings.summaryDepth * 0.5;
    final String title2 = TITLE + " Depth Histogram (normalised)";
    final Plot plot2 = new Plot(title2, "Depth (nm)", "Recall");
    plot2.setLimits(limits[0] + halfBinWidth, limits[1] + halfBinWidth, 0, MathUtils.min(1, MathUtils.max(v2)));
    plot2.setColor(Color.black);
    plot2.addLabel(0, 0, "Blue = Fitted; Red = Filtered");
    plot2.setColor(Color.blue);
    plot2.addPoints(points, v2, Plot.LINE);
    plot2.setColor(Color.red);
    plot2.addPoints(points, v3, Plot.LINE);
    plot2.setColor(Color.magenta);
    if (-halfSummaryDepth - halfBinWidth >= limits[0]) {
        plot2.drawLine(-halfSummaryDepth, 0, -halfSummaryDepth, getSplineValue(spline3, spline3b, -halfSummaryDepth - halfBinWidth) / getSplineValue(spline1, spline1b, -halfSummaryDepth - halfBinWidth));
    }
    if (halfSummaryDepth - halfBinWidth <= limits[1]) {
        plot2.drawLine(halfSummaryDepth, 0, halfSummaryDepth, getSplineValue(spline3, spline3b, halfSummaryDepth - halfBinWidth) / getSplineValue(spline1, spline1b, halfSummaryDepth - halfBinWidth));
    }
    ImageJUtils.display(title2, plot2, wo);
    return allAssignments;
}

Example 19 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class Fire method runQEstimation.

@SuppressWarnings("null")
private void runQEstimation() {
    IJ.showStatus(pluginTitle + " ...");
    if (!showQEstimationInputDialog()) {
        return;
    }
    MemoryPeakResults inputResults = ResultsManager.loadInputResults(settings.inputOption, false, null, null);
    if (MemoryPeakResults.isEmpty(inputResults)) {
        IJ.error(pluginTitle, "No results could be loaded");
        return;
    }
    if (inputResults.getCalibration() == null) {
        IJ.error(pluginTitle, "The results are not calibrated");
        return;
    }
    inputResults = cropToRoi(inputResults);
    if (inputResults.size() < 2) {
        IJ.error(pluginTitle, "No results within the crop region");
        return;
    }
    initialise(inputResults, null);
    // We need localisation precision.
    // Build a histogram of the localisation precision.
    // Get the initial mean and SD and plot as a Gaussian.
    final PrecisionHistogram histogram = calculatePrecisionHistogram();
    if (histogram == null) {
        IJ.error(pluginTitle, "No localisation precision available.\n \nPlease choose " + PrecisionMethod.FIXED + " and enter a precision mean and SD.");
        return;
    }
    final StoredDataStatistics precision = histogram.precision;
    final double fourierImageScale = Settings.scaleValues[settings.imageScaleIndex];
    final int imageSize = Settings.imageSizeValues[settings.imageSizeIndex];
    // Create the image and compute the numerator of FRC.
    // Do not use the signal so results.size() is the number of localisations.
    IJ.showStatus("Computing FRC curve ...");
    final FireImages images = createImages(fourierImageScale, imageSize, false);
    // DEBUGGING - Save the two images to disk. Load the images into the Matlab
    // code that calculates the Q-estimation and make this plugin match the functionality.
    // IJ.save(new ImagePlus("i1", images.ip1), "/scratch/i1.tif");
    // IJ.save(new ImagePlus("i2", images.ip2), "/scratch/i2.tif");
    final Frc frc = new Frc();
    frc.setTrackProgress(progress);
    frc.setFourierMethod(fourierMethod);
    frc.setSamplingMethod(samplingMethod);
    frc.setPerimeterSamplingFactor(settings.perimeterSamplingFactor);
    final FrcCurve frcCurve = frc.calculateFrcCurve(images.ip1, images.ip2, images.nmPerPixel);
    if (frcCurve == null) {
        IJ.error(pluginTitle, "Failed to compute FRC curve");
        return;
    }
    IJ.showStatus("Running Q-estimation ...");
    // Note:
    // The method implemented here is based on Matlab code provided by Bernd Rieger.
    // The idea is to compute the spurious correlation component of the FRC Numerator
    // using an initial estimate of distribution of the localisation precision (assumed
    // to be Gaussian). This component is the contribution of repeat localisations of
    // the same molecule to the numerator and is modelled as an exponential decay
    // (exp_decay). The component is scaled by the Q-value which
    // is the average number of times a molecule is seen in addition to the first time.
    // At large spatial frequencies the scaled component should match the numerator,
    // i.e. at high resolution (low FIRE number) the numerator is made up of repeat
    // localisations of the same molecule and not actual structure in the image.
    // The best fit is where the numerator equals the scaled component, i.e. num / (q*exp_decay) ==
    // 1.
    // The FRC Numerator is plotted and Q can be determined by
    // adjusting Q and the precision mean and SD to maximise the cost function.
    // This can be done interactively by the user with the effect on the FRC curve
    // dynamically updated and displayed.
    // Compute the scaled FRC numerator
    final double qNorm = (1 / frcCurve.mean1 + 1 / frcCurve.mean2);
    final double[] frcnum = new double[frcCurve.getSize()];
    for (int i = 0; i < frcnum.length; i++) {
        final FrcCurveResult r = frcCurve.get(i);
        frcnum[i] = qNorm * r.getNumerator() / r.getNumberOfSamples();
    }
    // Compute the spatial frequency and the region for curve fitting
    final double[] q = Frc.computeQ(frcCurve, false);
    int low = 0;
    int high = q.length;
    while (high > 0 && q[high - 1] > settings.maxQ) {
        high--;
    }
    while (low < q.length && q[low] < settings.minQ) {
        low++;
    }
    // Require we fit at least 10% of the curve
    if (high - low < q.length * 0.1) {
        IJ.error(pluginTitle, "Not enough points for Q estimation");
        return;
    }
    // Obtain initial estimate of Q plateau height and decay.
    // This can be done by fitting the precision histogram and then fixing the mean and sigma.
    // Or it can be done by allowing the precision to be sampled and the mean and sigma
    // become parameters for fitting.
    // Check if we can sample precision values
    final boolean sampleDecay = precision != null && settings.sampleDecay;
    double[] expDecay;
    if (sampleDecay) {
        // Random sample of precision values from the distribution is used to
        // construct the decay curve
        final int[] sample = RandomUtils.sample(10000, precision.getN(), UniformRandomProviders.create());
        final double four_pi2 = 4 * Math.PI * Math.PI;
        final double[] pre = new double[q.length];
        for (int i = 1; i < q.length; i++) {
            pre[i] = -four_pi2 * q[i] * q[i];
        }
        // Sample
        final int n = sample.length;
        final double[] hq = new double[n];
        for (int j = 0; j < n; j++) {
            // Scale to SR pixels
            double s2 = precision.getValue(sample[j]) / images.nmPerPixel;
            s2 *= s2;
            for (int i = 1; i < q.length; i++) {
                hq[i] += StdMath.exp(pre[i] * s2);
            }
        }
        for (int i = 1; i < q.length; i++) {
            hq[i] /= n;
        }
        expDecay = new double[q.length];
        expDecay[0] = 1;
        for (int i = 1; i < q.length; i++) {
            final double sinc_q = sinc(Math.PI * q[i]);
            expDecay[i] = sinc_q * sinc_q * hq[i];
        }
    } else {
        // Note: The sigma mean and std should be in the units of super-resolution
        // pixels so scale to SR pixels
        expDecay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
    }
    // Smoothing
    double[] smooth;
    if (settings.loessSmoothing) {
        // Note: This computes the log then smooths it
        final double bandwidth = 0.1;
        final int robustness = 0;
        final double[] l = new double[expDecay.length];
        for (int i = 0; i < l.length; i++) {
            // Original Matlab code computes the log for each array.
            // This is equivalent to a single log on the fraction of the two.
            // Perhaps the two log method is more numerically stable.
            // l[i] = Math.log(Math.abs(frcnum[i])) - Math.log(exp_decay[i]);
            l[i] = Math.log(Math.abs(frcnum[i] / expDecay[i]));
        }
        try {
            final LoessInterpolator loess = new LoessInterpolator(bandwidth, robustness);
            smooth = loess.smooth(q, l);
        } catch (final Exception ex) {
            IJ.error(pluginTitle, "LOESS smoothing failed");
            return;
        }
    } else {
        // Note: This smooths the curve before computing the log
        final double[] norm = new double[expDecay.length];
        for (int i = 0; i < norm.length; i++) {
            norm[i] = frcnum[i] / expDecay[i];
        }
        // Median window of 5 == radius of 2
        final DoubleMedianWindow mw = DoubleMedianWindow.wrap(norm, 2);
        smooth = new double[expDecay.length];
        for (int i = 0; i < norm.length; i++) {
            smooth[i] = Math.log(Math.abs(mw.getMedian()));
            mw.increment();
        }
    }
    // Fit with quadratic to find the initial guess.
    // Note: example Matlab code frc_Qcorrection7.m identifies regions of the
    // smoothed log curve with low derivative and only fits those. The fit is
    // used for the final estimate. Fitting a subset with low derivative is not
    // implemented here since the initial estimate is subsequently optimised
    // to maximise a cost function.
    final Quadratic curve = new Quadratic();
    final SimpleCurveFitter fit = SimpleCurveFitter.create(curve, new double[2]);
    final WeightedObservedPoints points = new WeightedObservedPoints();
    for (int i = low; i < high; i++) {
        points.add(q[i], smooth[i]);
    }
    final double[] estimate = fit.fit(points.toList());
    double qvalue = StdMath.exp(estimate[0]);
    // This could be made an option. Just use for debugging
    final boolean debug = false;
    if (debug) {
        // Plot the initial fit and the fit curve
        final double[] qScaled = Frc.computeQ(frcCurve, true);
        final double[] line = new double[q.length];
        for (int i = 0; i < q.length; i++) {
            line[i] = curve.value(q[i], estimate);
        }
        final String title = pluginTitle + " Initial fit";
        final Plot plot = new Plot(title, "Spatial Frequency (nm^-1)", "FRC Numerator");
        final String label = String.format("Q = %.3f", qvalue);
        plot.addPoints(qScaled, smooth, Plot.LINE);
        plot.setColor(Color.red);
        plot.addPoints(qScaled, line, Plot.LINE);
        plot.setColor(Color.black);
        plot.addLabel(0, 0, label);
        ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT);
    }
    if (settings.fitPrecision) {
        // Q - Should this be optional?
        if (sampleDecay) {
            // If a sample of the precision was used to construct the data for the initial fit
            // then update the estimate using the fit result since it will be a better start point.
            histogram.sigma = precision.getStandardDeviation();
            // Normalise sum-of-squares to the SR pixel size
            final double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
            histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
        }
        // Do a multivariate fit ...
        final SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
        PointValuePair pair = null;
        final MultiPlateauness f = new MultiPlateauness(frcnum, q, low, high);
        final double[] initial = new double[] { histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, qvalue };
        pair = findMin(pair, opt, f, scale(initial, 0.1));
        pair = findMin(pair, opt, f, scale(initial, 0.5));
        pair = findMin(pair, opt, f, initial);
        pair = findMin(pair, opt, f, scale(initial, 2));
        pair = findMin(pair, opt, f, scale(initial, 10));
        if (pair != null) {
            final double[] point = pair.getPointRef();
            histogram.mean = point[0] * images.nmPerPixel;
            histogram.sigma = point[1] * images.nmPerPixel;
            qvalue = point[2];
        }
    } else {
        // If so then this should be optional.
        if (sampleDecay) {
            if (precisionMethod != PrecisionMethod.FIXED) {
                histogram.sigma = precision.getStandardDeviation();
                // Normalise sum-of-squares to the SR pixel size
                final double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
                histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
            }
            expDecay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
        }
        // Estimate spurious component by promoting plateauness.
        // The Matlab code used random initial points for a Simplex optimiser.
        // A Brent line search should be pretty deterministic so do simple repeats.
        // However it will proceed downhill so if the initial point is wrong then
        // it will find a sub-optimal result.
        final UnivariateOptimizer o = new BrentOptimizer(1e-3, 1e-6);
        final Plateauness f = new Plateauness(frcnum, expDecay, low, high);
        UnivariatePointValuePair result = null;
        result = findMin(result, o, f, qvalue, 0.1);
        result = findMin(result, o, f, qvalue, 0.2);
        result = findMin(result, o, f, qvalue, 0.333);
        result = findMin(result, o, f, qvalue, 0.5);
        // Do some Simplex repeats as well
        final SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
        result = findMin(result, opt, f, qvalue * 0.1);
        result = findMin(result, opt, f, qvalue * 0.5);
        result = findMin(result, opt, f, qvalue);
        result = findMin(result, opt, f, qvalue * 2);
        result = findMin(result, opt, f, qvalue * 10);
        if (result != null) {
            qvalue = result.getPoint();
        }
    }
    final QPlot qplot = new QPlot(frcCurve, qvalue, low, high);
    // Interactive dialog to estimate Q (blinking events per flourophore) using
    // sliders for the mean and standard deviation of the localisation precision.
    showQEstimationDialog(histogram, qplot, images.nmPerPixel);
    IJ.showStatus(pluginTitle + " complete");
}

Example 20 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class EmGainAnalysis method fit.

/**
 * Fit the EM-gain distribution (Gaussian * Gamma).
 *
 * @param histogram The distribution
 */
private void fit(int[] histogram) {
    final int[] limits = limits(histogram);
    final double[] x = getX(limits);
    final double[] y = getY(histogram, limits);
    Plot plot = new Plot(TITLE, "ADU", "Frequency");
    double yMax = MathUtils.max(y);
    plot.setLimits(limits[0], limits[1], 0, yMax);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    ImageJUtils.display(TITLE, plot);
    // Estimate remaining parameters.
    // Assuming a gamma_distribution(shape,scale) then mean = shape * scale
    // scale = gain
    // shape = Photons = mean / gain
    double mean = getMean(histogram) - settings.bias;
    // Note: if the bias is too high then the mean will be negative. Just move the bias.
    while (mean < 0) {
        settings.bias -= 1;
        mean += 1;
    }
    double photons = mean / settings.gain;
    if (settings.settingSimulate) {
        ImageJUtils.log("Simulated bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.settingBias, MathUtils.rounded(settings.settingGain), MathUtils.rounded(settings.settingNoise), MathUtils.rounded(settings.settingPhotons));
    }
    ImageJUtils.log("Estimate bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    final int max = (int) x[x.length - 1];
    double[] pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot.setColor(Color.blue);
    plot.addPoints(x, pg, Plot.LINE);
    ImageJUtils.display(TITLE, plot);
    // Perform a fit
    final CustomPowellOptimizer o = new CustomPowellOptimizer(1e-6, 1e-16, 1e-6, 1e-16);
    final double[] startPoint = new double[] { photons, settings.gain, settings.noise, settings.bias };
    int maxEval = 3000;
    final String[] paramNames = { "Photons", "Gain", "Noise", "Bias" };
    // Set bounds
    final double[] lower = new double[] { 0, 0.5 * settings.gain, 0, settings.bias - settings.noise };
    final double[] upper = new double[] { 2 * photons, 2 * settings.gain, settings.gain, settings.bias + settings.noise };
    // Restart until converged.
    // TODO - Maybe fix this with a better optimiser. This needs to be tested on real data.
    PointValuePair solution = null;
    for (int iter = 0; iter < 3; iter++) {
        IJ.showStatus("Fitting histogram ... Iteration " + iter);
        try {
            // Basic Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(fun), GoalType.MINIMIZE, new InitialGuess((solution == null) ? startPoint : solution.getPointRef()));
            if (solution == null || optimum.getValue() < solution.getValue()) {
                final double[] point = optimum.getPointRef();
                // Check the bounds
                for (int i = 0; i < point.length; i++) {
                    if (point[i] < lower[i] || point[i] > upper[i]) {
                        throw new ComputationException(String.format("Fit out of of estimated range: %s %f", paramNames[i], point[i]));
                    }
                }
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
        try {
            // Bounded Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(fun, lower, upper);
            PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded((solution == null) ? startPoint : solution.getPointRef())));
            final double[] point = adapter.unboundedToBounded(optimum.getPointRef());
            optimum = new PointValuePair(point, optimum.getValue());
            if (solution == null || optimum.getValue() < solution.getValue()) {
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Bounded Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
    }
    ImageJUtils.finished();
    if (solution == null) {
        ImageJUtils.log("Failed to fit the distribution");
        return;
    }
    final double[] point = solution.getPointRef();
    photons = point[0];
    settings.gain = point[1];
    settings.noise = point[2];
    settings.bias = (int) Math.round(point[3]);
    final String label = String.format("Fitted bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    ImageJUtils.log(label);
    if (settings.settingSimulate) {
        ImageJUtils.log("Relative Error bias=%s, gain=%s, noise=%s, photons=%s", MathUtils.rounded(relativeError(settings.bias, settings.settingBias)), MathUtils.rounded(relativeError(settings.gain, settings.settingGain)), MathUtils.rounded(relativeError(settings.noise, settings.settingNoise)), MathUtils.rounded(relativeError(photons, settings.settingPhotons)));
    }
    // Show the PoissonGammaGaussian approximation
    double[] approxValues = null;
    if (settings.showApproximation) {
        approxValues = new double[x.length];
        final PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / settings.gain, settings.noise);
        final double expected = photons * settings.gain;
        for (int i = 0; i < approxValues.length; i++) {
            approxValues[i] = fun.likelihood(x[i] - settings.bias, expected);
        }
        yMax = MathUtils.maxDefault(yMax, approxValues);
    }
    // Replot
    pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot = new Plot(TITLE, "ADU", "Frequency");
    plot.setLimits(limits[0], limits[1], 0, yMax * 1.05);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    plot.setColor(Color.red);
    plot.addPoints(x, pg, Plot.LINE);
    plot.addLabel(0, 0, label);
    if (settings.showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, approxValues, Plot.LINE);
    }
    ImageJUtils.display(TITLE, plot);
}