SurfaceBSplineVolumePhantom.java

package edu.stanford.rsl.conrad.geometry.splines;

import java.io.IOException;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Iterator;

import ij.IJ;
import ij.gui.GenericDialog;
import ij.process.FloatProcessor;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.geometry.AbstractShape;
import edu.stanford.rsl.conrad.geometry.shapes.simple.Edge;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
import edu.stanford.rsl.conrad.geometry.shapes.simple.VectorPoint3D;
import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
import edu.stanford.rsl.conrad.geometry.transforms.Translation;
import edu.stanford.rsl.conrad.numerics.SimpleOperators;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.conrad.phantom.AnalyticPhantom4D;
import edu.stanford.rsl.conrad.phantom.workers.SliceWorker;
import edu.stanford.rsl.conrad.phantom.xcat.XCatScene;
import edu.stanford.rsl.conrad.physics.PhysicalObject;
import edu.stanford.rsl.conrad.physics.absorption.AbsorptionModel;
import edu.stanford.rsl.conrad.physics.absorption.SelectableEnergyMonochromaticAbsorptionModel;
import edu.stanford.rsl.conrad.physics.detector.XRayDetector;
import edu.stanford.rsl.conrad.physics.materials.Material;
import edu.stanford.rsl.conrad.physics.materials.utils.AttenuationType;
import edu.stanford.rsl.conrad.rendering.AbstractRayTracer;
import edu.stanford.rsl.conrad.rendering.AbstractScene;
import edu.stanford.rsl.conrad.rendering.PrioritizableScene;
import edu.stanford.rsl.conrad.rendering.Priority1DRayTracer;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.FileUtil;
import edu.stanford.rsl.conrad.utils.RegKeys;
import edu.stanford.rsl.conrad.utils.UserUtil;

public class SurfaceBSplineVolumePhantom extends SliceWorker {

	private ArrayList<SurfaceBSpline> list;
	private AbstractScene scene;
	private PointND [] hullPoints;
	private ArrayList<Iterator<VectorPoint3D>> voxelList;
	private SimpleVector bounds;
	private boolean preserveAspects = true;
	private ArrayList<SimpleVector> localBounds;
	private boolean showRaster = false;
	private boolean showVertices = false;
	private boolean renderSolid = true;

	protected double xrayEnergy = 80;
	protected AttenuationType attType = AttenuationType.TOTAL_WITH_COHERENT_ATTENUATION;
	boolean renderAttenuation = false;


