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).
*/