Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet.cpp @ 1000:50e5acf5553b
changed RTSTRUCT rendering from polygons to segments
author | Benjamin Golinvaux <bgo@osimis.io> |
---|---|
date | Fri, 20 Sep 2019 11:59:29 +0200 |
parents | 1f74bc3459ba |
children | 4f28d9459e31 |
comparison
equal
deleted
inserted
replaced
999:2d69b8bee484 | 1000:50e5acf5553b |
---|---|
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. | 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. |
19 **/ | 19 **/ |
20 | 20 |
21 | 21 |
22 #include "DicomStructureSet.h" | 22 #include "DicomStructureSet.h" |
23 #include "DicomStructureSetUtils.h" | |
23 | 24 |
24 #include "../Toolbox/GeometryToolbox.h" | 25 #include "../Toolbox/GeometryToolbox.h" |
25 | 26 |
26 #include <Core/Logging.h> | 27 #include <Core/Logging.h> |
27 #include <Core/OrthancException.h> | 28 #include <Core/OrthancException.h> |
28 #include <Core/Toolbox.h> | 29 #include <Core/Toolbox.h> |
29 #include <Plugins/Samples/Common/FullOrthancDataset.h> | 30 #include <Plugins/Samples/Common/FullOrthancDataset.h> |
30 #include <Plugins/Samples/Common/DicomDatasetReader.h> | 31 #include <Plugins/Samples/Common/DicomDatasetReader.h> |
32 | |
33 #if defined(_MSC_VER) | |
34 #pragma warning(push) | |
35 #pragma warning(disable:4244) | |
36 #endif | |
31 | 37 |
32 #include <limits> | 38 #include <limits> |
33 #include <stdio.h> | 39 #include <stdio.h> |
34 #include <boost/geometry.hpp> | 40 #include <boost/geometry.hpp> |
35 #include <boost/geometry/geometries/point_xy.hpp> | 41 #include <boost/geometry/geometries/point_xy.hpp> |
36 #include <boost/geometry/geometries/polygon.hpp> | 42 #include <boost/geometry/geometries/polygon.hpp> |
37 #include <boost/geometry/multi/geometries/multi_polygon.hpp> | 43 #include <boost/geometry/multi/geometries/multi_polygon.hpp> |
38 | 44 |
45 #if defined(_MSC_VER) | |
46 #pragma warning(pop) | |
47 #endif | |
48 | |
39 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; | 49 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; |
40 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; | 50 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; |
41 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; | 51 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; |
42 | 52 |
43 | 53 |
69 output = tmp; | 79 output = tmp; |
70 } | 80 } |
71 } | 81 } |
72 } | 82 } |
73 | 83 |
84 #ifdef USE_OLD_SJO_CUT_CODE | |
74 | 85 |
75 static BoostPolygon CreateRectangle(float x1, float y1, | 86 static BoostPolygon CreateRectangle(float x1, float y1, |
76 float x2, float y2) | 87 float x2, float y2) |
77 { | 88 { |
78 BoostPolygon r; | 89 BoostPolygon r; |
81 boost::geometry::append(r, BoostPoint(x2, y2)); | 92 boost::geometry::append(r, BoostPoint(x2, y2)); |
82 boost::geometry::append(r, BoostPoint(x2, y1)); | 93 boost::geometry::append(r, BoostPoint(x2, y1)); |
83 return r; | 94 return r; |
84 } | 95 } |
85 | 96 |
86 | 97 #else |
98 | |
99 namespace OrthancStone | |
100 { | |
101 static RtStructRectangleInSlab CreateRectangle(float x1, float y1, | |
102 float x2, float y2) | |
103 { | |
104 RtStructRectangleInSlab rect; | |
105 rect.xmin = std::min(x1, x2); | |
106 rect.xmax = std::max(x1, x2); | |
107 rect.ymin = std::min(y1, y2); | |
108 rect.ymax = std::max(y1, y2); | |
109 return rect; | |
110 } | |
111 | |
112 bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1, const std::pair<RtStructRectangleInSlab, double>& r2) | |
113 { | |
114 return r1.second < r2.second; | |
115 } | |
116 | |
117 bool CompareSlabsY(const RtStructRectanglesInSlab& r1, const RtStructRectanglesInSlab& r2) | |
118 { | |
119 if ((r1.size() == 0) || (r2.size() == 0)) | |
120 return false; | |
121 | |
122 return r1[0].ymax < r2[0].ymax; | |
123 } | |
124 } | |
125 | |
126 #endif | |
87 | 127 |
88 namespace OrthancStone | 128 namespace OrthancStone |
89 { | 129 { |
90 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); | 130 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); |
91 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); | 131 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); |
265 | 305 |
266 bool isOpposite; | 306 bool isOpposite; |
267 if (GeometryToolbox::IsParallelOrOpposite | 307 if (GeometryToolbox::IsParallelOrOpposite |
268 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | 308 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) |
269 { | 309 { |
310 // plane is constant Y | |
311 | |
270 if (y < extent_.GetY1() || | 312 if (y < extent_.GetY1() || |
271 y > extent_.GetY2()) | 313 y > extent_.GetY2()) |
272 { | 314 { |
273 // The polygon does not intersect the input slice | 315 // The polygon does not intersect the input slice |
274 return false; | 316 return false; |
285 { | 327 { |
286 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py | 328 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py |
287 double curX, curY; | 329 double curX, curY; |
288 geometry_.ProjectPoint(curX, curY, points_[i]); | 330 geometry_.ProjectPoint(curX, curY, points_[i]); |
289 | 331 |
332 // if prev* and cur* are on opposite sides of y, this means that the | |
333 // segment intersects the plane. | |
290 if ((prevY < y && curY > y) || | 334 if ((prevY < y && curY > y) || |
291 (prevY > y && curY < y)) | 335 (prevY > y && curY < y)) |
292 { | 336 { |
293 double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); | 337 double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); |
294 xmin = std::min(xmin, p); | 338 xmin = std::min(xmin, p); |
295 xmax = std::max(xmax, p); | 339 xmax = std::max(xmax, p); |
296 isFirst = false; | 340 isFirst = false; |
341 | |
342 // xmin and xmax represent the extent of the rectangle along the | |
343 // intersection between the plane and the polygon geometry | |
344 | |
297 } | 345 } |
298 | 346 |
299 prevX = curX; | 347 prevX = curX; |
300 prevY = curY; | 348 prevY = curY; |
301 } | 349 } |
302 | 350 |
351 // if NO segment intersects the plane | |
303 if (isFirst) | 352 if (isFirst) |
304 { | 353 { |
305 return false; | 354 return false; |
306 } | 355 } |
307 else | 356 else |
308 { | 357 { |
358 // y is the plane y coord in the polygon geometry | |
359 // xmin and xmax are ALSO expressed in the polygon geometry | |
360 | |
361 // let's convert them to 3D world geometry... | |
309 Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + | 362 Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + |
310 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 363 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
311 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - | 364 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - |
312 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 365 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
313 | 366 |
367 // then to the cutting plane geometry... | |
314 slice.ProjectPoint(x1, y1, p1); | 368 slice.ProjectPoint(x1, y1, p1); |
315 slice.ProjectPoint(x2, y2, p2); | 369 slice.ProjectPoint(x2, y2, p2); |
316 return true; | 370 return true; |
317 } | 371 } |
318 } | 372 } |
319 else if (GeometryToolbox::IsParallelOrOpposite | 373 else if (GeometryToolbox::IsParallelOrOpposite |
320 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | 374 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) |
321 { | 375 { |
376 // plane is constant X | |
377 | |
378 /* | |
379 Please read the comments in the section above, by taking into account | |
380 the fact that, in this case, the plane has a constant X, not Y (in | |
381 polygon geometry_ coordinates) | |
382 */ | |
383 | |
322 if (x < extent_.GetX1() || | 384 if (x < extent_.GetX1() || |
323 x > extent_.GetX2()) | 385 x > extent_.GetX2()) |
324 { | 386 { |
325 return false; | 387 return false; |
326 } | 388 } |
739 { | 801 { |
740 return referencedSlices_.begin()->second.geometry_.GetNormal(); | 802 return referencedSlices_.begin()->second.geometry_.GetNormal(); |
741 } | 803 } |
742 } | 804 } |
743 | 805 |
744 | 806 #ifdef USE_OLD_SJO_CUT_CODE |
745 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons, | 807 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<Point2D> >& polygons, |
746 const Structure& structure, | 808 const Structure& structure, |
747 const CoordinateSystem3D& slice) const | 809 const CoordinateSystem3D& slice) const |
748 { | 810 #else |
811 bool DicomStructureSet::ProjectStructure(std::vector< std::pair<Point2D, Point2D> >& segments, | |
812 const Structure& structure, | |
813 const CoordinateSystem3D& slice) const | |
814 #endif | |
815 { | |
816 #ifdef USE_OLD_SJO_CUT_CODE | |
749 polygons.clear(); | 817 polygons.clear(); |
818 #else | |
819 segments.clear(); | |
820 #endif | |
750 | 821 |
751 Vector normal = GetNormal(); | 822 Vector normal = GetNormal(); |
752 | 823 |
753 bool isOpposite; | 824 bool isOpposite; |
754 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) | 825 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) |
758 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 829 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
759 polygon != structure.polygons_.end(); ++polygon) | 830 polygon != structure.polygons_.end(); ++polygon) |
760 { | 831 { |
761 if (polygon->IsOnSlice(slice)) | 832 if (polygon->IsOnSlice(slice)) |
762 { | 833 { |
763 polygons.push_back(std::vector<PolygonPoint2D>()); | 834 #ifdef USE_OLD_SJO_CUT_CODE |
835 polygons.push_back(std::vector<Point2D>()); | |
764 | 836 |
765 for (Points::const_iterator p = polygon->GetPoints().begin(); | 837 for (Points::const_iterator p = polygon->GetPoints().begin(); |
766 p != polygon->GetPoints().end(); ++p) | 838 p != polygon->GetPoints().end(); ++p) |
767 { | 839 { |
768 double x, y; | 840 double x, y; |
769 slice.ProjectPoint(x, y, *p); | 841 slice.ProjectPoint(x, y, *p); |
770 polygons.back().push_back(std::make_pair(x, y)); | 842 polygons.back().push_back(Point2D(x, y)); |
771 } | 843 } |
844 #else | |
845 // we need to add all the segments corresponding to this polygon | |
846 const std::vector<Vector>& points3D = polygon->GetPoints(); | |
847 if (points3D.size() >= 3) | |
848 { | |
849 Point2D prev2D; | |
850 { | |
851 Vector prev = points3D[0]; | |
852 double prevX, prevY; | |
853 slice.ProjectPoint(prevX, prevY, prev); | |
854 prev2D = Point2D(prevX, prevY); | |
855 } | |
856 | |
857 size_t pointCount = points3D.size(); | |
858 for (size_t ipt = 1; ipt < pointCount; ++ipt) | |
859 { | |
860 Vector next = points3D[ipt]; | |
861 double nextX, nextY; | |
862 slice.ProjectPoint(nextX, nextY, next); | |
863 Point2D next2D(nextX, nextY); | |
864 segments.push_back(std::pair<Point2D, Point2D>(prev2D, next2D)); | |
865 prev2D = next2D; | |
866 } | |
867 } | |
868 else | |
869 { | |
870 LOG(ERROR) << "Contour with less than 3 points!"; | |
871 // !!! | |
872 } | |
873 #endif | |
772 } | 874 } |
773 } | 875 } |
774 | 876 |
775 return true; | 877 return true; |
776 } | 878 } |
777 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | 879 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || |
778 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 880 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) |
779 { | 881 { |
780 #if 1 | 882 #if 1 |
781 // Sagittal or coronal projection | 883 // Sagittal or coronal projection |
884 | |
885 #ifdef USE_OLD_SJO_CUT_CODE | |
782 std::vector<BoostPolygon> projected; | 886 std::vector<BoostPolygon> projected; |
783 | 887 #else |
888 // this will contain the intersection of the polygon slab with | |
889 // the cutting plane, projected on the cutting plane coord system | |
890 // (that yields a rectangle) + the Z coordinate of the polygon | |
891 // (this is required to group polygons with the same Z later) | |
892 std::vector<std::pair<RtStructRectangleInSlab, double> > projected; | |
893 #endif | |
784 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 894 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
785 polygon != structure.polygons_.end(); ++polygon) | 895 polygon != structure.polygons_.end(); ++polygon) |
786 { | 896 { |
787 double x1, y1, x2, y2; | 897 double x1, y1, x2, y2; |
898 | |
788 if (polygon->Project(x1, y1, x2, y2, slice)) | 899 if (polygon->Project(x1, y1, x2, y2, slice)) |
789 { | 900 { |
790 projected.push_back(CreateRectangle( | 901 double curZ = polygon->GetGeometryOrigin()[2]; |
902 | |
903 // x1,y1 and x2,y2 are in "slice" coordinates (the cutting plane | |
904 // geometry) | |
905 projected.push_back(std::make_pair(CreateRectangle( | |
791 static_cast<float>(x1), | 906 static_cast<float>(x1), |
792 static_cast<float>(y1), | 907 static_cast<float>(y1), |
793 static_cast<float>(x2), | 908 static_cast<float>(x2), |
794 static_cast<float>(y2))); | 909 static_cast<float>(y2)),curZ)); |
795 } | 910 } |
796 } | 911 } |
797 | 912 #ifndef USE_OLD_SJO_CUT_CODE |
913 // projected contains a set of rectangles specified by two opposite | |
914 // corners (x1,y1,x2,y2) | |
915 // we need to merge them | |
916 // each slab yields ONE polygon! | |
917 | |
918 // we need to sorted all the rectangles that originate from the same Z | |
919 // into lanes. To make sure they are grouped together in the array, we | |
920 // sort it. | |
921 std::sort(projected.begin(), projected.end(), CompareRectanglesForProjection); | |
922 | |
923 std::vector<RtStructRectanglesInSlab> rectanglesForEachSlab; | |
924 rectanglesForEachSlab.reserve(projected.size()); | |
925 | |
926 double curZ = 0; | |
927 for (size_t i = 0; i < projected.size(); ++i) | |
928 { | |
929 #if 0 | |
930 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
931 #else | |
932 if (i == 0) | |
933 { | |
934 curZ = projected[i].second; | |
935 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
936 } | |
937 else | |
938 { | |
939 // this check is needed to prevent creating a new slab if | |
940 // the new polygon is at the same Z coord than last one | |
941 if (!LinearAlgebra::IsNear(curZ, projected[i].second)) | |
942 { | |
943 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
944 curZ = projected[i].second; | |
945 } | |
946 } | |
947 #endif | |
948 | |
949 rectanglesForEachSlab.back().push_back(projected[i].first); | |
950 | |
951 // as long as they have the same y, we should put them into the same lane | |
952 // BUT in Sebastien's code, there is only one polygon per lane. | |
953 | |
954 //std::cout << "rect: xmin = " << rect.xmin << " xmax = " << rect.xmax << " ymin = " << rect.ymin << " ymax = " << rect.ymax << std::endl; | |
955 } | |
956 | |
957 // now we need to sort the slabs in increasing Y order (see ConvertListOfSlabsToSegments) | |
958 std::sort(rectanglesForEachSlab.begin(), rectanglesForEachSlab.end(), CompareSlabsY); | |
959 | |
960 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, projected.size()); | |
961 #else | |
798 BoostMultiPolygon merged; | 962 BoostMultiPolygon merged; |
799 Union(merged, projected); | 963 Union(merged, projected); |
800 | 964 |
801 polygons.resize(merged.size()); | 965 polygons.resize(merged.size()); |
802 for (size_t i = 0; i < merged.size(); i++) | 966 for (size_t i = 0; i < merged.size(); i++) |
804 const std::vector<BoostPoint>& outer = merged[i].outer(); | 968 const std::vector<BoostPoint>& outer = merged[i].outer(); |
805 | 969 |
806 polygons[i].resize(outer.size()); | 970 polygons[i].resize(outer.size()); |
807 for (size_t j = 0; j < outer.size(); j++) | 971 for (size_t j = 0; j < outer.size(); j++) |
808 { | 972 { |
809 polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); | 973 polygons[i][j] = Point2D(outer[j].x(), outer[j].y()); |
810 } | 974 } |
811 } | 975 } |
976 #endif | |
977 | |
812 #else | 978 #else |
813 for (Polygons::iterator polygon = structure.polygons_.begin(); | 979 for (Polygons::iterator polygon = structure.polygons_.begin(); |
814 polygon != structure.polygons_.end(); ++polygon) | 980 polygon != structure.polygons_.end(); ++polygon) |
815 { | 981 { |
816 double x1, y1, x2, y2; | 982 double x1, y1, x2, y2; |
817 if (polygon->Project(x1, y1, x2, y2, slice)) | 983 if (polygon->Project(x1, y1, x2, y2, slice)) |
818 { | 984 { |
819 std::vector<PolygonPoint2D> p(4); | 985 std::vector<Point2D> p(4); |
820 p[0] = std::make_pair(x1, y1); | 986 p[0] = std::make_pair(x1, y1); |
821 p[1] = std::make_pair(x2, y1); | 987 p[1] = std::make_pair(x2, y1); |
822 p[2] = std::make_pair(x2, y2); | 988 p[2] = std::make_pair(x2, y2); |
823 p[3] = std::make_pair(x1, y2); | 989 p[3] = std::make_pair(x1, y2); |
824 polygons.push_back(p); | 990 polygons.push_back(p); |