Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | List of all members
DicomFileCT Class Reference

#include <Doxymodules_medical.h>

Inheritance diagram for DicomFileCT:
DicomVFileImage DicomVFile

Public Member Functions

 DicomFileCT ()
 
 DicomFileCT (DcmDataset *dset)
 
 ~DicomFileCT ()
 
void BuildMaterials ()
 
void DumpMateIDsToTextFile (std::ofstream &fout)
 
void DumpDensitiesToTextFile (std::ofstream &fout)
 
void BuildStructureIDs ()
 
void DumpStructureIDsToTextFile (std::ofstream &fout)
 
- Public Member Functions inherited from DicomVFileImage
 DicomVFileImage ()
 
 DicomVFileImage (DcmDataset *dset)
 
 ~DicomVFileImage ()
 
virtual void ReadData ()
 
void operator+= (const DicomVFileImage &rhs)
 
DicomVFileImage operator+ (const DicomVFileImage &rhs)
 
void DumpHeaderToTextFile (std::ofstream &fout)
 
G4int GetNoVoxelsX () const
 
G4int GetNoVoxelsY () const
 
G4int GetNoVoxelsZ () const
 
G4int GetNoVoxels () const
 
G4double GetMinX () const
 
G4double GetMinY () const
 
G4double GetMinZ () const
 
G4double GetMaxX () const
 
G4double GetMaxY () const
 
G4double GetMaxZ () const
 
void SetNoVoxelsX (const G4int &val)
 
void SetNoVoxelsY (const G4int &val)
 
void SetNoVoxelsZ (const G4int &val)
 
void SetMinX (const G4double &val)
 
void SetMaxX (const G4double &val)
 
void SetMinY (const G4double &val)
 
void SetMaxY (const G4double &val)
 
void SetMinZ (const G4double &val)
 
void SetMaxZ (const G4double &val)
 
const G4double & GetLocation () const
 
void SetLocation (const G4double &val)
 
G4ThreeVector GetOrientationRows () const
 
G4ThreeVector GetOrientationColumns () const
 
- Public Member Functions inherited from DicomVFile
 DicomVFile ()
 
 DicomVFile (DcmDataset *dset)
 
 ~DicomVFile ()
 
void SetFileName (G4String fName)
 

Private Attributes

std::vector< size_t > fMateIDs
 
std::vector< G4double > fDensities
 
std::vector< G4int > fStructure
 

Additional Inherited Members

- Protected Member Functions inherited from DicomVFile
virtual std::vector< G4double > Read1Data (DcmDataset *dset, DcmTagKey tagKey, G4int nData)
 
virtual OFString Read1DataStr (DcmDataset *dset, DcmTagKey tagKey)
 
- Protected Attributes inherited from DicomVFileImage
G4double fLocation
 
G4double fBitAllocated
 
G4double fRescaleSlope
 
G4double fRescaleIntercept
 
G4int fNoVoxelsX
 
G4int fNoVoxelsY
 
G4int fNoVoxelsZ
 
G4double fMinX
 
G4double fMinY
 
G4double fMinZ
 
G4double fMaxX
 
G4double fMaxY
 
G4double fMaxZ
 
G4double fVoxelDimX
 
G4double fVoxelDimY
 
G4double fVoxelDimZ
 
G4ThreeVector fOrientationRows
 
G4ThreeVector fOrientationColumns
 
std::vector< int > fHounsfieldV
 
DicomFileMgrtheFileMgr
 
- Protected Attributes inherited from DicomVFile
DcmDataset * theDataset
 
G4String fFileName
 

Detailed Description

Definition at line 48 of file Doxymodules_medical.h.

Constructor & Destructor Documentation

◆ DicomFileCT() [1/2]

DicomFileCT::DicomFileCT ( )

Definition at line 42 of file DicomFileCT.cc.

43{
44}

◆ DicomFileCT() [2/2]

DicomFileCT::DicomFileCT ( DcmDataset *  dset)

Definition at line 47 of file DicomFileCT.cc.

47 : DicomVFileImage(dset)
48{
49}

◆ ~DicomFileCT()

DicomFileCT::~DicomFileCT ( )
inline

Definition at line 36 of file DicomFileCT.hh.

36{};

Member Function Documentation

◆ BuildMaterials()

void DicomFileCT::BuildMaterials ( )

Definition at line 52 of file DicomFileCT.cc.

