Examples with LocalisationModel gdsc.smlm.model.LocalisationModel used on opensource projects

Example 6 with LocalisationModel

use of gdsc.smlm.model.LocalisationModel in project GDSC-SMLM by aherbert.

the class CreateData method createImagePSF.

/**
	 * Create a PSF model from the image that contains all the z-slices needed to draw the given localisations
	 * 
	 * @param localisationSets
	 * @return
	 */
private ImagePSFModel createImagePSF(List<LocalisationModelSet> localisationSets) {
    ImagePlus imp = WindowManager.getImage(settings.psfImageName);
    if (imp == null) {
        IJ.error(TITLE, "Unable to create the PSF model from image: " + settings.psfImageName);
        return null;
    }
    try {
        Object o = XmlUtils.fromXML(imp.getProperty("Info").toString());
        if (!(o != null && o instanceof PSFSettings))
            throw new RuntimeException("Unknown PSF settings for image: " + imp.getTitle());
        PSFSettings psfSettings = (PSFSettings) o;
        // Check all the settings have values
        if (psfSettings.nmPerPixel <= 0)
            throw new RuntimeException("Missing nmPerPixel calibration settings for image: " + imp.getTitle());
        if (psfSettings.nmPerSlice <= 0)
            throw new RuntimeException("Missing nmPerSlice calibration settings for image: " + imp.getTitle());
        if (psfSettings.zCentre <= 0)
            throw new RuntimeException("Missing zCentre calibration settings for image: " + imp.getTitle());
        if (psfSettings.fwhm <= 0)
            throw new RuntimeException("Missing FWHM calibration settings for image: " + imp.getTitle());
        // To save memory construct the Image PSF using only the slices that are within 
        // the depth of field of the simulation
        double minZ = Double.POSITIVE_INFINITY, maxZ = Double.NEGATIVE_INFINITY;
        for (LocalisationModelSet l : localisationSets) {
            for (LocalisationModel m : l.getLocalisations()) {
                final double z = m.getZ();
                if (minZ > z)
                    minZ = z;
                if (maxZ < z)
                    maxZ = z;
            }
        }
        int nSlices = imp.getStackSize();
        // z-centre should be an index and not the ImageJ slice number so subtract 1
        int zCentre = psfSettings.zCentre - 1;
        // Calculate the start/end slices to cover the depth of field
        // This logic must match the ImagePSFModel.
        final double unitsPerSlice = psfSettings.nmPerSlice / settings.pixelPitch;
        // We assume the PSF was imaged axially with increasing z-stage position (moving the stage 
        // closer to the objective). Thus higher z-coordinate are for higher slice numbers.
        int lower = (int) Math.round(minZ / unitsPerSlice) + zCentre;
        int upper = (int) Math.round(maxZ / unitsPerSlice) + zCentre;
        upper = (upper < 0) ? 0 : (upper >= nSlices) ? nSlices - 1 : upper;
        lower = (lower < 0) ? 0 : (lower >= nSlices) ? nSlices - 1 : lower;
        // Image PSF requires the z-centre for normalisation
        if (!(lower <= zCentre && upper >= zCentre)) {
            // Ensure we include the zCentre
            lower = Math.min(lower, zCentre);
            upper = Math.max(upper, zCentre);
        }
        final double noiseFraction = 1e-3;
        final ImagePSFModel model = new ImagePSFModel(extractImageStack(imp, lower, upper), zCentre - lower, psfSettings.nmPerPixel / settings.pixelPitch, unitsPerSlice, psfSettings.fwhm, noiseFraction);
        // Add the calibrated centres
        if (psfSettings.offset != null) {
            final int sliceOffset = lower + 1;
            for (PSFOffset offset : psfSettings.offset) {
                model.setRelativeCentre(offset.slice - sliceOffset, offset.cx, offset.cy);
            }
        }
        // Initialise the HWHM table so that it can be cloned
        model.initialiseHWHM();
        return model;
    } catch (Exception e) {
        IJ.error(TITLE, "Unable to create the image PSF model:\n" + e.getMessage());
        return null;
    }
}

Example 7 with LocalisationModel

use of gdsc.smlm.model.LocalisationModel in project GDSC-SMLM by aherbert.

the class CreateData method removeFilteredFluorophores.

/**
	 * Remove all fluorophores which were not drawn
	 * 
	 * @param fluorophores
	 * @param localisations
	 * @return
	 */
