package edu.stanford.rsl.conrad.phantom.xcat; import java.io.IOException; import java.util.ArrayList; import java.util.Collections; import java.util.Iterator; import edu.stanford.rsl.conrad.geometry.AbstractShape; import edu.stanford.rsl.conrad.geometry.General; import edu.stanford.rsl.conrad.geometry.Rotations; import edu.stanford.rsl.conrad.geometry.motion.MotionField; import edu.stanford.rsl.conrad.geometry.motion.PointBasedMotionField; import edu.stanford.rsl.conrad.geometry.motion.timewarp.HarmonicTimeWarper; import edu.stanford.rsl.conrad.geometry.motion.timewarp.TimeWarper; import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND; import edu.stanford.rsl.conrad.geometry.splines.BSpline; import edu.stanford.rsl.conrad.geometry.splines.SurfaceBSpline; import edu.stanford.rsl.conrad.geometry.splines.TimeVariantSurfaceBSpline; import edu.stanford.rsl.conrad.geometry.transforms.ScaleRotate; import edu.stanford.rsl.conrad.geometry.transforms.Transform; import edu.stanford.rsl.conrad.geometry.transforms.Transformable; import edu.stanford.rsl.conrad.geometry.transforms.Translation; import edu.stanford.rsl.conrad.numerics.SimpleMatrix; import edu.stanford.rsl.conrad.numerics.SimpleOperators; import edu.stanford.rsl.conrad.numerics.SimpleVector; import edu.stanford.rsl.conrad.parallel.ParallelThreadExecutor; 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.DoubleArrayUtil; import edu.stanford.rsl.conrad.utils.RegKeys; import edu.stanford.rsl.jpop.utils.UserUtil; public class CoronaryScene extends XCatScene { /** * */ private static final long serialVersionUID = -8901686267381868405L; protected ArrayList<SurfaceBSpline> arteries; private boolean renderHeartBody = false; private MotionField heartMotionField; ScaleRotate heartRotation; public CoronaryScene(){ } @Override public void configure() throws Exception{ double heartCycles = UserUtil.queryDouble("Number of heart beats in Scene (acquisition time * heart rate):", 4.0); Configuration config = Configuration.getGlobalConfiguration(); int dimz = config.getGeometry().getProjectionStackSize() * config.getNumSweeps(); double [] heartStates = new double [dimz]; for (int i = 0; i < dimz; i++){ heartStates[i]= (i*heartCycles) / (dimz); } config.getRegistry().put(RegKeys.HEART_PHASES, DoubleArrayUtil.toString(heartStates)); warper = new HarmonicTimeWarper(heartCycles); init(); createArteryTree(this); createPhysicalObjects(); } public ArrayList<SurfaceBSpline> readHeartState(int state){ try { return SurfaceBSpline.readSplinesFromFile(XCatDirectory + "/heart_base" + state + ".nrb"); } catch (IOException e) { // TODO Auto-generated catch block e.printStackTrace(); return null; } } /** * * @param referenceMotion */ public void createArteryTree(MotionField referenceMotion){ try { arteries = SurfaceBSpline.readSplinesFromFile(XCatDirectory + "/artery_tree.nrb"); for (SurfaceBSpline spline:arteries){ applyTransform(spline, heartRotation); System.out.println("Converting spline " + spline.getTitle() + " to time-variant mode."); TimeVariantSurfaceBSpline tspline = new TimeVariantSurfaceBSpline(spline, referenceMotion, 27, true); variants.add(tspline); } } catch (IOException e) { // TODO Auto-generated catch block e.printStackTrace(); } } public void init(){ setName("Heart Scene"); SimpleMatrix rotation = SimpleOperators.multiplyMatrixProd(SimpleOperators.multiplyMatrixProd( Rotations.createBasicZRotationMatrix(General.toRadians(165)), SimpleOperators.multiplyMatrixProd( Rotations.createBasicXRotationMatrix(General.toRadians(45)), Rotations.createBasicYRotationMatrix(General.toRadians(45)))), Rotations.createBasicZRotationMatrix(General.toRadians(50))); rotation = SimpleOperators.multiplyMatrixProd( SimpleOperators.multiplyMatrixProd(SimpleOperators.multiplyMatrixProd( Rotations.createBasicZRotationMatrix(General.toRadians(165)), SimpleOperators.multiplyMatrixProd( Rotations.createBasicXRotationMatrix(General.toRadians(45)), Rotations.createBasicYRotationMatrix(General.toRadians(45)))), Rotations.createBasicZRotationMatrix(General.toRadians(50))), Rotations.createBasicZRotationMatrix(General.toRadians(180))); heartRotation = new ScaleRotate(rotation.multipliedBy(0.85)); String rotationString = Configuration.getGlobalConfiguration().getRegistryEntry(RegKeys.XCAT_HEART_ROTATION); if (rotationString != null){ String [] values = rotationString.split(", "); rotation = SimpleOperators.multiplyMatrixProd( Rotations.createBasicXRotationMatrix(General.toRadians(Double.parseDouble(values[0]))), SimpleOperators.multiplyMatrixProd( Rotations.createBasicYRotationMatrix(General.toRadians(Double.parseDouble(values[1]))), Rotations.createBasicZRotationMatrix(General.toRadians(Double.parseDouble(values[2]))))); heartRotation = new ScaleRotate(rotation.multipliedBy(Double.parseDouble(values[3]))); } ArrayList<ArrayList<SurfaceBSpline>> completeSplineList = new ArrayList<ArrayList<SurfaceBSpline>>(); int splineNum = 0; for (int k = 0; k < 27; k++){ ArrayList<SurfaceBSpline> list = readHeartState(k); splineNum = list.size(); completeSplineList.add(list); } max = new PointND(-Double.MAX_VALUE, -Double.MAX_VALUE, -Double.MAX_VALUE); min = new PointND(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE); for (int i=0; i<splineNum; i++){ ArrayList<SurfaceBSpline> surfaces = new ArrayList<SurfaceBSpline>(); for (int k = 0; k < 27; k++){ surfaces.add(completeSplineList.get(k).get(i)); } TimeVariantSurfaceBSpline variant = new TimeVariantSurfaceBSpline(surfaces); max.updateIfHigher(variant.getMax()); min.updateIfLower(variant.getMin()); variants.add(variant); } heartMotionField = new PointBasedMotionField(variants, 5); if (!renderHeartBody) { variants.clear(); } //applyTransform(heartRotation); } private void updateMinAndMax(){ max = new PointND(-Double.MAX_VALUE, -Double.MAX_VALUE, -Double.MAX_VALUE); min = new PointND(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE); for (int i=0; i <variants.size(); i++){ TimeVariantSurfaceBSpline variant = variants.get(i); max.updateIfHigher(variant.getMax()); min.updateIfLower(variant.getMin()); } } private void applyTransform(Transform t){ SimpleVector center = SimpleOperators.add(min.getAbstractVector(), max.getAbstractVector()).dividedBy(2); Translation backShift = new Translation(center.clone()); center.negate(); Translation centerShift = new Translation(center); for (int i=0; i <variants.size(); i++){ TimeVariantSurfaceBSpline variant = variants.get(i); variant.applyTransform(centerShift); variant.applyTransform(t); variant.applyTransform(backShift); } heartMotionField = new PointBasedMotionField(variants, 5); updateMinAndMax(); } private void applyTransform(Transformable shape, Transform t){ SimpleVector center = SimpleOperators.add(min.getAbstractVector(), max.getAbstractVector()).dividedBy(2); Translation backShift = new Translation(center.clone()); center.negate(); Translation centerShift = new Translation(center); shape.applyTransform(centerShift); shape.applyTransform(t); shape.applyTransform(backShift); } @Override public float[] getBinaryRepresentation() { int totalsize = 2; float [][] timeVariants = new float [variants.size()][]; for (int i = 0; i<variants.size();i++){ timeVariants[i] = variants.get(i).getBinaryRepresentation(); totalsize += timeVariants[i].length; } float [] binary = new float [totalsize + timeVariants.length + timeVariants.length]; binary[0] = BSpline.BSPLINECOLLECTION; binary[1] = totalsize; binary[2] = timeVariants.length; int index = 3; for (int i = 0; i<timeVariants.length;i++){ for (int j=0;j<timeVariants[i].length; j++){ binary[index]= timeVariants[i][j]; index ++; } } for (int i = 0; i<timeVariants.length;i++){ binary[index] = XCatScene.getSplinePriorityLUT().get(variants.get(i).getTitle()); } for (int i = 0; i<timeVariants.length;i++){ binary[index] = (float) XCatMaterialGenerator.generateFromMaterialName(XCatScene.getSplineNameMaterialNameLUT().get(variants.get(i).getTitle())).getDensity(); } return binary; } public PrioritizableScene tessellateScene(double time){ // remove previous state scene objects. PrioritizableScene scene = new PrioritizableScene(); clearObjectsOnly(); TessellationThread[] threads = new TessellationThread[CONRAD.getNumberOfThreads()]; Iterator<TimeVariantSurfaceBSpline> list = Collections.synchronizedList(variants).iterator(); for (int splineNum = 0; splineNum < CONRAD.getNumberOfThreads(); splineNum++) { threads[splineNum] = new TessellationThread(list, warper.warpTime(time)); } ParallelThreadExecutor executor = new ParallelThreadExecutor(threads); try { executor.execute(); } catch (InterruptedException e) { e.printStackTrace(); } for (int splineNum = 0; splineNum < CONRAD.getNumberOfThreads(); splineNum++) { for (AbstractShape s: threads[splineNum].getObjects()) { add(scene, s, s.getName()); } } return scene; } @Override public ArrayList<SurfaceBSpline> getSplines() { return readHeartState(0); } @Override public PointND getPosition(PointND initialPosition, double initialTime, double time) { return heartMotionField.getPosition(initialPosition, initialTime, time); } @Override public ArrayList<PointND> getPositions(PointND initialPosition, double initialTime, double... times) { return heartMotionField.getPositions(initialPosition, initialTime, times); } @Override public void setTimeWarper(TimeWarper warper){ super.setTimeWarper(warper); heartMotionField.setTimeWarper(warper); } @Override public String getName() { return "XCat Coronary"; } @Override public String getBibtexCitation() { return "@article{Segars08-RCS,\n" + " author={Segars WP, Mahesh M, Beck TJ, Frey EC, Tsui BMW},\n" + " title={Realistic CT simulation using the 4D XCAT phantom},\n" + " journal={Med. Phys.},\n" + " volume={35},\n" + " pages={3800-3808},\n" + " year={2008}\n" + "}"; } @Override public String getMedlineCitation() { return "Segars WP, Mahesh M, Beck TJ, Frey EC, Tsui BMW. Realistic CT simulation using the 4D XCAT phantom. Med. Phys. 35:3800-3808 (2008);"; } } /* * Copyright (C) 2010-2014 Andreas Maier * CONRAD is developed as an Open Source project under the GNU General Public License (GPL). */