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

Example 1 with LocalisationModelSet

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

the class CreateData method drawImage.

//StoredDataStatistics rawPhotons = new StoredDataStatistics();
//StoredDataStatistics drawPhotons = new StoredDataStatistics();
//	private synchronized void addRaw(double d)
//	{
//		//rawPhotons.add(d);
//	}
//
//	private synchronized void addDraw(double d)
//	{
//		//drawPhotons.add(d);
//	}
/**
	 * Create an image from the localisations using the configured PSF width. Draws a new stack
	 * image.
	 * <p>
	 * Note that the localisations are filtered using the signal. The input list of localisations will be updated.
	 * 
	 * @param localisationSets
	 * @return The localisations
	 */
private List<LocalisationModel> drawImage(final List<LocalisationModelSet> localisationSets) {
    if (localisationSets.isEmpty())
        return null;
    // Create a new list for all localisation that are drawn (i.e. pass the signal filters)
    List<LocalisationModelSet> newLocalisations = Collections.synchronizedList(new ArrayList<LocalisationModelSet>(localisationSets.size()));
    photonsRemoved = new AtomicInteger();
    t1Removed = new AtomicInteger();
    tNRemoved = new AtomicInteger();
    photonStats = new SummaryStatistics();
    // Add drawn spots to memory
    results = new MemoryPeakResults();
    Calibration c = new Calibration(settings.pixelPitch, settings.getTotalGain(), settings.exposureTime);
    c.setEmCCD((settings.getEmGain() > 1));
    c.setBias(settings.bias);
    c.setReadNoise(settings.readNoise * ((settings.getCameraGain() > 0) ? settings.getCameraGain() : 1));
    c.setAmplification(settings.getAmplification());
    results.setCalibration(c);
    results.setSortAfterEnd(true);
    results.begin();
    maxT = localisationSets.get(localisationSets.size() - 1).getTime();
    // Display image
    ImageStack stack = new ImageStack(settings.size, settings.size, maxT);
    final double psfSD = getPsfSD();
    if (psfSD <= 0)
        return null;
    ImagePSFModel imagePSFModel = null;
    if (imagePSF) {
        // Create one Image PSF model that can be copied
        imagePSFModel = createImagePSF(localisationSets);
        if (imagePSFModel == null)
            return null;
    }
    IJ.showStatus("Drawing image ...");
    // Multi-thread for speed
    // Note that the default Executors.newCachedThreadPool() will continue to make threads if
    // new tasks are added. We need to limit the tasks that can be added using a fixed size
    // blocking queue.
    // http://stackoverflow.com/questions/1800317/impossible-to-make-a-cached-thread-pool-with-a-size-limit
    // ExecutorService threadPool = Executors.newCachedThreadPool();
    ExecutorService threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
    List<Future<?>> futures = new LinkedList<Future<?>>();
    // Count all the frames to process
    frame = 0;
    totalFrames = maxT;
    // Collect statistics on the number of photons actually simulated
    // Process all frames
    int i = 0;
    int lastT = -1;
    for (LocalisationModelSet l : localisationSets) {
        if (Utils.isInterrupted())
            break;
        if (l.getTime() != lastT) {
            lastT = l.getTime();
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, i, lastT, createPSFModel(imagePSFModel), results, stack, poissonNoise, new RandomDataGenerator(createRandomGenerator()))));
        }
        i++;
    }
    // Finish processing data
    Utils.waitForCompletion(futures);
    futures.clear();
    if (Utils.isInterrupted()) {
        IJ.showProgress(1);
        return null;
    }
    // Do all the frames that had no localisations
    for (int t = 1; t <= maxT; t++) {
        if (Utils.isInterrupted())
            break;
        if (stack.getPixels(t) == null) {
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, maxT, t, null, results, stack, poissonNoise, new RandomDataGenerator(createRandomGenerator()))));
        }
    }
    // Finish
    Utils.waitForCompletion(futures);
    threadPool.shutdown();
    IJ.showProgress(1);
    if (Utils.isInterrupted()) {
        return null;
    }
    results.end();
    // Clear memory
    imagePSFModel = null;
    threadPool = null;
    futures.clear();
    futures = null;
    if (photonsRemoved.get() > 0)
        Utils.log("Removed %d localisations with less than %.1f rendered photons", photonsRemoved.get(), settings.minPhotons);
    if (t1Removed.get() > 0)
        Utils.log("Removed %d localisations with no neighbours @ SNR %.2f", t1Removed.get(), settings.minSNRt1);
    if (tNRemoved.get() > 0)
        Utils.log("Removed %d localisations with valid neighbours @ SNR %.2f", tNRemoved.get(), settings.minSNRtN);
    if (photonStats.getN() > 0)
        Utils.log("Average photons rendered = %s +/- %s", Utils.rounded(photonStats.getMean()), Utils.rounded(photonStats.getStandardDeviation()));
    //System.out.printf("rawPhotons = %f\n", rawPhotons.getMean());
    //System.out.printf("drawPhotons = %f\n", drawPhotons.getMean());
    //Utils.showHistogram("draw photons", drawPhotons, "photons", true, 0, 1000);
    // Update with all those localisation that have been drawn
    localisationSets.clear();
    localisationSets.addAll(newLocalisations);
    newLocalisations = null;
    IJ.showStatus("Displaying image ...");
    ImageStack newStack = stack;
    if (!settings.rawImage) {
        // Get the global limits and ensure all values can be represented
        Object[] imageArray = stack.getImageArray();
        float[] limits = Maths.limits((float[]) imageArray[0]);
        for (int j = 1; j < imageArray.length; j++) limits = Maths.limits(limits, (float[]) imageArray[j]);
        // Leave bias in place
        limits[0] = 0;
        // Check if the image will fit in a 16-bit range
        if ((limits[1] - limits[0]) < 65535) {
            // Convert to 16-bit
            newStack = new ImageStack(stack.getWidth(), stack.getHeight(), stack.getSize());
            // Account for rounding
            final float min = (float) (limits[0] - 0.5);
            for (int j = 0; j < imageArray.length; j++) {
                float[] image = (float[]) imageArray[j];
                short[] pixels = new short[image.length];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = (short) (image[k] - min);
                }
                newStack.setPixels(pixels, j + 1);
                // Free memory
                imageArray[j] = null;
                // Attempt to stay within memory (check vs 32MB)
                if (MemoryPeakResults.freeMemory() < 33554432L)
                    MemoryPeakResults.runGCOnce();
            }
        } else {
            // Keep as 32-bit but round to whole numbers
            for (int j = 0; j < imageArray.length; j++) {
                float[] pixels = (float[]) imageArray[j];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = Math.round(pixels[k]);
                }
            }
        }
    }
    // Show image
    ImagePlus imp = Utils.display(CREATE_DATA_IMAGE_TITLE, newStack);
    ij.measure.Calibration cal = new ij.measure.Calibration();
    String unit = "nm";
    double unitPerPixel = settings.pixelPitch;
    if (unitPerPixel > 100) {
        unit = "um";
        unitPerPixel /= 1000.0;
    }
    cal.setUnit(unit);
    cal.pixelHeight = cal.pixelWidth = unitPerPixel;
    imp.setCalibration(cal);
    imp.setDimensions(1, 1, newStack.getSize());
    imp.resetDisplayRange();
    imp.updateAndDraw();
    saveImage(imp);
    results.setSource(new IJImageSource(imp));
    results.setName(CREATE_DATA_IMAGE_TITLE + " (" + TITLE + ")");
    results.setConfiguration(createConfiguration((float) psfSD));
    results.setBounds(new Rectangle(0, 0, settings.size, settings.size));
    MemoryPeakResults.addResults(results);
    setBenchmarkResults(imp, results);
    if (benchmarkMode && benchmarkParameters != null)
        benchmarkParameters.setPhotons(results);
    List<LocalisationModel> localisations = toLocalisations(localisationSets);
    savePulses(localisations, results, CREATE_DATA_IMAGE_TITLE);
    // Saved the fixed and moving localisations into different datasets
    saveFixedAndMoving(results, CREATE_DATA_IMAGE_TITLE);
    return localisations;
}

