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

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 2157:917e40af6b45 default

fix rendering of RT-STRUCT

author	Sebastien Jodogne <s.jodogne@gmail.com>
date	Fri, 27 Sep 2024 22:29:51 +0200
parents	f8cbfd7175c3
children

comparison

equal deleted inserted replaced

-:340bde744884
+:917e40af6b45
 }
 }
 }
-static void AddSegmentIfIntersection(Extent2D& extent,
+void DicomStructureSet::Polygon::Project(std::list<Extent2D>& target,
-const CoordinateSystem3D& cuttingPlane,
-const Vector& p1,
-const Vector& p2,
-double originDistance)
-{
-// Does this segment intersects the cutting plane?
-double d1 = cuttingPlane.ProjectAlongNormal(p1);
-double d2 = cuttingPlane.ProjectAlongNormal(p2);
-if ((d1 < originDistance && d2 > originDistance) ||
-(d1 > originDistance && d2 < originDistance))
-{
-// This is an intersection: Add the segment
-double x, y;
-cuttingPlane.ProjectPoint2(x, y, p1);
-extent.AddPoint(x, y);
-cuttingPlane.ProjectPoint2(x, y, p2);
-extent.AddPoint(x, y);
-}
-}
-bool DicomStructureSet::Polygon::Project(double& x1,
-double& y1,
-double& x2,
-double& y2,
 const CoordinateSystem3D& cuttingPlane,
 const Vector& estimatedNormal,
 double estimatedSliceThickness) const
 {
-if (points_.size() <= 1)
+CoordinateSystem3D geometry;
-{
-return false;
-}
-Vector normal = estimatedNormal;
 double thickness = estimatedSliceThickness;
+/**
+* 1. Estimate the 3D plane associated with this polygon.
+**/
 if (hasSlice_)
 {
-normal = geometry_.GetNormal();
+// The exact geometry is known for this slice
+geometry = geometry_;
 thickness = sliceThickness_;
 }
+else if (points_.size() < 2)
+{
+return;
+}
+else
+{
+// Estimate the geometry
+Vector origin = points_[0];
+Vector axisX;
+bool found = false;
+for (size_t i = 1; i < points_.size(); i++)
+{
+axisX = points_[1] - origin;
+if (boost::numeric::ublas::norm_2(axisX) > 10.0 * std::numeric_limits<double>::epsilon())
+{
+found = true;
+break;
+}
+}
+if (!found)
+{
+return;  // The polygon is too small to extract a reliable geometry out of it
+}
+LinearAlgebra::NormalizeVector(axisX);
+Vector axisY;
+LinearAlgebra::CrossProduct(axisY, axisX, estimatedNormal);
+geometry = CoordinateSystem3D(origin, axisX, axisY);
+}
+/**
+* 2. Project the 3D cutting plane as a 2D line onto the polygon plane.
+**/
+double cuttingX1, cuttingY1, cuttingX2, cuttingY2;
+geometry.ProjectPoint(cuttingX1, cuttingY1, cuttingPlane.GetOrigin());
 bool isOpposite;
-if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) &&
+if (GeometryToolbox::IsParallelOrOpposite(isOpposite, geometry.GetNormal(), cuttingPlane.GetAxisX()))
-!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY()))
+{
-{
+geometry.ProjectPoint(cuttingX2, cuttingY2, cuttingPlane.GetOrigin() + cuttingPlane.GetAxisY());
-return false;
+}
-}
+else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, geometry.GetNormal(), cuttingPlane.GetAxisY()))
+{
-const double d = cuttingPlane.ProjectAlongNormal(cuttingPlane.GetOrigin());
+geometry.ProjectPoint(cuttingX2, cuttingY2, cuttingPlane.GetOrigin() + cuttingPlane.GetAxisX());
-Extent2D extent;
-AddSegmentIfIntersection(extent, cuttingPlane, points_[points_.size() - 1], points_[0], d);
-for (size_t i = 1; i < points_.size(); i++)
-{
-AddSegmentIfIntersection(extent, cuttingPlane, points_[i - 1], points_[i], d);
-}
-if (extent.GetWidth() > 0 ||
-extent.GetHeight() > 0)
-{
-// Let's convert them to 3D world geometry to add the slice thickness
-Vector p1 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX1(), extent.GetY1()) +
-thickness / 2.0 * normal);
-Vector p2 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX2(), extent.GetY2()) -
-thickness / 2.0 * normal);
-// Then back to the cutting plane geometry
-cuttingPlane.ProjectPoint2(x1, y1, p1);
-cuttingPlane.ProjectPoint2(x2, y2, p2);
-return true;
 }
 else
 {
-return false;
+return;
+}
+/**
+* 3. Compute the intersection of the 2D cutting line with the polygon.
+**/
+// Initialize the projection of a point onto a line:
+// https://stackoverflow.com/a/64330724
+const double abx = cuttingX2 - cuttingX1;
+const double aby = cuttingY2 - cuttingY1;
+const double denominator = abx * abx + aby * aby;
+if (LinearAlgebra::IsCloseToZero(denominator))
+{
+return;  // Should never happen
+}
+std::vector<double> intersections;
+intersections.reserve(points_.size());
+for (size_t i = 0; i < points_.size(); i++)
+{
+double segmentX1, segmentY1, segmentX2, segmentY2;
+geometry.ProjectPoint(segmentX1, segmentY1, points_[i]);
+geometry.ProjectPoint(segmentX2, segmentY2, points_[(i + 1) % points_.size()]);
+double x, y;
+if (GeometryToolbox::IntersectLineAndSegment(x, y, cuttingX1, cuttingY1, cuttingX2, cuttingY2,
+segmentX1, segmentY1, segmentX2, segmentY2))
+{
+// For each polygon segment that intersect the cutting line,
+// register its offset over the cutting line
+const double acx = x - cuttingX1;
+const double acy = y - cuttingY1;
+intersections.push_back((abx * acx + aby * acy) / denominator);
+}
+}
+/**
+* 4. Sort the intersection offsets, then generates one 2D rectangle on the
+* cutting plane from each pair of successive intersections.
+**/
+std::sort(intersections.begin(), intersections.end());
+if (intersections.size() % 2 == 1)
+{
+return;  // Should never happen
+}
+for (size_t i = 0; i < intersections.size(); i += 2)
+{
+Vector p1, p2;
+{
+// Let's convert them to 3D world geometry to add the slice thickness
+const double x1 = cuttingX1 + intersections[i] * abx;
+const double y1 = cuttingY1 + intersections[i] * aby;
+const double x2 = cuttingX1 + intersections[i + 1] * abx;
+const double y2 = cuttingY1 + intersections[i + 1] * aby;
+p1 = (geometry.MapSliceToWorldCoordinates(x1, y1) + thickness / 2.0 * geometry.GetNormal());
+p2 = (geometry.MapSliceToWorldCoordinates(x2, y2) - thickness / 2.0 * geometry.GetNormal());
+}
+{
+// Then back to the cutting plane geometry
+double x1, y1, x2, y2;
+cuttingPlane.ProjectPoint2(x1, y1, p1);
+cuttingPlane.ProjectPoint2(x2, y2, p2);
+target.push_back(Extent2D(x1, y1, x2, y2));
+}
 }
 }
 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const
 std::vector<BoostPolygon> projected;
 for (Polygons::const_iterator polygon = structure.polygons_.begin();
 polygon != structure.polygons_.end(); ++polygon)
 {
-double x1, y1, x2, y2;
+std::list<Extent2D> rectangles;
+polygon->Project(rectangles, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness());
-if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()))
-{
+for (std::list<Extent2D>::const_iterator it = rectangles.begin(); it != rectangles.end(); ++it)
-projected.push_back(CreateRectangle(x1, y1, x2, y2));
+{
+projected.push_back(CreateRectangle(it->GetX1(), it->GetY1(), it->GetX2(), it->GetY2()));
 }
 }
 BoostMultiPolygon merged;
 Union(merged, projected);
 std::list<Extent2D> rectangles;
 for (Polygons::const_iterator polygon = structure.polygons_.begin();
 polygon != structure.polygons_.end(); ++polygon)
 {
-double x1, y1, x2, y2;
+polygon->Project(rectangles, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness());
-if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()))
-{
-rectangles.push_back(Extent2D(x1, y1, x2, y2));
-}
 }
 typedef std::list< std::vector<ScenePoint2D> >  Contours;
 Contours contours;

Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 2157:917e40af6b45 default