Problem 990

Summary: Problems with G4PhantomParameterisation and large datasets
Product: Geant4 Reporter: marcus.h.mendenhall
Component: geometry/navigationAssignee: John Apostolakis <John.Apostolakis>
Status: RESOLVED FIXED    
Severity: normal CC: Gabriele.Cosmo, Pedro.Arce, robert.a.weller
Priority: P5    
Version: 9.1   
Hardware: All   
OS: All   

Description marcus.h.mendenhall 2008-01-03 20:35:56 CET
I am attempting to adapt a code I have for gigavoxel-level datasets to the new G4PhantomParameterisation and G4RegularNavigator systems in Geant4.9.1.  There are some issues with very large datasets. The code compresses the datasets dynamically in memory, so the number of voxels can typically hugely exceed the amount of memory on the machine.

First, the code starting at line 134 of G4Region.cc which looks like:

 * G4VPhysicalVolume* daughterPVol = lv->GetDaughter(0);
 * if (daughterPVol->IsParameterised())
 * {
 *   // Adopt special treatment in case of parameterised volumes,
 *   // where parameterisation involves a new material scan
 *   //
 *   G4VPVParameterisation* pParam = daughterPVol->GetParameterisation();
 *
 *   if (pParam->IsNested())
 *   {
 *     size_t matNo = pParam->GetMaterialScanner()->GetNumberOfMaterials();
 *     for (register size_t mat=0; mat<matNo; mat++)
 *    {
 *      volMat = pParam->GetMaterialScanner()->GetMaterial(mat);

does not use the the parameterisation's MaterialScanner unless the parameterisation is declared as nested.  The G4PhantomParameterisation is not nested, but if you have of order 10^9 voxels, scanning the slow way to fund the materials is painful.  I suggest decoupling this from the IsNested flag, and instead doing 

 *   if (pParam->GetMaterialScanner())
 *   {
 *     size_t matNo = pParam->GetMaterialScanner()->GetNumberOfMaterials();
 *     for (register size_t mat=0; mat<matNo; mat++)
 *    {
 *      volMat = pParam->GetMaterialScanner()->GetMaterial(mat);

so anyone who returns a valid material scanner pointer can use a potentially much more efficient method.  In my case, for example, I keep track of the materials as they are added, and can directly return the list of unique materials with no computational overhead. For now, I am declaring my un-nested parameterisation as nested, with no ill effect, but this makes for a bad coding style.

Next, the geometry optimizer does not recognize the fact that this parameterisation has its 
	param_pv->SetRegularStructureId(1); 
set to true, so it attempts to optimize 10^9 voxels, assuming it knows nothing about their organization.  This results in it allocating a vector of data with as many elements as the parameterisation, and therefore an assured out-of-memory condition. The geometry manager should completely bail out if the  structure is marked as regular, since it doesn't need to know anything else about what is inside, really.  

If these issues could be addressed, it would really open up the ability of Geant4 to work with full-resolution medical voxel volumes. 

Thanks.

Marcus Mendenhall
Comment 1 marcus.h.mendenhall 2008-01-03 23:44:09 CET
I have looked at the G4GeometryManager and would like to propose the following modification to fix this (if it has terrible, unintended side effects, I will be happy to entertain other methods).

This snippet is a modified section of G4GeometryManager.cc starting at line 170:

     volume->SetVoxelHeader(0);
     if (    (volume->IsToOptimise())
          && (volume->GetNoDaughters()>=kMinVoxelVolumesLevel1&&allOpts)
          || ( (volume->GetNoDaughters()==1)
            && (volume->GetDaughter(0)->IsReplicated()==true) && 
			(volume->GetDaughter(0)->GetRegularStructureId()!=1) ) ) 
     {

Note the final term in the conditional is all I changed.

Marcus

Comment 2 marcus.h.mendenhall 2008-01-09 21:54:59 CET
I would like to add some additional ideas here.  I know that, to some extent, this material belongs in the 'user requirements' tracker, but it is so closely related to the current issue that I am combining them.

It seems to me that it would make sense to add a couple of virtual methods to the G4VPVParameterisation base class:  
virtual G4bool hasMaterialBoundaryFinder(void) 

which returns a flag to indicate whether this parametrisation has a fast way to find the next material in a given direction, and 

virtual G4double DistanceToNextMaterial(G4ThreeVector &p, G4ThreeVector &v)

which returns the distance to the next material boundary starting from point p and going in direction v.  Then, the G4Navigator class would check the hasMaterialBoundaryFinder flag, and if it is set true, it will call DistanceToNextMaterial, instead of using its own boundary finder.  In the default G4VPVParameterisation, this flag would return false, allowing the the Navigator to use the current default behavior.

This would allow the voxel array to search efficiently, using its own internal information, for the next edge.  In particular, I have a voxel array which allows slightly fuzzy boundaries.  If the path clips the corner of a voxel without penetrating very far into it, it is OK to ignore the material transition.  This saves much of the delicate computational geometry associated with finding exact voxel boundaries.  Mine also interacts with the compression algorithm efficiently to look up the material in a huge voxel array appropriately.

This change would make the phantom navigation system much more extensible to the users' specific needs, with only a couple lines change in the G4 core routines.

Comment 3 John Apostolakis 2008-01-18 19:41:41 CET
Dear Marcus,

Thank you for your quick feedback to this 

I believe that I understand the essense of your feedback.  My first assessment is that the first ones are indeed corrections of our revised design and implementation.
  
In particular enabling users to create a material scanner is a simple correction or significant improvement, depending on the viewpoint.

If I have understood correctly, attempting to optimise volumes which respond have 
param_pv->SetRegularStructureId(1)
is simply a bug.  

I agree also that some aspects are primarily new user requirements. (Comment #3)
I have not yet had the opportunity to consider the proposed revision of the GeometryManager (Comment #2).  I will comment further on these in the future.

 
Comment 4 Gabriele Cosmo 2008-05-16 15:53:40 CEST
Issues #1 and #2 have been corrected and suggested fixes are included in tag "geommng-V09-01-04"
to be considered for the next release. Thanks!

Request #3 to be transferred to the URD tracker.