Mercurial > hg > orthanc-stone
changeset 998:38b6bb0bdd72
added a new set of classes that correctly handle non-convex polygons (not
used yet because of limitations in coordinates computing): DicomStructure2,
DicomStructureSet2, DicomStructurePolygon2, DicomStructureSetSlicer2. Too
many shortcuts have been taken when computing the actual position.
| author | Benjamin Golinvaux <bgo@osimis.io> |
|---|---|
| date | Fri, 20 Sep 2019 11:58:00 +0200 |
| parents | 9893fa8cd7a6 |
| children | 2d69b8bee484 |
| files | Framework/Loaders/DicomStructureSetLoader2.cpp Framework/Loaders/DicomStructureSetLoader2.h Framework/Toolbox/DicomStructure2.cpp Framework/Toolbox/DicomStructure2.h Framework/Toolbox/DicomStructurePolygon2.cpp Framework/Toolbox/DicomStructurePolygon2.h Framework/Toolbox/DicomStructureSet2.cpp Framework/Toolbox/DicomStructureSet2.h Framework/Toolbox/DicomStructureSetUtils.cpp Framework/Toolbox/DicomStructureSetUtils.h Framework/Volumes/DicomStructureSetSlicer2.cpp Framework/Volumes/DicomStructureSetSlicer2.h |
| diffstat | 12 files changed, 1774 insertions(+), 806 deletions(-) [+] |
line wrap: on
line diff
--- a/Framework/Loaders/DicomStructureSetLoader2.cpp Sat Sep 14 17:27:41 2019 +0200 +++ b/Framework/Loaders/DicomStructureSetLoader2.cpp Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,119 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#include "DicomStructureSetLoader2.h" + +#include "../Messages/IObservable.h" +#include "../Oracle/IOracle.h" +#include "../Oracle/OracleCommandExceptionMessage.h" + +namespace OrthancStone +{ + + DicomStructureSetLoader2::DicomStructureSetLoader2( + DicomStructureSet2& structureSet + , IOracle& oracle + , IObservable& oracleObservable) + : IObserver(oracleObservable.GetBroker()) + , IObservable(oracleObservable.GetBroker()) + , structureSet_(structureSet) + , oracle_(oracle) + , oracleObservable_(oracleObservable) + , structuresReady_(false) + { + LOG(TRACE) << "DicomStructureSetLoader2(" << std::hex << this << std::dec << ")::DicomStructureSetLoader2()"; + + oracleObservable.RegisterObserverCallback( + new Callable<DicomStructureSetLoader2, OrthancRestApiCommand::SuccessMessage> + (*this, &DicomStructureSetLoader2::HandleSuccessMessage)); + + oracleObservable.RegisterObserverCallback( + new Callable<DicomStructureSetLoader2, OracleCommandExceptionMessage> + (*this, &DicomStructureSetLoader2::HandleExceptionMessage)); + } + + DicomStructureSetLoader2::~DicomStructureSetLoader2() + { + LOG(TRACE) << "DicomStructureSetLoader2(" << std::hex << this << std::dec << ")::~DicomStructureSetLoader2()"; + oracleObservable_.Unregister(this); + } + + void DicomStructureSetLoader2::LoadInstanceFromString(const std::string& body) + { + OrthancPlugins::FullOrthancDataset dicom(body); + //loader.content_.reset(new DicomStructureSet(dicom)); + structureSet_.Clear(); + structureSet_.SetContents(dicom); + SetStructuresReady(); + } + + void DicomStructureSetLoader2::HandleSuccessMessage(const OrthancRestApiCommand::SuccessMessage& message) + { + const std::string& body = message.GetAnswer(); + LoadInstanceFromString(body); + } + + void DicomStructureSetLoader2::HandleExceptionMessage(const OracleCommandExceptionMessage& message) + { + LOG(ERROR) << "DicomStructureSetLoader2::HandleExceptionMessage: error when trying to load data. " + << "Error: " << message.GetException().What() << " Details: " + << message.GetException().GetDetails(); + } + + void DicomStructureSetLoader2::LoadInstance(const std::string& instanceId) + { + std::auto_ptr<OrthancRestApiCommand> command(new OrthancRestApiCommand); + command->SetHttpHeader("Accept-Encoding", "gzip"); + + std::string uri = "/instances/" + instanceId + "/tags?ignore-length=3006-0050"; + + command->SetUri(uri); + oracle_.Schedule(*this, command.release()); + } + + void DicomStructureSetLoader2::SetStructuresReady() + { + structuresReady_ = true; + } + + bool DicomStructureSetLoader2::AreStructuresReady() const + { + return structuresReady_; + } + + /* + + void LoaderStateMachine::HandleExceptionMessage(const OracleCommandExceptionMessage& message) + { + LOG(ERROR) << "LoaderStateMachine::HandleExceptionMessage: error in the state machine, stopping all processing"; + LOG(ERROR) << "Error: " << message.GetException().What() << " Details: " << + message.GetException().GetDetails(); + Clear(); + } + + LoaderStateMachine::~LoaderStateMachine() + { + Clear(); + } + + + */ + +}
--- a/Framework/Loaders/DicomStructureSetLoader2.h Sat Sep 14 17:27:41 2019 +0200 +++ b/Framework/Loaders/DicomStructureSetLoader2.h Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,80 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#pragma once + +#include "../Toolbox/DicomStructureSet2.h" +#include "../Messages/IMessage.h" +#include "../Messages/IObserver.h" +#include "../Messages/IObservable.h" +#include "../Oracle/OrthancRestApiCommand.h" + +#include <boost/noncopyable.hpp> + +namespace OrthancStone +{ + class IOracle; + class IObservable; + class OrthancRestApiCommand; + class OracleCommandExceptionMessage; + + class DicomStructureSetLoader2 : public IObserver, public IObservable + { + public: + ORTHANC_STONE_DEFINE_ORIGIN_MESSAGE(__FILE__, __LINE__, StructuresReady, DicomStructureSetLoader2); + + /** + Warning: the structureSet, oracle and oracleObservable objects must live + at least as long as this object (TODO: shared_ptr?) + */ + DicomStructureSetLoader2(DicomStructureSet2& structureSet, IOracle& oracle, IObservable& oracleObservable); + + ~DicomStructureSetLoader2(); + + void LoadInstance(const std::string& instanceId); + + /** Internal use */ + void LoadInstanceFromString(const std::string& body); + + void SetStructuresReady(); + bool AreStructuresReady() const; + + private: + /** + Called back by the oracle when data is ready! + */ + void HandleSuccessMessage(const OrthancRestApiCommand::SuccessMessage& message); + + /** + Called back by the oracle when shit hits the fan + */ + void HandleExceptionMessage(const OracleCommandExceptionMessage& message); + + /** + The structure set that will be (cleared and) filled with data from the + loader + */ + DicomStructureSet2& structureSet_; + + IOracle& oracle_; + IObservable& oracleObservable_; + bool structuresReady_; + }; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/DicomStructure2.cpp Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,287 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#include "DicomStructure2.h" + +#include "../Toolbox/GeometryToolbox.h" +#include "../Toolbox/DisjointDataSet.h" + +namespace OrthancStone +{ + // see header + //void DicomStructure2::ComputeNormal() + //{ + // try + // { + // if (polygons_.size() > 0) + // { + + // // TODO: check all polygons are OK + // const DicomStructurePolygon2 polygon = polygons_[0]; + // $$$$$$$$$$$$$$$$$ + // state_ = NormalComputed; + // } + // else + // { + // // bogus! no polygons. Let's assign a "nothing here" value + // LinearAlgebra::AssignVector(normal_, 0, 0, 0); + // state_ = Invalid; + // } + // } + // catch (const Orthanc::OrthancException& e) + // { + // state_ = Invalid; + // if (e.HasDetails()) + // { + // LOG(ERROR) << "OrthancException in ComputeNormal: " << e.What() << " Details: " << e.GetDetails(); + // } + // else + // { + // LOG(ERROR) << "OrthancException in ComputeNormal: " << e.What(); + // } + // throw; + // } + // catch (const std::exception& e) + // { + // state_ = Invalid; + // LOG(ERROR) << "std::exception in ComputeNormal: " << e.what(); + // throw; + // } + // catch (...) + // { + // state_ = Invalid; + // LOG(ERROR) << "Unknown exception in ComputeNormal"; + // throw; + // } + //} + + void DicomStructure2::ComputeSliceThickness() + { + if (state_ != NormalComputed) + { + LOG(ERROR) << "DicomStructure2::ComputeSliceThickness - state must be NormalComputed"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + if (polygons_.size() < 2) + { + // cannot compute thickness if there are not at least 2 slabs (contours) + sliceThickness_ = 1.0; + state_ = Invalid; + } + else + { + // normal can be (1,0,0), (0,1,0) or (0,0,1), nothing else. + // these can be compared with == (exact double representation) + if (normal_[0] == 1) + { + // in a single polygon, all the points have the same X + sliceThickness_ = fabs(polygons_[0].GetPoint(0)[0] - polygons_[1].GetPoint(0)[0]); + } + else if (normal_[1] == 1) + { + // in a single polygon, all the points have the same X + sliceThickness_ = fabs(polygons_[0].GetPoint(0)[1] - polygons_[1].GetPoint(0)[1]); + } + else if (normal_[2] == 1) + { + // in a single polygon, all the points have the same X + sliceThickness_ = fabs(polygons_[0].GetPoint(0)[2] - polygons_[1].GetPoint(0)[2]); + } + else + { + ORTHANC_ASSERT(false); + state_ = Invalid; + } + } + state_ = Valid; + } + + void DicomStructure2::AddPolygon(const DicomStructurePolygon2& polygon) + { + if (state_ != Building) + { + LOG(ERROR) << "DicomStructure2::AddPolygon - can only add polygon while building"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + polygons_.push_back(polygon); + } + + void DicomStructure2::ComputeDependentProperties() + { + if (state_ != Building) + { + LOG(ERROR) << "DicomStructure2::ComputeDependentProperties - can only be called once"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + for (size_t i = 0; i < polygons_.size(); ++i) + { + // "compute" the polygon normal + polygons_[i].ComputeDependentProperties(); + } + if (polygons_.size() > 0) + { + normal_ = polygons_[0].GetNormal(); + state_ = NormalComputed; + } + else + { + LinearAlgebra::AssignVector(normal_, 0, 0, 0); + state_ = Invalid; // THIS MAY HAPPEN !!! (for instance for instance 72c773ac-5059f2c4-2e6a9120-4fd4bca1-45701661 :) ) + } + if (polygons_.size() >= 2) + ComputeSliceThickness(); // this will change state_ from NormalComputed to Valid + } + + OrthancStone::Vector DicomStructure2::GetNormal() const + { + if (state_ != Valid && state_ != Invalid) + { + LOG(ERROR) << "DicomStructure2::GetNormal() -- please call ComputeDependentProperties first."; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + if (state_ == Invalid) + { + LOG(ERROR) << "DicomStructure2::GetNormal() -- The Dicom structure is invalid. The normal is set to 0,0,0"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + return normal_; + } + + const DicomStructurePolygon2* DicomStructure2::GetPolygonClosestToSlice( + const CoordinateSystem3D& plane) const + { + ORTHANC_ASSERT(state_ == Valid); + + // we assume 0,0,1 for now + ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0)); + ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0)); + + for (size_t i = 0; i < polygons_.size(); ++i) + { + const DicomStructurePolygon2& polygon = polygons_[i]; + + // "height" of cutting plane + double cutZ = plane.GetOrigin()[2]; + + if (LinearAlgebra::IsNear( + cutZ, polygon.GetZ(), + sliceThickness_ / 2.0 /* in mm */)) + return &polygon; + } + return NULL; + } + + + bool DicomStructure2::Project(std::vector< std::pair<Point2D, Point2D> > & segments, const CoordinateSystem3D & plane) const + { + segments.clear(); + + Vector normal = GetNormal(); + + size_t totalRectCount = 0; + + // dummy var + bool isOpposite = false; + + // This is an axial projection + if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, plane.GetNormal())) + { + const DicomStructurePolygon2* polygon = GetPolygonClosestToSlice(plane); + if (polygon) + { + polygon->ProjectOnParallelPlane(segments, plane); + } + } + else + { + // let's compute the dot product of the plane normal and the polygons + // normal. + double dot = LinearAlgebra::DotProduct(plane.GetNormal(), normal); + + if (LinearAlgebra::IsNear(dot, 0)) + { + // Coronal or sagittal projection + + // vector of vector of rectangles that will be merged in a single big contour: + + // each polygon slab cut by a perpendicular plane yields 0..* rectangles + std::vector< RtStructRectanglesInSlab > rectanglesForEachSlab; + + for (size_t i = 0; i < polygons_.size(); ++i) + { + // book an entry for this slab + rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); + + // let's compute the intersection between the polygon and the plane + // intersections are in plane coords + std::vector<Point2D> intersections; + + polygons_[i].ProjectOnConstantPlane(intersections, plane); + + // for each pair of intersections, we add a rectangle. + if ((intersections.size() % 2) != 0) + { + LOG(WARNING) << "Odd number of intersections between structure " + << name_ << ", polygon # " << i + << " and plane where X axis is parallel to polygon normal vector"; + } + + size_t numRects = intersections.size() / 2; + + // we keep count of the total number of rects for vector pre-allocations + totalRectCount += numRects; + + for (size_t iRect = 0; iRect < numRects; ++iRect) + { + RtStructRectangleInSlab rectangle; + ORTHANC_ASSERT(LinearAlgebra::IsNear(intersections[2 * iRect].y, intersections[2 * iRect + 1].y)); + ORTHANC_ASSERT((2 * iRect + 1) < intersections.size()); + double x1 = intersections[2 * iRect].x; + double x2 = intersections[2 * iRect + 1].x; + double y1 = intersections[2 * iRect].y - sliceThickness_ * 0.5; + double y2 = intersections[2 * iRect].y + sliceThickness_ * 0.5; + + rectangle.xmin = std::min(x1, x2); + rectangle.xmax = std::max(x1, x2); + rectangle.ymin = std::min(y1, y2); + rectangle.ymax = std::max(y1, y2); + + // TODO: keep them sorted!!!! + + rectanglesForEachSlab.back().push_back(rectangle); + } + } + // now we need to merge all the slabs into a set of polygons (1 or more) + ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, totalRectCount); + } + else + { + // plane is not perpendicular to the polygons + // 180.0 / [Math]::Pi = 57.2957795130823 + double acDot = 57.2957795130823 * acos(dot); + LOG(ERROR) << "DicomStructure2::Project -- cutting plane must be " + << "perpendicular to the contours, but dot product is: " + << dot << " and (180/pi)*acos(dot) = " << acDot; + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + return segments.size() != 0; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/DicomStructure2.h Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,155 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#pragma once + +#include "DicomStructurePolygon2.h" +#include "DicomStructureSetUtils.h" + +namespace OrthancStone +{ + + /* + A structure has a color, a name, a set of slices.. + + Each slice is a polygon. + */ + struct DicomStructure2 + { + DicomStructure2() : + red_(0), green_(0), blue_(0), sliceThickness_(0), state_(Building) {} + + void AddPolygon(const DicomStructurePolygon2& polygon); + + /** + Once all polygons have been added, this method will determine: + - the slice orientation (through the normal vector) + - the spacing between slices (slice thickness) + + it will also set up the info required to efficiently compute plane + intersections later on. + */ + void ComputeDependentProperties(); + + /** + Being given a plane that is PARALLEL to the set of polygon contours, this + returns a pointer to the polygon located at that position (if it is closer + than thickness/2) or NULL if there is none. + + TODO: use sorted vector to improve + + DO NOT STORE THE RETURNED POINTER! + */ + const DicomStructurePolygon2* GetPolygonClosestToSlice(const CoordinateSystem3D& plane) const; + + Vector GetNormal() const; + + Color GetColor() const + { + return Color(red_, green_, blue_); + } + + bool IsValid() const + { + return state_ == Valid; + } + + /** + This method is used to project the 3D structure on a 2D plane. + + A structure is a stack of polygons, representing a volume. + + We need to compute the intersection between this volume and the supplied + cutting plane (the "slice"). This is more than a cutting plane: it is also + a 2D-coordinate system (the plane has axes vectors) + + The cutting plane is always parallel to the plane defined by two of the + world coordinate system axes. + + The result is a set of closed polygons. + + If the cut is parallel to the polygons, we pick the polygon closest to + the slice, project it on the slice and return it in slice coordinates. + + If the cut is perpendicular to the polygons, for each polygon, we compute + the intersection between the cutting plane and the polygon slab (imaginary + volume created by extruding the polygon above and below its plane by + thickness/2) : + - each slab, intersected by the plane, gives a set of 0..* rectangles \ + (only one if the polygon is convex) + - when doing this for the whole stack of slabs, we get a set of rectangles: + To compute these rectangles, for each polygon, we compute the intersection + between : + - the line defined by the intersection of the polygon plane and the cutting + plane + - the polygon itself + This yields 0 or 2*K points along the line C. These are turned into K + rectangles by taking two consecutive points along the line and extruding + this segment by sliceThickness/2 in the orientation of the polygon normal, + in both directions. + + Then, once this list of rectangles is computed, we need to group the + connected rectangles together. Connected, here, means sharing at least part + of an edge --> union/find data structures and algorithm. + */ + bool Project(std::vector< std::pair<Point2D, Point2D> >& polygons, const CoordinateSystem3D& plane) const; + + std::string interpretation_; + std::string name_; + uint8_t red_; + uint8_t green_; + uint8_t blue_; + + /** Internal */ + const std::vector<DicomStructurePolygon2>& GetPolygons() const + { + return polygons_; + } + + /** Internal */ + double GetSliceThickness() const + { + return sliceThickness_; + } + + private: + enum State + { + Building, + NormalComputed, + Valid, // When normal components AND slice thickness are computed + Invalid + }; + + void ComputeNormal(); + void ComputeSliceThickness(); + + std::vector<DicomStructurePolygon2> polygons_; + Vector3D normal_; + double sliceThickness_; + + /* + After creation (and while polygons are added), state is Building. + After ComputeDependentProperties() is called, state can either be + Valid or Invalid. In any case, the object becomes immutable. + */ + State state_; + }; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/DicomStructurePolygon2.cpp Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,306 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#include "DicomStructurePolygon2.h" + +#include "../Toolbox/LinearAlgebra.h" + +namespace OrthancStone +{ + void DicomStructurePolygon2::ComputeDependentProperties() + { + ORTHANC_ASSERT(state_ == Building); + + for (size_t j = 0; j < points_.size(); ++j) + { + // TODO: move to AddPoint! + const Point3D& p = points_[j]; + if (p[0] < minX_) + minX_ = p[0]; + if (p[0] > maxX_) + maxX_ = p[0]; + + if (p[1] < minY_) + minY_ = p[1]; + if (p[1] > maxY_) + maxY_ = p[1]; + + if (p[2] < minZ_) + minZ_ = p[2]; + if (p[2] > maxZ_) + maxZ_ = p[2]; + } + + if (LinearAlgebra::IsNear(minX_, maxX_)) + { + LinearAlgebra::AssignVector(normal_, 1, 0, 0); + //ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX, maxX)); + ORTHANC_ASSERT(!LinearAlgebra::IsNear(minY_, maxY_)); + ORTHANC_ASSERT(!LinearAlgebra::IsNear(minZ_, maxZ_)); + } + else if (LinearAlgebra::IsNear(minY_, maxY_)) + { + LinearAlgebra::AssignVector(normal_, 0, 1, 0); + ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX_, maxX_)); + ORTHANC_ASSERT(!LinearAlgebra::IsNear(minZ_, maxZ_)); + } + else if (LinearAlgebra::IsNear(minZ_, maxZ_)) + { + LinearAlgebra::AssignVector(normal_, 0, 0, 1); + ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX_, maxX_)); + ORTHANC_ASSERT(!LinearAlgebra::IsNear(minY_, maxY_)); + } + else + { + LOG(ERROR) << "The contour is not coplanar and not parallel to any axis."; + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + state_ = Valid; + } + + + void DicomStructurePolygon2::ProjectOnConstantPlane( + std::vector<Point2D>& intersections, const CoordinateSystem3D& plane) const + { + // the plane can either have constant X, or constant Y. + // - for constant Z planes, use the ProjectOnParallelPlane method + // - other type of planes are not supported + + // V is the coordinate that is constant in the plane + double planeV = 0.0; + + // if true, then "u" in the code is "x" and "v" is "y". + // (v is constant in the plane) + bool uvxy = false; + // if true, then "u" in the code is "y" and "v" is "x" + // (v is constant in the plane) + bool uvyx = false; + + size_t uindex = static_cast<size_t>(-1); + size_t vindex = static_cast<size_t>(-1); + + ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[2], 0.0)); + + if (LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0)) + { + // normal is 1,0,0 (or -1,0,0). + // plane is constant X + uindex = 1; + vindex = 0; + + uvxy = false; + uvyx = true; + planeV = plane.GetOrigin()[0]; + if (planeV < minX_) + return; + if (planeV > maxX_) + return; + } + else if (LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0)) + { + // normal is 0,1,0 (or 0,-1,0). + // plane is constant Y + uindex = 0; + vindex = 1; + + uvxy = true; + uvyx = false; + planeV = plane.GetOrigin()[1]; + if (planeV < minY_) + return; + if (planeV > maxY_) + return; + } + else + { + // if the following assertion(s) fail(s), it means the plane is NOT a constant-X or constant-Y plane + LOG(ERROR) << "Plane normal must be (a,0,0) or (0,a,0), with a == -1 or a == 1"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + + size_t pointCount = GetPointCount(); + if (pointCount >= 3) + { + // this vector will contain the coordinates of the intersection points + // between the plane and the polygon. + // these are expressed in the U coordinate, that is either X or Y, + // depending upon the plane orientation + std::vector<double> uIntersections; + + // we loop on the segments of the polygon (TODO: optimize) + // and we compute the intersection between each segment and the cut + // cutting plane (slice) has a constant X + + for (size_t iPoint = 0; iPoint < pointCount; ++iPoint) + { + double u1 = points_[iPoint][uindex]; + double v1 = points_[iPoint][vindex]; + + double u2 = 0; + double v2 = 0; + + if (iPoint < pointCount - 1) + { + u2 = points_[iPoint + 1][uindex]; + v2 = points_[iPoint + 1][vindex]; + } + else + { + u2 = points_[0][uindex]; + v2 = points_[0][vindex]; + } + + // Check if the segment intersects the plane + if ((std::min(v1, v2) <= planeV) && (std::max(v1, v2) >= planeV)) + { + // special case: the segment is parallel to the plane but close to it + if (LinearAlgebra::IsNear(v1, v2)) + { + // in that case, we choose to label both points as an intersection + double x, y; + plane.ProjectPoint(x, y, points_[iPoint]); + intersections.push_back(Point2D(x, y)); + + plane.ProjectPoint(x, y, points_[iPoint + 1]); + intersections.push_back(Point2D(x, y)); + } + else + { + // we are looking for u so that (u,planeV) belongs to the segment + // let's define alpha = (u-u2)/(u1-u2) --> u = alpha*(u1-u2) + u2 + // alpha = (v2-planeV)/(v2-v1) + // because the following two triangles are similar + // [ (planeY,x) , (y2,x2), (planeY,x2) ] or + // [ (planeX,y) , (x2,y2), (planeX,y2) ] + // and + // [ (y1 ,x1) , (y2,x2), (y1 ,x2) ] or + // [ (x1 ,y1) , (x2,y2), (x1 ,y2) ] + + /* + void CoordinateSystem3D::ProjectPoint(double& offsetX, + double& offsetY, + const Vector& point) const + */ + double alpha = (v2 - planeV) / (v2 - v1); + + // get rid of numerical oddities + if (alpha < 0.0) + alpha = 0.0; + if (alpha > 1.0) + alpha = 1.0; + double u = alpha * (u1 - u2) + u2; + + // here is the intersection in world coordinates + Vector intersection; + if(uvxy) + LinearAlgebra::AssignVector(intersection, u, planeV, minZ_); + else + LinearAlgebra::AssignVector(intersection, planeV, u, minZ_); + + // and we convert it to plane coordinates + { + double xi, yi; + plane.ProjectPoint(xi, yi, intersection); + + // we consider that the x axis is always parallel to the polygons + // TODO: is this hypothesis safe?????? + uIntersections.insert(std::lower_bound(uIntersections.begin(), uIntersections.end(), xi), xi); + } + } + } + } // end of for (size_t iPoint = 0; iPoint < pointCount; ++iPoint) + + // now we convert the intersections to plane points + // we consider that the x axis is always parallel to the polygons + // TODO: same hypothesis as above: plane is perpendicular to polygons, + // plane is parallel to the XZ (constant Y) or YZ (constant X) 3D planes + for (size_t i = 0; i < uIntersections.size(); ++i) + { + double x = uIntersections[i]; + intersections.push_back(Point2D(x, minZ_)); + } + } // end of if (pointCount >= 3) + else + { + LOG(ERROR) << "This polygon has " << pointCount << " vertices, which is less than 3 --> skipping"; + } + } + +void OrthancStone::DicomStructurePolygon2::ProjectOnParallelPlane( + std::vector< std::pair<Point2D, Point2D> >& segments, + const CoordinateSystem3D& plane) const +{ + if (points_.size() < 3) + return; + + // the plane is horizontal + ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0)); + ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0)); + + size_t pointCount = GetPointCount(); + + segments.clear(); + segments.reserve(points_.size()); + // since the returned values need to be expressed in the supplied coordinate + // system, we need to subtract origin_ from the returned points + + double planeOriginX = plane.GetOrigin()[0]; + double planeOriginY = plane.GetOrigin()[1]; + + // precondition: points_.size() >= 3 + for (size_t j = 0; j < points_.size()-1; ++j) + { + // segment between point j and j+1 + + const Point3D& point0 = GetPoint(j); + // subtract plane origin x and y + Point2D p0(point0[0] - planeOriginX, point0[1] - planeOriginY); + + const Point3D& point1 = GetPoint(j+1); + // subtract plane origin x and y + Point2D p1(point1[0] - planeOriginX, point1[1] - planeOriginY); + + segments.push_back(std::pair<Point2D, Point2D>(p0,p1)); + } + + + // final segment + + const Point3D& point0 = GetPoint(points_.size() - 1); + // subtract plane origin x and y + Point2D p0(point0[0] - planeOriginX, point0[1] - planeOriginY); + + const Point3D& point1 = GetPoint(0); + // subtract plane origin x and y + Point2D p1(point1[0] - planeOriginX, point1[1] - planeOriginY); + + segments.push_back(std::pair<Point2D, Point2D>(p0, p1)); +} + +double OrthancStone::DicomStructurePolygon2::GetZ() const +{ + ORTHANC_ASSERT(LinearAlgebra::IsNear(normal_[0], 0.0)); + ORTHANC_ASSERT(LinearAlgebra::IsNear(normal_[1], 0.0)); + ORTHANC_ASSERT(LinearAlgebra::IsNear(minZ_, maxZ_)); + return minZ_; +} + + +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/DicomStructurePolygon2.h Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,155 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#pragma once + +#include "CoordinateSystem3D.h" +#include "DicomStructureSetUtils.h" +#include "Extent2D.h" + +#include "../Scene2D/Color.h" +#include "../StoneException.h" + +#include <Plugins/Samples/Common/FullOrthancDataset.h> + +#include <list> +#include <string> + +namespace OrthancStone +{ + + /** + Only polygons that are planar and parallel to either the X,Y or Z plane + ("X plane" == plane where X is equal to a constant for each point) are + supported. + */ + class DicomStructurePolygon2 + { + public: + enum Type + { + ClosedPlanar, + Unsupported + }; + + DicomStructurePolygon2(std::string referencedSopInstanceUid, const std::string& type) + : referencedSopInstanceUid_(referencedSopInstanceUid) + , state_(Building) + , minX_(std::numeric_limits<double>::max()) + , maxX_(-std::numeric_limits<double>::max()) + , minY_(std::numeric_limits<double>::max()) + , maxY_(-std::numeric_limits<double>::max()) + , minZ_(std::numeric_limits<double>::max()) + , maxZ_(-std::numeric_limits<double>::max()) + , type_(TypeFromString(type)) + { + ORTHANC_ASSERT(type_ == ClosedPlanar); + } + + void ComputeDependentProperties(); + + size_t GetPointCount() const + { + ORTHANC_ASSERT(state_ == Valid); + return points_.size(); + } + + const Point3D& GetPoint(size_t i) const + { + ORTHANC_ASSERT(state_ == Valid); + return points_.at(i); + } + + void AddPoint(const Point3D& v) + { + ORTHANC_ASSERT(state_ == Building); + points_.push_back(v); + } + + void Reserve(size_t n) + { + ORTHANC_ASSERT(state_ == Building); + points_.reserve(n); + } + + /** + This method takes a plane+coord system that is parallel to the polygon + and adds to polygons a new vector with the ordered set of points projected + on the plane, in the plane coordinate system. + */ + void ProjectOnParallelPlane( + std::vector< std::pair<Point2D,Point2D> >& segments, + const CoordinateSystem3D& plane) const; + + /** + Returns the coordinates of the intersection of the polygon and a plane + that is perpendicular to the polygons (plane has either constant X or + constant Y) + */ + void ProjectOnConstantPlane( + std::vector<Point2D>& intersections, + const CoordinateSystem3D& plane) const; + + /** + This method assumes polygon has a normal equal to 0,0,-1 and 0,0,1 (thus, + the polygon is parallel to the XY plane) and returns the Z coordinate of + all the polygon points + */ + double GetZ() const; + + /** + The normal sign is left undefined for now + */ + Vector3D GetNormal() const + { + return normal_; + } + + /** + This method will compute the intersection between a polygon and + a plane where either X, Y or Z is constant. + The plane is given with an origin and a normal. If the normal is + not parallel to an axis, an error is raised. + */ + void ComputeIntersectionWithPlane(const CoordinateSystem3D& plane); + + private: + static Type TypeFromString(const std::string& s) + { + if (s == "CLOSED_PLANAR") + return ClosedPlanar; + else + return Unsupported; + } + enum State + { + Building, + Valid + }; + std::string referencedSopInstanceUid_; + CoordinateSystem3D geometry_; + std::vector<Point3D> points_; + Vector3D normal_; // sign is irrelevant for now + State state_; + double minX_, maxX_, minY_, maxY_, minZ_, maxZ_; + Type type_; + }; +} +
--- a/Framework/Toolbox/DicomStructureSet2.cpp Sat Sep 14 17:27:41 2019 +0200 +++ b/Framework/Toolbox/DicomStructureSet2.cpp Fri Sep 20 11:58:00 2019 +0200 @@ -18,73 +18,19 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. **/ - #include "DicomStructureSet2.h" -#include "../Toolbox/GeometryToolbox.h" +#include "../Toolbox/LinearAlgebra.h" +#include "../StoneException.h" #include <Core/Logging.h> #include <Core/OrthancException.h> #include <Core/Toolbox.h> +#include <Core/DicomFormat/DicomTag.h> + #include <Plugins/Samples/Common/FullOrthancDataset.h> #include <Plugins/Samples/Common/DicomDatasetReader.h> -#include <limits> -#include <stdio.h> -#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> - -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) -{ - for (size_t i = 0; i < input.size(); i++) - { - boost::geometry::correct(input[i]); - } - - if (input.size() == 0) - { - output.clear(); - } - else if (input.size() == 1) - { - output.resize(1); - output[0] = input[0]; - } - else - { - boost::geometry::union_(input[0], input[1], output); - - for (size_t i = 0; i < input.size(); i++) - { - BoostMultiPolygon tmp; - boost::geometry::union_(output, input[i], tmp); - output = tmp; - } - } -} - - -static BoostPolygon CreateRectangle(float x1, float y1, - float x2, float y2) -{ - BoostPolygon r; - boost::geometry::append(r, BoostPoint(x1, y1)); - boost::geometry::append(r, BoostPoint(x1, y2)); - boost::geometry::append(r, BoostPoint(x2, y2)); - boost::geometry::append(r, BoostPoint(x2, y1)); - return r; -} - - - namespace OrthancStone { static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); @@ -100,8 +46,7 @@ static const OrthancPlugins::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080); static const OrthancPlugins::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020); - - static uint8_t ConvertColor(double v) + static inline uint8_t ConvertAndClipToByte(double v) { if (v < 0) { @@ -117,325 +62,98 @@ } } - - static bool ParseVector(Vector& target, - const OrthancPlugins::IDicomDataset& dataset, - const OrthancPlugins::DicomPath& tag) + static bool ReadDicomToVector(Vector& target, + const OrthancPlugins::IDicomDataset& dataset, + const OrthancPlugins::DicomPath& tag) { std::string value; return (dataset.GetStringValue(value, tag) && - LinearAlgebra::ParseVector(target, value)); + LinearAlgebra::ParseVector(target, value)); } - - void DicomStructureSet2::Polygon::CheckPointIsOnSlice(const Vector& v) const + + std::ostream& operator<<(std::ostream& s, const OrthancPlugins::DicomPath& dicomPath) { - if (hasSlice_) - { - double magnitude = - GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); - if(!LinearAlgebra::IsNear( - magnitude, - projectionAlongNormal_, - sliceThickness_ / 2.0 /* in mm */ )) - { - LOG(ERROR) << "This RT-STRUCT contains a point that is off the " - << "slice of its instance | " - << "magnitude = " << magnitude << " | " - << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " - << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); - - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - } - } - - bool DicomStructureSet2::Polygon::IsPointOnSliceIfAny(const Vector& v) const - { - if (hasSlice_) + std::stringstream tmp; + for (size_t i = 0; i < dicomPath.GetPrefixLength(); ++i) { - double magnitude = - GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); - bool onSlice = LinearAlgebra::IsNear( - magnitude, - projectionAlongNormal_, - sliceThickness_ / 2.0 /* in mm */); - if (!onSlice) - { - LOG(WARNING) << "This RT-STRUCT contains a point that is off the " - << "slice of its instance | " - << "magnitude = " << magnitude << " | " - << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " - << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); - } - return onSlice; + OrthancPlugins::DicomTag tag = dicomPath.GetPrefixTag(i); + + // We use this other object to be able to use GetMainTagsName + // and Format + Orthanc::DicomTag tag2(tag.GetGroup(), tag.GetElement()); + size_t index = dicomPath.GetPrefixIndex(i); + tmp << tag2.GetMainTagsName() << " (" << tag2.Format() << ") [" << index << "] / "; } - else - { - return false; - } - } - - void DicomStructureSet2::Polygon::AddPoint(const Vector& v) - { -#if 1 - // BGO 2019-09-03 - if (IsPointOnSliceIfAny(v)) - { - points_.push_back(v); - } -#else - CheckPoint(v); - points_.push_back(v); -#endif + const OrthancPlugins::DicomTag& tag = dicomPath.GetFinalTag(); + Orthanc::DicomTag tag2(tag.GetGroup(), tag.GetElement()); + tmp << tag2.GetMainTagsName() << " (" << tag2.Format() << ")"; + s << tmp.str(); + return s; } - bool DicomStructureSet2::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices) + void DicomStructureSet2::SetContents(const OrthancPlugins::FullOrthancDataset& tags) { - if (hasSlice_) - { - return true; - } - else + FillStructuresFromDataset(tags); + ComputeDependentProperties(); + } + + void DicomStructureSet2::ComputeDependentProperties() + { + for (size_t i = 0; i < structures_.size(); ++i) { - ReferencedSlices::const_iterator it = slices.find(sopInstanceUid_); - - if (it == slices.end()) - { - return false; - } - else - { - const CoordinateSystem3D& geometry = it->second.geometry_; - - hasSlice_ = true; - geometry_ = geometry; - projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); - sliceThickness_ = it->second.thickness_; - - extent_.Reset(); - - for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) - { - if (IsPointOnSliceIfAny(*it)) - { - double x, y; - geometry.ProjectPoint(x, y, *it); - extent_.AddPoint(x, y); - } - } - return true; - } + structures_[i].ComputeDependentProperties(); } } - bool DicomStructureSet2::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const - { - bool isOpposite = false; - - if (points_.empty() || - !hasSlice_ || - !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) - { - return false; - } - - double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); - - return (LinearAlgebra::IsNear(d, projectionAlongNormal_, - sliceThickness_ / 2.0)); - } - - bool DicomStructureSet2::Polygon::Project(double& x1, - double& y1, - double& x2, - double& y2, - const CoordinateSystem3D& slice) const - { - // TODO: optimize this method using a sweep-line algorithm for polygons - - if (!hasSlice_ || - points_.size() <= 1) - { - return false; - } - - double x, y; - geometry_.ProjectPoint(x, y, slice.GetOrigin()); - - bool isOpposite; - if (GeometryToolbox::IsParallelOrOpposite - (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) - { - 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<double>::infinity(); - double xmax = -std::numeric_limits<double>::infinity(); - - double prevX, prevY; - geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); - - for (size_t i = 0; i < points_.size(); i++) - { - // 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(xmax, y) - - sliceThickness_ / 2.0 * geometry_.GetNormal()); - - slice.ProjectPoint(x1, y1, p1); - slice.ProjectPoint(x2, y2, p2); - return true; - } - } - else if (GeometryToolbox::IsParallelOrOpposite - (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) - { - if (x < extent_.GetX1() || - x > extent_.GetX2()) - { - return false; - } - - bool isFirst = true; - double ymin = std::numeric_limits<double>::infinity(); - double ymax = -std::numeric_limits<double>::infinity(); - - double prevX, prevY; - geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); - - for (size_t i = 0; i < points_.size(); i++) - { - // 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, 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 - { - // Should not happen - return false; - } - } - - - const DicomStructureSet2::Structure& DicomStructureSet2::GetStructure(size_t index) const - { - if (index >= structures_.size()) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); - } - - return structures_[index]; - } - - - DicomStructureSet2::Structure& DicomStructureSet2::GetStructure(size_t index) - { - if (index >= structures_.size()) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); - } - - return structures_[index]; - } - - DicomStructureSet2::DicomStructureSet2(const OrthancPlugins::FullOrthancDataset& tags) + void DicomStructureSet2::FillStructuresFromDataset(const OrthancPlugins::FullOrthancDataset& tags) { OrthancPlugins::DicomDatasetReader reader(tags); - - size_t count, tmp; + + // a few sanity checks + size_t count = 0, tmp = 0; + + // DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE (0x3006, 0x0080); + // DICOM_TAG_ROI_CONTOUR_SEQUENCE (0x3006, 0x0039); + // DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE (0x3006, 0x0020); if (!tags.GetSequenceSize(count, DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE) || - !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) || - tmp != count || - !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) || - tmp != count) + !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) || + tmp != count || + !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) || + tmp != count) { throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); } + // let's now parse the structures stored in the dicom file + // DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE (0x3006, 0x0080) + // DICOM_TAG_RT_ROI_INTERPRETED_TYPE (0x3006, 0x00a4) + // DICOM_TAG_ROI_DISPLAY_COLOR (0x3006, 0x002a) + // DICOM_TAG_ROI_NAME (0x3006, 0x0026) structures_.resize(count); for (size_t i = 0; i < count; i++) { + // (0x3006, 0x0080)[i]/(0x3006, 0x00a4) structures_[i].interpretation_ = reader.GetStringValue - (OrthancPlugins::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, - DICOM_TAG_RT_ROI_INTERPRETED_TYPE), - "No interpretation"); + (OrthancPlugins::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, + DICOM_TAG_RT_ROI_INTERPRETED_TYPE), + "No interpretation"); + // (0x3006, 0x0020)[i]/(0x3006, 0x0026) structures_[i].name_ = reader.GetStringValue - (OrthancPlugins::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, - DICOM_TAG_ROI_NAME), - "No name"); + (OrthancPlugins::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, + DICOM_TAG_ROI_NAME), + "No name"); Vector color; - if (ParseVector(color, tags, OrthancPlugins::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_ROI_DISPLAY_COLOR)) && - color.size() == 3) + // (0x3006, 0x0039)[i]/(0x3006, 0x002a) + if (ReadDicomToVector(color, tags, OrthancPlugins::DicomPath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, DICOM_TAG_ROI_DISPLAY_COLOR)) + && color.size() == 3) { - structures_[i].red_ = ConvertColor(color[0]); - structures_[i].green_ = ConvertColor(color[1]); - structures_[i].blue_ = ConvertColor(color[2]); + structures_[i].red_ = ConvertAndClipToByte(color[0]); + structures_[i].green_ = ConvertAndClipToByte(color[1]); + structures_[i].blue_ = ConvertAndClipToByte(color[2]); } else { @@ -445,91 +163,103 @@ } size_t countSlices; - if (!tags.GetSequenceSize(countSlices, OrthancPlugins::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_CONTOUR_SEQUENCE))) + // DICOM_TAG_ROI_CONTOUR_SEQUENCE (0x3006, 0x0039); + // DICOM_TAG_CONTOUR_SEQUENCE (0x3006, 0x0040); + if (!tags.GetSequenceSize(countSlices, OrthancPlugins::DicomPath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, DICOM_TAG_CONTOUR_SEQUENCE))) { + LOG(WARNING) << "DicomStructureSet2::SetContents | structure \"" << structures_[i].name_ << "\" has no slices!"; countSlices = 0; } - LOG(INFO) << "New RT structure: \"" << structures_[i].name_ - << "\" with interpretation \"" << structures_[i].interpretation_ - << "\" containing " << countSlices << " slices (color: " - << static_cast<int>(structures_[i].red_) << "," - << static_cast<int>(structures_[i].green_) << "," - << static_cast<int>(structures_[i].blue_) << ")"; + LOG(INFO) << "New RT structure: \"" << structures_[i].name_ + << "\" with interpretation \"" << structures_[i].interpretation_ + << "\" containing " << countSlices << " slices (color: " + << static_cast<int>(structures_[i].red_) << "," + << static_cast<int>(structures_[i].green_) << "," + << static_cast<int>(structures_[i].blue_) << ")"; // These temporary variables avoid allocating many vectors in the loop below - OrthancPlugins::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_CONTOUR_SEQUENCE, 0, - DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); + + // (0x3006, 0x0039)[i]/(0x3006, 0x0040)[0]/(0x3006, 0x0046) + OrthancPlugins::DicomPath countPointsPath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, + DICOM_TAG_CONTOUR_SEQUENCE, 0, + DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); - OrthancPlugins::DicomPath geometricTypePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_CONTOUR_SEQUENCE, 0, - DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); - - OrthancPlugins::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_CONTOUR_SEQUENCE, 0, - DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); + OrthancPlugins::DicomPath geometricTypePath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, + DICOM_TAG_CONTOUR_SEQUENCE, 0, + DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); + + OrthancPlugins::DicomPath imageSequencePath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, + DICOM_TAG_CONTOUR_SEQUENCE, 0, + DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155) - OrthancPlugins::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_CONTOUR_SEQUENCE, 0, - DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, - DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); + OrthancPlugins::DicomPath referencedInstancePath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, + DICOM_TAG_CONTOUR_SEQUENCE, 0, + DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, + DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); - OrthancPlugins::DicomPath contourDataPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, - DICOM_TAG_CONTOUR_SEQUENCE, 0, - DICOM_TAG_CONTOUR_DATA); + OrthancPlugins::DicomPath contourDataPath( + DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, + DICOM_TAG_CONTOUR_SEQUENCE, 0, + DICOM_TAG_CONTOUR_DATA); for (size_t j = 0; j < countSlices; j++) { - unsigned int countPoints; + unsigned int countPoints = 0; countPointsPath.SetPrefixIndex(1, j); if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath)) { + LOG(ERROR) << "Dicom path " << countPointsPath << " is not valid (should contain an unsigned integer)"; throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); } - + //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices"; geometricTypePath.SetPrefixIndex(1, j); std::string type = reader.GetMandatoryStringValue(geometricTypePath); if (type != "CLOSED_PLANAR") { + // TODO: support points!! LOG(WARNING) << "Ignoring contour with geometry type: " << type; continue; } - size_t size; + size_t size = 0; imageSequencePath.SetPrefixIndex(1, j); if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1) { LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry."; - throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); } referencedInstancePath.SetPrefixIndex(1, j); std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); - contourDataPath.SetPrefixIndex(1, j); + contourDataPath.SetPrefixIndex(1, j); std::string slicesData = reader.GetMandatoryStringValue(contourDataPath); Vector points; if (!LinearAlgebra::ParseVector(points, slicesData) || - points.size() != 3 * countPoints) + points.size() != 3 * countPoints) { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); } // seen in real world - if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") + if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") { LOG(ERROR) << "WARNING. The following Dicom tag (Referenced SOP Instance UID) contains an empty value : // (3006,0039)[" << i << "] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)"; } - Polygon polygon(sopInstanceUid); + DicomStructurePolygon2 polygon(sopInstanceUid,type); polygon.Reserve(countPoints); for (size_t k = 0; k < countPoints; k++) @@ -540,297 +270,15 @@ v[2] = points[3 * k + 2]; polygon.AddPoint(v); } - - structures_[i].polygons_.push_back(polygon); - } - } - } - - - Vector DicomStructureSet2::GetStructureCenter(size_t index) const - { - const Structure& structure = GetStructure(index); - - Vector center; - LinearAlgebra::AssignVector(center, 0, 0, 0); - if (structure.polygons_.empty()) - { - return center; - } - - double n = static_cast<double>(structure.polygons_.size()); - - for (Polygons::const_iterator polygon = structure.polygons_.begin(); - polygon != structure.polygons_.end(); ++polygon) - { - if (!polygon->GetPoints().empty()) - { - center += polygon->GetPoints().front() / n; - } - } - - return center; - } - - - const std::string& DicomStructureSet2::GetStructureName(size_t index) const - { - return GetStructure(index).name_; - } - - - const std::string& DicomStructureSet2::GetStructureInterpretation(size_t index) const - { - return GetStructure(index).interpretation_; - } - - - Color DicomStructureSet2::GetStructureColor(size_t index) const - { - const Structure& s = GetStructure(index); - return Color(s.red_, s.green_, s.blue_); - } - - - void DicomStructureSet2::GetStructureColor(uint8_t& red, - uint8_t& green, - uint8_t& blue, - size_t index) const - { - const Structure& s = GetStructure(index); - red = s.red_; - green = s.green_; - blue = s.blue_; - } - - - void DicomStructureSet2::GetReferencedInstances(std::set<std::string>& instances) - { - for (Structures::const_iterator structure = structures_.begin(); - structure != structures_.end(); ++structure) - { - for (Polygons::const_iterator polygon = structure->polygons_.begin(); - polygon != structure->polygons_.end(); ++polygon) - { - instances.insert(polygon->GetSopInstanceUid()); - } - } - } - - - void DicomStructureSet2::AddReferencedSlice(const std::string& sopInstanceUid, - const std::string& seriesInstanceUid, - const CoordinateSystem3D& geometry, - double thickness) - { - if (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()) - { - // This geometry is already known - LOG(ERROR) << "DicomStructureSet2::AddReferencedSlice(): (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()). sopInstanceUid = " << sopInstanceUid; - - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); - } - else - { - if (thickness < 0) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); - } - - if (!referencedSlices_.empty()) - { - const ReferencedSlice& reference = referencedSlices_.begin()->second; - - if (reference.seriesInstanceUid_ != seriesInstanceUid) - { - LOG(ERROR) << "This RT-STRUCT refers to several different series"; - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - - if (!GeometryToolbox::IsParallel(reference.geometry_.GetNormal(), geometry.GetNormal())) - { - LOG(ERROR) << "The slices in this RT-STRUCT are not parallel"; - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - } - - referencedSlices_[sopInstanceUid] = ReferencedSlice(seriesInstanceUid, geometry, thickness); - - for (Structures::iterator structure = structures_.begin(); - structure != structures_.end(); ++structure) - { - for (Polygons::iterator polygon = structure->polygons_.begin(); - polygon != structure->polygons_.end(); ++polygon) - { - polygon->UpdateReferencedSlice(referencedSlices_); - } + structures_[i].AddPolygon(polygon); } } } - void DicomStructureSet2::AddReferencedSlice(const Orthanc::DicomMap& dataset) - { - CoordinateSystem3D slice(dataset); - - double thickness = 1; // 1 mm by default - - std::string s; - Vector v; - if (dataset.LookupStringValue(s, Orthanc::DICOM_TAG_SLICE_THICKNESS, false) && - LinearAlgebra::ParseVector(v, s) && - v.size() > 0) - { - thickness = v[0]; - } - - std::string instance, series; - if (dataset.LookupStringValue(instance, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false) && - dataset.LookupStringValue(series, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false)) - { - AddReferencedSlice(instance, series, slice, thickness); - } - else - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - } - - - void DicomStructureSet2::CheckReferencedSlices() + void DicomStructureSet2::Clear() { - for (Structures::iterator structure = structures_.begin(); - structure != structures_.end(); ++structure) - { - for (Polygons::iterator polygon = structure->polygons_.begin(); - polygon != structure->polygons_.end(); ++polygon) - { - if (!polygon->UpdateReferencedSlice(referencedSlices_)) - { - std::string sopInstanceUid = polygon->GetSopInstanceUid(); - if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") - { - LOG(ERROR) << "DicomStructureSet2::CheckReferencedSlices(): " - << " missing information about referenced instance " - << "(sopInstanceUid is empty!)"; - } - else - { - LOG(ERROR) << "DicomStructureSet2::CheckReferencedSlices(): " - << " missing information about referenced instance " - << "(sopInstanceUid = " << sopInstanceUid << ")"; - } - //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); - } - } - } - } - - - Vector DicomStructureSet2::GetNormal() const - { - if (referencedSlices_.empty()) - { - Vector v; - LinearAlgebra::AssignVector(v, 0, 0, 1); - return v; - } - else - { - return referencedSlices_.begin()->second.geometry_.GetNormal(); - } + structures_.clear(); } - - bool DicomStructureSet2::ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons, - const Structure& structure, - const CoordinateSystem3D& slice) const - { - polygons.clear(); - - Vector normal = GetNormal(); - - bool isOpposite; - if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) - { - // This is an axial projection - - for (Polygons::const_iterator polygon = structure.polygons_.begin(); - polygon != structure.polygons_.end(); ++polygon) - { - if (polygon->IsOnSlice(slice)) - { - polygons.push_back(std::vector<PolygonPoint2D>()); - - for (Points::const_iterator p = polygon->GetPoints().begin(); - p != polygon->GetPoints().end(); ++p) - { - double x, y; - slice.ProjectPoint(x, y, *p); - polygons.back().push_back(std::make_pair(x, y)); - } - } - } - - return true; - } - else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || - GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) - { -#if 1 - // Sagittal or coronal projection - std::vector<BoostPolygon> projected; - - for (Polygons::const_iterator polygon = structure.polygons_.begin(); - polygon != structure.polygons_.end(); ++polygon) - { - double x1, y1, x2, y2; - if (polygon->Project(x1, y1, x2, y2, slice)) - { - projected.push_back(CreateRectangle( - static_cast<float>(x1), - static_cast<float>(y1), - static_cast<float>(x2), - static_cast<float>(y2))); - } - } - - BoostMultiPolygon merged; - Union(merged, projected); - - polygons.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()); - for (size_t j = 0; j < outer.size(); j++) - { - 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<PolygonPoint2D> 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; - } - else - { - return false; - } - } }
--- a/Framework/Toolbox/DicomStructureSet2.h Sat Sep 14 17:27:41 2019 +0200 +++ b/Framework/Toolbox/DicomStructureSet2.h Fri Sep 20 11:58:00 2019 +0200 @@ -18,9 +18,9 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. **/ - #pragma once +#include "DicomStructure2.h" #include "CoordinateSystem3D.h" #include "Extent2D.h" #include "../Scene2D/Color.h" @@ -31,155 +31,34 @@ namespace OrthancStone { - class DicomStructureSet2 : public boost::noncopyable + class DicomStructureSet2 : public boost::noncopyable { public: - typedef std::pair<double, double> PolygonPoint2D; - - private: - struct ReferencedSlice - { - std::string seriesInstanceUid_; - CoordinateSystem3D geometry_; - double thickness_; - - ReferencedSlice() - { - } - - ReferencedSlice(const std::string& seriesInstanceUid, - const CoordinateSystem3D& geometry, - double thickness) : - seriesInstanceUid_(seriesInstanceUid), - geometry_(geometry), - thickness_(thickness) - { - } - }; - - typedef std::map<std::string, ReferencedSlice> ReferencedSlices; - - typedef std::vector<Vector> Points; - - class Polygon - { - private: - std::string sopInstanceUid_; - bool hasSlice_; - CoordinateSystem3D geometry_; - double projectionAlongNormal_; - double sliceThickness_; // In millimeters - Points points_; - Extent2D extent_; - - void CheckPointIsOnSlice(const Vector& v) const; - bool IsPointOnSliceIfAny(const Vector& v) const; - - public: - Polygon(const std::string& sopInstanceUid) : - sopInstanceUid_(sopInstanceUid), - hasSlice_(false) - { - } - - void Reserve(size_t n) - { - points_.reserve(n); - } - - void AddPoint(const Vector& v); - - bool UpdateReferencedSlice(const ReferencedSlices& slices); - - bool IsOnSlice(const CoordinateSystem3D& geometry) const; - - const std::string& GetSopInstanceUid() const - { - return sopInstanceUid_; - } - - const Points& GetPoints() const - { - return points_; - } - - double GetSliceThickness() const - { - return sliceThickness_; - } - - bool Project(double& x1, - double& y1, - double& x2, - double& y2, - const CoordinateSystem3D& slice) const; - }; - - typedef std::list<Polygon> Polygons; - - struct Structure - { - std::string name_; - std::string interpretation_; - Polygons polygons_; - uint8_t red_; - uint8_t green_; - uint8_t blue_; - }; - - typedef std::vector<Structure> Structures; - - Structures structures_; - ReferencedSlices referencedSlices_; - - const Structure& GetStructure(size_t index) const; - - Structure& GetStructure(size_t index); - - bool ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons, - const Structure& structure, - const CoordinateSystem3D& slice) const; - - public: - DicomStructureSet2(const OrthancPlugins::FullOrthancDataset& instance); + DicomStructureSet2(); + ~DicomStructureSet2(); + + void SetContents(const OrthancPlugins::FullOrthancDataset& tags); size_t GetStructuresCount() const { return structures_.size(); } - Vector GetStructureCenter(size_t index) const; - - const std::string& GetStructureName(size_t index) const; - - const std::string& GetStructureInterpretation(size_t index) const; - - Color GetStructureColor(size_t index) const; + void Clear(); - // TODO - remove - void GetStructureColor(uint8_t& red, - uint8_t& green, - uint8_t& blue, - size_t index) const; - - void GetReferencedInstances(std::set<std::string>& instances); + const DicomStructure2& GetStructure(size_t i) const + { + // at() is like []() but with range check + return structures_.at(i); + } - void AddReferencedSlice(const std::string& sopInstanceUid, - const std::string& seriesInstanceUid, - const CoordinateSystem3D& geometry, - double thickness); - - void AddReferencedSlice(const Orthanc::DicomMap& dataset); - - void CheckReferencedSlices(); + /** Internal use only */ + void FillStructuresFromDataset(const OrthancPlugins::FullOrthancDataset& tags); - Vector GetNormal() const; + /** Internal use only */ + void ComputeDependentProperties(); - bool ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons, - size_t index, - const CoordinateSystem3D& slice) const - { - return ProjectStructure(polygons, GetStructure(index), slice); - } + /** Internal use only */ + std::vector<DicomStructure2> structures_; }; }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/DicomStructureSetUtils.cpp Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,276 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#include "DicomStructureSetUtils.h" + +namespace OrthancStone +{ + +#if 0 + void DicomStructure2::PartitionRectangleList(std::vector< std::vector<size_t> > & sets, const std::vector<RtStructRectanglesInSlab> slabCuts) + { + // map position ( <slabIndex,rectIndex> )--> disjoint set index + std::map<std::pair<size_t, size_t>, size_t> posToIndex; + + // disjoint set index --> position + std::map<size_t, std::pair<size_t, size_t> > indexToPos; + + size_t nextIndex = 0; + for (size_t i = 0; i < slabCuts.size(); ++i) + { + for (size_t j = 0; j < slabCuts[i].size(); ++j) + { + std::pair<size_t, size_t> pos(i, j); + posToIndex<pos> = nextIndex; + indexToPos<nextIndex> = pos; + } + } + // nextIndex is now the total rectangle count + DisjointDataSet ds(nextIndex); + + // we loop on all slabs (except the last one) and we connect all rectangles + if (slabCuts.size() < 2) + { +#error write special case + } + else + { + for (size_t i = 0; i < slabCuts.size() - 1; ++i) + { + for (size_t j = 0; j < slabCuts[i].size(); ++j) + { + const RtStructRectangleInSlab& r1 = slabCuts[i][j]; + const size_t r1i = posToIndex(std::pair<size_t, size_t>(i, j)); + for (size_t k = 0; k < slabCuts[i + 1].size(); ++k) + { + const RtStructRectangleInSlab& r2 = slabCuts[i + 1][k]; + const size_t r2i = posToIndex(std::pair<size_t, size_t>(i, j)); + // rect.xmin <= rectBottom.xmax && rectBottom.xmin <= rect.xmax + if ((r1.xmin <= r2.xmax) && (r2.xmin <= r1.xmax)) + { +#error now go! + } + + } + } + } + } +#endif + + /* + + compute list of segments : + + numberOfRectsFromHereOn = 0 + possibleNext = {in_k,in_kplus1} + + for all boundaries: + - we create a vertical segment and we push it + - if boundary is a start, numberOfRectsFromHereOn += 1. + - if we switch from 0 to 1, we start a segment + - if we switch from 1 to 2, we end the current segment and we record it + - if boundary is an end, numberOfRectsFromHereOn -= 1. + - if we switch from 1 to 0, we end the current segment and we record it + - if we switch from 2 to 1, we start a segment + */ + + // static + void AddSlabBoundaries( + std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries, + const std::vector<RtStructRectanglesInSlab> & slabCuts, size_t iSlab) + { + if (iSlab < slabCuts.size()) + { + const RtStructRectanglesInSlab& slab = slabCuts[iSlab]; + for (size_t iRect = 0; iRect < slab.size(); ++iRect) + { + const RtStructRectangleInSlab& rect = slab[iRect]; + { + std::pair<double, RectangleBoundaryKind> boundary(rect.xmin, RectangleBoundaryKind_Start); + boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary); + } + { + std::pair<double, RectangleBoundaryKind> boundary(rect.xmax, RectangleBoundaryKind_End); + boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary); + } + } + } + } + + // static + void ProcessBoundaryList( + std::vector< std::pair<Point2D, Point2D> > & segments, + const std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries, + double y) + { + Point2D start; + Point2D end; + int curNumberOfSegments = 0; // we count the number of segments. we only draw if it is 1 (not 0 or 2) + for (size_t i = 0; i < boundaries.size(); ++i) + { + switch (boundaries[i].second) + { + case RectangleBoundaryKind_Start: + curNumberOfSegments += 1; + switch (curNumberOfSegments) + { + case 0: + assert(false); + break; + case 1: + // a new segment has begun! + start.x = boundaries[i].first; + start.y = y; + break; + case 2: + // an extra segment has begun : stop the current one (we don't draw overlaps) + end.x = boundaries[i].first; + end.y = y; + segments.push_back(std::pair<Point2D, Point2D>(start, end)); + break; + default: + //assert(false); // seen IRL ! + break; + } + break; + case RectangleBoundaryKind_End: + curNumberOfSegments -= 1; + switch (curNumberOfSegments) + { + case 0: + // a lone (thus active) segment has ended. + end.x = boundaries[i].first; + end.y = y; + segments.push_back(std::pair<Point2D, Point2D>(start, end)); + break; + case 1: + // an extra segment has ended : start a new one one + start.x = boundaries[i].first; + start.y = y; + break; + default: + // this should not happen! + //assert(false); + break; + } + break; + default: + assert(false); + break; + } + } + } + +#if 0 + void ConvertListOfSlabsToSegments( + std::vector< std::pair<Point2D, Point2D> >& segments, + const std::vector<RtStructRectanglesInSlab>& slabCuts, + const size_t totalRectCount) + { +#error to delete + } +#else + // See https://www.dropbox.com/s/bllco6q8aazxk44/2019-09-18-rtstruct-cut-algorithm-rect-merge.png + void ConvertListOfSlabsToSegments( + std::vector< std::pair<Point2D, Point2D> > & segments, + const std::vector<RtStructRectanglesInSlab> & slabCuts, + const size_t totalRectCount) + { + if (slabCuts.size() == 0) + return; + + if (totalRectCount > 0) + segments.reserve(4 * totalRectCount); // worst case, but common. + + /* + VERTICAL + */ + for (size_t iSlab = 0; iSlab < slabCuts.size(); ++iSlab) + { + for (size_t iRect = 0; iRect < slabCuts[iSlab].size(); ++iRect) + { + const RtStructRectangleInSlab& rect = slabCuts[iSlab][iRect]; + { + Point2D p1(rect.xmin, rect.ymin); + Point2D p2(rect.xmin, rect.ymax); + segments.push_back(std::pair<Point2D, Point2D>(p1, p2)); + } + { + Point2D p1(rect.xmax, rect.ymin); + Point2D p2(rect.xmax, rect.ymax); + segments.push_back(std::pair<Point2D, Point2D>(p1, p2)); + } + } + } + + /* + HORIZONTAL + */ + + // if we have N slabs, we have N+1 potential vertical positions for horizontal segments + // - one for top of slab 0 + // - N-1 for all positions between two slabs + // - one for bottom of slab N-1 + + // this adds all the horizontal segments for the tops of 3the rectangles + // in row 0 + if (slabCuts[0].size() > 0) + { + std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; + AddSlabBoundaries(boundaries, slabCuts, 0); + + ProcessBoundaryList(segments, boundaries, slabCuts[0][0].ymin); + } + + // this adds all the horizontal segments belonging to two slabs + for (size_t iSlab = 0; iSlab < slabCuts.size() - 1; ++iSlab) + { + std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; + AddSlabBoundaries(boundaries, slabCuts, iSlab); + AddSlabBoundaries(boundaries, slabCuts, iSlab + 1); + double curY = 0; + if (slabCuts[iSlab].size() > 0) + { + curY = slabCuts[iSlab][0].ymax; + ProcessBoundaryList(segments, boundaries, curY); + } + else if (slabCuts[iSlab + 1].size() > 0) + { + curY = slabCuts[iSlab + 1][0].ymin; + ProcessBoundaryList(segments, boundaries, curY); + } + else + { + // nothing to do!! : both slab lists are empty! + } + } + + // this adds all the horizontal segments for the BOTTOM of the rectangles + // on last row + if (slabCuts[slabCuts.size() - 1].size() > 0) + { + std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; + AddSlabBoundaries(boundaries, slabCuts, slabCuts.size() - 1); + + ProcessBoundaryList(segments, boundaries, slabCuts[slabCuts.size() - 1][0].ymax); + } + } +#endif + }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/DicomStructureSetUtils.h Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,84 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#pragma once + +#include <vector> +#include <utility> + +#include "../Toolbox/LinearAlgebra.h" + +namespace OrthancStone +{ +#if 0 + struct Point3D + { + Point3D(double x, double y, double z) : x(x), y(y), z(z) {} + Point3D() : x(0), y(0), z(0) {} + double x, y, z; + }; + + struct Vector3D + { + Vector3D(double x, double y, double z) : x(x), y(y), z(z) {} + Vector3D() : x(0), y(0), z(0) {} + double x, y, z; + }; +#else + typedef Vector Vector3D; + typedef Vector Point3D; +#endif + + struct Point2D + { + Point2D(double x, double y) : x(x), y(y) {} + Point2D() : x(0), y(0) {} + double x, y; + }; + + + /** Internal */ + struct RtStructRectangleInSlab + { + double xmin, xmax, ymin, ymax; + }; + typedef std::vector<RtStructRectangleInSlab> RtStructRectanglesInSlab; + + enum RectangleBoundaryKind + { + RectangleBoundaryKind_Start, + RectangleBoundaryKind_End + }; + +#if 0 + /** Internal */ + void PartitionRectangleList(std::vector< std::vector<size_t> > & sets, const std::vector<RtStructRectanglesInSlab>); +#endif + + /** Internal */ + void ConvertListOfSlabsToSegments(std::vector< std::pair<Point2D, Point2D> >& segments, const std::vector<RtStructRectanglesInSlab>& slabCuts, const size_t totalRectCount); + + /** Internal */ + void AddSlabBoundaries(std::vector<std::pair<double, RectangleBoundaryKind> >& boundaries, const std::vector<RtStructRectanglesInSlab>& slabCuts, size_t iSlab); + + /** Internal */ + void ProcessBoundaryList(std::vector< std::pair<Point2D, Point2D> >& segments, const std::vector<std::pair<double, RectangleBoundaryKind> >& boundaries, double y); + +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Volumes/DicomStructureSetSlicer2.cpp Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,109 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + +#include "DicomStructureSetSlicer2.h" + +#include "../Toolbox/GeometryToolbox.h" +#include "../Volumes/IVolumeSlicer.h" +#include "../Scene2D/PolylineSceneLayer.h" + +namespace OrthancStone +{ + DicomStructureSetSlicer2::DicomStructureSetSlicer2(boost::shared_ptr<DicomStructureSet2> structureSet) + : structureSet_(structureSet) + {} + + IVolumeSlicer::IExtractedSlice* DicomStructureSetSlicer2::ExtractSlice(const CoordinateSystem3D& cuttingPlane) + { + // revision is always the same, hence 0 + return new DicomStructureSetSlice2(structureSet_, 0, cuttingPlane); + } + + DicomStructureSetSlice2::DicomStructureSetSlice2( + boost::weak_ptr<DicomStructureSet2> structureSet, + uint64_t revision, + const CoordinateSystem3D& cuttingPlane) + : structureSet_(structureSet.lock()) + , isValid_(false) + { + bool opposite = false; + + if (structureSet_->GetStructuresCount() == 0) + { + isValid_ = false; + } + else + { + // some structures seen in real life have no polygons. We must be + // careful + bool found = false; + size_t curStructure = 0; + while (!found && curStructure < structureSet_->GetStructuresCount()) + { + if (structureSet_->GetStructure(curStructure).IsValid()) + { + found = true; + const Vector normal = structureSet_->GetStructure(0).GetNormal(); + isValid_ = ( + GeometryToolbox::IsParallelOrOpposite(opposite, normal, cuttingPlane.GetNormal()) || + GeometryToolbox::IsParallelOrOpposite(opposite, normal, cuttingPlane.GetAxisX()) || + GeometryToolbox::IsParallelOrOpposite(opposite, normal, cuttingPlane.GetAxisY())); + } + } + } + } + + ISceneLayer* DicomStructureSetSlice2::CreateSceneLayer( + const ILayerStyleConfigurator* configurator, + const CoordinateSystem3D& cuttingPlane) + { + assert(isValid_); + + std::auto_ptr<PolylineSceneLayer> layer(new PolylineSceneLayer); + layer->SetThickness(2); // thickness of the on-screen line + + for (size_t i = 0; i < structureSet_->GetStructuresCount(); i++) + { + const DicomStructure2& structure = structureSet_->GetStructure(i); + if (structure.IsValid()) + { + const Color& color = structure.GetColor(); + + std::vector< std::pair<Point2D, Point2D> > segments; + + if (structure.Project(segments, cuttingPlane)) + { + for (size_t j = 0; j < segments.size(); j++) + { + PolylineSceneLayer::Chain chain; + chain.resize(2); + + chain[0] = ScenePoint2D(segments[j].first.x, segments[j].first.y); + chain[1] = ScenePoint2D(segments[j].second.x, segments[j].second.y); + + layer->AddChain(chain, false /* NOT closed */, color); + } + } + } + } + return layer.release(); + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Volumes/DicomStructureSetSlicer2.h Fri Sep 20 11:58:00 2019 +0200 @@ -0,0 +1,70 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2019 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + + +#pragma once + +#include "../Toolbox/DicomStructureSet2.h" +#include "../Volumes/IVolumeSlicer.h" + +#include <boost/shared_ptr.hpp> +#include <boost/weak_ptr.hpp> + +namespace OrthancStone +{ + class DicomStructureSetSlice2 : public IVolumeSlicer::IExtractedSlice + { + public: + DicomStructureSetSlice2( + boost::weak_ptr<DicomStructureSet2> structureSet, + uint64_t revision, + const CoordinateSystem3D& cuttingPlane); + + virtual bool IsValid() ORTHANC_OVERRIDE + { + return isValid_; + } + + virtual uint64_t GetRevision() ORTHANC_OVERRIDE + { + return revision_; + } + + virtual ISceneLayer* CreateSceneLayer( + const ILayerStyleConfigurator* configurator, // possibly absent + const CoordinateSystem3D& cuttingPlane) ORTHANC_OVERRIDE; + + private: + boost::shared_ptr<DicomStructureSet2> structureSet_; + bool isValid_; + uint64_t revision_; + }; + + class DicomStructureSetSlicer2 : public IVolumeSlicer + { + public: + DicomStructureSetSlicer2(boost::shared_ptr<DicomStructureSet2> structureSet); + + /** IVolumeSlicer impl */ + virtual IExtractedSlice* ExtractSlice(const CoordinateSystem3D& cuttingPlane) ORTHANC_OVERRIDE; + private: + boost::weak_ptr<DicomStructureSet2> structureSet_; + }; +}