Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet.cpp @ 126:c9e88e7935a4 wasm
rt-struct projection
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 08 Nov 2017 17:27:38 +0100 |
parents | 44fc253d4876 |
children | 3e6163a53b16 |
comparison
equal
deleted
inserted
replaced
125:44fc253d4876 | 126:c9e88e7935a4 |
---|---|
26 #include <Core/Logging.h> | 26 #include <Core/Logging.h> |
27 #include <Core/OrthancException.h> | 27 #include <Core/OrthancException.h> |
28 #include <Plugins/Samples/Common/FullOrthancDataset.h> | 28 #include <Plugins/Samples/Common/FullOrthancDataset.h> |
29 #include <Plugins/Samples/Common/DicomDatasetReader.h> | 29 #include <Plugins/Samples/Common/DicomDatasetReader.h> |
30 | 30 |
31 #include <limits> | |
31 #include <stdio.h> | 32 #include <stdio.h> |
32 #include <boost/lexical_cast.hpp> | 33 #include <boost/lexical_cast.hpp> |
34 #include <boost/geometry.hpp> | |
35 #include <boost/geometry/geometries/point_xy.hpp> | |
36 #include <boost/geometry/geometries/polygon.hpp> | |
37 | |
38 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; | |
39 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; | |
40 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; | |
41 | |
42 | |
43 static void Union(BoostMultiPolygon& output, | |
44 std::vector<BoostPolygon>& input) | |
45 { | |
46 for (size_t i = 0; i < input.size(); i++) | |
47 { | |
48 boost::geometry::correct(input[i]); | |
49 } | |
50 | |
51 if (input.size() == 0) | |
52 { | |
53 output.clear(); | |
54 } | |
55 else if (input.size() == 1) | |
56 { | |
57 output.resize(1); | |
58 output[0] = input[0]; | |
59 } | |
60 else | |
61 { | |
62 boost::geometry::union_(input[0], input[1], output); | |
63 | |
64 for (size_t i = 0; i < input.size(); i++) | |
65 { | |
66 BoostMultiPolygon tmp; | |
67 boost::geometry::union_(output, input[i], tmp); | |
68 output = tmp; | |
69 } | |
70 } | |
71 } | |
72 | |
73 | |
74 static BoostPolygon CreateRectangle(float x1, float y1, | |
75 float x2, float y2) | |
76 { | |
77 BoostPolygon r; | |
78 boost::geometry::append(r, BoostPoint(x1, y1)); | |
79 boost::geometry::append(r, BoostPoint(x1, y2)); | |
80 boost::geometry::append(r, BoostPoint(x2, y2)); | |
81 boost::geometry::append(r, BoostPoint(x2, y1)); | |
82 return r; | |
83 } | |
84 | |
85 | |
33 | 86 |
34 namespace OrthancStone | 87 namespace OrthancStone |
35 { | 88 { |
36 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); | 89 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); |
37 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); | 90 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); |
447 } | 500 } |
448 | 501 |
449 | 502 |
450 Vector DicomStructureSet::GetNormal() const | 503 Vector DicomStructureSet::GetNormal() const |
451 { | 504 { |
452 if (!referencedSlices_.empty()) | 505 if (referencedSlices_.empty()) |
453 { | 506 { |
454 Vector v; | 507 Vector v; |
455 GeometryToolbox::AssignVector(v, 0, 0, 1); | 508 GeometryToolbox::AssignVector(v, 0, 0, 1); |
456 return v; | 509 return v; |
457 } | 510 } |
538 return true; | 591 return true; |
539 } | 592 } |
540 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | 593 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || |
541 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 594 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) |
542 { | 595 { |
543 // Sagittal or coronal projections | 596 // Sagittal or coronal projection |
597 | |
598 std::vector<BoostPolygon> projected; | |
599 | |
600 for (Polygons::iterator polygon = structure.polygons_.begin(); | |
601 polygon != structure.polygons_.end(); ++polygon) | |
602 { | |
603 polygon->UpdateReferencedSlice(referencedSlices_); | |
604 | |
605 // First, check that the polygon intersects the input slice | |
606 double zmin = std::numeric_limits<double>::infinity(); | |
607 double zmax = -std::numeric_limits<double>::infinity(); | |
608 | |
609 double zref = slice.ProjectAlongNormal(slice.GetOrigin()); | |
610 | |
611 for (Points::const_iterator p = polygon->GetPoints().begin(); | |
612 p != polygon->GetPoints().end(); ++p) | |
613 { | |
614 double z = slice.ProjectAlongNormal(*p); | |
615 zmin = std::min(zmin, z); | |
616 zmax = std::max(zmax, z); | |
617 } | |
618 | |
619 if (zmin < zref && zref < zmax) | |
620 { | |
621 // The polygon intersect the input slice, project the polygon | |
622 double xmin = std::numeric_limits<double>::infinity(); | |
623 double xmax = -std::numeric_limits<double>::infinity(); | |
624 double ymin = std::numeric_limits<double>::infinity(); | |
625 double ymax = -std::numeric_limits<double>::infinity(); | |
626 | |
627 for (Points::const_iterator p = polygon->GetPoints().begin(); | |
628 p != polygon->GetPoints().end(); ++p) | |
629 { | |
630 double x, y; | |
631 slice.ProjectPoint(x, y, *p); | |
632 xmin = std::min(xmin, x); | |
633 xmax = std::max(xmax, x); | |
634 ymin = std::min(ymin, y); | |
635 ymax = std::max(ymax, y); | |
636 } | |
637 | |
638 if (GeometryToolbox::IsNear(xmin, xmax)) | |
639 { | |
640 projected.push_back(CreateRectangle(xmin - polygon->GetSliceThickness() / 2.0, ymin, | |
641 xmax + polygon->GetSliceThickness() / 2.0, ymax)); | |
642 } | |
643 else if (GeometryToolbox::IsNear(ymin, ymax)) | |
644 { | |
645 projected.push_back(CreateRectangle(xmin, ymin - polygon->GetSliceThickness() / 2.0, | |
646 xmax, ymax + polygon->GetSliceThickness() / 2.0)); | |
647 } | |
648 else | |
649 { | |
650 // Should not happen | |
651 } | |
652 } | |
653 } | |
654 | |
655 BoostMultiPolygon merged; | |
656 Union(merged, projected); | |
657 | |
658 polygons.resize(merged.size()); | |
659 for (size_t i = 0; i < merged.size(); i++) | |
660 { | |
661 const std::vector<BoostPoint>& outer = merged[i].outer(); | |
662 | |
663 polygons[i].resize(outer.size()); | |
664 for (size_t j = 0; j < outer.size(); j++) | |
665 { | |
666 polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); | |
667 } | |
668 } | |
544 | 669 |
670 return true; | |
671 } | |
672 else | |
673 { | |
545 return false; | 674 return false; |
546 } | 675 } |
547 else | |
548 { | |
549 return false; | |
550 } | |
551 } | 676 } |
552 } | 677 } |