use of uk.ac.sussex.gdsc.smlm.model.LocalisationModelSet in project GDSC-SMLM by aherbert.
the class CreateData method run.
@Override
public void run(String arg) {
SmlmUsageTracker.recordPlugin(this.getClass(), arg);
extraOptions = ImageJUtils.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 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.setExposureTime(1000);
areaInUm = settings.getSize() * settings.getPixelPitch() * settings.getSize() * settings.getPixelPitch() / 1e6;
// Number of spots per frame
int count = 0;
int[] nextN = null;
SpatialDistribution dist;
if (benchmarkMode) {
// --------------------
// BENCHMARK SIMULATION
// --------------------
// Draw the same point on the image repeatedly
count = 1;
dist = createFixedDistribution();
try {
reportAndSaveFittingLimits(dist);
} catch (final Exception ex) {
// This will be from the computation of the CRLB
IJ.error(TITLE, ex.getMessage());
return;
}
} 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.getParticles() * 2];
Arrays.fill(nextN, 0, settings.getParticles(), 1);
RandomUtils.shuffle(nextN, UniformRandomProviders.create());
// Only put spots in the central part of the image
final double border = settings.getSize() / 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.getSamplePerFrame()) {
final double mean = areaInUm * settings.getDensity();
ImageJUtils.log("Mean samples = %f", mean);
if (mean < 0.5) {
final GenericDialog gd = new GenericDialog(TITLE);
gd.addMessage("The mean samples per frame is low: " + MathUtils.rounded(mean) + "\n \nContinue?");
gd.enableYesNoCancel();
gd.hideCancelButton();
gd.showDialog();
if (!gd.wasOKed()) {
return;
}
}
final PoissonSampler poisson = new PoissonSampler(createRandomGenerator(), mean);
final StoredDataStatistics samples = new StoredDataStatistics(settings.getParticles());
while (samples.getSum() < settings.getParticles()) {
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
count = (int) Math.max(1, Math.round(areaInUm * settings.getDensity()));
}
}
UniformRandomProvider rng = null;
localisationSets = new ArrayList<>(settings.getParticles());
final int minPhotons = (int) settings.getPhotonsPerSecond();
final int range = (int) settings.getPhotonsPerSecondMaximum() - minPhotons + 1;
if (range > 1) {
rng = createRandomGenerator();
}
// Add frames at the specified density until the number of particles has been reached
int id = 0;
int time = 0;
while (id < settings.getParticles()) {
// Allow the number per frame to be specified
if (nextN != null) {
if (time >= nextN.length) {
break;
}
count = nextN[time];
}
// Simulate random positions in the frame for the specified density
time++;
for (int j = 0; j < count; 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 + ((rng != null) ? rng.nextInt(range) : 0);
final LocalisationModel m = new LocalisationModel(id, time, xyz, intensity, LocalisationModel.CONTINUOUS);
// Each localisation can be a separate localisation set
final LocalisationModelSet set = new LocalisationModelSet(id, time);
set.add(m);
localisationSets.add(set);
id++;
}
}
} else {
if (!showDialog()) {
return;
}
resetMemory();
areaInUm = settings.getSize() * settings.getPixelPitch() * settings.getSize() * settings.getPixelPitch() / 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.setSeconds((int) Math.ceil(settings.getParticles() * (settings.getExposureTime() + settings.getTOn()) / 1000));
totalSteps = 0;
final double simulationStepsPerFrame = (settings.getStepsPerSecond() * settings.getExposureTime()) / 1000.0;
imageModel = new FixedLifetimeImageModel(settings.getStepsPerSecond() * settings.getTOn() / 1000.0, simulationStepsPerFrame, createRandomGenerator());
} 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.
final SpatialIllumination activationIllumination = createIllumination(settings.getPulseRatio(), settings.getPulseInterval());
// Generate additional frames so that each frame has the set number of simulation steps
totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
// Since we have an exponential decay of activations
// ensure half of the particles have activated by 30% of the frames.
final double eAct = totalSteps * 0.3 * activationIllumination.getAveragePhotons();
// Q. Does tOn/tOff change depending on the illumination strength?
imageModel = new ActivationEnergyImageModel(eAct, activationIllumination, settings.getStepsPerSecond() * settings.getTOn() / 1000.0, settings.getStepsPerSecond() * settings.getTOffShort() / 1000.0, settings.getStepsPerSecond() * settings.getTOffLong() / 1000.0, settings.getNBlinksShort(), settings.getNBlinksLong(), createRandomGenerator());
imageModel.setUseGeometricDistribution(settings.getNBlinksGeometricDistribution());
// Only use the correlation if selected for the distribution
if (PHOTON_DISTRIBUTION[PHOTON_CORRELATED].equals(settings.getPhotonDistribution())) {
correlation = settings.getCorrelation();
}
}
imageModel.setPhotonBudgetPerFrame(true);
imageModel.setDiffusion2D(settings.getDiffuse2D());
imageModel.setRotation2D(settings.getRotate2D());
IJ.showStatus("Creating molecules ...");
final SpatialDistribution distribution = createDistribution();
final List<CompoundMoleculeModel> compounds = createCompoundMolecules();
if (compounds == null) {
return;
}
final List<CompoundMoleculeModel> molecules = imageModel.createMolecules(compounds, settings.getParticles(), distribution, settings.getRotateInitialOrientation());
// 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;
}
// Map the fluorophore ID to the compound for mixtures
if (compounds.size() > 1) {
idToCompound = new TIntIntHashMap(fluorophores.size());
for (final FluorophoreSequenceModel l : fluorophores) {
idToCompound.put(l.getId(), l.getLabel());
}
}
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());
try {
imageModel.setConfinementDistribution(createConfinementDistribution());
} catch (final ConfigurationException ex) {
// We asked the user if it was OK to continue and they said no
return;
}
// This should be optimised
imageModel.setConfinementAttempts(10);
final List<LocalisationModel> localisations = imageModel.createImage(molecules, settings.getFixedFraction(), totalSteps, settings.getPhotonsPerSecond() / settings.getStepsPerSecond(), correlation, settings.getRotateDuringSimulation());
// Re-adjust the fluorophores to the correct time
if (settings.getStepsPerSecond() != 1) {
final double scale = 1.0 / settings.getStepsPerSecond();
for (final FluorophoreSequenceModel f : fluorophores) {
f.adjustTime(scale);
}
}
// Integrate the frames
localisationSets = combineSimulationSteps(localisations);
localisationSets = filterToImageBounds(localisationSets);
}
datasetNumber.getAndIncrement();
final List<LocalisationModel> localisations = drawImage(localisationSets);
if (localisations == null || localisations.isEmpty()) {
IJ.error(TITLE, "No localisations created");
return;
}
fluorophores = removeFilteredFluorophores(fluorophores, localisations);
final double signalPerFrame = showSummary(fluorophores, localisations);
if (!benchmarkMode) {
final boolean fullSimulation = (!(simpleMode || spotMode));
saveSimulationParameters(localisations.size(), fullSimulation, signalPerFrame);
}
IJ.showStatus("Saving data ...");
saveFluorophores(fluorophores);
saveImageResults(results);
saveLocalisations(localisations);
// The settings for the filenames may have changed
SettingsManager.writeSettings(settings.build());
IJ.showStatus("Done");
}
use of uk.ac.sussex.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 the localisation sets
* @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();
removedT1 = new AtomicInteger();
removedTn = new AtomicInteger();
photonStats = new SummaryStatistics();
// Add drawn spots to memory
results = new MemoryPeakResults();
final CalibrationWriter c = new CalibrationWriter();
c.setDistanceUnit(DistanceUnit.PIXEL);
c.setIntensityUnit(IntensityUnit.PHOTON);
c.setNmPerPixel(settings.getPixelPitch());
c.setExposureTime(settings.getExposureTime());
c.setCameraType(settings.getCameraType());
if (settings.getCameraType() == CameraType.SCMOS) {
c.setCameraModelName(settings.getCameraModelName());
} else {
final CreateDataSettingsHelper helper = new CreateDataSettingsHelper(settings);
c.setCountPerPhoton(helper.getTotalGainSafe());
c.setBias(settings.getBias());
c.setReadNoise(helper.getReadNoiseInCounts());
c.setQuantumEfficiency(helper.getQuantumEfficiency());
}
results.setCalibration(c.getCalibration());
results.setSortAfterEnd(true);
results.begin();
maxT = localisationSets.get(localisationSets.size() - 1).getTime();
// Display image
final ImageStack stack = new ImageStack(settings.getSize(), settings.getSize(), maxT);
final double psfSd = getPsfSd();
if (psfSd <= 0) {
return null;
}
final PsfModel psfModel = createPsfModel(localisationSets);
if (psfModel == null) {
return null;
}
// Create the camera noise model
createPerPixelCameraModelData(cameraModel);
createBackgroundPixels();
final int threadCount = Prefs.getThreads();
IJ.showStatus("Drawing image ...");
// Multi-thread for speed
final PeakResults syncResults = SynchronizedPeakResults.create(results, threadCount);
final ExecutorService threadPool = Executors.newFixedThreadPool(threadCount);
final List<Future<?>> futures = new LinkedList<>();
// Count all the frames to process
final Ticker ticker = ImageJUtils.createTicker(maxT, threadCount);
// Collect statistics on the number of photons actually simulated
// Process all frames
int index = 0;
int lastT = -1;
for (final LocalisationModelSet l : localisationSets) {
if (ImageJUtils.isInterrupted()) {
break;
}
if (l.getTime() != lastT) {
lastT = l.getTime();
futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, index, lastT, psfModel.copy(), syncResults, stack, poissonNoise, createRandomGenerator(), ticker)));
}
index++;
}
// Finish processing data
ConcurrencyUtils.waitForCompletionUnchecked(futures);
futures.clear();
if (ImageJUtils.isInterrupted()) {
IJ.showProgress(1);
return null;
}
// Do all the frames that had no localisations
float[] limits = null;
for (int t = 1; t <= maxT; t++) {
if (ImageJUtils.isInterrupted()) {
break;
}
final Object pixels = stack.getPixels(t);
if (pixels == null) {
futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, maxT, t, null, syncResults, stack, poissonNoise, createRandomGenerator(), ticker)));
} else if (limits == null) {
limits = MathUtils.limits((float[]) pixels);
}
}
// Finish
ConcurrencyUtils.waitForCompletionUnchecked(futures);
futures.clear();
threadPool.shutdown();
IJ.showProgress(1);
if (ImageJUtils.isInterrupted() || limits == null) {
return null;
}
results.end();
if (photonsRemoved.get() > 0) {
ImageJUtils.log("Removed %d localisations with less than %.1f rendered photons", photonsRemoved.get(), settings.getMinPhotons());
}
if (removedT1.get() > 0) {
ImageJUtils.log("Removed %d localisations with no neighbours @ SNR %.2f", removedT1.get(), settings.getMinSnrT1());
}
if (removedTn.get() > 0) {
ImageJUtils.log("Removed %d localisations with valid neighbours @ SNR %.2f", removedTn.get(), settings.getMinSnrTN());
}
if (photonStats.getN() > 0) {
ImageJUtils.log("Average photons rendered = %s +/- %s", MathUtils.rounded(photonStats.getMean()), MathUtils.rounded(photonStats.getStandardDeviation()));
}
// System.out.printf("rawPhotons = %f\n", rawPhotons.getMean());
// System.out.printf("drawPhotons = %f\n", drawPhotons.getMean());
// new HistogramPlotBuilder("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.getRawImage()) {
// Get the global limits and ensure all values can be represented
final Object[] imageArray = stack.getImageArray();
limits = MathUtils.limits((float[]) imageArray[0]);
for (int j = 1; j < imageArray.length; j++) {
limits = MathUtils.limits(limits, (float[]) imageArray[j]);
}
// float limits0 = limits[0];
// Leave bias in place
final float limits0 = 0;
// Check if the image will fit in a 16-bit range
if ((limits[1] - limits0) < 65535) {
// Convert to 16-bit
newStack = new ImageStack(stack.getWidth(), stack.getHeight(), stack.getSize());
// Account for rounding
final float min = (float) (limits0 - 0.5);
for (int j = 0; j < imageArray.length; j++) {
final float[] image = (float[]) imageArray[j];
final 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 (MemoryUtils.getFreeMemory() < 33554432L) {
MemoryUtils.runGarbageCollectorOnce();
}
}
for (int k = 2; k-- > 0; ) {
limits[k] = (float) Math.floor(limits[k] - min);
}
} else {
// Keep as 32-bit but round to whole numbers
for (int j = 0; j < imageArray.length; j++) {
final float[] pixels = (float[]) imageArray[j];
for (int k = 0; k < pixels.length; k++) {
pixels[k] = Math.round(pixels[k]);
}
}
for (int k = 2; k-- > 0; ) {
limits[k] = Math.round(limits[k]);
}
}
}
// Show image
final ImagePlus imp = ImageJUtils.display(CREATE_DATA_IMAGE_TITLE, newStack);
final ij.measure.Calibration cal = new ij.measure.Calibration();
String unit = "nm";
double unitPerPixel = settings.getPixelPitch();
if (unitPerPixel > 100) {
unit = "um";
unitPerPixel /= 1000.0;
}
cal.setUnit(unit);
cal.pixelHeight = cal.pixelWidth = unitPerPixel;
final Rectangle bounds = cameraModel.getBounds();
if (bounds != null) {
cal.xOrigin = -bounds.x;
cal.yOrigin = -bounds.y;
}
imp.setCalibration(cal);
imp.setDimensions(1, 1, newStack.getSize());
imp.setDisplayRange(limits[0], limits[1]);
imp.updateAndDraw();
saveImage(imp);
final IJImageSource imageSource = new IJImageSource(imp);
// Shift simulation image source to correct location
results.setSource(imageSource);
results.setName(CREATE_DATA_IMAGE_TITLE + " (" + TITLE + ")");
// Bounds are relative to the image source
results.setBounds(new Rectangle(settings.getSize(), settings.getSize()));
results.setPsf(createPsf(psfSd));
MemoryPeakResults.addResults(results);
setBenchmarkResults(imp, results);
if (benchmarkMode && benchmarkParameters != null) {
benchmarkParameters.setPhotons(results);
}
final List<LocalisationModel> localisations = toLocalisations(localisationSets);
savePulses(localisations, results);
// Saved the fixed and moving localisations into different datasets
saveFixedAndMoving(results);
saveCompoundMolecules(results);
return localisations;
}
use of uk.ac.sussex.gdsc.smlm.model.LocalisationModelSet in project GDSC-SMLM by aherbert.
the class CreateData method filterToImageBounds.
/**
* Filter those not in the bounds of the distribution.
*
* @param localisationSets the localisation sets
* @return the list
*/
private List<LocalisationModelSet> filterToImageBounds(List<LocalisationModelSet> localisationSets) {
final List<LocalisationModelSet> newLocalisations = new ArrayList<>(localisationSets.size());
final SpatialDistribution bounds = createUniformDistribution(0);
for (final LocalisationModelSet s : localisationSets) {
if (bounds.isWithinXy(s.toLocalisation().getCoordinates())) {
newLocalisations.add(s);
}
}
return newLocalisations;
}
use of uk.ac.sussex.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.getStepsPerSecond() * settings.getExposureTime()) / 1000.0;
final List<LocalisationModelSet> newLocalisations = new ArrayList<>((int) (localisations.size() / simulationStepsPerFrame));
// System.out.printf("combineSimulationSteps @ %f\n", simulationStepsPerFrame);
// final double gain = new CreateDataSettingsHelper(settings).getTotalGainSafe();
sortLocalisationsByIdThenTime(localisations);
final int[] idList = getIds(localisations);
movingMolecules = new TIntHashSet(idList.length);
int index = 0;
for (final int id : idList) {
final int fromIndex = findIndexById(localisations, index, id);
if (fromIndex > -1) {
final int toIndex = findLastIndexById(localisations, fromIndex, id);
final 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);
final LocalisationModelSet[] sets = new LocalisationModelSet[intLastFrame - intFirstFrame + 1];
// Process each step
for (final LocalisationModel l : subset) {
// Get the fractional start and end frames
final double startFrame = getStartFrame(l, simulationStepsPerFrame);
final double endFrame = getEndFrame(l, simulationStepsPerFrame);
// Round down to get the actual frames that are overlapped
final int start = (int) startFrame;
int end = (int) endFrame;
// 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
final 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
final 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.
// This is used later to filter based on SNR
sets[i].setData(new double[] { // * gain
0, // * gain
0, // * gain
0, // * gain
0, // * gain
sets[i].getIntensity() });
newLocalisations.add(sets[i]);
}
previous = sets[i];
}
}
}
// Sort by time
Collections.sort(newLocalisations, (r1, r2) -> Integer.compare(r1.getTime(), r2.getTime()));
return newLocalisations;
}
use of uk.ac.sussex.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 the localisation sets
* @return the image PSF model
*/
private ImagePsfModel createImagePsf(List<LocalisationModelSet> localisationSets) {
final ImagePlus imp = WindowManager.getImage(settings.getPsfImageName());
if (imp == null) {
IJ.error(TITLE, "Unable to create the PSF model from image: " + settings.getPsfImageName());
return null;
}
try {
final ImagePSF psfSettings = ImagePsfHelper.fromString(imp.getProperty("Info").toString());
if (psfSettings == null) {
throw new IllegalStateException("Unknown PSF settings for image: " + imp.getTitle());
}
// Check all the settings have values
if (psfSettings.getPixelSize() <= 0) {
throw new IllegalStateException("Missing nmPerPixel calibration settings for image: " + imp.getTitle());
}
if (psfSettings.getPixelDepth() <= 0) {
throw new IllegalStateException("Missing nmPerSlice calibration settings for image: " + imp.getTitle());
}
if (psfSettings.getCentreImage() <= 0) {
throw new IllegalStateException("Missing zCentre calibration settings for image: " + imp.getTitle());
}
if (psfSettings.getFwhm() <= 0) {
throw new IllegalStateException("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;
double maxZ = Double.NEGATIVE_INFINITY;
for (final LocalisationModelSet l : localisationSets) {
for (final LocalisationModel m : l.getLocalisations()) {
final double z = m.getZ();
if (minZ > z) {
minZ = z;
}
if (maxZ < z) {
maxZ = z;
}
}
}
final int nSlices = imp.getStackSize();
// z-centre should be an index and not the ImageJ slice number so subtract 1
final int zCentre = psfSettings.getCentreImage() - 1;
// Calculate the start/end slices to cover the depth of field
// This logic must match the ImagePSFModel.
final double unitsPerSlice = psfSettings.getPixelDepth() / settings.getPixelPitch();
// 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;
// Add extra to the range so that gradients can be computed.
lower--;
upper++;
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 float[][] image = extractImageStack(imp, lower, upper);
final ImagePsfModel model = new ImagePsfModel(image, zCentre - lower, psfSettings.getPixelSize() / settings.getPixelPitch(), unitsPerSlice, noiseFraction);
// Add the calibrated centres. The map will not be null
final Map<Integer, Offset> map = psfSettings.getOffsetsMap();
if (!map.isEmpty()) {
final int sliceOffset = lower + 1;
for (final Entry<Integer, Offset> entry : map.entrySet()) {
model.setRelativeCentre(entry.getKey() - sliceOffset, entry.getValue().getCx(), entry.getValue().getCy());
}
} else {
// Use the CoM if present
final double cx = psfSettings.getXCentre();
final double cy = psfSettings.getYCentre();
if (cx != 0 || cy != 0) {
for (int slice = 0; slice < image.length; slice++) {
model.setCentre(slice, cx, cy);
}
}
}
// Initialise the HWHM table so that it can be cloned
model.initialiseHwhm();
return model;
} catch (final Exception ex) {
IJ.error(TITLE, "Unable to create the image PSF model:\n" + ex.getMessage());
return null;
}
}
Aggregations