MotionCompensatedRecon.java

/* Copyright (C) 2015 Bernhard Stimpel
 * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
 */
package edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation;

import java.io.BufferedReader;
import java.io.DataInputStream;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.util.ArrayList;
import java.util.StringTokenizer;

import ij.IJ;
import ij.ImageJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.measure.ResultsTable;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import ij.process.StackStatistics;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.filtering.CosineWeightingTool;
import edu.stanford.rsl.conrad.filtering.ImageFilteringTool;
import edu.stanford.rsl.conrad.filtering.RampFilteringTool;
import edu.stanford.rsl.conrad.filtering.rampfilters.RampFilter;
import edu.stanford.rsl.conrad.filtering.rampfilters.SheppLoganRampFilter;
import edu.stanford.rsl.conrad.filtering.redundancy.ParkerWeightingTool;
import edu.stanford.rsl.conrad.geometry.Projection;
import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
import edu.stanford.rsl.conrad.io.FileProjectionSource;
import edu.stanford.rsl.conrad.io.ImagePlusDataSink;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.pipeline.BufferedProjectionSink;
import edu.stanford.rsl.conrad.pipeline.ParallelImageFilterPipeliner;
import edu.stanford.rsl.conrad.pipeline.ProjectionSource;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.ImageUtil;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.MotionCompensatedRecon;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.ECG.ECGGating;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.mip.OpenCLMaximumIntensityProjection;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.morphology.Morphology;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.morphology.StructuringElement;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.reconWithStreakReduction.ConeBeamBackprojectorStreakReduction;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.reconWithStreakReduction.ConeBeamBackprojectorStreakReductionWithMotionCompensation;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.MiscTools;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.Param;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.Transformation;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.bUnwarpJ_;



/**
 * Class to model Cardiac Vasculature Reconstruction with Residual Motion 
 * Compensation using ECG gating, Backprojection with Streak Reduction
 * Forward MIP and MotionCompensated Reconstruction
 * */
public class MotionCompensatedRecon{
	
	//	TODO Set ROI dimensions if needed
	//	x and y define the top left corner of the ROI.  		
	private int xROI = 0;
	private int yROI = 0;
	private int widthROI = 960;
	private int heightROI = 960;	
	
	public Grid3D transformationVectorsX = null;
	public Grid3D transformationVectorsY = null;
		
	// Projection geometry parameters
	private double sizeU = 960;
	private double sizeV = 960;
	private double spacingU = 0.32;
	private double spacingV = 0.32;
	private double sourceDetectorDist = 1200;
	private double sourceIsocenterDist = 800;
	
	private int gat_ign = 3;	
	
	// ECG Parameters
	private String ecgFile = "C:/workspace/CavarevData/phasedata_card.txt";
	static double REF_HEART_PHASE = 0.90;	
	private double gatingWidth = 4.0;
	private double gatingSlope = 4.0;
	
	Configuration cnfg = null;

	// Path to the ecg and pmat files
	private String projMatricesFile = "C:/workspace/CavarevData/cavarev.matrices.bin";
	static String projectionFile = "C:/workspace/CavarevData/cavarev_noBreathing.tif";
	static String registrationPath = "C:/workspace/CavarevData/Registration/";

	public MotionCompensatedRecon(){	
	}
	
	public static void main(String args[]){
		
		new ImageJ();
		
		//	Computing the initial reconstruction
		MotionCompensatedRecon moCoRecon = new MotionCompensatedRecon();		
		Grid3D initialReco = null;		
		try {
			System.out.println("Starting the initial reconstruction.");
			initialReco = moCoRecon.reconstructConeBeam(false, moCoRecon.getGat_ign());		
			} catch (Exception e) {
			e.printStackTrace();
		}

		Grid3D resultingReco = moCoRecon.IterativeSteps(initialReco);
		resultingReco.show();
	}
			
	
	/**
	 * Steps that can be repeated iteratively during Cardiac Vasculature Reconstruction
	 * @param resultRecon: result of reconstruction (in first iteration result of initial reconstruction in next iteration result of MotionCompReconstruction)
	 * @param projectionsPath: Path where the projections are saved
	 * @param iterationCounter: used to save image with corresponding Number of Iteration
	 * @param referenceHeartPhase: ECG gating parameter
	 * @param numHeartBeats: ECG gating parameter
	 * */
	public Grid3D IterativeSteps(Grid3D resultRecon){		
		
		Grid3D source = null;
		Grid3D target = null;
		configure();
		
		// Preparing the source and target images for the registration process
		System.out.println("Preparing the initial reconstruction.");
		target = prepReconstruction(resultRecon);

		System.out.println("Pre-processing the original projections.");
		source = preProcessProj();
		
		transformationVectorsX = new Grid3D(source.getSize()[0], source.getSize()[1], source.getSize()[2]);
		transformationVectorsY = new Grid3D(target.getSize()[0], target.getSize()[1], target.getSize()[2]);
		
		// Cropping the images to save computation time.
		source = cropImage(source, xROI, yROI, widthROI, heightROI);	
		target = cropImage(target, xROI, yROI, widthROI, heightROI);
		
		bUnwarpJ_ bUnwarpJ = new bUnwarpJ_();

		//	Starting the registration process. The resulting deformationfields are saved in the transformationVectors. 
		performRegistration(target, source, bUnwarpJ);
		
		//	Performing the reconstruction.
		resultRecon = reconstructConeBeam(true, getGat_ign());
		
		return resultRecon;
	}
	
