Mercurial > hg > orthanc-stone
diff 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 |
line wrap: on
line diff
--- a/OrthancStone/Sources/Toolbox/DicomStructureSet.cpp Wed Jan 19 14:51:55 2022 +0100 +++ b/OrthancStone/Sources/Toolbox/DicomStructureSet.cpp Wed Jan 26 17:19:37 2022 +0100 @@ -21,6 +21,8 @@ **/ +#define USE_BOOST_UNION_FOR_POLYGONS 1 + #include "DicomStructureSet.h" #include "DicomStructureSetUtils.h" // TODO REMOVE @@ -42,13 +44,20 @@ # include <boost/date_time/posix_time/posix_time.hpp> #endif +#if !defined(USE_BOOST_UNION_FOR_POLYGONS) +# error Macro USE_BOOST_UNION_FOR_POLYGONS must be defined +#endif + #include <limits> #include <stdio.h> #include <boost/math/constants/constants.hpp> -#include <boost/geometry.hpp> -#include <boost/geometry/geometries/point_xy.hpp> -#include <boost/geometry/geometries/polygon.hpp> -#include <boost/geometry/multi/geometries/multi_polygon.hpp> + +#if USE_BOOST_UNION_FOR_POLYGONS == 1 +# include <boost/geometry.hpp> +# include <boost/geometry/geometries/point_xy.hpp> +# include <boost/geometry/geometries/polygon.hpp> +# include <boost/geometry/multi/geometries/multi_polygon.hpp> +#endif #if defined(_MSC_VER) # pragma warning(pop) @@ -59,11 +68,12 @@ #endif +#if USE_BOOST_UNION_FOR_POLYGONS == 1 + typedef boost::geometry::model::d2::point_xy<double> BoostPoint; typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; - static void Union(BoostMultiPolygon& output, std::vector<BoostPolygon>& input) { @@ -94,8 +104,6 @@ } } -#if USE_BOOST_UNION_FOR_POLYGONS == 1 - static BoostPolygon CreateRectangle(float x1, float y1, float x2, float y2) { @@ -122,12 +130,14 @@ return rect; } - bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1, const std::pair<RtStructRectangleInSlab, double>& r2) + bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1, + const std::pair<RtStructRectangleInSlab, double>& r2) { return r1.second < r2.second; } - bool CompareSlabsY(const RtStructRectanglesInSlab& r1, const RtStructRectanglesInSlab& r2) + bool CompareSlabsY(const RtStructRectanglesInSlab& r1, + const RtStructRectanglesInSlab& r2) { if ((r1.size() == 0) || (r2.size() == 0)) return false; @@ -829,22 +839,13 @@ } } - bool DicomStructureSet::ProjectStructure( -#if USE_BOOST_UNION_FOR_POLYGONS == 1 - std::vector< std::vector<ScenePoint2D> >& polygons, -#else - std::vector< std::pair<ScenePoint2D, ScenePoint2D> >& segments, -#endif - const Structure& structure, - const CoordinateSystem3D& sourceSlice) const + bool DicomStructureSet::ProjectStructure(std::vector< std::vector<ScenePoint2D> >& chains, + const Structure& structure, + const CoordinateSystem3D& sourceSlice) const { const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); -#if USE_BOOST_UNION_FOR_POLYGONS == 1 - polygons.clear(); -#else - segments.clear(); -#endif + chains.clear(); Vector normal = GetNormal(); @@ -853,51 +854,30 @@ { // This is an axial projection + chains.reserve(structure.polygons_.size()); + for (Polygons::const_iterator polygon = structure.polygons_.begin(); polygon != structure.polygons_.end(); ++polygon) { - if (polygon->IsOnSlice(slice)) + const Points& points = polygon->GetPoints(); + + if (polygon->IsOnSlice(slice) && + !points.empty()) { -#if USE_BOOST_UNION_FOR_POLYGONS == 1 - polygons.push_back(std::vector<ScenePoint2D>()); + chains.push_back(std::vector<ScenePoint2D>()); + chains.back().reserve(points.size() + 1); - for (Points::const_iterator p = polygon->GetPoints().begin(); - p != polygon->GetPoints().end(); ++p) + for (Points::const_iterator p = points.begin(); + p != points.end(); ++p) { double x, y; slice.ProjectPoint2(x, y, *p); - polygons.back().push_back(ScenePoint2D(x, y)); + chains.back().push_back(ScenePoint2D(x, y)); } -#else - // we need to add all the segments corresponding to this polygon - const std::vector<Vector>& points3D = polygon->GetPoints(); - if (points3D.size() >= 3) - { - ScenePoint2D prev2D; - { - Vector prev = points3D[0]; - double prevX, prevY; - slice.ProjectPoint2(prevX, prevY, prev); - prev2D = ScenePoint2D(prevX, prevY); - } - size_t pointCount = points3D.size(); - for (size_t ipt = 1; ipt < pointCount; ++ipt) - { - Vector next = points3D[ipt]; - double nextX, nextY; - slice.ProjectPoint2(nextX, nextY, next); - ScenePoint2D next2D(nextX, nextY); - segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(prev2D, next2D)); - prev2D = next2D; - } - } - else - { - LOG(ERROR) << "Contour with less than 3 points!"; - // !!! - } -#endif + double x0, y0; + slice.ProjectPoint2(x0, y0, points.front()); + chains.back().push_back(ScenePoint2D(x0, y0)); } } @@ -906,7 +886,6 @@ else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) { -#if 1 // Sagittal or coronal projection #if USE_BOOST_UNION_FOR_POLYGONS == 1 @@ -997,40 +976,33 @@ // now we need to sort the slabs in increasing Y order (see ConvertListOfSlabsToSegments) std::sort(rectanglesForEachSlab.begin(), rectanglesForEachSlab.end(), CompareSlabsY); + std::vector< std::pair<ScenePoint2D, ScenePoint2D> > segments; ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, projected.size()); + + chains.resize(segments.size()); + for (size_t i = 0; i < segments.size(); i++) + { + chains[i].resize(2); + chains[i][0] = segments[i].first; + chains[i][1] = segments[i].second; + } + #else BoostMultiPolygon merged; Union(merged, projected); - polygons.resize(merged.size()); + chains.resize(merged.size()); for (size_t i = 0; i < merged.size(); i++) { const std::vector<BoostPoint>& outer = merged[i].outer(); - polygons[i].resize(outer.size()); + chains[i].resize(outer.size()); for (size_t j = 0; j < outer.size(); j++) { - polygons[i][j] = ScenePoint2D(outer[j].x(), outer[j].y()); + chains[i][j] = ScenePoint2D(outer[j].x(), outer[j].y()); } } #endif - -#else - for (Polygons::iterator polygon = structure.polygons_.begin(); - polygon != structure.polygons_.end(); ++polygon) - { - double x1, y1, x2, y2; - if (polygon->Project(x1, y1, x2, y2, slice)) - { - std::vector<ScenePoint2D> p(4); - p[0] = std::make_pair(x1, y1); - p[1] = std::make_pair(x2, y1); - p[2] = std::make_pair(x2, y2); - p[3] = std::make_pair(x1, y2); - polygons.push_back(p); - } - } -#endif return true; } @@ -1046,38 +1018,15 @@ size_t structureIndex, const Color& color) const { -#if USE_BOOST_UNION_FOR_POLYGONS == 1 - std::vector< std::vector<ScenePoint2D> > polygons; - if (ProjectStructure(polygons, structureIndex, plane)) + std::vector< std::vector<ScenePoint2D> > chains; + + if (ProjectStructure(chains, structureIndex, plane)) { - for (size_t j = 0; j < polygons.size(); j++) + for (size_t j = 0; j < chains.size(); j++) { - std::vector<ScenePoint2D> chain; - chain.reserve(polygons[j].size()); - - for (size_t k = 0; k < polygons[j].size(); k++) - { - chain.push_back(ScenePoint2D(polygons[j][k].x, polygons[j][k].y)); - } - - layer.AddChain(chain, true, color.GetRed(), color.GetGreen(), color.GetBlue()); + layer.AddChain(chains[j], false, color.GetRed(), color.GetGreen(), color.GetBlue()); } } - -#else - std::vector< std::pair<ScenePoint2D, ScenePoint2D> > segments; - - if (ProjectStructure(segments, structureIndex, plane)) - { - for (size_t j = 0; j < segments.size(); j++) - { - std::vector<ScenePoint2D> chain(2); - chain[0] = ScenePoint2D(segments[j].first.GetX(), segments[j].first.GetY()); - chain[1] = ScenePoint2D(segments[j].second.GetX(), segments[j].second.GetY()); - layer.AddChain(chain, false, color.GetRed(), color.GetGreen(), color.GetBlue()); - } - } -#endif }