Example 2 with LocalisationModelSet

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

the class CreateData method combineSimulationSteps.

private List<LocalisationModelSet> combineSimulationSteps(List<LocalisationModel> localisations) {
    // Allow fractional integration steps
    final double simulationStepsPerFrame = (settings.stepsPerSecond * settings.exposureTime) / 1000.0;
    List<LocalisationModelSet> newLocalisations = new ArrayList<LocalisationModelSet>((int) (localisations.size() / simulationStepsPerFrame));
    //System.out.printf("combineSimulationSteps @ %f\n", simulationStepsPerFrame);
    final double gain = (settings.getTotalGain() > 0) ? settings.getTotalGain() : 1;
    sortLocalisationsByIdThenTime(localisations);
    int[] idList = getIds(localisations);
    movingMolecules = new TIntHashSet(idList.length);
    int index = 0;
    for (int id : idList) {
        int fromIndex = findIndexById(localisations, index, id);
        if (fromIndex > -1) {
            int toIndex = findLastIndexById(localisations, fromIndex, id);
            List<LocalisationModel> subset = localisations.subList(fromIndex, toIndex + 1);
            index = toIndex;
            // Store the IDs of any moving molecules
            if (isMoving(subset))
                movingMolecules.add(id);
            // The frames may be longer or shorter than the simulation steps. Allocate the step
            // proportionately to each frame it overlaps:
            //
            // Steps:  |-- 0 --|-- 1 --|-- 2 --|--
            // Frames: |--- 0 ---|--- 1 ---|--- 2 ---|
            //
            //         ^       ^
            //         |       |
            //         |       End frame
            //         |
            //         Start frame
            final double firstFrame = getStartFrame(subset.get(0), simulationStepsPerFrame);
            final double lastFrame = getEndFrame(subset.get(subset.size() - 1), simulationStepsPerFrame);
            // Get the first frame offset and allocate space to store all potential frames  
            final int intFirstFrame = (int) firstFrame;
            final int intLastFrame = (int) Math.ceil(lastFrame);
            LocalisationModelSet[] sets = new LocalisationModelSet[intLastFrame - intFirstFrame + 1];
            // Process each step
            for (LocalisationModel l : subset) {
                // Get the fractional start and end frames 
                double startFrame = getStartFrame(l, simulationStepsPerFrame);
                double endFrame = getEndFrame(l, simulationStepsPerFrame);
                // Round down to get the actual frames that are overlapped
                int start = (int) startFrame;
                int end = (int) endFrame;
                // Check if the span covers a fraction of the end frame, otherwise decrement to ignore that frame
                if (end > start && endFrame == end) {
                    // E.g. convert 
                    // Steps:      |-- 0 --|
                    // Frames: |- 0 -|- 1 -|- 2 -|
                    // to
                    // Steps:      |-- 0 --|
                    // Frames: |- 0 -|- 1 -|
                    end--;
                }
                if (start == end) {
                    // If the step falls within one frame then add it to the set
                    int tIndex = start - intFirstFrame;
                    if (sets[tIndex] == null)
                        sets[tIndex] = new LocalisationModelSet(id, start);
                    sets[tIndex].add(l);
                } else {
                    // Add the localisation to all the frames that the step spans
                    final double total = endFrame - startFrame;
                    //double t = 0;
                    for (int frame = start; frame <= end; frame++) {
                        // Get the fraction to allocate to this frame
                        double fraction;
                        int state = (l.isContinuous() ? LocalisationModel.CONTINUOUS : 0);
                        if (frame == start) {
                            state |= LocalisationModel.NEXT | (l.hasPrevious() ? LocalisationModel.PREVIOUS : 0);
                            // start
                            if (startFrame == start)
                                fraction = 1;
                            else
                                fraction = (Math.ceil(startFrame) - startFrame);
                        } else if (frame == end) {
                            state |= LocalisationModel.PREVIOUS | (l.hasNext() ? LocalisationModel.NEXT : 0);
                            // |=====|----|
                            // |     |     
                            // |     endFrame
                            // end
                            fraction = (endFrame - end);
                        } else {
                            state |= LocalisationModel.CONTINUOUS;
                            fraction = 1;
                        }
                        //t += fraction;
                        // Add to the set
                        int tIndex = frame - intFirstFrame;
                        if (sets[tIndex] == null)
                            sets[tIndex] = new LocalisationModelSet(id, frame);
                        sets[tIndex].add(getFraction(l, fraction / total, state));
                    }
                //if (t < total * 0.98)
                //{
                //	System.out.printf("Total error %g < %g : %f (%d) -> %f (%d)\n", t, total, startFrame,
                //			start, endFrame, end);
                //}
                }
            }
            LocalisationModelSet previous = null;
            for (int i = 0; i < sets.length; i++) {
                if (sets[i] != null) {
                    sets[i].setPrevious(previous);
                    // Create a data array and store the current intensity after gain. 
                    // This is used later to filter based on SNR
                    sets[i].setData(new double[] { 0, 0, 0, 0, sets[i].getIntensity() * gain });
                    newLocalisations.add(sets[i]);
                }
                previous = sets[i];
            }
        }
    }
    // Sort by time
    Collections.sort(newLocalisations);
    return newLocalisations;
}