	/**
	 * This method performs a percentile histogram thresholding followed by a Maximum Intensity Forward Projection.
	 * @param initReconstruction: the initial reconstruction
	 * @return the Maximum Intensity Forward Projection as Grid3D
	 */
	public Grid3D prepReconstruction(Grid3D initReconstruction){
		ImagePlus reconIP = ImageUtil.wrapGrid3D(initReconstruction, "");
		
		//TODO Adjusting the thresholding values in this method may be necessary 
		
		//	Adapt thresholding value to current image intensity values via percentile thresholding based on histograms
		double percentile = 0.001;	
		Grid3D reconGrid = percentileHistogramThreshold(reconIP, percentile);
		
		//	Perform the maximum intensity projection
		Grid3D mipProjGrid = mipProjection(reconGrid);
		
		//	Second thresholding operation
		Grid3D fwpCopy = threshold3D(mipProjGrid, 0.011);
		return fwpCopy;
	}
	
	/** This methods performs the pre-processing of the original projections. 
	 * In the first step a top-hat filter is applied to the projections. The necessary parameters
	 * are defined in the CoronaryConfig. In the second step the achieved result is binarized.  
	 * @return the filtered and binarized projections as Grid3D
 	 */	
	public Grid3D preProcessProj(){		
		//TODO	Adjusting the filter size and binarization threshold in this method may be necessary
		
		//Performing the morphological filtering
		ImagePlus preProjIP = morphFilterProjs(projectionFile, 10, "topHat");
		
		//Wrapping and binarizing the morph-filtered projections
		Grid3D preProjGrid = ImageUtil.wrapImagePlus(preProjIP);			
		preProjGrid = binarize3D(preProjGrid, 3.0);		
		
		return preProjGrid;
	}
	
	/**
	 * ConeBeamReconstruction Method using ConeBeamBackprojector from Tutorial Code
	 * including: ConeBeamCosineFilter, Parker Weights and RamLakFilter
	 * */
	public Grid3D reconstructConeBeam(boolean motionCompensation, int gat_ign){
		
		configure();		
		Grid3D sino = null;
		Grid3D resultRecon = null;		

		// Preparing the projections in advance to the reconstruction. 
		try {
			sino = prepProjections();
		} catch (Exception e) {
			e.printStackTrace();
		}		
		
		//	Performing the actual backprojection, either with or without motion compensation. 
		if(motionCompensation == true){
			ConeBeamBackprojectorStreakReductionWithMotionCompensation cbbp = new ConeBeamBackprojectorStreakReductionWithMotionCompensation();
			cbbp.setGat_ign(gat_ign);
			resultRecon = cbbp.backprojectPixelDrivenCL(sino, transformationVectorsX, transformationVectorsY);
		}
		else{
			ConeBeamBackprojectorStreakReduction cbbp = new ConeBeamBackprojectorStreakReduction();
			cbbp.setGat_ign(gat_ign);		
			resultRecon = cbbp.backprojectPixelDrivenCL(sino);
		}
		return resultRecon;
	}
	
