Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1908:affde38b84de

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
moved tentative bgo reimplementation of rt-struct into graveyard
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 01 Feb 2022 08:38:32 +0100
parents e318b524ad3f
children 782ba9eb6f22
comparison
equal deleted inserted replaced
1907:0208f99b8bde 1908:affde38b84de
55 # include <boost/geometry.hpp> 55 # include <boost/geometry.hpp>
56 # include <boost/geometry/geometries/point_xy.hpp> 56 # include <boost/geometry/geometries/point_xy.hpp>
57 # include <boost/geometry/geometries/polygon.hpp> 57 # include <boost/geometry/geometries/polygon.hpp>
58 # include <boost/geometry/multi/geometries/multi_polygon.hpp> 58 # include <boost/geometry/multi/geometries/multi_polygon.hpp>
59 #else 59 #else
60 # include "DicomStructureSetUtils.h" // TODO REMOVE
61 # include "UnionOfRectangles.h" 60 # include "UnionOfRectangles.h"
62 #endif 61 #endif
63 62
64 #if defined(_MSC_VER) 63 #if defined(_MSC_VER)
65 # pragma warning(pop) 64 # pragma warning(pop)
115 boost::geometry::append(r, BoostPoint(x2, y2)); 114 boost::geometry::append(r, BoostPoint(x2, y2));
116 boost::geometry::append(r, BoostPoint(x2, y1)); 115 boost::geometry::append(r, BoostPoint(x2, y1));
117 return r; 116 return r;
118 } 117 }
119 118
120 #else
121
122 namespace OrthancStone
123 {
124 static RtStructRectangleInSlab CreateRectangle(float x1, float y1,
125 float x2, float y2)
126 {
127 RtStructRectangleInSlab rect;
128 rect.xmin = std::min(x1, x2);
129 rect.xmax = std::max(x1, x2);
130 rect.ymin = std::min(y1, y2);
131 rect.ymax = std::max(y1, y2);
132 return rect;
133 }
134
135 bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1,
136 const std::pair<RtStructRectangleInSlab, double>& r2)
137 {
138 return r1.second < r2.second;
139 }
140
141 bool CompareSlabsY(const RtStructRectanglesInSlab& r1,
142 const RtStructRectanglesInSlab& r2)
143 {
144 if ((r1.size() == 0) || (r2.size() == 0))
145 return false;
146
147 return r1[0].ymax < r2[0].ymax;
148 }
149 }
150
151 #endif 119 #endif
120
152 121
153 namespace OrthancStone 122 namespace OrthancStone
154 { 123 {
155 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); 124 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042);
156 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); 125 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016);
919 { 888 {
920 chains[i][j] = ScenePoint2D(outer[j].x(), outer[j].y()); 889 chains[i][j] = ScenePoint2D(outer[j].x(), outer[j].y());
921 } 890 }
922 } 891 }
923 892
924 #elif 1 893 #else
925 894
926 std::list<Extent2D> rectangles; 895 std::list<Extent2D> rectangles;
927 896
928 for (Polygons::const_iterator polygon = structure.polygons_.begin(); 897 for (Polygons::const_iterator polygon = structure.polygons_.begin();
929 polygon != structure.polygons_.end(); ++polygon) 898 polygon != structure.polygons_.end(); ++polygon)
946 for (Contours::const_iterator it = contours.begin(); it != contours.end(); ++it) 915 for (Contours::const_iterator it = contours.begin(); it != contours.end(); ++it)
947 { 916 {
948 chains.push_back(*it); 917 chains.push_back(*it);
949 } 918 }
950 919
951 #else
952 // this will contain the intersection of the polygon slab with
953 // the cutting plane, projected on the cutting plane coord system
954 // (that yields a rectangle) + the Z coordinate of the polygon
955 // (this is required to group polygons with the same Z later)
956 std::vector<std::pair<RtStructRectangleInSlab, double> > projected;
957
958 for (Polygons::const_iterator polygon = structure.polygons_.begin();
959 polygon != structure.polygons_.end(); ++polygon)
960 {
961 double x1, y1, x2, y2;
962
963 if (polygon->Project(x1, y1, x2, y2, slice, GetEstimatedNormal(), GetEstimatedSliceThickness()))
964 {
965 double curZ = polygon->GetGeometryOrigin()[2];
966
967 // x1,y1 and x2,y2 are in "slice" coordinates (the cutting plane
968 // geometry)
969 projected.push_back(std::make_pair(CreateRectangle(
970 static_cast<float>(x1),
971 static_cast<float>(y1),
972 static_cast<float>(x2),
973 static_cast<float>(y2)),curZ));
974 }
975 }
976
977 // projected contains a set of rectangles specified by two opposite
978 // corners (x1,y1,x2,y2)
979 // we need to merge them
980 // each slab yields ONE polygon!
981
982 // we need to sorted all the rectangles that originate from the same Z
983 // into lanes. To make sure they are grouped together in the array, we
984 // sort it.
985 std::sort(projected.begin(), projected.end(), CompareRectanglesForProjection);
986
987 std::vector<RtStructRectanglesInSlab> rectanglesForEachSlab;
988 rectanglesForEachSlab.reserve(projected.size());
989
990 double curZ = 0;
991 for (size_t i = 0; i < projected.size(); ++i)
992 {
993 #if 0
994 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab());
995 #else
996 if (i == 0)
997 {
998 curZ = projected[i].second;
999 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab());
1000 }
1001 else
1002 {
1003 // this check is needed to prevent creating a new slab if
1004 // the new polygon is at the same Z coord than last one
1005 if (!LinearAlgebra::IsNear(curZ, projected[i].second))
1006 {
1007 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab());
1008 curZ = projected[i].second;
1009 }
1010 }
1011 #endif 920 #endif
1012 921
1013 rectanglesForEachSlab.back().push_back(projected[i].first);
1014
1015 // as long as they have the same y, we should put them into the same lane
1016 // BUT in Sebastien's code, there is only one polygon per lane.
1017
1018 //std::cout << "rect: xmin = " << rect.xmin << " xmax = " << rect.xmax << " ymin = " << rect.ymin << " ymax = " << rect.ymax << std::endl;
1019 }
1020
1021 // now we need to sort the slabs in increasing Y order (see ConvertListOfSlabsToSegments)
1022 std::sort(rectanglesForEachSlab.begin(), rectanglesForEachSlab.end(), CompareSlabsY);
1023
1024 std::vector< std::pair<ScenePoint2D, ScenePoint2D> > segments;
1025 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, projected.size());
1026
1027 chains.resize(segments.size());
1028 for (size_t i = 0; i < segments.size(); i++)
1029 {
1030 chains[i].resize(2);
1031 chains[i][0] = segments[i].first;
1032 chains[i][1] = segments[i].second;
1033 }
1034 #endif
1035
1036 return true; 922 return true;
1037 } 923 }
1038 else 924 else
1039 { 925 {
1040 return false; 926 return false;