	@Override
	public void workOnSlice(int sliceNumber) {
		Trajectory geom = Configuration.getGlobalConfiguration().getGeometry();
		double originIndexX = geom.getOriginInPixelsX();
		double originIndexY = geom.getOriginInPixelsY();
		double originIndexZ = geom.getOriginInPixelsZ();
		double voxelSizeX = geom.getVoxelSpacingX();
		double voxelSizeY = geom.getVoxelSpacingY();
		double voxelSizeZ = geom.getVoxelSpacingZ();
		FloatProcessor slice = new FloatProcessor(geom.getReconDimensionX(), geom.getReconDimensionY());
		if (renderSolid){
			if (renderAttenuation && xrayEnergy > 0)
				slice.setValue(scene.getBackgroundMaterial().getAttenuation(xrayEnergy,attType));
			else
				slice.setValue(scene.getBackgroundMaterial().getDensity());
			slice.fill();
		}

		AbstractRayTracer tracer = new Priority1DRayTracer();
		tracer.setScene(scene);

		if (showVertices) {
			Iterator<VectorPoint3D> iterator = voxelList.get(sliceNumber);
			while(iterator.hasNext()){
				VectorPoint3D voxel = iterator.next();
				slice.putPixelValue((int)voxel.getX(), (int)voxel.getY(), voxel.getVector().getElement(0));
			}
		}

		double z  = (sliceNumber - originIndexZ) * voxelSizeZ;
		double xFirst = (-originIndexX) * voxelSizeX;
		double xLast = (slice.getWidth()-originIndexX) * voxelSizeX;
		if (renderSolid) {
			//long basetimeSlice = System.currentTimeMillis();
			for (int i = 0; i < slice.getHeight(); i ++){
				double y  = (i - originIndexY) * voxelSizeY;
				PointND pointLeft = new PointND(xFirst, y, z);
				PointND pointRight = new PointND(xLast, y, z);
				//long basetime = System.currentTimeMillis();
				StraightLine line = new StraightLine(pointLeft, pointRight);
				line.normalize();
				//System.out.println(line.getDirection());
				ArrayList<PhysicalObject> segments = tracer.castRay(line);
				//long castTime = System.currentTimeMillis() - basetime;
				//basetime = System.currentTimeMillis();
				if (segments != null) {
					for (PhysicalObject o: segments){
						Edge lineSegment = (Edge) o.getShape(); 
						PointND p1 = lineSegment.getPoint();
						PointND p2 = lineSegment.getEnd();
						int ix1 = (int) Math.round((p1.get(0) / voxelSizeX) + originIndexX);
						int ix2 = (int) Math.round((p2.get(0) / voxelSizeX) + originIndexX);
						int iy = (int) Math.round((y / voxelSizeY) + originIndexY);
						if (renderAttenuation && xrayEnergy > 0)
							slice.setValue(o.getMaterial().getAttenuation(xrayEnergy, attType)*o.getMaterial().getDensity());
						else
							slice.setValue(o.getMaterial().getDensity());
						slice.drawLine(ix1, iy, ix2, iy);
					}
					//long renderTime = System.currentTimeMillis() - basetime;
					//if (sliceNumber == 50) System.out.println("Cast time " + castTime + " render time " + renderTime + " " + segments.size());
				}

			}
			//long slicetime = System.currentTimeMillis() - basetimeSlice;
			//if (sliceNumber == 123) System.out.println("Slice time: " + slicetime);

		}
		if (showRaster) {
			for (double splineNum = 0; splineNum < list.size(); splineNum++) {
				for  (PointND p: hullPoints){
					if (Math.round(((p.get(2) / voxelSizeZ) + originIndexZ)) == sliceNumber){
						slice.putPixelValue((int)((p.get(0) / voxelSizeX) + originIndexX), (int)((p.get(1) / voxelSizeY) + originIndexY), splineNum+1);
					}
				}
			}
		}
		Grid2D grid = new Grid2D((float[]) slice.getPixels(), slice.getWidth(), slice.getHeight());
		imageBuffer.add(grid, sliceNumber);
	}

	@Override
	public String getProcessName() {
		return "Surface BSpline Phantom";
	}

	@Override
	public String getBibtexCitation() {
		return "see medline";
	}

	@Override
	public String getMedlineCitation() {
		return "L. Piegl and W. Tiller. The NURBS Book. Second Edition. Springer, Berlin, Heidelberg, New York. 1997.";
	}

	private void initBounds(){
		localBounds = new ArrayList<SimpleVector>();
		for (int i =0 ; i < list.size(); i++){
			SurfaceBSpline spline = list.get(i);
			SimpleVector bounds = SimpleOperators.concatenateVertically(spline.getMin().getAbstractVector(), spline.getMax().getAbstractVector());
			//System.out.println(bounds);
			localBounds.add(bounds);
		}
		SimpleVector max = SimpleOperators.max((SimpleVector []) localBounds.toArray(new SimpleVector [localBounds.size()])).getSubVec(3, 3);
		SimpleVector min = SimpleOperators.min((SimpleVector []) localBounds.toArray(new SimpleVector [localBounds.size()])).getSubVec(0, 3);

		bounds = SimpleOperators.concatenateVertically(min, max);		
	}

	public void readSplineListFromFile(String filename) throws IOException{
		list = SurfaceBSpline.readSplinesFromFile(filename);
		initBounds();
	}

	public void setSplineList(ArrayList<SurfaceBSpline> list){
		this.list = list;
		initBounds();
	}