	/**
	 * Method performing the Maximum Intensity Projection
	 * @param inputProjection: The pre-processed inputProjection as Grid3D
	 * @return the resulting Maximum Intensity Forward Projection as Grid3D
	 */
	public Grid3D mipProjection(Grid3D inputProjection){
		ImagePlus resultingMIP = null;
		
		//perform the actual Maximum Intensity Projection
		OpenCLMaximumIntensityProjection clForwardProjector = new OpenCLMaximumIntensityProjection();		
		try {
			clForwardProjector.configure();
		} catch (Exception e) {
			e.printStackTrace();
		}				
		clForwardProjector.setTex3D(ImageUtil.wrapGrid3D(inputProjection, "Maximum Intensity Forward Projection"));
		resultingMIP = clForwardProjector.project();
		
		//Wraping and returning the result
		Grid3D resultingMIPGrid = ImageUtil.wrapImagePlus(resultingMIP);
		return resultingMIPGrid;
	}
	
	/**
	 * Preparing the original projections for reconstruction.
	 * @param projections : Original projections (currently not used; projectionstream as source)
	 * @return The prepared projections 
	 * @throws Exception
	 */
	private Grid3D prepProjections() throws Exception{
		
		Trajectory geo  = cnfg.getGeometry();
		
		//	Initializing and configuring the tools. 
		CosineWeightingTool cbTool = new CosineWeightingTool();
		cbTool.configure();		
		ParkerWeightingTool redundancyTool = new ParkerWeightingTool(geo);
		redundancyTool.configure();		
		RampFilteringTool rampFiltTool = new RampFilteringTool();		
		RampFilter ramp = new SheppLoganRampFilter();		
		rampFiltTool.setRamp(ramp);
		rampFiltTool.setConfigured(true);
		
		ImageFilteringTool[] filters = new ImageFilteringTool[]{cbTool,redundancyTool,rampFiltTool};
		
		//	Retrieving the projection source as projection stream.
		ProjectionSource pSource = FileProjectionSource.openProjectionStream(projectionFile);
		BufferedProjectionSink sink = new ImagePlusDataSink();
		
		//	Creating the filter pipeline. 
		ParallelImageFilterPipeliner filteringPipeline = 
				new ParallelImageFilterPipeliner(pSource, filters, sink);
				
		//	Executing the previously created pipeline.
		Thread thread = new Thread(new Runnable(){
			public void run(){
				try {
					filteringPipeline.project();
				} catch (Exception e) {
					e.printStackTrace();
				}
			} 
		});
		thread.start();
		
		//	Applying the ECG-gating to the projection source 
		Grid3D preppedProj = sink.getResult();
		double[] ecg = readEcg(ecgFile);
		preppedProj = applyGating(preppedProj, ecg);
		
		return preppedProj;
	}
	
