content: package edu.stanford.rsl.tutorial.truncation;
  
  import ij.ImageJ;
  import edu.stanford.rsl.apps.gui.Citeable;
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.tutorial.filters.RamLakKernel;
  import edu.stanford.rsl.tutorial.parallel.ParallelBackprojector2D;
  import edu.stanford.rsl.tutorial.parallel.ParallelProjector2D;
  import edu.stanford.rsl.tutorial.phantoms.Phantom;
  
  /**
   * Test implementation to show that a constant offset inside the ROI in the artifact image is not possible.
   * Enforcing the 0 attenuation inside the FOV and a constant value in the artifact image converges to 0.
   * @author akmaier
   *
   */
  public class NonZeroArtifactImage extends Phantom implements Citeable {
  
  	public String getBibtexCitation() {
  		return "@article{Yu2009, \n"+
  		"  author = {Yu, Hengyong and Yang, Jiansheng and Jiang, Ming and Wang, Ge},\n"+
  		"  doi = {10.1088/0031-9155/54/18/N04.Supplemental},\n"+
  		"  file = {:C$\backslash$:/Users/z002xsyz/Desktop/supplemental.pdf:pdf},\n"+
  		"  institution = {CT Laboratory, Biomedical Imaging Division, VT-WFU School of Biomedical Engineering, Virginia Tech, Blacksburg, VA 24061, USA. hengyong-yu@ieee.org},\n"+
  		"  journal = {Physics in Medicine and Biology},\n"+
  		"  keywords = {artifact image,compressed sensing,computed tomography,ct,interior tomography,minimization,total variation},\n"+
  		"  number = {18},\n"+
  		"  pages = {N425--N432},\n"+
  		"  publisher = {IOP Publishing},\n"+
  		"  title = {{Supplemental analysis on compressed sensing based interior tomography.}},\n"+
  		"  volume = {54},\n"+
  		"  year = {2009}\n"+
  		"  }";
  	}
  
  	public String getMedlineCitation() {
  		return "Yu, H., Yang, J., Jiang, M., & Wang, G. (2009). Supplemental analysis on compressed sensing based interior tomography. Physics in Medicine and Biology, 54(18), N425–N432. doi:10.1088/0031-9155/54/18/N04.Supplemental";
  	}
  
  	
  	public NonZeroArtifactImage(int x, int y) {
  		super(x, y, "NonZeroArtifactImage");
  		boolean showProgress = false;
  
  		// Create circle in the grid.
  		double radius = 0.40 * x;
  		double radius2 = 0.50 * x;
  		int xcenter = x / 2;
  		int ycenter = y / 2;
  		float val = 1.f;
  
  		for (int i = 0; i < x; i++) {
  			for (int j = 0; j < y; j++) {
  				if (Math.pow(i - xcenter, 2) + Math.pow(j - ycenter, 2) <= (radius * radius)) {
  					super.setAtIndex(i, j, val);
  				} else {
  					if (Math.pow(i - xcenter, 2) + Math.pow(j - ycenter, 2) <= (radius2 * radius2)){
  						super.setAtIndex(i, j, -4.f);
  					}
  				}
  			}
  		}
  		int iterations = 1001;
  		int truncationSize = x/4;
  		radius = truncationSize;
  		int numS = x;
  		double deltaS = x / (double)(x+1);
  		int numTheta = 360;
  
  		ParallelProjector2D projector2d = new ParallelProjector2D(Math.PI, Math.PI/(numTheta-1), numS, deltaS);
  		Grid2D sinogram = projector2d.projectRayDrivenCL(this);
  		if (showProgress) sinogram.show("Cylinder Sinogram");
  		Grid2D truncatedSinogram = new Grid2D(sinogram);
  
  		for (int iter = 1; iter < iterations; iter++){
  			if (showProgress) if (iter%50 == 0) truncatedSinogram.show("Artifact Sinogram " + iter);
  			truncatedSinogram = new Grid2D(truncatedSinogram);
  			for (int i=0; i < numTheta; i++) {
  				for (int j = truncationSize; j < numS - truncationSize; j++){
  					truncatedSinogram.setAtIndex(i, j, 0);
  				}
  			}
  			RamLakKernel ramLakFilter = new RamLakKernel((int) (numS / deltaS), deltaS);
  			for (int theta = 0; theta < sinogram.getSize()[0]; ++theta) {
  				ramLakFilter.applyToGrid(truncatedSinogram.getSubGrid(theta));
  			}
  			ParallelBackprojector2D backprojector2d = new ParallelBackprojector2D(numS, numS, 1, 1);
  			Grid2D truncatedImage = backprojector2d.backprojectPixelDriven(truncatedSinogram);
  			if (showProgress) if (iter%50 == 0) truncatedImage.show("Artifact Image" + iter);
  			truncatedImage = new Grid2D(truncatedImage);
  			double sum = 0;
  			long number = 0;
  			for (int i = 0; i < x; i++) {
  				for (int j = 0; j < y; j++) {
  					if (Math.pow(i - xcenter, 2) + Math.pow(j - ycenter, 2) <= (radius * radius)) {
  						sum += truncatedImage.getAtIndex(i, j);
  						number++;
  					}
  					if (Math.pow(i - xcenter, 2) + Math.pow(j - ycenter, 2) > (radius2 * radius2)){
  						truncatedImage.setAtIndex(i, j, 0);
  					}
  				}
  			}
  			double mean = sum / number;
  			for (int i = 0; i < x; i++) {
  				for (int j = 0; j < y; j++) {
  					if (Math.pow(i - xcenter, 2) + Math.pow(j - ycenter, 2) <= (radius * radius)) {
  						truncatedImage.setAtIndex(i, j, (float)mean);
  					} 
  				}
  			}
  			truncatedSinogram = projector2d.projectRayDrivenCL(truncatedImage);
  		}
  	}
  
  	/**
  	 * @param args
  	 */
  	public static void main(String[] args) {
  		new ImageJ();
  		Grid2D phantom = new NonZeroArtifactImage(200, 200);
  		phantom.show("Phantom");
  	}
  
  }
  /*
   * Copyright (C) 2010-2014 Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */