content: /*
   * Copyright (C) 2014 Michael Manhart, Kerstin Mueller
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  
  package edu.stanford.rsl.tutorial.phantoms;
  
  import java.io.IOException;
  import java.util.ArrayList;
  
  import ij.ImageJ;
  import ij.io.FileInfo;
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.Grid3D;
  import edu.stanford.rsl.conrad.data.numeric.NumericGridOperator;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.conrad.data.numeric.iterators.NumericPointwiseIteratorND;
  import edu.stanford.rsl.conrad.geometry.Projection;
  import edu.stanford.rsl.conrad.geometry.trajectories.ConfigFileBasedTrajectory;
  import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
  import edu.stanford.rsl.conrad.io.GridRawIOUtil;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.numerics.SimpleOperators;
  import edu.stanford.rsl.conrad.opencl.OpenCLForwardProjectorDynamicVolume;
  import edu.stanford.rsl.conrad.physics.materials.Material;
  import edu.stanford.rsl.conrad.physics.materials.database.MaterialsDB;
  import edu.stanford.rsl.conrad.physics.materials.utils.AttenuationType;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.Configuration;
  import edu.stanford.rsl.conrad.utils.XmlUtils;
  
  
  /**
   * @author Kerstin Mueller
   * BrainPerfusionPhantom allows to generate noise-free 2D projection images of a brain perfusion scan
   */
  public class BrainPerfusionPhantom {
  	BrainPerfusionPhantomConfig config;
  
  	// forward and backward sweeps, calibrated with projection matrices
  	protected Trajectory gFwd;
  	protected Trajectory gBwd;
  
  	// pre-generated data from the MATLAB Script which can be downloaded here: http://www5.cs.fau.de/en/research/data/digital-brain-perfusion-phantom/
  	protected Grid3D volSkull = null;
  	protected Grid3D volTissue = null;
  
  	// variables for the contrast flow simulation
  	protected Grid3D  volContrastPrev = null;
  	protected Grid3D  volContrastNext = null;
  	protected int     volContrastPrevIdx = -1;
  	protected int     volContrastNextIdx = -1;
  
  	OpenCLForwardProjectorDynamicVolume openclFwdProjector;
  
  	protected static NumericGridOperator gop = NumericGridOperator.getInstance();
  	
  	public static void main(String[] args) {
  		CONRAD.setup();
  
  		
  		
  		// Create new config file for perfusion phantom
  		BrainPerfusionPhantomConfig perfConfig = new BrainPerfusionPhantomConfig();
  		if (args.length >= 5)
  		{
  			// the directory containing, skull.raw and all generated volumes with MATLAB
  			perfConfig.phantomDirectory = args[0];
  			// complete path with filename to projection matrix for forward run
  			perfConfig.calibrationFwdMatFile = args[1];
  			// complete path with filename to projection matrix for backward run
  			perfConfig.calibrationBwdMatFile = args[2];
  			// the output directory for all projections
  			perfConfig.phantomOutDirectory = args[3];
  			// if motion should be applied
  			perfConfig.motion = Boolean.parseBoolean(args[4]);
  			// sampling of the generated volumes
  			perfConfig.phantomSampling = 1.0f;
  			if (args.length==6) {
  				if (!args[5].isEmpty())
  				{
  					// If 3D Markers should be added to the skull phantom
  					perfConfig.markerFile = args[5];
  				}
  			}
  		}  else {
  			System.err.println("Usage -directory of the phantom - fwd projection matrix file [.txt]"
  					+ "- bwd projection matrix file [.txt] - output directory -motion[true/false]");
  			System.exit(1);
  		}
  
  		BrainPerfusionPhantom perfPhantom = new BrainPerfusionPhantom(perfConfig);
  
  		if ( !perfConfig.markerFile.isEmpty())
  		{
  			perfPhantom.addMarkers();
  		}
  		
  		try {
  			perfPhantom.createProjectionData();
  		}catch(IOException io){
  			System.err.println("PerfusionBrainPhantom.createProjectionData(): IO error!");
  			io.printStackTrace();
  		}
  
  	}
  
  	/**
  	 * Adds 3D markers to the skull volume, which are predefined in an xml file
  	 * 
  	 */
  	public void addMarkers()
  	{
  		// Add markers
  		// load our 3d marker positions
  		ArrayList<double[]> threeDMarkerPositions=null;
  		try {
  			threeDMarkerPositions = (ArrayList<double[]>) XmlUtils.importFromXML(config.markerFile);
  		} catch (Exception e) {
  			// TODO Auto-generated catch block
  			e.printStackTrace();
  			System.out.println("Exception while reading xml file: " + e);
  		}	
  		
  		Grid3D tempSkullVolume = getSkullVolume();
  		Grid3D tempTissueVolume = getTissueVolume();
  		
  		// Create the tantalum markers with the correct attenuation
  		float radius_mm = 0.0f;// mm
  		Material beads = MaterialsDB.getMaterial("tantalum");
  		float attenutation_beads = (float)(beads.getAttenuation(80, AttenuationType.TOTAL_WITH_COHERENT_ATTENUATION))/10;
  		
  		int radius_px = (int) (radius_mm / BrainPerfusionPhantomConfig.voxelSize);
  		for ( int i = 0; i < threeDMarkerPositions.size(); i++)
  		{
  			int indexX = (int)(threeDMarkerPositions.get(i)[0]);
  			int indexY = (int)(threeDMarkerPositions.get(i)[1]);
  			int indexZ = (int)(threeDMarkerPositions.get(i)[2]);
  			
  			for (int rx =-radius_px; rx <= radius_px; rx++)
  			{
  				for (int ry =-radius_px; ry <= radius_px; ry++)
  				{
  					for (int rz =-radius_px; rz <= radius_px; rz++)
  					{
  						tempSkullVolume.addAtIndex(indexX+rx, indexY+ry, indexZ+rz, attenutation_beads);
  					}
  				}
  			}
  		}
  		setSkullVolume(tempSkullVolume);
  		gop.addBySave(tempSkullVolume,tempTissueVolume);
  		tempSkullVolume.show("Skull with attached markers");
  	}
  	
  	/**
  	 * Creates the 2D projection data
  	 * @throws IOException
  	 */
  	public void createProjectionData() throws IOException
  	{
  		// -----------------------------------------------------
  		// material decomposed dynamic forward projection 
  		// -----------------------------------------------------		
  		dynamicDensityForwardProjection();
  	}
  
  	/**
  	 * Performs the dynamic density forward projection in forward, backward runs divided for tissue, skull and contrast agent
  	 * @throws IOException
  	 */
  	protected void dynamicDensityForwardProjection() throws IOException
  	{
  		openclFwdProjector = new OpenCLForwardProjectorDynamicVolume();
  		int[] volumeSize  = new int[3];
  		float[] voxelSize = new float[3];
  		volumeSize[0]=BrainPerfusionPhantomConfig.sizeX;
  		volumeSize[1]=BrainPerfusionPhantomConfig.sizeY;
  		volumeSize[2]=BrainPerfusionPhantomConfig.sizeZ;
  		voxelSize[0]=BrainPerfusionPhantomConfig.voxelSize;
  		voxelSize[1]=BrainPerfusionPhantomConfig.voxelSize;
  		voxelSize[2]=BrainPerfusionPhantomConfig.voxelSize;
  		openclFwdProjector.configure(gFwd, gBwd, volumeSize, voxelSize);
  
  		Grid2D[] projections = new Grid2D[gFwd.getNumProjectionMatrices()];
  		// skull // already sets the projection matrices to the graphics card
  		setVolumeToProject(volSkull);
  		
  		float rot = 30.0f;
  		float [] trans = {0,0,0};			
  		for(int projection_id = 0; projection_id < gFwd.getNumProjectionMatrices(); projection_id++)
  		{
  			if (config.motion && projection_id > 120 )
  			{
  				Projection proj = gFwd.getProjectionMatrix(projection_id);
  				Projection projMat = applyRigidMotion(proj, rot, trans);
  				gFwd.setProjectionMatrix(projection_id, projMat);
  			}				
  			System.out.println("Forward projection: skull forward rotation projection " + (projection_id+1));	
  			projections[projection_id] = applyForwardProjection(projection_id, true, config.motion);
  		}
  		saveProjectionData("skull_fwd", projections);
  		for(int projection_id = 0; projection_id < gBwd.getNumProjectionMatrices(); projection_id++)
  		{
  			if (config.motion && projection_id > 120 )
  			{
  				Projection proj = gFwd.getProjectionMatrix(projection_id);
  				Projection projMat = applyRigidMotion(proj, rot, trans);
  				gFwd.setProjectionMatrix(projection_id, projMat);
  			}	
  			System.out.println("Forward projection: skull backward rotation projection " + (projection_id+1));
  			projections[projection_id] = applyForwardProjection(projection_id, false, config.motion);
  		}				
  		saveProjectionData("skull_bwd", projections);
  
  		//tissue
  		setVolumeToProject(volTissue);
  		for(int projection_id = 0; projection_id < gFwd.getNumProjectionMatrices(); projection_id++)
  		{
  			if (config.motion && projection_id > 120 )
  			{
  				Projection proj = gFwd.getProjectionMatrix(projection_id);
  				Projection projMat = applyRigidMotion(proj, rot, trans);
  				gFwd.setProjectionMatrix(projection_id, projMat);
  			}	
  			System.out.println("Forward projection: tissue forward rotation projection " + (projection_id+1));
  			projections[projection_id] = applyForwardProjection(projection_id, true, config.motion);
  		}
  		saveProjectionData("tissue_fwd", projections);
  		for(int projection_id = 0; projection_id < gBwd.getNumProjectionMatrices(); projection_id++)
  		{
  			if (config.motion && projection_id > 120 )
  			{
  				Projection proj = gFwd.getProjectionMatrix(projection_id);
  				Projection projMat = applyRigidMotion(proj, rot, trans);
  				gFwd.setProjectionMatrix(projection_id, projMat);
  			}	
  			System.out.println("Forward projection: tissue backward rotation projection " + (projection_id+1));
  			projections[projection_id] = applyForwardProjection( projection_id, false, config.motion);
  		}
  		saveProjectionData("tissue_bwd", projections);
  
  		// contrast agent
  		double time_increment_projection = config.tRot/(gFwd.getNumProjectionMatrices()-1);
  		for(int sweep = 0; sweep < config.numSweeps; sweep++)
  		{
  			boolean fwd = (sweep%2 == 0);
  			float time = config.tStart+sweep*(config.tRot+config.tPause);
  			for(int projection_id = 0; projection_id < gFwd.getNumProjectionMatrices(); projection_id++)
  			{
  				if (config.motion && projection_id > 120 )
  				{
  					Projection proj = gFwd.getProjectionMatrix(projection_id);
  					Projection projMat = applyRigidMotion(proj, rot, trans);
  					gFwd.setProjectionMatrix(projection_id, projMat);
  				}	
  				System.out.println("Forward projection contrast agent: sweep " + (sweep+1) + " projection " + (projection_id+1));
  				setVolumeToProject(getDynamicContrastVolume(time));
  				projections[projection_id] = applyForwardProjection(projection_id, fwd, config.motion);
  				time += time_increment_projection;
  			}
  			String fn_proj_out = String.format("sweep%03d_contrast", sweep);
  			saveProjectionData(fn_proj_out, projections);
  		}	
  
  	}
  
  	/**
  	 * Applies 3D rotation and translation to the projection matrix 
  	 * y = P * [R|t] * x; P(3x4) and R|t (4x4)
  	 * @param proj projection where to apply the motion
  	 * @param angle_degree
  	 * @param translation float[]   
  	 * @return projection with changed projection matrix
  	 */
  	protected Projection applyRigidMotion(Projection proj, float angle_degree, float [] translation)
  	{
  		double angle_radian = 0.0;
  		angle_radian = angle_degree*(Math.PI/180.0);
  			
  		SimpleMatrix P = proj.computeP();
  		SimpleMatrix Motion = new SimpleMatrix(4,4);
  		Motion.fill(0.0);
  		Motion.setElementValue(0, 0, Math.cos(angle_radian)); // angle in radians
  		Motion.setElementValue(0, 1, Math.sin(angle_radian)); // angle in radians
  		Motion.setElementValue(1, 0, -Math.sin(angle_radian)); // angle in radians
  		Motion.setElementValue(1, 1, Math.cos(angle_radian)); // angle in radians
  		Motion.setElementValue(2, 2, 1); 
  		Motion.setElementValue(0, 3, translation[0]);
  		Motion.setElementValue(1, 3, translation[1]);
  		Motion.setElementValue(2, 3, translation[2]);
  		Motion.setElementValue(3, 3, 1);
  		SimpleMatrix newProj = SimpleOperators.multiplyMatrixProd(P, Motion);
  		Projection projNew = new Projection();
  		projNew.initFromP(newProj);
  		return projNew; 
  	}
  
  	/**
  	 * Sets the 3D tissue volume for projection 
  	 * @param 
  	 */
  	protected void setTissueVolume(Grid3D tissueVolume)
  	{
  		volTissue = tissueVolume;
  	}
  	
  	
  	/**
  	 * Gets the 3D skull volume for projection 
  	 * @param 
  	 */
  	protected Grid3D getTissueVolume()
  	{
  		return volTissue;
  	}
  	
  	/**
  	 * Sets the 3D skull volume for projection 
  	 * @param 
  	 */
  	protected void setSkullVolume(Grid3D skullVolume)
  	{
  		volSkull = skullVolume;
  	}
  	
  	
  	/**
  	 * Gets the 3D skull volume for projection 
  	 * @param 
  	 */
  	protected Grid3D getSkullVolume()
  	{
  		return volSkull;
  	}
  	
  	/**
  	 * Sets the 3D volume for projection 
  	 * @param volumeBuffer Grid3D containing the information for forward projection
  	 */
  	protected void setVolumeToProject(Grid3D volumeBuffer)
  	{
  		float[] linearBuffer = new float[volumeBuffer.getSize()[0]*volumeBuffer.getSize()[1]*volumeBuffer.getSize()[2]];
  		int i = 0;
  		for(Grid2D grid : volumeBuffer.getBuffer()){
  			for(int j = 0; j < grid.getBuffer().length; j++){
  				linearBuffer[j+i] = grid.getBuffer()[j];
  			}
  			i += grid.getBuffer().length;
  		}
  		openclFwdProjector.setVolume(linearBuffer);
  	}
  
  	/**
  	 * Performs the forward projection onto a given projection image
  	 * @param projectionId the number for the projection image and matrix
  	 * @param fwd defines if it is a forward or backward run
  	 * @return 2D projection image 
  	 */
  	protected Grid2D applyForwardProjection(int projectionId, boolean fwd, boolean motion )
  	{	
  		Trajectory g = fwd?gFwd:gBwd;
  		Grid2D projection = new Grid2D(g.getDetectorWidth(),g.getDetectorHeight());
  		openclFwdProjector.applyForwardProjection(projectionId, fwd, projection.getBuffer(), motion);
  		return projection;	
  	}
  
  	/**
  	 * Saves the projection image stack to file
  	 * @param filename the name of the file (note: not complete path!)
  	 * @param grid contains an array of all projection images
  	 */
  	public void saveProjectionData(String filename, Grid2D[] grid)
  	{
  		FileInfo fI = GridRawIOUtil.getDefaultFloat32BigEndianFileInfo(grid[0]);
  		String sep = System.getProperty("file.separator");
  		String complFilename = config.phantomOutDirectory+sep+filename;
  		try {
  			GridRawIOUtil.saveRawDataGrid(grid, fI, complFilename);
  		} catch (IOException e) {
  			System.err.println("Could not write projection data to file "+ filename);
  			e.printStackTrace();
  		}
  	}
  
  	/**
  	 * Computes the dynamic contrasted volume inside the brain tissue
  	 * @param time the time point in s 
  	 * @return the 3D volume of the contrasted brain
  	 * @throws IOException
  	 */
  	protected Grid3D getDynamicContrastVolume(float time) throws IOException
  	{
  		Grid3D volume = new Grid3D(BrainPerfusionPhantomConfig.sizeX,BrainPerfusionPhantomConfig.sizeY,BrainPerfusionPhantomConfig.sizeZ);
  		FileInfo fI = GridRawIOUtil.getDefaultFloat32BigEndianFileInfo();
  		fI.nImages = BrainPerfusionPhantomConfig.sizeZ;
  
  		int prev_idx = (int)(time/config.phantomSampling)+1;
  		float prev_time = (float)((prev_idx-1)*config.phantomSampling);
  		int next_idx = (int)(time/config.phantomSampling)+2;
  		float next_time = (float)((next_idx-1)*config.phantomSampling);
  
  		String sep = System.getProperty("file.separator");
  
  		if(prev_idx != volContrastPrevIdx) {
  			String complFilename = config.phantomDirectory+sep+Integer.toString(prev_idx);
  			volContrastPrev = (Grid3D) GridRawIOUtil.loadFromRawData(fI, complFilename);
  			HUToAttenuationValues(volContrastPrev);
  			volContrastPrevIdx = prev_idx;
  		}
  		if(next_idx != volContrastNextIdx) {
  			String complFilename = config.phantomDirectory+sep+Integer.toString(next_idx);
  			volContrastNext = (Grid3D) GridRawIOUtil.loadFromRawData(fI, complFilename);
  			HUToAttenuationValues(volContrastNext);
  			volContrastNextIdx = next_idx;
  		}
  
  		float weight = (float)(time-prev_time)/(next_time-prev_time);
  
  		NumericPointwiseIteratorND pIter = new NumericPointwiseIteratorND(volume);
  		NumericPointwiseIteratorND pIterNext = new NumericPointwiseIteratorND(volContrastNext);
  		NumericPointwiseIteratorND pIterPrev = new NumericPointwiseIteratorND(volContrastPrev);
  
  		while( pIter.hasNext() && pIterPrev.hasNext() && pIterNext.hasNext())
  		{
  			float perfusionValue = (1-weight)*pIterPrev.getNext() + weight*pIterNext.getNext();
  			pIter.setNext(perfusionValue);
  		}
  		return volume;
  	}
  
  	
  	public void HUToAttenuationValues(Grid3D volume)
  	{
  		NumericPointwiseIteratorND pIter = new NumericPointwiseIteratorND(volume);
  		while (pIter.hasNext())
  		{
  			// :1000 + 1: scale from HU values to density (0 HU = 1 mg/mm^3; voxel volume is 1 mm^3)
  			pIter.setNext(pIter.get()/1000 + 1);
  		}
  	}
  	
  	/** 
  	 * Reads the skull file generated with the MATLAB Tool (crate_phantom_material_decomposition.m)
  	 */
  	public void readSkullInfo()  
  	{
  		FileInfo fi = GridRawIOUtil.getDefaultFloat32BigEndianFileInfo();
  		fi.width = BrainPerfusionPhantomConfig.sizeX;
  		fi.height = BrainPerfusionPhantomConfig.sizeY;
  		fi.nImages = BrainPerfusionPhantomConfig.sizeZ;
  		String sep = System.getProperty("file.separator");
  		String complFilename = config.phantomDirectory+sep+"skull";
  		volSkull = (Grid3D) GridRawIOUtil.loadFromRawData(fi, complFilename);
  		HUToAttenuationValues(volSkull);
  		
  	}
  
  	/**
  	 * Reads the tissue file generated with the MATLAB Tool (crate_phantom_material_decomposition.m)
  	 */
  	public void readTissueInfo()  
  	{
  		FileInfo fi = GridRawIOUtil.getDefaultFloat32BigEndianFileInfo();
  		fi.width = BrainPerfusionPhantomConfig.sizeX;
  		fi.height = BrainPerfusionPhantomConfig.sizeY;
  		fi.nImages = BrainPerfusionPhantomConfig.sizeZ;
  		String sep = System.getProperty("file.separator");
  		String complFilename = config.phantomDirectory+sep+"tissue";
  		volTissue = (Grid3D) GridRawIOUtil.loadFromRawData(fi, complFilename);
  		HUToAttenuationValues(volTissue);
  	}
  
  	/**
  	 * Constructor
  	 * @param config contains all information regarding the phantom generation process
  	 */
  	public BrainPerfusionPhantom(BrainPerfusionPhantomConfig config)
  	{
  		this.config = config;
  
  		gFwd = ConfigFileBasedTrajectory.openAsGeometrySource(config.calibrationFwdMatFile, Configuration.getGlobalConfiguration().getGeometry());
  		gBwd = ConfigFileBasedTrajectory.openAsGeometrySource(config.calibrationFwdMatFile, Configuration.getGlobalConfiguration().getGeometry());
  
  
  		readSkullInfo();
  		readTissueInfo();
  		if ( volSkull==null || volTissue==null) 
  		{
  			System.err.println("BrainPerfusionPhantom: Error reading phantom input data!");
  		}
  		new ImageJ();
  		
  	}
  
  }