Search in sources :

Example 1 with ActivationEnergyImageModel

use of gdsc.smlm.model.ActivationEnergyImageModel 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 2 with ActivationEnergyImageModel

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

the class BlinkEstimatorTest method estimateBlinking.

private TIntHashSet estimateBlinking(double nBlinks, double tOn, double tOff, int particles, double fixedFraction, boolean timeAtLowerBound, boolean doAssert) {
    SpatialIllumination activationIllumination = new UniformIllumination(100);
    int totalSteps = 100;
    double eAct = totalSteps * 0.3 * activationIllumination.getAveragePhotons();
    ImageModel imageModel = new ActivationEnergyImageModel(eAct, activationIllumination, tOn, 0, tOff, 0, nBlinks);
    imageModel.setRandomGenerator(rand);
    double[] max = new double[] { 256, 256, 32 };
    double[] min = new double[3];
    SpatialDistribution distribution = new UniformDistribution(min, max, rand.nextInt());
    List<CompoundMoleculeModel> compounds = new ArrayList<CompoundMoleculeModel>(1);
    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);
    List<CompoundMoleculeModel> molecules = imageModel.createMolecules(compounds, particles, distribution, false);
    // Activate fluorophores
    List<? extends FluorophoreSequenceModel> fluorophores = imageModel.createFluorophores(molecules, totalSteps);
    totalSteps = checkTotalSteps(totalSteps, fluorophores);
    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;
    //		System.out.printf("Photons = %f\n", photonStats.getMean());
    //
    //		List<FluorophoreSequenceModel> newFluorophores = new ArrayList<FluorophoreSequenceModel>(fluorophores.size());
    //		for (FluorophoreSequenceModel f : fluorophores)
    //		{
    //			if (id[f.getId()])
    //				newFluorophores.add(f);
    //		}
    //		fluorophores = newFluorophores;
    MemoryPeakResults results = new MemoryPeakResults();
    results.setCalibration(new Calibration(pixelPitch, 1, msPerFrame));
    for (LocalisationModel l : localisations) {
        // Remove by intensity threshold and optionally at random.
        if (l.getIntensity() < minPhotons || rand.nextDouble() < pDelete)
            continue;
        float[] params = new float[7];
        params[Gaussian2DFunction.X_POSITION] = (float) l.getX();
        params[Gaussian2DFunction.Y_POSITION] = (float) l.getY();
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        params[Gaussian2DFunction.SIGNAL] = (float) (l.getIntensity());
        results.addf(l.getTime(), 0, 0, 0, 0, 0, params, null);
    }
    // Add random localisations
    for (int i = (int) (localisations.size() * pAdd); i-- > 0; ) {
        float[] params = new float[7];
        params[Gaussian2DFunction.X_POSITION] = (float) (rand.nextDouble() * max[0]);
        params[Gaussian2DFunction.Y_POSITION] = (float) (rand.nextDouble() * max[1]);
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        // Intensity doesn't matter at the moment for tracing
        params[Gaussian2DFunction.SIGNAL] = (float) (photons);
        results.addf(1 + rand.nextInt(totalSteps), 0, 0, 0, 0, 0, params, null);
    }
    // Get actual simulated stats ...
    Statistics statsNBlinks = new Statistics();
    Statistics statsTOn = new Statistics();
    Statistics statsTOff = new Statistics();
    Statistics statsSampledNBlinks = new Statistics();
    Statistics statsSampledTOn = new Statistics();
    StoredDataStatistics statsSampledTOff = new StoredDataStatistics();
    for (FluorophoreSequenceModel f : fluorophores) {
        statsNBlinks.add(f.getNumberOfBlinks());
        statsTOn.add(f.getOnTimes());
        statsTOff.add(f.getOffTimes());
        int[] on = f.getSampledOnTimes();
        statsSampledNBlinks.add(on.length);
        statsSampledTOn.add(on);
        statsSampledTOff.add(f.getSampledOffTimes());
    }
    System.out.printf("N = %d (%d), N-blinks = %f, tOn = %f, tOff = %f, Fixed = %f\n", fluorophores.size(), localisations.size(), nBlinks, tOn, tOff, fixedFraction);
    System.out.printf("Actual N-blinks = %f (%f), tOn = %f (%f), tOff = %f (%f), 95%% = %f, max = %f\n", statsNBlinks.getMean(), statsSampledNBlinks.getMean(), statsTOn.getMean(), statsSampledTOn.getMean(), statsTOff.getMean(), statsSampledTOff.getMean(), statsSampledTOff.getStatistics().getPercentile(95), statsSampledTOff.getStatistics().getMax());
    System.out.printf("-=-=--=-\n");
    BlinkEstimator be = new BlinkEstimator();
    be.maxDarkTime = (int) (tOff * 10);
    be.msPerFrame = msPerFrame;
    be.relativeDistance = false;
    double d = ImageModel.getRandomMoveDistance(diffusionRate);
    be.searchDistance = (fixedFraction < 1) ? Math.sqrt(2 * d * d) * 3 : 0;
    be.timeAtLowerBound = timeAtLowerBound;
    be.showPlots = false;
    //Assert.assertTrue("Max dark time must exceed the dark time of the data (otherwise no plateau)",
    //		be.maxDarkTime > statsSampledTOff.getStatistics().getMax());
    int nMolecules = fluorophores.size();
    if (usePopulationStatistics) {
        nBlinks = statsNBlinks.getMean();
        tOff = statsTOff.getMean();
    } else {
        nBlinks = statsSampledNBlinks.getMean();
        tOff = statsSampledTOff.getMean();
    }
    // See if any fitting regime gets a correct answer
    TIntHashSet ok = new TIntHashSet();
    for (int nFittedPoints = MIN_FITTED_POINTS; nFittedPoints <= MAX_FITTED_POINTS; nFittedPoints++) {
        be.nFittedPoints = nFittedPoints;
        be.computeBlinkingRate(results, true);
        double moleculesError = DoubleEquality.relativeError(nMolecules, be.getNMolecules());
        double blinksError = DoubleEquality.relativeError(nBlinks, be.getNBlinks());
        double offError = DoubleEquality.relativeError(tOff * msPerFrame, be.getTOff());
        System.out.printf("Error %d: N = %f, blinks = %f, tOff = %f : %f\n", nFittedPoints, moleculesError, blinksError, offError, (moleculesError + blinksError + offError) / 3);
        if (moleculesError < relativeError && blinksError < relativeError && offError < relativeError) {
            ok.add(nFittedPoints);
            System.out.printf("-=-=--=-\n");
            System.out.printf("*** Correct at %d fitted points ***\n", nFittedPoints);
            if (doAssert)
                break;
        }
    //if (!be.isIncreaseNFittedPoints())
    //	break;
    }
    System.out.printf("-=-=--=-\n");
    if (doAssert)
        Assert.assertFalse(ok.isEmpty());
    //Assert.assertEquals("Invalid t-off", tOff * msPerFrame, be.getTOff(), tOff * msPerFrame * relativeError);
    return ok;
}
Also used : ActivationEnergyImageModel(gdsc.smlm.model.ActivationEnergyImageModel) CompoundMoleculeModel(gdsc.smlm.model.CompoundMoleculeModel) ArrayList(java.util.ArrayList) TIntHashSet(gnu.trove.set.hash.TIntHashSet) CompoundMoleculeModel(gdsc.smlm.model.CompoundMoleculeModel) MoleculeModel(gdsc.smlm.model.MoleculeModel) SpatialIllumination(gdsc.smlm.model.SpatialIllumination) MemoryPeakResults(gdsc.smlm.results.MemoryPeakResults) SpatialDistribution(gdsc.smlm.model.SpatialDistribution) UniformDistribution(gdsc.smlm.model.UniformDistribution) StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) Calibration(gdsc.smlm.results.Calibration) StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) Statistics(gdsc.core.utils.Statistics) UniformIllumination(gdsc.smlm.model.UniformIllumination) LocalisationModel(gdsc.smlm.model.LocalisationModel) FluorophoreSequenceModel(gdsc.smlm.model.FluorophoreSequenceModel) ImageModel(gdsc.smlm.model.ImageModel) ActivationEnergyImageModel(gdsc.smlm.model.ActivationEnergyImageModel)

