# HG changeset patch # User Sebastien Jodogne # Date 1510938091 -3600 # Node ID 35c2b85836ce3b1cde4495db5cea279ab1462c52 # Parent 3e6163a53b1660c4c72c622d2935a6d128f0acdc fix rtstruct projections diff -r 3e6163a53b16 -r 35c2b85836ce Framework/Layers/DicomStructureSetRendererFactory.cpp --- a/Framework/Layers/DicomStructureSetRendererFactory.cpp Thu Nov 16 15:00:45 2017 +0100 +++ b/Framework/Layers/DicomStructureSetRendererFactory.cpp Fri Nov 17 18:01:31 2017 +0100 @@ -26,57 +26,85 @@ class DicomStructureSetRendererFactory::Renderer : public ILayerRenderer { private: - DicomStructureSet& structureSet_; + class Structure + { + private: + bool visible_; + uint8_t red_; + uint8_t green_; + uint8_t blue_; + std::string name_; + std::vector< std::vector > polygons_; + + public: + Structure(DicomStructureSet& structureSet, + const CoordinateSystem3D& slice, + size_t index) : + name_(structureSet.GetStructureName(index)) + { + structureSet.GetStructureColor(red_, green_, blue_, index); + visible_ = structureSet.ProjectStructure(polygons_, index, slice); + } + + void Render(CairoContext& context) + { + if (visible_) + { + cairo_t* cr = context.GetObject(); + + context.SetSourceColor(red_, green_, blue_); + + for (size_t i = 0; i < polygons_.size(); i++) + { + cairo_move_to(cr, polygons_[i][0].first, polygons_[i][0].second); + + for (size_t j = 1; j < polygons_[i].size(); j++) + { + cairo_line_to(cr, polygons_[i][j].first, polygons_[i][j].second); + } + + cairo_line_to(cr, polygons_[i][0].first, polygons_[i][0].second); + cairo_stroke(cr); + } + } + } + }; + + typedef std::list Structures; + CoordinateSystem3D slice_; - bool visible_; - + Structures structures_; + public: Renderer(DicomStructureSet& structureSet, const CoordinateSystem3D& slice) : - structureSet_(structureSet), - slice_(slice), - visible_(true) + slice_(slice) { + for (size_t k = 0; k < structureSet.GetStructureCount(); k++) + { + structures_.push_back(new Structure(structureSet, slice, k)); + } + } + + virtual ~Renderer() + { + for (Structures::iterator it = structures_.begin(); + it != structures_.end(); ++it) + { + delete *it; + } } virtual bool RenderLayer(CairoContext& context, const ViewportGeometry& view) { - if (visible_) - { - cairo_set_line_width(context.GetObject(), 3.0f / view.GetZoom()); - - cairo_t* cr = context.GetObject(); - - for (size_t k = 0; k < structureSet_.GetStructureCount(); k++) - { - /*if (structureSet_.GetStructureName(k) != "CORD") - { - continue; - }*/ - - std::vector< std::vector > polygons; + cairo_set_line_width(context.GetObject(), 2.0f / view.GetZoom()); - if (structureSet_.ProjectStructure(polygons, k, slice_)) - { - uint8_t red, green, blue; - structureSet_.GetStructureColor(red, green, blue, k); - context.SetSourceColor(red, green, blue); - - for (size_t i = 0; i < polygons.size(); i++) - { - cairo_move_to(cr, polygons[i][0].first, polygons[i][0].second); - - for (size_t j = 1; j < polygons[i].size(); j++) - { - cairo_line_to(cr, polygons[i][j].first, polygons[i][j].second); - } - - cairo_line_to(cr, polygons[i][0].first, polygons[i][0].second); - cairo_stroke(cr); - } - } - } + for (Structures::const_iterator it = structures_.begin(); + it != structures_.end(); ++it) + { + assert(*it != NULL); + (*it)->Render(context); } return true; @@ -89,9 +117,8 @@ virtual void SetLayerStyle(const RenderStyle& style) { - visible_ = style.visible_; } - + virtual bool IsFullQuality() { return true; diff -r 3e6163a53b16 -r 35c2b85836ce Framework/Toolbox/DicomStructureSet.cpp --- a/Framework/Toolbox/DicomStructureSet.cpp Thu Nov 16 15:00:45 2017 +0100 +++ b/Framework/Toolbox/DicomStructureSet.cpp Fri Nov 17 18:01:31 2017 +0100 @@ -213,8 +213,10 @@ double& y2, const CoordinateSystem3D& slice) const { + // TODO Optimize this method using a sweep-line algorithm for polygons + if (!hasSlice_ || - points_.empty()) + points_.size() <= 1) { return false; } @@ -226,42 +228,110 @@ if (GeometryToolbox::IsParallelOrOpposite (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) { - if (y >= extent_.GetY1() && - y <= extent_.GetY2()) + if (y < extent_.GetY1() || + y > extent_.GetY2()) + { + // The polygon does not intersect the input slice + return false; + } + + bool isFirst = true; + double xmin = std::numeric_limits::infinity(); + double xmax = -std::numeric_limits::infinity(); + + double prevX, prevY; + geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); + + for (size_t i = 0; i < points_.size(); i++) { - Vector p1 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX1(), y) + + // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py + double curX, curY; + geometry_.ProjectPoint(curX, curY, points_[i]); + + if ((prevY < y && curY > y) || + (prevY > y && curY < y)) + { + double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); + xmin = std::min(xmin, p); + xmax = std::max(xmax, p); + isFirst = false; + } + + prevX = curX; + prevY = curY; + } + + if (isFirst) + { + return false; + } + else + { + Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + sliceThickness_ / 2.0 * geometry_.GetNormal()); - Vector p2 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX2(), y) - + Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - sliceThickness_ / 2.0 * geometry_.GetNormal()); - + slice.ProjectPoint(x1, y1, p1); slice.ProjectPoint(x2, y2, p2); return true; } - else - { - return false; - } } else if (GeometryToolbox::IsParallelOrOpposite (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) { - if (x >= extent_.GetX1() && - x <= extent_.GetX2()) + if (x < extent_.GetX1() || + x > extent_.GetX2()) + { + return false; + } + + bool isFirst = true; + double ymin = std::numeric_limits::infinity(); + double ymax = -std::numeric_limits::infinity(); + + double prevX, prevY; + geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); + + for (size_t i = 0; i < points_.size(); i++) { - Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY1()) + + // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py + double curX, curY; + geometry_.ProjectPoint(curX, curY, points_[i]); + + if ((prevX < x && curX > x) || + (prevX > x && curX < x)) + { + double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); + ymin = std::min(ymin, p); + ymax = std::max(ymax, p); + isFirst = false; + } + + prevX = curX; + prevY = curY; + } + + if (isFirst) + { + return false; + } + else + { + Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + sliceThickness_ / 2.0 * geometry_.GetNormal()); - Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY2()) - + Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - sliceThickness_ / 2.0 * geometry_.GetNormal()); slice.ProjectPoint(x1, y1, p1); slice.ProjectPoint(x2, y2, p2); + + // TODO WHY THIS??? + y1 = -y1; + y2 = -y2; + return true; } - else - { - return false; - } } else { @@ -404,6 +474,7 @@ } Polygon polygon(sopInstanceUid); + polygon.Reserve(countPoints); for (size_t k = 0; k < countPoints; k++) { @@ -671,6 +742,7 @@ else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) { +#if 1 // Sagittal or coronal projection std::vector projected; @@ -698,6 +770,22 @@ polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); } } +#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 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; } diff -r 3e6163a53b16 -r 35c2b85836ce Framework/Toolbox/DicomStructureSet.h --- a/Framework/Toolbox/DicomStructureSet.h Thu Nov 16 15:00:45 2017 +0100 +++ b/Framework/Toolbox/DicomStructureSet.h Fri Nov 17 18:01:31 2017 +0100 @@ -57,7 +57,7 @@ typedef std::map ReferencedSlices; - typedef std::list Points; + typedef std::vector Points; class Polygon { @@ -79,6 +79,11 @@ { } + void Reserve(size_t n) + { + points_.reserve(n); + } + void AddPoint(const Vector& v); bool UpdateReferencedSlice(const ReferencedSlices& slices); diff -r 3e6163a53b16 -r 35c2b85836ce Resources/Computations/IntersectSegmentAndHorizontalLine.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Resources/Computations/IntersectSegmentAndHorizontalLine.py Fri Nov 17 18:01:31 2017 +0100 @@ -0,0 +1,21 @@ +#!/usr/bin/env python + +from sympy import * + +# Intersection between the 2D line segment (prevX,prevY)-(curX,curY) and the +# horizontal line "y = y0" using homogeneous coordinates + +prevX, prevY, curX, curY, y0 = symbols('prevX prevY curX curY y0') + +p1 = Matrix([prevX, prevY, 1]) +p2 = Matrix([curX, curY, 1]) +l1 = p1.cross(p2) + +h1 = Matrix([0, y0, 1]) +h2 = Matrix([1, y0, 1]) +l2 = h1.cross(h2) + +a = l1.cross(l2) + +#pprint(cse(a/a[2], symbols = symbols('a b'))) +pprint(a / a[2]) diff -r 3e6163a53b16 -r 35c2b85836ce Resources/Computations/IntersectSegmentAndVerticalLine.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Resources/Computations/IntersectSegmentAndVerticalLine.py Fri Nov 17 18:01:31 2017 +0100 @@ -0,0 +1,20 @@ +#!/usr/bin/env python + +from sympy import * + +# Intersection between the 2D line segment (prevX,prevY)-(curX,curY) and the +# vertical line "x = x0" using homogeneous coordinates + +prevX, prevY, curX, curY, x0 = symbols('prevX prevY curX curY x0') + +p1 = Matrix([prevX, prevY, 1]) +p2 = Matrix([curX, curY, 1]) +l1 = p1.cross(p2) + +h1 = Matrix([x0, 0, 1]) +h2 = Matrix([x0, 1, 1]) +l2 = h1.cross(h2) + +a = l1.cross(l2) + +pprint(a / a[2])