/*
 * Copyright (C) 2014 Madgalena Herbst
 * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
package edu.stanford.rsl.tutorial.noncircularfov;

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


/**
 * Short example of how to use the Correspondences class.
 * 
 * @author Magdalena
 *
 */
public class TrajectoryExample {
	
	public static void main(String[] args) {
		
		// required parameters
		int imgSzXMM = 256, // [mm]
		imgSzYMM = imgSzXMM; // [mm]
		float pxSzXMM = 1.0f, // [mm]
		pxSzYMM = pxSzXMM; // [mm]
		int imgSzXGU = (int) Math.floor(imgSzXMM / pxSzXMM), // [GU]
		imgSzYGU = (int) Math.floor(imgSzYMM / pxSzYMM); // [GU]
		
		double gammaM =2* 11.768288932020647 * Math.PI / 180, maxT = 501, deltaT = 1.0, focalLength = (maxT / 2.0 - 0.5)
				* deltaT / Math.tan(gammaM), maxBeta = 360 * Math.PI / 180, // +gammaM*2,
		deltaBeta = maxBeta / 360;
		
		// create Phantom
		System.out.println("Creating phantom...");
		Phantom phantom = new TwoCirclesGrid(imgSzXMM,imgSzYMM);
		phantom.setSpacing(pxSzXMM, pxSzYMM);
		// origin is given in (negative) world coordinates
		phantom.setOrigin(-(imgSzXGU * phantom.getSpacing()[0]) / 2,
				-(imgSzYGU * phantom.getSpacing()[1]) / 2);
		
		new ImageJ();
		
		phantom.show("Phantom");
		
		// make projections
		System.out.println("projecting...");
		FanBeamProjector2D fanBeamProjector = new FanBeamProjector2D(focalLength, maxBeta, deltaBeta, maxT, deltaT);
		
		Grid2D projectionP = new Grid2D(phantom);

		Grid2D sino = fanBeamProjector.projectRayDriven(projectionP);
		sino.clone().show("Sinogram");
		
		// computing the minimal set
		System.out.println("Computing the minimal set...");
		
		double detSize = 130;	// [mm], pixelsize deltaT
		int detSizeGU = (int) (detSize/deltaT);
		
		Correspondences c = new Correspondences(focalLength, sino);
		
		int[] minimalSet = c.evaluateNumer(sino, detSizeGU);
		
		System.out.println("Minimal required rotation: "  + minimalSet[0]*deltaBeta  + "degress, corresp. start angle: " + minimalSet[1]*deltaBeta);
		
		// displaying the minimal set
		System.out.println("Display minimal set");
		
		int width = sino.getWidth();
		int height = sino.getHeight();
		
		Grid2D mask = new Grid2D(width, height);
		Grid2D data = new Grid2D(width, height);
		
		int start = minimalSet[1];
		int rot = minimalSet[0];
		
		for (int y = start; y < start + rot; y++) {
			for (int x = 0; x < width; x++) {
				//if (sino.getAtIndex(x, y%height) > 0) {
				if((x>150) && (x < 150+detSize)){
					for (int i = 0; i < detSize; i++) {
						if (x + i < width) {
							mask.setAtIndex(x + i, y%height, 1.f);
							data.setAtIndex(x + i, y%height, sino.getAtIndex(x + i, y%height));
						}
					}
					break;
				}
			}
		}
		
		Grid2D res = c.findCorrespondences(mask, sino);
		res.setSpacing(deltaT, deltaBeta);	// required for reconstruction
		Grid2D dou = c.findRedundancies(mask, sino);
		Grid2D diff = new Grid2D(width, height);
		
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				
				// count missing data
				if(sino.getAtIndex(x, y) > 0 && res.getAtIndex(x, y) == 0) {
					diff.setAtIndex(x, y, 1.f);
				}
			}
		}
		
		data.show("Acquired Data");
		res.show("Reconstructed Sinogram");
		dou.show("redundant areas");
		diff.show("Missing Pixels");
		mask.show("Mask");
		
		// resonstruct the completed sinogram
		System.out.println("Reconstruct the completed sinogram...");
		
		RamLakKernel ramLak = new RamLakKernel((int) (maxT / deltaT), deltaT);
		CosineFilter cKern = new CosineFilter(focalLength, maxT, deltaT);
		// Apply filtering
		for (int theta = 0; theta < res.getSize()[1]; ++theta) {
			cKern.applyToGrid(res.getSubGrid(theta));

		}

		for (int theta = 0; theta < res.getSize()[1]; ++theta) {
			ramLak.applyToGrid(res.getSubGrid(theta));
		}
		
		res.show("After Filtering");
		
		// Do the backprojection
		FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength,
				deltaT, deltaBeta, imgSzXMM, imgSzYMM);

		Grid2D reco = fbp.backprojectPixelDriven(res);
		reco.show("Reconstructed result");
		
		//finished
		System.out.println("Everything done!");
	}

}