Search in sources :

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;
}
Also used : Rectangle(java.awt.Rectangle) MemoryPeakResults(gdsc.smlm.results.MemoryPeakResults) ImagePSFModel(gdsc.smlm.model.ImagePSFModel) ImageStack(ij.ImageStack) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) SummaryStatistics(org.apache.commons.math3.stat.descriptive.SummaryStatistics) Calibration(gdsc.smlm.results.Calibration) ImagePlus(ij.ImagePlus) LinkedList(java.util.LinkedList) IJImageSource(gdsc.smlm.ij.IJImageSource) LocalisationModel(gdsc.smlm.model.LocalisationModel) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future) LocalisationModelSet(gdsc.smlm.model.LocalisationModelSet)

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;
}
Also used : LocalisationModel(gdsc.smlm.model.LocalisationModel) ArrayList(java.util.ArrayList) LocalisationModelSet(gdsc.smlm.model.LocalisationModelSet) TIntHashSet(gnu.trove.set.hash.TIntHashSet)

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;
}
Also used : SpatialDistribution(gdsc.smlm.model.SpatialDistribution) ArrayList(java.util.ArrayList) LocalisationModelSet(gdsc.smlm.model.LocalisationModelSet)

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");
}
Also used : PoissonDistribution(org.apache.commons.math3.distribution.PoissonDistribution) ActivationEnergyImageModel(gdsc.smlm.model.ActivationEnergyImageModel) CompoundMoleculeModel(gdsc.smlm.model.CompoundMoleculeModel) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) Random(gdsc.core.utils.Random) GenericDialog(ij.gui.GenericDialog) SpatialIllumination(gdsc.smlm.model.SpatialIllumination) SpatialDistribution(gdsc.smlm.model.SpatialDistribution) StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) LocalisationModel(gdsc.smlm.model.LocalisationModel) FluorophoreSequenceModel(gdsc.smlm.model.FluorophoreSequenceModel) FixedLifetimeImageModel(gdsc.smlm.model.FixedLifetimeImageModel) LocalisationModelSet(gdsc.smlm.model.LocalisationModelSet) ImageModel(gdsc.smlm.model.ImageModel) ActivationEnergyImageModel(gdsc.smlm.model.ActivationEnergyImageModel) FixedLifetimeImageModel(gdsc.smlm.model.FixedLifetimeImageModel)

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;
    }
}
Also used : LocalisationModel(gdsc.smlm.model.LocalisationModel) PSFOffset(gdsc.smlm.ij.settings.PSFOffset) LocalisationModelSet(gdsc.smlm.model.LocalisationModelSet) ImagePSFModel(gdsc.smlm.model.ImagePSFModel) ImagePlus(ij.ImagePlus) PSFSettings(gdsc.smlm.ij.settings.PSFSettings) IOException(java.io.IOException) NullArgumentException(org.apache.commons.math3.exception.NullArgumentException)

Aggregations

LocalisationModelSet (gdsc.smlm.model.LocalisationModelSet)5 LocalisationModel (gdsc.smlm.model.LocalisationModel)4 ImagePSFModel (gdsc.smlm.model.ImagePSFModel)2 SpatialDistribution (gdsc.smlm.model.SpatialDistribution)2 ImagePlus (ij.ImagePlus)2 ArrayList (java.util.ArrayList)2 Random (gdsc.core.utils.Random)1 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)1 IJImageSource (gdsc.smlm.ij.IJImageSource)1 PSFOffset (gdsc.smlm.ij.settings.PSFOffset)1 PSFSettings (gdsc.smlm.ij.settings.PSFSettings)1 ActivationEnergyImageModel (gdsc.smlm.model.ActivationEnergyImageModel)1 CompoundMoleculeModel (gdsc.smlm.model.CompoundMoleculeModel)1 FixedLifetimeImageModel (gdsc.smlm.model.FixedLifetimeImageModel)1 FluorophoreSequenceModel (gdsc.smlm.model.FluorophoreSequenceModel)1 ImageModel (gdsc.smlm.model.ImageModel)1 SpatialIllumination (gdsc.smlm.model.SpatialIllumination)1 Calibration (gdsc.smlm.results.Calibration)1 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)1 TIntHashSet (gnu.trove.set.hash.TIntHashSet)1