content: package edu.stanford.rsl.tutorial.scalespace;
  
  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.NumericGrid;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.tutorial.fan.CosineFilter;
  import edu.stanford.rsl.tutorial.fan.FanBeamBackprojector2D;
  import edu.stanford.rsl.tutorial.fan.FanBeamProjector2D;
  import edu.stanford.rsl.tutorial.filters.RamLakKernel;
  import edu.stanford.rsl.tutorial.phantoms.SheppLogan;
  import edu.stanford.rsl.tutorial.phantoms.UniformCircleGrid2D;
  import edu.stanford.rsl.tutorial.atract.LaplaceKernel2D;
  import edu.stanford.rsl.tutorial.atract.LaplaceKernel1D;
  import ij.ImageJ;
  
  /**
   * Class to convolve projection images of a phantom with either a Gaussian, Laplacian or a Laplacian of Gaussian.
   * Different values for sigma are used to show the impact of their size. Moreover, it is possible to compare the
   * Laplacian of Gaussian with the successive application of filtering an image with the Gaussian and then with the
   * Laplacian.
   * 
   * @author Markus Wolf
   */
  
  public class ScaleSpaceStudies extends Grid2D {
  	
  	
  	public static void main(String args[]) {
  		
  		new ImageJ();
  		
  		// phantom respectively image parameters
  		int phantomType = 4; // 0 = custom phantom
  							 // 1 = SheppLogan
  							 // 2 = UniformCircleGrid2D
  							 // 3 = Ring
  							 // 4 = Dartboard
  		
  		int phantomWidth  = 256,
  			phantomHeight = 256;
  		
  		int imgSizeX = phantomWidth,
  			imgSizeY = phantomHeight;
  		
  		// parameters for Gaussian/Laplacian of Gaussian
  		int kernelsize = 30;
  		
  		int method = 2; // 0 = Gaussian
  						// 1 = Laplace
  						// 2 = Laplacian of Gaussian (LoG)
  						// 3 = Comparison between LoG referred to method 2 and applying Laplacian on a Gaussian
  		
  		double[] sigmaValue = {0.1, 1.0, 3.0, 5.0};
  		
  		// fan beam bp parameters
  		double gammaM      = 11.768288932020647*Math.PI/180,
  			   maxT        = (int)Math.round(Math.sqrt((phantomWidth*phantomWidth) + (phantomHeight*phantomHeight))),
  			   deltaT      = 1.0, 
  			   focalLength = (maxT/2.0-0.5)*deltaT/Math.tan(gammaM),
  		   	   maxBeta     = 360*Math.PI/180,
  		       deltaBeta   = maxBeta / 180;
  		
  		// create phantom		
  		Grid2D phantom;
  		
  		switch(phantomType) {
  		case 0:
  			phantom = new ScaleSpaceStudies(phantomWidth, phantomHeight);
  			break;
  		case 1:
  			phantom = new SheppLogan(phantomWidth,false);
  			break;
  		case 2:
  			phantom = new UniformCircleGrid2D(phantomWidth, phantomHeight);
  			break;
  		case 3:
  			Grid2D outerCircle = new UniformCircleGrid2D(phantomWidth, phantomHeight, 0.4);
  			Grid2D innerCircle = new UniformCircleGrid2D(phantomWidth, phantomHeight, 0.37);
  			phantom = new Grid2D(outerCircle);
  			NumericPointwiseOperators.subtractBy(phantom, innerCircle);
  			break;
  		case 4:
  			phantom = phantomDartboard(phantomWidth, phantomHeight, 0.4);
  			break;
  		default:
  			phantom = new SheppLogan(phantomWidth,false);
  		}
  
  		phantom.show("phantom");
  
  		// create projections (sinogram) out of the phantom
  		FanBeamProjector2D fanBeamProjector = new FanBeamProjector2D(focalLength, maxBeta, deltaBeta, maxT, deltaT);
  				
  		Grid2D projectionP = new Grid2D(phantom);
  		Grid2D fanBeamSinoRay = fanBeamProjector.projectRayDrivenCL(projectionP);	
  		fanBeamSinoRay.show("sinogram");
  		
  		
  		// Laplacian image (method 1)
  		if (method == 1) {
  			// create sinogram
  			Grid2D sinoLaplace = new Grid2D(fanBeamSinoRay);
  		
  			// apply Laplacian
  			LaplaceKernel2D laplacianKernelSino = new LaplaceKernel2D();
  			laplacianKernelSino.applyToGrid(sinoLaplace);
  			sinoLaplace.show("laplacian sinogram");
  		
  			// correct the sinogram	
  			RamLakKernel ramLakLaplace = new RamLakKernel((int) (maxT / deltaT), deltaT);
  			CosineFilter cKernLaplace = new CosineFilter(focalLength, maxT, deltaT);
  							
  			// apply filtering
  			for (int theta = 0; theta < sinoLaplace.getSize()[1]; ++theta) {
  				cKernLaplace.applyToGrid(sinoLaplace.getSubGrid(theta));
  			}
  						
  			for (int theta = 0; theta < sinoLaplace.getSize()[1]; ++theta) {
  				ramLakLaplace.applyToGrid(sinoLaplace.getSubGrid(theta));
  			}
  				
  			// do the backprojection
  			FanBeamBackprojector2D fbpLaplace = new FanBeamBackprojector2D(focalLength, deltaT, deltaBeta, imgSizeX, imgSizeY);
  			Grid2D recoWorkingGridLaplace;
  			recoWorkingGridLaplace = fbpLaplace.backprojectPixelDrivenCL(sinoLaplace);
  			recoWorkingGridLaplace.show("laplace reconstruction");
  		}
  		
  
  		// calculate Gaussian (method 0) or Laplacian of Gaussian (method 2) on projection images for different values of sigma
  		if (method == 0 || method == 2) {
  			
  			for (int i = 0; i < sigmaValue.length; i++) {
  				
  				Grid2D workingGrid;
  				
  				switch (method){
  				case 0:
  					workingGrid = gaussianSinogram(fanBeamSinoRay, sigmaValue[i]);
  					workingGrid.show("gaussian (sigma=" + sigmaValue[i] + ") sinogram");
  					break;
  				case 2:
  					workingGrid = laplacianOfGaussianSinogram(fanBeamSinoRay, kernelsize, sigmaValue[i]);
  					workingGrid.show("laplacian of gaussian (sigma=" + sigmaValue[i] + ") sinogram");
  					break;
  				default:
  					workingGrid = laplacianOfGaussianSinogram(fanBeamSinoRay, kernelsize, sigmaValue[i]);
  					workingGrid.show("laplacian of gaussian (sigma=" + sigmaValue[i] + ") sinogram");
  				}
  				
  				// correct the sinograms	
  				RamLakKernel ramLak = new RamLakKernel((int) (maxT / deltaT), deltaT);
  				CosineFilter cKern = new CosineFilter(focalLength, maxT, deltaT);
  				
  				// apply filtering
  				for (int theta = 0; theta < workingGrid.getSize()[1]; ++theta) {
  					cKern.applyToGrid(workingGrid.getSubGrid(theta));
  				}
  			
  				for (int theta = 0; theta < workingGrid.getSize()[1]; ++theta) {
  					ramLak.applyToGrid(workingGrid.getSubGrid(theta));
  				}
  							
  				// do the backprojection
  				FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength, deltaT, deltaBeta, imgSizeX, imgSizeY);
  				
  				Grid2D recoWorkingGrid;
  				
  				switch (method) {
  				case 0:
  					recoWorkingGrid = fbp.backprojectPixelDrivenCL(workingGrid);
  					recoWorkingGrid.show("gaussian (sigma=" + sigmaValue[i] + ") reconstruction");
  					break;
  				case 2:
  					recoWorkingGrid = fbp.backprojectPixelDrivenCL(workingGrid);
  					recoWorkingGrid.show("laplacian of gaussian (sigma=" + sigmaValue[i] + ") reconstruction");
  					break;
  				default:
  					recoWorkingGrid = fbp.backprojectPixelDrivenCL(workingGrid);
  					recoWorkingGrid.show("laplacian of gaussian (sigma=" + sigmaValue[i] + ") reconstruction");
  				}
  							
  				// difference between phantom and reconstruction
  				// HINT: Only makes sense for Gaussian.
  				if (method == 0) {
  					NumericGrid recoDiff = NumericPointwiseOperators.subtractedBy(phantom, recoWorkingGrid);
  					recoDiff.show("gaussian (sigma=" + sigmaValue[i] + ") RecoDiff");
  				}
  			}
  		}
  				
  		// difference between LoG and Laplacian on Gaussian (method 3)
  		if (method == 3) {
  				
  			for (int i = 0; i < sigmaValue.length; i++) {
  					
  				Grid2D workingGridLoG;
  				Grid2D workingGridLaplace;
  				
  				// sinogram LoG
  				workingGridLoG = laplacianOfGaussianSinogram(fanBeamSinoRay, kernelsize, sigmaValue[i]);
  				workingGridLoG.show("laplacian of gaussian (sigma=" + sigmaValue[i] + ") sinogram");
  				
  				// applying Laplace on gaussian sinogram
  				workingGridLaplace = gaussianSinogram(fanBeamSinoRay, sigmaValue[i]);
  //				LaplaceKernel2D laplacianKernel = new LaplaceKernel2D();
  				LaplaceKernel1D laplacianKernel = new LaplaceKernel1D(workingGridLaplace.getSize()[0]);
  				laplacianKernel.applyToGrid(workingGridLaplace);
  				workingGridLaplace.show("laplacian on gaussian (sigma=" + sigmaValue[i] + ") sinogram");
  				
  				// correct the sinograms	
  				RamLakKernel ramLak = new RamLakKernel((int) (maxT / deltaT), deltaT);
  				CosineFilter cKern = new CosineFilter(focalLength, maxT, deltaT);
  				
  				// apply filtering
  				for (int theta = 0; theta < workingGridLoG.getSize()[1]; ++theta) {
  					cKern.applyToGrid(workingGridLoG.getSubGrid(theta));
  					cKern.applyToGrid(workingGridLaplace.getSubGrid(theta));
  				}
  			
  				for (int theta = 0; theta < workingGridLoG.getSize()[1]; ++theta) {
  					ramLak.applyToGrid(workingGridLoG.getSubGrid(theta));
  					ramLak.applyToGrid(workingGridLaplace.getSubGrid(theta));
  				}
  				
  				// do the backprojection
  				FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength, deltaT, deltaBeta, imgSizeX, imgSizeY);
  				
  				Grid2D recoWorkingGridLoG;
  				Grid2D recoWorkingGridLaplace;
  				
  				recoWorkingGridLoG = fbp.backprojectPixelDrivenCL(workingGridLoG);
  				recoWorkingGridLoG.show("laplacian of gaussian (sigma=" + sigmaValue[i] + ") reconstruction");
  				
  				recoWorkingGridLaplace = fbp.backprojectPixelDrivenCL(workingGridLaplace);
  				recoWorkingGridLaplace.show("laplacian on gaussian (sigma=" + sigmaValue[i] + ") reconstruction");
  				
  				// Difference between LoG and applying Laplacian on a Gaussian
  				NumericGrid diff = NumericPointwiseOperators.subtractedBy(recoWorkingGridLoG, recoWorkingGridLaplace);
  				diff.show("LoG - laplacian on gaussian (sigma=" + sigmaValue[i] + ")");	
  			}
  					
  		}	
  	}
  	
  	/**
  	 * create own phantom consisting of
  	 * <ul><li>a big circle with the intensity 1.0,
  	 * <li>a small circle with the intensity 1.5,
  	 * <li>an ellipse with the intensity 0.9,
  	 * <li>a rectangle with the intensity 1.5
  	 */
  	public ScaleSpaceStudies(int width, int height){
  		super(width, height);
  		
  		// center of grid
  		final double ctr_w = width/2;
  		final double ctr_h = height/2;	
  		
  		// put some concrete shapes in the phantom
  		for (int i = 0; i < width; i++){
  			for (int j = 0; j < height; j++){
  				// big circle
  				if(((i-ctr_w-10)*(i-ctr_w-10)+(j-ctr_h-10)*(j-ctr_h-10)) < Math.pow(19, 2)){
  				putPixelValue(i,j,1.0);
  				}
  				// small circle
  				if(((i-ctr_w+10)*(i-ctr_w+10)+(j-ctr_h+10)*(j-ctr_h+10)) < Math.pow(1.5, 2)){
  					putPixelValue(i,j,1.5);
  				}
  				// left ellipse
  				if((((i-ctr_w+35)*(i-ctr_w+35))/26+((j-ctr_h+20)*(j-ctr_h+20)))/32 < 1){
  					putPixelValue(i,j,0.9);
  				}
  			}
  		}
  		
  		// rectangle
  		for (int i = 10; i < ctr_w; i++){
  			for (int j = (int)ctr_h; j < height-15; j++){
  				if (getPixelValue(i,j) == 0){
  					putPixelValue(i,j,1.5);
  				} else {
  					putPixelValue(i,j, getPixelValue(i,j)+1.5);
  				}		
  			}
  		}
  	}
  	
  	/**
  	 * create phantom that looks like a dartboard with 9 rings
  	 */
  	public static Grid2D phantomDartboard(int width, int height, double r){
  			
  		Grid2D dartboard = new Grid2D(width, height);
  		
  		float[] val = {0.05f, 0.06f, 0.08f, 0.09f, 0.1f, 0.13f, 0.14f, 0.15f, 0.2f};
  		double[] radius = {r * width, (r-0.03)*width, (r-0.06)*width, (r-0.1)*width, (r-0.2)*width, (r-0.25)*width, (r-0.30)*width, (r-0.35)*width, (r-0.38)*width};
  		int xcenter = width / 2;
  		int ycenter = height / 2;
  
  		for (int k = 0; k < val.length; k++) {			
  			for (int i = 0; i < width; i++) {
  				for (int j = 0; j < height; j++) {
  					if (Math.pow(i - xcenter, 2) + Math.pow(j - ycenter, 2) <= (radius[k] * radius[k])) {
  						dartboard.addAtIndex(i, j, val[k]);
  					}
  				}
  			}
  		}
  		return dartboard; 
  	}
  	
  	/**
  	 * convolution of sinogram and Gaussian -> multiplication in Fourier Space
  	 */
  	public static Grid2D gaussianSinogram(Grid2D sinogram, double sigma) {
  		Grid2D gausSinogram = new Grid2D(sinogram);
  		
  		for(int i = 0; i < sinogram.getHeight(); i++){
  			Grid1DComplex sinoRows = new Grid1DComplex(sinogram.getSubGrid(i));
  			
  			// Fourier Transformation
  			sinoRows.transformForward();
  			
  			for(int j = 0; j < sinoRows.getSize()[0]/2; j++){
  				
  				double omega = j*(2*Math.PI/sinoRows.getSize()[0]);
  				
  				float gausValueReal = sinoRows.getRealAtIndex(j)*(float)Math.exp(-0.5*Math.pow(sigma, 2)*Math.pow(omega, 2));
  				float gausValueImag = sinoRows.getImagAtIndex(j)*(float)Math.exp(-0.5*Math.pow(sigma, 2)*Math.pow(omega, 2));
  				
  				sinoRows.setRealAtIndex(j, gausValueReal);
  				sinoRows.setImagAtIndex(j, gausValueImag);
  			}
  			
  			for(int j = sinoRows.getSize()[0]/2; j < sinoRows.getSize()[0]; j++){
  				
  				double omega = (j - sinoRows.getSize()[0])*(2*Math.PI/sinoRows.getSize()[0]);
  				
  				float gausValueReal = sinoRows.getRealAtIndex(j)*(float)Math.exp(-0.5*Math.pow(sigma, 2)*Math.pow(omega, 2));
  				float gausValueImag = sinoRows.getImagAtIndex(j)*(float)Math.exp(-0.5*Math.pow(sigma, 2)*Math.pow(omega, 2));
  				
  				sinoRows.setRealAtIndex(j, gausValueReal);
  				sinoRows.setImagAtIndex(j, gausValueImag);
  			}
  			
  			// inverse Fourier Tranformation
  			sinoRows.transformInverse();
  			
  			// copy again into sinogram
  			for(int k = 0; k < sinogram.getWidth(); k++){
  				float val = (float)Math.sqrt((float)Math.pow(sinoRows.getRealAtIndex(k), 2) + (float)Math.pow(sinoRows.getImagAtIndex(k), 2));
  				gausSinogram.setAtIndex(k, i, val);
  			}
  		}
  		
  		return gausSinogram;
  	}
  
  	/**
  	 * convolution of sinogram and Laplacian of Gaussian
  	 */
  	public static Grid2D laplacianOfGaussianSinogram(Grid2D sinogram, int kernelsize, double sigma) {
  		
  		Grid2D sinogramLoG = new Grid2D(sinogram);
  		
  		for(int i = 0; i < sinogram.getHeight(); i++) {
  			Grid1D sinoRow = new Grid1D(sinogram.getSubGrid(i));
  
  			Grid1D result = convolutionLoG(sinoRow, kernelsize, sigma);
  		
  			// copy again into sinogram
  			for(int k = 0; k < sinogram.getWidth(); k++){
  				sinogramLoG.setAtIndex(k, i, result.getAtIndex(k));
  			}
  		}	
  		
  		return sinogramLoG;	
  	}
  
  	/**
  	 * compute convolution with Laplacian of Gaussian
  	 */
  	public static Grid1D convolutionLoG(Grid1D signal, int kernelsize, double sigma) {
  
  		int N = signal.getSize()[0];
  		Grid1D result = new Grid1D(N);
  
  		for (int n = 0; n < N; n++) { // n is number of signal element
  
  			float valueConvolved = 0.0f;			
  			result.setAtIndex(n, 0);
  
  			// borders for kernel window
  			int kmin, kmax;
  			int kernelrad = kernelsize/2;
  			kmin = (n > kernelrad) ? n - kernelrad : 0;
  			kmax = (n < N - kernelrad) ? n + kernelrad : N - 1;
  				
  			// local sampling
  			float deltaU;			// additional factor (between 0.0 and 0.9) for sampling rate
  			float v;				// sampling rate
  			float interpolSigVal;	// interpolated signal value
  			
  			for (int k = kmin; k < kmax; k++) {
  				for (int u = 0; u < 10; u++) {
  					
  					deltaU = 0.1f * u;
  					v = k + deltaU;
  					
  					interpolSigVal = (1-deltaU)*signal.getAtIndex(k) + deltaU*signal.getAtIndex(k+1);	// interpolation due to local sampling
  					valueConvolved += interpolSigVal * laplacianOfGaussian(n-v, sigma);					// actual convolution
  				}
  			}
  			
  			result.setAtIndex(n, valueConvolved);
  		}
  		
  		return result;
  	}
  	
  	/**
  	 * Laplacian of Gaussian -> second derivative of Gaussian distribution
  	 */
  	public static double laplacianOfGaussian(double index, double sigma) {
  		
  		double value = (-1/(Math.pow(sigma, 3) * Math.sqrt(2*Math.PI))) * (1 - (Math.pow(index, 2) / Math.pow(sigma, 2))) * Math.exp(-(Math.pow(index, 2) / (2*Math.pow(sigma, 2))));
  	
  		return value;
  	}
  
  }
  
  /*
   * Copyright (C) 2010-2015 Markus Wolf
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  