Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1896:b3c08e607d9f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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,