ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomFileCT.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DicomFileCT.cc
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......
42 DicomFileCT::DicomFileCT()
43 {
44 }
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 DicomFileCT::DicomFileCT(DcmDataset* dset) : DicomVFileImage(dset)
48 {
49 }
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 void DicomFileCT::BuildMaterials()
53 {
54  G4int fCompress = theFileMgr->GetCompression();
55  if( fNoVoxelX%fCompress != 0 || fNoVoxelY%fCompress != 0 ) {
56  G4Exception("DicompFileMgr.:BuildMaterials",
57  "DFC004",
58  FatalException,
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;
83  if( theFileMgr->IsMaterialsDensity() ) {
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......
130 void DicomFileCT::BuildStructureIDs()
131 {
132  G4int fCompress = theFileMgr->GetCompression();
133  std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
134  if( dfs.size() == 0 ) return;
135 
136  G4int NMAXROI = DicomFileMgr::GetInstance()->GetStructureNMaxROI();
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));
201  if( DicomFileMgr::verbose >= debugVerb )
202  G4cout << " minXc " << minXc << " < " << fMinX
203  << " maxXc " << maxXc << " > " << fMaxX
204  << " minYc " << minYc << " < " << fMinY
205  << " maxYc " << maxYc << " > " << fMaxY << G4endl;
206  if( DicomFileMgr::verbose >= debugVerb )
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",
299  FatalException,
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;
315  if( DicomFileMgr::verbose >= testVerb ){
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......
377 void DicomFileCT::DumpStructureIDsToTextFile(std::ofstream& fout)
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