Merge pull request #219 from djboersma/allow_some_HU_out_of_table

add option to allow out-of-bounds HU values
OpenGATE · Feb 13, 2019 · e80adc6 · e80adc6
2 parents eccc25a + 50adde4
commit e80adc6
Show file tree

Hide file tree

Showing 4 changed files with 75 additions and 12 deletions.
diff --git a/source/geometry/include/GateVImageVolume.hh b/source/geometry/include/GateVImageVolume.hh
@@ -165,6 +165,7 @@ public:
   void SetDensityImageFilename(G4String filename);
   void SetMassImageFilename   (G4String filename) {mMassImageFilename = filename;}
   void EnableBoundingBoxOnly(bool b);
+  void SetMaxOutOfRangeFraction(double f);
 
 protected:
 
@@ -250,6 +251,9 @@ protected:
 
   //-----------------------------------------------------------------------------
   bool mIsBoundingBoxOnlyModeEnabled;
+  unsigned int mUnderflow;
+  unsigned int mOverflow;
+  double mMaxOutOfRangeFraction;
 };
 // EO class GateVImageVolume
 //-----------------------------------------------------------------------------

diff --git a/source/geometry/include/GateVImageVolumeMessenger.hh b/source/geometry/include/GateVImageVolumeMessenger.hh
@@ -26,6 +26,7 @@ class GateVImageVolume;
 class G4UIcmdWithAString;
 class G4UIcmdWith3VectorAndUnit;
 class G4UIcmdWithABool;
+class G4UIcmdWithADouble;
 
 //-----------------------------------------------------------------------------
 /// \brief Messenger of GateVImageVolume
@@ -53,6 +54,7 @@ private:
   G4UIcmdWithAString        * pBuildDensityImageCmd;
   G4UIcmdWithAString        * pBuildMassImageCmd;
   G4UIcmdWithABool          * pDoNotBuildVoxelsCmd;
+  G4UIcmdWithADouble        * pSetMaxOutOfRangeFractionCmd;
 };
 //-----------------------------------------------------------------------------
 

