Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 2157:917e40af6b45 default
fix rendering of RT-STRUCT
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Fri, 27 Sep 2024 22:29:51 +0200 |
parents | f8cbfd7175c3 |
children |
comparison
equal
deleted
inserted
replaced
2156:340bde744884 | 2157:917e40af6b45 |
---|---|
281 } | 281 } |
282 } | 282 } |
283 } | 283 } |
284 | 284 |
285 | 285 |
286 static void AddSegmentIfIntersection(Extent2D& extent, | 286 void DicomStructureSet::Polygon::Project(std::list<Extent2D>& target, |
287 const CoordinateSystem3D& cuttingPlane, | |
288 const Vector& p1, | |
289 const Vector& p2, | |
290 double originDistance) | |
291 { | |
292 // Does this segment intersects the cutting plane? | |
293 double d1 = cuttingPlane.ProjectAlongNormal(p1); | |
294 double d2 = cuttingPlane.ProjectAlongNormal(p2); | |
295 | |
296 if ((d1 < originDistance && d2 > originDistance) || | |
297 (d1 > originDistance && d2 < originDistance)) | |
298 { | |
299 // This is an intersection: Add the segment | |
300 double x, y; | |
301 cuttingPlane.ProjectPoint2(x, y, p1); | |
302 extent.AddPoint(x, y); | |
303 cuttingPlane.ProjectPoint2(x, y, p2); | |
304 extent.AddPoint(x, y); | |
305 } | |
306 } | |
307 | |
308 bool DicomStructureSet::Polygon::Project(double& x1, | |
309 double& y1, | |
310 double& x2, | |
311 double& y2, | |
312 const CoordinateSystem3D& cuttingPlane, | 287 const CoordinateSystem3D& cuttingPlane, |
313 const Vector& estimatedNormal, | 288 const Vector& estimatedNormal, |
314 double estimatedSliceThickness) const | 289 double estimatedSliceThickness) const |
315 { | 290 { |
316 if (points_.size() <= 1) | 291 CoordinateSystem3D geometry; |
317 { | |
318 return false; | |
319 } | |
320 | |
321 Vector normal = estimatedNormal; | |
322 double thickness = estimatedSliceThickness; | 292 double thickness = estimatedSliceThickness; |
293 | |
294 /** | |
295 * 1. Estimate the 3D plane associated with this polygon. | |
296 **/ | |
297 | |
323 if (hasSlice_) | 298 if (hasSlice_) |
324 { | 299 { |
325 normal = geometry_.GetNormal(); | 300 // The exact geometry is known for this slice |
301 geometry = geometry_; | |
326 thickness = sliceThickness_; | 302 thickness = sliceThickness_; |
327 } | 303 } |
304 else if (points_.size() < 2) | |
305 { | |
306 return; | |
307 } | |
308 else | |
309 { | |
310 // Estimate the geometry | |
311 Vector origin = points_[0]; | |
312 | |
313 Vector axisX; | |
314 bool found = false; | |
315 for (size_t i = 1; i < points_.size(); i++) | |
316 { | |
317 axisX = points_[1] - origin; | |
318 if (boost::numeric::ublas::norm_2(axisX) > 10.0 * std::numeric_limits<double>::epsilon()) | |
319 { | |
320 found = true; | |
321 break; | |
322 } | |
323 } | |
324 | |
325 if (!found) | |
326 { | |
327 return; // The polygon is too small to extract a reliable geometry out of it | |
328 } | |
329 | |
330 LinearAlgebra::NormalizeVector(axisX); | |
331 | |
332 Vector axisY; | |
333 LinearAlgebra::CrossProduct(axisY, axisX, estimatedNormal); | |
334 geometry = CoordinateSystem3D(origin, axisX, axisY); | |
335 } | |
336 | |
337 | |
338 /** | |
339 * 2. Project the 3D cutting plane as a 2D line onto the polygon plane. | |
340 **/ | |
341 | |
342 double cuttingX1, cuttingY1, cuttingX2, cuttingY2; | |
343 geometry.ProjectPoint(cuttingX1, cuttingY1, cuttingPlane.GetOrigin()); | |
328 | 344 |
329 bool isOpposite; | 345 bool isOpposite; |
330 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) && | 346 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, geometry.GetNormal(), cuttingPlane.GetAxisX())) |
331 !GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY())) | 347 { |
332 { | 348 geometry.ProjectPoint(cuttingX2, cuttingY2, cuttingPlane.GetOrigin() + cuttingPlane.GetAxisY()); |
333 return false; | 349 } |
334 } | 350 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, geometry.GetNormal(), cuttingPlane.GetAxisY())) |
335 | 351 { |
336 const double d = cuttingPlane.ProjectAlongNormal(cuttingPlane.GetOrigin()); | 352 geometry.ProjectPoint(cuttingX2, cuttingY2, cuttingPlane.GetOrigin() + cuttingPlane.GetAxisX()); |
337 | |
338 Extent2D extent; | |
339 AddSegmentIfIntersection(extent, cuttingPlane, points_[points_.size() - 1], points_[0], d); | |
340 for (size_t i = 1; i < points_.size(); i++) | |
341 { | |
342 AddSegmentIfIntersection(extent, cuttingPlane, points_[i - 1], points_[i], d); | |
343 } | |
344 | |
345 if (extent.GetWidth() > 0 || | |
346 extent.GetHeight() > 0) | |
347 { | |
348 // Let's convert them to 3D world geometry to add the slice thickness | |
349 Vector p1 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX1(), extent.GetY1()) + | |
350 thickness / 2.0 * normal); | |
351 Vector p2 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX2(), extent.GetY2()) - | |
352 thickness / 2.0 * normal); | |
353 | |
354 // Then back to the cutting plane geometry | |
355 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
356 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
357 return true; | |
358 } | 353 } |
359 else | 354 else |
360 { | 355 { |
361 return false; | 356 return; |
357 } | |
358 | |
359 | |
360 /** | |
361 * 3. Compute the intersection of the 2D cutting line with the polygon. | |
362 **/ | |
363 | |
364 // Initialize the projection of a point onto a line: | |
365 // https://stackoverflow.com/a/64330724 | |
366 const double abx = cuttingX2 - cuttingX1; | |
367 const double aby = cuttingY2 - cuttingY1; | |
368 const double denominator = abx * abx + aby * aby; | |
369 if (LinearAlgebra::IsCloseToZero(denominator)) | |
370 { | |
371 return; // Should never happen | |
372 } | |
373 | |
374 std::vector<double> intersections; | |
375 intersections.reserve(points_.size()); | |
376 | |
377 for (size_t i = 0; i < points_.size(); i++) | |
378 { | |
379 double segmentX1, segmentY1, segmentX2, segmentY2; | |
380 geometry.ProjectPoint(segmentX1, segmentY1, points_[i]); | |
381 geometry.ProjectPoint(segmentX2, segmentY2, points_[(i + 1) % points_.size()]); | |
382 | |
383 double x, y; | |
384 if (GeometryToolbox::IntersectLineAndSegment(x, y, cuttingX1, cuttingY1, cuttingX2, cuttingY2, | |
385 segmentX1, segmentY1, segmentX2, segmentY2)) | |
386 { | |
387 // For each polygon segment that intersect the cutting line, | |
388 // register its offset over the cutting line | |
389 const double acx = x - cuttingX1; | |
390 const double acy = y - cuttingY1; | |
391 intersections.push_back((abx * acx + aby * acy) / denominator); | |
392 } | |
393 } | |
394 | |
395 | |
396 /** | |
397 * 4. Sort the intersection offsets, then generates one 2D rectangle on the | |
398 * cutting plane from each pair of successive intersections. | |
399 **/ | |
400 | |
401 std::sort(intersections.begin(), intersections.end()); | |
402 | |
403 if (intersections.size() % 2 == 1) | |
404 { | |
405 return; // Should never happen | |
406 } | |
407 | |
408 for (size_t i = 0; i < intersections.size(); i += 2) | |
409 { | |
410 Vector p1, p2; | |
411 | |
412 { | |
413 // Let's convert them to 3D world geometry to add the slice thickness | |
414 const double x1 = cuttingX1 + intersections[i] * abx; | |
415 const double y1 = cuttingY1 + intersections[i] * aby; | |
416 const double x2 = cuttingX1 + intersections[i + 1] * abx; | |
417 const double y2 = cuttingY1 + intersections[i + 1] * aby; | |
418 | |
419 p1 = (geometry.MapSliceToWorldCoordinates(x1, y1) + thickness / 2.0 * geometry.GetNormal()); | |
420 p2 = (geometry.MapSliceToWorldCoordinates(x2, y2) - thickness / 2.0 * geometry.GetNormal()); | |
421 } | |
422 | |
423 { | |
424 // Then back to the cutting plane geometry | |
425 double x1, y1, x2, y2; | |
426 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
427 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
428 | |
429 target.push_back(Extent2D(x1, y1, x2, y2)); | |
430 } | |
362 } | 431 } |
363 } | 432 } |
364 | 433 |
365 | 434 |
366 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const | 435 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const |
888 std::vector<BoostPolygon> projected; | 957 std::vector<BoostPolygon> projected; |
889 | 958 |
890 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 959 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
891 polygon != structure.polygons_.end(); ++polygon) | 960 polygon != structure.polygons_.end(); ++polygon) |
892 { | 961 { |
893 double x1, y1, x2, y2; | 962 std::list<Extent2D> rectangles; |
894 | 963 polygon->Project(rectangles, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()); |
895 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness())) | 964 |
896 { | 965 for (std::list<Extent2D>::const_iterator it = rectangles.begin(); it != rectangles.end(); ++it) |
897 projected.push_back(CreateRectangle(x1, y1, x2, y2)); | 966 { |
967 projected.push_back(CreateRectangle(it->GetX1(), it->GetY1(), it->GetX2(), it->GetY2())); | |
898 } | 968 } |
899 } | 969 } |
900 | 970 |
901 BoostMultiPolygon merged; | 971 BoostMultiPolygon merged; |
902 Union(merged, projected); | 972 Union(merged, projected); |
918 std::list<Extent2D> rectangles; | 988 std::list<Extent2D> rectangles; |
919 | 989 |
920 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 990 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
921 polygon != structure.polygons_.end(); ++polygon) | 991 polygon != structure.polygons_.end(); ++polygon) |
922 { | 992 { |
923 double x1, y1, x2, y2; | 993 polygon->Project(rectangles, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()); |
924 | |
925 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness())) | |
926 { | |
927 rectangles.push_back(Extent2D(x1, y1, x2, y2)); | |
928 } | |
929 } | 994 } |
930 | 995 |
931 typedef std::list< std::vector<ScenePoint2D> > Contours; | 996 typedef std::list< std::vector<ScenePoint2D> > Contours; |
932 | 997 |
933 Contours contours; | 998 Contours contours; |