Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1909:782ba9eb6f22
improved variable names
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 01 Feb 2022 08:56:36 +0100 |
parents | affde38b84de |
children | f81cdf283859 |
comparison
equal
deleted
inserted
replaced
1908:affde38b84de | 1909:782ba9eb6f22 |
---|---|
241 return true; | 241 return true; |
242 } | 242 } |
243 } | 243 } |
244 } | 244 } |
245 | 245 |
246 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& slice, | 246 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& cuttingPlane, |
247 const Vector& estimatedNormal, | 247 const Vector& estimatedNormal, |
248 double estimatedSliceThickness) const | 248 double estimatedSliceThickness) const |
249 { | 249 { |
250 bool isOpposite = false; | 250 bool isOpposite = false; |
251 | 251 |
254 return false; | 254 return false; |
255 } | 255 } |
256 else if (hasSlice_) | 256 else if (hasSlice_) |
257 { | 257 { |
258 // Use the actual geometry of this specific slice | 258 // Use the actual geometry of this specific slice |
259 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) | 259 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, cuttingPlane.GetNormal(), geometry_.GetNormal())) |
260 { | 260 { |
261 return false; | 261 return false; |
262 } | 262 } |
263 else | 263 else |
264 { | 264 { |
265 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); | 265 double d = GeometryToolbox::ProjectAlongNormal(cuttingPlane.GetOrigin(), geometry_.GetNormal()); |
266 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, sliceThickness_ / 2.0)); | 266 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, sliceThickness_ / 2.0)); |
267 } | 267 } |
268 } | 268 } |
269 else | 269 else |
270 { | 270 { |
271 // Use the estimated geometry for the global RT-STRUCT volume | 271 // Use the estimated geometry for the global RT-STRUCT volume |
272 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), estimatedNormal)) | 272 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, cuttingPlane.GetNormal(), estimatedNormal)) |
273 { | 273 { |
274 return false; | 274 return false; |
275 } | 275 } |
276 else | 276 else |
277 { | 277 { |
278 double d1 = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), estimatedNormal); | 278 double d1 = GeometryToolbox::ProjectAlongNormal(cuttingPlane.GetOrigin(), estimatedNormal); |
279 double d2 = GeometryToolbox::ProjectAlongNormal(points_.front(), estimatedNormal); | 279 double d2 = GeometryToolbox::ProjectAlongNormal(points_.front(), estimatedNormal); |
280 return (LinearAlgebra::IsNear(d1, d2, estimatedSliceThickness / 2.0)); | 280 return (LinearAlgebra::IsNear(d1, d2, estimatedSliceThickness / 2.0)); |
281 } | 281 } |
282 } | 282 } |
283 } | 283 } |
284 | 284 |
285 bool DicomStructureSet::Polygon::Project(double& x1, | 285 bool DicomStructureSet::Polygon::Project(double& x1, |
286 double& y1, | 286 double& y1, |
287 double& x2, | 287 double& x2, |
288 double& y2, | 288 double& y2, |
289 const CoordinateSystem3D& slice, | 289 const CoordinateSystem3D& cuttingPlane, |
290 const Vector& estimatedNormal, | 290 const Vector& estimatedNormal, |
291 double estimatedSliceThickness) const | 291 double estimatedSliceThickness) const |
292 { | 292 { |
293 #if 0 | |
294 if (points_.size() <= 1) | |
295 { | |
296 return false; | |
297 } | |
298 | |
299 Vector normal = estimatedNormal; | |
300 double thickness = estimatedSliceThickness; | |
301 if (hasSlice_) | |
302 { | |
303 normal = geometry_.GetNormal(); | |
304 thickness = sliceThickness_; | |
305 } | |
306 | |
307 bool isOpposite; | |
308 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) && | |
309 !GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY())) | |
310 { | |
311 printf("UUUU\n"); | |
312 return false; | |
313 } | |
314 | |
315 return false; | |
316 | |
317 #else | |
293 if (!hasSlice_ || | 318 if (!hasSlice_ || |
294 points_.size() <= 1) | 319 points_.size() <= 1) |
295 { | 320 { |
296 return false; | 321 return false; |
297 } | 322 } |
298 | 323 |
299 double x, y; | 324 double x, y; |
300 geometry_.ProjectPoint2(x, y, slice.GetOrigin()); | 325 geometry_.ProjectPoint2(x, y, cuttingPlane.GetOrigin()); |
301 | 326 |
302 bool isOpposite; | 327 bool isOpposite; |
303 if (GeometryToolbox::IsParallelOrOpposite | 328 if (GeometryToolbox::IsParallelOrOpposite |
304 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | 329 (isOpposite, cuttingPlane.GetNormal(), geometry_.GetAxisY())) |
305 { | 330 { |
306 // plane is constant Y | 331 // plane is constant Y |
307 | 332 |
308 if (y < extent_.GetY1() || | 333 if (y < extent_.GetY1() || |
309 y > extent_.GetY2()) | 334 y > extent_.GetY2()) |
359 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 384 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
360 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - | 385 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - |
361 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 386 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
362 | 387 |
363 // then to the cutting plane geometry... | 388 // then to the cutting plane geometry... |
364 slice.ProjectPoint2(x1, y1, p1); | 389 cuttingPlane.ProjectPoint2(x1, y1, p1); |
365 slice.ProjectPoint2(x2, y2, p2); | 390 cuttingPlane.ProjectPoint2(x2, y2, p2); |
366 return true; | 391 return true; |
367 } | 392 } |
368 } | 393 } |
369 else if (GeometryToolbox::IsParallelOrOpposite | 394 else if (GeometryToolbox::IsParallelOrOpposite |
370 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | 395 (isOpposite, cuttingPlane.GetNormal(), geometry_.GetAxisX())) |
371 { | 396 { |
372 // plane is constant X => Sagittal view (remember that in the | 397 // plane is constant X => Sagittal view (remember that in the |
373 // sagittal projection, the normal must be swapped) | 398 // sagittal projection, the normal must be swapped) |
374 | 399 |
375 | 400 |
420 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + | 445 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + |
421 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 446 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
422 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - | 447 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - |
423 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 448 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
424 | 449 |
425 slice.ProjectPoint2(x1, y1, p1); | 450 cuttingPlane.ProjectPoint2(x1, y1, p1); |
426 slice.ProjectPoint2(x2, y2, p2); | 451 cuttingPlane.ProjectPoint2(x2, y2, p2); |
427 | 452 |
428 return true; | 453 return true; |
429 } | 454 } |
430 } | 455 } |
431 else | 456 else |
432 { | 457 { |
433 // Should not happen | 458 // Should not happen |
434 return false; | 459 return false; |
435 } | 460 } |
461 #endif | |
436 } | 462 } |
437 | 463 |
438 | 464 |
439 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const | 465 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const |
440 { | 466 { |
812 } | 838 } |
813 } | 839 } |
814 | 840 |
815 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<ScenePoint2D> >& chains, | 841 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<ScenePoint2D> >& chains, |
816 const Structure& structure, | 842 const Structure& structure, |
817 const CoordinateSystem3D& sourceSlice) const | 843 const CoordinateSystem3D& cuttingPlane) const |
818 { | 844 { |
819 const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); | 845 const CoordinateSystem3D cutting = CoordinateSystem3D::NormalizeCuttingPlane(cuttingPlane); |
820 | 846 |
821 chains.clear(); | 847 chains.clear(); |
822 | 848 |
823 Vector normal = GetNormal(); | 849 Vector normal = GetNormal(); |
824 | 850 |
825 bool isOpposite; | 851 bool isOpposite; |
826 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) | 852 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cutting.GetNormal())) |
827 { | 853 { |
828 // This is an axial projection | 854 // This is an axial projection |
829 | 855 |
830 chains.reserve(structure.polygons_.size()); | 856 chains.reserve(structure.polygons_.size()); |
831 | 857 |
832 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 858 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
833 polygon != structure.polygons_.end(); ++polygon) | 859 polygon != structure.polygons_.end(); ++polygon) |
834 { | 860 { |
835 const Points& points = polygon->GetPoints(); | 861 const Points& points = polygon->GetPoints(); |
836 | 862 |
837 if (polygon->IsOnSlice(slice, GetEstimatedNormal(), GetEstimatedSliceThickness()) && | 863 if (polygon->IsOnSlice(cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()) && |
838 !points.empty()) | 864 !points.empty()) |
839 { | 865 { |
840 chains.push_back(std::vector<ScenePoint2D>()); | 866 chains.push_back(std::vector<ScenePoint2D>()); |
841 chains.back().reserve(points.size() + 1); | 867 chains.back().reserve(points.size() + 1); |
842 | 868 |
843 for (Points::const_iterator p = points.begin(); | 869 for (Points::const_iterator p = points.begin(); |
844 p != points.end(); ++p) | 870 p != points.end(); ++p) |
845 { | 871 { |
846 double x, y; | 872 double x, y; |
847 slice.ProjectPoint2(x, y, *p); | 873 cutting.ProjectPoint2(x, y, *p); |
848 chains.back().push_back(ScenePoint2D(x, y)); | 874 chains.back().push_back(ScenePoint2D(x, y)); |
849 } | 875 } |
850 | 876 |
851 double x0, y0; | 877 double x0, y0; |
852 slice.ProjectPoint2(x0, y0, points.front()); | 878 cutting.ProjectPoint2(x0, y0, points.front()); |
853 chains.back().push_back(ScenePoint2D(x0, y0)); | 879 chains.back().push_back(ScenePoint2D(x0, y0)); |
854 } | 880 } |
855 } | 881 } |
856 | 882 |
857 return true; | 883 return true; |
858 } | 884 } |
859 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | 885 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cutting.GetAxisX()) || |
860 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 886 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cutting.GetAxisY())) |
861 { | 887 { |
862 // Sagittal or coronal projection | 888 // Sagittal or coronal projection |
863 | 889 |
864 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | 890 #if USE_BOOST_UNION_FOR_POLYGONS == 1 |
865 std::vector<BoostPolygon> projected; | 891 std::vector<BoostPolygon> projected; |
867 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 893 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
868 polygon != structure.polygons_.end(); ++polygon) | 894 polygon != structure.polygons_.end(); ++polygon) |
869 { | 895 { |
870 double x1, y1, x2, y2; | 896 double x1, y1, x2, y2; |
871 | 897 |
872 if (polygon->Project(x1, y1, x2, y2, slice, GetEstimatedNormal(), GetEstimatedSliceThickness())) | 898 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness())) |
873 { | 899 { |
874 projected.push_back(CreateRectangle(x1, y1, x2, y2)); | 900 projected.push_back(CreateRectangle(x1, y1, x2, y2)); |
875 } | 901 } |
876 } | 902 } |
877 | 903 |
897 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 923 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
898 polygon != structure.polygons_.end(); ++polygon) | 924 polygon != structure.polygons_.end(); ++polygon) |
899 { | 925 { |
900 double x1, y1, x2, y2; | 926 double x1, y1, x2, y2; |
901 | 927 |
902 if (polygon->Project(x1, y1, x2, y2, slice, GetEstimatedNormal(), GetEstimatedSliceThickness())) | 928 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness())) |
903 { | 929 { |
904 rectangles.push_back(Extent2D(x1, y1, x2, y2)); | 930 rectangles.push_back(Extent2D(x1, y1, x2, y2)); |
905 } | 931 } |
906 } | 932 } |
907 | 933 |
927 } | 953 } |
928 } | 954 } |
929 | 955 |
930 | 956 |
931 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, | 957 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, |
932 const CoordinateSystem3D& plane, | 958 const CoordinateSystem3D& cuttingPlane, |
933 size_t structureIndex, | 959 size_t structureIndex, |
934 const Color& color) const | 960 const Color& color) const |
935 { | 961 { |
936 std::vector< std::vector<ScenePoint2D> > chains; | 962 std::vector< std::vector<ScenePoint2D> > chains; |
937 | 963 |
938 if (ProjectStructure(chains, structureIndex, plane)) | 964 if (ProjectStructure(chains, structureIndex, cuttingPlane)) |
939 { | 965 { |
940 for (size_t j = 0; j < chains.size(); j++) | 966 for (size_t j = 0; j < chains.size(); j++) |
941 { | 967 { |
942 layer.AddChain(chains[j], false, color.GetRed(), color.GetGreen(), color.GetBlue()); | 968 layer.AddChain(chains[j], false, color.GetRed(), color.GetGreen(), color.GetBlue()); |
943 } | 969 } |