131{
134 if( dfs.size() == 0 ) return;
135
137 G4int NMAXROI_real = std::log(INT_MAX)/std::log(NMAXROI);
138
139
140 for(
int ir = 0; ir <
fNoVoxelsY/fCompress; ir++ ) {
141 for(
int ic = 0; ic <
fNoVoxelsX/fCompress; ic++ ) {
142
144 }
145 }
146
147 std::set<double> distInters;
148
149
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++ ){
155 << rois[jj]->GetNumber() << " " << rois[jj]->GetName() << G4endl;
156 std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
157
158 G4int roiID = rois[jj]->GetNumber();
159 for( size_t kk = 0; kk < contours.size(); kk++ ){
160
162
163
165 <<
" INTERS CONTOUR " << roic->
GetZ() <<
" SLICE Z " <<
fMinZ <<
" " <<
fMaxZ << G4endl;
166
167
172
173 std::vector<G4ThreeVector> points = roic->
GetPoints();
175 << points.size() << G4endl;
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 }
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
212
213 for( int ix = idMinX; ix <= idMaxX; ix++ ) {
214 for( int iy = idMinY; iy <= idMaxY; iy++ ) {
215 G4bool bOK = false;
216 G4int bOKs;
220
221 for( int icx = 0; icx <= 1; icx++ ){
222 for( int icy = 0; icy <= 1; icy++ ){
223 bOKs = 0;
224 if( bOK ) continue;
227 double v0x = 1.;
228 if( icx == 1 ) v0x = -1.;
231 << iy << " + " << icy << " CORNER (" << p0x << "," << p0y << ") "
232 << " DIR= (" << v0x << "," << v0y << ") " << G4endl;
234 for( int ip = 0; ip < NTRIES; ip++) {
235 distInters.clear();
236 v0y -= ip*0.1*v0y;
237 G4double halfDiagonal = 0.5*
fVoxelDimX*fCompress;
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 ) {
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
255 double lambdap = fac3/fac1;
256 if( lambdap > 0. ) {
257 distInters.insert(lambdap);
259 <<lambdaq << " (" << d0x+p0x+lambdaq*w0x << "," << d0y+p0y+lambdaq*w0y
260 << ") = (" << p0x+lambdap*v0x << "," << p0y+lambdap*v0y << ") "
261 << " N " << distInters.size() << G4endl;
262 }
264 << lambdaq << " P " << lambdap << G4endl;
265
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
273 bOKs += 1;
275 << " N INTERS " << distInters.size() << " " << *(distInters.begin())
276 << " > " << halfDiagonal << G4endl;
277 }
278 }
279 }
280 if(bOKs == NTRIES ) {
281 bOK = true;
285 }
286
287 }
288 }
289 if( bOK ) {
290
292
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 ) {
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 }
316 G4cout << " WITH PREVIOUS ROI IN VOXEL " << roival << G4endl;
317 }
318 } else {
320 }
321 }
322
323 }
324 }
325 }
326 }
327 }
328 }
329
331
332
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++ ){
341
345 << kk <<
" INTERS CONTOUR " << roic->
GetZ() <<
" SLICE Z "
348 std::vector<G4ThreeVector> points = roic->
GetPoints();
351 << std::endl;
352 for( size_t ll = 0; ll < points.size(); 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
362 for(
int ir = 0; ir <
fNoVoxelsY/fCompress; ir++ ) {
363 for(
int ic = 0; ic <
fNoVoxelsX/fCompress; ic++ ) {
366 <<
" " << ir <<
" STRUCTURE VOXEL (" <<
fMinX +
fVoxelDimX*fCompress * (ic+0.5)
369 }
370 }
371 }
372 }
373
374}
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