Examples with MemoryPeakResults uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults used on opensource projects

Example 6 with MemoryPeakResults

use of uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class CreateData method savePulses.

/**
 * Create a set of results that represent the molecule continuous on-times (pulses).
 *
 * @param localisations the localisations
 * @param results the results
 */
@SuppressWarnings("null")
private void savePulses(List<LocalisationModel> localisations, MemoryPeakResults results) {
    sortLocalisationsByIdThenTime(localisations);
    final MemoryPeakResults traceResults = copyMemoryPeakResults("Pulses");
    LocalisationModel start = null;
    int currentId = -1;
    int count = 0;
    float[] params = Gaussian2DPeakResultHelper.createTwoAxisParams(0, 0, 0, 0, 0, 0, 0);
    final int isx = Gaussian2DPeakResultHelper.INDEX_SX;
    final int isy = Gaussian2DPeakResultHelper.INDEX_SY;
    double noise = 0;
    int lastT = -1;
    for (final LocalisationModel localisation : localisations) {
        if (currentId != localisation.getId() || lastT + 1 != localisation.getTime()) {
            if (count > 0) {
                params[PeakResult.BACKGROUND] /= count;
                params[PeakResult.X] /= count;
                params[PeakResult.Y] /= count;
                params[isx] /= count;
                params[isy] /= count;
                final ExtendedPeakResult p = new ExtendedPeakResult(start.getTime(), (int) Math.round(start.getX()), (int) Math.round(start.getY()), 0, 0, (float) (Math.sqrt(noise)), 0, params, null, lastT, currentId);
                // if (p.getPrecision(107, 1) > 2000)
                // {
                // System.out.printf("Weird precision = %g (%d)\n", p.getPrecision(107, 1), n);
                // }
                traceResults.add(p);
            }
            start = localisation;
            currentId = localisation.getId();
            count = 0;
            params = new float[7];
            noise = 0;
        }
        final double[] data = localisation.getData();
        params[PeakResult.BACKGROUND] += data[0];
        params[PeakResult.X] += localisation.getX();
        params[PeakResult.Y] += localisation.getY();
        params[PeakResult.INTENSITY] += localisation.getIntensity();
        noise += data[1] * data[1];
        params[isx] += data[2];
        params[isy] += data[3];
        count++;
        lastT = localisation.getTime();
    }
    // Final pulse
    if (count > 0) {
        params[PeakResult.BACKGROUND] /= count;
        params[PeakResult.X] /= count;
        params[PeakResult.Y] /= count;
        params[isx] /= count;
        params[isy] /= count;
        traceResults.add(new ExtendedPeakResult(start.getTime(), (int) Math.round(start.getX()), (int) Math.round(start.getY()), 0, 0, (float) (Math.sqrt(noise)), 0, params, null, lastT, currentId));
    }
    traceResults.end();
}

Example 7 with MemoryPeakResults

use of uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class PcPalmClusters method doClustering.

/**
 * Extract the results from the PCPALM molecules using the area ROI and then do clustering to
 * obtain the histogram of molecules per cluster.
 *
 * @return the histogram data
 */
