Message: Re: Execution Speed and Parameterization Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Feedback Re: Execution Speed and Parameterization 

Forum: Medical Applications
Re: Question Execution Speed and Parameterization (Ken Sutherland)
Re: None Re: Execution Speed and Parameterization (Hongyu Jiang)
Re: Sad Re: Execution Speed and Parameterization
Date: 10 Dec, 2004
From: Ken Sutherland <ken@cme.hokudai.ac.jp>

Dr. Sasaki,

This is what I have been able to come up following Dr. Jiang's paper. 
Note that this solution is completely untested, I've only got it 
working today! If anyone can point out any problems, please do.

First, get rid of PatientParameterisation::patientPlacementX, Y, and Z.
Storing the position of each voxel requires 512*512*8 = 2 MB of memory
per slice (in the worst case). The x,y,z position can easily be 
calculated in ComputeTransformation based on copyNo. Here is what mine 
looks like:

void PatientParameterisation::ComputeTransformation(const G4int copyNo, G4VPhysicalVolume* physVol) const
{
//	G4double originX = patientPlacementX[copyNo];
//	G4double originY = patientPlacementY[copyNo];
//	G4double originZ = patientPlacementZ[copyNo]*mm-middleLocationValue*mm-sliceThickness/2;

	// Compute the voxel position from the copy nummber.
	G4int nVoxelX = copyNo % lenColumns;
	G4int nVoxelY = (copyNo / lenColumns) % lenRows;
	G4int nVoxelZ = copyNo / (lenRows * lenColumns);

	// Compute the world position.
	G4double originX = compression*(nVoxelX - lenColumns/2)*pixelSpacingX;
	G4double originY = compression*(lenRows/2 - nVoxelY)*pixelSpacingY;
	G4double originZ = (((G4double)nVoxelZ) - (totalNumberOfFile - 1)/2.0)*sliceThickness;

	G4ThreeVector origin( 
		originX, 
		originY, 
		originZ );

	physVol->SetTranslation(origin);

You will also have to modify ComputeMaterial to return "Air" if the
density is less than 0.207 because the original DICOM sample does
not store voxels for air. 

Next, I modified G4VPhysicalVolume.hh to hold the six adjacent voxel
numbers.

    public:
	void SetVoxelLeft(G4int n)	{ nAdjacentVoxel[0] = n;}
	void SetVoxelRight(G4int n)	{ nAdjacentVoxel[1] = n;}
	void SetVoxelUp(G4int n)	{ nAdjacentVoxel[2] = n;}
	void SetVoxelDown(G4int n)	{ nAdjacentVoxel[3] = n;}
	void SetVoxelFront(G4int n)	{ nAdjacentVoxel[4] = n;}
	void SetVoxelBack(G4int n)	{ nAdjacentVoxel[5] = n;}

	G4int GetVoxelLeft()		{ return nAdjacentVoxel[0];}
	G4int GetVoxelRight()		{ return nAdjacentVoxel[1];}
	G4int GetVoxelUp()		{ return nAdjacentVoxel[2];}
	G4int GetVoxelDown()		{ return nAdjacentVoxel[3];}
	G4int GetVoxelFront()		{ return nAdjacentVoxel[4];}
	G4int GetVoxelBack()		{ return nAdjacentVoxel[5];}

	G4int* GetAdjacentVoxel() { return nAdjacentVoxel;}

    private:
	G4int nAdjacentVoxel[6];

And I added code to ComputeTransformation to compute the adjacent
voxel numbers:

	if (nVoxelX > 0)
		physVol->SetVoxelLeft(copyNo - 1);
	else
		physVol->SetVoxelLeft(-1);

	if (nVoxelX < lenColumns - 1)
		physVol->SetVoxelRight(copyNo + 1);
	else
		physVol->SetVoxelRight(-1);

	if (nVoxelY > 0)
		physVol->SetVoxelUp(copyNo - lenColumns);
	else
		physVol->SetVoxelUp(-1);

	if (nVoxelY < lenRows - 1)
		physVol->SetVoxelDown(copyNo + lenColumns);
	else
		physVol->SetVoxelDown(-1);

	if (nVoxelZ > 0)
		physVol->SetVoxelBack(copyNo - (lenColumns * lenRows));
	else
		physVol->SetVoxelBack(-1);

	if (nVoxelZ < totalNumberOfFile - 1)
		physVol->SetVoxelFront(copyNo + (lenColumns * lenRows));
	else
		physVol->SetVoxelFront(-1);

Now the hard part. I derived a new class from G4ParameterisedNavigation 
and overloaded the function LevelLocate.

class G4ParameterisedNavigationOpt : public G4ParameterisedNavigation
{
  public:  

    G4bool LevelLocate( G4NavigationHistory& history,
                  const G4VPhysicalVolume* blockedVol,
                  const G4int blockedNum,
                  const G4ThreeVector& globalPoint,
                  const G4ThreeVector* globalDirection,
                  const G4bool pLocatedOnEdge, 
                        G4ThreeVector& localPoint );
};

I modified G4Navigator.hh to use my new class.

//#include "G4ParameterisedNavigation.hh"	// Original
#include "G4ParameterisedNavigationOpt.hh"	// Optimized

//G4ParameterisedNavigation fparamNav;		// Original
  G4ParameterisedNavigationOpt fparamNav;	// Optimized

I copied the original LevelLocate function from 
G4ParameterisedNavigation.icc and started messing with it. I ended up
with this:

G4bool
G4ParameterisedNavigationOpt::LevelLocate( G4NavigationHistory& history,
                                  const G4VPhysicalVolume* blockedVol,
                                  const G4int blockedNum,
                                  const G4ThreeVector& globalPoint,
                                  const G4ThreeVector* globalDirection,
                                  const G4bool pLocatedOnEdge, 
                                        G4ThreeVector& localPoint )
{
	G4SmartVoxelHeader *motherVoxelHeader;
	G4SmartVoxelNode *motherVoxelNode;
	G4VPhysicalVolume *motherPhysical, *pPhysical;
	G4VPVParameterisationOpt *pParam;
	G4LogicalVolume *motherLogical;
	G4VSolid *pSolid;
	G4ThreeVector samplePoint;
	G4int voxelNoDaughters, replicaNo;

	motherPhysical = history.GetTopVolume();
	motherLogical = motherPhysical->GetLogicalVolume();
	motherVoxelHeader = motherLogical->GetVoxelHeader();

	// Find the voxel containing the point
	//
	motherVoxelNode = ParamVoxelLocate(motherVoxelHeader,localPoint);
	  
	voxelNoDaughters = 1;		  
	pPhysical = motherLogical->GetDaughter(0);
	pParam = (G4VPVParameterisationOpt*)pPhysical->GetParameterisation();

	bool bSearchAdjacent = false;
	G4int pAdjacentVoxel[6];

	if (blockedNum >= 0) {
		G4EventManager* pEventManager = G4EventManager::GetEventManager();
		const G4TrackingManager* pTrackingManager = pEventManager->GetTrackingManager();
		const G4Track* pTrack = pTrackingManager->GetTrack();
		int nStepNumber = pTrack->GetCurrentStepNumber();
		if (nStepNumber > 0) {
			// Turn on the flag!
			bSearchAdjacent = true;
			voxelNoDaughters = 6;

			// NOTE: I'd rather get the adjacent info from blockedVol, but can't because it's "const".
			pParam->ComputeTransformation(blockedNum, pPhysical);
			pAdjacentVoxel[0] = pPhysical->GetVoxelLeft();
			pAdjacentVoxel[1] = pPhysical->GetVoxelRight();
			pAdjacentVoxel[2] = pPhysical->GetVoxelUp();
			pAdjacentVoxel[3] = pPhysical->GetVoxelDown();
			pAdjacentVoxel[4] = pPhysical->GetVoxelFront();
			pAdjacentVoxel[5] = pPhysical->GetVoxelBack();
		}
	}


	// Search replicated daughter volume
	//
	for ( register int sampleNo = 0; sampleNo < voxelNoDaughters; sampleNo++ )
	{
		if (bSearchAdjacent)
			replicaNo = pAdjacentVoxel[sampleNo];
		else
			//replicaNo = motherVoxelNode->GetVolume(sampleNo);
			replicaNo = pParam->ComputeCopyNo(globalPoint);

		if ( replicaNo != -1 )
		{
			// Obtain solid (as it can vary) and obtain its parameters
			//
			pSolid = pParam->ComputeSolid(replicaNo, pPhysical); 
			pSolid->ComputeDimensions(pParam, replicaNo, pPhysical);
			pParam->ComputeTransformation(replicaNo, pPhysical);

			// Setup history
			//
			history.NewLevel(pPhysical, kParameterised, replicaNo);
			samplePoint = history.GetTopTransform().TransformPoint(globalPoint);
			if ( bSearchAdjacent && !G4AuxiliaryNavServices::CheckPointOnSurface(pSolid,
										samplePoint, globalDirection, 
										history.GetTopTransform(), pLocatedOnEdge) )
			{
				history.BackLevel();
			}
			else
			{
				// Enter this daughter
				//
				localPoint = samplePoint;
			    
				// Set the correct copy number in physical
				//
				pPhysical->SetCopyNo(replicaNo);
			    
				// Set the correct solid and material in Logical Volume
				//
				G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
				pLogical->SetSolid(pSolid);
				pLogical->SetMaterial(pParam->ComputeMaterial(replicaNo, pPhysical));

				return true;
			}
		}
	}

	//G4cerr << "Oops! G4ParameterisedNavigationOpt::LevelLocate failed!\n";
	return LevelLocateOriginal(
		history,
		blockedVol,
		blockedNum,
		globalPoint,
		globalDirection,
		pLocatedOnEdge, 
		localPoint );
}

The first time LevelLocate is called with blockedNum == -1. In this 
case, the original code always seemed to return the node at the center
of the volume. To handle this case, I derived a class from 
G4VPVParameterisation and added function called ComputeCopyNo.

class G4VPVParameterisationOpt : public G4VPVParameterisation
{
	public:

	    virtual G4int ComputeCopyNo(const G4ThreeVector& globalPoint) const = 0;

};

Then my implementation of this function is:

G4int PatientParameterisation::ComputeCopyNo(const G4ThreeVector& point) const
{
	G4int x = (int)((columns/2. + (point.x() / pixelSpacingX))/compression);
	G4int y = (int)((rows/2. - (point.y() / pixelSpacingY))/compression);
	G4int z = (int)(totalNumberOfFile/2. + (point.z() / sliceThickness));

	if ((x < 0) || (x > lenColumns-1) || (y < 0) || (y > lenRows-1) || (z < 0) || (z > totalNumberOfFile-1))
		return -1;

	int nCopyNo = (z*lenColumns*lenRows) + (y*lenColumns) + x;
	return nCopyNo;
}

This function will used again during the second call to LevelLocate,
when step number == 1. This time globalPoint will contain the initial 
point.

Finally, if all else fails, I fall back on a copy of the original 
function. This happens when a particle enters a diagonal voxel. Or it
may happen when a particle "skips over" a voxel, if that is possible.
Anyway, the original function seems to be called only very rarely. I'm
not sure if adding the code to handle diagonal voxels is worth the 
overhead. You would have to consider at least 8 additional adjacent
voxels. I tried it, but gave up.

Anyway, I hope that this helps anyone out there get started with 
this.

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

 Add Message Add Message
to: "Re: Execution Speed and Parameterization"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]