Warning, /geant4/examples/extended/medical/DICOM/DICOM1/README.md is written in an unsupported language. File is not indexed.
0001 \page ExampleDICOM1 Example DICOM1
0002
0003 # DICOM1 Example
0004
0005 The DICOM application has been originally developed by the Geant4 users:
0006 Louis Archambault,(1)Luc Beaulieu, (2)Vincent Hubert-Tremblay.
0007
0008 - (1) Centre Hospitalier Universitaire de Quebec (CHUQ),
0009 Hotel-Dieu de Quebec, departement de Radio-oncologie
0010 11 cote du palais. Quebec, QC, Canada, G1R 2J6
0011 tel (418) 525-4444 #6720
0012 fax (418) 691 5268
0013 web : thomson.phy.ulaval.ca/phys_med
0014 - (2) Université Laval, Québec (QC) Canada
0015
0016
0017 And it has been deeply reviewed by Pedro Arce in December 2007. \n
0018 Very small changes by Stephane Chauvie in January 2008. \n
0019 Stephane Chauvie, Oct 2009: changed Physics list; changes in DICOM read. \n
0020 Stephane Chauvie and Andrea Armando; June 2010 adapted for reading whatever DICOM file \n
0021 Jonathan Madsen, Nov 2013: updated DICOM to utilize multithreading now available in Geant4.10 \n
0022 Alexander Howard, Dec 2023: updated scorer to give correct voxel index for both regular and nested geometries
0023
0024 ## Introduction
0025
0026 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.
0027
0028 You can find the phantom reproduced in the image PhantomCT.jpg.
0029 In the application the phantom is placed on a table.
0030
0031 ## Run the example:
0032
0033 - batch mode:
0034 ```
0035 ./dicom1 run.mac
0036 ```
0037 - interactive mode:
0038 ```
0039 ./dicom1
0040 ```
0041 the file vis.mac is read in order to visualise the phantom with OpenGL, DAWN or VRML
0042
0043 ## Metadata file:
0044
0045 The old version of "Data.dat" is found in "Data.dat.old", when the project is configured with DICOM_USE_DCMTK=OFF,
0046 "Data.dat.old" is copied into the binary directory at "Data.dat".
0047 - i.e. cp ${PROJECT_SOURCE_DIR}/Data.dat.old ${PROJECT_BINARY_DIR}/Data.dat
0048
0049 The new version of "Data.dat" is found in "Data.dat.new", when the project is configured with DICOM_USE_DCMTK=ON,
0050 "Data.dat.new" is copied into the binary directory as "Data.dat".
0051 - i.e. cp ${PROJECT_SOURCE_DIR}/Data.dat.new ${PROJECT_BINARY_DIR}/Data.dat
0052
0053 ### Metadata file, OLD version:
0054
0055 The file Data.dat has the following information
0056 - A line with the compression value (used only to create the .g4dcm and .g4dcmb, not to read it)
0057 - A line with the number of files
0058 - 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)
0059
0060 In case you want to convert DICOM files to text files, it must have the following lines:
0061 - The number of materials you want to use
0062 - 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.
0063
0064 ### Metadata file, NEW version (based on DCMTK):
0065
0066 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).
0067 The following tags should be used:
0068
0069 :COMPRESSION level
0070
0071 Where "level" is the number of voxels that will be merged into one in the X and Y dimen-
0072 sions. The Hounsfield numbers of the voxels merged are averaged to give the
0073 resulting value for the new voxel.
0074 Example:
0075 :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
0076
0077 :FILE file_name
0078 These are the list of files (one line per file) in DICOM format that will be treated.
0079 They can be of modality CT, RTSTRUCT or RTPLAN (the code will automatically
0080 detect its modality and treat it correspondingly).
0081 Example:
0082 :FILE 1.dcm
0083 :FILE 2.dcm
0084 :FILE 3.dcm
0085
0086 :CT2D Hounsfield_number density
0087 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.
0088 Example:
0089 :CT2D -5000 0.
0090 :CT2D -1000 0.01
0091 :CT2D -400 0.602
0092 :CT2D 300 1.145
0093 :CT2D 2000 1.856
0094
0095 :MATE material_name upper_bound_of_material_Hounsfield_number_interval
0096 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.
0097 Example:
0098 :MATE G4_AIR -800
0099 :MATE G4_LUNG_ICRP -145
0100 :MATE G4_ADIPOSE_TISSUE_ICRP -60
0101 :MATE G4_WATER 0
0102
0103 Alternatively to the use of :MATE, you can use the :MATE_DENS
0104 :MATE_DENS material_name upper_bound_of_material_density_interval
0105 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.
0106 Example:
0107 :MATE_DENS G4_AIR 0.207
0108 :MATE_DENS G4_LUNG_ICRP 0.919
0109 :MATE_DENS G4_ADIPOSE_TISSUE_ICRP 0.979
0110 :MATE_DENS G4_WATER 1.01
0111
0112 We recommend the use of :MATE instead of :MATE_DENS as this is the way is used more often in the literature.
0113
0114 :FILE_OUT file_name
0115 Name of output file containing the DICOM information in ASCII format
0116
0117
0118 ## Conversion of Hounsfield numbers to materials:
0119
0120 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.
0121 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.
0122
0123 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.
0124 The file CT2Density.dat contains the calibration curve to convert CT (Hounsfield) number to physical density
0125 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:
0126
0127 #####################################################
0128 # Density Range Material #
0129 #---------------------------------------------------#
0130 # mg/cm3 - #
0131 #---------------------------------------------------#
0132 # [ 0. , 0.207 ) Air #
0133 # [ 0.207 , 0.481 ) Lungs (inhale) #
0134 # [ 0.481 , 0.919 ) Lungs (exhale) #
0135 # [ 0.919 , 0.979 ) Adipose #
0136 # [ 0.979 , 1.004 ) Breast #
0137 # [ 1.004 , 1.043 ) Phantom #
0138 # [ 1.043 , 1.109 ) Liver #
0139 # [ 1.109 , 1.113 ) Muscle #
0140 # [ 1.113 , 1.496 ) Trabecular Bone#
0141 # [ 1.496 , 1.654 ] Dense Bone #
0142 #####################################################
0143
0144 Data taken from the International Commission on Radiation Units and measurements (ICRU) report 46 was used to build the materials (lung, liver, breast, bones, ...).
0145
0146 When using the Digital Head Phantom, the CT2Density.dat is not used. The conversion is performed directly in the Dicom Handler.cc
0147
0148 ## Splitting materials in density intervals:
0149
0150 In the class DicomDetectorConstruction, it is defined a density interval
0151
0152 G4double densityDiff = 0.1;
0153
0154 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.
0155
0156 ## Voxel colouring:
0157
0158 The file Colormap.dat defines the colour that will be assigned to the voxels of each material.
0159
0160 ## DICOM file formats:
0161
0162 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):
0163
0164 - A line with the number of materials
0165 - A line for each material with its index and name (the same name of materials that you construct as G4Material's)
0166 - A line with the number of voxels in X, Y and Z
0167 - A line with the minimum and maximum extension in X (mm)
0168 - A line with the minimum and maximum extension in Y (mm)
0169 - A line with the minimum and maximum extension in Z (mm)
0170 - A number of lines containing the nVoxelX*nVoxelY*nVoxelZ material indices (one per voxel)
0171 - A number of lines containing the nVoxelX*nVoxelY*nVoxelZ material densities (one per voxel)
0172
0173 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.
0174
0175 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 .
0176
0177 ## Choosing different parameterisation/navigation options:
0178
0179 There are four possible ways in GEANT4 to treat the navigation in regular voxelised volumes:
0180
0181 - 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
0182 - Using G4NestedParameterisation. The search is done hierarchically in X, Y and Z. It is fast and does not require big memory
0183 - 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).
0184
0185 Obsolete option:
0186 - 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.
0187
0188 You can select amongst the four options in the following way:
0189
0190 - By default the example will run with G4RegularNavigation
0191
0192 - To use the first option at RegularDicomDetectorConstruction.cc you must set
0193 ```cpp
0194 patient_phys->SetRegularStructureId(0);
0195 ```
0196 - To use the second option (Nested Parameterisation) you must set the enviromental variable DICOM_NESTED_PARAM to 1
0197
0198 - To use the final, obsolete 1D-option, apart from the change above at RegularDicomDetectorConstruction.cc you need to replace (i.e. use kUndefined)
0199 ```cpp
0200 G4PVParameterised * patient_phys = new G4PVParameterised("Patient",voxel_logic,container_logic,
0201 kUndefined, nVoxelX*nVoxelY*nVoxelZ, param);
0202 ```
0203 by
0204 ```cpp
0205 G4PVParameterised * patient_phys = new G4PVParameterised("Patient",voxel_logic,container_logic,
0206 kXAxis, nVoxelX*nVoxelY*nVoxelZ, param);
0207 ```
0208 Note also you must *not* set the enviromental variable DICOM_NESTED_PARAM.
0209
0210
0211 ## Calculating dose in phantom voxels for regular navigation
0212
0213 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.
0214
0215 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.
0216 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.
0217 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.
0218
0219 This scorer is implemented in this examples if the regular navigation option is
0220 chosen. It is triggered at the method RegularDicomDetectorConstruction::ConstructPhantom() by the call
0221 ```cpp
0222 SetScorer(voxel_logic);
0223 ```
0224 ## Calculating dose in phantom voxels for nested parameterisation
0225
0226 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
0227 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.
0228
0229 ## Output
0230 dicom.out is produced running the macro file run.mac. It has 2 columns: the first is the number of
0231 voxel (ordered in x,y,z) and the second the dose there deposited (in Gy)
0232 It is produced, as an example, with a compression value of 32
0233
0234
0235 ## Partial phantom
0236 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
0237 ```
0238 /dicom/intersectWithUserVolume 0. 0. 0. 45.*deg 0. 0. TUBE 0. 150. 100.
0239 ```
0240 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
0241
0242 Alternatively you can intersect the phantom with an existing Geant4 volume with the command
0243 ```
0244 /dicom/intersectWithG4Volume VOLUME_NAME
0245 ```
0246 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
0247
0248 ## Visualisation
0249
0250 The Geant4 drivers are not meant for visualizing millions of voxel and visualising the DICOM geometries can be very computationally demanding.
0251 The users may want to visualise each DICOM slice separately or use higher compression values when visualising a part of DICOM project.
0252
0253
0254