Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet.cpp @ 131:3e6163a53b16 wasm
optimization
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Thu, 16 Nov 2017 15:00:45 +0100 |
parents | c9e88e7935a4 |
children | 35c2b85836ce |
comparison
equal
deleted
inserted
replaced
130:1982d6c1d2ff | 131:3e6163a53b16 |
---|---|
129 | 129 |
130 void DicomStructureSet::Polygon::CheckPoint(const Vector& v) | 130 void DicomStructureSet::Polygon::CheckPoint(const Vector& v) |
131 { | 131 { |
132 if (hasSlice_) | 132 if (hasSlice_) |
133 { | 133 { |
134 if (!GeometryToolbox::IsNear(GeometryToolbox::ProjectAlongNormal(v, normal_), | 134 if (!GeometryToolbox::IsNear(GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()), |
135 projectionAlongNormal_, | 135 projectionAlongNormal_, |
136 sliceThickness_ / 2.0 /* in mm */)) | 136 sliceThickness_ / 2.0 /* in mm */)) |
137 { | 137 { |
138 LOG(ERROR) << "This RT-STRUCT contains a point that is off the slice of its instance"; | 138 LOG(ERROR) << "This RT-STRUCT contains a point that is off the slice of its instance"; |
139 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 139 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
166 else | 166 else |
167 { | 167 { |
168 const CoordinateSystem3D& geometry = it->second.geometry_; | 168 const CoordinateSystem3D& geometry = it->second.geometry_; |
169 | 169 |
170 hasSlice_ = true; | 170 hasSlice_ = true; |
171 normal_ = geometry.GetNormal(); | 171 geometry_ = geometry; |
172 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), normal_); | 172 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); |
173 sliceThickness_ = it->second.thickness_; | 173 sliceThickness_ = it->second.thickness_; |
174 | 174 |
175 extent_.Reset(); | |
176 | |
175 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) | 177 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) |
176 { | 178 { |
177 CheckPoint(*it); | 179 CheckPoint(*it); |
180 | |
181 double x, y; | |
182 geometry.ProjectPoint(x, y, *it); | |
183 extent_.AddPoint(x, y); | |
178 } | 184 } |
179 | 185 |
180 return true; | 186 return true; |
181 } | 187 } |
182 } | 188 } |
187 { | 193 { |
188 bool isOpposite; | 194 bool isOpposite; |
189 | 195 |
190 if (points_.empty() || | 196 if (points_.empty() || |
191 !hasSlice_ || | 197 !hasSlice_ || |
192 !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), normal_)) | 198 !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) |
193 { | 199 { |
194 return false; | 200 return false; |
195 } | 201 } |
196 | 202 |
197 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), normal_); | 203 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); |
198 | 204 |
199 return (GeometryToolbox::IsNear(d, projectionAlongNormal_, | 205 return (GeometryToolbox::IsNear(d, projectionAlongNormal_, |
200 sliceThickness_ / 2.0)); | 206 sliceThickness_ / 2.0)); |
207 } | |
208 | |
209 | |
210 bool DicomStructureSet::Polygon::Project(double& x1, | |
211 double& y1, | |
212 double& x2, | |
213 double& y2, | |
214 const CoordinateSystem3D& slice) const | |
215 { | |
216 if (!hasSlice_ || | |
217 points_.empty()) | |
218 { | |
219 return false; | |
220 } | |
221 | |
222 double x, y; | |
223 geometry_.ProjectPoint(x, y, slice.GetOrigin()); | |
224 | |
225 bool isOpposite; | |
226 if (GeometryToolbox::IsParallelOrOpposite | |
227 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | |
228 { | |
229 if (y >= extent_.GetY1() && | |
230 y <= extent_.GetY2()) | |
231 { | |
232 Vector p1 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX1(), y) + | |
233 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
234 Vector p2 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX2(), y) - | |
235 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
236 | |
237 slice.ProjectPoint(x1, y1, p1); | |
238 slice.ProjectPoint(x2, y2, p2); | |
239 return true; | |
240 } | |
241 else | |
242 { | |
243 return false; | |
244 } | |
245 } | |
246 else if (GeometryToolbox::IsParallelOrOpposite | |
247 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | |
248 { | |
249 if (x >= extent_.GetX1() && | |
250 x <= extent_.GetX2()) | |
251 { | |
252 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY1()) + | |
253 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
254 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY2()) - | |
255 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
256 | |
257 slice.ProjectPoint(x1, y1, p1); | |
258 slice.ProjectPoint(x2, y2, p2); | |
259 return true; | |
260 } | |
261 else | |
262 { | |
263 return false; | |
264 } | |
265 } | |
266 else | |
267 { | |
268 // Should not happen | |
269 return false; | |
270 } | |
201 } | 271 } |
202 | 272 |
203 | 273 |
204 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const | 274 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const |
205 { | 275 { |
447 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 517 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
448 } | 518 } |
449 } | 519 } |
450 | 520 |
451 referencedSlices_[sopInstanceUid] = ReferencedSlice(seriesInstanceUid, geometry, thickness); | 521 referencedSlices_[sopInstanceUid] = ReferencedSlice(seriesInstanceUid, geometry, thickness); |
522 | |
523 for (Structures::iterator structure = structures_.begin(); | |
524 structure != structures_.end(); ++structure) | |
525 { | |
526 for (Polygons::iterator polygon = structure->polygons_.begin(); | |
527 polygon != structure->polygons_.end(); ++polygon) | |
528 { | |
529 polygon->UpdateReferencedSlice(referencedSlices_); | |
530 } | |
531 } | |
452 } | 532 } |
453 } | 533 } |
454 | 534 |
455 | 535 |
456 void DicomStructureSet::AddReferencedSlice(const Orthanc::DicomMap& dataset) | 536 void DicomStructureSet::AddReferencedSlice(const Orthanc::DicomMap& dataset) |
570 // This is an axial projection | 650 // This is an axial projection |
571 | 651 |
572 for (Polygons::iterator polygon = structure.polygons_.begin(); | 652 for (Polygons::iterator polygon = structure.polygons_.begin(); |
573 polygon != structure.polygons_.end(); ++polygon) | 653 polygon != structure.polygons_.end(); ++polygon) |
574 { | 654 { |
575 polygon->UpdateReferencedSlice(referencedSlices_); | |
576 | |
577 if (polygon->IsOnSlice(slice)) | 655 if (polygon->IsOnSlice(slice)) |
578 { | 656 { |
579 polygons.push_back(std::vector<PolygonPoint>()); | 657 polygons.push_back(std::vector<PolygonPoint>()); |
580 | 658 |
581 for (Points::const_iterator p = polygon->GetPoints().begin(); | 659 for (Points::const_iterator p = polygon->GetPoints().begin(); |
592 } | 670 } |
593 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | 671 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || |
594 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 672 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) |
595 { | 673 { |
596 // Sagittal or coronal projection | 674 // Sagittal or coronal projection |
597 | |
598 std::vector<BoostPolygon> projected; | 675 std::vector<BoostPolygon> projected; |
599 | 676 |
600 for (Polygons::iterator polygon = structure.polygons_.begin(); | 677 for (Polygons::iterator polygon = structure.polygons_.begin(); |
601 polygon != structure.polygons_.end(); ++polygon) | 678 polygon != structure.polygons_.end(); ++polygon) |
602 { | 679 { |
603 polygon->UpdateReferencedSlice(referencedSlices_); | 680 double x1, y1, x2, y2; |
604 | 681 if (polygon->Project(x1, y1, x2, y2, slice)) |
605 // First, check that the polygon intersects the input slice | 682 { |
606 double zmin = std::numeric_limits<double>::infinity(); | 683 projected.push_back(CreateRectangle(x1, y1, x2, y2)); |
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 } | 684 } |
653 } | 685 } |
654 | 686 |
655 BoostMultiPolygon merged; | 687 BoostMultiPolygon merged; |
656 Union(merged, projected); | 688 Union(merged, projected); |