	/**
	 * Method to perform image registration using bUnwarpJ
	 * Gets the targetGrid, sourceGrid and an instance of bUnwarpJ
	 * the transformationVectors in 2D in x and y direction are saved 
	 * and can be accessed via getters, the registration windows are
	 * closed after each registration of mapping source to target
	 * */
	public void performRegistration(Grid3D targetGrid, Grid3D sourceGrid, bUnwarpJ_ bUnwarpJ){
		
		bUnwarpJ_.setRichOutput(true);
		bUnwarpJ_.setSaveTransformation(true);	
		
		// Iterate over all projection slices
		for(int i = 0; i < targetGrid.getSize()[2]; i++){
			
				System.out.println("Perform registration for slice: " + i);
			
				// Take the current slice from the whole stack and convert it to the ImagePlus format
				Grid2D targetGridSlice = targetGrid.getSubGrid(i);
				Grid2D sourceGridSlice = sourceGrid.getSubGrid(i);				
				
				ImageProcessor target = ImageUtil.wrapGrid2D(targetGridSlice);
				ImageProcessor source = ImageUtil.wrapGrid2D(sourceGridSlice);
				
				ImagePlus ipTarget = new ImagePlus("Target", target);
				ImagePlus ipSource = new ImagePlus("Source", source);				
				
				//TODO Defining the Registration Parameters
				// - accuracy mode (0 - Fast, 1 - Accurate, 2 - Mono)
				// - minimum scale deformation (0 - Very Coarse, 1 - Coarse, 2 - Fine, 3 - Very Fine)
				// - maximum scale deformation (0 - Very Coarse, 1 - Coarse, 2 - Fine, 3 - Very Fine, 4 - Super Fine)
				//  - divergence weight
				// - curl weight
				// - landmark weight
				// - similarity weight
				// - consistency weight
				// - stopping threshold			
				Param parameters = new Param(2, 0 ,3, 4, 0.5, 0.5, 0, 1, 10, 0.01);
				
				//	Starting the actual registration process
				bUnwarpJ_.alignImagesBatch(ipTarget, ipSource, null, null, parameters);
				
				//	Get transformation vectors from GrayscaleResultTileMaker
				double[][] transformation_x = Transformation.getTrafo_x();
				double[][] transformation_y = Transformation.getTrafo_y();				
				
				if(i < targetGrid.getSize()[2]){
					final double transformedImage [][] = new double [transformation_x.length][transformation_x[0].length];
					int stepv = Math.min(Math.max(10, transformation_x.length/15),30);
					int stepu = Math.min(Math.max(10, transformation_x[0].length/15),30);
	
					for (int v=0; v<transformation_x.length; v++){
						for (int u=0; u<transformation_x[0].length; u++){
							transformedImage[v][u]=255;
						}
					}
					
					for (int v=0; v< transformation_x.length; v+=stepv){
						for (int u=0; u< transformation_x[0].length; u+=stepu){
								final double x = transformation_x[v][u];
								final double y = transformation_y[v][u];
									MiscTools.drawArrow(
											transformedImage,
											u,v,(int)Math.round(x),(int)Math.round(y),0,2);
						}
					}
	
					// Set it to the image stack
					FloatProcessor fp=new FloatProcessor(transformation_x[0].length,transformation_x.length);
					for (int v=0; v<transformation_x.length; v++){
						for (int u=0; u<transformation_x[0].length; u++){
							fp.putPixelValue(u, v, transformedImage[v][u]);
						}
					}					

					//Grid2D grid = ImageUtil.wrapFloatProcessor(fp);
					//grid.show("DeformationVectors");
				}				
				
				//	Compute real deformation vectors				
				for(int v = 0; v < transformation_x.length; v ++){
					for(int u = 0; u < transformation_x[0].length; u++){
						
						final double x = transformation_x[v][u];
						final double y = transformation_y[v][u];
						transformation_x[v][u] = (x - u);
						transformation_y[v][u] = (y - v);						
					}
				}
				Grid2D grid2D_x = convertToGrid2D(transformation_x);
				Grid2D grid2D_y = convertToGrid2D(transformation_y);				
			
				this.transformationVectorsX.setSubGrid(i, grid2D_x);
				this.transformationVectorsY.setSubGrid(i, grid2D_y);
								
				IJ.run("Close All");
		}		
	 }

	 /**
	  * Method to configure the geometry
	  */
	public void configure(){
		
		//	Loading the CONRAD config. 
		Configuration.loadConfiguration();
		Configuration cnfg = Configuration.getGlobalConfiguration();
		
		//	Reading the projection matrices.
		Trajectory geo = cnfg.getGeometry();
		Projection[] pMat = readProjections(projMatricesFile);
		double[] angles = new double[pMat.length];
		for(int i = 0; i < pMat.length; i++){
			angles[i] = 200.0d/(pMat.length-1)*i;
		}
		//	Loading the actual geometry. 
		geo.setProjectionMatrices(pMat);
		geo.setNumProjectionMatrices(pMat.length);
		geo.setProjectionStackSize(pMat.length);
		geo.setPrimaryAngleArray(angles);
		geo.setSourceToAxisDistance(sourceIsocenterDist);
		geo.setSourceToDetectorDistance(sourceDetectorDist);
		geo.setDetectorWidth((int)sizeU);
		geo.setDetectorHeight((int)sizeV);
		geo.setPixelDimensionX(spacingU);
		geo.setPixelDimensionY(spacingV);
		cnfg.setGeometry(geo);
		Configuration.setGlobalConfiguration(cnfg);
		this.cnfg = cnfg;
	}	
		
	/**convert array[][] to Grid2D*/
	public static Grid2D convertToGrid2D(double[][] array){

		int height = array.length;
		int width = array[0].length;
		
		Grid2D grid2D = new Grid2D(height, width);
		for(int i = 0; i < height; i++){
			for(int j = 0; j < width; j++){
				double val = array[i][j];
				grid2D.setAtIndex(i, j,(float) val);
			}
		}
		return grid2D;
	}
	
	
	/**Thresholding Method for a Grid3D.
	 * Values above threshold are kept, values below are set to zero.*/
	public Grid3D threshold3D(Grid3D grid, double threshold){
			int[] size = grid.getSize();
			int width = size[0];
			int height = size[1];
			int depth = size[2];
			Grid3D result = new Grid3D(width, height, depth);
			
			for(int i = 0; i < width; i ++){
				for(int j = 0; j < height; j++){
					for(int k = 0; k < depth; k++){
						if(grid.getAtIndex(i, j, k) < threshold){
							result.setAtIndex(i, j, k, 0.0f);
						}
						else{
							result.setAtIndex(i,j,k,grid.getAtIndex(i, j, k));
						}
					}
				}				
			}
			return result;
	}
	

