Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructureSet.cpp @ 981:c20dbaab360c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Ability to cope with empty "Referenced SOP Instance UID" (dicom path (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)) + better logs + code formating
author Benjamin Golinvaux <bgo@osimis.io>
date Fri, 06 Sep 2019 09:38:18 +0200
parents 13e078adfb94
children 4c9b4c4de814
comparison
equal deleted inserted replaced
978:0e21ecafcc23 981:c20dbaab360c
23 23
24 #include "../Toolbox/GeometryToolbox.h" 24 #include "../Toolbox/GeometryToolbox.h"
25 25
26 #include <Core/Logging.h> 26 #include <Core/Logging.h>
27 #include <Core/OrthancException.h> 27 #include <Core/OrthancException.h>
28 #include <Core/Toolbox.h>
28 #include <Plugins/Samples/Common/FullOrthancDataset.h> 29 #include <Plugins/Samples/Common/FullOrthancDataset.h>
29 #include <Plugins/Samples/Common/DicomDatasetReader.h> 30 #include <Plugins/Samples/Common/DicomDatasetReader.h>
30 31
31 #include <limits> 32 #include <limits>
32 #include <stdio.h> 33 #include <stdio.h>
124 std::string value; 125 std::string value;
125 return (dataset.GetStringValue(value, tag) && 126 return (dataset.GetStringValue(value, tag) &&
126 LinearAlgebra::ParseVector(target, value)); 127 LinearAlgebra::ParseVector(target, value));
127 } 128 }
128 129
129 130 void DicomStructureSet::Polygon::CheckPointIsOnSlice(const Vector& v) const
130 void DicomStructureSet::Polygon::CheckPoint(const Vector& v)
131 { 131 {
132 if (hasSlice_) 132 if (hasSlice_)
133 { 133 {
134 if (!LinearAlgebra::IsNear(GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()), 134 double magnitude =
135 projectionAlongNormal_, 135 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal());
136 sliceThickness_ / 2.0 /* in mm */)) 136 if(!LinearAlgebra::IsNear(
137 { 137 magnitude,
138 LOG(ERROR) << "This RT-STRUCT contains a point that is off the slice of its instance"; 138 projectionAlongNormal_,
139 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); 139 sliceThickness_ / 2.0 /* in mm */ ))
140 } 140 {
141 } 141 LOG(ERROR) << "This RT-STRUCT contains a point that is off the "
142 } 142 << "slice of its instance | "
143 143 << "magnitude = " << magnitude << " | "
144 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | "
145 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0);
146
147 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
148 }
149 }
150 }
151
152 bool DicomStructureSet::Polygon::IsPointOnSlice(const Vector& v) const
153 {
154 if (hasSlice_)
155 {
156 double magnitude =
157 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal());
158 bool onSlice = LinearAlgebra::IsNear(
159 magnitude,
160 projectionAlongNormal_,
161 sliceThickness_ / 2.0 /* in mm */);
162 if (!onSlice)
163 {
164 LOG(WARNING) << "This RT-STRUCT contains a point that is off the "
165 << "slice of its instance | "
166 << "magnitude = " << magnitude << " | "
167 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | "
168 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0);
169 }
170 return onSlice;
171 }
172 else
173 {
174 return false;
175 }
176 }
144 177
145 void DicomStructureSet::Polygon::AddPoint(const Vector& v) 178 void DicomStructureSet::Polygon::AddPoint(const Vector& v)
146 { 179 {
180 #if 1
181 // BGO 2019-09-03
182 if (IsPointOnSlice(v))
183 {
184 points_.push_back(v);
185 }
186 #else
147 CheckPoint(v); 187 CheckPoint(v);
148 points_.push_back(v); 188 points_.push_back(v);
189 #endif
149 } 190 }
150 191
151 192
152 bool DicomStructureSet::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices) 193 bool DicomStructureSet::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices)
153 { 194 {
174 215
175 extent_.Reset(); 216 extent_.Reset();
176 217
177 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) 218 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it)
178 { 219 {
179 CheckPoint(*it); 220 if (IsPointOnSlice(*it))
180 221 {
181 double x, y; 222 double x, y;
182 geometry.ProjectPoint(x, y, *it); 223 geometry.ProjectPoint(x, y, *it);
183 extent_.AddPoint(x, y); 224 extent_.AddPoint(x, y);
184 } 225 }
185 226 }
186 return true; 227 return true;
187 } 228 }
188 } 229 }
189 } 230 }
190 231
191
192 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const 232 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const
193 { 233 {
194 bool isOpposite; 234 bool isOpposite = false;
195 235
196 if (points_.empty() || 236 if (points_.empty() ||
197 !hasSlice_ || 237 !hasSlice_ ||
198 !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) 238 !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal()))
199 { 239 {
203 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); 243 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal());
204 244
205 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, 245 return (LinearAlgebra::IsNear(d, projectionAlongNormal_,
206 sliceThickness_ / 2.0)); 246 sliceThickness_ / 2.0));
207 } 247 }
208 248
209
210 bool DicomStructureSet::Polygon::Project(double& x1, 249 bool DicomStructureSet::Polygon::Project(double& x1,
211 double& y1, 250 double& y1,
212 double& x2, 251 double& x2,
213 double& y2, 252 double& y2,
214 const CoordinateSystem3D& slice) const 253 const CoordinateSystem3D& slice) const
215 { 254 {
216 // TODO Optimize this method using a sweep-line algorithm for polygons 255 // TODO: optimize this method using a sweep-line algorithm for polygons
217 256
218 if (!hasSlice_ || 257 if (!hasSlice_ ||
219 points_.size() <= 1) 258 points_.size() <= 1)
220 { 259 {
221 return false; 260 return false;
430 469
431 OrthancPlugins::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, 470 OrthancPlugins::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
432 DICOM_TAG_CONTOUR_SEQUENCE, 0, 471 DICOM_TAG_CONTOUR_SEQUENCE, 0,
433 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); 472 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE);
434 473
474 // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)
435 OrthancPlugins::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, 475 OrthancPlugins::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
436 DICOM_TAG_CONTOUR_SEQUENCE, 0, 476 DICOM_TAG_CONTOUR_SEQUENCE, 0,
437 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, 477 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0,
438 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); 478 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID);
439 479
462 } 502 }
463 503
464 size_t size; 504 size_t size;
465 505
466 imageSequencePath.SetPrefixIndex(1, j); 506 imageSequencePath.SetPrefixIndex(1, j);
467 if (!tags.GetSequenceSize(size, imageSequencePath) || 507 if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1)
468 size != 1) 508 {
469 { 509 LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry.";
470 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); 510 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
471 } 511 }
472 512
473 referencedInstancePath.SetPrefixIndex(1, j); 513 referencedInstancePath.SetPrefixIndex(1, j);
474 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); 514 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath);
479 Vector points; 519 Vector points;
480 if (!LinearAlgebra::ParseVector(points, slicesData) || 520 if (!LinearAlgebra::ParseVector(points, slicesData) ||
481 points.size() != 3 * countPoints) 521 points.size() != 3 * countPoints)
482 { 522 {
483 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); 523 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
524 }
525
526 // seen in real world
527 if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "")
528 {
529 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)";
484 } 530 }
485 531
486 Polygon polygon(sopInstanceUid); 532 Polygon polygon(sopInstanceUid);
487 polygon.Reserve(countPoints); 533 polygon.Reserve(countPoints);
488 534
579 { 625 {
580 if (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()) 626 if (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end())
581 { 627 {
582 // This geometry is already known 628 // This geometry is already known
583 LOG(ERROR) << "DicomStructureSet::AddReferencedSlice(): (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()). sopInstanceUid = " << sopInstanceUid; 629 LOG(ERROR) << "DicomStructureSet::AddReferencedSlice(): (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()). sopInstanceUid = " << sopInstanceUid;
584 630
585 // TODO: the following assertion has been disabled on 20190822 by BGO 631 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
586 // because it occurred from time to time. Since it wrecked havoc on the
587
588 //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
589 } 632 }
590 else 633 else
591 { 634 {
592 if (thickness < 0) 635 if (thickness < 0)
593 { 636 {
662 for (Polygons::iterator polygon = structure->polygons_.begin(); 705 for (Polygons::iterator polygon = structure->polygons_.begin();
663 polygon != structure->polygons_.end(); ++polygon) 706 polygon != structure->polygons_.end(); ++polygon)
664 { 707 {
665 if (!polygon->UpdateReferencedSlice(referencedSlices_)) 708 if (!polygon->UpdateReferencedSlice(referencedSlices_))
666 { 709 {
667 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): missing information about referenced instance: " 710 std::string sopInstanceUid = polygon->GetSopInstanceUid();
668 << polygon->GetSopInstanceUid(); 711 if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "")
669 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); 712 {
713 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): "
714 << " missing information about referenced instance "
715 << "(sopInstanceUid is empty!)";
716 }
717 else
718 {
719 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): "
720 << " missing information about referenced instance "
721 << "(sopInstanceUid = " << sopInstanceUid << ")";
722 }
723 //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
670 } 724 }
671 } 725 }
672 } 726 }
673 } 727 }
674 728
686 return referencedSlices_.begin()->second.geometry_.GetNormal(); 740 return referencedSlices_.begin()->second.geometry_.GetNormal();
687 } 741 }
688 } 742 }
689 743
690 744
691 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<PolygonPoint> >& polygons, 745 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons,
692 const Structure& structure, 746 const Structure& structure,
693 const CoordinateSystem3D& slice) const 747 const CoordinateSystem3D& slice) const
694 { 748 {
695 polygons.clear(); 749 polygons.clear();
696 750
704 for (Polygons::const_iterator polygon = structure.polygons_.begin(); 758 for (Polygons::const_iterator polygon = structure.polygons_.begin();
705 polygon != structure.polygons_.end(); ++polygon) 759 polygon != structure.polygons_.end(); ++polygon)
706 { 760 {
707 if (polygon->IsOnSlice(slice)) 761 if (polygon->IsOnSlice(slice))
708 { 762 {
709 polygons.push_back(std::vector<PolygonPoint>()); 763 polygons.push_back(std::vector<PolygonPoint2D>());
710 764
711 for (Points::const_iterator p = polygon->GetPoints().begin(); 765 for (Points::const_iterator p = polygon->GetPoints().begin();
712 p != polygon->GetPoints().end(); ++p) 766 p != polygon->GetPoints().end(); ++p)
713 { 767 {
714 double x, y; 768 double x, y;
760 polygon != structure.polygons_.end(); ++polygon) 814 polygon != structure.polygons_.end(); ++polygon)
761 { 815 {
762 double x1, y1, x2, y2; 816 double x1, y1, x2, y2;
763 if (polygon->Project(x1, y1, x2, y2, slice)) 817 if (polygon->Project(x1, y1, x2, y2, slice))
764 { 818 {
765 std::vector<PolygonPoint> p(4); 819 std::vector<PolygonPoint2D> p(4);
766 p[0] = std::make_pair(x1, y1); 820 p[0] = std::make_pair(x1, y1);
767 p[1] = std::make_pair(x2, y1); 821 p[1] = std::make_pair(x2, y1);
768 p[2] = std::make_pair(x2, y2); 822 p[2] = std::make_pair(x2, y2);
769 p[3] = std::make_pair(x1, y2); 823 p[3] = std::make_pair(x1, y2);
770 polygons.push_back(p); 824 polygons.push_back(p);