NumericalSheppLogan3D.java

/**
 * 
 */
package edu.stanford.rsl.conrad.phantom;

import edu.stanford.rsl.conrad.data.numeric.Grid3D;

/**
 * @author febmeng
 *
 */
public class NumericalSheppLogan3D {

	boolean init = false;
	int nx = 0;
	int ny = 0;
	int nz = 0;
	double wx, wy, wz;
	double dx, dy, dz;
	Grid3D image;
	SheppLogan3D phan;

	public NumericalSheppLogan3D( int maxI, int maxJ, int maxK ){
		
		initSheppLogan3D(  maxI, maxJ, maxK);
		
		if (nx == 0 || ny == 0 || nz == 0 ){
			System.out.println("Errpr: Wrong volume image size!");			
		}
		this.image = new Grid3D(nx, ny, nz);
		
	}
	
	
	public NumericalSheppLogan3D( Grid3D grid ){
		this.image = grid;
		initSheppLogan3D( image.getSize()[0], image.getSize()[1], image.getSize()[2] );
	}
	
	
	void initSheppLogan3D( int maxI, int maxJ, int maxK){

		this.nx = maxI;
		this.ny = maxJ;
		this.nz = maxK;
		wx = ( nx - 1 ) / 2.0 ;
		wy = ( ny - 1 ) / 2.0 ;
		wz = ( nz - 1 ) / 2.0 ;
		dx = 2.0 / ( nx - 1 );
		dy = 2.0 / ( ny - 1 );
		dz = 2.0 / ( nz - 1 );
		
		phan = new SheppLogan3D( ellipsoidsModifid );
		init = true;
	}

	void computeNumericalPhantom( ){
		
		//System.out.print("nx = " + nx + ", ny = " + ny + ", nz = " + nz + " .\n");
		//System.out.print("dx = " + dx + ", dy = " + dy + ", dz = " + dz + " .\n");
		//System.out.print("wx = " + wx + ", wy = " + wy + ", wz = " + wz + " .\n");
		
		for( int ix = 0; ix < nx; ix++ ){
			
			double x = ( ix - wx ) * dx ;
			
			for( int iy = 0; iy < ny; iy++ ){
				
				double y = ( iy - wy ) * dy ;
				
				for( int iz = 0; iz < nz; iz++ ){
					
					double z = ( iz - wz ) * dz ;
					
					image.setAtIndex(ix, iy, iz, (float)phan.ImageDomainSignal(x, y, z));
					
				}
			}
		}
		
	}
	
	public Grid3D getNumericalSheppLoganPhantom(){
		if( init) 
			computeNumericalPhantom();
		return image;
	}
	
	
	private static double[][] ellipsoidsModifid =
	// { delta_x, delta_y, delta_z,        a,       b,       c,            phi,  theta,  psi,     rho }
	 {  {       0,       0,       0,     0.69,    0.92,    0.81,              0,      0,    0,     1 },
		{       0, -0.0184,       0,   0.6624,   0.874,    0.78,              0,      0,    0,   -0.8 },
		{    0.22,       0,       0,     0.11,    0.31,    0.22, -(Math.PI)/10.,      0,    0,   -0.2 },
		{   -0.22,       0,       0,     0.16,    0.41,    0.28,  (Math.PI)/10.,      0,    0,   -0.2 },
		{       0,    0.35,   -0.15,     0.21,    0.25,    0.41,              0,      0,    0,    0.1 },
		{       0,     0.1,    0.25,    0.046,   0.046,    0.05,              0,      0,    0,    0.1 },
		{       0,    -0.1,    0.25,    0.046,   0.046,    0.05,              0,      0,    0,    0.1 },
		{   -0.08,  -0.605,       0,    0.046,   0.023,    0.05,              0,      0,    0,    0.1 },
		{       0,  -0.605,       0,    0.023,   0.023,    0.02,              0,      0,    0,    0.1 },
		{    0.06,  -0.605,       0,    0.023,   0.046,    0.02,              0,      0,    0,    0.1 }};

	
	private static double[][] sphere =
	// { delta_x, delta_y, delta_z,        a,       b,       c,            phi,  theta,  psi,     rho }
	 {  {       0,       0,       0,     0.69,    0.69,    0.69,              0,      0,    0,     1 }};
	
	private static double[][] point =
	// { delta_x, delta_y, delta_z,        a,       b,       c,            phi,  theta,  psi,     rho }
	 {  {       0.6,      0.4,       0,     0.19,    0.19,    0.19,              0,      0,    0,     1 }};
	
	private static double[][] cylinder =
			// { delta_x, delta_y, delta_z,        a,       b,       c,            phi,  theta,  psi,     rho }
	{  {       0,      0,       0,     0.69,    0.69,      10,        0,      0,    0,     1 }};
			

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