	/**Binarization Method for a Grid3D. 
	 * values above threshold are set to one, values below to zero.*/
	public Grid3D binarize3D(Grid3D grid, double threshold){
			int[] size = grid.getSize();
			int width = size[0];
			int height = size[1];
			int depth = size[2];
			Grid3D result = new Grid3D(width, height, depth);
			
			for(int i = 0; i < width; i ++){
				for(int j = 0; j < height; j++){
					for(int k = 0; k < depth; k++){
						if(grid.getAtIndex(i, j, k) < threshold){
							result.setAtIndex(i, j, k, 0.0f);
						}
						else{
							result.setAtIndex(i, j, k, 1.0f);
						}
					}
				}				
			}
			return result;
	}
	
	/**threshold with percentile threshold*/
	public Grid3D percentileHistogramThreshold(ImagePlus imp, double percentile){
		
		//	Computing the threshold.
		long percentileNum = getPercentileNum(imp, percentile);
		double threshold = findThreshold(imp, percentileNum);
		
		// Performing the actual thresholding with the previously computed threshold value. 
		Grid3D grid = ImageUtil.wrapImagePlus(imp);
		grid = threshold3D(grid, threshold);
		
		return grid;
	}
	
	public long getPercentileNum(ImagePlus imp, double percentile){
		ResultsTable rt = new ResultsTable();
		long[] histogram = null; //is.histogram;
		long totalCount = 0;
		double value = 0.0;
		double binWidth = 0.0;
		double min = imp.getProcessor().getMin();
		StackStatistics ss1 = new StackStatistics(imp);
		min = ss1.histMin;
		if(imp.getBitDepth() == 8 || imp.getBitDepth() == 24){
			for(int i = 0; i < histogram.length; i++){
				rt.setValue("Value", i, i);
			    rt.setValue("Count", i, histogram[i]);
			}
		}else{
			value = min;
			binWidth = ss1.binSize;
			histogram = ss1.getHistogram();
			
			for (int i=0; i<histogram.length; i++) {
		     rt.setValue("Value", i, value);
		     rt.setValue("Count", i, histogram[i]);
		     value += binWidth;
		     totalCount += histogram[i];
			}
	
		}
		long percentileNum = (long) (percentile*totalCount);
		
		return percentileNum;
	}
	
	public double findThreshold(ImagePlus imp, long percentileNum){
		ResultsTable rt = new ResultsTable();
		long[] histogram = null; //is.histogram;
		long totalCount = 0;
		double value = 0.0;
		double binWidth = 0.0;
		StackStatistics ss1 = new StackStatistics(imp);
		double max = ss1.histMax;
		histogram = ss1.getHistogram();
		if(imp.getBitDepth() == 8 || imp.getBitDepth() == 24){
			for(int i = histogram.length - 1; i > - 1; i--){
				rt.setValue("Value", i, i);
			    rt.setValue("Count", i, histogram[i]);
			     if(percentileNum < totalCount){
			    	 return value;
			     }
			}
		}else{
			value = max;
			binWidth = ss1.binSize;
			histogram = ss1.getHistogram();
			
			for(int i = histogram.length - 1; i > - 1; i--){
		     value -= binWidth;
		     totalCount += histogram[i];
		     if(percentileNum < totalCount){
		    	 return value;
		     }
			}
			rt.show("Threshold Results");
		}
		return 0.0;
	}
	/**
	 * A method to perform a morphological filtering of the projections.  
	 * @param projFile : location of the projection file on the disc
	 * @param structElementSize : radius of the structuring element
	 * @param filter : choose which filter should be used. Currently "topHat" or "dilation" are implemented. 
	 * @return ImagePlus with the filtered result.
	 */
	public static ImagePlus morphFilterProjs(String projFile, int structElementSize, String filter){
		
		//	Creating the structuring element.
		StructuringElement s = new StructuringElement("Circle", structElementSize, false);		
		
		// Loading the projection source and creating an ImageStack.
		Grid3D projGrid = ImageUtil.wrapImagePlus(IJ.openImage(projFile));		
		ImagePlus result = null;		
		ImageStack stack = new ImageStack(projGrid.getSize()[0],projGrid.getSize()[1]);
		
		// Applying the morphological filter slice by slice. 
		System.out.println("Applying the morphological filter. This may take a while...");
		for(int i = 0; i < projGrid.getSize()[2]; i++)
		{			
			ImagePlus workPlus = new ImagePlus("Img", ImageUtil.wrapGrid2D(projGrid.getSubGrid(i)));
			if(filter == "topHat") result = Morphology.topHat(workPlus, s, true);
			if(filter == "dilation") result = Morphology.dilate(workPlus, s);
			stack.addSlice(result.getProcessor());			
		}
		
		ImagePlus projIP = new ImagePlus("ProjectionsFiltered", stack);
		return projIP;
	}
	
	
	public int getGat_ign() {
		return gat_ign;
	}

