Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructureSet.cpp @ 126:c9e88e7935a4 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
rt-struct projection
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 08 Nov 2017 17:27:38 +0100
parents 44fc253d4876
children 3e6163a53b16
comparison
equal deleted inserted replaced
125:44fc253d4876 126:c9e88e7935a4
26 #include <Core/Logging.h> 26 #include <Core/Logging.h>
27 #include <Core/OrthancException.h> 27 #include <Core/OrthancException.h>
28 #include <Plugins/Samples/Common/FullOrthancDataset.h> 28 #include <Plugins/Samples/Common/FullOrthancDataset.h>
29 #include <Plugins/Samples/Common/DicomDatasetReader.h> 29 #include <Plugins/Samples/Common/DicomDatasetReader.h>
30 30
31 #include <limits>
31 #include <stdio.h> 32 #include <stdio.h>
32 #include <boost/lexical_cast.hpp> 33 #include <boost/lexical_cast.hpp>
34 #include <boost/geometry.hpp>
35 #include <boost/geometry/geometries/point_xy.hpp>
36 #include <boost/geometry/geometries/polygon.hpp>
37
38 typedef boost::geometry::model::d2::point_xy<double> BoostPoint;
39 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon;
40 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon;
41
42
43 static void Union(BoostMultiPolygon& output,
44 std::vector<BoostPolygon>& input)
45 {
46 for (size_t i = 0; i < input.size(); i++)
47 {
48 boost::geometry::correct(input[i]);
49 }
50
51 if (input.size() == 0)
52 {
53 output.clear();
54 }
55 else if (input.size() == 1)
56 {
57 output.resize(1);
58 output[0] = input[0];
59 }
60 else
61 {
62 boost::geometry::union_(input[0], input[1], output);
63
64 for (size_t i = 0; i < input.size(); i++)
65 {
66 BoostMultiPolygon tmp;
67 boost::geometry::union_(output, input[i], tmp);
68 output = tmp;
69 }
70 }
71 }
72
73
74 static BoostPolygon CreateRectangle(float x1, float y1,
75 float x2, float y2)
76 {
77 BoostPolygon r;
78 boost::geometry::append(r, BoostPoint(x1, y1));
79 boost::geometry::append(r, BoostPoint(x1, y2));
80 boost::geometry::append(r, BoostPoint(x2, y2));
81 boost::geometry::append(r, BoostPoint(x2, y1));
82 return r;
83 }
84
85
33 86
34 namespace OrthancStone 87 namespace OrthancStone
35 { 88 {
36 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); 89 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042);
37 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); 90 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016);
447 } 500 }
448 501
449 502
450 Vector DicomStructureSet::GetNormal() const 503 Vector DicomStructureSet::GetNormal() const
451 { 504 {
452 if (!referencedSlices_.empty()) 505 if (referencedSlices_.empty())
453 { 506 {
454 Vector v; 507 Vector v;
455 GeometryToolbox::AssignVector(v, 0, 0, 1); 508 GeometryToolbox::AssignVector(v, 0, 0, 1);
456 return v; 509 return v;
457 } 510 }
538 return true; 591 return true;
539 } 592 }
540 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || 593 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) ||
541 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) 594 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY()))
542 { 595 {
543 // Sagittal or coronal projections 596 // Sagittal or coronal projection
597
598 std::vector<BoostPolygon> projected;
599
600 for (Polygons::iterator polygon = structure.polygons_.begin();
601 polygon != structure.polygons_.end(); ++polygon)
602 {
603 polygon->UpdateReferencedSlice(referencedSlices_);
604
605 // First, check that the polygon intersects the input slice
606 double zmin = std::numeric_limits<double>::infinity();
607 double zmax = -std::numeric_limits<double>::infinity();
608
609 double zref = slice.ProjectAlongNormal(slice.GetOrigin());
610
611 for (Points::const_iterator p = polygon->GetPoints().begin();
612 p != polygon->GetPoints().end(); ++p)
613 {
614 double z = slice.ProjectAlongNormal(*p);
615 zmin = std::min(zmin, z);
616 zmax = std::max(zmax, z);
617 }
618
619 if (zmin < zref && zref < zmax)
620 {
621 // The polygon intersect the input slice, project the polygon
622 double xmin = std::numeric_limits<double>::infinity();
623 double xmax = -std::numeric_limits<double>::infinity();
624 double ymin = std::numeric_limits<double>::infinity();
625 double ymax = -std::numeric_limits<double>::infinity();
626
627 for (Points::const_iterator p = polygon->GetPoints().begin();
628 p != polygon->GetPoints().end(); ++p)
629 {
630 double x, y;
631 slice.ProjectPoint(x, y, *p);
632 xmin = std::min(xmin, x);
633 xmax = std::max(xmax, x);
634 ymin = std::min(ymin, y);
635 ymax = std::max(ymax, y);
636 }
637
638 if (GeometryToolbox::IsNear(xmin, xmax))
639 {
640 projected.push_back(CreateRectangle(xmin - polygon->GetSliceThickness() / 2.0, ymin,
641 xmax + polygon->GetSliceThickness() / 2.0, ymax));
642 }
643 else if (GeometryToolbox::IsNear(ymin, ymax))
644 {
645 projected.push_back(CreateRectangle(xmin, ymin - polygon->GetSliceThickness() / 2.0,
646 xmax, ymax + polygon->GetSliceThickness() / 2.0));
647 }
648 else
649 {
650 // Should not happen
651 }
652 }
653 }
654
655 BoostMultiPolygon merged;
656 Union(merged, projected);
657
658 polygons.resize(merged.size());
659 for (size_t i = 0; i < merged.size(); i++)
660 {
661 const std::vector<BoostPoint>& outer = merged[i].outer();
662
663 polygons[i].resize(outer.size());
664 for (size_t j = 0; j < outer.size(); j++)
665 {
666 polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y());
667 }
668 }
544 669
670 return true;
671 }
672 else
673 {
545 return false; 674 return false;
546 } 675 }
547 else
548 {
549 return false;
550 }
551 } 676 }
552 } 677 }