package edu.stanford.rsl.tutorial.truncation;

import ij.ImageJ;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
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.UniformCircleGrid2D;


/**
 * Example Class to show the linearity of the truncation artifacts.
 * If we split the projection image into the truncated and the non-truncated part,
 * we can see that the filtered backprojection of both combined in image space
 * gives a correct reconstruction. Furthermore, the forward projection of the artifact
 * image gives nearly 0 in the area where the truncation took place. Note that this
 * is only perfect in continuous domain. With discrete sampling a slight residual artifact
 * remains.
 * 
 * @author akmaier
 *
 */
public class WaterCylinderFBPExample {


	public static void main (String [] args){
		new ImageJ();
		int numS = 200;
		double deltaS = 200.0 / 201.0;
		int numTheta = 360;
		int truncationSize = 50;
		UniformCircleGrid2D cylinder = new UniformCircleGrid2D(numS, numS);
		ParallelProjector2D projector2d = new ParallelProjector2D(Math.PI, Math.PI/(numTheta-1), numS, deltaS);
		Grid2D sinogram = projector2d.projectRayDrivenCL(cylinder);
		sinogram.show("Cylinder Sinogram");
		Grid2D truncatedSinogram = new Grid2D(sinogram);
		for (int i=0; i < numTheta; i++) {
			for (int j = 0; j < truncationSize; j++){
				truncatedSinogram.setAtIndex(i, j, 0);
			}
			for (int j = numS - truncationSize+1; j < numS; j++){
				truncatedSinogram.setAtIndex(i, j, 0);
			}
		}
		truncatedSinogram.show("Truncated Image");
		Grid2D artifactSinogram = new Grid2D(sinogram);
		NumericPointwiseOperators.multiplyBy(sinogram, -1);
		artifactSinogram = (Grid2D)NumericPointwiseOperators.addedBy(truncatedSinogram, sinogram);
		NumericPointwiseOperators.multiplyBy(sinogram, -1);
		NumericPointwiseOperators.multiplyBy(artifactSinogram, -1);
		artifactSinogram.show("Artifact Image");
		RamLakKernel ramLakFilter = new RamLakKernel((int) (numS / deltaS), deltaS);
		for (int theta = 0; theta < sinogram.getSize()[0]; ++theta) {
			ramLakFilter.applyToGrid(artifactSinogram.getSubGrid(theta));
			ramLakFilter.applyToGrid(truncatedSinogram.getSubGrid(theta));
		}
		ParallelBackprojector2D backprojector2d = new ParallelBackprojector2D(numS, numS, 1, 1);
		Grid2D truncatedImage = backprojector2d.backprojectPixelDriven(truncatedSinogram);
		truncatedImage.show("Truncated Image");
		Grid2D artifactImage = backprojector2d.backprojectPixelDriven(artifactSinogram);
		artifactImage.show("Artifact Image");
		cylinder = (UniformCircleGrid2D)NumericPointwiseOperators.addedBy(artifactImage, truncatedImage);
		cylinder.show("Combination of both");
		Grid2D FPofArtifactImage = projector2d.projectRayDrivenCL(artifactImage);
		FPofArtifactImage.show("FP of Artifact Image");
		Grid2D FPofTruncImage = projector2d.projectRayDrivenCL(truncatedImage);
		FPofTruncImage.show("FP of Truncation Image");
		Grid2D combinedSinogram = new Grid2D(sinogram);
		combinedSinogram = (Grid2D)NumericPointwiseOperators.addedBy(FPofArtifactImage, FPofTruncImage);
		combinedSinogram.show("Combination of FP");
	}
}
/*
 * Copyright (C) 2010-2014 Andreas Maier
 * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/