	public void resizeVolumeToMatchBounds(PointND min, PointND max){
		bounds = SimpleOperators.concatenateVertically(min.getAbstractVector(), max.getAbstractVector());
		double rangeX = (max.getAbstractVector().getElement(0) - min.getAbstractVector().getElement(0))*1.05;
		double rangeY = (max.getAbstractVector().getElement(1) - min.getAbstractVector().getElement(1))*1.05;
		double rangeZ = (max.getAbstractVector().getElement(2) - min.getAbstractVector().getElement(2))*1.05;
		//System.out.println("Range Z: " + rangeZ);
		Trajectory traj = Configuration.getGlobalConfiguration().getGeometry();
		double voxelSizeX = rangeX / ((double) traj.getReconDimensionX());
		double voxelSizeY = rangeY / ((double) traj.getReconDimensionY());
		double voxelSizeZ = rangeZ / ((double) traj.getReconDimensionZ());
		double shiftX = 0;
		double shiftY = 0;
		double shiftZ = 0;
		if (preserveAspects) {
			double maxRatio = voxelSizeX;
			if (maxRatio < voxelSizeY) maxRatio = voxelSizeY;
			if (maxRatio < voxelSizeZ) maxRatio = voxelSizeZ;
			// size of object in voxel space
			double sizeX = rangeX / maxRatio;
			double sizeY = rangeY / maxRatio;
			double sizeZ = rangeZ / maxRatio;
			voxelSizeX = maxRatio;
			voxelSizeY = maxRatio;
			voxelSizeZ = maxRatio;
			shiftX = (traj.getReconDimensionX() - sizeX) / 2;
			shiftY = (traj.getReconDimensionY() - sizeY) / 2;
			shiftZ = (traj.getReconDimensionZ() - sizeZ) / 2;
		}
		traj.setOriginInPixelsX(shiftX-(min.getAbstractVector().getElement(0) / voxelSizeX));
		traj.setOriginInPixelsY(shiftY-(min.getAbstractVector().getElement(1) / voxelSizeY));
		traj.setOriginInPixelsZ(shiftZ-(min.getAbstractVector().getElement(2) / voxelSizeZ));
		traj.setVoxelSpacingX(voxelSizeZ);
		traj.setVoxelSpacingY(voxelSizeY);
		traj.setVoxelSpacingZ(voxelSizeX);
	}

	public void setBoundsFromGeometry(AnalyticPhantom4D phantom4d){
		Trajectory traj = Configuration.getGlobalConfiguration().getGeometry();
		SimpleVector minVec = new SimpleVector(traj.getOriginX(), traj.getOriginY(), traj.getOriginZ());
		String key = Configuration.getGlobalConfiguration().getRegistry().get(RegKeys.RENDER_PHANTOM_VOLUME_AUTO_CENTER);
		boolean set= false;
		if (key != null){
			if (key.equals("true")) {
				traj.setOriginToPhantomCenter(phantom4d);
				set=true;	
			}
		} 
		if (!set) {
			SimpleVector max = new SimpleVector(traj.getVoxelSpacingX()*traj.getReconDimensionX(),traj.getVoxelSpacingY()*traj.getReconDimensionY(), traj.getVoxelSpacingZ()*traj.getReconDimensionZ());
			max.add(minVec);
			bounds = SimpleOperators.concatenateVertically(minVec, max);
		}
	}