diff --git a/source/geometry/src/GateVImageVolume.cc b/source/geometry/src/GateVImageVolume.cc
@@ -56,6 +56,9 @@ GateVImageVolume::GateVImageVolume( const G4String& name,G4bool acceptsChildren,
   mMassImageFilename    = "none";
   mIsBoundingBoxOnlyModeEnabled = false;
   mImageMaterialsFromHounsfieldTableDone = false;
+  mUnderflow = 0;
+  mOverflow = 0;
+  mMaxOutOfRangeFraction = 0.0;
   GateMessageDec("Volume",5,"End GateVImageVolume("<<name<<")\n");
 
   // do not display all voxels, only bounding box
@@ -84,6 +87,17 @@ void GateVImageVolume::EnableBoundingBoxOnly(bool b) {
 }
 //--------------------------------------------------------------------
 
+//--------------------------------------------------------------------
+void GateVImageVolume::SetMaxOutOfRangeFraction(double f){
+  GateMessage("Volume",3,"setting a new threshold for the fraction of voxels with out-of-range HU values: " << f);
+  if ( (f>=0.) && (f<=1.) ){
+    mMaxOutOfRangeFraction=f;
+    return;
+  }
+  GateError("Fraction should be given as a number between 0.0 and 1.0!!");
+}
+//--------------------------------------------------------------------
+
 
 //--------------------------------------------------------------------
 void GateVImageVolume::SetIsoCenter(const G4ThreeVector & i)
@@ -347,8 +361,8 @@ void GateVImageVolume::LoadImageMaterialsFromHounsfieldTable() {
 
   //FIXME: remove these two lines, we should load the HU-file as-is. It's up to the user to make sure it works.
   // GetOutsideValue returns the lowest value found in the image - 1. NOT the lowest value in mHounsfieldToImageMaterialTableFilename
-  G4String parentMat = GetParentVolume()->GetMaterialName();
-  mHounsfieldMaterialTable.AddMaterial(pImage->GetOutsideValue(),pImage->GetOutsideValue(),parentMat);
+  // G4String parentMat = GetParentVolume()->GetMaterialName();
+  // mHounsfieldMaterialTable.AddMaterial(pImage->GetOutsideValue(),pImage->GetOutsideValue(),parentMat);
 
   double low  =  1e6; //must start oppositely for the comparisons to work.
   double high = -1e6;
@@ -390,14 +404,13 @@ void GateVImageVolume::LoadImageMaterialsFromHounsfieldTable() {
   GateMessage("Volume",5,"HUMinValue   : " << low << ", HUMaxValue: " << high << Gateendl);
 
   if (pImage->GetMinValue() < low || pImage->GetMaxValue() > high) {
-    GateError("The image contains HU indices out of range of the HU range found in " <<
-              mHounsfieldToImageMaterialTableFilename << Gateendl <<
-              "HU    min, max: " << low << ", " << high << Gateendl <<
-              "Image min, max: " << pImage->GetMinValue() << ", " << pImage->GetMaxValue() << Gateendl <<
-              "Abort." << Gateendl);
-  }
-  //  if (mHounsfieldMaterialTable.GetNumberOfMaterials() == 0) {
-  if (mHounsfieldMaterialTable.GetNumberOfMaterials() == 1 ) {//there is a default mat = worldDefaultAir
+    GateWarning( "The image contains HU indices out of range of the HU range found in " <<
+                 mHounsfieldToImageMaterialTableFilename << Gateendl <<
+                 "HU    min, max: " << low << ", " << high << Gateendl <<
+                 "Image min, max: " << pImage->GetMinValue() << ", " << pImage->GetMaxValue() << Gateendl );
+    // GateError( "Abort." << Gateendl);
+  }
+  if (mHounsfieldMaterialTable.GetNumberOfMaterials() == 0 ) {
     GateError("No Hounsfield material defined in the file "
               << mHounsfieldToImageMaterialTableFilename << ". Abort.\n");
   }
@@ -411,21 +424,48 @@ void GateVImageVolume::LoadImageMaterialsFromHounsfieldTable() {
   while (iter != pImage->end()) {
     double label = mHounsfieldMaterialTable.GetLabelFromH(*iter);
     if (label<0) {
-      GateError(" I find H=" << *iter
+      GateMessage("Volume",1," I find H=" << *iter
                 << " in the image, while Hounsfield range start at "
                 << mHounsfieldMaterialTable[0].mH1 << Gateendl);
+      label = 0;
+      ++mUnderflow;
     }
     if (label>=mHounsfieldMaterialTable.GetNumberOfMaterials()) {
-      GateError(" I find H=" << *iter
+      GateMessage("Volume",1," I find H=" << *iter
                 << " in the image, while Hounsfield range stop at "
                 << mHounsfieldMaterialTable[mHounsfieldMaterialTable.GetNumberOfMaterials()-1].mH2
                 << Gateendl);
+      label = mHounsfieldMaterialTable.GetNumberOfMaterials() - 1;
+      ++mOverflow;
     }
     //GateMessage("Core", 0, " pix = " << (*iter) << " lab = " << label << Gateendl);
     (*iter) = label;
     ++iter;
   }
 
+  assert( pImage->GetNumberOfValues() > 0 );
+  // double out_of_range_fraction = double(mUnderflow+mOverflow)/pImage->GetNumberOfValues(); // not yet
+  double out_of_range_fraction = double(mOverflow)/pImage->GetNumberOfValues();
+  if ( out_of_range_fraction > mMaxOutOfRangeFraction ){
+    int n_bad_max(std::floor(mMaxOutOfRangeFraction * pImage->GetNumberOfValues()));
+    GateMessage( "Volume", 0, "ERROR: too many HU values are out of the range of the materials table: " << Gateendl
+                 << "******** " << mUnderflow << " underflows (HU<" << low << ") ******** " << Gateendl
+                 << "******** " << mOverflow << " overflows (HU>" << high << ") ********" << Gateendl);
+    if ( (mUnderflow > 0) && ( std::abs( mHalfSize.x() - pImage->GetHalfSize().x()) > 0.1*pImage->GetVoxelSize().x() ) ){
+        const G4ThreeVector & size = pImage->GetResolution();
+        int n_margin = pImage->GetNumberOfValues() - (((int)lrint(size.x())-2)*((int)lrint(size.y())-2)*((int)lrint(size.z())-2));
+        GateMessage( "Volume", 0, "(" << n_margin << " of the " << mUnderflow << " underflows are in the 1 voxel extra margin.)" << Gateendl);
+    }
+    GateMessage( "Volume", 0,
+                 "(The maximum 'bad' fraction is " << mMaxOutOfRangeFraction
+              << " and the image has " << pImage->GetNumberOfValues() << " voxels, " << Gateendl
+              << "so maximum " << n_bad_max << " out-of-range HU values are allowed." << Gateendl
+              << "Possible solutions: fix your input image, "
+              << "or expand the HU range of the HU-to-materials table. "  << Gateendl
+              << "Or, if you really know what you are doing, "
+              << "you can set the 'setMaxOutOfRangeFraction' option to a nonzero value larger than " << out_of_range_fraction << " .)" << Gateendl );
+    GateError( "ABORT" );
+  }
   // Debug
   // for(uint i=0; i<mHounsfieldMaterialTable.GetH1Vector().size(); i++) {
   //     double h = mHounsfieldMaterialTable.GetH1Vector()[i];
@@ -846,6 +886,14 @@ void GateVImageVolume::DestroyOwnSolidAndLogicalVolume()
   pBoxLog = 0;
   if (pBoxSolid) delete pBoxSolid;
   pBoxSolid = 0;
+  // If I put these warnings in the destructor, then they will NEVER be displayed.
+  if (mUnderflow>0) GateMessage("Volume",0,"There were " << mUnderflow << " voxels with HU values less than the minimum HU value in the HU-to-materials table." << Gateendl);
+  if ( (mUnderflow > 0) && ( std::abs( mHalfSize.x() - pImage->GetHalfSize().x()) > 0.1*pImage->GetVoxelSize().x() ) ){
+    const G4ThreeVector & size = pImage->GetResolution();
+    int n_margin = pImage->GetNumberOfValues() - (((int)lrint(size.x())-2)*((int)lrint(size.y())-2)*((int)lrint(size.z())-2));
+    GateMessage( "Volume", 0, "(" << n_margin << " of the " << mUnderflow << " underflows are in the 1 voxel extra margin.)" << Gateendl);
+  }
+  if (mOverflow>0) GateMessage("Volume",0,"There were " << mOverflow << " voxels with HU values greater than the maximum HU value in the HU-to-materials table." << Gateendl);
 }
 //--------------------------------------------------------------------
 

diff --git a/source/geometry/src/GateVImageVolumeMessenger.cc b/source/geometry/src/GateVImageVolumeMessenger.cc
@@ -17,6 +17,7 @@
 #include "G4UIcmdWithAString.hh"
 #include "G4UIcmdWith3VectorAndUnit.hh"
 #include "G4UIcmdWithABool.hh"
+#include "G4UIcmdWithADouble.hh"
 
 //---------------------------------------------------------------------------
 GateVImageVolumeMessenger::GateVImageVolumeMessenger(GateVImageVolume* volume)
@@ -98,6 +99,10 @@ GateVImageVolumeMessenger::GateVImageVolumeMessenger(GateVImageVolume* volume)
   pDoNotBuildVoxelsCmd = 0;
   pDoNotBuildVoxelsCmd = new G4UIcmdWithABool(n,this);
   pDoNotBuildVoxelsCmd->SetGuidance("Only build the bounding box (no voxels !), for visualization purpose only.");
+
+  n = dir +"/setMaxOutOfRangeFraction";
+  pSetMaxOutOfRangeFractionCmd = new G4UIcmdWithADouble(n,this);
+  pSetMaxOutOfRangeFractionCmd->SetGuidance("Maximum fraction (number between 0.0 and 1.0) of voxels that have a HU value out of the range of the materials table.");
 }
 //---------------------------------------------------------------------------
 
@@ -119,6 +124,7 @@ GateVImageVolumeMessenger::~GateVImageVolumeMessenger()
   delete pBuildMassImageCmd;
   delete pDoNotBuildVoxelsCmd;
   delete pIsoCenterRotationFlagCmd;
+  delete pSetMaxOutOfRangeFractionCmd;
 }
 //---------------------------------------------------------------------------
 
@@ -165,6 +171,9 @@ void GateVImageVolumeMessenger::SetNewValue(G4UIcommand* command,
   else if (command == pIsoCenterRotationFlagCmd) {
     pVImageVolume->SetIsoCenterRotationFlag(pIsoCenterRotationFlagCmd->GetNewBoolValue(newValue));
   }
+  else if ( command == pSetMaxOutOfRangeFractionCmd) {
+    pVImageVolume->SetMaxOutOfRangeFraction(pSetMaxOutOfRangeFractionCmd->GetNewDoubleValue(newValue));
+  }
   // It is necessary to call GateVolumeMessenger::SetNewValue if the command
   // is not recognized
   else {