use of uk.ac.sussex.gdsc.smlm.model.MoleculeModel in project GDSC-SMLM by aherbert.
the class CreateData method createCompoundMolecules.
private List<CompoundMoleculeModel> createCompoundMolecules() {
// Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
final double diffusionFactor = (1000000.0 / (settings.getPixelPitch() * settings.getPixelPitch())) / settings.getStepsPerSecond();
List<CompoundMoleculeModel> compounds;
if (settings.getCompoundMolecules()) {
// Try and load the compounds from the text specification
try {
// Convert from the serialised objects to the compound model
final String text = settings.getCompoundText();
final Mixture.Builder builder = Mixture.newBuilder();
TextFormat.merge(text, builder);
compounds = new ArrayList<>(builder.getMoleculeCount());
int id = 1;
compoundNames = new ArrayList<>(builder.getMoleculeCount());
for (final Molecule m : builder.getMoleculeList()) {
final MoleculeModel[] molecules = new MoleculeModel[m.getAtomCount()];
for (int i = 0; i < molecules.length; i++) {
final AtomOrBuilder a = m.getAtomOrBuilder(i);
molecules[i] = new MoleculeModel(a.getMass(), a.getX(), a.getY(), a.getZ());
}
final CompoundMoleculeModel cm = new CompoundMoleculeModel(id++, 0, 0, 0, Arrays.asList(molecules));
cm.setFraction(m.getFraction());
cm.setDiffusionRate(m.getDiffusionRate() * diffusionFactor);
cm.setDiffusionType(DiffusionType.fromString(m.getDiffusionType()));
compounds.add(cm);
compoundNames.add(String.format("Fraction=%s, D=%s um^2/s", MathUtils.rounded(cm.getFraction()), MathUtils.rounded(m.getDiffusionRate())));
}
// Convert coordinates from nm to pixels
final double scaleFactor = 1.0 / settings.getPixelPitch();
for (final CompoundMoleculeModel c : compounds) {
c.scale(scaleFactor);
}
} catch (final IOException ex) {
IJ.error(TITLE, "Unable to create compound molecules");
return null;
}
} else {
// Create a simple compound with one molecule at the origin
compounds = new ArrayList<>(1);
final CompoundMoleculeModel m = new CompoundMoleculeModel(1, 0, 0, 0, Arrays.asList(new MoleculeModel(0, 0, 0, 0)));
m.setDiffusionRate(settings.getDiffusionRate() * diffusionFactor);
m.setDiffusionType(CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType()));
compounds.add(m);
}
return compounds;
}
use of uk.ac.sussex.gdsc.smlm.model.MoleculeModel in project GDSC-SMLM by aherbert.
the class CreateData method convertRelativeToAbsolute.
/**
* Update the fluorophores relative coordinates to absolute.
*
* @param molecules the molecules
*/
@SuppressWarnings("unused")
private static void convertRelativeToAbsolute(List<CompoundMoleculeModel> molecules) {
for (final CompoundMoleculeModel c : molecules) {
final double[] xyz = c.getCoordinates();
for (int n = c.getSize(); n-- > 0; ) {
final MoleculeModel m = c.getMolecule(n);
final double[] xyz2 = m.getCoordinates();
for (int i = 0; i < 3; i++) {
xyz2[i] += xyz[i];
}
}
}
}
use of uk.ac.sussex.gdsc.smlm.model.MoleculeModel in project GDSC-SMLM by aherbert.
the class BlinkEstimatorTest method estimateBlinking.
private TIntHashSet estimateBlinking(UniformRandomProvider rg, double blinkingRate, double ton, double toff, int particles, double fixedFraction, boolean timeAtLowerBound, boolean doAssert) {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MAXIMUM));
final SpatialIllumination activationIllumination = new UniformIllumination(100);
int totalSteps = 100;
final double eAct = totalSteps * 0.3 * activationIllumination.getAveragePhotons();
final ImageModel imageModel = new ActivationEnergyImageModel(eAct, activationIllumination, ton, 0, toff, 0, blinkingRate, rg);
final double[] max = new double[] { 256, 256, 32 };
final double[] min = new double[3];
final SpatialDistribution distribution = new UniformDistribution(min, max, rg.nextInt());
final List<CompoundMoleculeModel> compounds = new ArrayList<>(1);
final CompoundMoleculeModel c = new CompoundMoleculeModel(1, 0, 0, 0, Arrays.asList(new MoleculeModel(0, 0, 0, 0)));
c.setDiffusionRate(diffusionRate);
c.setDiffusionType(DiffusionType.RANDOM_WALK);
compounds.add(c);
final List<CompoundMoleculeModel> molecules = imageModel.createMolecules(compounds, particles, distribution, false);
// Activate fluorophores
final List<? extends FluorophoreSequenceModel> fluorophores = imageModel.createFluorophores(molecules, totalSteps);
totalSteps = checkTotalSteps(totalSteps, fluorophores);
final List<LocalisationModel> localisations = imageModel.createImage(molecules, fixedFraction, totalSteps, photons, 0.5, false);
// // Remove localisations to simulate missed counts.
// List<LocalisationModel> newLocalisations = new
// ArrayList<LocalisationModel>(localisations.size());
// boolean[] id = new boolean[fluorophores.size() + 1];
// Statistics photonStats = new Statistics();
// for (LocalisationModel l : localisations)
// {
// photonStats.add(l.getIntensity());
// // Remove by intensity threshold and optionally at random.
// if (l.getIntensity() < minPhotons || rand.nextDouble() < pDelete)
// continue;
// newLocalisations.add(l);
// id[l.getId()] = true;
// }
// localisations = newLocalisations;
// logger.info("Photons = %f", photonStats.getMean());
//
// List<FluorophoreSequenceModel> newFluorophores = new
// ArrayList<FluorophoreSequenceModel>(fluorophores.size());
// for (FluorophoreSequenceModel f : fluorophores)
// {
// if (id[f.getId()])
// newFluorophores.add(f);
// }
// fluorophores = newFluorophores;
final MemoryPeakResults results = new MemoryPeakResults();
final CalibrationWriter calibration = new CalibrationWriter();
calibration.setNmPerPixel(pixelPitch);
calibration.setExposureTime(msPerFrame);
calibration.setCountPerPhoton(1);
results.setCalibration(calibration.getCalibration());
results.setPsf(PsfHelper.create(PSFType.ONE_AXIS_GAUSSIAN_2D));
final float b = 0;
float intensity;
final float z = 0;
for (final LocalisationModel l : localisations) {
// Remove by intensity threshold and optionally at random.
if (l.getIntensity() < minPhotons || rg.nextDouble() < probabilityDelete) {
continue;
}
final int frame = l.getTime();
intensity = (float) l.getIntensity();
final float x = (float) l.getX();
final float y = (float) l.getY();
final float[] params = Gaussian2DPeakResultHelper.createParams(b, intensity, x, y, z, psfWidth);
results.add(frame, 0, 0, 0, 0, 0, 0, params, null);
}
// Add random localisations
// Intensity doesn't matter at the moment for tracing
intensity = (float) photons;
for (int i = (int) (localisations.size() * probabilityAdd); i-- > 0; ) {
final int frame = 1 + rg.nextInt(totalSteps);
final float x = (float) (rg.nextDouble() * max[0]);
final float y = (float) (rg.nextDouble() * max[1]);
final float[] params = Gaussian2DPeakResultHelper.createParams(b, intensity, x, y, z, psfWidth);
results.add(frame, 0, 0, 0, 0, 0, 0, params, null);
}
// Get actual simulated stats ...
final Statistics statsNBlinks = new Statistics();
final Statistics statsTOn = new Statistics();
final Statistics statsTOff = new Statistics();
final Statistics statsSampledNBlinks = new Statistics();
final Statistics statsSampledTOn = new Statistics();
final StoredDataStatistics statsSampledTOff = new StoredDataStatistics();
for (final FluorophoreSequenceModel f : fluorophores) {
statsNBlinks.add(f.getNumberOfBlinks());
statsTOn.add(f.getOnTimes());
statsTOff.add(f.getOffTimes());
final int[] on = f.getSampledOnTimes();
statsSampledNBlinks.add(on.length);
statsSampledTOn.add(on);
statsSampledTOff.add(f.getSampledOffTimes());
}
logger.info(FunctionUtils.getSupplier("N = %d (%d), N-blinks = %f, tOn = %f, tOff = %f, Fixed = %f", fluorophores.size(), localisations.size(), blinkingRate, ton, toff, fixedFraction));
logger.info(FunctionUtils.getSupplier("Actual N-blinks = %f (%f), tOn = %f (%f), tOff = %f (%f), 95%% = %f, max = %f", statsNBlinks.getMean(), statsSampledNBlinks.getMean(), statsTOn.getMean(), statsSampledTOn.getMean(), statsTOff.getMean(), statsSampledTOff.getMean(), statsSampledTOff.getStatistics().getPercentile(95), statsSampledTOff.getStatistics().getMax()));
logger.info("-=-=--=-");
final BlinkEstimator be = new BlinkEstimator();
be.setMaxDarkTime((int) (toff * 10));
be.setMsPerFrame(msPerFrame);
be.setRelativeDistance(false);
final double d = ImageModel.getRandomMoveDistance(diffusionRate);
be.setSearchDistance((fixedFraction < 1) ? Math.sqrt(2 * d * d) * 3 : 0);
be.setTimeAtLowerBound(timeAtLowerBound);
// Assertions.assertTrue("Max dark time must exceed the dark time of the data (otherwise no
// plateau)",
// be.maxDarkTime > statsSampledTOff.getStatistics().getMax());
final int nMolecules = fluorophores.size();
if (usePopulationStatistics) {
blinkingRate = statsNBlinks.getMean();
toff = statsTOff.getMean();
} else {
blinkingRate = statsSampledNBlinks.getMean();
toff = statsSampledTOff.getMean();
}
// See if any fitting regime gets a correct answer
final TIntHashSet ok = new TIntHashSet();
for (int numberOfFittedPoints = MIN_FITTED_POINTS; numberOfFittedPoints <= MAX_FITTED_POINTS; numberOfFittedPoints++) {
be.setNumberOfFittedPoints(numberOfFittedPoints);
be.computeBlinkingRate(results, true);
final double moleculesError = DoubleEquality.relativeError(nMolecules, be.getNMolecules());
final double blinksError = DoubleEquality.relativeError(blinkingRate, be.getNBlinks());
final double offError = DoubleEquality.relativeError(toff * msPerFrame, be.getTOff());
logger.info(FunctionUtils.getSupplier("Error %d: N = %f, blinks = %f, tOff = %f : %f", numberOfFittedPoints, moleculesError, blinksError, offError, (moleculesError + blinksError + offError) / 3));
if (moleculesError < relativeError && blinksError < relativeError && offError < relativeError) {
ok.add(numberOfFittedPoints);
logger.info("-=-=--=-");
logger.info(FunctionUtils.getSupplier("*** Correct at %d fitted points ***", numberOfFittedPoints));
if (doAssert) {
break;
}
}
// if (!be.isIncreaseNFittedPoints())
// break;
}
logger.info("-=-=--=-");
if (doAssert) {
Assertions.assertFalse(ok.isEmpty());
}
// relativeError);
return ok;
}
use of uk.ac.sussex.gdsc.smlm.model.MoleculeModel in project GDSC-SMLM by aherbert.
the class DiffusionRateTest method showExample.
private void showExample(int totalSteps, double diffusionSigma, UniformRandomProvider rng) {
final MoleculeModel m = new MoleculeModel(0, new double[3]);
final float[] xValues = new float[totalSteps];
final float[] x = new float[totalSteps];
final float[] y = new float[totalSteps];
final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
double[] axis;
if (diffusionType == DiffusionType.LINEAR_WALK) {
axis = nextVector(SamplerUtils.createNormalizedGaussianSampler(rng));
} else {
axis = null;
}
for (int j = 0; j < totalSteps; j++) {
if (diffusionType == DiffusionType.GRID_WALK) {
m.walk(diffusionSigma, rng);
} else if (diffusionType == DiffusionType.LINEAR_WALK) {
m.slide(diffusionSigma, axis, rng);
} else {
m.move(diffusionSigma, rng);
}
x[j] = (float) (m.getX());
y[j] = (float) (m.getY());
xValues[j] = (float) ((j + 1) / settings.getStepsPerSecond());
}
// Plot x and y coords on a timeline
final String title = TITLE + " example coordinates";
final Plot plot = new Plot(title, "Time (seconds)", "Distance (um)");
final float[] xUm = convertToUm(x);
final float[] yUm = convertToUm(y);
float[] limits = MathUtils.limits(xUm);
limits = MathUtils.limits(limits, yUm);
plot.setLimits(0, totalSteps / settings.getStepsPerSecond(), limits[0], limits[1]);
plot.setColor(Color.red);
plot.addPoints(xValues, xUm, Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(xValues, yUm, Plot.LINE);
ImageJUtils.display(title, plot);
// Scale up and draw 2D position
for (int j = 0; j < totalSteps; j++) {
x[j] *= pluginSettings.magnification;
y[j] *= pluginSettings.magnification;
}
final float[] limitsx = getLimits(x);
final float[] limitsy = getLimits(y);
int width = (int) (limitsx[1] - limitsx[0]);
int height = (int) (limitsy[1] - limitsy[0]);
// Ensure we draw something, even it is a simple dot at the centre for no diffusion
if (width == 0) {
width = (int) (32 * pluginSettings.magnification);
limitsx[0] = -width / 2.0f;
}
if (height == 0) {
height = (int) (32 * pluginSettings.magnification);
limitsy[0] = -height / 2.0f;
}
final ImageProcessor ip = new ByteProcessor(width, height);
// Adjust x and y using the minimum to centre
x[0] -= limitsx[0];
y[0] -= limitsy[0];
for (int j = 1; j < totalSteps; j++) {
// Adjust x and y using the minimum to centre
x[j] -= limitsx[0];
y[j] -= limitsy[0];
// Draw a line
ip.setColor(32 + (223 * j) / (totalSteps - 1));
ip.drawLine(round(x[j - 1]), round(y[j - 1]), round(x[j]), round(y[j]));
}
// Draw the final position
ip.putPixel(round(x[totalSteps - 1]), round(y[totalSteps - 1]), 255);
final ImagePlus imp = ImageJUtils.display(TITLE + " example", ip);
// Apply the fire lookup table
WindowManager.setTempCurrentImage(imp);
final LutLoader lut = new LutLoader();
lut.run("fire");
WindowManager.setTempCurrentImage(null);
}
use of uk.ac.sussex.gdsc.smlm.model.MoleculeModel in project GDSC-SMLM by aherbert.
the class DiffusionRateTest method run.
@Override
public void run(String arg) {
SmlmUsageTracker.recordPlugin(this.getClass(), arg);
pluginSettings = Settings.load();
pluginSettings.save();
if (IJ.controlKeyDown()) {
simpleTest();
return;
}
extraOptions = ImageJUtils.isExtraOptions();
if (!showDialog()) {
return;
}
lastSimulation.set(null);
final int totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
conversionFactor = 1000000.0 / (settings.getPixelPitch() * settings.getPixelPitch());
// Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
final double diffusionRateInPixelsPerSecond = settings.getDiffusionRate() * conversionFactor;
final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.getStepsPerSecond();
final double precisionInPixels = myPrecision / settings.getPixelPitch();
final boolean addError = myPrecision != 0;
ImageJUtils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", MathUtils.rounded(settings.getDiffusionRate(), 4), MathUtils.rounded(myPrecision, 4));
ImageJUtils.log("Mean-displacement per dimension = %s nm/sec", MathUtils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.getDiffusionRate()), 4));
if (extraOptions) {
ImageJUtils.log("Step size = %s, precision = %s", MathUtils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), MathUtils.rounded(precisionInPixels));
}
// Convert diffusion co-efficient into the standard deviation for the random walk
final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
final double diffusionSigma = ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
ImageJUtils.log("Simulation step-size = %s nm", MathUtils.rounded(settings.getPixelPitch() * diffusionSigma, 4));
// Move the molecules and get the diffusion rate
IJ.showStatus("Simulating ...");
final long start = System.nanoTime();
final UniformRandomProvider random = UniformRandomProviders.create();
final Statistics[] stats2D = new Statistics[totalSteps];
final Statistics[] stats3D = new Statistics[totalSteps];
final StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
final StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
for (int j = 0; j < totalSteps; j++) {
stats2D[j] = new Statistics();
stats3D[j] = new Statistics();
}
final SphericalDistribution dist = new SphericalDistribution(settings.getConfinementRadius() / settings.getPixelPitch());
final Statistics asymptote = new Statistics();
// Save results to memory
final MemoryPeakResults results = new MemoryPeakResults(totalSteps);
results.setCalibration(CalibrationHelper.create(settings.getPixelPitch(), 1, 1000.0 / settings.getStepsPerSecond()));
results.setName(TITLE);
results.setPsf(PsfHelper.create(PSFType.CUSTOM));
int peak = 0;
// Store raw coordinates
final ArrayList<Point> points = new ArrayList<>(totalSteps);
final StoredData totalJumpDistances1D = new StoredData(settings.getParticles());
final StoredData totalJumpDistances2D = new StoredData(settings.getParticles());
final StoredData totalJumpDistances3D = new StoredData(settings.getParticles());
final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
for (int i = 0; i < settings.getParticles(); i++) {
if (i % 16 == 0) {
IJ.showProgress(i, settings.getParticles());
if (ImageJUtils.isInterrupted()) {
return;
}
}
// Increment the frame so that tracing analysis can distinguish traces
peak++;
double[] origin = new double[3];
final int id = i + 1;
final MoleculeModel m = new MoleculeModel(id, origin.clone());
if (addError) {
origin = addError(origin, precisionInPixels, gauss);
}
if (pluginSettings.useConfinement) {
// Note: When using confinement the average displacement should asymptote
// at the average distance of a point from the centre of a ball. This is 3r/4.
// See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM
// The equivalent in 2D is 2r/3. However although we are plotting 2D distance
// this is a projection of the 3D position onto the plane and so the particles
// will not be evenly spread (there will be clustering at centre caused by the
// poles)
final double[] axis = (diffusionType == DiffusionType.LINEAR_WALK) ? nextVector(gauss) : null;
for (int j = 0; j < totalSteps; j++) {
double[] xyz = m.getCoordinates();
final double[] originalXyz = xyz.clone();
for (int n = pluginSettings.confinementAttempts; n-- > 0; ) {
if (diffusionType == DiffusionType.GRID_WALK) {
m.walk(diffusionSigma, random);
} else if (diffusionType == DiffusionType.LINEAR_WALK) {
m.slide(diffusionSigma, axis, random);
} else {
m.move(diffusionSigma, random);
}
if (!dist.isWithin(m.getCoordinates())) {
// Reset position
for (int k = 0; k < 3; k++) {
xyz[k] = originalXyz[k];
}
} else {
// The move was allowed
break;
}
}
points.add(new Point(id, xyz));
if (addError) {
xyz = addError(xyz, precisionInPixels, gauss);
}
peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
}
asymptote.add(distance(m.getCoordinates()));
} else if (diffusionType == DiffusionType.GRID_WALK) {
for (int j = 0; j < totalSteps; j++) {
m.walk(diffusionSigma, random);
double[] xyz = m.getCoordinates();
points.add(new Point(id, xyz));
if (addError) {
xyz = addError(xyz, precisionInPixels, gauss);
}
peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
}
} else if (diffusionType == DiffusionType.LINEAR_WALK) {
final double[] axis = nextVector(gauss);
for (int j = 0; j < totalSteps; j++) {
m.slide(diffusionSigma, axis, random);
double[] xyz = m.getCoordinates();
points.add(new Point(id, xyz));
if (addError) {
xyz = addError(xyz, precisionInPixels, gauss);
}
peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
}
} else {
for (int j = 0; j < totalSteps; j++) {
m.move(diffusionSigma, random);
double[] xyz = m.getCoordinates();
points.add(new Point(id, xyz));
if (addError) {
xyz = addError(xyz, precisionInPixels, gauss);
}
peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
}
}
// Debug: record all the particles so they can be analysed
// System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ());
final double[] xyz = m.getCoordinates();
double d2 = 0;
totalJumpDistances1D.add(d2 = xyz[0] * xyz[0]);
totalJumpDistances2D.add(d2 += xyz[1] * xyz[1]);
totalJumpDistances3D.add(d2 += xyz[2] * xyz[2]);
}
final long nanoseconds = System.nanoTime() - start;
IJ.showProgress(1);
MemoryPeakResults.addResults(results);
simulation = new SimulationData(results.getName(), myPrecision);
// Convert pixels^2/step to um^2/sec
final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
ImageJUtils.log("Raw data D=%s um^2/s, Precision = %s nm, N=%d, step=%s s, mean2D=%s um^2, " + "MSD 2D = %s um^2/s, mean3D=%s um^2, MSD 3D = %s um^2/s", MathUtils.rounded(settings.getDiffusionRate()), MathUtils.rounded(myPrecision), jumpDistances2D.getN(), MathUtils.rounded(results.getCalibrationReader().getExposureTime() / 1000), MathUtils.rounded(jumpDistances2D.getMean() / conversionFactor), MathUtils.rounded(msd2D), MathUtils.rounded(jumpDistances3D.getMean() / conversionFactor), MathUtils.rounded(msd3D));
aggregateIntoFrames(points, addError, precisionInPixels, gauss);
IJ.showStatus("Analysing results ...");
if (pluginSettings.showDiffusionExample) {
showExample(totalSteps, diffusionSigma, random);
}
// Plot a graph of mean squared distance
final double[] xValues = new double[stats2D.length];
final double[] yValues2D = new double[stats2D.length];
final double[] yValues3D = new double[stats3D.length];
final double[] upper2D = new double[stats2D.length];
final double[] lower2D = new double[stats2D.length];
final double[] upper3D = new double[stats3D.length];
final double[] lower3D = new double[stats3D.length];
final SimpleRegression r2D = new SimpleRegression(false);
final SimpleRegression r3D = new SimpleRegression(false);
final int firstN = (pluginSettings.useConfinement) ? pluginSettings.fitN : totalSteps;
for (int j = 0; j < totalSteps; j++) {
// Convert steps to seconds
xValues[j] = (j + 1) / settings.getStepsPerSecond();
// Convert values in pixels^2 to um^2
final double mean2D = stats2D[j].getMean() / conversionFactor;
final double mean3D = stats3D[j].getMean() / conversionFactor;
final double sd2D = stats2D[j].getStandardDeviation() / conversionFactor;
final double sd3D = stats3D[j].getStandardDeviation() / conversionFactor;
yValues2D[j] = mean2D;
yValues3D[j] = mean3D;
upper2D[j] = mean2D + sd2D;
lower2D[j] = mean2D - sd2D;
upper3D[j] = mean3D + sd3D;
lower3D[j] = mean3D - sd3D;
if (j < firstN) {
r2D.addData(xValues[j], yValues2D[j]);
r3D.addData(xValues[j], yValues3D[j]);
}
}
// TODO - Fit using the equation for 2D confined diffusion:
// MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2)
// s = localisation precision
// R = confinement radius
// D = 2D diffusion coefficient
// t = time
final PolynomialFunction fitted2D;
final PolynomialFunction fitted3D;
if (r2D.getN() > 0) {
// Do linear regression to get diffusion rate
final double[] best2D = new double[] { r2D.getIntercept(), r2D.getSlope() };
fitted2D = new PolynomialFunction(best2D);
final double[] best3D = new double[] { r3D.getIntercept(), r3D.getSlope() };
fitted3D = new PolynomialFunction(best3D);
// For 2D diffusion: d^2 = 4D
// where: d^2 = mean-square displacement
double diffCoeff = best2D[1] / 4.0;
final String msg = "2D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")";
IJ.showStatus(msg);
ImageJUtils.log(msg);
diffCoeff = best3D[1] / 6.0;
ImageJUtils.log("3D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")");
} else {
fitted2D = fitted3D = null;
}
// Create plots
plotMsd(totalSteps, xValues, yValues2D, lower2D, upper2D, fitted2D, 2);
plotMsd(totalSteps, xValues, yValues3D, lower3D, upper3D, fitted3D, 3);
plotJumpDistances(TITLE, jumpDistances2D, 2, 1);
plotJumpDistances(TITLE, jumpDistances3D, 3, 1);
// Show the total jump length for debugging
// plotJumpDistances(TITLE + " total", totalJumpDistances1D, 1, totalSteps);
// plotJumpDistances(TITLE + " total", totalJumpDistances2D, 2, totalSteps);
// plotJumpDistances(TITLE + " total", totalJumpDistances3D, 3, totalSteps);
windowOrganiser.tile();
if (pluginSettings.useConfinement) {
ImageJUtils.log("3D asymptote distance = %s nm (expected %.2f)", MathUtils.rounded(asymptote.getMean() * settings.getPixelPitch(), 4), 3 * settings.getConfinementRadius() / 4);
}
}
Aggregations