53{
54 G4int fCompress = theFileMgr->GetCompression();
55 if( fNoVoxelsX%fCompress != 0 || fNoVoxelsY%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(fNoVoxelsX)
61 + " and Y " + std::to_string(fNoVoxelsY)).c_str());
62 }
63
64 // if( DicomVerb(debugVerb) ) G4cout << " BuildMaterials " << fFileName << G4endl;
65 double meanHV = 0.;
66 for( int ir = 0; ir < fNoVoxelsY; ir += fCompress ) {
67 for( int ic = 0; ic < fNoVoxelsX; ic += fCompress ) {
68 meanHV = 0.;
69 int isumrMax = std::min(ir+fCompress,fNoVoxelsY);
70 int isumcMax = std::min(ic+fCompress,fNoVoxelsX);
71 for( int isumr = ir; isumr < isumrMax; isumr ++ ) {
72 for( int isumc = ic; isumc < isumcMax; isumc ++ ) {
73 meanHV += fHounsfieldV[isumc+isumr*fNoVoxelsX];
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}
std::vector< size_t > fMateIDs
std::vector< G4double > fDensities
size_t GetMaterialIndex(G4double Hval)
G4int GetCompression() const
G4double Hounsfield2density(Uint32 Hval)
G4bool IsMaterialsDensity() const
size_t GetMaterialIndexByDensity(G4double density)
DicomFileMgr * theFileMgr
std::vector< int > fHounsfieldV

◆ DumpMateIDsToTextFile()

void DicomFileCT::DumpMateIDsToTextFile ( std::ofstream &  fout)

Definition at line 95 of file DicomFileCT.cc.

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 < fNoVoxelsY/fCompress; ir++ ) {
101 for( int ic = 0; ic < fNoVoxelsX/fCompress; ic++ ) {
102 fout << fMateIDs[ic+ir*fNoVoxelsX/fCompress];
103 if( ic != fNoVoxelsX/fCompress-1) fout << " ";
104 }
105 fout << G4endl;
106 }
107}
@ warningVerb
static G4int verbose
G4String fFileName
Definition DicomVFile.hh:55

◆ DumpDensitiesToTextFile()

void DicomFileCT::DumpDensitiesToTextFile ( std::ofstream &  fout)

Definition at line 110 of file DicomFileCT.cc.

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 < fNoVoxelsY/fCompress; ir++ ) {
118 for( int ic = 0; ic < fNoVoxelsX/fCompress; ic++ ) {
119 fout << fDensities[ic+ir*fNoVoxelsX/fCompress];
120 if( ic != fNoVoxelsX/fCompress-1) fout << " ";
121 if( copyNo%8 == 7 ) fout << G4endl;
122 copyNo++;
123 }
124 if( copyNo%8 != 0 ) fout << G4endl;
125 }
126
127}

◆ BuildStructureIDs()

void DicomFileCT::BuildStructureIDs ( )

Definition at line 130 of file DicomFileCT.cc.

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(fNoVoxelsX*fNoVoxelsY);
140 for( int ir = 0; ir < fNoVoxelsY/fCompress; ir++ ) {
141 for( int ic = 0; ic < fNoVoxelsX/fCompress; ic++ ) {
142 // fStructure[ic+ir*fNoVoxelsX] = 0;
143 fStructure.push_back(0);
144 }
145 }
146
147 std::set<double> distInters;
148
149 // std::fill_n(fStructure,fNoVoxelsX*fNoVoxelsY,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(fNoVoxelsX/fCompress-1,int((maxXc-fMinX)/fVoxelDimX/fCompress+1));
199 int idMinY = std::max(0,int((minYc-fMinY)/fVoxelDimY/fCompress));
200 int idMaxY = std::min(fNoVoxelsY/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*fNoVoxelsX/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\
302GING -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*fNoVoxelsX/fCompress] = roival;
316 G4cout << " WITH PREVIOUS ROI IN VOXEL " << roival << G4endl;
317 }
318 } else {
319 fStructure[ix+iy*fNoVoxelsX/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 < fNoVoxelsY/fCompress; ir++ ) {
363 for( int ic = 0; ic < fNoVoxelsX/fCompress; ic++ ) {
364 if( fStructure[ic+ir*fNoVoxelsX/fCompress] != 0 ) {
365 if( DicomFileMgr::verbose >= 0 ) G4cout << ic+ir*fNoVoxelsX/fCompress << " = " << ic
366 << " " << ir << " STRUCTURE VOXEL (" << fMinX + fVoxelDimX*fCompress * (ic+0.5)
367 << "," << fMinY + fVoxelDimY*fCompress * (ir+0.5) << ") = "
368 << fStructure[ic+ir*fNoVoxelsX/fCompress] << G4endl;
369 }
370 }
371 }
372 }
373
374}
@ testVerb
@ debugVerb
std::vector< G4int > fStructure
std::vector< DicomFileStructure * > GetStructFiles() const
G4int GetStructureNMaxROI() const
G4int GetStructureNCheck() const
static DicomFileMgr * GetInstance()
std::vector< G4ThreeVector > GetDirections()
std::vector< G4ThreeVector > GetPoints() const
OFString GetGeomType() const

◆ DumpStructureIDsToTextFile()

void DicomFileCT::DumpStructureIDsToTextFile ( std::ofstream &  fout)

Definition at line 377 of file DicomFileCT.cc.

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 < fNoVoxelsY/fCompress; ir++ ) {
386 for( int ic = 0; ic < fNoVoxelsX/fCompress; ic++ ) {
387 fout << fStructure[ic+ir*fNoVoxelsX/fCompress];
388 if( ic != fNoVoxelsX/fCompress-1) fout << " ";
389 }
390 fout << G4endl;
391 }
392}

Member Data Documentation

◆ fMateIDs

std::vector<size_t> DicomFileCT::fMateIDs
private

Definition at line 46 of file DicomFileCT.hh.

◆ fDensities

std::vector<G4double> DicomFileCT::fDensities
private

Definition at line 47 of file DicomFileCT.hh.

◆ fStructure

std::vector<G4int> DicomFileCT::fStructure
private

Definition at line 48 of file DicomFileCT.hh.


The documentation for this class was generated from the following files:

Applications | User Support | Publications | Collaboration