private List<? extends FluorophoreSequenceModel> removeFilteredFluorophores(List<? extends FluorophoreSequenceModel> fluorophores, List<LocalisationModel> localisations) {
    if (fluorophores == null)
        return null;
    // movingMolecules will be created with an initial capacity to hold all the unique IDs
    TIntHashSet idSet = new TIntHashSet((movingMolecules != null) ? movingMolecules.capacity() : 0);
    for (LocalisationModel l : localisations) idSet.add(l.getId());
    List<FluorophoreSequenceModel> newFluorophores = new ArrayList<FluorophoreSequenceModel>(idSet.size());
    for (FluorophoreSequenceModel f : fluorophores) {
        if (idSet.contains(f.getId()))
            newFluorophores.add(f);
    }
    return newFluorophores;
}

Example 8 with LocalisationModel

use of gdsc.smlm.model.LocalisationModel in project GDSC-SMLM by aherbert.

the class CreateData method showSummary.

private double showSummary(List<? extends FluorophoreSequenceModel> fluorophores, List<LocalisationModel> localisations) {
    IJ.showStatus("Calculating statistics ...");
    createSummaryTable();
    Statistics[] stats = new Statistics[NAMES.length];
    for (int i = 0; i < stats.length; i++) {
        stats[i] = (settings.showHistograms || alwaysRemoveOutliers[i]) ? new StoredDataStatistics() : new Statistics();
    }
    // Find the largest timepoint
    ImagePlus outputImp = WindowManager.getImage(benchmarkImageId);
    int nFrames;
    if (outputImp == null) {
        sortLocalisationsByTime(localisations);
        nFrames = localisations.get(localisations.size() - 1).getTime();
    } else {
        nFrames = outputImp.getStackSize();
    }
    int[] countHistogram = new int[nFrames + 1];
    // Use the localisations that were drawn to create the sampled on/off times
    rebuildNeighbours(localisations);
    // Assume that there is at least one localisation
    LocalisationModel first = localisations.get(0);
    // The current localisation
    int currentId = first.getId();
    // The last time this localisation was on
    int lastT = first.getTime();
    // Number of blinks
    int blinks = 0;
    // On-time of current pulse
    int currentT = 0;
    double signal = 0;
    final double centreOffset = settings.size * 0.5;
    // Used to convert the sampled times in frames into seconds
    final double framesPerSecond = 1000.0 / settings.exposureTime;
    final double gain = (settings.getTotalGain() > 0) ? settings.getTotalGain() : 1;
    for (LocalisationModel l : localisations) {
        if (l.getData() == null)
            System.out.println("No localisation data. This should not happen!");
        final double noise = (l.getData() != null) ? l.getData()[1] : 1;
        final double intensity = (l.getData() != null) ? l.getData()[4] : l.getIntensity();
        final double intensityInPhotons = intensity / gain;
        // Q. What if the noise is zero, i.e. no background photon / read noise?
        // Just ignore it at current.
        final double snr = intensity / noise;
        stats[SIGNAL].add(intensityInPhotons);
        stats[NOISE].add(noise / gain);
        if (noise != 0)
            stats[SNR].add(snr);
        //if (l.isContinuous())
        if (l.getNext() != null && l.getPrevious() != null) {
            stats[SIGNAL_CONTINUOUS].add(intensityInPhotons);
            if (noise != 0)
                stats[SNR_CONTINUOUS].add(snr);
        }
        int id = l.getId();
        // Check if this a new fluorophore
        if (currentId != id) {
            // Add previous fluorophore
            stats[SAMPLED_BLINKS].add(blinks);
            stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
            stats[TOTAL_SIGNAL].add(signal);
            // Reset
            blinks = 0;
            currentT = 1;
            currentId = id;
            signal = intensityInPhotons;
        } else {
            signal += intensityInPhotons;
            // Check if the current fluorophore pulse is broken (i.e. a blink)
            if (l.getTime() - 1 > lastT) {
                blinks++;
                stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
                currentT = 1;
                stats[SAMPLED_T_OFF].add(((l.getTime() - 1) - lastT) / framesPerSecond);
            } else {
                // Continuous on-time
                currentT++;
            }
        }
        lastT = l.getTime();
        countHistogram[lastT]++;
        stats[X].add((l.getX() - centreOffset) * settings.pixelPitch);
        stats[Y].add((l.getY() - centreOffset) * settings.pixelPitch);
        stats[Z].add(l.getZ() * settings.pixelPitch);
    }
    // Final fluorophore
    stats[SAMPLED_BLINKS].add(blinks);
    stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
    stats[TOTAL_SIGNAL].add(signal);
    // Samples per frame
    for (int t = 1; t < countHistogram.length; t++) stats[SAMPLES].add(countHistogram[t]);
    if (fluorophores != null) {
        for (FluorophoreSequenceModel f : fluorophores) {
            stats[BLINKS].add(f.getNumberOfBlinks());
            // On-time
            for (double t : f.getOnTimes()) stats[T_ON].add(t);
            // Off-time
            for (double t : f.getOffTimes()) stats[T_OFF].add(t);
        }
    } else {
        // show no blinks
        stats[BLINKS].add(0);
        stats[T_ON].add(1);
    //stats[T_OFF].add(0);
    }
    if (results != null) {
        final boolean emCCD = (settings.getEmGain() > 1);
        // Convert depth-of-field to pixels
        final double depth = settings.depthOfField / settings.pixelPitch;
        for (PeakResult r : results.getResults()) {
            final double precision = r.getPrecision(settings.pixelPitch, gain, emCCD);
            stats[PRECISION].add(precision);
            // The error stores the z-depth in pixels
            if (Math.abs(r.error) < depth)
                stats[PRECISION_IN_FOCUS].add(precision);
            stats[WIDTH].add(r.getSD());
        }
        // Compute density per frame. Multithread for speed
        if (settings.densityRadius > 0) {
            IJ.showStatus("Calculating density ...");
            ExecutorService threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
            List<Future<?>> futures = new LinkedList<Future<?>>();
            final ArrayList<float[]> coords = new ArrayList<float[]>();
            int t = results.getHead().getFrame();
            final Statistics densityStats = stats[DENSITY];
            final float radius = (float) (settings.densityRadius * getHWHM());
            final Rectangle bounds = results.getBounds();
            currentIndex = 0;
            finalIndex = results.getTail().getFrame();
            // Store the density for each result.
            int[] allDensity = new int[results.size()];
            int allIndex = 0;
            for (PeakResult r : results.getResults()) {
                if (t != r.getFrame()) {
                    allIndex += runDensityCalculation(threadPool, futures, coords, densityStats, radius, bounds, allDensity, allIndex);
                }
                coords.add(new float[] { r.getXPosition(), r.getYPosition() });
                t = r.getFrame();
            }
            runDensityCalculation(threadPool, futures, coords, densityStats, radius, bounds, allDensity, allIndex);
            Utils.waitForCompletion(futures);
            threadPool.shutdownNow();
            threadPool = null;
            IJ.showProgress(1);
            // Split results into singles (density = 0) and clustered (density > 0)
            MemoryPeakResults singles = copyMemoryPeakResults("No Density");
            MemoryPeakResults clustered = copyMemoryPeakResults("Density");
            int i = 0;
            for (PeakResult r : results.getResults()) {
                // Store density in the original value field
                r.origValue = allDensity[i];
                if (allDensity[i++] == 0)
                    singles.add(r);
                else
                    clustered.add(r);
            }
        }
    }
    StringBuilder sb = new StringBuilder();
    sb.append(datasetNumber).append("\t");
    sb.append((fluorophores == null) ? localisations.size() : fluorophores.size()).append("\t");
    sb.append(stats[SAMPLED_BLINKS].getN() + (int) stats[SAMPLED_BLINKS].getSum()).append("\t");
    sb.append(localisations.size()).append("\t");
    sb.append(nFrames).append("\t");
    sb.append(Utils.rounded(areaInUm)).append("\t");
    sb.append(Utils.rounded(localisations.size() / (areaInUm * nFrames), 4)).append("\t");
    sb.append(Utils.rounded(getHWHM(), 4)).append("\t");
    double s = getPsfSD();
    sb.append(Utils.rounded(s, 4)).append("\t");
    s *= settings.pixelPitch;
    final double sa = PSFCalculator.squarePixelAdjustment(s, settings.pixelPitch) / settings.pixelPitch;
    sb.append(Utils.rounded(sa, 4)).append("\t");
    // Width not valid for the Image PSF
    int nStats = (imagePSF) ? stats.length - 1 : stats.length;
    for (int i = 0; i < nStats; i++) {
        double centre = (alwaysRemoveOutliers[i]) ? ((StoredDataStatistics) stats[i]).getStatistics().getPercentile(50) : stats[i].getMean();
        sb.append(Utils.rounded(centre, 4)).append("\t");
    }
    if (java.awt.GraphicsEnvironment.isHeadless()) {
        IJ.log(sb.toString());
        return stats[SIGNAL].getMean();
    } else {
        summaryTable.append(sb.toString());
    }
    // Show histograms
    if (settings.showHistograms) {
        IJ.showStatus("Calculating histograms ...");
        boolean[] chosenHistograms = getChoosenHistograms();
        WindowOrganiser wo = new WindowOrganiser();
        boolean requireRetile = false;
        for (int i = 0; i < NAMES.length; i++) {
            if (chosenHistograms[i]) {
                wo.add(Utils.showHistogram(TITLE, (StoredDataStatistics) stats[i], NAMES[i], (integerDisplay[i]) ? 1 : 0, (settings.removeOutliers || alwaysRemoveOutliers[i]) ? 2 : 0, settings.histogramBins * ((integerDisplay[i]) ? 100 : 1)));
                requireRetile = requireRetile || Utils.isNewWindow();
            }
        }
        wo.tile();
    }
    IJ.showStatus("");
    return stats[SIGNAL].getMean();
}

