orthanc-stone: OrthancStone/Sources/Toolbox/DicomStructureSet.cpp comparison

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1892:cdf91ad891a5

estimated geometry of rt-struct

author	Sebastien Jodogne <s.jodogne@gmail.com>
date	Wed, 19 Jan 2022 13:50:28 +0100
parents	6ce81914f7e4
children	90b5e116a5f9

comparison

equal deleted inserted replaced

-:3716d72161d2
+:cdf91ad891a5
 #include "GeometryToolbox.h"
 #include "GenericToolbox.h"
 #include "OrthancDatasets/DicomDatasetReader.h"
+#include "BucketAccumulator2D.h"
 #include <Logging.h>
 #include <OrthancException.h>
 #include <Toolbox.h>
 #if defined(_MSC_VER)
 # include <boost/date_time/posix_time/posix_time.hpp>
 #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 ORTHANC_ENABLE_DCMTK == 1
 #  include "ParsedDicomDataset.h"
 #endif
 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,
 }
 structures_[i].polygons_.push_back(polygon);
 }
 }
+EstimateGeometry();
 #if STONE_TIME_BLOCKING_OPS
 boost::posix_time::ptime timerEnd = boost::posix_time::microsec_clock::universal_time();
 boost::posix_time::time_duration duration = timerEnd - timerStart;
 int64_t durationMs = duration.total_milliseconds();
 LOG(WARNING) << "DicomStructureSet::Setup took " << durationMs << " ms";
 }
 }
 }
-void DicomStructureSet::Test()
+void DicomStructureSet::EstimateGeometry()
 {
-printf("OK\n");
+static const double PI = boost::math::constants::pi<double>();
+BucketAccumulator2D accumulator(0, PI, 9,   /* for acos() */
+-PI, PI, 9, /* for atan() */
+true /* store values */);
+unsigned int countPolygons = 0;
 for (size_t i = 0; i < structures_.size(); i++)
 {
 const Polygons& polygons = structures_[i].polygons_;
 for (Polygons::const_iterator it = polygons.begin(); it != polygons.end(); ++it)
 {
+countPolygons++;
 const Points& points = it->GetPoints();
 if (points.size() >= 3)
 {
-const Vector& a = points[0];
+// Compute the normal of the polygon using 3 successive points
-const Vector& b = points[1];
+Vector normal;
-const Vector& c = points[2];
+bool valid = false;
-Vector n;
-LinearAlgebra::CrossProduct(n, b - a, c - a);
+for (size_t i = 0; i + 2 < points.size(); i++)
-LinearAlgebra::NormalizeVector(n);
+{
-if (n[2] < 0)
+const Vector& a = points[i];
-n = -n;
+const Vector& b = points[i + 1];
-//LinearAlgebra::Print(n);
+const Vector& c = points[i + 2];
+LinearAlgebra::CrossProduct(normal, b - a, c - a);  // (*)
-// https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates#Vector_fields_2
+LinearAlgebra::NormalizeVector(normal);
-double r = 1.0;
-double theta = acos(n[2]);
+if (LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal), 1.0))  // (**)
-double phi = atan(n[1]);
+{
-printf("%.02f %.02f %.02f => ", n[0], n[1], n[2]);
+valid = true;
-printf("%.02f %.02f =>", theta, phi);
+break;
-printf("%.02f %.02f %.02f\n", r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta));
+}
 }
-}
+if (valid)
+{
+// Check that all the points of the polygon lie in the plane defined by the normal
+double d1 = GeometryToolbox::ProjectAlongNormal(points[0], normal);
+for (size_t i = 1; i < points.size(); i++)
+{
+double d2 = GeometryToolbox::ProjectAlongNormal(points[i], normal);
+if (!LinearAlgebra::IsNear(d1, d2))
+{
+valid = false;
+break;
+}
+}
+}
+if (valid)
+{
+if (normal[2] < 0)
+{
+normal = -normal;
+}
+// The normal is a non-zero unit vector because of (*) and (**), so "r == 1"
+// https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates#Vector_fields_2
+const double theta = acos(normal[2]);
+const double phi = atan(normal[1]);
+accumulator.AddValue(theta, phi);
+}
+}
+}
+}
+size_t bestX, bestY;
+accumulator.FindBestBucket(bestX, bestY);
+if (accumulator.GetBucketContentSize(bestX, bestY) > 0)
+{
+double normalTheta, normalPhi;
+accumulator.ComputeBestMedian(normalTheta, normalPhi);
+// Back to (x,y,z) normal, taking "r == 1"
+// https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates#Vector_fields_2
+double sinTheta = sin(normalTheta);
+estimatedNormal_ = LinearAlgebra::CreateVector(sinTheta * cos(normalPhi), sinTheta * sin(normalPhi), cos(normalTheta));
+}
+std::vector<double> polygonsProjection;
+polygonsProjection.reserve(countPolygons);
+for (size_t i = 0; i < structures_.size(); i++)
+{
+const Polygons& polygons = structures_[i].polygons_;
+for (Polygons::const_iterator it = polygons.begin(); it != polygons.end(); ++it)
+{
+const Points& points = it->GetPoints();
+polygonsProjection.push_back(GeometryToolbox::ProjectAlongNormal(points[0], estimatedNormal_));
+}
+}
+std::sort(polygonsProjection.begin(), polygonsProjection.end());
+std::vector<double> deltas;
+deltas.reserve(polygonsProjection.size());
+for (size_t i = 0; i + 1 < polygonsProjection.size(); i++)
+{
+if (!LinearAlgebra::IsNear(polygonsProjection[i], polygonsProjection[i + 1]))
+{
+assert(polygonsProjection[i + 1] > polygonsProjection[i]);
+deltas.push_back(polygonsProjection[i + 1] - polygonsProjection[i]);
+}
+}
+if (deltas.empty())
+{
+estimatedSliceThickness_ = 1;
+}
+else
+{
+estimatedSliceThickness_ = LinearAlgebra::ComputeMedian(deltas);
 }
 }
 }

Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1892:cdf91ad891a5