content: package edu.stanford.rsl.tutorial.dmip;
  
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.conrad.utils.ImageUtil;
  import edu.stanford.rsl.conrad.utils.VisualizationUtil;
  import ij.IJ;
  import ij.ImageJ;
  import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.numerics.SimpleOperators;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  import edu.stanford.rsl.conrad.numerics.DecompositionSVD;
  
  public class ImageUndistortion{
  
  	
  	public void doImageUndstortion(){
  
  		/////////////////////
  		// 1.Preprocessing //
  		/////////////////////
  		// In the preprocessing part, an artificial distortion field is generated.
  		// With this field, a distorted image is generated.
  		
  		
  		// 1. Load undistorted image
  		
  		// TODO : adapt the paths
  		int caseNo = 0;
  		String filename = "C:/StanfordRepo/CONRAD/src/edu/stanford/rsl/tutorial/dmip/frame32.jpg";
  		
  		if(caseNo == 0)
  		{
  			filename = "C:/StanfordRepo/CONRAD/src/edu/stanford/rsl/tutorial/dmip/frame32.jpg";
  		}
  		else if(caseNo == 1)
  		{
  			filename = "C:/StanfordRepo/CONRAD/src/edu/stanford/rsl/tutorial/dmip/undistorted.jpg";
  		}else if(caseNo == 2)
  		{
  			filename = "C:/StanfordRepo/CONRAD/src/edu/stanford/rsl/tutorial/dmip/frame90.jpg";
  		}
  				
  		Grid2D image = ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0);
  		
  		
  		// 2. Normalize intensity values to [0,1]
  		float max = NumericPointwiseOperators.max(image);
  		float min = NumericPointwiseOperators.min(image);
  		
  		for(int i = 0; i < image.getWidth(); i++)
  		{
  			for(int j = 0; j < image.getHeight(); j++)
  			{
  				float scaledValue = (image.getAtIndex(i, j) - min) / (max-min);
  				image.setAtIndex(i, j, scaledValue);
  			}
  		}
  		image.show("Input Image");
  		
  		
  		// 3. Make the image quadratic
  		int h = image.getHeight();
  		int w = image.getWidth();
  		
  		int imSize = Math.min(h, w);
  		
  		Grid2D quadraticImage = new Grid2D(imSize, imSize); 
  		
  		for(int i  = 0; i < image.getWidth(); i++)
  		{
  			for(int j = 0; j < image.getHeight(); j++)
  			{
  				quadraticImage.setAtIndex(i, j, image.getPixelValue(i, j));
  			}
  		}
  		
  		quadraticImage.show("Quadratic Input Image");
  		
  		
  		// 4. Generate a grid to sample the image
  		// undistorted coordinates X, Y
  		Grid2D X = new Grid2D(imSize, imSize);
  		Grid2D Y = new Grid2D(imSize, imSize);
  		
  		for(int i = 0; i < X.getWidth(); i++)
  		{
  			for(int j = 0; j < X.getHeight(); j++)
  			{
  				X.setAtIndex(i, j, i);
  				Y.setAtIndex(i, j, j);
  			}
  		}
  		
  		// 5. Create an artificial elliptical distortion field (R)
  		// a: spread among x-direction
  		// b: spread among y-direction
  		// d: (d/2) = maximal value at the radius boundary
  		// NumericPointwiseOperators.subtractBy(R, (float) (maxR * 0.5)):
  		// 		shifting the positive range to half positive/half negative range
  		Grid2D R = new Grid2D(imSize, imSize);
  		float a = 3;
  		float b = 9; 
  		float d = 12;
  		
  		int half = imSize / 2;
  		
  		for(int i = 0; i < R.getWidth(); i++)
  		{
  			for(int j = 0; j < R.getHeight(); j++)
  			{
  				float r = (float) (d * Math.sqrt( a * Math.pow(((X.getAtIndex(i, j) - half) / imSize) , 2)  + b * Math.pow(((Y.getAtIndex(i, j) - half) / imSize) , 2) ) ) ;
  				R.setAtIndex(i, j, r);
  			}
  		}
  		
  		R.show("Distortion Field");
  		Grid2D shiftedR = new Grid2D(R);
  		float maxR = NumericPointwiseOperators.max(R);
  		NumericPointwiseOperators.subtractBy(shiftedR, (float) (maxR * 0.5));
  		
  		shiftedR.show("Shifted Distortion Field");
  		
  		// 5. Create the distorted image coordinates:
  		// distorted image coordinates = undistorted image points + artificial distortion field
  		
  		// distorted coordinates Xd, Yd
  		Grid2D Xd = new Grid2D(X);
  		Grid2D Yd = new Grid2D(Y);
  		
  		NumericPointwiseOperators.addBy(Xd, shiftedR);
  		NumericPointwiseOperators.addBy(Yd, shiftedR);
  		
  		// 6. Create the distorted image
  		// Resample the original input image at the distorted image points (XD,YD) to create artificial distortion
  		
  		Grid2D distortedImage = new Grid2D(imSize, imSize);
  		for(int i = 0; i < R.getWidth(); i++)
  		{
  			for(int j = 0; j < R.getHeight(); j++)
  			{
  				float distortedValue = InterpolationOperators.interpolateLinear(image, Xd.getAtIndex(i, j), Yd.getAtIndex(i, j));
  				distortedImage.setAtIndex(i, j, distortedValue);
  			}
  		}
  		
  		distortedImage.show("Distorted Image");
  				
  		///////////////////////////////////
  		// Image Undistortion - Workflow //
  		///////////////////////////////////
  		
  		
  		// 1. Number of lattice points (this only works for symmetric images).
  		// nx, ny feature points: usually provided by tracking point
  		// correspondences in the phantom during the calibration step.
  		// Here, the distorted and undistorted coordinates from the preprocessing
  		// can be used.
  		
  		// Number of lattice points
  		// TODO: define the number of lattice points
  		// change the value of nx, ny
  		int nx = 0;
  		int ny = 0;
  		
  		// step size
  		// TODO: calculate the stepsize of the lattice points 
  		float fx = 0;
  		float fy = 0;
  		
  		// Fill the distorted and undistorted lattice points with the 
  		// grid coordinates from the preprocessing part.
  		SimpleMatrix Xu2 = new SimpleMatrix(ny,nx);
  		SimpleMatrix Yu2 = new SimpleMatrix(ny,nx);
  		SimpleMatrix Xd2 = new SimpleMatrix(ny,nx);
  		SimpleMatrix Yd2 = new SimpleMatrix(ny,nx);
  		
  		for(int i = 0; i < ny; i++)
  		{
  			for(int j = 0; j < nx; j++)
  			{
  				//TODO: sample the distorted and undistorted grid points at the lattice points
  				// TODO
  				// TODO
  				// TODO
  				// TODO
  			}
  		}
  		
  		// Compute the distorted points: be aware of the fact, that the artificial deformation takes
  		// place from the distorted to undistorted! 
  		// In the Preprocessing: X + distortion = Xd
  		// Now: We correct the distorted image to get the undistorted one!
  		// Thus we have the flip the influence of the distortion field.
  		
  		
  		
  		// Compute the distorted points:
  		// XD2 = XU2 + (XU2 - XD2)
  		// TODO:
  		// TODO:
  		// TODO:
  		// TODO:
  		
  		
  		// 2. Polynom of degree d
  		// Polynom of degree d -> (d-1): extrema
  		// d=0: constant (horizontal line with y-intercept a_0 -> f(x)=a_0)
  		// d=1: oblique line with y-intercept a_0 & slope a_1 -> f(x)=a_0 + a_1 x
  		// d=2: parabola
  		// d>=2: continuous non-linear curve 
  		// E.g. d=5: 4 extrema
  		// d = 10 -> NumKoeff: 66 -> but only 64 lattice points are known
  		int degree = 5; //Polynomial's degree: 2,...,10
  		
  		// Number of Coefficients
  		// TODO:
  		int numCoeff = 0;
  		
  		// Number of Correspondences
  		// TODO:
  		int numCorresp = 0;
  		
  		// Print out of the used parameters
  		System.out.println("Polynom of degree: " + degree);
  		System.out.println("Number of Coefficients: " + numCoeff);
  		System.out.println("Number of Correspondences: " + numCorresp);
  		
  		// 3.Create the matrix A
  		SimpleMatrix A = new SimpleMatrix(numCorresp, numCoeff);
  		A.zeros();
  		
  		// Realign the grid matrix into a vector
  		// Easier access in the next step
  		SimpleVector Xu2_vec = new SimpleVector(numCorresp);
  		SimpleVector Yu2_vec = new SimpleVector(numCorresp);
  		SimpleVector Xd2_vec = new SimpleVector(numCorresp);
  		SimpleVector Yd2_vec = new SimpleVector(numCorresp);
  		
  		for(int i = 0; i < ny; i++)
  		{
  			for(int j = 0; j < nx; j++)
  			{
  				Xu2_vec.setElementValue(i * ny + j, Xu2.getElement(j, i));
  				Yu2_vec.setElementValue(i * ny + j, Yu2.getElement(j, i));
  				Xd2_vec.setElementValue(i * ny + j, Xd2.getElement(j, i));
  				Yd2_vec.setElementValue(i * ny + j, Yd2.getElement(j, i));
  			}
  		}
  		
  		// Compute matrix A
  		for(int r = 0; r < numCorresp; r++)
  		{
  			int cc = 0;
  			for(int i = 0; i <= degree; i++)
  			{
  				for(int j = 0; j <= (degree-i); j++)
  				{
  					// TODO:
  					
  				}
  			}
  		}
  		
  		// Compute the pseudo-inverse of A with the help of the SVD (class: DecompositionSVD)
  		// TODO
  		// TODO
  		
  		
  		// Compute the distortion coefficients
  		// TODO
  		// TODO
  		
  		
  		// 4. Compute the distorted grid points (xDist, yDist) which are used to sample the
  		// distorted image to get the undistorted image
  		// (x,y) is the position in the undistorted image and (XDist,YDist) the
  		// position in the distorted (observed) X-ray image. 
  		
  		Grid2D xDist = new Grid2D(imSize, imSize);  
  		Grid2D yDist = new Grid2D(imSize, imSize); 
  		
  		for(int x = 0; x < imSize; x++)
  		{
  			for(int y = 0; y < imSize; y++)
  			{
  				int cc = 0;
  				for(int k = 0; k <= degree; k++)
  				{
  					for(int l = 0; l <= degree - k; l++)
  					{
  						// TODO
  						// TODO
  						// TODO
  					}
  				}
  			}
  		}
  		
  		Grid2D undistortedImage = new Grid2D(imSize, imSize);
  		
  		for(int i = 0; i < imSize; i++)
  		{
  			for(int j = 0; j < imSize; j++)
  			{
  				// TODO
  				// TODO
  			}
  		}
  		undistortedImage.show("Undistorted Image");
  		
  		Grid2D differenceImage = (Grid2D) NumericPointwiseOperators.subtractedBy(quadraticImage, undistortedImage);
  		differenceImage.show("diffImage");
  	}
  	
  	
  	
  	public static void main(String[] args) {
  		ImageJ ij = new ImageJ();
  		ImageUndistortion iu = new ImageUndistortion();
  		iu.doImageUndstortion();
  	}
  }
  