Example 9 with LocalisationModel

use of gdsc.smlm.model.LocalisationModel 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
	 * @param results
	 * @param title
	 */
private void savePulses(List<LocalisationModel> localisations, MemoryPeakResults results, String title) {
    sortLocalisationsByIdThenTime(localisations);
    MemoryPeakResults traceResults = copyMemoryPeakResults("Pulses");
    LocalisationModel start = null;
    int currentId = -1;
    int n = 0;
    float[] params = new float[7];
    double noise = 0;
    int lastT = -1;
    for (LocalisationModel localisation : localisations) {
        if (currentId != localisation.getId() || lastT + 1 != localisation.getTime()) {
            if (n > 0) {
                params[Gaussian2DFunction.BACKGROUND] /= n;
                params[Gaussian2DFunction.X_POSITION] /= n;
                params[Gaussian2DFunction.Y_POSITION] /= n;
                params[Gaussian2DFunction.X_SD] /= n;
                params[Gaussian2DFunction.Y_SD] /= n;
                ExtendedPeakResult p = new ExtendedPeakResult(start.getTime(), (int) Math.round(start.getX()), (int) Math.round(start.getY()), 0, 0, (float) (Math.sqrt(noise)), 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();
            n = 0;
            params = new float[7];
            noise = 0;
        }
        final double[] data = localisation.getData();
        params[Gaussian2DFunction.BACKGROUND] += data[0];
        params[Gaussian2DFunction.X_POSITION] += localisation.getX();
        params[Gaussian2DFunction.Y_POSITION] += localisation.getY();
        params[Gaussian2DFunction.SIGNAL] += localisation.getIntensity();
        noise += data[1] * data[1];
        params[Gaussian2DFunction.X_SD] += data[2];
        params[Gaussian2DFunction.Y_SD] += data[3];
        n++;
        lastT = localisation.getTime();
    }
    // Final pulse
    if (n > 0) {
        params[Gaussian2DFunction.BACKGROUND] /= n;
        params[Gaussian2DFunction.X_POSITION] /= n;
        params[Gaussian2DFunction.Y_POSITION] /= n;
        params[Gaussian2DFunction.X_SD] /= n;
        params[Gaussian2DFunction.Y_SD] /= n;
        traceResults.add(new ExtendedPeakResult(start.getTime(), (int) Math.round(start.getX()), (int) Math.round(start.getY()), 0, 0, (float) (Math.sqrt(noise)), params, null, lastT, currentId));
    }
    traceResults.end();
}

Example 10 with LocalisationModel

use of gdsc.smlm.model.LocalisationModel in project GDSC-SMLM by aherbert.

the class CreateData method getIds.

private int[] getIds(List<LocalisationModel> localisations) {
    int[] ids = new int[localisations.size()];
    if (localisations.isEmpty())
        return ids;
    int count = 0;
    int id = localisations.get(0).getId();
    ids[count++] = id;
    for (LocalisationModel l : localisations) {
        if (id != l.getId()) {
            id = l.getId();
            ids[count++] = id;
        }
    }
    return Arrays.copyOf(ids, count);
}