content: package edu.stanford.rsl.tutorial.dmip;
  
  import java.util.ArrayList;
  
  import edu.stanford.rsl.conrad.data.numeric.Grid1D;
  import edu.stanford.rsl.conrad.data.numeric.Grid1DComplex;
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.Box;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
  import edu.stanford.rsl.conrad.geometry.transforms.Transform;
  import edu.stanford.rsl.conrad.geometry.transforms.Translation;
  import edu.stanford.rsl.conrad.numerics.SimpleOperators;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  import edu.stanford.rsl.tutorial.phantoms.SheppLogan;
  import ij.ImageJ;;
  
  /**
   * Exercise 5 of Diagnostic Medical Image Processing (DMIP)
   * @author Bastian Bier
   *
   */
  
  public class DMIP_ParallelBeam {
  	
  	
  	public enum RampFilterType {NONE, RAMLAK, SHEPPLOGAN};
  	
  	
  	/**
  	 * Forward projection of the phantom onto the detector
  	 * Rule of thumb: Always sample in the domain where you expect the output!
  	 * Thus, we sample at the detector pixel positions and sum up the informations along one ray
  	 * 
  	 * @param grid the image
  	 * @param maxTheta the angular range in radians
  	 * @param deltaTheta the angular step size in radians
  	 * @param maxS the detector size in [mm]
  	 * @param deltaS the detector element size in [mm]
  	 */
  	
  	public Grid2D projectRayDriven(Grid2D grid, double maxTheta, double deltaTheta, double maxS, double deltaS) {
  		
  		int maxSIndex = (int) (maxS / deltaS + 1);
  		int maxThetaIndex = (int) (maxTheta / deltaTheta + 1);
  		
  		final double samplingRate = 3.d; // # of samples per pixel
  		Grid2D sino = new Grid2D(new float[maxThetaIndex*maxSIndex], maxSIndex, maxThetaIndex);
  		sino.setSpacing(deltaS, deltaTheta);
  
  		// set up image bounding box in WC
  		Translation trans = new Translation(
  				-(grid.getSize()[0] * grid.getSpacing()[0])/2, -(grid.getSize()[1] * grid.getSpacing()[1])/2, -1);
  		Transform inverse = trans.inverse();
  
  		Box b = new Box((grid.getSize()[0] * grid.getSpacing()[0]), (grid.getSize()[1] * grid.getSpacing()[1]), 2);
  		b.applyTransform(trans);
  
  		for(int e=0; e<maxThetaIndex; ++e){
  			// compute theta [rad] and angular functions.
  			double theta = deltaTheta * e;
  			double cosTheta = Math.cos(theta);
  			double sinTheta = Math.sin(theta);
  
  			for (int i = 0; i < maxSIndex; ++i) {
  				// compute s, the distance from the detector edge in WC [mm]
  				double s = deltaS * i - maxS / 2;
  				// compute two points on the line through s and theta
  				// We use PointND for Points in 3D space and SimpleVector for directions.
  				PointND p1 = new PointND(s * cosTheta, s * sinTheta, .0d);
  				PointND p2 = new PointND(-sinTheta + (s * cosTheta),
  						(s * sinTheta) + cosTheta, .0d);
  				// set up line equation
  				StraightLine line = new StraightLine(p1, p2);
  				// compute intersections between bounding box and intersection line.
  				ArrayList<PointND> points = b.intersect(line);
  
  				// only if we have intersections
  				if (2 != points.size()){
  					if(points.size() == 0) {
  						line.getDirection().multiplyBy(-1.d);
  						points = b.intersect(line);
  					}
  					if(points.size() == 0)
  						continue;
  				}
  
  				PointND start = points.get(0); // [mm]
  				PointND end = points.get(1);   // [mm]
  
  				// get the normalized increment
  				SimpleVector increment = new SimpleVector(
  						end.getAbstractVector());
  				increment.subtract(start.getAbstractVector());
  				double distance = increment.normL2();
  				increment.divideBy(distance * samplingRate);
  
  				double sum = .0;
  				start = inverse.transform(start);
  
  				// compute the integral along the line.
  				for (double t = 0.0; t < distance * samplingRate; ++t) {
  					PointND current = new PointND(start);
  					current.getAbstractVector().add(increment.multipliedBy(t));
  
  					double x = current.get(0) / grid.getSpacing()[0],
  							y = current.get(1) / grid.getSpacing()[1];
  
  					if (grid.getSize()[0] <= x + 1
  							|| grid.getSize()[1] <= y + 1
  							|| x < 0 || y < 0)
  						continue;
  
  					sum += InterpolationOperators.interpolateLinear(grid, x, y);
  				}
  
  				// normalize by the number of interpolation points
  				sum /= samplingRate;
  				// write integral value into the sinogram.
  				sino.setAtIndex(i, e, (float)sum);
  			}
  		}
  		return sino;
  	}
  	
  	/**
  	 * Sampling of projections is defined in the constructor.
  	 * Backprojection of the projections/sinogram
  	 * Rule of thumb: Always sample in the domain where you expect the output!
  	 * Here, we want to reconstruct the volume, thus we sample in the reconstructed grid!
  	 * The projections are created pixel driven!
  	 * 
  	 * @param sino  the sinogram
  	 * @param imageSizeX
  	 * @param imageSizeY
  	 * @param pxSzXMM
  	 * @param pxSzYMM	 
  	 */
  	
  	public Grid2D backprojectPixelDriven(Grid2D sino, int imageSizeX, int imageSizeY, float pxSzXMM, float pxSzYMM) {
  		
  		int maxThetaIndex = sino.getSize()[1];
  		double deltaTheta = sino.getSpacing()[1];
  		int maxSIndex = sino.getSize()[0];
  		double deltaS = sino.getSpacing()[0];
  		double maxS = (maxSIndex-1) * deltaS;
  		
  		Grid2D grid = new Grid2D(imageSizeX, imageSizeY);
  		grid.setSpacing(pxSzXMM, pxSzYMM);
  		grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0])/2, -(grid.getSize()[1]*grid.getSpacing()[1])/2);
  
  		// loop over the projection angles
  		for (int i = 0; i < maxThetaIndex; i++) {
  			// compute actual value for theta
  			double theta = deltaTheta * i;
  			// precompute sine and cosines for faster computation
  			double cosTheta = Math.cos(theta);
  			double sinTheta = Math.sin(theta);
  			// get detector direction vector
  			SimpleVector dirDetector = new SimpleVector(sinTheta,cosTheta);
  			// loops over the image grid
  			for (int x = 0; x < grid.getSize()[0]; x++) {
  				for (int y = 0; y < grid.getSize()[1]; y++) {
  					// compute world coordinate of current pixel
  					double[] w = grid.indexToPhysical(x, y);
  					// wrap into vector
  					SimpleVector pixel = new SimpleVector(w[0], w[1]);
  					//  project pixel onto detector
  					double s = SimpleOperators.multiplyInnerProd(pixel, dirDetector);
  					// compute detector element index from world coordinates
  					s += maxS/2; // [mm]
  					s /= deltaS; // [GU]
  					// get detector grid
  					Grid1D subgrid = sino.getSubGrid(i);
  					// check detector bounds, continue if out of array
  					if (subgrid.getSize()[0] <= s + 1
  							||  s < 0)
  						continue;
  					// get interpolated value
  					float val = InterpolationOperators.interpolateLinear(subgrid, s);
  					// sum value to sinogram
  					grid.addAtIndex(x, y, val);
  				}
  
  			}
  		}
  		// apply correct scaling
  		NumericPointwiseOperators.divideBy(grid, (float) (maxThetaIndex / Math.PI));
  		return grid;
  	}
  	
  	
  	/**
  	 * Filtering the sinogram with a high pass filter
  	 * 
  	 * The ramp filters are defined in the spatial domain but 
  	 * they are applied in the frequency domain.
  	 * Remember: a convolution in the spatial domain corresponds to a 
  	 * multiplication in the frequency domain.
  	 * 
  	 * Both, the sinogram and the ramp filter are transformed into
  	 * the frequency domain and multiplied there.
  	 * 
  	 * @param sinogram  a line of the sinogram
  	 *  
  	 */
  	public Grid1D rampFiltering(Grid1D sinogram, RampFilterType filter){
  		
  		double deltaS = 1;
  		
  		// Initialize the ramp filter
  		// Define the filter in the spatial domain on the full padded size!
  		Grid1DComplex ramp = new Grid1DComplex(sinogram.getSize()[0]);
  		
  		int paddedSize = ramp.getSize()[0];
  		
  		if(filter == RampFilterType.RAMLAK)
  		{
  			// TODO: implement the ram-lak filter in the spatial domain 
  			
  		}
  		else if(filter == RampFilterType.SHEPPLOGAN)
  		{
  			// TODO: implement the Shepp-Logan filter in the spatial domain
  			
  		}
  		else
  		{
  			// if no filtering is used
  			return sinogram;
  		}
  		
  		// TODO: Transform ramp filter into frequency domain
  		
  		
  		Grid1DComplex sinogramF = new Grid1DComplex(sinogram,true);
  		// TODO: Transform the input sinogram signal into the frequency domain
  		
  		
  		// TODO: Multiply the ramp filter with the transformed sinogram
  		
  		
  		// TODO: Backtransformation
  
  		
  		// Crop the image to its initial size
  		Grid1D ret = new Grid1D(sinogram);
  		ret = sinogramF.getRealSubGrid(0, sinogram.getSize()[0]);
  		
  		return ret;
  	}
  	
  	public static void main(String[] args)
  	{
  		ImageJ ij = new ImageJ();
  		
  		DMIP_ParallelBeam parallel = new DMIP_ParallelBeam();
  		
  		// 0. Parameters
  		
  		// size of the phantom
  		int phantomSize = 128; 
  		// projection image range 	
  		double angularRange = Math.PI; 	
  		// number of projection images	
  		int projectionNumber = 180;	
  		// angle in between adjacent projections
  		double angularStepSize 	= angularRange / projectionNumber;
  		// detector size in [mm]
  		float detectorSize = 200; 
  		// size of a detector Element [mm]
  		float detectorSpacing = 1.0f;	
  		// filterType: NONE, RAMLAK, SHEPPLOGAN
  		RampFilterType filter = RampFilterType.NONE;				
  		
  		// 1. Create the Shepp Logan Phantom
  		SheppLogan sheppLoganPhantom = new SheppLogan(phantomSize);
  		sheppLoganPhantom.show();
  		
  		// 2. Acquire forward projection images with a parallel projector
  		Grid2D sinogram = parallel.projectRayDriven(sheppLoganPhantom, angularRange, angularStepSize, detectorSize, detectorSpacing);
  		sinogram.show("The Sinogram");
  		Grid2D filteredSinogram = new Grid2D(sinogram);
  		
  		// 4. Ramp Filtering
  		for (int theta = 0; theta < sinogram.getSize()[1]; ++theta) 
  		{
  			// Filter each line of the sinogram independently
  			Grid1D tmp = parallel.rampFiltering(sinogram.getSubGrid(theta), filter);
  			
  			for(int i = 0; i < tmp.getSize()[0]; i++)
  			{
  				filteredSinogram.putPixelValue(i, theta, tmp.getAtIndex(i));
  			}
  		}
  		
  		filteredSinogram.show("Filtered Sinogram");
  		
  		// 5. Reconstruct the object with the information in the sinogram	
  		Grid2D reco = parallel.backprojectPixelDriven(filteredSinogram, phantomSize, phantomSize, detectorSpacing, detectorSpacing);
  		reco.show("Reconstruction");
  	}
  	
  }
  