use of uk.ac.sussex.gdsc.smlm.model.camera.CcdCameraModel in project GDSC-SMLM by aherbert.
the class PsfCreator method runUsingAlignment.
private void runUsingAlignment() {
if (!showAlignmentDialog()) {
return;
}
boxRadius = (int) Math.ceil(settings.getRadius());
final CalibrationReader calibration = new CalibrationReader(settings.getCalibration());
// Limit this
final int halfBoxRadius = boxRadius / 2;
settings.setAnalysisWindow(Math.min(settings.getAnalysisWindow(), halfBoxRadius));
// Find the selected PSF spots x,y,z centre
// We offset the centre to the middle of pixel.
BasePoint[] centres = getSpots(0.5f, false);
if (centres.length == 0) {
IJ.error(TITLE, "No PSFs");
return;
}
CameraModel cameraModel = null;
if (calibration.isScmos()) {
cameraModel = CameraModelManager.load(calibration.getCameraModelName());
if (cameraModel == null) {
IJ.error(TITLE, "No camera model");
return;
}
cameraModel = PeakFit.cropCameraModel(cameraModel, IJImageSource.getBounds(imp), null, true);
} else {
cameraModel = new CcdCameraModel(calibration.getBias(), 1);
}
// Extract the image data for processing as float
final float[][] image = CreateData.extractImageStack(imp, 0, imp.getStackSize() - 1);
for (final float[] data : image) {
cameraModel.removeBiasAndGain(data);
}
zSelector = new PsfCentreSelector();
// Relocate the initial centres
ImageJUtils.showStatus("Relocating initial centres");
centres = relocateCentres(image, centres);
if (centres == null) {
return;
}
zRadius = (int) Math.ceil(settings.getAlignmentZRadius());
// Check the region overlap in 3D and exclude overlapping PSFs
boolean[] bad = findSpotOverlap(centres, null);
centres = getNonBadSpots(centres, bad);
if (centres.length == 0) {
IJ.error(TITLE, "No PSFs without neighbours within the box region");
return;
}
// Multi-thread for speed
if (threadPool == null) {
threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
}
// Extract each PSF into a scaled PSF
ImageJUtils.showStatus(String.format("[%d] Extracting PSFs", 0));
ExtractedPsf[] psfs = extractPsfs(image, centres);
Point location = null;
// Iterate until centres have converged
boolean converged = false;
for (int iter = 0; !converged && iter < settings.getMaxIterations(); iter++) {
if (ImageJUtils.isInterrupted()) {
return;
}
// Combine all PSFs
ImageJUtils.showStatus(String.format("[%d] Aligning PSFs", iter + 1));
final ExtractedPsf combined = combine(psfs);
combined.createProjections();
// Get the current combined z-centre.
// This is used to get the centre of mass for repositioning.
// It also effects the alignment so do it for the first iteration.
zSelector.setPsf(combined);
if (iter == 0) {
// TODO - check if the z-centre should be guessed here.
// We assume that the combined PSF may be easier to guess if the initial
// guess for each individual PSF was OK. It may not be necessary since all
// the PSFs are combined around their z-centres. Once alignment has
// started we skip this step.
zSelector.analyse();
zSelector.guessZCentre();
}
if (settings.getInteractiveMode()) {
if (iter != 0) {
zSelector.analyse();
}
// zSelector.guessZCentre();
final double dz = zSelector.run("Update combined PSF z-centre", true, false, false, null);
if (dz < 0) {
return;
}
}
// Align each to the combined PSF
final float[][] translation = align(combined, psfs);
if (ImageJUtils.isInterrupted()) {
return;
}
// Find the new centre using the old centre plus the alignment shift
for (int j = 0; j < psfs.length; j++) {
centres[j] = psfs[j].updateCentre(translation[j]);
// Update to get the correct scale
translation[j][0] = centres[j].getX() - psfs[j].centre.getX();
translation[j][1] = centres[j].getY() - psfs[j].centre.getY();
translation[j][2] = centres[j].getZ() - psfs[j].centre.getZ();
ImageJUtils.log("[%d] Centre %d : Shift X = %s : Shift Y = %s : Shift Z = %s", iter, j + 1, rounder.toString(translation[j][0]), rounder.toString(translation[j][1]), rounder.toString(translation[j][2]));
}
final boolean[] excluded = new boolean[psfs.length];
if (checkAlignments) {
combined.show(TITLE_PSF);
// Ask about each centre in turn.
// Update Point ROI using float coordinates and set image slice to
// correct z-centre.
// imp.saveRoi();
imp.killRoi();
final ImageCanvas ic = imp.getCanvas();
// ic.setMagnification(16);
int reject = 0;
final float box = boxRadius + 0.5f;
final int n = imp.getStackSize();
for (int j = 0; j < centres.length; j++) {
psfs[j].show(TITLE_SPOT_PSF);
final Overlay o = new Overlay();
o.add(createRoi(psfs[j].centre.getX(), psfs[j].centre.getY(), Color.RED));
final float cx = centres[j].getX();
final float cy = centres[j].getY();
o.add(createRoi(cx, cy, Color.GREEN));
final Roi roi = new Roi(cx - box, cy - box, 2 * box, 2 * box);
o.add(roi);
// The centre is absolute within the original stack
imp.setSlice(MathUtils.clip(1, n, Math.round(centres[j].getZ())));
final Rectangle r = ic.getSrcRect();
final int x = centres[j].getXint();
final int y = centres[j].getYint();
if (!r.contains(x, y)) {
r.x = x - r.width / 2;
r.y = y - r.height / 2;
ic.setSourceRect(r);
}
imp.setOverlay(o);
imp.updateAndDraw();
final NonBlockingExtendedGenericDialog gd = new NonBlockingExtendedGenericDialog(TITLE);
ImageJUtils.addMessage(gd, "Shift X = %s\nShift Y = %s\nShift Z = %s", rounder.toString(translation[j][0]), rounder.toString(translation[j][1]), rounder.toString(translation[j][2]));
final int spotIndex = j;
gd.addAndGetButton("Exclude spot", event -> {
if (excluded[spotIndex]) {
ImageJUtils.log("Included spot %d", spotIndex + 1);
excluded[spotIndex] = false;
} else {
ImageJUtils.log("Excluded spot %d", spotIndex + 1);
excluded[spotIndex] = true;
}
});
gd.enableYesNoCancel("Accept", "Reject");
if (location != null) {
gd.setLocation(location.x, location.y);
}
gd.showDialog();
if (gd.wasCanceled()) {
resetImp();
return;
}
final boolean failed = excluded[spotIndex] || !gd.wasOKed();
if (failed) {
reject++;
centres[j] = psfs[j].centre;
// For RMSD computation
Arrays.fill(translation[j], 0f);
}
location = gd.getLocation();
}
resetImp();
if (reject == psfs.length) {
IJ.error(TITLE, "No PSF translations were accepted");
return;
}
}
bad = findSpotOverlap(centres, excluded);
final int badCount = count(bad);
final int excludedCount = count(excluded);
int ok = bad.length - badCount - excludedCount;
if (ok < bad.length) {
if (badCount != 0 && settings.getInteractiveMode()) {
final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Warning: Regions now overlap!");
gd.addMessage("OK = " + TextUtils.pleural(ok, "PSF"));
gd.addMessage("Overlapping = " + TextUtils.pleural(badCount, "PSF"));
// gd.addMessage("Excluded = " + TextUtils.pleural(excludedCount, "PSF"));
gd.enableYesNoCancel("Exclude", "Include");
if (location != null) {
gd.setLocation(location.x, location.y);
}
gd.showDialog();
if (gd.wasCanceled()) {
resetImp();
return;
}
if (!gd.wasOKed()) {
// allow bad spots
Arrays.fill(bad, false);
ok = bad.length;
}
location = gd.getLocation();
}
if (ok == 0) {
IJ.error(TITLE, "No PSFs remaining");
resetImp();
return;
}
}
// Merge bad and excluded to get new centres
for (int i = 0; i < bad.length; i++) {
if (excluded[i]) {
bad[i] = true;
}
}
ok = bad.length - count(bad);
final BasePoint[] newCentres = getNonBadSpots(centres, bad);
// Find the change in centres
final double[] rmsd = new double[2];
for (int j = 0; j < psfs.length; j++) {
if (bad[j]) {
continue;
}
rmsd[0] += MathUtils.pow2(translation[j][0]) + MathUtils.pow2(translation[j][1]);
rmsd[1] += MathUtils.pow2(translation[j][2]);
}
for (int j = 0; j < 2; j++) {
rmsd[j] = Math.sqrt(rmsd[j] / ok);
}
ImageJUtils.showStatus(String.format("[%d] Checking combined PSF", iter + 1));
// Compute CoM shift using the current z-centre and z-window
final double[] shift = combined.getCentreOfMassXyShift(zSelector.getCentreSlice());
final double shiftd = Math.sqrt(shift[0] * shift[0] + shift[1] * shift[1]);
ImageJUtils.log("[%d] RMSD XY = %s : RMSD Z = %s : Combined CoM shift = %s,%s (%s)", iter, rounder.toString(rmsd[0]), rounder.toString(rmsd[1]), rounder.toString(shift[0]), rounder.toString(shift[1]), rounder.toString(shiftd));
if (settings.getInteractiveMode()) {
// Ask if OK to continue?
final GenericDialog gd = new GenericDialog(TITLE);
ImageJUtils.addMessage(gd, "RMSD XY = %s\nRMSD Z = %s\nCombined CoM shift = %s,%s (%s)", rounder.toString(rmsd[0]), rounder.toString(rmsd[1]), rounder.toString(shift[0]), rounder.toString(shift[1]), rounder.toString(shiftd));
// Check if we can do more iterations
if (iter + 1 < settings.getMaxIterations()) {
gd.enableYesNoCancel("Continue", "Converged");
} else {
gd.setOKLabel("Converged");
}
gd.showDialog();
if (gd.wasCanceled()) {
return;
}
converged = !gd.wasOKed();
} else {
// Check convergence thresholds
converged = rmsd[0] < settings.getRmsdXyThreshold() && rmsd[1] < settings.getRmsdZThreshold() && shiftd < settings.getComShiftThreshold();
}
// For the next round we move to the non-overlapping spots
centres = newCentres;
// Update the centres using the centre-of-mass of the combined PSF
centres = updateUsingCentreOfMassXyShift(shift, shiftd, combined, centres);
// Extract each PSF into a scaled PSF
ImageJUtils.showStatus(String.format("[%d] Extracting PSFs", iter + 1));
psfs = extractPsfs(image, centres);
}
// Combine all
ExtractedPsf combined = combine(psfs);
// Show an interactive dialog for cropping the PSF and choosing the
// final output
final PsfOutputSelector cropSelector = new PsfOutputSelector(combined);
combined = cropSelector.run();
if (combined == null) {
return;
}
if (settings.getUpdateRoi()) {
final float[] ox = new float[centres.length];
final float[] oy = new float[centres.length];
for (int i = 0; i < centres.length; i++) {
ox[i] = centres[i].getX();
oy[i] = centres[i].getY();
}
imp.setRoi(new OffsetPointRoi(ox, oy));
}
// For an image PSF we can just enlarge the PSF and window.
// For a CSpline then we already have the 3D cubic spline function.
// However we want to post-process the function to allow windowing and
// normalisation. So we enlarge by 3 in each dimension.
// The CSpline can be created by solving the coefficients for the
// 4x4x4 (64) sampled points on each node.
int magnification;
if (settings.getOutputType() == OUTPUT_TYPE_IMAGE_PSF) {
magnification = settings.getPsfMagnification();
} else {
magnification = 3;
}
// Enlarge the combined PSF for final processing
ExtractedPsf finalPsf = combined.enlarge(magnification, threadPool);
// Show a dialog to collect final z-centre interactively
ImageJUtils.showStatus("Analysing PSF");
zSelector.setPsf(finalPsf);
zSelector.analyse();
// zSelector.guessZCentre(); // No need to guess the centre
final double dz = zSelector.run("Finalise PSF", true, true, true, null);
if (dz < 0) {
return;
}
zCentre = zSelector.getCentreSlice();
if (settings.getCropToZCentre()) {
finalPsf = finalPsf.cropToZCentre(zCentre);
// Back to 1-based index
zCentre = finalPsf.stackZCentre + 1;
}
// When click ok the background is subtracted from the PSF
// All pixels below the background are set to zero
// Apply a Tukey window to roll-off to zero at the outer pixels
ImageJUtils.showStatus("Windowing PSF");
final double[] wx = ImageWindow.tukeyEdge(finalPsf.maxx, settings.getWindow());
final double[] wz = ImageWindow.tukeyEdge(finalPsf.psf.length, settings.getWindow());
// Normalisation so the max intensity frame is one
final float[][] psf = finalPsf.psf;
final int maxz = psf.length;
final double[] sum = new double[maxz];
for (int z = 0; z < maxz; z++) {
sum[z] = applyWindow(psf[z], z, wx, wz, zSelector.background);
}
// Smooth the intensity
ImageJUtils.showStatus("Normalising PSF");
final Smoother smoother = zSelector.ssmoother;
final double[] ssum = smoother.smooth(sum).getDSmooth();
// Compute normalisation and apply.
SimpleArrayUtils.multiply(ssum, 1.0 / MathUtils.max(ssum));
for (int z = 0; z < psf.length; z++) {
if (sum[z] != 0) {
SimpleArrayUtils.multiply(psf[z], ssum[z] / sum[z]);
}
sum[z] = MathUtils.sum(psf[z]);
}
// Show the final intensity profile
final double[] slice = SimpleArrayUtils.newArray(maxz, 1, 1.0);
final Plot plot = new Plot(TITLE_SIGNAL, "Slice", "Signal");
final double[] range = MathUtils.limits(sum);
plot.setLimits(1, maxz, range[0], range[1]);
plot.setColor(Color.black);
plot.addPoints(slice, sum, Plot.LINE);
ImageJUtils.display(TITLE_SIGNAL, plot);
// Create a new extracted PSF and show
ImageJUtils.showStatus("Displaying PSF");
magnification = finalPsf.magnification;
finalPsf = new ExtractedPsf(psf, finalPsf.maxx, finalPsf.centre, magnification);
finalPsf.createProjections();
psfOut = finalPsf.show(TITLE_PSF, zCentre);
psfImp = psfOut[0];
// Add image info
final int imageCount = centres.length;
final ImagePSF.Builder imagePsf = ImagePsfHelper.create(zCentre, nmPerPixel / magnification, settings.getNmPerSlice() / magnification, imageCount, 0, createNote()).toBuilder();
// Add the CoM
// Find the XY centre around the z centre
final double[] com = getCentreOfMassXy(finalPsf.psf, finalPsf.maxx, finalPsf.maxy, zCentre - 1, settings.getComWindow(), getComXyBorder(finalPsf.maxx, finalPsf.maxy));
imagePsf.setXCentre(com[0]);
imagePsf.setYCentre(com[1]);
imagePsf.setZCentre(zCentre - 1);
psfImp.setProperty("Info", ImagePsfHelper.toString(imagePsf));
psfImp.setRoi(new OffsetPointRoi(com[0], com[1]));
psfImp.setSlice(zCentre);
psfImp.resetDisplayRange();
psfImp.updateAndDraw();
ImageJUtils.log("Final Centre-of-mass = %s,%s\n", rounder.toString(com[0]), rounder.toString(com[1]));
ImageJUtils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images\n", psfImp.getTitle(), zCentre, MathUtils.rounded(imagePsf.getPixelSize(), 3), MathUtils.rounded(imagePsf.getPixelDepth(), 3), imageCount);
if (settings.getOutputType() == OUTPUT_TYPE_CSPLINE) {
// Ask this again as it is important
// if (TextUtils.isNullOrEmpty(settings.getSplineFilename()))
// {
final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addFilenameField("Spline_filename", settings.getSplineFilename());
gd.showDialog(true);
if (gd.wasCanceled()) {
return;
}
settings.setSplineFilename(gd.getNextString());
// }
if (!TextUtils.isNullOrEmpty(settings.getSplineFilename())) {
// Save the result ...
IJ.showStatus("Creating cubic spline");
final CubicSplinePsf cubicSplinePsf = CubicSplineManager.createCubicSpline(imagePsf, psfImp.getImageStack(), settings.getSinglePrecision());
IJ.showStatus("Saving cubic spline");
CubicSplineManager.save(cubicSplinePsf, settings.getSplineFilename());
final String msg = "Spline saved to " + settings.getSplineFilename();
IJ.showStatus(msg);
IJ.log(msg);
// To leave the status message
return;
}
}
IJ.showStatus("");
}
use of uk.ac.sussex.gdsc.smlm.model.camera.CcdCameraModel in project GDSC-SMLM by aherbert.
the class Noise method createCameraModel.
private void createCameraModel() {
yAxisTitle = Y_AXIS_PHOTON;
switch(calibration.getCameraType()) {
case CCD:
cameraModel = new CcdCameraModel(calibration.getBias(), calibration.getCountPerPhoton());
break;
case EMCCD:
cameraModel = new EmCcdCameraModel(calibration.getBias(), calibration.getCountPerPhoton());
break;
case SCMOS:
cameraModel = CameraModelManager.load(calibration.getCameraModelName());
if (cameraModel == null) {
throw new IllegalStateException("No camera model for camera type: " + calibration.getCameraType());
}
cameraModel = PeakFit.cropCameraModel(cameraModel, IJImageSource.getBounds(imp), null, false);
// Store for next time
final Rectangle bounds = cameraModel.getBounds();
final int ox = bounds.x;
final int oy = bounds.y;
// Reset origin for filtering
if (ox != 0 || oy != 0) {
cameraModel = cameraModel.copy();
cameraModel.setOrigin(0, 0);
}
break;
case CAMERA_TYPE_NA:
case UNRECOGNIZED:
default:
cameraModel = new NullCameraModel();
yAxisTitle = Y_AXIS_COUNT;
break;
}
}
use of uk.ac.sussex.gdsc.smlm.model.camera.CcdCameraModel in project GDSC-SMLM by aherbert.
the class FitConfiguration method getCameraModel.
/**
* Gets the camera model.
*
* @return the camera model
* @throws IllegalStateException if no camera model exists for the camera type
*/
public CameraModel getCameraModel() {
if (cameraModel == null) {
final int value = getCameraTypeValue();
switch(value) {
case CameraType.CAMERA_TYPE_NA_VALUE:
// We can support this by doing nothing to pixels values
cameraModel = new NullCameraModel();
break;
case CameraType.CCD_VALUE:
case CameraType.EMCCD_VALUE:
final double bias = calibration.getBias();
final double gain = calibration.getCountPerPhoton();
final double variance = MathUtils.pow2(calibration.getReadNoise());
// This will throw an exception if the calibration is invalid.
cameraModel = (value == CameraType.EMCCD_VALUE) ? new EmCcdCameraModel(bias, gain, variance) : new CcdCameraModel(bias, gain, variance);
break;
case CameraType.SCMOS_VALUE:
default:
throw new IllegalStateException("No camera model for camera type: " + getCameraType());
}
}
return cameraModel;
}
Aggregations