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