package edu.stanford.rsl.conrad.calibration;
import ij.ImagePlus;
import ij.gui.EllipseRoi;
import ij.gui.Overlay;
import ij.gui.PointRoi;
import ij.gui.Roi;
import ij.gui.TextRoi;
import ij.process.FloatProcessor;
import ij.process.FloatStatistics;
import java.awt.Color;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import Jama.EigenvalueDecomposition;
import Jama.Matrix;
import edu.mines.jtk.util.Array;
import edu.stanford.rsl.apps.gui.blobdetection.AutomaticMarkerDetectionWorker;
import edu.stanford.rsl.apps.gui.blobdetection.MarkerDetection;
import edu.stanford.rsl.apps.gui.blobdetection.MarkerDetectionWorker;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.filtering.FastRadialSymmetryTool;
import edu.stanford.rsl.conrad.filtering.NumericalDerivativeComputationTool;
import edu.stanford.rsl.conrad.geometry.General;
import edu.stanford.rsl.conrad.geometry.Projection;
import edu.stanford.rsl.conrad.geometry.Rotations;
import edu.stanford.rsl.conrad.geometry.Rotations.BasicAxis;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.trajectories.CircularTrajectory;
import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
import edu.stanford.rsl.conrad.numerics.DecompositionQR;
import edu.stanford.rsl.conrad.numerics.DecompositionSVD;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix.InversionType;
import edu.stanford.rsl.conrad.numerics.SimpleOperators;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.conrad.numerics.Solvers;
import edu.stanford.rsl.conrad.phantom.AbstractCalibrationPhantom;
import edu.stanford.rsl.conrad.phantom.AnalyticPhantom;
import edu.stanford.rsl.conrad.phantom.MathematicalPhantom;
import edu.stanford.rsl.conrad.phantom.RandomDistributionPhantom;
import edu.stanford.rsl.conrad.phantom.RandomizedHelixPhantom;
import edu.stanford.rsl.conrad.utils.CONRAD;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.ImageUtil;
import edu.stanford.rsl.conrad.utils.UserUtil;
import edu.stanford.rsl.hough.LineHoughSpace;
public class GeometricCalibration {
/**
* configuration, storing the ideal trajectory
*/
Configuration config;
/**
* storing ideal projection matrices derived from config
*
* @see config
*/
Projection[] pMatricesIdeal;
/**
* angle between principal components of ideal and detected beads
*/
double angle = 0.0;
/**
* rotation matrix to eliminate the angle between ideal and detected beads
*/
SimpleMatrix rotation;
/**
* detected beads
*/
ArrayList<PointND> beads2DUnsorted;
/**
* ideally projected beads
*/
ArrayList<PointND> beads2DIdeal;
/**
* formerly used for bead identification
*
* @deprecated
*/
ArrayList<PointND> beads3Dbackprojected;
/**
* formerly used for bead identification
*
* @deprecated
*/
ArrayList<PointND> beads3DReal;
/**
* correspondences between 2D and 3D coordinates of the beads
*/
ArrayList<CalibrationBead> cBeads;
/**
* bead IDs stored according to cBeads
*
* @see cBeads
*/
int[] ids;
/**
* storing back projection error for each projection
*/
ArrayList<Double> errors;
/**
* bead radius in px
*/
int beadRadius = 7;
/**
* value to adjust the bead detection
*/
double threshold = 0.000;
double largeBeadPercentileD = 0.999;
double smallBeadPercentileD = 0.9998;
double thresh1PercentileD = 0.000;
double thresh2PercentileD = 0.05;
/**
* minimal distance between to beads to be detected in px
*/
private double minimalDistanceBetweenTwoBeadsPx = 20;
/**
* number of projection
*/
int slice = 0;
/**
* calibration phantom name, inialized with RDP
*/
String phantomName = "Random Distribution Phantom";
/**
* calibration phantom
*/
AbstractCalibrationPhantom phantom;
/**
* storing the projection matrix for each projection
*/
ArrayList<SimpleMatrix> projectionMatrices = new ArrayList<SimpleMatrix>();
/**
* constructor, initializes phantom, config, rotation, erros and projection
* matrices
*
* @see phantom
* @see config
* @see rotation
* @see errors
* @see projectionMatrices
*/
public GeometricCalibration() {
initPhantom();
config = new Configuration();
rotation = SimpleMatrix.I_4;
errors = new ArrayList<Double>();
projectionMatrices = new ArrayList<SimpleMatrix>();
}
/**
* method, to change phantom in the GUI
*
* @param phantom
* , name of the phantom
*/
public void setPhantom(String phantom) {
this.phantomName = new String(phantom);
initPhantom();
}
/**
* method, to adjust slice in the GUI
*
* @param slice
*/
public void setSlice(int slice) {
this.slice = slice;
}
/**
* initializes phantom according to phantomName
*
* @see phantom
* @see phantomName
*/
private void initPhantom() {
if (phantomName.equals("Random Distribution Phantom")) {
phantom = new RandomDistributionPhantom();
} else if (phantomName.equals("Randomized Helix Phantom")) {
phantom = new RandomizedHelixPhantom();
} else {
phantom = new MathematicalPhantom();
}
phantom.init();
}
/**
* method to adjust detection in the GUI
*
* @param rad
* @param thresh
*/
public void setDetection(int rad, double thresh) {
beadRadius = rad;
threshold = thresh;
}
/**
*
* @param id
* @return 3D PointND containing x, y and z coordinates of the bead with id
* id
*/
private PointND get3DCoordinates(int id) {
return new PointND(phantom.getX(id), phantom.getY(id), phantom.getZ(id));
}
/**
* method to set config in the GUI
*
* @param config
* @see config
* @deprecated
*/
public void setConfig(String config) {
this.config = Configuration.loadConfiguration("d:\\Desktop\\" + config);
}
public void configure() {
config = new Configuration();
// pMatricesIdeal = ((CircularTrajectory) config.getGeometry())
// .getProjectionMatrices();
initPhantom();
rotation = SimpleMatrix.I_4;
}
/**
* method to initialize pMatricesIdeal after loading a configuration in the
* GUI
*/
public void reConfigure() {
pMatricesIdeal = ((CircularTrajectory) config.getGeometry())
.getProjectionMatrices();
}
/**
* Computes the projection matrices using a list of correspondences
*
* @param listCoords
* @return SimpleMatrix that maps world to image points
* @author Andreas Maier, Philipp Roser
*/
public static SimpleMatrix computePMatrix(
ArrayList<CalibrationBead> listCoords) {
Collections.sort(listCoords);
int N = listCoords.size();
// initialize measurement & dVector matrix
SimpleMatrix measurements = new SimpleMatrix(2 * N, 11);
SimpleMatrix dVector = new SimpleMatrix(2 * N, 1);
for (int i = 0; i < N; i++) {
CalibrationBead coords = listCoords.get(i);
// 1st row of M
measurements.setElementValue(i * 2, 0, coords.getX());
measurements.setElementValue(i * 2, 1, coords.getY());
measurements.setElementValue(i * 2, 2, coords.getZ());
measurements.setElementValue(i * 2, 3, 1);
measurements.setElementValue(i * 2, 4, 0);
measurements.setElementValue(i * 2, 5, 0);
measurements.setElementValue(i * 2, 6, 0);
measurements.setElementValue(i * 2, 7, 0);
measurements.setElementValue(i * 2, 8,
-coords.getX() * coords.getU());
measurements.setElementValue(i * 2, 9,
-coords.getY() * coords.getU());
measurements.setElementValue(i * 2, 10,
-coords.getZ() * coords.getU());
// 2nd row of M
measurements.setElementValue(i * 2 + 1, 0, 0);
measurements.setElementValue(i * 2 + 1, 1, 0);
measurements.setElementValue(i * 2 + 1, 2, 0);
measurements.setElementValue(i * 2 + 1, 3, 0);
measurements.setElementValue(i * 2 + 1, 4, coords.getX());
measurements.setElementValue(i * 2 + 1, 5, coords.getY());
measurements.setElementValue(i * 2 + 1, 6, coords.getZ());
measurements.setElementValue(i * 2 + 1, 7, 1);
measurements.setElementValue(i * 2 + 1, 8,
-coords.getX() * coords.getV());
measurements.setElementValue(i * 2 + 1, 9,
-coords.getY() * coords.getV());
measurements.setElementValue(i * 2 + 1, 10,
-coords.getZ() * coords.getV());
// 1st & 2nd row of vector d
dVector.setElementValue(i * 2, 0, coords.getU());
dVector.setElementValue(i * 2 + 1, 0, coords.getV());
}
// compute pVector (11*1)
DecompositionSVD svd = new DecompositionSVD(measurements);
SimpleMatrix pVector = SimpleOperators.multiplyMatrixProd(
svd.inverse(true), dVector);
// rewrite to Matrix format 3 x 4 (p23 == 1)
SimpleMatrix pMatrix = new SimpleMatrix(3, 4);
pMatrix.setElementValue(0, 0, 700 * pVector.getElement(0, 0));
pMatrix.setElementValue(0, 1, 700 * pVector.getElement(1, 0));
pMatrix.setElementValue(0, 2, 700 * pVector.getElement(2, 0));
pMatrix.setElementValue(0, 3, 700 * pVector.getElement(3, 0));
pMatrix.setElementValue(1, 0, 700 * pVector.getElement(4, 0));
pMatrix.setElementValue(1, 1, 700 * pVector.getElement(5, 0));
pMatrix.setElementValue(1, 2, 700 * pVector.getElement(6, 0));
pMatrix.setElementValue(1, 3, 700 * pVector.getElement(7, 0));
pMatrix.setElementValue(2, 0, 700 * pVector.getElement(8, 0));
pMatrix.setElementValue(2, 1, 700 * pVector.getElement(9, 0));
pMatrix.setElementValue(2, 2, 700 * pVector.getElement(10, 0));
pMatrix.setElementValue(2, 3, 700 * 1);
return pMatrix;
}
/**
* estimates threshold
*
* @param histogram
* @param stats
* @param percentile
* @return threshold
* @author Andreas Maier
*/
private float getHistogramThreshold(long[] histogram,
FloatStatistics stats, double percentile) {
int threshold = 0;
long count = 0;
for (threshold = 0; count < percentile * stats.pixelCount; threshold++) {
count += histogram[threshold];
}
return (float) (threshold * stats.binSize + stats.histMin);
}
/**
*
* @param linearHoughSpace
* @return
* @author Andreas Maier
*/
private int lineExtractBestLineCandidate(LineHoughSpace linearHoughSpace) {
FloatProcessor houghspaceAsFloatProcessor = (FloatProcessor) linearHoughSpace
.getImagePlus().getChannelProcessor();
FloatStatistics stats = new FloatStatistics(houghspaceAsFloatProcessor);
long[] histogram = stats.getHistogram();
float edgeThreshold = getHistogramThreshold(histogram, stats, 0.9999);
ArrayList<PointND> lineCandidate = General.extractCandidatePoints(
houghspaceAsFloatProcessor, edgeThreshold); // 245
ArrayList<PointND> lines = General.extractClusterCenter(lineCandidate,
20); // distance
return (int) lines.get(0).get(0);
}
/**
*
* @param currentImage
* @param derivativeTool
* @param rightBorder
* @param leftBorder
* @return
* @author Andreas Maier
*/
protected Roi[] detectHoughLines(ImagePlus currentImage,
NumericalDerivativeComputationTool derivativeTool, Roi rightBorder,
Roi leftBorder) {
// Compute derivative
double houghLineThresh1 = 100;
double houghLineThresh2 = 360;
Grid2D before = ImageUtil.wrapImageProcessor(currentImage
.getProcessor().toFloat(0, null));
// /VisualizationUtil.showGrid2D(before, "before");
// Configuration.getGlobalConfiguration().getGeometry()
// .setPixelDimensionX(0.3);
Grid2D derivative = derivativeTool.applyToolToImage(new Grid2D(before));
for (int j = 0; j < derivative.getHeight(); j++) {
for (int i = 0; i < derivative.getWidth(); i++) {
derivative.setAtIndex(i, j,
Math.abs(derivative.getPixelValue(i, j)));
}
}
// VisualizationUtil.showGrid2D(derivative, "Deriv");
FloatProcessor edgeImage = ImageUtil.wrapGrid2D(derivative);
FloatStatistics stats = new FloatStatistics(edgeImage);
long[] histogram = stats.getHistogram();
houghLineThresh1 = getHistogramThreshold(histogram, stats, 0.95);
houghLineThresh2 = getHistogramThreshold(histogram, stats, 0.99);
LineHoughSpace houghLineLeft = new LineHoughSpace(1.0, 3.0,
derivative.getWidth(), derivative.getHeight());
LineHoughSpace houghLineRight = new LineHoughSpace(1.0, 3.0,
derivative.getWidth(), derivative.getHeight());
for (int j = 0; j < derivative.getHeight(); j++) {
for (int i = 0; i < derivative.getWidth(); i++) {
double value = Math.abs(derivative.getPixelValue(i, j));
if (value > houghLineThresh1 && value < houghLineThresh2) { // thresholds
// for
// line
if ((i > leftBorder.getBounds().x)
&& (i < rightBorder.getBounds().x
+ rightBorder.getBounds().width)) // set
// x-axis
// boundaries
// for
// outer
// cylinder
// line
if (i < leftBorder.getBounds().x
+ leftBorder.getBounds().width) {
houghLineLeft.fill(i, j, 1.0); // cylinder boundary
// line
} else if (i > rightBorder.getBounds().x) {
houghLineRight.fill(i, j, 1.0); // cylinder boundary
// line
}
}
}
}
// houghLineLeft.getImagePlus().show();
// houghLineRight.getImagePlus().show();
Roi leftLineRoi;
Roi rightLineRoi;
try {
int leftLine = lineExtractBestLineCandidate(houghLineLeft);
leftLineRoi = new Roi(leftLine, 0, 1, currentImage.getHeight());
leftLineRoi.setName("LeftLineRoi");
int rightLine = lineExtractBestLineCandidate(houghLineRight);
rightLineRoi = new Roi(rightLine, 0, 1, currentImage.getHeight());
rightLineRoi.setName("RightLineRoi");
} catch (Exception e) {
leftLineRoi = null;
rightLineRoi = null;
}
return new Roi[] { leftLineRoi, rightLineRoi };
// this happens in the GUI
// overlayLines();
// this.revalidate();
// this.repaint();
}
/**
*
* @param inputGrid
* @param radius
* @param percentile
* @param distance
* @param leftLineRoi
* @param rightLineRoi
* @return
* @author Andreas Maier
*/
protected ArrayList<PointND> extractBeads(Grid2D inputGrid, double radius,
double percentile, double distance, Roi leftLineRoi,
Roi rightLineRoi) {
FastRadialSymmetryTool frst = new FastRadialSymmetryTool(new double[] {
radius, radius + 2 }, 3, 0, null, 0);
inputGrid = frst.applyToolToImage(new Grid2D(inputGrid));
FloatProcessor floatP = ImageUtil.wrapGrid2D(inputGrid);
FloatStatistics stats = new FloatStatistics(floatP);
long[] histogram = stats.getHistogram();
double thresh1 = getHistogramThreshold(histogram, stats, 0.000);
double thresh2 = getHistogramThreshold(histogram, stats, 0.05);
for (int j = 0; j < floatP.getHeight(); j++) {
for (int i = 0; i < floatP.getWidth(); i++) {
double value = (floatP.getPixelValue(i, j));
if (value > thresh1 && value < thresh2) { // thresholds for line
if ((i > leftLineRoi.getBounds().x)
&& (i < rightLineRoi.getBounds().x)) // set x-axis
// boundaries
// for outer
// cylinder
// line
floatP.putPixelValue(i, j, -floatP.getPixelValue(i, j));
else {
floatP.putPixelValue(i, j, 0.0);
}
} else
floatP.putPixelValue(i, j, 0.0);
}
}
stats = new FloatStatistics(floatP);
histogram = stats.getHistogram();
thresh1 = getHistogramThreshold(histogram, stats, percentile);
ArrayList<PointND> candidate = General.extractCandidatePoints(floatP,
thresh1);
// filter with min distance
return General.extractClusterCenter(candidate, distance, false); // distance
}
/**
*
* @param tmpBead
* @param pMatrix
* @return
* @author Andreas Maier
*/
public static double[] compute2DCoordinates(CalibrationBead tmpBead,
SimpleMatrix pMatrix) {
SimpleMatrix Q3D = new SimpleMatrix(4, 1);
Q3D.setElementValue(0, 0, tmpBead.getX());
Q3D.setElementValue(1, 0, tmpBead.getY());
Q3D.setElementValue(2, 0, tmpBead.getZ());
Q3D.setElementValue(3, 0, 1);
SimpleMatrix Q2D = SimpleOperators.multiplyMatrixProd(pMatrix, Q3D);
return new double[] { Q2D.getElement(0, 0) / Q2D.getElement(2, 0),
Q2D.getElement(1, 0) / Q2D.getElement(2, 0) };
}
/**
*
* @param point3D
* @param pMatrix
* @param slice
* @return
* @author Martin Berger
*/
private double[] compute2Dfrom3D(double[] point3D, Projection pMatrix,
int slice) {
// Compute coordinates in projection data.
SimpleVector homogeneousPoint = SimpleOperators.multiply(pMatrix
.computeP(), new SimpleVector(point3D[0], point3D[1],
point3D[2], 1));
// Do forward projection to 2D coordinates
double coordU = homogeneousPoint.getElement(0)
/ homogeneousPoint.getElement(2);
double coordV = homogeneousPoint.getElement(1)
/ homogeneousPoint.getElement(2);
return new double[] { coordU, coordV, (double) slice };
}
/**
*
* @param point2D
* @param p
* @return
* @author Martin Berger
*/
private PointND compute3Dfrom2D(PointND point2D, Projection p) {
SimpleVector homP = SimpleOperators.multiply(
p.computeP().inverse(InversionType.INVERT_SVD),
new SimpleVector(point2D.get(0), point2D.get(1), 1));
double x = homP.getElement(0) / homP.getElement(3);
double y = homP.getElement(1) / homP.getElement(3);
double z = homP.getElement(2) / homP.getElement(3);
return new PointND(x, y, z);
}
protected double measureDistance(int i, PointND p1, PointND p2) {
if (i < 1 || i > 3) {
throw new IllegalArgumentException();
} else if (i == 1) {
return Math.min(
p1.euclideanDistance(p2),
Math.abs(p1.get(p1.getDimension() - 1)
- p2.get(p2.getDimension() - 1)));
} else {
return p1.euclideanDistance(p2);
}
}
/**
* method only for test issues
*
* @deprecated
*/
protected void computeBeadIDsBack() {
// initialize array storing all euclidean distances between real and
// real 3D positions
double[][] distMeasure = new double[beads3Dbackprojected.size()][beads3DReal
.size()];
int r = 0;
int i = 0;
// actual measurements
for (PointND real : beads3Dbackprojected) {
for (PointND ideal : beads3DReal) {
double dist = measureDistance(2, real, ideal);
distMeasure[r][i] = dist;
i++;
}
i = 0;
r++;
}
// initalize arrays for potential minima and corresponding bead ids
double[] potentialMinima = new double[beads3Dbackprojected.size()];
int[] potentialIds = new int[beads3Dbackprojected.size()];
// find minimum for each real 2D position detected
for (r = 0; r < beads3Dbackprojected.size(); r++) {
double min = Array.max(distMeasure[r]);
for (i = 0; i < beads3DReal.size(); i++) {
if (distMeasure[r][i] < min) {
min = distMeasure[r][i];
potentialIds[r] = i;
}
}
potentialMinima[r] = min;
}
// initialize list and arrays for actual calibration
cBeads = new ArrayList<CalibrationBead>();
CalibrationBead[] minima = new CalibrationBead[beads3Dbackprojected
.size()];
ids = new int[beads3Dbackprojected.size()];
// find n = 8 correspondences with smallest euclidean distance
double max = Array.max(potentialMinima);
for (int n = 0; n < beads3Dbackprojected.size(); n++) {
double min = Array.max(potentialMinima);
int id = 0;
int pos = 0;
for (r = 0; r < potentialMinima.length; r++) {
if (potentialMinima[r] <= min) {
min = potentialMinima[r];
id = potentialIds[r];
pos = r;
}
}
potentialMinima[pos] = max;
ids[n] = id;
minima[n] = new CalibrationBead(beads2DUnsorted.get(pos).get(0),
beads2DUnsorted.get(pos).get(1));
phantom.setBeadCoordinates(minima[n], ids[n]);
cBeads.add(minima[n]);
}
}
// only for test issues
/**
* method only for test issues
*
* @deprecated
*/
protected void computeBeadIDsIdeal() {
// initialize array storing all euclidean distances between real and
// ideal 2D positions
double[][] distMeasure = new double[beads2DIdeal.size()][beads2DIdeal
.size()];
int r = 0;
int i = 0;
// actual measurements
for (PointND real : beads2DIdeal) {
for (PointND ideal : beads2DIdeal) {
double dist = measureDistance(2, real, ideal);
distMeasure[r][i] = dist;
i++;
}
i = 0;
r++;
}
// initalize arrays for potential minima and corresponding bead ids
double[] potentialMinima = new double[beads2DIdeal.size()];
int[] potentialIds = new int[beads2DIdeal.size()];
// find minimum for each real 2D position detected
for (r = 0; r < beads2DIdeal.size(); r++) {
double min = 1000;
for (i = 0; i < beads2DIdeal.size(); i++) {
if (distMeasure[r][i] < min) {
min = distMeasure[r][i];
potentialIds[r] = i;
}
}
potentialMinima[r] = min;
}
// initialize list and arrays for actual calibration
cBeads = new ArrayList<CalibrationBead>();
CalibrationBead[] minima = new CalibrationBead[beads2DIdeal.size()];
ids = new int[beads2DIdeal.size()];
// find n = 8 correspondences with smallest euclidean distance
double max = 1000;
for (int n = 0; n < beads2DIdeal.size(); n++) {
double min = 1000;
int id = 0;
int pos = 0;
for (r = 0; r < potentialMinima.length; r++) {
if (potentialMinima[r] <= min) {
min = potentialMinima[r];
id = potentialIds[r];
pos = r;
}
}
potentialMinima[pos] = max;
ids[n] = id;
minima[n] = new CalibrationBead(beads2DIdeal.get(pos).get(0),
beads2DIdeal.get(pos).get(1));
phantom.setBeadCoordinates(minima[n], ids[n]);
cBeads.add(minima[n]);
}
}
/**
* estimates rotation between ideally projected and actually detected beads
* using PCA
*
* @see PrincipalComponentAnalysis2D
*/
protected void estimateRotation() {
PrincipalComponentAnalysis2D real = new PrincipalComponentAnalysis2D(
beads2DUnsorted);
PrincipalComponentAnalysis2D ideal = new PrincipalComponentAnalysis2D(
beads2DIdeal);
System.out.println("Real EV:" + real.stringEigenvectors());
System.out.println("Ideal EV:" + ideal.stringEigenvectors());
Matrix realV = real.getEigenvectors()[0]
.plus(real.getEigenvectors()[1]);
Matrix idealV = ideal.getEigenvectors()[0]
.plus(ideal.getEigenvectors()[1]);
double scalar = Math.abs(realV.get(0, 0) * idealV.get(0, 0))
+ Math.abs(realV.get(1, 0) * idealV.get(1, 0));
angle = Math.acos(scalar / (realV.norm2() * idealV.norm2()));
System.out.println("Angle = " + angle * 180.0 / Math.PI);
double s = Math.sin(angle);
double c = Math.cos(angle);
rotation = new SimpleMatrix(
new double[][] { { 1.0, 0.0, 0.0, 0.0 }, { 0.0, c, -s, 0.0 },
{ 0.0, s, c, 0.0 }, { 0.0, 0.0, 0.0, 1.0 } });
}
// compares real and ideal 2D positions to establish most probable
// correspondences
/**
* establishes correspondences between world and image points by linking
* ideal and detected beads via the euclidean distance
*/
protected void computeBeadIDsForward() {
beads2DIdeal = new ArrayList<PointND>();
for (int i = 0; i < numberOfBeads(); i++) {
double[] coord = compute2Dfrom3D(
new double[] { get3DCoordinates(i).get(0),
get3DCoordinates(i).get(1),
get3DCoordinates(i).get(2) },
new Projection(SimpleOperators.multiplyMatrixProd(
pMatricesIdeal[slice].computeP(), rotation)), slice);
beads2DIdeal.add(i, new PointND(coord[0], coord[1]));
}
// initialize array storing all euclidean distances between real and
// ideal 2D positions
double[][] distMeasure = new double[beads2DUnsorted.size()][beads2DIdeal
.size()];
int r = 0;
int i = 0;
// actual measurements
for (PointND real : beads2DUnsorted) {
for (PointND ideal : beads2DIdeal) {
double dist = real.euclideanDistance(ideal);
distMeasure[r][i] = dist;
i++;
}
i = 0;
r++;
}
// initalize arrays for potential minima and corresponding bead ids
double[] potentialMinima = new double[beads2DUnsorted.size()];
int[] potentialIds = new int[beads2DUnsorted.size()];
// find minimum for each real 2D position detected
for (r = 0; r < beads2DUnsorted.size(); r++) {
double min = Array.max(distMeasure[r]);
for (i = 0; i < beads2DIdeal.size(); i++) {
if (distMeasure[r][i] < min) {
min = distMeasure[r][i];
potentialIds[r] = i;
}
}
potentialMinima[r] = min;
}
// initialize list and arrays for actual calibration
cBeads = new ArrayList<CalibrationBead>();
CalibrationBead[] minima = new CalibrationBead[beads2DUnsorted.size()];
ids = new int[beads2DUnsorted.size()];
// find n = 8 correspondences with smallest euclidean distance
double max = Array.max(potentialMinima) + 1;
for (int n = 0; n < beads2DUnsorted.size(); n++) {
double min = Array.max(potentialMinima) + 1;
int id = 0;
int pos = 0;
for (r = 0; r < potentialMinima.length; r++) {
if (potentialMinima[r] <= min) {
min = potentialMinima[r];
id = potentialIds[r];
pos = r;
}
}
potentialMinima[pos] = max;
ids[n] = id;
CalibrationBead cb = new CalibrationBead(beads2DUnsorted.get(pos)
.get(0), beads2DUnsorted.get(pos).get(1));
minima[n] = new CalibrationBead(beads2DUnsorted.get(pos).get(0),
beads2DUnsorted.get(pos).get(1));
phantom.setBeadCoordinates(minima[n], ids[n]);
phantom.setBeadCoordinates(cb, id);
cBeads.add(cb);
}
}
/**
* contains all image points of all frames
*/
ArrayList<ArrayList<PointND>> imagePoints;
/**
* reference points to sort newly detected beads
*/
ArrayList<PointND> referencePoints;
/**
* potential bead ids to be removed
*/
ArrayList<Integer> remove;
int loc = 0;
PointND[][] image;
PointND[] reference;
/**
* Finds initial points for Factorization. Currently configured for
* debugging.
*
* @param pointSelector
* , if true initial points are estimated manually
* @param im
* @param leftLineRoi
* @param rightLineRoi
*/
protected void getInitialPoints(boolean pointSelector, ImagePlus im,
Roi leftLineRoi, Roi rightLineRoi) {
referencePoints = new ArrayList<PointND>();
imagePoints = new ArrayList<ArrayList<PointND>>();
remove = new ArrayList<Integer>();
if (pointSelector) {
} else {
// detectBeads(im, leftLineRoi, rightLineRoi);
// referencePoints = (ArrayList<PointND>) beads2DUnsorted.clone();
// imagePoints.add(0, (ArrayList<PointND>) referencePoints.clone());
idealBeads();
imagePoints.add(new ArrayList<PointND>(beads2DIdeal));
}
}
/**
* method for test issues
*
* @param pointSelector
* @param im
* @param leftLineRoi
* @param rightLineRoi
* @deprecated
*/
protected void getInitialPointsArray(boolean pointSelector, ImagePlus im,
Roi leftLineRoi, Roi rightLineRoi) {
reference = new PointND[numberOfBeads()];
image = new PointND[im.getImageStackSize()][numberOfBeads()];
remove = new ArrayList<Integer>();
if (pointSelector) {
} else {
detectBeads(im, leftLineRoi, rightLineRoi);
for (int i = 0; i < reference.length; i++) {
reference[i] = beads2DUnsorted.get(i);
}
image[0] = reference;
}
}
/**
* Establishes correspondences via linking detected beads to former
* established correspondences
*
* @param im
* @param leftLineRoi
* @param rightLineRoi
*/
protected void findPointsNaive(ImagePlus im, Roi leftLineRoi,
Roi rightLineRoi) {
ArrayList<CalibrationBead> newC = new ArrayList<CalibrationBead>();
for (CalibrationBead cb : cBeads) {
boolean found = false;
int pos = -1;
PointND ref = new PointND(cb.getU(), cb.getV());
double thresh = 15.0;
for (PointND det : beads2DUnsorted) {
CalibrationBead newBead = new CalibrationBead(0.0, 0.0);
newC.add(newBead);
double dist = ref.euclideanDistance(det);
if (dist < thresh) {
thresh = dist;
newBead.setU(det.get(0));
newBead.setV(det.get(1));
newBead.setX(cb.getX());
newBead.setY(cb.getY());
newBead.setZ(cb.getZ());
found = true;
newC.add(newBead);
pos = beads2DUnsorted.indexOf(det);
}
}
if (!found) {
newC.add(cb);
} else {
beads2DUnsorted.remove(pos);
}
}
cBeads = new ArrayList<CalibrationBead>(newC);
}
/**
* establishes point sets for Factorization. Configured for debugging
*
* @param im
* @param leftLineRoi
* @param rightLineRoi
*/
protected void findPoints(ImagePlus im, Roi leftLineRoi, Roi rightLineRoi) {
beads2DIdeal = new ArrayList<PointND>();
idealBeads();
ArrayList<PointND> newL = new ArrayList<PointND>(beads2DIdeal);
imagePoints.add(newL);
/*
* beads2DUnsorted = new ArrayList<PointND>(); detectBeads(im,
* leftLineRoi, rightLineRoi); ArrayList<PointND> list = new
* ArrayList<PointND>(numberOfBeads()); int size = 0;
*
* for (PointND ref : referencePoints) {
*
* boolean found = false; int pos = referencePoints.indexOf(ref);
*
* if (ref == null) { continue; }
*
* for (PointND det : beads2DUnsorted) {
*
* if (Math.abs(ref.get(1) - det.get(1)) < 15.0) { if
* (list.contains(det)) { if (found &&
* list.get(pos).euclideanDistance(ref) > det .euclideanDistance(ref)) {
* list.add(pos, det); } } if (!found && !list.contains(det)) {
* list.add(pos, det); found = true; size++; } } }
*
* if (!found) { list.add(pos, referencePoints.get(pos)); remove.add(new
* Integer(pos)); }
*
* } System.out.println("Slice: " + slice + ", beads detected: " +
* beads2DUnsorted.size() + ", beads found: " + size);
* System.out.println(list); imagePoints.add(list); referencePoints =
* list;
*/
}
/**
* still testing
*/
protected void removePoints() {
remove = new ArrayList<Integer>();
for (Integer i : remove) {
for (ArrayList<PointND> list : imagePoints) {
list.remove(remove.get(i).intValue());
}
}
}
/**
* computes Factorization
*/
protected void factorize() {
Factorization fac = new Factorization(imagePoints);
projectionMatrices = fac.getProjectionMatrices();
}
/**
*
* @return number of beads of the used phantom
*/
private int numberOfBeads() {
return phantom.getNumberOfBeads();
}
// sorted IDs
/**
* computes ideal beads
*/
protected void idealBeads() {
beads2DIdeal = new ArrayList<PointND>();
for (int i = 0; i < numberOfBeads(); i++) {
double[] coord = compute2Dfrom3D(new double[] {
get3DCoordinates(i).get(0), get3DCoordinates(i).get(1),
get3DCoordinates(i).get(2) }, pMatricesIdeal[slice], slice);
beads2DIdeal.add(i, new PointND(coord[0], coord[1]));
}
}
// unsorted IDs
/**
* detect beads using FRST
*
* @param im
* @param leftLineRoi
* @param rightLineRoi
*/
protected void detectBeads(ImagePlus im, Roi leftLineRoi, Roi rightLineRoi) {
beads2DUnsorted = new ArrayList<PointND>();
beads2DUnsorted = extractBeads(ImageUtil.wrapImageProcessor(im
.getProcessor().toFloat(0, null)), beadRadius,
smallBeadPercentileD, minimalDistanceBetweenTwoBeadsPx,
leftLineRoi, rightLineRoi);
}
/**
* computes backprojected beads
*
* @deprecated
*/
protected void backprojectBeads() {
beads3Dbackprojected = new ArrayList<PointND>();
int i = 0;
for (PointND p2D : beads2DUnsorted) {
beads3Dbackprojected.add(i,
compute3Dfrom2D(p2D, pMatricesIdeal[slice]));
i++;
}
}
/**
* provides 3D coordinates of the beads
*/
protected void real3DBeads() {
beads3DReal = new ArrayList<PointND>();
for (int i = 0; i < numberOfBeads(); i++) {
beads3DReal.add(i, get3DCoordinates(i));
}
}
// created only for test issues
/**
* only for test issues
*
* @deprecated
* @return
*/
protected SimpleMatrix computePMatrixIdeal() {
cBeads = new ArrayList<CalibrationBead>();
for (int i = 0; i < beads3DReal.size(); i++) {
CalibrationBead cb = new CalibrationBead(
beads2DIdeal.get(i).get(0), beads2DIdeal.get(i).get(1));
cb.setX(beads3DReal.get(i).get(0));
cb.setY(beads3DReal.get(i).get(1));
cb.setZ(beads3DReal.get(i).get(2));
cBeads.add(cb);
}
return computePMatrix(cBeads);
}
/**
* calibration using factorization, not working completely yet
*
* @param sliceNumber
* @param overlay
* @return back projection error
* @author Andreas Maier, Philipp Roser
*/
public double calibrateF(int sliceNumber, Overlay overlay) {
SimpleMatrix pMatrix = projectionMatrices.get(sliceNumber);
System.out.println("slice" + sliceNumber + ": " + pMatrix);
double backprojectionError = 0;
for (int i = 0; i < numberOfBeads(); i++) {
CalibrationBead tmpBead = new CalibrationBead(0, 0);
phantom.setBeadCoordinates(tmpBead, i);
double[] coord2d = GeometricCalibration.compute2DCoordinates(
tmpBead, pMatrix);
for (PointND b : beads2DIdeal) {
PointND p = new PointND(coord2d[0], coord2d[1]);
if (b.euclideanDistance(p) < CONRAD.FLOAT_EPSILON) {
backprojectionError += b.euclideanDistance(p);
}
}
// mark projected bead in image.
overlay.add(new PointRoi(coord2d[0], coord2d[1]));
NumberFormat nf = NumberFormat.getInstance();
nf.setMaximumFractionDigits(2);
nf.setMinimumFractionDigits(2);
TextRoi text = new TextRoi(coord2d[0], coord2d[1] - 20, "Bead " + i);
overlay.add(text);
}
backprojectionError /= beads2DUnsorted.size();
// if(backprojectionError > 1.0){
// System.out.println("Error: "+ backprojectionError);
// }
// currentImage.setOverlay(overlay);
// config.getGeometry().getProjectionMatrices()[sliceNumber] = new
// Projection(
// new SimpleMatrix(pMatrix));
/**
* Focal spot Sw = -inverse of M * P3
*/
SimpleMatrix M = new SimpleMatrix(3, 3);
SimpleMatrix Sw = new SimpleMatrix(3, 1);
SimpleMatrix P3 = new SimpleMatrix(3, 1);
M.setElementValue(0, 0, pMatrix.getElement(0, 0));
M.setElementValue(0, 1, pMatrix.getElement(0, 1));
M.setElementValue(0, 2, pMatrix.getElement(0, 2));
M.setElementValue(1, 0, pMatrix.getElement(1, 0));
M.setElementValue(1, 1, pMatrix.getElement(1, 1));
M.setElementValue(1, 2, pMatrix.getElement(1, 2));
M.setElementValue(2, 0, pMatrix.getElement(2, 0));
M.setElementValue(2, 1, pMatrix.getElement(2, 1));
M.setElementValue(2, 2, pMatrix.getElement(2, 2));
P3.setElementValue(0, 0, pMatrix.getElement(0, 3));
P3.setElementValue(1, 0, pMatrix.getElement(1, 3));
P3.setElementValue(2, 0, pMatrix.getElement(2, 3));
Sw = SimpleOperators.multiplyMatrixProd(
M.inverse(SimpleMatrix.InversionType.INVERT_QR), P3)
.multipliedBy(-1.0);
// primary angle
double priAngle = Math.atan(Math.abs(Sw.getElement(1, 0)
/ Sw.getElement(0, 0)));
if (Sw.getElement(0, 0) < 0 && Sw.getElement(1, 0) >= 0)
priAngle = Math.PI - priAngle;// in 2nd quadrant
else if (Sw.getElement(0, 0) < 0 && Sw.getElement(1, 0) < 0)
priAngle = Math.PI + priAngle;// in 3rd quadrant
else if (Sw.getElement(0, 0) >= 0 && Sw.getElement(1, 0) < 0)
priAngle = 2 * Math.PI - priAngle;// in 4th quadrant
if (priAngle >= Math.PI)
priAngle = -(Math.PI * 2 - priAngle);
priAngle *= 180.0 / Math.PI;
// secondary angle
double sndAngle = Math.atan(Sw.getElement(2, 0)
/ Math.sqrt(Math.pow(Sw.getElement(0, 0), 2)
+ Math.pow(Sw.getElement(1, 0), 2)));
sndAngle *= 180.0 / Math.PI;
// config.getGeometry().getPrimaryAngles()[sliceNumber] = priAngle;
// config.getGeometry().getSecondaryAngles()[sliceNumber] = sndAngle;
// this.revalidate();
// this.repaint();
return backprojectionError;
}
/**
* calibration using correspondences
*
* @param sliceNumber
* @param overlay
* @return back projection error
* @author Andreas Maier, Philipp Roser
*/
public double calibrate(int sliceNumber, Overlay overlay) {
SimpleMatrix pMatrix = GeometricCalibration.computePMatrix(cBeads);
projectionMatrices.add(pMatrix);
double backprojectionError = 0;
for (int i = 0; i < numberOfBeads(); i++) {
CalibrationBead tmpBead = new CalibrationBead(0, 0);
phantom.setBeadCoordinates(tmpBead, i);
double[] coord2d = GeometricCalibration.compute2DCoordinates(
tmpBead, pMatrix);
// backprojectionError += beads2DIdeal.get(i).euclideanDistance(new
// PointND(coord2d[0], coord2d[1]));
for (CalibrationBead b : cBeads) {
PointND threeDPoint = new PointND(b.getX(), b.getY(), b.getZ());
PointND threeDPoint2 = new PointND(tmpBead.getX(),
tmpBead.getY(), tmpBead.getZ());
if (threeDPoint.euclideanDistance(threeDPoint2) < CONRAD.FLOAT_EPSILON) {
PointND twoDPoint = new PointND(b.getU(), b.getV());
PointND twoDPoint2 = new PointND(coord2d[0], coord2d[1]);
backprojectionError += twoDPoint
.euclideanDistance(twoDPoint2);
}
}
// mark projected bead in image.
overlay.add(new PointRoi(coord2d[0], coord2d[1]));
NumberFormat nf = NumberFormat.getInstance();
nf.setMaximumFractionDigits(2);
nf.setMinimumFractionDigits(2);
TextRoi text = new TextRoi(coord2d[0], coord2d[1] - 20, "Bead " + i);
overlay.add(text);
}
backprojectionError /= cBeads.size();
// if(backprojectionError > 1.0){
// System.out.println("Error: "+ backprojectionError);
// }
// currentImage.setOverlay(overlay);
// config.getGeometry().getProjectionMatrices()[sliceNumber] = new
// Projection(
// new SimpleMatrix(pMatrix));
/**
* Focal spot Sw = -inverse of M * P3
*/
SimpleMatrix M = new SimpleMatrix(3, 3);
SimpleMatrix Sw = new SimpleMatrix(3, 1);
SimpleMatrix P3 = new SimpleMatrix(3, 1);
M.setElementValue(0, 0, pMatrix.getElement(0, 0));
M.setElementValue(0, 1, pMatrix.getElement(0, 1));
M.setElementValue(0, 2, pMatrix.getElement(0, 2));
M.setElementValue(1, 0, pMatrix.getElement(1, 0));
M.setElementValue(1, 1, pMatrix.getElement(1, 1));
M.setElementValue(1, 2, pMatrix.getElement(1, 2));
M.setElementValue(2, 0, pMatrix.getElement(2, 0));
M.setElementValue(2, 1, pMatrix.getElement(2, 1));
M.setElementValue(2, 2, pMatrix.getElement(2, 2));
P3.setElementValue(0, 0, pMatrix.getElement(0, 3));
P3.setElementValue(1, 0, pMatrix.getElement(1, 3));
P3.setElementValue(2, 0, pMatrix.getElement(2, 3));
Sw = SimpleOperators.multiplyMatrixProd(
M.inverse(SimpleMatrix.InversionType.INVERT_QR), P3)
.multipliedBy(-1.0);
// primary angle
double priAngle = Math.atan(Math.abs(Sw.getElement(1, 0)
/ Sw.getElement(0, 0)));
if (Sw.getElement(0, 0) < 0 && Sw.getElement(1, 0) >= 0)
priAngle = Math.PI - priAngle;// in 2nd quadrant
else if (Sw.getElement(0, 0) < 0 && Sw.getElement(1, 0) < 0)
priAngle = Math.PI + priAngle;// in 3rd quadrant
else if (Sw.getElement(0, 0) >= 0 && Sw.getElement(1, 0) < 0)
priAngle = 2 * Math.PI - priAngle;// in 4th quadrant
if (priAngle >= Math.PI)
priAngle = -(Math.PI * 2 - priAngle);
priAngle *= 180.0 / Math.PI;
// secondary angle
double sndAngle = Math.atan(Sw.getElement(2, 0)
/ Math.sqrt(Math.pow(Sw.getElement(0, 0), 2)
+ Math.pow(Sw.getElement(1, 0), 2)));
sndAngle *= 180.0 / Math.PI;
// config.getGeometry().getPrimaryAngles()[sliceNumber] = priAngle;
// config.getGeometry().getSecondaryAngles()[sliceNumber] = sndAngle;
// this.revalidate();
// this.repaint();
if (slice < errors.size()) {
if (errors.get(slice) > backprojectionError) {
errors.remove(slice);
errors.add(new Double(backprojectionError));
}
} else {
errors.add(new Double(backprojectionError));
}
return backprojectionError;
}
}