	public void setGat_ign(int gat_ign) {
		this.gat_ign = gat_ign;
	}
	
	private Projection[] readProjections(String filename){
		 
		Projection[] pMat = null;
		try {
			FileInputStream fStream = new FileInputStream(filename);
			// Number of matrices is given as the total size of the file
			// divided by 4 bytes per float, divided by 12 floats per projection matrix
			int nMat = (int) (fStream.getChannel().size() / 4 / (4*3));
			DataInputStream in = new DataInputStream(fStream);
			
			pMat = new Projection[nMat];
			
			for(int m = 0; m < nMat; m++){
				SimpleMatrix mat = new SimpleMatrix(3,4);
				for(int i = 0; i < mat.getRows(); i++){
					for(int j = 0; j < mat.getCols(); j++){
						byte[] buffer = new byte[4];
					    int bytesRead = in.read(buffer);
					    float val = 0;
					    if(bytesRead == 4){
					    	val = ByteBuffer.wrap(buffer).order(ByteOrder.LITTLE_ENDIAN).getFloat();
					    }
					    mat.setElementValue(i, j, val);
					}
				}
				pMat[m] = new Projection(mat);
			}
			
			in.close();
			fStream.close();
			
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		} catch (IOException e) {
			e.printStackTrace();
		}
		
		return pMat;
	}
	
	/**
	 * Reads the ECG from a file. It is assumed that there are as many heart phases in a separate line
 	 * as there are projections.
	 * @param filename
	 * @return
	 */
	private double[] readEcg(String filename){
		ArrayList<Double> e = new ArrayList<Double>();
		FileReader fr;
		try {
			fr = new FileReader(filename);		
			BufferedReader br = new BufferedReader(fr);			
			String line = br.readLine();
			while(line != null){
				StringTokenizer tok = new StringTokenizer(line);
				e.add(Double.parseDouble(tok.nextToken()));
				line = br.readLine();
			}
			br.close();
			fr.close();
			
		} catch (FileNotFoundException e1) {
			e1.printStackTrace();
		} catch (IOException e1) {
			e1.printStackTrace();
		}
		
		double[] ecg = new double[e.size()];
		for(int i = 0; i < e.size(); i++){
			ecg[i] = e.get(i);
		}
		return ecg;
	}	
	/**
	 * Uses the ECGGating class to perform cosine-function based gating of the acquisition.
	 * @param img
	 * @param ecg
	 * @return
	 */
	private Grid3D applyGating(Grid3D img, double[] ecg){
		ECGGating gating = new ECGGating(gatingWidth, gatingSlope, REF_HEART_PHASE);
		return gating.weightProjections(img, ecg);
	}
	/**
	 * Returns the registration save path.
	 */
	public static String getRegistrationPath(){
		return registrationPath;
	}
	
	/**
	 * This method "crops" the image to the specified roi.
	 * The image doesn't get cropped actually, but all values outside the roi are set to zero.
	 */
	public Grid3D cropImage(Grid3D inputGrid, int x, int y, int width, int height ){
		
		for(int k = 0; k < inputGrid.getSize()[2]; k++){
			for(int i = 0; i < inputGrid.getSize()[0]; i++){
				for(int j = 0; j < inputGrid.getSize()[1]; j++){
					if(i < y || i > y + height || j < x || j > x+width){						
							inputGrid.setAtIndex(i, j, k, 0.0f);
					}
				}
			}
		}		
		return inputGrid;
	}
	
}