Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1896:b3c08e607d9f
simplified signature of DicomStructureSet::ProjectStructure()
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 26 Jan 2022 17:19:37 +0100 |
parents | 14c8f339d480 |
children | 144f8f82c15a |
comparison
equal
deleted
inserted
replaced
1895:14c8f339d480 | 1896:b3c08e607d9f |
---|---|
19 * License along with this program. If not, see | 19 * License along with this program. If not, see |
20 * <http://www.gnu.org/licenses/>. | 20 * <http://www.gnu.org/licenses/>. |
21 **/ | 21 **/ |
22 | 22 |
23 | 23 |
24 #define USE_BOOST_UNION_FOR_POLYGONS 1 | |
25 | |
24 #include "DicomStructureSet.h" | 26 #include "DicomStructureSet.h" |
25 #include "DicomStructureSetUtils.h" // TODO REMOVE | 27 #include "DicomStructureSetUtils.h" // TODO REMOVE |
26 | 28 |
27 #include "BucketAccumulator2D.h" | 29 #include "BucketAccumulator2D.h" |
28 #include "GenericToolbox.h" | 30 #include "GenericToolbox.h" |
40 | 42 |
41 #if STONE_TIME_BLOCKING_OPS | 43 #if STONE_TIME_BLOCKING_OPS |
42 # include <boost/date_time/posix_time/posix_time.hpp> | 44 # include <boost/date_time/posix_time/posix_time.hpp> |
43 #endif | 45 #endif |
44 | 46 |
47 #if !defined(USE_BOOST_UNION_FOR_POLYGONS) | |
48 # error Macro USE_BOOST_UNION_FOR_POLYGONS must be defined | |
49 #endif | |
50 | |
45 #include <limits> | 51 #include <limits> |
46 #include <stdio.h> | 52 #include <stdio.h> |
47 #include <boost/math/constants/constants.hpp> | 53 #include <boost/math/constants/constants.hpp> |
48 #include <boost/geometry.hpp> | 54 |
49 #include <boost/geometry/geometries/point_xy.hpp> | 55 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
50 #include <boost/geometry/geometries/polygon.hpp> | 56 # include <boost/geometry.hpp> |
51 #include <boost/geometry/multi/geometries/multi_polygon.hpp> | 57 # include <boost/geometry/geometries/point_xy.hpp> |
58 # include <boost/geometry/geometries/polygon.hpp> | |
59 # include <boost/geometry/multi/geometries/multi_polygon.hpp> | |
60 #endif | |
52 | 61 |
53 #if defined(_MSC_VER) | 62 #if defined(_MSC_VER) |
54 # pragma warning(pop) | 63 # pragma warning(pop) |
55 #endif | 64 #endif |
56 | 65 |
57 #if ORTHANC_ENABLE_DCMTK == 1 | 66 #if ORTHANC_ENABLE_DCMTK == 1 |
58 # include "ParsedDicomDataset.h" | 67 # include "ParsedDicomDataset.h" |
59 #endif | 68 #endif |
60 | 69 |
70 | |
71 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
61 | 72 |
62 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; | 73 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; |
63 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; | 74 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; |
64 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; | 75 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; |
65 | |
66 | 76 |
67 static void Union(BoostMultiPolygon& output, | 77 static void Union(BoostMultiPolygon& output, |
68 std::vector<BoostPolygon>& input) | 78 std::vector<BoostPolygon>& input) |
69 { | 79 { |
70 for (size_t i = 0; i < input.size(); i++) | 80 for (size_t i = 0; i < input.size(); i++) |
91 boost::geometry::union_(output, input[i], tmp); | 101 boost::geometry::union_(output, input[i], tmp); |
92 output = tmp; | 102 output = tmp; |
93 } | 103 } |
94 } | 104 } |
95 } | 105 } |
96 | |
97 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
98 | 106 |
99 static BoostPolygon CreateRectangle(float x1, float y1, | 107 static BoostPolygon CreateRectangle(float x1, float y1, |
100 float x2, float y2) | 108 float x2, float y2) |
101 { | 109 { |
102 BoostPolygon r; | 110 BoostPolygon r; |
120 rect.ymin = std::min(y1, y2); | 128 rect.ymin = std::min(y1, y2); |
121 rect.ymax = std::max(y1, y2); | 129 rect.ymax = std::max(y1, y2); |
122 return rect; | 130 return rect; |
123 } | 131 } |
124 | 132 |
125 bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1, const std::pair<RtStructRectangleInSlab, double>& r2) | 133 bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1, |
134 const std::pair<RtStructRectangleInSlab, double>& r2) | |
126 { | 135 { |
127 return r1.second < r2.second; | 136 return r1.second < r2.second; |
128 } | 137 } |
129 | 138 |
130 bool CompareSlabsY(const RtStructRectanglesInSlab& r1, const RtStructRectanglesInSlab& r2) | 139 bool CompareSlabsY(const RtStructRectanglesInSlab& r1, |
140 const RtStructRectanglesInSlab& r2) | |
131 { | 141 { |
132 if ((r1.size() == 0) || (r2.size() == 0)) | 142 if ((r1.size() == 0) || (r2.size() == 0)) |
133 return false; | 143 return false; |
134 | 144 |
135 return r1[0].ymax < r2[0].ymax; | 145 return r1[0].ymax < r2[0].ymax; |
827 { | 837 { |
828 return referencedSlices_.begin()->second.geometry_.GetNormal(); | 838 return referencedSlices_.begin()->second.geometry_.GetNormal(); |
829 } | 839 } |
830 } | 840 } |
831 | 841 |
832 bool DicomStructureSet::ProjectStructure( | 842 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<ScenePoint2D> >& chains, |
833 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | 843 const Structure& structure, |
834 std::vector< std::vector<ScenePoint2D> >& polygons, | 844 const CoordinateSystem3D& sourceSlice) const |
835 #else | |
836 std::vector< std::pair<ScenePoint2D, ScenePoint2D> >& segments, | |
837 #endif | |
838 const Structure& structure, | |
839 const CoordinateSystem3D& sourceSlice) const | |
840 { | 845 { |
841 const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); | 846 const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); |
842 | 847 |
843 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | 848 chains.clear(); |
844 polygons.clear(); | |
845 #else | |
846 segments.clear(); | |
847 #endif | |
848 | 849 |
849 Vector normal = GetNormal(); | 850 Vector normal = GetNormal(); |
850 | 851 |
851 bool isOpposite; | 852 bool isOpposite; |
852 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) | 853 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) |
853 { | 854 { |
854 // This is an axial projection | 855 // This is an axial projection |
855 | 856 |
857 chains.reserve(structure.polygons_.size()); | |
858 | |
856 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 859 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
857 polygon != structure.polygons_.end(); ++polygon) | 860 polygon != structure.polygons_.end(); ++polygon) |
858 { | 861 { |
859 if (polygon->IsOnSlice(slice)) | 862 const Points& points = polygon->GetPoints(); |
860 { | 863 |
861 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | 864 if (polygon->IsOnSlice(slice) && |
862 polygons.push_back(std::vector<ScenePoint2D>()); | 865 !points.empty()) |
866 { | |
867 chains.push_back(std::vector<ScenePoint2D>()); | |
868 chains.back().reserve(points.size() + 1); | |
863 | 869 |
864 for (Points::const_iterator p = polygon->GetPoints().begin(); | 870 for (Points::const_iterator p = points.begin(); |
865 p != polygon->GetPoints().end(); ++p) | 871 p != points.end(); ++p) |
866 { | 872 { |
867 double x, y; | 873 double x, y; |
868 slice.ProjectPoint2(x, y, *p); | 874 slice.ProjectPoint2(x, y, *p); |
869 polygons.back().push_back(ScenePoint2D(x, y)); | 875 chains.back().push_back(ScenePoint2D(x, y)); |
870 } | 876 } |
871 #else | 877 |
872 // we need to add all the segments corresponding to this polygon | 878 double x0, y0; |
873 const std::vector<Vector>& points3D = polygon->GetPoints(); | 879 slice.ProjectPoint2(x0, y0, points.front()); |
874 if (points3D.size() >= 3) | 880 chains.back().push_back(ScenePoint2D(x0, y0)); |
875 { | |
876 ScenePoint2D prev2D; | |
877 { | |
878 Vector prev = points3D[0]; | |
879 double prevX, prevY; | |
880 slice.ProjectPoint2(prevX, prevY, prev); | |
881 prev2D = ScenePoint2D(prevX, prevY); | |
882 } | |
883 | |
884 size_t pointCount = points3D.size(); | |
885 for (size_t ipt = 1; ipt < pointCount; ++ipt) | |
886 { | |
887 Vector next = points3D[ipt]; | |
888 double nextX, nextY; | |
889 slice.ProjectPoint2(nextX, nextY, next); | |
890 ScenePoint2D next2D(nextX, nextY); | |
891 segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(prev2D, next2D)); | |
892 prev2D = next2D; | |
893 } | |
894 } | |
895 else | |
896 { | |
897 LOG(ERROR) << "Contour with less than 3 points!"; | |
898 // !!! | |
899 } | |
900 #endif | |
901 } | 881 } |
902 } | 882 } |
903 | 883 |
904 return true; | 884 return true; |
905 } | 885 } |
906 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | 886 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || |
907 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 887 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) |
908 { | 888 { |
909 #if 1 | |
910 // Sagittal or coronal projection | 889 // Sagittal or coronal projection |
911 | 890 |
912 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | 891 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
913 std::vector<BoostPolygon> projected; | 892 std::vector<BoostPolygon> projected; |
914 | 893 |
995 } | 974 } |
996 | 975 |
997 // now we need to sort the slabs in increasing Y order (see ConvertListOfSlabsToSegments) | 976 // now we need to sort the slabs in increasing Y order (see ConvertListOfSlabsToSegments) |
998 std::sort(rectanglesForEachSlab.begin(), rectanglesForEachSlab.end(), CompareSlabsY); | 977 std::sort(rectanglesForEachSlab.begin(), rectanglesForEachSlab.end(), CompareSlabsY); |
999 | 978 |
979 std::vector< std::pair<ScenePoint2D, ScenePoint2D> > segments; | |
1000 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, projected.size()); | 980 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, projected.size()); |
981 | |
982 chains.resize(segments.size()); | |
983 for (size_t i = 0; i < segments.size(); i++) | |
984 { | |
985 chains[i].resize(2); | |
986 chains[i][0] = segments[i].first; | |
987 chains[i][1] = segments[i].second; | |
988 } | |
989 | |
1001 #else | 990 #else |
1002 BoostMultiPolygon merged; | 991 BoostMultiPolygon merged; |
1003 Union(merged, projected); | 992 Union(merged, projected); |
1004 | 993 |
1005 polygons.resize(merged.size()); | 994 chains.resize(merged.size()); |
1006 for (size_t i = 0; i < merged.size(); i++) | 995 for (size_t i = 0; i < merged.size(); i++) |
1007 { | 996 { |
1008 const std::vector<BoostPoint>& outer = merged[i].outer(); | 997 const std::vector<BoostPoint>& outer = merged[i].outer(); |
1009 | 998 |
1010 polygons[i].resize(outer.size()); | 999 chains[i].resize(outer.size()); |
1011 for (size_t j = 0; j < outer.size(); j++) | 1000 for (size_t j = 0; j < outer.size(); j++) |
1012 { | 1001 { |
1013 polygons[i][j] = ScenePoint2D(outer[j].x(), outer[j].y()); | 1002 chains[i][j] = ScenePoint2D(outer[j].x(), outer[j].y()); |
1014 } | 1003 } |
1015 } | 1004 } |
1016 #endif | |
1017 | |
1018 #else | |
1019 for (Polygons::iterator polygon = structure.polygons_.begin(); | |
1020 polygon != structure.polygons_.end(); ++polygon) | |
1021 { | |
1022 double x1, y1, x2, y2; | |
1023 if (polygon->Project(x1, y1, x2, y2, slice)) | |
1024 { | |
1025 std::vector<ScenePoint2D> p(4); | |
1026 p[0] = std::make_pair(x1, y1); | |
1027 p[1] = std::make_pair(x2, y1); | |
1028 p[2] = std::make_pair(x2, y2); | |
1029 p[3] = std::make_pair(x1, y2); | |
1030 polygons.push_back(p); | |
1031 } | |
1032 } | |
1033 #endif | 1005 #endif |
1034 | 1006 |
1035 return true; | 1007 return true; |
1036 } | 1008 } |
1037 else | 1009 else |
1044 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, | 1016 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, |
1045 const CoordinateSystem3D& plane, | 1017 const CoordinateSystem3D& plane, |
1046 size_t structureIndex, | 1018 size_t structureIndex, |
1047 const Color& color) const | 1019 const Color& color) const |
1048 { | 1020 { |
1049 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | 1021 std::vector< std::vector<ScenePoint2D> > chains; |
1050 std::vector< std::vector<ScenePoint2D> > polygons; | |
1051 if (ProjectStructure(polygons, structureIndex, plane)) | |
1052 { | |
1053 for (size_t j = 0; j < polygons.size(); j++) | |
1054 { | |
1055 std::vector<ScenePoint2D> chain; | |
1056 chain.reserve(polygons[j].size()); | |
1057 | |
1058 for (size_t k = 0; k < polygons[j].size(); k++) | |
1059 { | |
1060 chain.push_back(ScenePoint2D(polygons[j][k].x, polygons[j][k].y)); | |
1061 } | |
1062 | |
1063 layer.AddChain(chain, true, color.GetRed(), color.GetGreen(), color.GetBlue()); | |
1064 } | |
1065 } | |
1066 | 1022 |
1067 #else | 1023 if (ProjectStructure(chains, structureIndex, plane)) |
1068 std::vector< std::pair<ScenePoint2D, ScenePoint2D> > segments; | 1024 { |
1069 | 1025 for (size_t j = 0; j < chains.size(); j++) |
1070 if (ProjectStructure(segments, structureIndex, plane)) | 1026 { |
1071 { | 1027 layer.AddChain(chains[j], false, color.GetRed(), color.GetGreen(), color.GetBlue()); |
1072 for (size_t j = 0; j < segments.size(); j++) | 1028 } |
1073 { | 1029 } |
1074 std::vector<ScenePoint2D> chain(2); | |
1075 chain[0] = ScenePoint2D(segments[j].first.GetX(), segments[j].first.GetY()); | |
1076 chain[1] = ScenePoint2D(segments[j].second.GetX(), segments[j].second.GetY()); | |
1077 layer.AddChain(chain, false, color.GetRed(), color.GetGreen(), color.GetBlue()); | |
1078 } | |
1079 } | |
1080 #endif | |
1081 } | 1030 } |
1082 | 1031 |
1083 | 1032 |
1084 void DicomStructureSet::GetStructurePoints(std::list< std::vector<Vector> >& target, | 1033 void DicomStructureSet::GetStructurePoints(std::list< std::vector<Vector> >& target, |
1085 size_t structureIndex, | 1034 size_t structureIndex, |