private HistogramData doClustering() {
    // Perform clustering analysis to generate the histogram of cluster sizes
    final PcPalmAnalysis analysis = new PcPalmAnalysis();
    final List<Molecule> molecules = analysis.cropToRoi(WindowManager.getCurrentImage(), moleculesResults);
    if (molecules.size() < 2) {
        error("No results within the crop region");
        return null;
    }
    ImageJUtils.log("Using %d molecules (Density = %s um^-2) @ %s nm", molecules.size(), MathUtils.rounded(molecules.size() / analysis.croppedArea), MathUtils.rounded(settings.distance));
    final long s1 = System.nanoTime();
    final ClusteringEngine engine = new ClusteringEngine(1, clusteringAlgorithm, SimpleImageJTrackProgress.getInstance());
    if (settings.multiThread) {
        engine.setThreadCount(Prefs.getThreads());
    }
    engine.setTracker(SimpleImageJTrackProgress.getInstance());
    IJ.showStatus("Clustering ...");
    final List<Cluster> clusters = engine.findClusters(convertToPoint(molecules), settings.distance);
    IJ.showStatus("");
    if (clusters == null) {
        ImageJUtils.log("Aborted");
        return null;
    }
    numberOfMolecules = molecules.size();
    ImageJUtils.log("Finished : %d total clusters (%s ms)", clusters.size(), MathUtils.rounded((System.nanoTime() - s1) / 1e6));
    // Save cluster centroids to a results set in memory. Then they can be plotted.
    final MemoryPeakResults results = new MemoryPeakResults(clusters.size());
    results.setName(TITLE);
    // Set an arbitrary calibration so that the lifetime of the results is stored in the exposure
    // time
    // The results will be handled as a single mega-frame containing all localisation.
    results.setCalibration(CalibrationHelper.create(100, 1, moleculesResults.seconds * 1000));
    int id = 0;
    for (final Cluster c : clusters) {
        results.add(new ExtendedPeakResult((float) c.getX(), (float) c.getY(), c.getSize(), ++id));
    }
    MemoryPeakResults.addResults(results);
    // Get the data for fitting
    final float[] values = new float[clusters.size()];
    for (int i = 0; i < values.length; i++) {
        values[i] = clusters.get(i).getSize();
    }
    final float yMax = (int) Math.ceil(MathUtils.max(values));
    final int nBins = (int) (yMax + 1);
    final float[][] hist = HistogramPlot.calcHistogram(values, 0, yMax, nBins);
    final HistogramData histogramData = (settings.calibrateHistogram) ? new HistogramData(hist, settings.frames, settings.area, Settings.UNITS[settings.units]) : new HistogramData(hist);
    saveHistogram(histogramData);
    return histogramData;
}

Example 8 with MemoryPeakResults

use of uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method traceMolecules.

/**
 * Trace localisations.
 *
 * @param results The results
 * @param precisions the precisions
 * @param distance The distance threshold (nm)
 * @param time The time threshold (frames)
 * @param singles a list of the singles (not grouped into molecules)
 * @return a list of molecules
 */
private static ArrayList<Molecule> traceMolecules(MemoryPeakResults results, double[] precisions, double distance, int time, ArrayList<Molecule> singles) {
    // These plugins are not really supported so just leave them to throw an exception if
    // the data cannot be handled
    final TypeConverter<IntensityUnit> ic = results.getCalibrationReader().getIntensityConverter(IntensityUnit.PHOTON);
    final TypeConverter<DistanceUnit> dc = results.getCalibrationReader().getDistanceConverter(DistanceUnit.NM);
    // Create a new dataset with the precision
    final MemoryPeakResults results2 = new MemoryPeakResults(results.size());
    for (int i = 0, size = results.size(); i < size; i++) {
        final AttributePeakResult peak2 = new AttributePeakResult(results.get(i));
        peak2.setPrecision(precisions[i]);
        results2.add(peak2);
    }
    final TraceManager tm = new TraceManager(results2);
    final double distanceThreshold = dc.convertBack(distance);
    tm.traceMolecules(distanceThreshold, time);
    final Trace[] traces = tm.getTraces();
    final ArrayList<Molecule> molecules = new ArrayList<>(traces.length);
    for (final Trace t : traces) {
        final double p = t.getLocalisationPrecision(dc);
        final float[] centroid = t.getCentroid();
        final List<Molecule> list = t.size() == 1 ? singles : molecules;
        list.add(new Molecule(dc.convert(centroid[0]), dc.convert(centroid[1]), p, ic.convert(t.getSignal())));
    }
    log("  %d localisations traced to %d molecules (%d singles, %d traces) using d=%.2f nm," + " t=%d frames (%s s)", results.size(), molecules.size() + singles.size(), singles.size(), molecules.size(), distance, time, MathUtils.rounded(time * results.getCalibrationReader().getExposureTime() / 1000.0));
    return molecules;
}

