content: /*
   * Copyright (C) 2010-2014 Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  package edu.stanford.rsl.apps.gui;
  
  import ij.IJ;
  import ij.ImagePlus;
  import ij.ImageStack;
  import ij.measure.Calibration;
  
  import java.io.FileOutputStream;
  import java.io.ObjectOutputStream;
  import java.util.ArrayList;
  import java.util.Collections;
  
  import edu.stanford.rsl.conrad.geometry.General;
  import edu.stanford.rsl.conrad.geometry.motion.MotionUtil;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
  import edu.stanford.rsl.conrad.geometry.splines.SurfaceBSpline;
  import edu.stanford.rsl.conrad.geometry.splines.SurfaceBSplineVolumePhantom;
  import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  import edu.stanford.rsl.conrad.phantom.AnalyticPhantom4D;
  import edu.stanford.rsl.conrad.phantom.renderer.MetricPhantomRenderer;
  import edu.stanford.rsl.conrad.phantom.xcat.SquatScene;
  import edu.stanford.rsl.conrad.phantom.xcat.XCatScene;
  import edu.stanford.rsl.conrad.rendering.PrioritizableScene;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.Configuration;
  import edu.stanford.rsl.conrad.utils.ImageGridBuffer;
  import edu.stanford.rsl.conrad.utils.RegKeys;
  import edu.stanford.rsl.conrad.utils.UserUtil;
  
  public class XCatMetricPhantomCreator {
  
  	private int steps = 2;
  
  	/**
  	 * Method to render 3-D versions of the phantoms.
  	 * @param args the arguments (none)
  	 */
  	public static void main(String [] args){
  
  		CONRAD.setup();
  		XCatMetricPhantomCreator creator = new XCatMetricPhantomCreator();
  		AnalyticPhantom4D scene = creator.instantiateScene();
  		ImagePlus hyper = creator.renderMetricVolumePhantom(scene);
  		hyper.show();
  	}
  
  	public void renderMetricVolumePhantom(){
  		renderMetricVolumePhantom(null);
  	}
  
  	public AnalyticPhantom4D instantiateScene(){
  
  		try {
  			steps = UserUtil.queryInt("Enter number of steps for sampling:", steps);
  		} catch (Exception e2) {
  			// TODO Auto-generated catch block
  			e2.printStackTrace();
  		}
  		
  		AnalyticPhantom4D scene = null;
  		boolean filenameSet = Configuration.getGlobalConfiguration().getRegistry().get(RegKeys.SPLINE_4D_LOCATION) != null;
  		//select the scene
  		if (filenameSet){
  			scene = (AnalyticPhantom4D) MotionUtil.get4DSpline();
  		}
  
  
  		if(scene == null){
  			try {
  				scene = (AnalyticPhantom4D) UserUtil.queryPhantom("Scene Selection", "Please select a 4D phantom:");
  			} catch (Exception e1) {
  				// TODO Auto-generated catch block
  				e1.printStackTrace();
  			}
  			try {
  				scene.configure();
  				boolean write = true;
  				if (Configuration.getGlobalConfiguration().getRegistry().get(RegKeys.SPLINE_4D_LOCATION) == null) write = false;
  				if (write) {
  					FileOutputStream fos = new FileOutputStream(Configuration.getGlobalConfiguration().getRegistry().get(RegKeys.SPLINE_4D_LOCATION));
  					ObjectOutputStream oos = new ObjectOutputStream(fos);
  					oos.writeObject(scene);
  					oos.flush();
  					oos.close();
  				}
  			} catch (Exception e) {
  				// TODO Auto-generated catch block
  				e.printStackTrace();
  			}
  		}
  		return scene;
  	}
  
  	public ImagePlus renderMetricVolumePhantom(AnalyticPhantom4D scene){
  		try {
  			// read global configuration
  			Trajectory geom = Configuration.getGlobalConfiguration().getGeometry();
  			int dimz = geom.getReconDimensionZ();
  			int dimx = geom.getReconDimensionX();
  			int dimy = geom.getReconDimensionY();
  			// create hyperstack (4D volume), 3D volume that has time line
  			ImageStack hyperStack = new ImageStack(dimx, dimy);
  
  
  			//scene.setTimeWarper(new HarmonicTimeWarper(1));
  			ArrayList<SurfaceBSpline> list = null;//scene.getSplines();
  			if (scene instanceof XCatScene) {
  				list = ((XCatScene)scene).getSplines();
  			}
  			Calibration calibration = null;
  
  						// render the scene
  			// the number of steps gives the number of volumes rendered. (k < 40)
  
  			double[] timeList = new double[getSteps()];
  			for (int i = 0; i < timeList.length; i++) {
  				double steps = (getSteps()==1)?1.0:((double)(getSteps()-1));
  				timeList[i]=((double)i)/steps;
  			}
  
  
  			for (int k = 0; k < getSteps(); k++){
  				ImageGridBuffer buffer = new ImageGridBuffer();
  				IJ.showStatus("Rendering State " + k);
  				IJ.showProgress(((double)k)/getSteps());
  
  				// configure worker
  				// worker: renders the scene, it basically cast a ray thru the scene for every pixel for every projection
  				SurfaceBSplineVolumePhantom phantomWorker = new SurfaceBSplineVolumePhantom();
  				phantomWorker.setImageProcessorBuffer(buffer);
  				phantomWorker.setShowStatus(true);
  				boolean resize = true;
  				if (Configuration.getGlobalConfiguration().getRegistry().get(RegKeys.XCAT_AUTO_RESIZE) != null){
  					if (Configuration.getGlobalConfiguration().getRegistry().get(RegKeys.XCAT_AUTO_RESIZE).equals("false")) {
  						resize = false;
  					} else {
  						// set scene boundaries to recon volume.
  						// Thus we have to create a new origin location that centers the recon volume around the scene volume without changing the voxel spacing.
  						PointND origin = new PointND(scene.getMax());
  						origin.getAbstractVector().subtract(scene.getMin().getAbstractVector());
  						origin.getAbstractVector().divideBy(2.0);
  						origin.getAbstractVector().add(scene.getMin().getAbstractVector());
  						PointND maximumCoordinate = new PointND(Configuration.getGlobalConfiguration().getGeometry().getReconVoxelSizes()[0]*Configuration.getGlobalConfiguration().getGeometry().getReconDimensionX(),
  								Configuration.getGlobalConfiguration().getGeometry().getReconVoxelSizes()[1]*Configuration.getGlobalConfiguration().getGeometry().getReconDimensionY(),
  								Configuration.getGlobalConfiguration().getGeometry().getReconVoxelSizes()[2]*Configuration.getGlobalConfiguration().getGeometry().getReconDimensionZ());
  						SimpleVector reg = maximumCoordinate.getAbstractVector().dividedBy(2.0);
  						origin.getAbstractVector().subtract(reg);
  						maximumCoordinate.getAbstractVector().add(origin.getAbstractVector());
  						Configuration.getGlobalConfiguration().getGeometry().setOriginInWorld(origin);
  						scene.setMin(origin);
  						scene.setMax(maximumCoordinate);
  					}
  				}
  				if (list != null && list.size() > 0) {
  					phantomWorker.setSplineList(list);
  				}
  
  				/**
  				 * Time calculation
  				 * 
  				 * * .4 => Inspriration only
  				 */
  				double time = 0;
  				if (timeList == null){
  					time= (k/(double)getSteps());
  				}
  				else{
  					time = timeList[k];
  				}
  				System.out.println("Current timepoint: "+ time);
  
  
  				if (scene instanceof SquatScene){
  					// Update Origin according to the knee center point
  					/*
  					PointND centeredOrigin = new PointND(269.3, 287.2, -638.3);
  					PointND origin = new PointND(0,0,0);
  					for (int i = 0; i < geom.getReconDimensions().length; i++) {
  						double size = geom.getReconDimensions()[i]*geom.getReconVoxelSizes()[i];
  						origin.set(i, centeredOrigin.get(i)-(size/2.0-0.5));
  					}
  
  
  
  					PointND leftTibiaTop = new PointND(7.6, 287.2,-638.3);
  					PointND rightTibiaTop = new PointND(269.3, 287.2, -638.3);
  					SimpleVector centerOffset = SimpleOperators.subtract(rightTibiaTop.getAbstractVector(), leftTibiaTop.getAbstractVector()).dividedBy(2.0);
  					PointND initialOrigin = new PointND(SimpleOperators.add(leftTibiaTop.getAbstractVector(), centerOffset));
  					 */
  					double yoffset = 50;
  					PointND rightTibiaTop = new PointND(269.3, 287.2, -638.3);
  					PointND initialOrigin = rightTibiaTop;
  
  					PointND currentOrigin = scene.getPosition(initialOrigin, 0, time);
  					// update the origin, i.e. the volume position to the current knee position
  					System.out.println("XCat centered origin: " + currentOrigin);
  
  					PointND currentOriginCorner = currentOrigin.clone();
  					for (int i = 0; i < geom.getReconDimensions().length; i++) {
  						double size = geom.getReconDimensions()[i]*geom.getReconVoxelSizes()[i];
  						currentOriginCorner.set(i, currentOrigin.get(i)-(size/2.0-0.5));
  					}
  					currentOriginCorner.set(1, currentOriginCorner.get(1)+yoffset);
  					geom.setOriginInWorld(currentOriginCorner);
  
  
  
  					System.out.println("Voxel Coordinate System centered origin: " +  new SimpleVector(General.worldToVoxel(currentOrigin.getCoordinates(), geom.getReconVoxelSizes(), currentOriginCorner.getCoordinates())));
  					System.out.println("Voxel Coordinate System moving origin: " +  new SimpleVector(General.worldToVoxel(rightTibiaTop.getCoordinates(), geom.getReconVoxelSizes(), currentOriginCorner.getCoordinates())));
  					System.out.println(" ");
  					System.out.println(" ");
  					//SimpleOperators.subtract(currentOrigin.getAbstractVector(),initialOrigin.getAbstractVector()));
  
  
  					/*
  					PointND leftFemurTop = new PointND(16.6, 273.6,-207.42);
  					/*PointND rightFemurTop = new PointND(269.3, 273.6,-207.42);
  					SimpleVector rotCenterOffset = SimpleOperators.subtract(rightFemurTop.getAbstractVector(), leftFemurTop.getAbstractVector()).dividedBy(2.0);
  					PointND rotCenter = new PointND(SimpleOperators.add(leftFemurTop.getAbstractVector(), rotCenterOffset));
  
  					double lowerZoffset = 0;
  					PointND leftTibiaTop = new PointND(7.6, 287.2,-638.3+lowerZoffset);
  					PointND rightTibiaTop = new PointND(269.3, 287.2, -638.3+lowerZoffset);
  					SimpleVector centerOffset = SimpleOperators.subtract(rightTibiaTop.getAbstractVector(), leftTibiaTop.getAbstractVector()).dividedBy(2.0);
  					PointND initialOrigin = new PointND(SimpleOperators.add(leftTibiaTop.getAbstractVector(), centerOffset));
  
  
  
  					MotionField rotMot = new RotationMotionField(leftFemurTop , new SimpleVector(1,0,0), General.toRadians(0));
  					PointND movingOrigin = rotMot.getPosition(initialOrigin, 0, time);
  					PointND currentOrigin = new PointND(new double[movingOrigin.getDimension()]);
  
  					// determine the current coordinate system origin by the moving center point
  					for (int i = 0; i < geom.getReconDimensions().length; i++) {
  						double size = geom.getReconDimensions()[i]*geom.getReconVoxelSizes()[i];
  						currentOrigin.set(i, movingOrigin.get(i)-(size/2.0-0.5));
  					}
  					System.out.println(currentOrigin);
  					// update the origin, i.e. the volume position to the current knee position
  					geom.setOriginInWorld(currentOrigin);*/
  
  				}
  
  
  
  				if (resize){
  					if (scene.getMin() == null){
  						phantomWorker.resizeVolumeToMatchSplineSpace();
  					} else {
  						phantomWorker.resizeVolumeToMatchBounds(scene.getMin(), scene.getMax());
  					}
  				} else {
  					phantomWorker.setBoundsFromGeometry(scene);
  				}
  
  
  				// tessellate scene -> convert to triangles -> render the triangles
  				PrioritizableScene current = scene.getScene(time);
  
  				phantomWorker.setScene(current);
  
  				// create renderer
  				MetricPhantomRenderer phantom = new MetricPhantomRenderer();
  				ArrayList<Integer> processors = new ArrayList<Integer>();
  				for (int i = 0; i < dimz; i++){
  					processors.add(new Integer(i));
  				}
  				phantomWorker.setSliceList(Collections.synchronizedList(processors).iterator());
  
  				// render
  				phantom.setModelWorker(phantomWorker);
  				phantom.createPhantom();
  				ImagePlus renderedBSpline = buffer.toImagePlus("State "+k);
  				calibration = renderedBSpline.getCalibration();
  
  				// the volume to the hyperstack
  				for (int i=1; i<= dimz; i++){
  					hyperStack.addSlice("Slice z = " +(i-1) + " t = " + k, renderedBSpline.getStack().getProcessor(i));
  				}
  
  			}
  
  
  
  			calibration.xOrigin = Configuration.getGlobalConfiguration().getGeometry().getOriginInPixelsX();
  			calibration.yOrigin = Configuration.getGlobalConfiguration().getGeometry().getOriginInPixelsY();
  			calibration.zOrigin = Configuration.getGlobalConfiguration().getGeometry().getOriginInPixelsZ();
  			calibration.pixelWidth = Configuration.getGlobalConfiguration().getGeometry().getVoxelSpacingX();
  			calibration.pixelHeight = Configuration.getGlobalConfiguration().getGeometry().getVoxelSpacingY();
  			calibration.pixelDepth = Configuration.getGlobalConfiguration().getGeometry().getVoxelSpacingZ();
  			// finalize the hyperstack
  			ImagePlus hyper = new ImagePlus();
  			hyper.setCalibration(calibration);
  			hyper.setStack(scene.getName(), hyperStack);
  			hyper.setDimensions(1, dimz, getSteps());
  			hyper.setOpenAsHyperStack(true);
  			IJ.showProgress(1.0);
  
  			//geom.setOriginInWorld(saveWorldOrigin);
  
  			// display the result.
  			return hyper;
  
  		} catch (Exception e) {
  			// TODO Auto-generated catch block
  			e.printStackTrace();
  			return null;
  		}
  	}
  
  	/**
  	 * @return the steps
  	 */
  	public int getSteps() {
  		return steps;
  	}
  
  	/**
  	 * @param steps the steps to set
  	 */
  	public void setSteps(int steps) {
  		this.steps = steps;
  	}
  }
  