	public void resizeVolumeToMatchSplineSpace(){
		double rangeX = (bounds.getElement(3) - bounds.getElement(0))*1.05;
		double rangeY = (bounds.getElement(4) - bounds.getElement(1))*1.05;
		double rangeZ = (bounds.getElement(5) - bounds.getElement(2))*1.05;
		//System.out.println("Range Z: " + rangeZ);
		Trajectory traj = Configuration.getGlobalConfiguration().getGeometry();
		double voxelSizeX = rangeX / ((double) traj.getReconDimensionX());
		double voxelSizeY = rangeY / ((double) traj.getReconDimensionY());
		double voxelSizeZ = rangeZ / ((double) traj.getReconDimensionZ());
		double shiftX = 0;
		double shiftY = 0;
		double shiftZ = 0;
		if (preserveAspects) {
			double maxRatio = voxelSizeX;
			if (maxRatio < voxelSizeY) maxRatio = voxelSizeY;
			if (maxRatio < voxelSizeZ) maxRatio = voxelSizeZ;
			// size of object in voxel space
			double sizeX = rangeX / maxRatio;
			double sizeY = rangeY / maxRatio;
			double sizeZ = rangeZ / maxRatio;
			voxelSizeX = maxRatio;
			voxelSizeY = maxRatio;
			voxelSizeZ = maxRatio;
			shiftX = (traj.getReconDimensionX() - sizeX) / 2;
			shiftY = (traj.getReconDimensionY() - sizeY) / 2;
			shiftZ = (traj.getReconDimensionZ() - sizeZ) / 2;
		}
		traj.setOriginInPixelsX(shiftX - (bounds.getElement(0) / voxelSizeX));
		traj.setOriginInPixelsY(shiftY - (bounds.getElement(1) / voxelSizeY));
		traj.setOriginInPixelsZ(shiftZ - (bounds.getElement(2) / voxelSizeZ));
		traj.setVoxelSpacingX(voxelSizeZ);
		traj.setVoxelSpacingY(voxelSizeY);
		traj.setVoxelSpacingZ(voxelSizeX);
	}

	public ArrayList<AbstractShape> tesselateSplines(double samplingU, double samplingV){
		Trajectory traj = Configuration.getGlobalConfiguration().getGeometry();
		ArrayList<AbstractShape> meshList = new ArrayList<AbstractShape>();
		int numSlices = traj.getReconDimensionZ();
		voxelList = new ArrayList<Iterator<VectorPoint3D>>();
		ArrayList<ArrayList<VectorPoint3D>> pointList = new ArrayList<ArrayList<VectorPoint3D>>();
		for (int i = 0; i < numSlices; i++){
			pointList.add(new ArrayList<VectorPoint3D>());
		}	
		double voxelSizeX = traj.getVoxelSpacingX();
		double voxelSizeY = traj.getVoxelSpacingY();
		double voxelSizeZ = traj.getVoxelSpacingZ();
		double originIndexX = traj.getOriginInPixelsX();
		double originIndexY = traj.getOriginInPixelsY();
		double originIndexZ = traj.getOriginInPixelsZ();
		SimpleVector scaling = new SimpleVector(1.0 / voxelSizeX, 1.0 / voxelSizeY, 1.0 / voxelSizeZ);
		SimpleVector shift = new SimpleVector(originIndexX, originIndexY, originIndexZ);
		// Voxelize surface
		double samplingFactorU = samplingU / traj.getReconDimensionX();
		double samplingFactorV = samplingV / traj.getReconDimensionY();
		for (double splineNum = 0; splineNum < list.size(); splineNum++) {
			SimpleVector bound = localBounds.get((int) splineNum);
			double rangeX = (bound.getElement(3) - bound.getElement(0)) / voxelSizeX;
			double rangeY = (bound.getElement(4) - bound.getElement(1)) / voxelSizeY;
			double rangeZ = (bound.getElement(5) - bound.getElement(2)) / voxelSizeZ;
			int maxRange = (int) Math.ceil(Math.max(Math.max(rangeX, rangeY), rangeZ));
			samplingU = ((int)(maxRange * samplingFactorU));
			samplingV = ((int)(maxRange * samplingFactorV));
			// Old rendering (Sampline of the surface at a fixed grid):

			if (showStatus) IJ.showStatus("Sampling Shape " + list.get((int) splineNum).getTitle() + " (dim "+ (int)rangeX + "x"+ (int)rangeY + "x" + (int)rangeZ +") with " + samplingU + "x" +samplingV + " grid");
			//System.out.println("Sampling Shape " + list.get((int) splineNum).getTitle());
			double u = splineNum/list.size();
			if (showStatus) IJ.showProgress(u);

			if (showVertices) {
				for(double i =0 ; i < samplingU; i++){
					for (double j = 0; j < samplingV; j++){
						PointND p = list.get((int)splineNum).evaluate(i/samplingU, j/ samplingV);
						p.getAbstractVector().multiplyElementWiseBy(scaling);
						p.getAbstractVector().add(shift);
						int zPosition = (int) p.getAbstractVector().getElement(2);
						if (zPosition >= 0 && zPosition < numSlices){
							ArrayList<VectorPoint3D> list = pointList.get(zPosition);
							list.add(new VectorPoint3D(p, splineNum+1));
						}
					}
				}
			}
			AbstractShape mesh = list.get((int)splineNum).tessellateMesh(samplingU, samplingV);
			//System.out.println("Building done");
			meshList.add(mesh);
		}
		if (showStatus) IJ.showProgress(1.0);
		// prepare iterators for rendering.
		for (int i = 0; i < numSlices; i++){
			voxelList.add(Collections.synchronizedList(pointList.get(i)).iterator());
		}
		return meshList;
	}

