Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomFileCT.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 #include "DicomFileCT.hh"
27 #include "DicomFileStructure.hh"
28 #include "DicomROI.hh"
29 
30 #include "G4GeometryTolerance.hh"
31 
32 #include "dcmtk/dcmdata/dcfilefo.h"
33 #include "dcmtk/dcmdata/dcdeftag.h"
34 #include "dcmtk/dcmdata/dcpixel.h"
35 #include "dcmtk/dcmdata/dcpxitem.h"
36 #include "dcmtk/dcmdata/dcpixseq.h"
37 #include "dcmtk/dcmrt/drtimage.h"
38 
39 #include <set>
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 {
44 }
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 DicomFileCT::DicomFileCT(DcmDataset* dset) : DicomVFileImage(dset)
48 {
49 }
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 {
54  G4int fCompress = theFileMgr->GetCompression();
55  if( fNoVoxelX%fCompress != 0 || fNoVoxelY%fCompress != 0 ) {
56  G4Exception("DicompFileMgr.:BuildMaterials",
57  "DFC004",
59  ("Compression factor = " + std::to_string(fCompress)
60  + " has to be a divisor of Number of voxels X = " + std::to_string(fNoVoxelX)
61  + " and Y " + std::to_string(fNoVoxelY)).c_str());
62  }
63 
64  // if( DicomVerb(debugVerb) ) G4cout << " BuildMaterials " << fFileName << G4endl;
65  double meanHV = 0.;
66  for( int ir = 0; ir < fNoVoxelY; ir += fCompress ) {
67  for( int ic = 0; ic < fNoVoxelX; ic += fCompress ) {
68  meanHV = 0.;
69  int isumrMax = std::min(ir+fCompress,fNoVoxelY);
70  int isumcMax = std::min(ic+fCompress,fNoVoxelX);
71  for( int isumr = ir; isumr < isumrMax; isumr ++ ) {
72  for( int isumc = ic; isumc < isumcMax; isumc ++ ) {
73  meanHV += fHounsfieldV[isumc+isumr*fNoVoxelX];
74  // G4cout << isumr << " " << isumc << " GET mean " << meanHV << G4endl;
75  }
76  }
77  meanHV /= (isumrMax-ir)*(isumcMax-ic);
78  G4double meanDens = theFileMgr->Hounsfield2density(meanHV);
79  // G4cout << ir << " " << ic << " FINAL mean " << meanDens << G4endl;
80 
81  fDensities.push_back(meanDens);
82  size_t mateID;
84  mateID = theFileMgr->GetMaterialIndexByDensity(meanDens);
85  } else {
86  mateID = theFileMgr->GetMaterialIndex(meanHV);
87  }
88  fMateIDs.push_back(mateID);
89  }
90  }
91 
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 void DicomFileCT::DumpMateIDsToTextFile(std::ofstream& fout)
96 {
97  G4int fCompress = theFileMgr->GetCompression();
98  if( DicomFileMgr::verbose >= warningVerb ) G4cout << fLocation << " DumpMateIDsToTextFile "
99  << fFileName << " " << fMateIDs.size() << G4endl;
100  for( int ir = 0; ir < fNoVoxelY/fCompress; ir++ ) {
101  for( int ic = 0; ic < fNoVoxelX/fCompress; ic++ ) {
102  fout << fMateIDs[ic+ir*fNoVoxelX/fCompress];
103  if( ic != fNoVoxelX/fCompress-1) fout << " ";
104  }
105  fout << G4endl;
106  }
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 void DicomFileCT::DumpDensitiesToTextFile(std::ofstream& fout)
111 {
112  G4int fCompress = theFileMgr->GetCompression();
113  if( DicomFileMgr::verbose >= warningVerb ) G4cout << fLocation << " DumpDensitiesToTextFile "
114  << fFileName << " " << fDensities.size() << G4endl;
115 
116  G4int copyNo = 0;
117  for( int ir = 0; ir < fNoVoxelY/fCompress; ir++ ) {
118  for( int ic = 0; ic < fNoVoxelX/fCompress; ic++ ) {
119  fout << fDensities[ic+ir*fNoVoxelX/fCompress];
120  if( ic != fNoVoxelX/fCompress-1) fout << " ";
121  if( copyNo%8 == 7 ) fout << G4endl;
122  copyNo++;
123  }
124  if( copyNo%8 != 0 ) fout << G4endl;
125  }
126 
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 {
132  G4int fCompress = theFileMgr->GetCompression();
133  std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
134  if( dfs.size() == 0 ) return;
135 
137  G4int NMAXROI_real = std::log(INT_MAX)/std::log(NMAXROI);
138 
139  // fStructure = new G4int(fNoVoxelX*fNoVoxelY);
140  for( int ir = 0; ir < fNoVoxelY/fCompress; ir++ ) {
141  for( int ic = 0; ic < fNoVoxelX/fCompress; ic++ ) {
142  // fStructure[ic+ir*fNoVoxelX] = 0;
143  fStructure.push_back(0);
144  }
145  }
146 
147  std::set<double> distInters;
148 
149  // std::fill_n(fStructure,fNoVoxelX*fNoVoxelY,0);
150  //
151  for( size_t ii = 0; ii < dfs.size(); ii++ ){
152  std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
153  for( size_t jj = 0; jj < rois.size(); jj++ ){
154  if( DicomFileMgr::verbose >= debugVerb ) std::cout << " BuildStructureIDs checking ROI "
155  << rois[jj]->GetNumber() << " " << rois[jj]->GetName() << G4endl;
156  std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
157  // G4int roi = jj+1;
158  G4int roiID = rois[jj]->GetNumber();
159  for( size_t kk = 0; kk < contours.size(); kk++ ){
160  // check contour corresponds to this CT slice
161  DicomROIContour* roic = contours[kk];
162  // if( DicomVerb(-debugVerb) ) G4cout << jj << " " << kk << " INTERS CONTOUR " << roic
163  // << " " << fLocation << G4endl;
164  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " " << roiID << " " << kk
165  << " INTERS CONTOUR " << roic->GetZ() << " SLICE Z " << fMinZ << " " << fMaxZ << G4endl;
166  // Check Contour correspond to slice
167 
168  if( roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ ) continue;
169  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " INTERS CONTOUR WITH Z SLIZE "
170  << fMinZ << " < " << roic->GetZ() << " < " << fMaxZ << G4endl;
171  if( roic->GetGeomType() == "CLOSED_PLANAR" ){
172  // Get min and max X and Y, and corresponding slice indexes
173  std::vector<G4ThreeVector> points = roic->GetPoints();
174  if( DicomFileMgr::verbose >= debugVerb ) G4cout << jj << " " << kk << " NPOINTS "
175  << points.size() << G4endl;
176  std::vector<G4ThreeVector> dirs = roic->GetDirections();
177  double minXc = DBL_MAX;
178  double maxXc = -DBL_MAX;
179  double minYc = DBL_MAX;
180  double maxYc = -DBL_MAX;
181  for( size_t ll = 0; ll < points.size(); ll++ ){
182  minXc = std::min(minXc,points[ll].x());
183  maxXc = std::max(maxXc,points[ll].x());
184  minYc = std::min(minYc,points[ll].y());
185  maxYc = std::max(maxYc,points[ll].y());
186  }
187  if( minXc < fMinX || maxXc > fMaxX || minYc < fMinY || maxYc > fMaxY ){
188  G4cerr << " minXc " << minXc << " < " << fMinX
189  << " maxXc " << maxXc << " > " << fMaxX
190  << " minYc " << minYc << " < " << fMinY
191  << " maxYc " << maxYc << " > " << fMaxY << G4endl;
192  G4Exception("DicomFileCT::BuildStructureIDs",
193  "DFCT001",
194  JustWarning,
195  "Contour limits exceed Z slice extent");
196  }
197  int idMinX = std::max(0,int((minXc-fMinX)/fVoxelDimX/fCompress));
198  int idMaxX = std::min(fNoVoxelX/fCompress-1,int((maxXc-fMinX)/fVoxelDimX/fCompress+1));
199  int idMinY = std::max(0,int((minYc-fMinY)/fVoxelDimY/fCompress));
200  int idMaxY = std::min(fNoVoxelY/fCompress-1,int((maxYc-fMinY)/fVoxelDimY/fCompress+1));
202  G4cout << " minXc " << minXc << " < " << fMinX
203  << " maxXc " << maxXc << " > " << fMaxX
204  << " minYc " << minYc << " < " << fMinY
205  << " maxYc " << maxYc << " > " << fMaxY << G4endl;
207  G4cout << " idMinX " << idMinX
208  << " idMaxX " << idMaxX
209  << " idMinY " << idMinY
210  << " idMaxY " << idMaxY << G4endl;
211  //for each voxel: build 4 lines from the corner towards the center
212  // and check how many contour segments it crosses, and the minimum distance to a segment
213  for( int ix = idMinX; ix <= idMaxX; ix++ ) {
214  for( int iy = idMinY; iy <= idMaxY; iy++ ) {
215  G4bool bOK = false;
216  G4int bOKs;
217  if( DicomFileMgr::verbose >= debugVerb ) G4cout << "@@@@@ TRYING POINT ("
218  << fMinX + fVoxelDimX*fCompress*(ix+0.5) << ","
219  << fMinY + fVoxelDimY*fCompress*(iy+0.5) << ")" << G4endl;
220  // four corners
221  for( int icx = 0; icx <= 1; icx++ ){
222  for( int icy = 0; icy <= 1; icy++ ){
223  bOKs = 0;
224  if( bOK ) continue;
225  double p0x = fMinX + fVoxelDimX*fCompress * (ix+icx);
226  double p0y = fMinY + fVoxelDimY*fCompress * (iy+icy);
227  double v0x = 1.;
228  if( icx == 1 ) v0x = -1.;
229  double v0y = 0.99*fVoxelDimY/fVoxelDimX*std::pow(-1.,icy);
230  if( DicomFileMgr::verbose >= testVerb ) G4cout << ix << " + " << icx << " "
231  << iy << " + " << icy << " CORNER (" << p0x << "," << p0y << ") "
232  << " DIR= (" << v0x << "," << v0y << ") " << G4endl;
233  int NTRIES = theFileMgr->GetStructureNCheck();
234  for( int ip = 0; ip < NTRIES; ip++) {
235  distInters.clear();
236  v0y -= ip*0.1*v0y;
237  G4double halfDiagonal = 0.5*fVoxelDimX*fCompress;
238  if( DicomFileMgr::verbose >= testVerb ) G4cout << ip
239  << " TRYING WITH DIRECTION (" << " DIR= (" << v0x << ","
240  << v0y << ") " << bOKs << G4endl;
241  for( size_t ll = 0; ll < points.size(); ll++ ){
242  double d0x = points[ll].x() - p0x;
243  double d0y = points[ll].y() - p0y;
244  double w0x = dirs[ll].x();
245  double w0y = dirs[ll].y();
246  double fac1 = w0x*v0y - w0y*v0x;
247  if( fac1 == 0 ) { // parallel lines
248  continue;
249  }
250  double fac2 = d0x*v0y - d0y*v0x;
251  double fac3 = d0y*w0x - d0x*w0y;
252  double lambdaq = -fac2/fac1;
253  if( lambdaq < 0. || lambdaq >= 1. ) continue;
254  // intersection further than segment length
255  double lambdap = fac3/fac1;
256  if( lambdap > 0. ) {
257  distInters.insert(lambdap);
258  if( DicomFileMgr::verbose >= testVerb ) G4cout << " !! GOOD INTERS "
259  <<lambdaq << " (" << d0x+p0x+lambdaq*w0x << "," << d0y+p0y+lambdaq*w0y
260  << ") = (" << p0x+lambdap*v0x << "," << p0y+lambdap*v0y << ") "
261  << " N " << distInters.size() << G4endl;
262  }
263  if( DicomFileMgr::verbose >= testVerb ) G4cout << " INTERS LAMBDAQ "
264  << lambdaq << " P " << lambdap << G4endl;
265 
266  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " INTERS POINT ("
267  << d0x+p0x+lambdaq*w0x << "," << d0y+p0y+lambdaq*w0y << ") = ("
268  << p0x+lambdap*v0x << "," << p0y+lambdap*v0y << ") " << G4endl;
269  }
270  if( distInters.size() % 2 == 1 ) {
271  if( *(distInters.begin() ) > halfDiagonal ) {
272  // bOK = true;
273  bOKs += 1;
274  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " OK= " << bOKs
275  << " N INTERS " << distInters.size() << " " << *(distInters.begin())
276  << " > " << halfDiagonal << G4endl;
277  }
278  }
279  }
280  if(bOKs == NTRIES ) {
281  bOK = true;
282  if( DicomFileMgr::verbose >= debugVerb ) G4cout << "@@@@@ POINT OK ("
283  << fMinX + fVoxelDimX*fCompress*(ix+0.5) << ","
284  << fMinY + fVoxelDimY*fCompress*(iy+0.5) << ")" << G4endl;
285  }
286 
287  }
288  } // loop to four corners
289  if( bOK ) {
290  // extract previous ROI value
291  int roival = fStructure[ix+iy*fNoVoxelX/fCompress];
292  // roival = 2 + NMAXROI*3 + NMAXROI*NMAXROI*15;
293  if(roival != 0 && roival != int(roiID) ) {
294  std::set<G4int> roisVoxel;
295  int nRois = std::log10(roival)/std::log10(NMAXROI)+1;
296  if( nRois > NMAXROI_real ) { // another one cannot be added
297  G4Exception("DicomFileCT:BuildStructureIDs",
298  "DFC0004",
300  G4String("Too many ROIs associated to a voxel: \
301 " + std::to_string(nRois) + " > " + std::to_string(NMAXROI_real) + " TRY CHAN\
302 GING -NStructureNMaxROI argument to a lower value").c_str());
303  }
304  for( int inr = 0; inr < nRois; inr++ ) {
305  roisVoxel.insert( int(roival/std::pow(NMAXROI,inr))%NMAXROI );
306  }
307  roisVoxel.insert(roiID);
308  roival = 0;
309  size_t inr = 0;
310  for( std::set<G4int>::const_iterator ite = roisVoxel.begin();
311  ite != roisVoxel.end(); ite++, inr++ ) {
312  roival += (*ite)*std::pow(NMAXROI,inr);
313  }
314  fStructure[ix+iy*fNoVoxelX/fCompress] = roival;
316  G4cout << " WITH PREVIOUS ROI IN VOXEL " << roival << G4endl;
317  }
318  } else {
319  fStructure[ix+iy*fNoVoxelX/fCompress] = roiID;
320  }
321  }
322 
323  }
324  }
325  }
326  }
327  }
328  }
329 
330  if( DicomFileMgr::verbose >= 1 ) {
331  //@@@@ PRINT structures
332  //@@@ PRINT points of structures in this Z slice
333  if( DicomFileMgr::verbose >= 0 ) G4cout << " STRUCTURES FOR SLICE " << fLocation << G4endl;
334  for( size_t ii = 0; ii < dfs.size(); ii++ ){
335  std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
336  for( size_t jj = 0; jj < rois.size(); jj++ ){
337  int roi = rois[jj]->GetNumber();
338  std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
339  for( size_t kk = 0; kk < contours.size(); kk++ ){
340  DicomROIContour* roic = contours[kk];
341  // check contour corresponds to this CT slice
342  if( roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ ) continue;
343  if( roic->GetGeomType() == "CLOSED_PLANAR" ){
344  if( DicomFileMgr::verbose >= testVerb ) G4cout << " PRINTING CONTOUR? " << roi << " "
345  << kk << " INTERS CONTOUR " << roic->GetZ() << " SLICE Z "
346  << fMinZ << " " << fMaxZ << G4endl;
347  if( roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ ) continue;
348  std::vector<G4ThreeVector> points = roic->GetPoints();
349  std::vector<G4ThreeVector> dirs = roic->GetDirections();
350  if( DicomFileMgr::verbose >= 0 ) std::cout << " STRUCTURE Z " << roic->GetZ()
351  << std::endl;
352  for( size_t ll = 0; ll < points.size(); ll++ ){
353  if( DicomFileMgr::verbose >= 0 ) G4cout << roi << " : " << ll
354  << " STRUCTURE POINT (" << points[ll].x() << "," << points[ll].y() << ") "
355  << " (" << dirs[ll].x() << "," << dirs[ll].y() << ") = " << roi << G4endl;
356  }
357  }
358  }
359  }
360  }
361  //@@@ PRINT points in slice inside structure
362  for( int ir = 0; ir < fNoVoxelY/fCompress; ir++ ) {
363  for( int ic = 0; ic < fNoVoxelX/fCompress; ic++ ) {
364  if( fStructure[ic+ir*fNoVoxelX/fCompress] != 0 ) {
365  if( DicomFileMgr::verbose >= 0 ) G4cout << ic+ir*fNoVoxelX/fCompress << " = " << ic
366  << " " << ir << " STRUCTURE VOXEL (" << fMinX + fVoxelDimX*fCompress * (ic+0.5)
367  << "," << fMinY + fVoxelDimY*fCompress * (ir+0.5) << ") = "
368  << fStructure[ic+ir*fNoVoxelX/fCompress] << G4endl;
369  }
370  }
371  }
372  }
373 
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
378 {
379  G4int fCompress = theFileMgr->GetCompression();
380  if( DicomFileMgr::verbose >= 0 ) G4cout << fLocation << " DumpStructureIDsToTextFile "
381  << fFileName << " " << fStructure.size() << G4endl;
382  std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
383  if( dfs.size() == 0 ) return;
384 
385  for( int ir = 0; ir < fNoVoxelY/fCompress; ir++ ) {
386  for( int ic = 0; ic < fNoVoxelX/fCompress; ic++ ) {
387  fout << fStructure[ic+ir*fNoVoxelX/fCompress];
388  if( ic != fNoVoxelX/fCompress-1) fout << " ";
389  }
390  fout << G4endl;
391  }
392 }
393 
G4bool IsMaterialsDensity() const
Definition: DicomFileMgr.hh:99
void DumpStructureIDsToTextFile(std::ofstream &fout)
Definition: DicomFileCT.cc:377
G4double Hounsfield2density(Uint32 Hval)
void DumpMateIDsToTextFile(std::ofstream &fout)
Definition: DicomFileCT.cc:95
G4int GetCompression() const
Definition: DicomFileMgr.hh:91
tuple x
Definition: test.py:50
std::vector< DicomFileStructure * > GetStructFiles() const
Definition: DicomFileMgr.hh:56
int G4int
Definition: G4Types.hh:78
G4String fFileName
Definition: DicomVFile.hh:55
void BuildMaterials()
Definition: DicomFileCT.cc:52
void BuildStructureIDs()
Definition: DicomFileCT.cc:130
G4GLOB_DLL std::ostream G4cout
std::vector< G4ThreeVector > GetDirections()
bool G4bool
Definition: G4Types.hh:79
size_t GetMaterialIndexByDensity(G4double density)
std::vector< int > fHounsfieldV
static int verbose
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define INT_MAX
Definition: templates.hh:111
DicomFileMgr * theFileMgr
OFString GetGeomType() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4ThreeVector > GetPoints() const
size_t GetMaterialIndex(G4double Hval)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
G4int GetStructureNCheck() const
Definition: DicomFileMgr.hh:82
void DumpDensitiesToTextFile(std::ofstream &fout)
Definition: DicomFileCT.cc:110
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4int GetStructureNMaxROI() const
Definition: DicomFileMgr.hh:88
static DicomFileMgr * GetInstance()
Definition: DicomFileMgr.cc:41
G4GLOB_DLL std::ostream G4cerr