ApplyDenseMotionFieldExample.java example

Explorer
CONRAD-master
- src
package edu.stanford.rsl.tutorial.motion;

import java.io.IOException;

import com.jogamp.opencl.CLContext;
import com.jogamp.opencl.CLDevice;

import ij.ImagePlus;
import edu.stanford.rsl.apps.gui.XCatMetricPhantomCreator;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
import edu.stanford.rsl.conrad.data.numeric.MultiChannelGrid2D;
import edu.stanford.rsl.conrad.geometry.motion.OpenCLParzenWindowMotionField;
import edu.stanford.rsl.conrad.geometry.motion.ParzenWindowMotionField;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.conrad.opencl.OpenCLUtil;
import edu.stanford.rsl.conrad.phantom.AnalyticPhantom4D;
import edu.stanford.rsl.conrad.utils.CONRAD;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.ImageUtil;

public class ApplyDenseMotionFieldExample {

	public static void main(String[] args) {
		int timesteps = 10; 
		CONRAD.setup();
		Configuration config = Configuration.getGlobalConfiguration();
		XCatMetricPhantomCreator creator = new XCatMetricPhantomCreator();
		creator.setSteps(timesteps);
		AnalyticPhantom4D scene = creator.instantiateScene();
		ImagePlus hyper = creator.renderMetricVolumePhantom(scene);
		Grid3D hyperGrid = ImageUtil.wrapImagePlus(hyper);
		hyperGrid.show("Sampled from Phantom");
		Grid3D deformed = new Grid3D(hyperGrid, false);
		float motionfield [] [] = new float [timesteps][];
		PointND splineSurfacePoints [] [] = new PointND [timesteps][];
		CLContext context = OpenCLUtil.createContext();
		CLDevice device = context.getMaxFlopsDevice();
		ParzenWindowMotionField motion = (ParzenWindowMotionField) scene.getMotionField();
		OpenCLParzenWindowMotionField cl;
		int factor = 16; 
		double originx = config.getGeometry().getOriginX();
		try {
			cl = new OpenCLParzenWindowMotionField(motion, context, device);
			for (int i = 0; i < timesteps; i++){
				double time = ((double)i)/(timesteps-1);
				motionfield[i] =cl.getMotionFieldAsArrayReduceZGridXY(time, 0, config.getGeometry().getReconDimensionX()/factor, config.getGeometry().getReconDimensionY()/factor); 	
				splineSurfacePoints [i] = cl.getRasterPoints(time);
			}
			
			for (int z = 0; z < deformed.getSize()[2]; z++){
				MultiChannelGrid2D slice = new MultiChannelGrid2D(deformed.getSize()[0], deformed.getSize()[1], timesteps);
				for (int i = 0; i < timesteps; i++){
					Grid2D slicet = new Grid2D(deformed.getSize()[0], deformed.getSize()[1]);
					double zCoord = edu.stanford.rsl.conrad.geometry.General.voxelToWorld(z, config.getGeometry().getVoxelSpacingZ(), config.getGeometry().getOriginZ());
					int idxz = z*(deformed.getSize()[0]*deformed.getSize()[1]);
					for (int y = 0; y < deformed.getSize()[1]; y++){
						double yCoord = edu.stanford.rsl.conrad.geometry.General.voxelToWorld(y, config.getGeometry().getVoxelSpacingY(), config.getGeometry().getOriginY());
						int idxy = y*(deformed.getSize()[0]);
						for (int x = 0; x < deformed.getSize()[0]; x++){
							double xCoord = edu.stanford.rsl.conrad.geometry.General.voxelToWorld(x, config.getGeometry().getVoxelSpacingX(), originx);
							int idx = idxz + idxy + x;
							double xCoordTemp = xCoord + (motionfield[i][(idx*3)]);
							double yCoordTemp = yCoord + (motionfield[i][(idx*3)+1]);
							double zCoordTemp = zCoord + (motionfield[i][(idx*3)+2]);
							double xCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(xCoordTemp, config.getGeometry().getVoxelSpacingX(), originx);
							double yCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(yCoordTemp, config.getGeometry().getVoxelSpacingY(), config.getGeometry().getOriginY());
							double zCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(zCoordTemp, config.getGeometry().getVoxelSpacingZ(), config.getGeometry().getOriginZ());
							float value = InterpolationOperators.interpolateLinear(hyperGrid, zCoordVoxel, xCoordVoxel, yCoordVoxel);
							slicet.setAtIndex(x, y, value);
						}
					}
					slice.setChannel(i, slicet);
				}
				deformed.setSubGrid(z, slice);
			}
			for (int i = 0; i < timesteps; i++){
				for (int j = 0; j < splineSurfacePoints[0].length; j++){
					PointND point = splineSurfacePoints [0][j].clone();
					PointND from = splineSurfacePoints [i][j].clone();
					
					double xCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(from.get(0), config.getGeometry().getVoxelSpacingX(), originx);
					double yCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(from.get(1), config.getGeometry().getVoxelSpacingY(), config.getGeometry().getOriginY());
					double zCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(from.get(2), config.getGeometry().getVoxelSpacingZ(), config.getGeometry().getOriginZ());
					int z = (int) Math.round(zCoordVoxel);
					int y = (int) Math.round(yCoordVoxel);
					int x = (int) Math.round(xCoordVoxel);
					int idx = z*(deformed.getSize()[0]*deformed.getSize()[1]) + y*(deformed.getSize()[0]) + x;
					
					double xCoordTemp = from.get(0) + (motionfield[i][(idx*3)]);
					double yCoordTemp = from.get(1) + (motionfield[i][(idx*3)+1]);
					double zCoordTemp = from.get(2) + (motionfield[i][(idx*3)+2]);
					PointND interp = new PointND(xCoordTemp, yCoordTemp, zCoordTemp);
					
					xCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(xCoordTemp, config.getGeometry().getVoxelSpacingX(), originx);
					yCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(yCoordTemp, config.getGeometry().getVoxelSpacingY(), config.getGeometry().getOriginY());
					zCoordVoxel = edu.stanford.rsl.conrad.geometry.General.worldToVoxel(zCoordTemp, config.getGeometry().getVoxelSpacingZ(), config.getGeometry().getOriginZ());
					//z = (int) Math.round(zCoordVoxel);
					//y = (int) Math.round(yCoordVoxel);
					//x = (int) Math.round(xCoordVoxel);
					
					MultiChannelGrid2D multi = (MultiChannelGrid2D) deformed.getSubGrid(z);
					if (multi !=null) {
						//multi.getChannel(i).setAtIndex(x, y, j);
					}
					point.getAbstractVector().subtract(from.getAbstractVector());
					interp.getAbstractVector().subtract(from.getAbstractVector());
					SimpleVector diff = interp.getAbstractVector().clone();
					diff.subtract(point.getAbstractVector());
					if (Math.abs(diff.normL2()) > 0.01) {
						//System.out.println(j+": error: " + diff.normL2() + " expected point:" + point + " interp:" +interp);
					}
				}
			}
			
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		deformed.show("Phantom created with deformation");
	}

}