	public void setScene(AbstractScene scene){
		this.scene = scene;
	}


	public void generateDefaultScene(double samplingU, double samplingV) throws IOException{
		scene = new PrioritizableScene();
		ArrayList<AbstractShape> meshList = tesselateSplines(samplingU, samplingV);
		for (int i = 0; i<meshList.size(); i++){
			AbstractShape mesh = meshList.get(i);
			PhysicalObject obj = new PhysicalObject();
			obj.setNameString(list.get((int)i).getTitle());
			obj.setMaterial(new Material(i+1));
			obj.setShape(mesh);
			((PrioritizableScene)scene).add(obj, PrioritizableScene.ADD_LOWEST_PRIORITY);
			if (showRaster) {
				hullPoints = mesh.getRasterPoints(200000);
			}
		}
	}

	public void configure() throws Exception{
		Trajectory traj = Configuration.getGlobalConfiguration().getGeometry();
		renderAttenuation = UserUtil.queryBoolean("Render energy dependent attenuation?");
		if(renderAttenuation)
			xrayEnergy = UserUtil.queryDouble("Monochromatic Xray energy [keV]", xrayEnergy);
		GenericDialog gd = new GenericDialog("Configure Surface BSpline Volume");
		int width = traj.getReconDimensionX();
		int height = traj.getReconDimensionY();
		gd.addSlider("Sampling in u direction", 10, width, width/2);
		gd.addSlider("Sampling in v direction", 10, height, height/2);
		gd.showDialog();
		if (gd.wasCanceled()){
			throw new RuntimeException("User cancelled selection.");
		}
		double samplingU = gd.getNextNumber();
		double samplingV = gd.getNextNumber();
		String filename = FileUtil.myFileChoose(".nrb", false);
		readSplineListFromFile(filename);
		generateDefaultScene(samplingU, samplingV);
		super.configure();
	}

	@Override
	public SliceWorker clone() {
		SurfaceBSplineVolumePhantom clone = new SurfaceBSplineVolumePhantom();
		clone.list = list;
		clone.voxelList = voxelList;
		clone.scene = scene;
		clone.hullPoints = hullPoints;
		clone.showStatus = showStatus;
		clone.renderAttenuation = renderAttenuation;
		clone.xrayEnergy = xrayEnergy;
		return clone;
	}

	/**
	 * Returns an SimpleVector that specifies the bounding box of the BSpline Phantom with six entries:<BR>
	 * <li>Minimum X Coordinate</li>
	 * <li>Minimum Y Coordinate</li>
	 * <li>Minimum Z Coordinate</li>
	 * <li>Maximum X Coordinate</li>
	 * <li>Maximum Y Coordinate</li>
	 * <li>Maximum Z Coordinate</li>
	 * @return the bounding box vector
	 */
	public SimpleVector getBounds(){
		return bounds;
	}

	
	public void setXrayEnergy(double xrayEnergy) {
		this.xrayEnergy = xrayEnergy;
	}
	
	public void setRenderAttenuation(boolean renderAttenuation) {
		this.renderAttenuation = renderAttenuation;
	}

}
/*
 * Copyright (C) 2010-2014 Andreas Maier
 * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
 */