Example 3 with LocalisationModelSet

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

the class CreateData method filterToImageBounds.

/**
	 * Filter those not in the distribution
	 * 
	 * @param localisationSets
	 * @return
	 */
private List<LocalisationModelSet> filterToImageBounds(List<LocalisationModelSet> localisationSets) {
    List<LocalisationModelSet> newLocalisations = new ArrayList<LocalisationModelSet>(localisationSets.size());
    SpatialDistribution bounds = createUniformDistribution(0);
    for (LocalisationModelSet s : localisationSets) {
        if (bounds.isWithinXY(s.toLocalisation().getCoordinates()))
            newLocalisations.add(s);
    }
    return newLocalisations;
}

Example 4 with LocalisationModelSet

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

the class CreateData method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    extraOptions = Utils.isExtraOptions();
    simpleMode = (arg != null && arg.contains("simple"));
    benchmarkMode = (arg != null && arg.contains("benchmark"));
    spotMode = (arg != null && arg.contains("spot"));
    trackMode = (arg != null && arg.contains("track"));
    if ("load".equals(arg)) {
        loadBenchmarkData();
        return;
    }
    // Each localisation is a simulated emission of light from a point in space and time
    List<LocalisationModel> localisations = null;
    // Each localisation set is a collection of localisations that represent all localisations
    // with the same ID that are on in the same image time frame (Note: the simulation
    // can create many localisations per fluorophore per time frame which is useful when 
    // modelling moving particles)
    List<LocalisationModelSet> localisationSets = null;
    // Each fluorophore contains the on and off times when light was emitted 
    List<? extends FluorophoreSequenceModel> fluorophores = null;
    if (simpleMode || benchmarkMode || spotMode) {
        if (!showSimpleDialog())
            return;
        resetMemory();
        // 1 second frames
        settings.exposureTime = 1000;
        areaInUm = settings.size * settings.pixelPitch * settings.size * settings.pixelPitch / 1e6;
        // Number of spots per frame
        int n = 0;
        int[] nextN = null;
        SpatialDistribution dist;
        if (benchmarkMode) {
            // --------------------
            // BENCHMARK SIMULATION
            // --------------------
            // Draw the same point on the image repeatedly
            n = 1;
            dist = createFixedDistribution();
            reportAndSaveFittingLimits(dist);
        } else if (spotMode) {
            // ---------------
            // SPOT SIMULATION
            // ---------------
            // The spot simulation draws 0 or 1 random point per frame. 
            // Ensure we have 50% of the frames with a spot.
            nextN = new int[settings.particles * 2];
            Arrays.fill(nextN, 0, settings.particles, 1);
            Random rand = new Random();
            rand.shuffle(nextN);
            // Only put spots in the central part of the image
            double border = settings.size / 4.0;
            dist = createUniformDistribution(border);
        } else {
            // -----------------
            // SIMPLE SIMULATION
            // -----------------
            // The simple simulation draws n random points per frame to achieve a specified density.
            // No points will appear in multiple frames.
            // Each point has a random number of photons sampled from a range.
            // We can optionally use a mask. Create his first as it updates the areaInUm
            dist = createDistribution();
            // Randomly sample (i.e. not uniform density in all frames)
            if (settings.samplePerFrame) {
                final double mean = areaInUm * settings.density;
                System.out.printf("Mean samples = %f\n", mean);
                if (mean < 0.5) {
                    GenericDialog gd = new GenericDialog(TITLE);
                    gd.addMessage("The mean samples per frame is low: " + Utils.rounded(mean) + "\n \nContinue?");
                    gd.enableYesNoCancel();
                    gd.hideCancelButton();
                    gd.showDialog();
                    if (!gd.wasOKed())
                        return;
                }
                PoissonDistribution poisson = new PoissonDistribution(createRandomGenerator(), mean, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
                StoredDataStatistics samples = new StoredDataStatistics(settings.particles);
                while (samples.getSum() < settings.particles) {
                    samples.add(poisson.sample());
                }
                nextN = new int[samples.getN()];
                for (int i = 0; i < nextN.length; i++) nextN[i] = (int) samples.getValue(i);
            } else {
                // Use the density to get the number per frame
                n = (int) FastMath.max(1, Math.round(areaInUm * settings.density));
            }
        }
        RandomGenerator random = null;
        localisations = new ArrayList<LocalisationModel>(settings.particles);
        localisationSets = new ArrayList<LocalisationModelSet>(settings.particles);
        final int minPhotons = (int) settings.photonsPerSecond;
        final int range = (int) settings.photonsPerSecondMaximum - minPhotons + 1;
        if (range > 1)
            random = createRandomGenerator();
        // Add frames at the specified density until the number of particles has been reached
        int id = 0;
        int t = 0;
        while (id < settings.particles) {
            // Allow the number per frame to be specified
            if (nextN != null) {
                if (t >= nextN.length)
                    break;
                n = nextN[t];
            }
            // Simulate random positions in the frame for the specified density
            t++;
            for (int j = 0; j < n; j++) {
                final double[] xyz = dist.next();
                // Ignore within border. We do not want to draw things we cannot fit.
                //if (!distBorder.isWithinXY(xyz))
                //	continue;
                // Simulate random photons
                final int intensity = minPhotons + ((random != null) ? random.nextInt(range) : 0);
                LocalisationModel m = new LocalisationModel(id, t, xyz, intensity, LocalisationModel.CONTINUOUS);
                localisations.add(m);
                // Each localisation can be a separate localisation set
                LocalisationModelSet set = new LocalisationModelSet(id, t);
                set.add(m);
                localisationSets.add(set);
                id++;
            }
        }
    } else {
        if (!showDialog())
            return;
        resetMemory();
        areaInUm = settings.size * settings.pixelPitch * settings.size * settings.pixelPitch / 1e6;
        int totalSteps;
        double correlation = 0;
        ImageModel imageModel;
        if (trackMode) {
            // ----------------
            // TRACK SIMULATION
            // ----------------
            // In track mode we create fixed lifetime fluorophores that do not overlap in time.
            // This is the simplest simulation to test moving molecules.
            settings.seconds = (int) Math.ceil(settings.particles * (settings.exposureTime + settings.tOn) / 1000);
            totalSteps = 0;
            final double simulationStepsPerFrame = (settings.stepsPerSecond * settings.exposureTime) / 1000.0;
            imageModel = new FixedLifetimeImageModel(settings.stepsPerSecond * settings.tOn / 1000.0, simulationStepsPerFrame);
        } else {
            // ---------------
            // FULL SIMULATION
            // ---------------
            // The full simulation draws n random points in space.
            // The same molecule may appear in multiple frames, move and blink.
            //
            // Points are modelled as fluorophores that must be activated and then will 
            // blink and photo-bleach. The molecules may diffuse and this can be simulated 
            // with many steps per image frame. All steps from a frame are collected
            // into a localisation set which can be drawn on the output image.
            SpatialIllumination activationIllumination = createIllumination(settings.pulseRatio, settings.pulseInterval);
            // Generate additional frames so that each frame has the set number of simulation steps
            totalSteps = (int) Math.ceil(settings.seconds * settings.stepsPerSecond);
            // Since we have an exponential decay of activations
            // ensure half of the particles have activated by 30% of the frames.
            double eAct = totalSteps * 0.3 * activationIllumination.getAveragePhotons();
            // Q. Does tOn/tOff change depending on the illumination strength?
            imageModel = new ActivationEnergyImageModel(eAct, activationIllumination, settings.stepsPerSecond * settings.tOn / 1000.0, settings.stepsPerSecond * settings.tOffShort / 1000.0, settings.stepsPerSecond * settings.tOffLong / 1000.0, settings.nBlinksShort, settings.nBlinksLong);
            imageModel.setUseGeometricDistribution(settings.nBlinksGeometricDistribution);
            // Only use the correlation if selected for the distribution
            if (PHOTON_DISTRIBUTION[PHOTON_CORRELATED].equals(settings.photonDistribution))
                correlation = settings.correlation;
        }
        imageModel.setRandomGenerator(createRandomGenerator());
        imageModel.setPhotonBudgetPerFrame(true);
        imageModel.setDiffusion2D(settings.diffuse2D);
        imageModel.setRotation2D(settings.rotate2D);
        IJ.showStatus("Creating molecules ...");
        SpatialDistribution distribution = createDistribution();
        List<CompoundMoleculeModel> compounds = createCompoundMolecules();
        if (compounds == null)
            return;
        List<CompoundMoleculeModel> molecules = imageModel.createMolecules(compounds, settings.particles, distribution, settings.rotateInitialOrientation);
        // Activate fluorophores
        IJ.showStatus("Creating fluorophores ...");
        // Note: molecules list will be converted to compounds containing fluorophores
        fluorophores = imageModel.createFluorophores(molecules, totalSteps);
        if (fluorophores.isEmpty()) {
            IJ.error(TITLE, "No fluorophores created");
            return;
        }
        IJ.showStatus("Creating localisations ...");
        // TODO - Output a molecule Id for each fluorophore if using compound molecules. This allows analysis
        // of the ratio of trimers, dimers, monomers, etc that could be detected.
        totalSteps = checkTotalSteps(totalSteps, fluorophores);
        if (totalSteps == 0)
            return;
        imageModel.setPhotonDistribution(createPhotonDistribution());
        imageModel.setConfinementDistribution(createConfinementDistribution());
        // This should be optimised
        imageModel.setConfinementAttempts(10);
        localisations = imageModel.createImage(molecules, settings.fixedFraction, totalSteps, (double) settings.photonsPerSecond / settings.stepsPerSecond, correlation, settings.rotateDuringSimulation);
        // Re-adjust the fluorophores to the correct time
        if (settings.stepsPerSecond != 1) {
            final double scale = 1.0 / settings.stepsPerSecond;
            for (FluorophoreSequenceModel f : fluorophores) f.adjustTime(scale);
        }
        // Integrate the frames
        localisationSets = combineSimulationSteps(localisations);
        localisationSets = filterToImageBounds(localisationSets);
    }
    datasetNumber++;
    localisations = drawImage(localisationSets);
    if (localisations == null || localisations.isEmpty()) {
        IJ.error(TITLE, "No localisations created");
        return;
    }
    fluorophores = removeFilteredFluorophores(fluorophores, localisations);
    double signalPerFrame = showSummary(fluorophores, localisations);
    if (!benchmarkMode) {
        boolean fullSimulation = (!(simpleMode || spotMode));
        saveSimulationParameters(localisations.size(), fullSimulation, signalPerFrame);
    }
    IJ.showStatus("Saving data ...");
    //convertRelativeToAbsolute(molecules);
    saveFluorophores(fluorophores);
    saveImageResults(results);
    saveLocalisations(localisations);
    // The settings for the filenames may have changed
    SettingsManager.saveSettings(globalSettings);
    IJ.showStatus("Done");
}

Example 5 with LocalisationModelSet

use of gdsc.smlm.model.LocalisationModelSet 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;
    }
}