Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1910:f81cdf283859
display RT-STRUCT before the referenced slices get loaded
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 01 Feb 2022 10:40:35 +0100 |
parents | 782ba9eb6f22 |
children | bd527bbc34df |
comparison
equal
deleted
inserted
replaced
1909:782ba9eb6f22 | 1910:f81cdf283859 |
---|---|
216 { | 216 { |
217 return false; | 217 return false; |
218 } | 218 } |
219 else | 219 else |
220 { | 220 { |
221 // return true; // TODO - TEST | |
222 | |
223 const CoordinateSystem3D& geometry = it->second.geometry_; | 221 const CoordinateSystem3D& geometry = it->second.geometry_; |
224 | 222 |
225 hasSlice_ = true; | 223 hasSlice_ = true; |
226 geometry_ = geometry; | 224 geometry_ = geometry; |
227 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); | 225 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); |
279 double d2 = GeometryToolbox::ProjectAlongNormal(points_.front(), estimatedNormal); | 277 double d2 = GeometryToolbox::ProjectAlongNormal(points_.front(), estimatedNormal); |
280 return (LinearAlgebra::IsNear(d1, d2, estimatedSliceThickness / 2.0)); | 278 return (LinearAlgebra::IsNear(d1, d2, estimatedSliceThickness / 2.0)); |
281 } | 279 } |
282 } | 280 } |
283 } | 281 } |
284 | 282 |
283 | |
284 static void AddSegmentIfIntersection(Extent2D& extent, | |
285 const CoordinateSystem3D& cuttingPlane, | |
286 const Vector& p1, | |
287 const Vector& p2, | |
288 double originDistance) | |
289 { | |
290 // Does this segment intersects the cutting plane? | |
291 double d1 = cuttingPlane.ProjectAlongNormal(p1); | |
292 double d2 = cuttingPlane.ProjectAlongNormal(p2); | |
293 | |
294 if ((d1 < originDistance && d2 > originDistance) || | |
295 (d1 > originDistance && d2 < originDistance)) | |
296 { | |
297 // This is an intersection: Add the segment | |
298 double x, y; | |
299 cuttingPlane.ProjectPoint2(x, y, p1); | |
300 extent.AddPoint(x, y); | |
301 cuttingPlane.ProjectPoint2(x, y, p2); | |
302 extent.AddPoint(x, y); | |
303 } | |
304 } | |
305 | |
285 bool DicomStructureSet::Polygon::Project(double& x1, | 306 bool DicomStructureSet::Polygon::Project(double& x1, |
286 double& y1, | 307 double& y1, |
287 double& x2, | 308 double& x2, |
288 double& y2, | 309 double& y2, |
289 const CoordinateSystem3D& cuttingPlane, | 310 const CoordinateSystem3D& cuttingPlane, |
290 const Vector& estimatedNormal, | 311 const Vector& estimatedNormal, |
291 double estimatedSliceThickness) const | 312 double estimatedSliceThickness) const |
292 { | 313 { |
293 #if 0 | |
294 if (points_.size() <= 1) | 314 if (points_.size() <= 1) |
295 { | 315 { |
296 return false; | 316 return false; |
297 } | 317 } |
298 | 318 |
301 if (hasSlice_) | 321 if (hasSlice_) |
302 { | 322 { |
303 normal = geometry_.GetNormal(); | 323 normal = geometry_.GetNormal(); |
304 thickness = sliceThickness_; | 324 thickness = sliceThickness_; |
305 } | 325 } |
306 | 326 |
307 bool isOpposite; | 327 bool isOpposite; |
308 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) && | 328 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) && |
309 !GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY())) | 329 !GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY())) |
310 { | 330 { |
311 printf("UUUU\n"); | |
312 return false; | 331 return false; |
313 } | 332 } |
314 | 333 |
315 return false; | 334 const double d = cuttingPlane.ProjectAlongNormal(cuttingPlane.GetOrigin()); |
316 | 335 |
317 #else | 336 Extent2D extent; |
318 if (!hasSlice_ || | 337 AddSegmentIfIntersection(extent, cuttingPlane, points_[points_.size() - 1], points_[0], d); |
319 points_.size() <= 1) | 338 for (size_t i = 1; i < points_.size(); i++) |
339 { | |
340 AddSegmentIfIntersection(extent, cuttingPlane, points_[i - 1], points_[i], d); | |
341 } | |
342 | |
343 if (extent.GetWidth() > 0 || | |
344 extent.GetHeight() > 0) | |
345 { | |
346 // Let's convert them to 3D world geometry to add the slice thickness | |
347 Vector p1 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX1(), extent.GetY1()) + | |
348 thickness / 2.0 * normal); | |
349 Vector p2 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX2(), extent.GetY2()) - | |
350 thickness / 2.0 * normal); | |
351 | |
352 // Then back to the cutting plane geometry | |
353 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
354 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
355 return true; | |
356 } | |
357 else | |
320 { | 358 { |
321 return false; | 359 return false; |
322 } | 360 } |
323 | |
324 double x, y; | |
325 geometry_.ProjectPoint2(x, y, cuttingPlane.GetOrigin()); | |
326 | |
327 bool isOpposite; | |
328 if (GeometryToolbox::IsParallelOrOpposite | |
329 (isOpposite, cuttingPlane.GetNormal(), geometry_.GetAxisY())) | |
330 { | |
331 // plane is constant Y | |
332 | |
333 if (y < extent_.GetY1() || | |
334 y > extent_.GetY2()) | |
335 { | |
336 // The polygon does not intersect the input slice | |
337 return false; | |
338 } | |
339 | |
340 bool isFirst = true; | |
341 double xmin = std::numeric_limits<double>::infinity(); | |
342 double xmax = -std::numeric_limits<double>::infinity(); | |
343 | |
344 double prevX, prevY; | |
345 geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); | |
346 | |
347 for (size_t i = 0; i < points_.size(); i++) | |
348 { | |
349 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py | |
350 double curX, curY; | |
351 geometry_.ProjectPoint2(curX, curY, points_[i]); | |
352 | |
353 // if prev* and cur* are on opposite sides of y, this means that the | |
354 // segment intersects the plane. | |
355 if ((prevY <= y && curY >= y) || | |
356 (prevY >= y && curY <= y)) | |
357 { | |
358 double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); | |
359 xmin = std::min(xmin, p); | |
360 xmax = std::max(xmax, p); | |
361 isFirst = false; | |
362 | |
363 // xmin and xmax represent the extent of the rectangle along the | |
364 // intersection between the plane and the polygon geometry | |
365 | |
366 } | |
367 | |
368 prevX = curX; | |
369 prevY = curY; | |
370 } | |
371 | |
372 // if NO segment intersects the plane | |
373 if (isFirst) | |
374 { | |
375 return false; | |
376 } | |
377 else | |
378 { | |
379 // y is the plane y coord in the polygon geometry | |
380 // xmin and xmax are ALSO expressed in the polygon geometry | |
381 | |
382 // let's convert them to 3D world geometry... | |
383 Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + | |
384 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
385 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - | |
386 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
387 | |
388 // then to the cutting plane geometry... | |
389 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
390 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
391 return true; | |
392 } | |
393 } | |
394 else if (GeometryToolbox::IsParallelOrOpposite | |
395 (isOpposite, cuttingPlane.GetNormal(), geometry_.GetAxisX())) | |
396 { | |
397 // plane is constant X => Sagittal view (remember that in the | |
398 // sagittal projection, the normal must be swapped) | |
399 | |
400 | |
401 /* | |
402 Please read the comments in the section above, by taking into account | |
403 the fact that, in this case, the plane has a constant X, not Y (in | |
404 polygon geometry_ coordinates) | |
405 */ | |
406 | |
407 if (x < extent_.GetX1() || | |
408 x > extent_.GetX2()) | |
409 { | |
410 return false; | |
411 } | |
412 | |
413 bool isFirst = true; | |
414 double ymin = std::numeric_limits<double>::infinity(); | |
415 double ymax = -std::numeric_limits<double>::infinity(); | |
416 | |
417 double prevX, prevY; | |
418 geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); | |
419 | |
420 for (size_t i = 0; i < points_.size(); i++) | |
421 { | |
422 // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py | |
423 double curX, curY; | |
424 geometry_.ProjectPoint2(curX, curY, points_[i]); | |
425 | |
426 if ((prevX <= x && curX >= x) || | |
427 (prevX >= x && curX <= x)) | |
428 { | |
429 double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); | |
430 ymin = std::min(ymin, p); | |
431 ymax = std::max(ymax, p); | |
432 isFirst = false; | |
433 } | |
434 | |
435 prevX = curX; | |
436 prevY = curY; | |
437 } | |
438 | |
439 if (isFirst) | |
440 { | |
441 return false; | |
442 } | |
443 else | |
444 { | |
445 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + | |
446 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
447 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - | |
448 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
449 | |
450 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
451 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
452 | |
453 return true; | |
454 } | |
455 } | |
456 else | |
457 { | |
458 // Should not happen | |
459 return false; | |
460 } | |
461 #endif | |
462 } | 361 } |
463 | 362 |
464 | 363 |
465 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const | 364 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const |
466 { | 365 { |