Example 9 with MemoryPeakResults

use of uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class BlinkEstimatorTest method estimateBlinking.

private TIntHashSet estimateBlinking(UniformRandomProvider rg, double blinkingRate, double ton, double toff, int particles, double fixedFraction, boolean timeAtLowerBound, boolean doAssert) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MAXIMUM));
    final SpatialIllumination activationIllumination = new UniformIllumination(100);
    int totalSteps = 100;
    final double eAct = totalSteps * 0.3 * activationIllumination.getAveragePhotons();
    final ImageModel imageModel = new ActivationEnergyImageModel(eAct, activationIllumination, ton, 0, toff, 0, blinkingRate, rg);
    final double[] max = new double[] { 256, 256, 32 };
    final double[] min = new double[3];
    final SpatialDistribution distribution = new UniformDistribution(min, max, rg.nextInt());
    final List<CompoundMoleculeModel> compounds = new ArrayList<>(1);
    final CompoundMoleculeModel c = new CompoundMoleculeModel(1, 0, 0, 0, Arrays.asList(new MoleculeModel(0, 0, 0, 0)));
    c.setDiffusionRate(diffusionRate);
    c.setDiffusionType(DiffusionType.RANDOM_WALK);
    compounds.add(c);
    final List<CompoundMoleculeModel> molecules = imageModel.createMolecules(compounds, particles, distribution, false);
    // Activate fluorophores
    final List<? extends FluorophoreSequenceModel> fluorophores = imageModel.createFluorophores(molecules, totalSteps);
    totalSteps = checkTotalSteps(totalSteps, fluorophores);
    final List<LocalisationModel> localisations = imageModel.createImage(molecules, fixedFraction, totalSteps, photons, 0.5, false);
    // // Remove localisations to simulate missed counts.
    // List<LocalisationModel> newLocalisations = new
    // ArrayList<LocalisationModel>(localisations.size());
    // boolean[] id = new boolean[fluorophores.size() + 1];
    // Statistics photonStats = new Statistics();
    // for (LocalisationModel l : localisations)
    // {
    // photonStats.add(l.getIntensity());
    // // Remove by intensity threshold and optionally at random.
    // if (l.getIntensity() < minPhotons || rand.nextDouble() < pDelete)
    // continue;
    // newLocalisations.add(l);
    // id[l.getId()] = true;
    // }
    // localisations = newLocalisations;
    // logger.info("Photons = %f", photonStats.getMean());
    // 
    // List<FluorophoreSequenceModel> newFluorophores = new
    // ArrayList<FluorophoreSequenceModel>(fluorophores.size());
    // for (FluorophoreSequenceModel f : fluorophores)
    // {
    // if (id[f.getId()])
    // newFluorophores.add(f);
    // }
    // fluorophores = newFluorophores;
    final MemoryPeakResults results = new MemoryPeakResults();
    final CalibrationWriter calibration = new CalibrationWriter();
    calibration.setNmPerPixel(pixelPitch);
    calibration.setExposureTime(msPerFrame);
    calibration.setCountPerPhoton(1);
    results.setCalibration(calibration.getCalibration());
    results.setPsf(PsfHelper.create(PSFType.ONE_AXIS_GAUSSIAN_2D));
    final float b = 0;
    float intensity;
    final float z = 0;
    for (final LocalisationModel l : localisations) {
        // Remove by intensity threshold and optionally at random.
        if (l.getIntensity() < minPhotons || rg.nextDouble() < probabilityDelete) {
            continue;
        }
        final int frame = l.getTime();
        intensity = (float) l.getIntensity();
        final float x = (float) l.getX();
        final float y = (float) l.getY();
        final float[] params = Gaussian2DPeakResultHelper.createParams(b, intensity, x, y, z, psfWidth);
        results.add(frame, 0, 0, 0, 0, 0, 0, params, null);
    }
    // Add random localisations
    // Intensity doesn't matter at the moment for tracing
    intensity = (float) photons;
    for (int i = (int) (localisations.size() * probabilityAdd); i-- > 0; ) {
        final int frame = 1 + rg.nextInt(totalSteps);
        final float x = (float) (rg.nextDouble() * max[0]);
        final float y = (float) (rg.nextDouble() * max[1]);
        final float[] params = Gaussian2DPeakResultHelper.createParams(b, intensity, x, y, z, psfWidth);
        results.add(frame, 0, 0, 0, 0, 0, 0, params, null);
    }
    // Get actual simulated stats ...
    final Statistics statsNBlinks = new Statistics();
    final Statistics statsTOn = new Statistics();
    final Statistics statsTOff = new Statistics();
    final Statistics statsSampledNBlinks = new Statistics();
    final Statistics statsSampledTOn = new Statistics();
    final StoredDataStatistics statsSampledTOff = new StoredDataStatistics();
    for (final FluorophoreSequenceModel f : fluorophores) {
        statsNBlinks.add(f.getNumberOfBlinks());
        statsTOn.add(f.getOnTimes());
        statsTOff.add(f.getOffTimes());
        final int[] on = f.getSampledOnTimes();
        statsSampledNBlinks.add(on.length);
        statsSampledTOn.add(on);
        statsSampledTOff.add(f.getSampledOffTimes());
    }
    logger.info(FunctionUtils.getSupplier("N = %d (%d), N-blinks = %f, tOn = %f, tOff = %f, Fixed = %f", fluorophores.size(), localisations.size(), blinkingRate, ton, toff, fixedFraction));
    logger.info(FunctionUtils.getSupplier("Actual N-blinks = %f (%f), tOn = %f (%f), tOff = %f (%f), 95%% = %f, max = %f", statsNBlinks.getMean(), statsSampledNBlinks.getMean(), statsTOn.getMean(), statsSampledTOn.getMean(), statsTOff.getMean(), statsSampledTOff.getMean(), statsSampledTOff.getStatistics().getPercentile(95), statsSampledTOff.getStatistics().getMax()));
    logger.info("-=-=--=-");
    final BlinkEstimator be = new BlinkEstimator();
    be.setMaxDarkTime((int) (toff * 10));
    be.setMsPerFrame(msPerFrame);
    be.setRelativeDistance(false);
    final double d = ImageModel.getRandomMoveDistance(diffusionRate);
    be.setSearchDistance((fixedFraction < 1) ? Math.sqrt(2 * d * d) * 3 : 0);
    be.setTimeAtLowerBound(timeAtLowerBound);
    // Assertions.assertTrue("Max dark time must exceed the dark time of the data (otherwise no
    // plateau)",
    // be.maxDarkTime > statsSampledTOff.getStatistics().getMax());
    final int nMolecules = fluorophores.size();
    if (usePopulationStatistics) {
        blinkingRate = statsNBlinks.getMean();
        toff = statsTOff.getMean();
    } else {
        blinkingRate = statsSampledNBlinks.getMean();
        toff = statsSampledTOff.getMean();
    }
    // See if any fitting regime gets a correct answer
    final TIntHashSet ok = new TIntHashSet();
    for (int numberOfFittedPoints = MIN_FITTED_POINTS; numberOfFittedPoints <= MAX_FITTED_POINTS; numberOfFittedPoints++) {
        be.setNumberOfFittedPoints(numberOfFittedPoints);
        be.computeBlinkingRate(results, true);
        final double moleculesError = DoubleEquality.relativeError(nMolecules, be.getNMolecules());
        final double blinksError = DoubleEquality.relativeError(blinkingRate, be.getNBlinks());
        final double offError = DoubleEquality.relativeError(toff * msPerFrame, be.getTOff());
        logger.info(FunctionUtils.getSupplier("Error %d: N = %f, blinks = %f, tOff = %f : %f", numberOfFittedPoints, moleculesError, blinksError, offError, (moleculesError + blinksError + offError) / 3));
        if (moleculesError < relativeError && blinksError < relativeError && offError < relativeError) {
            ok.add(numberOfFittedPoints);
            logger.info("-=-=--=-");
            logger.info(FunctionUtils.getSupplier("*** Correct at %d fitted points ***", numberOfFittedPoints));
            if (doAssert) {
                break;
            }
        }
    // if (!be.isIncreaseNFittedPoints())
    // break;
    }
    logger.info("-=-=--=-");
    if (doAssert) {
        Assertions.assertFalse(ok.isEmpty());
    }
    // relativeError);
    return ok;
}

