Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructureSet.cpp @ 1178:3c7cdbf32e2a broker

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 }