Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/DICOM/src/DicomIntersectVolume.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/medical/DICOM/src/DicomIntersectVolume.cc (Version 11.3.0) and /examples/extended/medical/DICOM/src/DicomIntersectVolume.cc (Version 11.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 /// \file medical/DICOM/src/DicomIntersectVolu     27 /// \file medical/DICOM/src/DicomIntersectVolume.cc
 28 /// \brief Implementation of the DicomIntersec     28 /// \brief Implementation of the DicomIntersectVolume class
 29 //                                                 29 //
 30                                                    30 
 31 #include "DicomIntersectVolume.hh"                 31 #include "DicomIntersectVolume.hh"
 32                                                    32 
                                                   >>  33 #include "G4UIcmdWithAString.hh"
 33 #include "G4LogicalVolume.hh"                      34 #include "G4LogicalVolume.hh"
 34 #include "G4LogicalVolumeStore.hh"                 35 #include "G4LogicalVolumeStore.hh"
 35 #include "G4Material.hh"                       << 
 36 #include "G4PVParameterised.hh"                << 
 37 #include "G4PhantomParameterisation.hh"        << 
 38 #include "G4PhysicalVolumeStore.hh"            << 
 39 #include "G4SystemOfUnits.hh"                  << 
 40 #include "G4UIcmdWithAString.hh"               << 
 41 #include "G4UIcommand.hh"                      << 
 42 #include "G4VPhysicalVolume.hh"                    36 #include "G4VPhysicalVolume.hh"
                                                   >>  37 #include "G4PhysicalVolumeStore.hh"
 43 #include "G4VSolid.hh"                             38 #include "G4VSolid.hh"
 44 #include "G4tgbVolume.hh"                      << 
 45 #include "G4tgrSolid.hh"                           39 #include "G4tgrSolid.hh"
                                                   >>  40 #include "G4tgbVolume.hh"
                                                   >>  41 #include "G4Material.hh"
                                                   >>  42 #include "G4PhantomParameterisation.hh"
                                                   >>  43 #include "G4PVParameterised.hh"
                                                   >>  44 #include "G4UIcommand.hh"
                                                   >>  45 #include "G4SystemOfUnits.hh"
 46                                                    46 
 47 //....oooOO0OOooo........oooOO0OOooo........oo     47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 48 DicomIntersectVolume::DicomIntersectVolume()       48 DicomIntersectVolume::DicomIntersectVolume()
 49   : G4UImessenger(), fG4VolumeCmd(0), fSolid(0 <<  49  : G4UImessenger(),
                                                   >>  50    fG4VolumeCmd(0),
                                                   >>  51    fSolid(0),
                                                   >>  52    fPhantomMinusCorner(),
                                                   >>  53    fVoxelIsInside(0)
 50 {                                                  54 {
 51   fUserVolumeCmd = new G4UIcmdWithAString("/di <<  55   fUserVolumeCmd= new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
 52   fUserVolumeCmd->SetGuidance(                 <<  56   fUserVolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume "
 53     "Intersects a phantom with a user-defined  <<  57      "and outputs the voxels that are totally inside the intersection as"
 54     "and outputs the voxels that are totally i <<  58      " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
 55     " a new phantom file. It must have the par <<  59      "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
 56     "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_ <<  60   fUserVolumeCmd->SetParameterName("choice",true);
 57   fUserVolumeCmd->SetParameterName("choice", t << 
 58   fUserVolumeCmd->AvailableForStates(G4State_I     61   fUserVolumeCmd->AvailableForStates(G4State_Idle);
 59                                                    62 
 60   fG4VolumeCmd = new G4UIcmdWithAString("/dico <<  63   fG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
 61   fG4VolumeCmd->SetGuidance(                   <<  64   fG4VolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume"
 62     "Intersects a phantom with a user-defined  <<  65      " and outputs the voxels that are totally inside the intersection as "
 63     " and outputs the voxels that are totally  <<  66      " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
 64     " a new phantom file. It must have the par <<  67      "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
 65     "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_ <<  68   fG4VolumeCmd->SetParameterName("choice",true);
 66   fG4VolumeCmd->SetParameterName("choice", tru << 
 67   fG4VolumeCmd->AvailableForStates(G4State_Idl     69   fG4VolumeCmd->AvailableForStates(G4State_Idle);
 68 }                                                  70 }
 69                                                    71 
 70 //....oooOO0OOooo........oooOO0OOooo........oo     72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 71 DicomIntersectVolume::~DicomIntersectVolume()      73 DicomIntersectVolume::~DicomIntersectVolume()
 72 {                                                  74 {
 73   if (fUserVolumeCmd) delete fUserVolumeCmd;   <<  75  if (fUserVolumeCmd) delete fUserVolumeCmd;
 74   if (fG4VolumeCmd) delete fG4VolumeCmd;       <<  76  if (fG4VolumeCmd) delete fG4VolumeCmd;
 75 }                                                  77 }
 76                                                    78 
 77 //....oooOO0OOooo........oooOO0OOooo........oo     79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 78 void DicomIntersectVolume::SetNewValue(G4UIcom <<  80 void DicomIntersectVolume::SetNewValue(G4UIcommand * command,
 79 {                                              <<  81                                              G4String newValues)    
                                                   >>  82 { 
 80   G4AffineTransform theVolumeTransform;            83   G4AffineTransform theVolumeTransform;
 81                                                    84 
 82   if (command == fUserVolumeCmd) {             <<  85   if (command == fUserVolumeCmd)
 83     std::vector<G4String> params = GetWordsInS <<  86   {
 84     if (params.size() < 8) {                   <<  87     std::vector<G4String> params = GetWordsInString( newValues );
                                                   >>  88     if( params.size() < 8 ) {
 85       G4Exception("DicomIntersectVolume::SetNe     89       G4Exception("DicomIntersectVolume::SetNewValue",
 86                   " There must be at least 8 p <<  90           " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
 87                   "ANG_X ANG_Y ANG_Z SOLID_PAR <<  91           "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
 88                   FatalErrorInArgument,            92                   FatalErrorInArgument,
 89                   G4String("Number of paramete <<  93                   G4String("Number of parameters given = " + 
 90                            + G4UIcommand::Conv <<  94                G4UIcommand::ConvertToString( G4int(params.size()) )).c_str());
 91                     .c_str());                 <<  95       
 92     }                                              96     }
 93                                                    97 
 94     //----- Build G4VSolid                     <<  98     //----- Build G4VSolid 
 95     BuildUserSolid(params);                        99     BuildUserSolid(params);
 96                                                   100 
 97     //----- Calculate volume inverse 3D transf    101     //----- Calculate volume inverse 3D transform
 98     G4ThreeVector pos = G4ThreeVector(G4UIcomm << 102     G4ThreeVector pos = G4ThreeVector( G4UIcommand::ConvertToDouble(params[0]), 
 99                                       G4UIcomm << 103                        G4UIcommand::ConvertToDouble(params[1]),
100                                       G4UIcomm << 104                        G4UIcommand::ConvertToDouble(params[2]) );
101     G4RotationMatrix* rotmat = new G4RotationM    105     G4RotationMatrix* rotmat = new G4RotationMatrix;
102     std::vector<G4double> angles;              << 106     std::vector<G4double> angles; 
103     rotmat->rotateX(G4UIcommand::ConvertToDoub << 107     rotmat->rotateX( G4UIcommand::ConvertToDouble(params[3]) );
104     rotmat->rotateY(G4UIcommand::ConvertToDoub << 108     rotmat->rotateY( G4UIcommand::ConvertToDouble(params[4]) );
105     rotmat->rotateY(G4UIcommand::ConvertToDoub << 109     rotmat->rotateY( G4UIcommand::ConvertToDouble(params[5]) );
106     theVolumeTransform = G4AffineTransform(rot << 110     theVolumeTransform = G4AffineTransform( rotmat, pos ).Invert();
107   }                                            << 111 
108   else if (command == fG4VolumeCmd) {          << 112   }
109     std::vector<G4String> params = GetWordsInS << 113   else if (command == fG4VolumeCmd)
110     if (params.size() != 1)                    << 114   {
111       G4Exception("DicomIntersectVolume::SetNe << 115     std::vector<G4String> params = GetWordsInString( newValues );
112                   G4String("Command: " + comma << 116     if( params.size() !=1 ) 
113                            + " " + newValues + << 117       G4Exception("DicomIntersectVolume::SetNewValue",
114                     .c_str());                 << 118                   "",
                                                   >> 119                   FatalErrorInArgument,
                                                   >> 120                   G4String("Command: "+ command->GetCommandPath() + 
                                                   >> 121                            "/" + command->GetCommandName() + 
                                                   >> 122                            " " + newValues + 
                                                   >> 123                            "  needs 1 argument: VOLUME_NAME").c_str()); 
115                                                   124 
116     //----- Build G4VSolid                     << 125     //----- Build G4VSolid 
117     BuildG4Solid(params);                         126     BuildG4Solid(params);
118                                                   127 
119     //----- Calculate volume inverse 3D transf    128     //----- Calculate volume inverse 3D transform
120     G4VPhysicalVolume* pv = GetPhysicalVolumes << 129     G4VPhysicalVolume* pv = GetPhysicalVolumes( params[0], 1, 1)[0];
121                                                << 130     
122     theVolumeTransform = G4AffineTransform(pv- << 131     theVolumeTransform = 
                                                   >> 132       G4AffineTransform( pv->GetFrameRotation(), pv->GetFrameTranslation() );
123   }                                               133   }
124                                                   134 
125   //----- Calculate relative phantom - volume  << 135   //----- Calculate relative phantom - volume 3D transform 
126   G4PhantomParameterisation* thePhantomParam =    136   G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
127                                                << 137   
128   G4RotationMatrix* rotph = new G4RotationMatr << 138   G4RotationMatrix* rotph = new G4RotationMatrix(); 
129   // assumes the phantom mother is not rotated << 139   // assumes the phantom mother is not rotated !!! 
130   G4AffineTransform thePhantomTransform(rotph, << 140   G4AffineTransform thePhantomTransform( rotph, G4ThreeVector() ); 
131   // assumes the phantom mother is not transla << 141   // assumes the phantom mother is not translated !!! 
132                                                << 142 
133   G4AffineTransform theTransform = theVolumeTr << 143   G4AffineTransform theTransform = theVolumeTransform*thePhantomTransform;
134                                                << 144 
135   G4ThreeVector axisX(1., 0., 0.);             << 145   G4ThreeVector axisX( 1., 0., 0. );
136   G4ThreeVector axisY(0., 1., 0.);             << 146   G4ThreeVector axisY( 0., 1., 0. );
137   G4ThreeVector axisZ(0., 0., 1.);             << 147   G4ThreeVector axisZ( 0., 0., 1. );
138   theTransform.ApplyAxisTransform(axisX);         148   theTransform.ApplyAxisTransform(axisX);
139   theTransform.ApplyAxisTransform(axisY);         149   theTransform.ApplyAxisTransform(axisY);
140   theTransform.ApplyAxisTransform(axisZ);         150   theTransform.ApplyAxisTransform(axisZ);
141                                                << 151   
142   //----- Write phantom header                    152   //----- Write phantom header
143   G4String thePhantomFileName = "phantom.g4pdc    153   G4String thePhantomFileName = "phantom.g4pdcm";
144   fout.open(thePhantomFileName);                  154   fout.open(thePhantomFileName);
145   std::vector<G4Material*> materials = thePhan    155   std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
146   fout << materials.size() << G4endl;             156   fout << materials.size() << G4endl;
147   for (unsigned int ii = 0; ii < materials.siz << 157   for( unsigned int ii = 0; ii < materials.size(); ++ii )
                                                   >> 158   {
148     fout << ii << " " << materials[ii]->GetNam    159     fout << ii << " " << materials[ii]->GetName() << G4endl;
149   }                                               160   }
150                                                   161 
151   //----- Loop to pantom voxels                   162   //----- Loop to pantom voxels
152   G4int nx = G4int(thePhantomParam->GetNoVoxel    163   G4int nx = G4int(thePhantomParam->GetNoVoxelsX());
153   G4int ny = G4int(thePhantomParam->GetNoVoxel    164   G4int ny = G4int(thePhantomParam->GetNoVoxelsY());
154   G4int nz = G4int(thePhantomParam->GetNoVoxel    165   G4int nz = G4int(thePhantomParam->GetNoVoxelsZ());
155   G4int nxy = nx * ny;                         << 166   G4int nxy = nx*ny;
156   fVoxelIsInside = new G4bool[nx * ny * nz];   << 167   fVoxelIsInside = new G4bool[nx*ny*nz];
157   G4double voxelHalfWidthX = thePhantomParam->    168   G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
158   G4double voxelHalfWidthY = thePhantomParam->    169   G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
159   G4double voxelHalfWidthZ = thePhantomParam->    170   G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
160                                                   171 
161   //----- Write phantom dimensions and limits     172   //----- Write phantom dimensions and limits
162   fout << nx << " " << ny << " " << nz << G4en    173   fout << nx << " " << ny << " " << nz << G4endl;
163   fout << -voxelHalfWidthX * nx + thePhantomTr << 174   fout << -voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << " "  
164        << voxelHalfWidthX * nx + thePhantomTra << 175        << voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << G4endl;
165   fout << -voxelHalfWidthY * ny + thePhantomTr << 176   fout << -voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << " " 
166        << voxelHalfWidthY * ny + thePhantomTra << 177        << voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << G4endl;
167   fout << -voxelHalfWidthZ * nz + thePhantomTr << 178   fout << -voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << " " 
168        << voxelHalfWidthZ * nz + thePhantomTra << 179        << voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << G4endl;
                                                   >> 180 
                                                   >> 181   for( G4int iz = 0; iz < nz; ++iz ){
                                                   >> 182     for( G4int iy = 0; iy < ny; ++iy) {
169                                                   183 
170   for (G4int iz = 0; iz < nz; ++iz) {          << 
171     for (G4int iy = 0; iy < ny; ++iy) {        << 
172       G4bool bPrevVoxelInside = true;             184       G4bool bPrevVoxelInside = true;
173       G4bool b1VoxelFoundInside = false;          185       G4bool b1VoxelFoundInside = false;
174       G4int firstVoxel = -1;                      186       G4int firstVoxel = -1;
175       G4int lastVoxel = -1;                       187       G4int lastVoxel = -1;
176       for (G4int ix = 0; ix < nx; ++ix) {      << 188       for(G4int ix = 0; ix < nx; ++ix ){
177         G4ThreeVector voxelCentre((-nx + ix *  << 189         G4ThreeVector voxelCentre( (-nx+ix*2+1)*voxelHalfWidthX, 
178                                   (-ny + iy *  << 190                    (-ny+iy*2+1)*voxelHalfWidthY, (-nz+iz*2+1)*voxelHalfWidthZ);
179                                   (-nz + iz *  << 
180         theTransform.ApplyPointTransform(voxel    191         theTransform.ApplyPointTransform(voxelCentre);
181         G4bool bVoxelIsInside = true;             192         G4bool bVoxelIsInside = true;
182         for (G4int ivx = -1; ivx <= 1; ivx +=  << 193         for( G4int ivx = -1; ivx <= 1; ivx+=2 ) {
183           for (G4int ivy = -1; ivy <= 1; ivy + << 194           for( G4int ivy = -1; ivy <= 1; ivy+=2 ){
184             for (G4int ivz = -1; ivz <= 1; ivz << 195             for( G4int ivz = -1; ivz <= 1; ivz+=2 ) {
185               G4ThreeVector voxelPoint = voxel << 196               G4ThreeVector voxelPoint = voxelCentre 
186                                          + ivy << 197                 + ivx*voxelHalfWidthX*axisX +
187                                          + ivz << 198                 ivy*voxelHalfWidthY*axisY + ivz*voxelHalfWidthZ*axisZ;
188               if (fSolid->Inside(voxelPoint) = << 199               if( fSolid->Inside( voxelPoint ) == kOutside ) {
189                 bVoxelIsInside = false;           200                 bVoxelIsInside = false;
190                 break;                            201                 break;
191               }                                << 202               } else {
192               else {                           << 
193               }                                   203               }
194             }                                     204             }
195             if (!bVoxelIsInside) break;        << 205             if( !bVoxelIsInside ) break;
196           }                                       206           }
197           if (!bVoxelIsInside) break;          << 207           if( !bVoxelIsInside ) break;
198         }                                         208         }
199                                                   209 
200         G4int copyNo = ix + nx * iy + nxy * iz << 210         G4int copyNo = ix + nx*iy + nxy*iz;
201         if (bVoxelIsInside) {                  << 211         if( bVoxelIsInside ) {
202           if (!bPrevVoxelInside) {             << 212           if( !bPrevVoxelInside ) {
203             G4Exception("DicomIntersectVolume: << 213             G4Exception("DicomIntersectVolume::SetNewValue",
204                         "Volume cannot interse << 214                         "",
205                         "please use other voxe << 215                         FatalException,
                                                   >> 216                  "Volume cannot intersect phantom in discontiguous voxels, "
                                                   >> 217                  "please use other voxel");
206           }                                       218           }
207           if (!b1VoxelFoundInside) {           << 219           if( !b1VoxelFoundInside ) {
208             firstVoxel = ix;                      220             firstVoxel = ix;
209             b1VoxelFoundInside = true;            221             b1VoxelFoundInside = true;
210           }                                       222           }
211           lastVoxel = ix;                         223           lastVoxel = ix;
212           fVoxelIsInside[copyNo] = true;          224           fVoxelIsInside[copyNo] = true;
213         }                                      << 225         } else {
214         else {                                 << 
215           fVoxelIsInside[copyNo] = false;         226           fVoxelIsInside[copyNo] = false;
216         }                                         227         }
                                                   >> 228            
217       }                                           229       }
218       fout << firstVoxel << " " << lastVoxel <    230       fout << firstVoxel << " " << lastVoxel << G4endl;
219     }                                             231     }
220   }                                               232   }
221                                                   233 
222   //-----  Now write the materials                234   //-----  Now write the materials
223   for (G4int iz = 0; iz < nz; ++iz) {          << 235   for( G4int iz = 0; iz < nz; ++iz ){
224     for (G4int iy = 0; iy < ny; ++iy) {        << 236     for( G4int iy = 0; iy < ny; ++iy) {
225       G4bool b1xFound = false;                 << 237       G4bool b1xFound = false; 
226       for (G4int ix = 0; ix < nx; ++ix) {      << 238       for(G4int ix = 0; ix < nx; ++ix ){
227         size_t copyNo = ix + ny * iy + nxy * i << 239         size_t copyNo = ix + ny*iy + nxy*iz;
228         if (fVoxelIsInside[copyNo]) {          << 240         if( fVoxelIsInside[copyNo] ) {
229           fout << thePhantomParam->GetMaterial << 241           fout << thePhantomParam->GetMaterialIndex(copyNo)<< " ";
230           b1xFound = true;                        242           b1xFound = true;
231         }                                         243         }
232       }                                           244       }
233       if (b1xFound) fout << G4endl;            << 245       if(b1xFound ) fout << G4endl;
234     }                                             246     }
235   }                                               247   }
236                                                   248 
237   // Now write densities                          249   // Now write densities
238   for (G4int iz = 0; iz < nz; ++iz) {          << 250   for( G4int iz = 0; iz < nz; ++iz ){
239     for (G4int iy = 0; iy < ny; ++iy) {        << 251     for( G4int iy = 0; iy < ny; ++iy) {
240       G4bool b1xFound = false;                 << 252       G4bool b1xFound = false; 
241       for (G4int ix = 0; ix < nx; ++ix) {      << 253       for(G4int ix = 0; ix < nx; ++ix ){
242         size_t copyNo = ix + ny * iy + nxy * i << 254         size_t copyNo = ix + ny*iy + nxy*iz;
243         if (fVoxelIsInside[copyNo]) {          << 255         if( fVoxelIsInside[copyNo] ) {
244           fout << thePhantomParam->GetMaterial << 256           fout <<thePhantomParam->GetMaterial(copyNo)->GetDensity()/g*cm3<< " ";
245           b1xFound = true;                        257           b1xFound = true;
246         }                                         258         }
247       }                                           259       }
248       if (b1xFound) fout << G4endl;            << 260       if(b1xFound ) fout << G4endl;
249     }                                             261     }
250   }                                               262   }
251 }                                              << 263 
                                                   >> 264 } 
252                                                   265 
253 //....oooOO0OOooo........oooOO0OOooo........oo    266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 void DicomIntersectVolume::BuildUserSolid(std: << 267 void DicomIntersectVolume::BuildUserSolid( std::vector<G4String> params )
255 {                                                 268 {
256   for (G4int ii = 0; ii < 6; ++ii)             << 269   for( G4int ii = 0; ii < 6; ++ii ) params.erase( params.begin() ); 
257     params.erase(params.begin());              << 
258   // take otu position and rotation angles        270   // take otu position and rotation angles
259   params.insert(params.begin(), ":SOLID");     << 271   params.insert( params.begin(), ":SOLID");
260   params.insert(params.begin(), params[1]);    << 272   params.insert( params.begin(), params[1] );
261   G4tgrSolid* tgrSolid = new G4tgrSolid(params    273   G4tgrSolid* tgrSolid = new G4tgrSolid(params);
262   G4tgbVolume* tgbVolume = new G4tgbVolume();  << 274   G4tgbVolume* tgbVolume = new G4tgbVolume();  
263   fSolid = tgbVolume->FindOrConstructG4Solid(t << 275   fSolid = tgbVolume->FindOrConstructG4Solid( tgrSolid );
                                                   >> 276 
264 }                                                 277 }
265                                                   278 
266 //....oooOO0OOooo........oooOO0OOooo........oo    279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 void DicomIntersectVolume::BuildG4Solid(std::v << 280 void DicomIntersectVolume::BuildG4Solid( std::vector<G4String> params )
268 {                                                 281 {
269   fSolid = GetLogicalVolumes(params[0], 1, 1)[ << 282   fSolid = GetLogicalVolumes( params[0], 1, 1)[0]->GetSolid();
                                                   >> 283  
270 }                                                 284 }
271                                                   285 
272 //....oooOO0OOooo........oooOO0OOooo........oo    286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
273 G4PhantomParameterisation* DicomIntersectVolum << 287 G4PhantomParameterisation*
                                                   >> 288 DicomIntersectVolume::GetPhantomParam(G4bool bMustExist)
274 {                                                 289 {
275   G4PhantomParameterisation* paramreg = 0;        290   G4PhantomParameterisation* paramreg = 0;
276                                                   291 
277   G4PhysicalVolumeStore* pvs = G4PhysicalVolum    292   G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
278   for (auto cite = pvs->cbegin(); cite != pvs- << 293   for( auto cite = pvs->cbegin(); cite != pvs->cend(); ++cite )
279     if (IsPhantomVolume(*cite)) {              << 294   {
280       const G4PVParameterised* pvparam = stati << 295     if( IsPhantomVolume( *cite ) ) {
                                                   >> 296       const G4PVParameterised* pvparam = 
                                                   >> 297         static_cast<const G4PVParameterised*>(*cite);
281       G4VPVParameterisation* param = pvparam->    298       G4VPVParameterisation* param = pvparam->GetParameterisation();
282       //    if( static_cast<const G4PhantomPar    299       //    if( static_cast<const G4PhantomParameterisation*>(param) ){
283       //    if( static_cast<const G4PhantomPar    300       //    if( static_cast<const G4PhantomParameterisation*>(param) ){
284       //      G4cout << "G4PhantomParameterisa    301       //      G4cout << "G4PhantomParameterisation volume found  "
285       // << (*cite)->GetName() << G4endl;         302       // << (*cite)->GetName() << G4endl;
286       paramreg = static_cast<G4PhantomParamete    303       paramreg = static_cast<G4PhantomParameterisation*>(param);
287     }                                             304     }
288   }                                               305   }
289                                                << 306   
290   if (!paramreg && bMustExist)                 << 307   if( !paramreg && bMustExist ) 
291     G4Exception("DicomIntersectVolume::GetPhan << 308     G4Exception("DicomIntersectVolume::GetPhantomParam",
                                                   >> 309                 "",
                                                   >> 310                 FatalErrorInArgument,
292                 " No G4PhantomParameterisation    311                 " No G4PhantomParameterisation found ");
293                                                << 312   
294   return paramreg;                                313   return paramreg;
                                                   >> 314   
295 }                                                 315 }
296                                                   316 
297 //....oooOO0OOooo........oooOO0OOooo........oo    317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
298 std::vector<G4VPhysicalVolume*> DicomIntersect << 318 std::vector<G4VPhysicalVolume*> DicomIntersectVolume::GetPhysicalVolumes(
299                                                << 319                  const G4String& name, G4bool exists, G4int nVols )
300 {                                                 320 {
301   std::vector<G4VPhysicalVolume*> vvolu;          321   std::vector<G4VPhysicalVolume*> vvolu;
302   std::string::size_type ial = name.rfind(":")    322   std::string::size_type ial = name.rfind(":");
303   G4String volname = "";                          323   G4String volname = "";
304   G4int volcopy = 0;                              324   G4int volcopy = 0;
305   if (ial != std::string::npos) {              << 325   if( ial != std::string::npos )
306     std::string::size_type ial2 = name.rfind(" << 326   {
307     if (ial2 != std::string::npos) {           << 327     std::string::size_type ial2 = name.rfind(":",ial-1);
308       G4Exception("DicomIntersectVolume::GetPh << 328     if( ial2 != std::string::npos )
                                                   >> 329     {
                                                   >> 330       G4Exception("DicomIntersectVolume::GetPhysicalVolumes",
                                                   >> 331                   "",
                                                   >> 332                   FatalErrorInArgument,
309                   G4String("Name corresponds t    333                   G4String("Name corresponds to a touchable " + name).c_str());
310     }                                             334     }
311     else {                                     << 335     else
312       volname = name.substr(0, ial);           << 336     { 
313       volcopy = G4UIcommand::ConvertToInt(name << 337       volname = name.substr( 0, ial );
                                                   >> 338       volcopy = G4UIcommand::ConvertToInt( name.substr(ial+1, name.length()).c_str() );
314     }                                             339     }
315   }                                               340   }
316   else {                                       << 341   else
                                                   >> 342   {
317     volname = name;                               343     volname = name;
318     volcopy = -1;                                 344     volcopy = -1;
319   }                                               345   }
320                                                   346 
321   G4PhysicalVolumeStore* pvs = G4PhysicalVolum    347   G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
322   for (auto citepv = pvs->cbegin(); citepv !=  << 348   for( auto citepv = pvs->cbegin(); citepv != pvs->cend(); ++citepv )
323     if (volname == (*citepv)->GetName() && ((* << 349   {
324       vvolu.push_back(*citepv);                << 350     if( volname == (*citepv)->GetName() 
                                                   >> 351         && ( (*citepv)->GetCopyNo() == volcopy || -1 == volcopy ) )
                                                   >> 352     {
                                                   >> 353       vvolu.push_back( *citepv );
325     }                                             354     }
326   }                                               355   }
327                                                << 356   
328   //----- Check that at least one volume was f    357   //----- Check that at least one volume was found
329   if (vvolu.size() == 0) {                     << 358   if( vvolu.size() == 0 ) {
330     if (exists) {                              << 359     if(exists) {
331       G4Exception(" DicomIntersectVolume::GetP << 360       G4Exception(" DicomIntersectVolume::GetPhysicalVolumes",
332                   G4String("No physical volume << 361                "",
333     }                                          << 362                FatalErrorInArgument,
334     else {                                     << 363                G4String("No physical volume found with name " + name).c_str());
335       G4cerr << "!!WARNING: DicomIntersectVolu << 364     } else {
336              << " no physical volume found wit << 365       G4cerr << "!!WARNING: DicomIntersectVolume::GetPhysicalVolumes: "<<
                                                   >> 366         " no physical volume found with name " << name << G4endl;
337     }                                             367     }
338   }                                               368   }
339                                                   369 
340   if (nVols != -1 && G4int(vvolu.size()) != nV << 370   if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
341     G4Exception("DicomIntersectVolume::GetLogi    371     G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
342                 "Wrong number of physical volu << 372                 "Wrong number of physical volumes found",
343                 ("Number of physical volumes " << 373                 FatalErrorInArgument,
344                  + ", requesting " + G4UIcomma << 374                 ("Number of physical volumes " + 
345                   .c_str());                   << 375          G4UIcommand::ConvertToString(G4int(vvolu.size())) + 
346   }                                            << 376          ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
                                                   >> 377   } 
347                                                   378 
348   return vvolu;                                   379   return vvolu;
349 }                                                 380 }
350                                                   381 
351 //....oooOO0OOooo........oooOO0OOooo........oo    382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352 G4bool DicomIntersectVolume::IsPhantomVolume(G << 383 G4bool DicomIntersectVolume::IsPhantomVolume( G4VPhysicalVolume* pv )
353 {                                                 384 {
354   EAxis axis;                                     385   EAxis axis;
355   G4int nReplicas;                                386   G4int nReplicas;
356   G4double width, offset;                      << 387   G4double width,offset;
357   G4bool consuming;                               388   G4bool consuming;
358   pv->GetReplicationData(axis, nReplicas, widt << 389   pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
359   EVolume type = (consuming) ? kReplica : kPar    390   EVolume type = (consuming) ? kReplica : kParameterised;
360   if (type == kParameterised && pv->GetRegular << 391   if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {  
361     return TRUE;                                  392     return TRUE;
362   }                                            << 393   } else {
363   else {                                       << 
364     return FALSE;                                 394     return FALSE;
365   }                                               395   }
366 }                                              << 396 
                                                   >> 397 } 
367                                                   398 
368 //....oooOO0OOooo........oooOO0OOooo........oo    399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
369 std::vector<G4LogicalVolume*> DicomIntersectVo << 400 std::vector<G4LogicalVolume*> DicomIntersectVolume::GetLogicalVolumes( 
370                                                << 401 const G4String& name, G4bool exists, G4int nVols )
371 {                                                 402 {
372   //  G4cout << "GetLogicalVolumes " << name <    403   //  G4cout << "GetLogicalVolumes " << name << " " << exists << G4endl;
373   std::vector<G4LogicalVolume*> vvolu;            404   std::vector<G4LogicalVolume*> vvolu;
374   G4int ial = G4int(name.rfind(":"));             405   G4int ial = G4int(name.rfind(":"));
375   if (ial != -1) {                             << 406   if( ial != -1 ) {
376     G4Exception("DicomIntersectVolume::GetLogi << 407     G4Exception("DicomIntersectVolume::GetLogicalVolumes",
377                 G4String("Name corresponds to  << 408                 "",
                                                   >> 409                 FatalErrorInArgument,
                                                   >> 410                 G4String("Name corresponds to a touchable or physcal volume" 
                                                   >> 411                          + name).c_str());
378   }                                               412   }
379                                                   413 
380   G4LogicalVolumeStore* lvs = G4LogicalVolumeS    414   G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
381   for (auto citelv = lvs->cbegin(); citelv !=  << 415   for( auto citelv = lvs->cbegin(); citelv != lvs->cend(); ++citelv )
382     if (name == (*citelv)->GetName()) {        << 416   {
383       vvolu.push_back(*citelv);                << 417     if( name == (*citelv)->GetName() ) {
                                                   >> 418       vvolu.push_back( *citelv );
384     }                                             419     }
385   }                                               420   }
386                                                << 421   
387   //----- Check that at least one volume was f    422   //----- Check that at least one volume was found
388   if (vvolu.size() == 0) {                     << 423   if( vvolu.size() == 0 ) {
389     if (exists) {                              << 424     if(exists) {
390       G4Exception("DicomIntersectVolume::GetLo << 425       G4Exception("DicomIntersectVolume::GetLogicalVolumes:","",
391                   ("no logical volume found wi << 426           FatalErrorInArgument,("no logical volume found with name " 
392     }                                          << 427                                 + name).c_str());
393     else {                                     << 428     } else {
394       G4Exception("DicomIntersectVolume::GetLo << 429       G4Exception("DicomIntersectVolume::GetLogicalVolumes:","",
395                   ("no  logical volume found w << 430           JustWarning,("no  logical volume found with name " + name).c_str());
396     }                                             431     }
397   }                                               432   }
398                                                   433 
399   if (nVols != -1 && G4int(vvolu.size()) != nV << 434   if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
400     G4Exception("DicomIntersectVolume::GetLogi << 435     G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
                                                   >> 436                 "Wrong number of logical volumes found",
401                 FatalErrorInArgument,             437                 FatalErrorInArgument,
402                 ("Number of logical volumes "  << 438                 ("Number of logical volumes " +
403                  + ", requesting " + G4UIcomma << 439          G4UIcommand::ConvertToString(G4int(vvolu.size())) + 
404                   .c_str());                   << 440          ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
405   }                                            << 441   } 
406                                                   442 
407   return vvolu;                                   443   return vvolu;
                                                   >> 444     
408 }                                                 445 }
409                                                   446 
410 //....oooOO0OOooo........oooOO0OOooo........oo    447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411 std::vector<G4String> DicomIntersectVolume::Ge << 448 std::vector<G4String> DicomIntersectVolume::GetWordsInString( 
                                                   >> 449                                                 const G4String& stemp)
412 {                                                 450 {
413   std::vector<G4String> wordlist;                 451   std::vector<G4String> wordlist;
414                                                   452 
415   //---------- Read a line of file:               453   //---------- Read a line of file:
416   //----- Clear wordlist                          454   //----- Clear wordlist
417   G4int ii;                                       455   G4int ii;
418   const char* cstr = stemp.c_str();               456   const char* cstr = stemp.c_str();
419   G4int siz = G4int(stemp.length());              457   G4int siz = G4int(stemp.length());
420   G4int istart = 0;                               458   G4int istart = 0;
421   G4int nQuotes = 0;                              459   G4int nQuotes = 0;
422   G4bool lastIsBlank = false;                     460   G4bool lastIsBlank = false;
423   G4bool lastIsQuote = false;                     461   G4bool lastIsQuote = false;
424   for (ii = 0; ii < siz; ++ii) {               << 462   for( ii = 0; ii < siz; ++ii )
425     if (cstr[ii] == '\"') {                    << 463   {
426       if (lastIsQuote) {                       << 464     if(cstr[ii] == '\"' ){
427         G4Exception("GmGenUtils:GetWordsFromSt << 465       if( lastIsQuote ){
428                     ("There cannot be two quot << 466         G4Exception("GmGenUtils:GetWordsFromString","",FatalException,
                                                   >> 467             ("There cannot be two quotes together " + stemp).c_str() );
429       }                                           468       }
430       if (nQuotes % 2 == 1) {                  << 469       if( nQuotes%2 == 1 ){
431         // close word                          << 470         //close word 
432         wordlist.push_back(stemp.substr(istart << 471         wordlist.push_back( stemp.substr(istart,ii-istart) );
433         //        G4cout << "GetWordsInString  << 472         //        G4cout << "GetWordsInString new word0 " 
434         //<< wordlist[wordlist.size()-1] << "  << 473     //<< wordlist[wordlist.size()-1] << " istart " << istart 
435         //<< " ii " << ii << G4endl;              474         //<< " ii " << ii << G4endl;
436       }                                        << 475       } else {
437       else {                                   << 476         istart = ii+1;
438         istart = ii + 1;                       << 
439       }                                           477       }
440       ++nQuotes;                                  478       ++nQuotes;
441       lastIsQuote = true;                         479       lastIsQuote = true;
442       lastIsBlank = false;                        480       lastIsBlank = false;
443     }                                          << 481     } else if(cstr[ii] == ' ' ){
444     else if (cstr[ii] == ' ') {                << 482       //            G4cout << "GetWordsInString blank nQuotes " 
445       //            G4cout << "GetWordsInStrin << 
446       //<< nQuotes << " lastIsBlank " << lastI    483       //<< nQuotes << " lastIsBlank " << lastIsBlank << G4endl;
447       if (nQuotes % 2 == 0) {                  << 484       if( nQuotes%2 == 0 ){
448         if (!lastIsBlank && !lastIsQuote) {    << 485         if( !lastIsBlank && !lastIsQuote ) {
449           wordlist.push_back(stemp.substr(ista << 486           wordlist.push_back( stemp.substr(istart,ii-istart) );
450           //          G4cout << "GetWordsInStr << 487           //          G4cout << "GetWordsInString new word1 " 
451           //<< wordlist[wordlist.size()-1] <<  << 488       //<< wordlist[wordlist.size()-1] << " lastIsBlank " 
452           //<< lastIsBlank << G4endl;             489           //<< lastIsBlank << G4endl;
453         }                                         490         }
454                                                   491 
455         istart = ii + 1;                       << 492         istart = ii+1;
456         lastIsQuote = false;                      493         lastIsQuote = false;
457         lastIsBlank = true;                       494         lastIsBlank = true;
458       }                                           495       }
459     }                                          << 496     } else {
460     else {                                     << 497       if( ii == siz-1 ) {
461       if (ii == siz - 1) {                     << 498         wordlist.push_back( stemp.substr(istart,ii-istart+1) );
462         wordlist.push_back(stemp.substr(istart << 499         //                G4cout << "GetWordsInString new word2 " 
463         //                G4cout << "GetWordsI << 500     //<< wordlist[wordlist.size()-1] << " istart " << istart << G4endl;
464         //<< wordlist[wordlist.size()-1] << "  << 
465       }                                           501       }
466       lastIsQuote = false;                        502       lastIsQuote = false;
467       lastIsBlank = false;                        503       lastIsBlank = false;
468     }                                             504     }
469   }                                               505   }
470   if (nQuotes % 2 == 1) {                      << 506   if( nQuotes%2 == 1 ) {
471     G4Exception("GmGenUtils:GetWordsFromString << 507     G4Exception("GmGenUtils:GetWordsFromString","",FatalException, 
472                 ("unbalanced quotes in line "  << 508         ("unbalanced quotes in line " + stemp).c_str() );
473   }                                               509   }
474                                                   510 
475   return wordlist;                                511   return wordlist;
476 }                                                 512 }
477                                                   513