Example 10 with MemoryPeakResults

use of uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class PsfCreator method fitPsf.

/**
 * Fit the new PSF image and show a graph of the amplitude/width.
 *
 * @param psfStack the psf stack
 * @param loess the loess
 * @param cz the cz
 * @param averageRange the average range
 * @param fitCom the fit com
 * @return The width of the PSF in the z-centre
 */
private double fitPsf(ImageStack psfStack, LoessInterpolator loess, int cz, double averageRange, final double[][] fitCom) {
    IJ.showStatus("Fitting final PSF");
    // is not appropriate for a normalised PSF.
    if (fitConfig.getFitSolver() != FitSolver.LVM_LSE) {
        ImageJUtils.log("  " + FitProtosHelper.getName(fitConfig.getFitSolver()) + " is not appropriate for final PSF fitting.");
        ImageJUtils.log("  Switching to Least Square Estimation");
        fitConfig.setFitSolver(FitSolver.LVM_LSE);
        if (settings.getInteractiveMode()) {
            // This assumes the LVM does not need the calibration
            PeakFit.configureFitSolver(config, null, null, 0);
        }
    }
    // Update the box radius since this is used in the fitSpot method.
    boxRadius = psfStack.getWidth() / 2;
    final int x = boxRadius;
    final int y = boxRadius;
    final double shift = fitConfig.getCoordinateShiftFactor();
    // Scale the PSF
    final PSF.Builder localPsf = fitConfig.getPsf().toBuilder();
    for (int i = 0; i < localPsf.getParametersCount(); i++) {
        final PSFParameter param = localPsf.getParameters(i);
        if (param.getUnit() == PSFParameterUnit.DISTANCE) {
            final PSFParameter.Builder b = param.toBuilder();
            b.setValue(b.getValue() * settings.getMagnification());
            localPsf.setParameters(i, b);
        }
    }
    fitConfig.setPsf(localPsf.build());
    // Need to be updated after the widths have been set
    fitConfig.setCoordinateShiftFactor(shift);
    fitConfig.setBackgroundFitting(false);
    // Since the PSF will be normalised remove the camera calibration
    fitConfig.setCameraType(CameraType.CAMERA_TYPE_NA);
    fitConfig.setMinPhotons(0);
    fitConfig.setBias(0);
    fitConfig.setGain(1);
    // No complex filtering so we get a fit. It should be easy to fit anyway.
    fitConfig.setPrecisionThreshold(0);
    fitConfig.setDirectFilter(null);
    // fitConfig.setDisableSimpleFilter(true);
    // Use this for debugging the fit
    // fitConfig.setLog(uk.ac.sussex.gdsc.core.ij.ImageJPluginLoggerHelper.getDefaultLogger());
    final MemoryPeakResults results = fitSpot(psfStack, psfStack.getWidth(), psfStack.getHeight(), x, y);
    if (results.size() < 5) {
        ImageJUtils.log("  Final PSF: Not enough fit results %d", results.size());
        return 0;
    }
    // Get the results for the spot centre and width
    final double[] z = new double[results.size()];
    final double[] xCoord = new double[z.length];
    final double[] yCoord = new double[z.length];
    final double[] sd = new double[z.length];
    final double[] a = new double[z.length];
    // Set limits for the fit
    final float maxWidth = (float) (Math.max(fitConfig.getInitialXSd(), fitConfig.getInitialYSd()) * settings.getMagnification() * 4);
    // PSF is normalised to 1
    final float maxSignal = 2;
    final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
    wp.getWxWy();
    final HeightResultProcedure hp = new HeightResultProcedure(results, IntensityUnit.COUNT);
    hp.getH();
    final Counter counter = new Counter();
    final Counter counterOk = new Counter();
    // We have fit the results so they will be in the preferred units
    results.forEach((PeakResultProcedure) peak -> {
        int index = counter.getAndIncrement();
        final float w = Math.max(wp.wx[index], wp.wy[index]);
        if (peak.getIntensity() > maxSignal || w > maxWidth) {
            return;
        }
        index = counterOk.getAndIncrement();
        z[index] = peak.getFrame();
        fitCom[0][peak.getFrame() - 1] = xCoord[index] = peak.getXPosition() - x;
        fitCom[1][peak.getFrame() - 1] = yCoord[index] = peak.getYPosition() - y;
        sd[index] = w;
        a[index] = hp.heights[index];
    });
    // Truncate
    final double[] z2 = Arrays.copyOf(z, counter.getCount());
    final double[] xCoord2 = Arrays.copyOf(xCoord, z2.length);
    final double[] yCoord2 = Arrays.copyOf(yCoord, z2.length);
    final double[] sd2 = Arrays.copyOf(sd, z2.length);
    final double[] a2 = Arrays.copyOf(a, z2.length);
    // Extract the average smoothed range from the individual fits
    final int r = (int) Math.ceil(averageRange / 2);
    int start = 0;
    int stop = z2.length - 1;
    for (int j = 0; j < z2.length; j++) {
        if (z2[j] > cz - r) {
            start = j;
            break;
        }
    }
    for (int j = z2.length; j-- > 0; ) {
        if (z2[j] < cz + r) {
            stop = j;
            break;
        }
    }
    // Extract xy centre coords and smooth
    double[] smoothX = new double[stop - start + 1];
    double[] smoothY = new double[smoothX.length];
    double[] smoothSd = new double[smoothX.length];
    double[] smoothA = new double[smoothX.length];
    final double[] newZ = new double[smoothX.length];
    int smoothCzIndex = 0;
    for (int j = start, k = 0; j <= stop; j++, k++) {
        smoothX[k] = xCoord2[j];
        smoothY[k] = yCoord2[j];
        smoothSd[k] = sd2[j];
        smoothA[k] = a2[j];
        newZ[k] = z2[j];
        if (newZ[k] == cz) {
            smoothCzIndex = k;
        }
    }
    smoothX = loess.smooth(newZ, smoothX);
    smoothY = loess.smooth(newZ, smoothY);
    smoothSd = loess.smooth(newZ, smoothSd);
    smoothA = loess.smooth(newZ, smoothA);
    // Update the widths and positions using the magnification
    final double scale = 1.0 / settings.getMagnification();
    for (int j = 0; j < xCoord2.length; j++) {
        xCoord2[j] *= scale;
        yCoord2[j] *= scale;
        sd2[j] *= scale;
    }
    for (int j = 0; j < smoothX.length; j++) {
        smoothX[j] *= scale;
        smoothY[j] *= scale;
        smoothSd[j] *= scale;
    }
    showPlots(z2, a2, newZ, smoothA, xCoord2, yCoord2, sd2, newZ, smoothX, smoothY, smoothSd, cz);
    // Store the data for replotting
    this.z = z2;
    this.amplitude = a2;
    this.smoothAz = newZ;
    this.smoothA = smoothA;
    this.xCoord = xCoord2;
    this.yCoord = yCoord2;
    this.sd = sd2;
    this.newZ = newZ;
    this.smoothX = smoothX;
    this.smoothY = smoothY;
    this.smoothSd = smoothSd;
    // maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
    return smoothSd[smoothCzIndex];
}