gt: 		FloatProcessor floatP = ImageUtil.wrapGrid2D(inputGrid);
context: 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));
  