Aggregations

StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)2 ActivationEnergyImageModel (gdsc.smlm.model.ActivationEnergyImageModel)2 CompoundMoleculeModel (gdsc.smlm.model.CompoundMoleculeModel)2 FluorophoreSequenceModel (gdsc.smlm.model.FluorophoreSequenceModel)2 ImageModel (gdsc.smlm.model.ImageModel)2 LocalisationModel (gdsc.smlm.model.LocalisationModel)2 SpatialDistribution (gdsc.smlm.model.SpatialDistribution)2 SpatialIllumination (gdsc.smlm.model.SpatialIllumination)2 Random (gdsc.core.utils.Random)1 Statistics (gdsc.core.utils.Statistics)1 FixedLifetimeImageModel (gdsc.smlm.model.FixedLifetimeImageModel)1 LocalisationModelSet (gdsc.smlm.model.LocalisationModelSet)1 MoleculeModel (gdsc.smlm.model.MoleculeModel)1 UniformDistribution (gdsc.smlm.model.UniformDistribution)1 UniformIllumination (gdsc.smlm.model.UniformIllumination)1 Calibration (gdsc.smlm.results.Calibration)1 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)1 TIntHashSet (gnu.trove.set.hash.TIntHashSet)1 GenericDialog (ij.gui.GenericDialog)1 ArrayList (java.util.ArrayList)1