Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1898:a5e54bd87b25

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
try to use UnionOfRectangles in DicomStructureSet
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 26 Jan 2022 19:42:04 +0100
parents 144f8f82c15a
children 738814c24574
comparison
equal deleted inserted replaced
1897:144f8f82c15a 1898:a5e54bd87b25
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 24 #define USE_BOOST_UNION_FOR_POLYGONS 0
25 25
26 #include "DicomStructureSet.h" 26 #include "DicomStructureSet.h"
27 #include "DicomStructureSetUtils.h" // TODO REMOVE
28 27
29 #include "BucketAccumulator2D.h" 28 #include "BucketAccumulator2D.h"
30 #include "GenericToolbox.h" 29 #include "GenericToolbox.h"
31 #include "GeometryToolbox.h" 30 #include "GeometryToolbox.h"
32 #include "OrthancDatasets/DicomDatasetReader.h" 31 #include "OrthancDatasets/DicomDatasetReader.h"
55 #if USE_BOOST_UNION_FOR_POLYGONS == 1 54 #if USE_BOOST_UNION_FOR_POLYGONS == 1
56 # include <boost/geometry.hpp> 55 # include <boost/geometry.hpp>
57 # include <boost/geometry/geometries/point_xy.hpp> 56 # include <boost/geometry/geometries/point_xy.hpp>
58 # include <boost/geometry/geometries/polygon.hpp> 57 # include <boost/geometry/geometries/polygon.hpp>
59 # include <boost/geometry/multi/geometries/multi_polygon.hpp> 58 # include <boost/geometry/multi/geometries/multi_polygon.hpp>
59 #else
60 # include "DicomStructureSetUtils.h" // TODO REMOVE
61 # include "UnionOfRectangles.h"
60 #endif 62 #endif
61 63
62 #if defined(_MSC_VER) 64 #if defined(_MSC_VER)
63 # pragma warning(pop) 65 # pragma warning(pop)
64 #endif 66 #endif
188 std::string value; 190 std::string value;
189 return (dataset.GetStringValue(value, tag) && 191 return (dataset.GetStringValue(value, tag) &&
190 GenericToolbox::FastParseVector(target, value)); 192 GenericToolbox::FastParseVector(target, value));
191 } 193 }
192 194
193 void DicomStructureSet::Polygon::CheckPointIsOnSlice(const Vector& v) const
194 {
195 if (hasSlice_)
196 {
197 double magnitude =
198 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal());
199 if(!LinearAlgebra::IsNear(
200 magnitude,
201 projectionAlongNormal_,
202 sliceThickness_ / 2.0 /* in mm */ ))
203 {
204 LOG(ERROR) << "This RT-STRUCT contains a point that is off the "
205 << "slice of its instance | "
206 << "magnitude = " << magnitude << " | "
207 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | "
208 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0);
209
210 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
211 }
212 }
213 }
214
215 bool DicomStructureSet::Polygon::IsPointOnSliceIfAny(const Vector& v) const 195 bool DicomStructureSet::Polygon::IsPointOnSliceIfAny(const Vector& v) const
216 { 196 {
217 if (hasSlice_) 197 if (hasSlice_)
218 { 198 {
219 double magnitude = 199 double magnitude =
267 { 247 {
268 return false; 248 return false;
269 } 249 }
270 else 250 else
271 { 251 {
252 // return true; // TODO - TEST
253
272 const CoordinateSystem3D& geometry = it->second.geometry_; 254 const CoordinateSystem3D& geometry = it->second.geometry_;
273 255
274 hasSlice_ = true; 256 hasSlice_ = true;
275 geometry_ = geometry; 257 geometry_ = geometry;
276 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); 258 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal());
313 double& y1, 295 double& y1,
314 double& x2, 296 double& x2,
315 double& y2, 297 double& y2,
316 const CoordinateSystem3D& slice) const 298 const CoordinateSystem3D& slice) const
317 { 299 {
318 // TODO: optimize this method using a sweep-line algorithm for polygons
319
320 if (!hasSlice_ || 300 if (!hasSlice_ ||
321 points_.size() <= 1) 301 points_.size() <= 1)
322 { 302 {
323 return false; 303 return false;
324 } 304 }
914 for (size_t j = 0; j < outer.size(); j++) 894 for (size_t j = 0; j < outer.size(); j++)
915 { 895 {
916 chains[i][j] = ScenePoint2D(outer[j].x(), outer[j].y()); 896 chains[i][j] = ScenePoint2D(outer[j].x(), outer[j].y());
917 } 897 }
918 } 898 }
899
900 #elif 0
901
902 // TODO - Fix possible infinite loop in UnionOfRectangles
903
904 std::list<Extent2D> rectangles;
905
906 for (Polygons::const_iterator polygon = structure.polygons_.begin();
907 polygon != structure.polygons_.end(); ++polygon)
908 {
909 double x1, y1, x2, y2;
910
911 if (polygon->Project(x1, y1, x2, y2, slice))
912 {
913 rectangles.push_back(Extent2D(x1, y1, x2, y2));
914 }
915 }
916
917 typedef std::list< std::vector<ScenePoint2D> > Contours;
918
919 Contours contours;
920 UnionOfRectangles::Apply(contours, rectangles);
921
922 chains.reserve(contours.size());
923
924 for (Contours::const_iterator it = contours.begin(); it != contours.end(); ++it)
925 {
926 chains.push_back(*it);
927 }
928
919 #else 929 #else
920 // this will contain the intersection of the polygon slab with 930 // this will contain the intersection of the polygon slab with
921 // the cutting plane, projected on the cutting plane coord system 931 // the cutting plane, projected on the cutting plane coord system
922 // (that yields a rectangle) + the Z coordinate of the polygon 932 // (that yields a rectangle) + the Z coordinate of the polygon
923 // (this is required to group polygons with the same Z later) 933 // (this is required to group polygons with the same Z later)