Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet.cpp @ 1178:3c7cdbf32e2a broker
DicomStructureSet::ProjectOntoLayer()
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 20 Nov 2019 14:12:11 +0100 |
parents | 1644de437a7b |
children | 9c8f557ea799 |
comparison
equal
deleted
inserted
replaced
1177:c3d219b6266b | 1178:3c7cdbf32e2a |
---|---|
29 #include <Core/Toolbox.h> | 29 #include <Core/Toolbox.h> |
30 #include <Plugins/Samples/Common/FullOrthancDataset.h> | 30 #include <Plugins/Samples/Common/FullOrthancDataset.h> |
31 #include <Plugins/Samples/Common/DicomDatasetReader.h> | 31 #include <Plugins/Samples/Common/DicomDatasetReader.h> |
32 | 32 |
33 #if defined(_MSC_VER) | 33 #if defined(_MSC_VER) |
34 #pragma warning(push) | 34 # pragma warning(push) |
35 #pragma warning(disable:4244) | 35 # pragma warning(disable:4244) |
36 #endif | 36 #endif |
37 | 37 |
38 #include <limits> | 38 #include <limits> |
39 #include <stdio.h> | 39 #include <stdio.h> |
40 #include <boost/geometry.hpp> | 40 #include <boost/geometry.hpp> |
41 #include <boost/geometry/geometries/point_xy.hpp> | 41 #include <boost/geometry/geometries/point_xy.hpp> |
42 #include <boost/geometry/geometries/polygon.hpp> | 42 #include <boost/geometry/geometries/polygon.hpp> |
43 #include <boost/geometry/multi/geometries/multi_polygon.hpp> | 43 #include <boost/geometry/multi/geometries/multi_polygon.hpp> |
44 | 44 |
45 #if defined(_MSC_VER) | 45 #if defined(_MSC_VER) |
46 #pragma warning(pop) | 46 # pragma warning(pop) |
47 #endif | 47 #endif |
48 | 48 |
49 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; | 49 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; |
50 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; | 50 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; |
51 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; | 51 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; |
79 output = tmp; | 79 output = tmp; |
80 } | 80 } |
81 } | 81 } |
82 } | 82 } |
83 | 83 |
84 #ifdef USE_BOOST_UNION_FOR_POLYGONS | 84 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
85 | 85 |
86 static BoostPolygon CreateRectangle(float x1, float y1, | 86 static BoostPolygon CreateRectangle(float x1, float y1, |
87 float x2, float y2) | 87 float x2, float y2) |
88 { | 88 { |
89 BoostPolygon r; | 89 BoostPolygon r; |
97 #else | 97 #else |
98 | 98 |
99 namespace OrthancStone | 99 namespace OrthancStone |
100 { | 100 { |
101 static RtStructRectangleInSlab CreateRectangle(float x1, float y1, | 101 static RtStructRectangleInSlab CreateRectangle(float x1, float y1, |
102 float x2, float y2) | 102 float x2, float y2) |
103 { | 103 { |
104 RtStructRectangleInSlab rect; | 104 RtStructRectangleInSlab rect; |
105 rect.xmin = std::min(x1, x2); | 105 rect.xmin = std::min(x1, x2); |
106 rect.xmax = std::max(x1, x2); | 106 rect.xmax = std::max(x1, x2); |
107 rect.ymin = std::min(y1, y2); | 107 rect.ymin = std::min(y1, y2); |
172 if (hasSlice_) | 172 if (hasSlice_) |
173 { | 173 { |
174 double magnitude = | 174 double magnitude = |
175 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); | 175 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); |
176 if(!LinearAlgebra::IsNear( | 176 if(!LinearAlgebra::IsNear( |
177 magnitude, | 177 magnitude, |
178 projectionAlongNormal_, | 178 projectionAlongNormal_, |
179 sliceThickness_ / 2.0 /* in mm */ )) | 179 sliceThickness_ / 2.0 /* in mm */ )) |
180 { | 180 { |
181 LOG(ERROR) << "This RT-STRUCT contains a point that is off the " | 181 LOG(ERROR) << "This RT-STRUCT contains a point that is off the " |
182 << "slice of its instance | " | 182 << "slice of its instance | " |
183 << "magnitude = " << magnitude << " | " | 183 << "magnitude = " << magnitude << " | " |
184 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " | 184 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " |
185 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); | 185 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); |
186 | 186 |
187 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 187 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
188 } | 188 } |
189 } | 189 } |
190 } | 190 } |
200 projectionAlongNormal_, | 200 projectionAlongNormal_, |
201 sliceThickness_ / 2.0 /* in mm */); | 201 sliceThickness_ / 2.0 /* in mm */); |
202 if (!onSlice) | 202 if (!onSlice) |
203 { | 203 { |
204 LOG(WARNING) << "This RT-STRUCT contains a point that is off the " | 204 LOG(WARNING) << "This RT-STRUCT contains a point that is off the " |
205 << "slice of its instance | " | 205 << "slice of its instance | " |
206 << "magnitude = " << magnitude << " | " | 206 << "magnitude = " << magnitude << " | " |
207 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " | 207 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " |
208 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); | 208 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); |
209 } | 209 } |
210 return onSlice; | 210 return onSlice; |
211 } | 211 } |
212 else | 212 else |
213 { | 213 { |
376 // plane is constant X => Sagittal view (remember that in the | 376 // plane is constant X => Sagittal view (remember that in the |
377 // sagittal projection, the normal must be swapped) | 377 // sagittal projection, the normal must be swapped) |
378 | 378 |
379 | 379 |
380 /* | 380 /* |
381 Please read the comments in the section above, by taking into account | 381 Please read the comments in the section above, by taking into account |
382 the fact that, in this case, the plane has a constant X, not Y (in | 382 the fact that, in this case, the plane has a constant X, not Y (in |
383 polygon geometry_ coordinates) | 383 polygon geometry_ coordinates) |
384 */ | 384 */ |
385 | 385 |
386 if (x < extent_.GetX1() || | 386 if (x < extent_.GetX1() || |
387 x > extent_.GetX2()) | 387 x > extent_.GetX2()) |
388 { | 388 { |
510 { | 510 { |
511 countSlices = 0; | 511 countSlices = 0; |
512 } | 512 } |
513 | 513 |
514 LOG(INFO) << "New RT structure: \"" << structures_[i].name_ | 514 LOG(INFO) << "New RT structure: \"" << structures_[i].name_ |
515 << "\" with interpretation \"" << structures_[i].interpretation_ | 515 << "\" with interpretation \"" << structures_[i].interpretation_ |
516 << "\" containing " << countSlices << " slices (color: " | 516 << "\" containing " << countSlices << " slices (color: " |
517 << static_cast<int>(structures_[i].red_) << "," | 517 << static_cast<int>(structures_[i].red_) << "," |
518 << static_cast<int>(structures_[i].green_) << "," | 518 << static_cast<int>(structures_[i].green_) << "," |
519 << static_cast<int>(structures_[i].blue_) << ")"; | 519 << static_cast<int>(structures_[i].blue_) << ")"; |
520 | 520 |
521 // These temporary variables avoid allocating many vectors in the loop below | 521 // These temporary variables avoid allocating many vectors in the loop below |
522 OrthancPlugins::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 522 OrthancPlugins::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, |
523 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 523 DICOM_TAG_CONTOUR_SEQUENCE, 0, |
524 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); | 524 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); |
769 { | 769 { |
770 std::string sopInstanceUid = polygon->GetSopInstanceUid(); | 770 std::string sopInstanceUid = polygon->GetSopInstanceUid(); |
771 if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") | 771 if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") |
772 { | 772 { |
773 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): " | 773 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): " |
774 << " missing information about referenced instance " | 774 << " missing information about referenced instance " |
775 << "(sopInstanceUid is empty!)"; | 775 << "(sopInstanceUid is empty!)"; |
776 } | 776 } |
777 else | 777 else |
778 { | 778 { |
779 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): " | 779 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): " |
780 << " missing information about referenced instance " | 780 << " missing information about referenced instance " |
781 << "(sopInstanceUid = " << sopInstanceUid << ")"; | 781 << "(sopInstanceUid = " << sopInstanceUid << ")"; |
782 } | 782 } |
783 //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | 783 //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); |
784 } | 784 } |
785 } | 785 } |
786 } | 786 } |
799 { | 799 { |
800 return referencedSlices_.begin()->second.geometry_.GetNormal(); | 800 return referencedSlices_.begin()->second.geometry_.GetNormal(); |
801 } | 801 } |
802 } | 802 } |
803 | 803 |
804 #ifdef USE_BOOST_UNION_FOR_POLYGONS | 804 bool DicomStructureSet::ProjectStructure( |
805 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<Point2D> >& polygons, | 805 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
806 const Structure& structure, | 806 std::vector< std::vector<Point2D> >& polygons, |
807 const CoordinateSystem3D& slice) const | |
808 #else | 807 #else |
809 bool DicomStructureSet::ProjectStructure(std::vector< std::pair<Point2D, Point2D> >& segments, | 808 std::vector< std::pair<Point2D, Point2D> >& segments, |
809 #endif | |
810 const Structure& structure, | 810 const Structure& structure, |
811 const CoordinateSystem3D& slice) const | 811 const CoordinateSystem3D& sourceSlice) const |
812 #endif | 812 { |
813 { | 813 const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); |
814 #ifdef USE_BOOST_UNION_FOR_POLYGONS | 814 |
815 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
815 polygons.clear(); | 816 polygons.clear(); |
816 #else | 817 #else |
817 segments.clear(); | 818 segments.clear(); |
818 #endif | 819 #endif |
819 | 820 |
827 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 828 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
828 polygon != structure.polygons_.end(); ++polygon) | 829 polygon != structure.polygons_.end(); ++polygon) |
829 { | 830 { |
830 if (polygon->IsOnSlice(slice)) | 831 if (polygon->IsOnSlice(slice)) |
831 { | 832 { |
832 #ifdef USE_BOOST_UNION_FOR_POLYGONS | 833 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
833 polygons.push_back(std::vector<Point2D>()); | 834 polygons.push_back(std::vector<Point2D>()); |
834 | 835 |
835 for (Points::const_iterator p = polygon->GetPoints().begin(); | 836 for (Points::const_iterator p = polygon->GetPoints().begin(); |
836 p != polygon->GetPoints().end(); ++p) | 837 p != polygon->GetPoints().end(); ++p) |
837 { | 838 { |
878 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 879 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) |
879 { | 880 { |
880 #if 1 | 881 #if 1 |
881 // Sagittal or coronal projection | 882 // Sagittal or coronal projection |
882 | 883 |
883 #ifdef USE_BOOST_UNION_FOR_POLYGONS | 884 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
884 std::vector<BoostPolygon> projected; | 885 std::vector<BoostPolygon> projected; |
885 | 886 |
886 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 887 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
887 polygon != structure.polygons_.end(); ++polygon) | 888 polygon != structure.polygons_.end(); ++polygon) |
888 { | 889 { |
910 double curZ = polygon->GetGeometryOrigin()[2]; | 911 double curZ = polygon->GetGeometryOrigin()[2]; |
911 | 912 |
912 // x1,y1 and x2,y2 are in "slice" coordinates (the cutting plane | 913 // x1,y1 and x2,y2 are in "slice" coordinates (the cutting plane |
913 // geometry) | 914 // geometry) |
914 projected.push_back(std::make_pair(CreateRectangle( | 915 projected.push_back(std::make_pair(CreateRectangle( |
915 static_cast<float>(x1), | 916 static_cast<float>(x1), |
916 static_cast<float>(y1), | 917 static_cast<float>(y1), |
917 static_cast<float>(x2), | 918 static_cast<float>(x2), |
918 static_cast<float>(y2)),curZ)); | 919 static_cast<float>(y2)),curZ)); |
919 } | 920 } |
920 } | 921 } |
921 #endif | 922 #endif |
922 | 923 |
923 #ifndef USE_BOOST_UNION_FOR_POLYGONS | 924 #if USE_BOOST_UNION_FOR_POLYGONS != 1 |
924 // projected contains a set of rectangles specified by two opposite | 925 // projected contains a set of rectangles specified by two opposite |
925 // corners (x1,y1,x2,y2) | 926 // corners (x1,y1,x2,y2) |
926 // we need to merge them | 927 // we need to merge them |
927 // each slab yields ONE polygon! | 928 // each slab yields ONE polygon! |
928 | 929 |
1008 else | 1009 else |
1009 { | 1010 { |
1010 return false; | 1011 return false; |
1011 } | 1012 } |
1012 } | 1013 } |
1014 | |
1015 | |
1016 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, | |
1017 const CoordinateSystem3D& plane, | |
1018 size_t structureIndex, | |
1019 const Color& color) const | |
1020 { | |
1021 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
1022 std::vector< std::vector<Point2D> > polygons; | |
1023 if (ProjectStructure(polygons, structureIndex, plane)) | |
1024 { | |
1025 for (size_t j = 0; j < polygons.size(); j++) | |
1026 { | |
1027 std::vector<ScenePoint2D> chain; | |
1028 chain.reserve(polygons[j].size()); | |
1029 | |
1030 for (size_t k = 0; k < polygons[j].size(); k++) | |
1031 { | |
1032 chain.push_back(ScenePoint2D(polygons[j][k].x, polygons[j][k].y)); | |
1033 } | |
1034 | |
1035 layer.AddChain(chain, true, color.GetRed(), color.GetGreen(), color.GetBlue()); | |
1036 } | |
1037 } | |
1038 | |
1039 #else | |
1040 std::vector< std::pair<Point2D, Point2D> > segments; | |
1041 | |
1042 if (ProjectStructure(segments, structureIndex, plane)) | |
1043 { | |
1044 for (size_t j = 0; j < segments.size(); j++) | |
1045 { | |
1046 std::vector<ScenePoint2D> chain(2); | |
1047 chain[0] = ScenePoint2D(segments[j].first.x, segments[j].first.y); | |
1048 chain[1] = ScenePoint2D(segments[j].second.x, segments[j].second.y); | |
1049 layer.AddChain(chain, false, color.GetRed(), color.GetGreen(), color.GetBlue()); | |
1050 } | |
1051 } | |
1052 #endif | |
1053 } | |
1013 } | 1054 } |