Back to home page

EIC code displayed by LXR

 
 

    


Warning, /geant4/examples/extended/medical/DICOM/DICOM1/README is written in an unsupported language. File is not indexed.

0001 # DICOM1 Example
0002 
0003 The DICOM application has been originally developed by the Geant4 users:
0004 Louis Archambault,+Luc Beaulieu, ++Vincent Hubert-Tremblay.
0005 
0006 + Centre Hospitalier Universitaire de Quebec (CHUQ),
0007 Hotel-Dieu de Quebec, departement de Radio-oncologie
0008 11 cote du palais. Quebec, QC, Canada, G1R 2J6
0009 tel (418) 525-4444 #6720
0010 fax (418) 691 5268
0011 web : thomson.phy.ulaval.ca/phys_med 
0012 
0013 ++ Université Laval, Québec (QC) Canada
0014 
0015 
0016 And it has been deeply reviewed by Pedro Arce in December 2007.
0017 Very small changes by Stephane Chauvie in January 2008.
0018 Stephane Chauvie, Oct 2009: changed Physics list; changes in DICOM read. 
0019 Stephane Chauvie and Andrea Armando; June 2010 adapted for reading  whatever DICOM file
0020 Jonathan Madsen, Nov 2013: updated DICOM to utilize multithreading now available in Geant4.10
0021 Alexander Howard, Dec 2023: updated scorer to give correct voxel index for both regular and nested geometries
0022 
0023 ## Introduction
0024 
0025 This example serves first to convert a DICOM file to a simple ASCII file, where the Hounsfield numbers are converted to materials and densities so that it can be used by GEANT4. It serves also to create a GEANT4 geometry based on the DICOM file information using the G4PhantomParameterisation. 
0026 
0027 You can find the phantom reproduced in the image PhantomCT.jpg.
0028 In the application the phantom is placed on a table.
0029  
0030 ## Run the example:
0031 
0032  - batch mode:
0033    - DICOM run.mac 
0034 
0035  - interactive mode:
0036    - DICOM
0037    the file vis.mac is read in order to visualise the phantom with OpenGL, DAWN or VRML
0038  
0039 ## Metadata file:
0040 
0041 The old version of "Data.dat" is found in "Data.dat.old", when the project is configured with DICOM_USE_DCMTK=OFF,
0042 "Data.dat.old" is copied into the binary directory at "Data.dat".
0043   - i.e. cp ${PROJECT_SOURCE_DIR}/Data.dat.old ${PROJECT_BINARY_DIR}/Data.dat
0044 
0045 The new version of "Data.dat" is found in "Data.dat.new", when the project is configured with DICOM_USE_DCMTK=ON,
0046 "Data.dat.new" is copied into the binary directory as "Data.dat".
0047   - i.e. cp ${PROJECT_SOURCE_DIR}/Data.dat.new ${PROJECT_BINARY_DIR}/Data.dat
0048 
0049 ### Metadata file, OLD version:
0050 
0051  The file Data.dat has the following information
0052   - A line with the compression value (used only to create the .g4dcm and .g4dcmb, not to read it)
0053   - A line with the number of files 
0054   - A line for each file name (to these names it will be added the suffix .dcm to read the DICOM files in their original format, and the suffix .g4dcm to read the text files that contain the DICOM information where the Hounsfield numbers have been converted to material and densities)
0055 
0056 In case you want to convert DICOM files to text files, it must have the following lines:
0057   - The number of materials you want to use
0058   - A line for each material describing its name and the upper bound of the density interval. The materials should be described in increasing order of density. The voxels with a density between 0. and the first upper bound will be assigned to the first material, those with a density between the first upper bound and the second upper bound will be assigned to the second material, etc.
0059 
0060 ### Metadata file, NEW version (based on DCMTK):
0061 
0062 As for the previous version, a Data.dat file has to be defined to manage the conversion options. The format of this file is though quite different from the previous version. The format of this file is based on tags (similary to the ASCII geometry files).
0063 The following tags should be used:
0064 
0065 :COMPRESSION level 
0066 
0067 Where "level" is the number of voxels that will be merged into one in the X and Y dimen-
0068 sions. The Hounsfield numbers of the voxels merged are averaged to give the
0069 resulting value for the new voxel.
0070 Example:
0071 :COMPRESSION 4 // 4 X 4 voxels will be merged, so that the number of voxels in X and Y dimensions will be reduced by a factor 4
0072 
0073 :FILE file_name
0074 These are the list of files (one line per file) in DICOM format that will be treated.
0075 They can be of modality CT, RTSTRUCT or RTPLAN (the code will automatically
0076 detect its modality and treat it correspondingly).
0077 Example:
0078 :FILE 1.dcm
0079 :FILE 2.dcm
0080 :FILE 3.dcm
0081 
0082 :CT2D  Hounsfield_number density
0083 These sets of value pairs build the calibration curve (linearly interpolating between them). In other words, each Hounsfield number is given a material density using a function that is built interpolating between this list of value pairs.
0084 Example:
0085 :CT2D -5000 0. 
0086 :CT2D -1000 0.01 
0087 :CT2D -400 0.602 
0088 :CT2D  300 1.145 
0089 :CT2D 2000 1.856
0090 
0091 :MATE material_name upper_bound_of_material_Hounsfield_number_interval
0092 This serves for the Hounsfield number to material name conversion. The voxels with a Hounsfield number between 0. and the first upper bound will be assigned to the first material, those with a Hounsfiled number between the first upper bound and the second upper bound will be assigned to the second material, etc.
0093 Example:
0094 :MATE G4_AIR -800
0095 :MATE G4_LUNG_ICRP -145
0096 :MATE G4_ADIPOSE_TISSUE_ICRP -60
0097 :MATE G4_WATER 0
0098 
0099 Alternatively to the use of :MATE, you can use the :MATE_DENS
0100 :MATE_DENS material_name upper_bound_of_material_density_interval
0101 This serves for the material density to material name conversion. The voxels with a density between 0. and the first upper bound will be assigned to the first material, those with a density between the first upper bound and the second upper bound will be assigned to the second material, etc.
0102 Example:
0103 :MATE_DENS G4_AIR 0.207
0104 :MATE_DENS G4_LUNG_ICRP 0.919
0105 :MATE_DENS G4_ADIPOSE_TISSUE_ICRP 0.979
0106 :MATE_DENS G4_WATER 1.01
0107 
0108 We recommend the use of :MATE instead of :MATE_DENS as this is the way is used more often in the literature.
0109 
0110 :FILE_OUT file_name
0111 Name of output file containing the DICOM information in ASCII format
0112 
0113 
0114 ## Conversion of Hounsfield numbers to materials:
0115 
0116 After reading the name of files from Data.dat, if a file .dcm is found, then it looks for the corresponding .g4dcm file and if not found creates it. 
0117 Each file corresponds to a Z slice. The Z slices will be merged at runtime to form a unique patient volume; therefore the different slices have to be contiguous in Z.
0118 
0119 The DICOM images pixel values represent CT (Hounsfield) numbers and they should be converted, first, to a given density and then to a material type.  The relation between CT number and density is more or less linear. 
0120 The file CT2Density.dat contains the calibration curve to convert CT (Hounsfield) number to physical density
0121 The assignment of material densities to materials is done following the information from the file Data.dat (see below). In this case we have used:
0122 
0123            #####################################################
0124            #     Density Range                   Material      #
0125            #---------------------------------------------------#
0126            #       mg/cm3                            -         #
0127            #---------------------------------------------------#
0128            #  [ 0.    , 0.207 )                 Air            #
0129            #  [ 0.207 , 0.481 )                 Lungs (inhale) #
0130            #  [ 0.481 , 0.919 )                 Lungs (exhale) #
0131            #  [ 0.919 , 0.979 )                 Adipose        #
0132            #  [ 0.979 , 1.004 )                 Breast         #
0133            #  [ 1.004 , 1.043 )                 Phantom        #
0134            #  [ 1.043 , 1.109 )                 Liver          #
0135            #  [ 1.109 , 1.113 )                 Muscle         #
0136            #  [ 1.113 , 1.496 )                 Trabecular Bone#
0137            #  [ 1.496 , 1.654 ]                 Dense Bone     #
0138            #####################################################
0139 
0140 Data taken from the International Commission on Radiation Units and measurements (ICRU) report 46 was used to build the materials (lung, liver, breast, bones, ...).
0141 
0142 When using the Digital Head Phantom, the CT2Density.dat is not used. The conversion is performed directly in the Dicom Handler.cc
0143 
0144 ## Splitting materials in density intervals:
0145 
0146 In the class DicomDetectorConstruction, it is defined a density interval
0147 
0148   G4double densityDiff = 0.1;
0149 
0150 This means that the voxels of each material will be grouped in density intervals of 0.1 g/cm3 and a new material will be created for each group of voxels.
0151 
0152 ## Voxel colouring:
0153 
0154 The file Colormap.dat defines the colour that will be assigned to the voxels of each material.
0155 
0156 ## DICOM file formats:
0157 
0158 The DICOM files are converted to a simple text format. You may create your own file with the following format (see e.g. 14196616.g4dcm):
0159 
0160 - A line with the number of materials
0161 - A line for each material with its index and name (the same name of materials that you construct as G4Material's)
0162 - A line with the number of voxels in X, Y and Z
0163 - A line with the minimum and maximum extension in X (mm)
0164 - A line with the minimum and maximum extension in Y (mm)
0165 - A line with the minimum and maximum extension in Z (mm)
0166 - A number of lines containing the nVoxelX*nVoxelY*nVoxelZ material indices (one per voxel)
0167 - A number of lines containing the nVoxelX*nVoxelY*nVoxelZ material densities (one per voxel)
0168 
0169 As commented before the DICOM files (.dcm) are assumed to describe one Z slice per file, and therefore the GEANT4 text files (.g4dcm) created from them have also one unique Z slice per file. Nevertheless if you create your own .g4dcm file you may include as many Z slices as desired. In any case you have to respect the rule that the Z slices must be contiguous.
0170 
0171 The same information is also used to fill a file in binary format, that contains the same information as the text format. Its name ends in .g4dcmb, instead of .g4dcm .
0172 
0173 ## Choosing different parameterisation/navigation options:
0174 
0175 There are four possible ways in GEANT4 to treat the navigation in regular voxelised volumes:
0176 
0177 - The 3D optimisation with G4SmartVoxel: a 3D grid is built, so that the location of voxels is fast, but it requires a lot of memory
0178 - Using G4NestedParameterisation. The search is done hierarchically in X, Y and Z. It is fast and does not require big memory
0179 - Using G4PhantomParameterisation/G4RegularNavigation: an special algorithm to navigate in regular voxelised geometries (see GEANT4 doc). This is the fastest way without any extra memory requirement (and it is the default in this example). It includes an option (default) to skip frontiers between voxels when they have the same material. When using this option at each step the energy is all deposited in the last voxel; for properly distribution of the dose (=energy/volume) the G4PSDoseDeposit scorer can be used for regular (see (10) below) and G4PSDoseDeposit3D for nested parameterisation (see (11) below).
0180 
0181 Obsolete option:
0182 - Use 1D optimisation in replica. It will be very slow because each time a track exits a voxel it has to loop to all other voxels in a 2D slide in order to locate which one it will enter.
0183 
0184 You can select amongst the four options in the following way:
0185 
0186 - By default the example will run with G4RegularNavigation
0187 
0188 - To use the first option at RegularDicomDetectorConstruction.cc you must set
0189 
0190   patient_phys->SetRegularStructureId(0);
0191 
0192 - To use the second option (Nested Parameterisation) you must set the enviromental variable DICOM_NESTED_PARAM to 1
0193 
0194 - To use the final, obsolete 1D-option, apart from the change above at RegularDicomDetectorConstruction.cc you need to replace (i.e. use kUndefined)
0195 
0196   G4PVParameterised * patient_phys = new G4PVParameterised("Patient",voxel_logic,container_logic,
0197                         kUndefined, nVoxelX*nVoxelY*nVoxelZ, param);
0198 by 
0199   G4PVParameterised * patient_phys = new G4PVParameterised("Patient",voxel_logic,container_logic,
0200                         kXAxis, nVoxelX*nVoxelY*nVoxelZ, param);
0201 
0202   Note also you must *not* set the enviromental variable DICOM_NESTED_PARAM.
0203 
0204 
0205 ## Calculating dose in phantom voxels for regular navigation
0206 
0207 As mentioned above the regular navigation has the option to keip voxel frontiers when two voxels share the same material, what can make the CPU time several times smaller. But this option makes that all energy deposited is computed in the last voxel, instead of distributing it along the voxels traversed. To properly calculate the dose in each voxel the G4PSDoseDeposit scorer can be used.
0208 
0209 It takes into account the fact that,  when the particle travels through the voxels it looses energy and therefore the energy lost per length (dEdx) is bigger and also the effect of the multiple scattering is bigger. 
0210 The algorithm to make this correction is an iterative one, as the step length increase due multiple scattering (that converts the geometrical step length in what we will call the true step length) and the energy loss are correlated.
0211 It works in the folloing way: first the total true step length is distributed among the voxels proportionally to their geometrical step length; with these values it is calculated one voxel after another the value of dEdx and then the value of the kinetic energy at the entrance of each voxel; with these values it is calculated the geometrical to true step corrections due to multiple scattering for each voxel; finally these new values are used to recalculate the energy lost in each voxel. It has been demonstrated for dose in a water phantom and in a real phantom that the two-step iteration described is enough to reproduce the dose calcualted when no skipping of voxel frontiers is done.
0212 
0213 This scorer is implemented in this examples if the regular navigation option is 
0214 chosen. It is triggered at the method RegularDicomDetectorConstruction::ConstructPhantom() by the call
0215 
0216 SetScorer(voxel_logic);
0217 
0218 ## Calculating dose in phantom voxels for nested parameterisation
0219 
0220 For the nested parameterisation the geometry comprises replicas in X and Y which are then parameterised in Z. This means that to get the correct voxel idendification the replica depth has to be taken into account. The G4PSDoseDeposit3D scorers uses a fixed algorithm to calculate the voxel ID, according to the number
0221 of voxels in each axes and the associated replica depth. G4PSDoseDeposit3D("DoseDeposit", fNoVoxelsZ, fNoVoxelsY, fNoVoxelsX, 0, 2, 1) contains the number of voxels at the top level (0) and then two daughter levels down for the Y-voxels and one depth down for X.
0222 
0223 ## Output
0224 dicom.out is produced running the macro file run.mac. It has 2 columns: the first is the number of 
0225 voxel (ordered in x,y,z) and the second the dose there deposited (in Gy)
0226 It is produced, as an example, with a compression value of 32
0227 
0228 
0229 ## Partial phantom 
0230 It is possible to create a partial phantom, that is the intersection of a phantom with a volume. You may define the volume with the command 
0231 
0232 /dicom/intersectWithUserVolume 0. 0. 0. 45.*deg 0. 0. TUBE 0. 150. 100.
0233 
0234 where the first three arguments are its position, its second three arguments are the rotation around the global X, Y and Z axis and the rest of the parameters are the same that you use to build a solid using the ASCII geometry format
0235 
0236 Alternatively you can intersect the phantom with an existing Geant4 volume with the command
0237 
0238 /dicom/intersectWithG4Volume VOLUME_NAME
0239 
0240 The job will create an ASCII file names "phantom.g4pdcm" containing the partial phantom. To read this file all what is needed is to set the enviromental variable DICOM_PARTIAL_PARAM to 1
0241 
0242 ## Visualisation
0243 
0244 The Geant4 drivers are not meant for visualizing millions of voxel and visualising the DICOM geometries can be very computationally demanding.
0245 The users may want to visualise each DICOM slice separately or use higher compression values when visualising a part of DICOM project.
0246 
0247 
0248