Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1512:244ad1e4e76a
reorganization of folders
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 07 Jul 2020 16:21:02 +0200 |
parents | Framework/Toolbox/DicomStructureSet.cpp@d8af188ab545 |
children | 85e117739eca |
comparison
equal
deleted
inserted
replaced
1511:9dfeee74c1e6 | 1512:244ad1e4e76a |
---|---|
1 /** | |
2 * Stone of Orthanc | |
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics | |
4 * Department, University Hospital of Liege, Belgium | |
5 * Copyright (C) 2017-2020 Osimis S.A., Belgium | |
6 * | |
7 * This program is free software: you can redistribute it and/or | |
8 * modify it under the terms of the GNU Affero General Public License | |
9 * as published by the Free Software Foundation, either version 3 of | |
10 * the License, or (at your option) any later version. | |
11 * | |
12 * This program is distributed in the hope that it will be useful, but | |
13 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 * Affero General Public License for more details. | |
16 * | |
17 * You should have received a copy of the GNU Affero General Public License | |
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. | |
19 **/ | |
20 | |
21 | |
22 #include "DicomStructureSet.h" | |
23 #include "DicomStructureSetUtils.h" | |
24 | |
25 #include "../Toolbox/GeometryToolbox.h" | |
26 #include "OrthancDatasets/DicomDatasetReader.h" | |
27 | |
28 #include <Logging.h> | |
29 #include <OrthancException.h> | |
30 #include <Toolbox.h> | |
31 | |
32 #if defined(_MSC_VER) | |
33 # pragma warning(push) | |
34 # pragma warning(disable:4244) | |
35 #endif | |
36 | |
37 #include <limits> | |
38 #include <stdio.h> | |
39 #include <boost/geometry.hpp> | |
40 #include <boost/geometry/geometries/point_xy.hpp> | |
41 #include <boost/geometry/geometries/polygon.hpp> | |
42 #include <boost/geometry/multi/geometries/multi_polygon.hpp> | |
43 | |
44 #if defined(_MSC_VER) | |
45 # pragma warning(pop) | |
46 #endif | |
47 | |
48 #if ORTHANC_ENABLE_DCMTK == 1 | |
49 # include "ParsedDicomDataset.h" | |
50 #endif | |
51 | |
52 | |
53 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; | |
54 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; | |
55 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; | |
56 | |
57 | |
58 static void Union(BoostMultiPolygon& output, | |
59 std::vector<BoostPolygon>& input) | |
60 { | |
61 for (size_t i = 0; i < input.size(); i++) | |
62 { | |
63 boost::geometry::correct(input[i]); | |
64 } | |
65 | |
66 if (input.size() == 0) | |
67 { | |
68 output.clear(); | |
69 } | |
70 else if (input.size() == 1) | |
71 { | |
72 output.resize(1); | |
73 output[0] = input[0]; | |
74 } | |
75 else | |
76 { | |
77 boost::geometry::union_(input[0], input[1], output); | |
78 | |
79 for (size_t i = 0; i < input.size(); i++) | |
80 { | |
81 BoostMultiPolygon tmp; | |
82 boost::geometry::union_(output, input[i], tmp); | |
83 output = tmp; | |
84 } | |
85 } | |
86 } | |
87 | |
88 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
89 | |
90 static BoostPolygon CreateRectangle(float x1, float y1, | |
91 float x2, float y2) | |
92 { | |
93 BoostPolygon r; | |
94 boost::geometry::append(r, BoostPoint(x1, y1)); | |
95 boost::geometry::append(r, BoostPoint(x1, y2)); | |
96 boost::geometry::append(r, BoostPoint(x2, y2)); | |
97 boost::geometry::append(r, BoostPoint(x2, y1)); | |
98 return r; | |
99 } | |
100 | |
101 #else | |
102 | |
103 namespace OrthancStone | |
104 { | |
105 static RtStructRectangleInSlab CreateRectangle(float x1, float y1, | |
106 float x2, float y2) | |
107 { | |
108 RtStructRectangleInSlab rect; | |
109 rect.xmin = std::min(x1, x2); | |
110 rect.xmax = std::max(x1, x2); | |
111 rect.ymin = std::min(y1, y2); | |
112 rect.ymax = std::max(y1, y2); | |
113 return rect; | |
114 } | |
115 | |
116 bool CompareRectanglesForProjection(const std::pair<RtStructRectangleInSlab,double>& r1, const std::pair<RtStructRectangleInSlab, double>& r2) | |
117 { | |
118 return r1.second < r2.second; | |
119 } | |
120 | |
121 bool CompareSlabsY(const RtStructRectanglesInSlab& r1, const RtStructRectanglesInSlab& r2) | |
122 { | |
123 if ((r1.size() == 0) || (r2.size() == 0)) | |
124 return false; | |
125 | |
126 return r1[0].ymax < r2[0].ymax; | |
127 } | |
128 } | |
129 | |
130 #endif | |
131 | |
132 namespace OrthancStone | |
133 { | |
134 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); | |
135 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); | |
136 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_SEQUENCE(0x3006, 0x0040); | |
137 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_DATA(0x3006, 0x0050); | |
138 static const Orthanc::DicomTag DICOM_TAG_NUMBER_OF_CONTOUR_POINTS(0x3006, 0x0046); | |
139 static const Orthanc::DicomTag DICOM_TAG_REFERENCED_SOP_INSTANCE_UID(0x0008, 0x1155); | |
140 static const Orthanc::DicomTag DICOM_TAG_ROI_CONTOUR_SEQUENCE(0x3006, 0x0039); | |
141 static const Orthanc::DicomTag DICOM_TAG_ROI_DISPLAY_COLOR(0x3006, 0x002a); | |
142 static const Orthanc::DicomTag DICOM_TAG_ROI_NAME(0x3006, 0x0026); | |
143 static const Orthanc::DicomTag DICOM_TAG_RT_ROI_INTERPRETED_TYPE(0x3006, 0x00a4); | |
144 static const Orthanc::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080); | |
145 static const Orthanc::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020); | |
146 | |
147 | |
148 static uint8_t ConvertColor(double v) | |
149 { | |
150 if (v < 0) | |
151 { | |
152 return 0; | |
153 } | |
154 else if (v >= 255) | |
155 { | |
156 return 255; | |
157 } | |
158 else | |
159 { | |
160 return static_cast<uint8_t>(v); | |
161 } | |
162 } | |
163 | |
164 | |
165 static bool ParseVector(Vector& target, | |
166 const IDicomDataset& dataset, | |
167 const DicomPath& tag) | |
168 { | |
169 std::string value; | |
170 return (dataset.GetStringValue(value, tag) && | |
171 LinearAlgebra::ParseVector(target, value)); | |
172 } | |
173 | |
174 void DicomStructureSet::Polygon::CheckPointIsOnSlice(const Vector& v) const | |
175 { | |
176 if (hasSlice_) | |
177 { | |
178 double magnitude = | |
179 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); | |
180 if(!LinearAlgebra::IsNear( | |
181 magnitude, | |
182 projectionAlongNormal_, | |
183 sliceThickness_ / 2.0 /* in mm */ )) | |
184 { | |
185 LOG(ERROR) << "This RT-STRUCT contains a point that is off the " | |
186 << "slice of its instance | " | |
187 << "magnitude = " << magnitude << " | " | |
188 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " | |
189 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); | |
190 | |
191 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
192 } | |
193 } | |
194 } | |
195 | |
196 bool DicomStructureSet::Polygon::IsPointOnSliceIfAny(const Vector& v) const | |
197 { | |
198 if (hasSlice_) | |
199 { | |
200 double magnitude = | |
201 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); | |
202 bool onSlice = LinearAlgebra::IsNear( | |
203 magnitude, | |
204 projectionAlongNormal_, | |
205 sliceThickness_ / 2.0 /* in mm */); | |
206 if (!onSlice) | |
207 { | |
208 LOG(WARNING) << "This RT-STRUCT contains a point that is off the " | |
209 << "slice of its instance | " | |
210 << "magnitude = " << magnitude << " | " | |
211 << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | " | |
212 << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0); | |
213 } | |
214 return onSlice; | |
215 } | |
216 else | |
217 { | |
218 return true; | |
219 } | |
220 } | |
221 | |
222 void DicomStructureSet::Polygon::AddPoint(const Vector& v) | |
223 { | |
224 #if 1 | |
225 // BGO 2019-09-03 | |
226 if (IsPointOnSliceIfAny(v)) | |
227 { | |
228 points_.push_back(v); | |
229 } | |
230 #else | |
231 CheckPoint(v); | |
232 points_.push_back(v); | |
233 #endif | |
234 } | |
235 | |
236 | |
237 bool DicomStructureSet::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices) | |
238 { | |
239 if (hasSlice_) | |
240 { | |
241 return true; | |
242 } | |
243 else | |
244 { | |
245 ReferencedSlices::const_iterator it = slices.find(sopInstanceUid_); | |
246 | |
247 if (it == slices.end()) | |
248 { | |
249 return false; | |
250 } | |
251 else | |
252 { | |
253 const CoordinateSystem3D& geometry = it->second.geometry_; | |
254 | |
255 hasSlice_ = true; | |
256 geometry_ = geometry; | |
257 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); | |
258 sliceThickness_ = it->second.thickness_; | |
259 | |
260 extent_.Reset(); | |
261 | |
262 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) | |
263 { | |
264 if (IsPointOnSliceIfAny(*it)) | |
265 { | |
266 double x, y; | |
267 geometry.ProjectPoint2(x, y, *it); | |
268 extent_.AddPoint(x, y); | |
269 } | |
270 } | |
271 return true; | |
272 } | |
273 } | |
274 } | |
275 | |
276 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const | |
277 { | |
278 bool isOpposite = false; | |
279 | |
280 if (points_.empty() || | |
281 !hasSlice_ || | |
282 !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) | |
283 { | |
284 return false; | |
285 } | |
286 | |
287 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); | |
288 | |
289 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, | |
290 sliceThickness_ / 2.0)); | |
291 } | |
292 | |
293 bool DicomStructureSet::Polygon::Project(double& x1, | |
294 double& y1, | |
295 double& x2, | |
296 double& y2, | |
297 const CoordinateSystem3D& slice) const | |
298 { | |
299 // TODO: optimize this method using a sweep-line algorithm for polygons | |
300 | |
301 if (!hasSlice_ || | |
302 points_.size() <= 1) | |
303 { | |
304 return false; | |
305 } | |
306 | |
307 double x, y; | |
308 geometry_.ProjectPoint2(x, y, slice.GetOrigin()); | |
309 | |
310 bool isOpposite; | |
311 if (GeometryToolbox::IsParallelOrOpposite | |
312 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | |
313 { | |
314 // plane is constant Y | |
315 | |
316 if (y < extent_.GetY1() || | |
317 y > extent_.GetY2()) | |
318 { | |
319 // The polygon does not intersect the input slice | |
320 return false; | |
321 } | |
322 | |
323 bool isFirst = true; | |
324 double xmin = std::numeric_limits<double>::infinity(); | |
325 double xmax = -std::numeric_limits<double>::infinity(); | |
326 | |
327 double prevX, prevY; | |
328 geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); | |
329 | |
330 for (size_t i = 0; i < points_.size(); i++) | |
331 { | |
332 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py | |
333 double curX, curY; | |
334 geometry_.ProjectPoint2(curX, curY, points_[i]); | |
335 | |
336 // if prev* and cur* are on opposite sides of y, this means that the | |
337 // segment intersects the plane. | |
338 if ((prevY < y && curY > y) || | |
339 (prevY > y && curY < y)) | |
340 { | |
341 double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); | |
342 xmin = std::min(xmin, p); | |
343 xmax = std::max(xmax, p); | |
344 isFirst = false; | |
345 | |
346 // xmin and xmax represent the extent of the rectangle along the | |
347 // intersection between the plane and the polygon geometry | |
348 | |
349 } | |
350 | |
351 prevX = curX; | |
352 prevY = curY; | |
353 } | |
354 | |
355 // if NO segment intersects the plane | |
356 if (isFirst) | |
357 { | |
358 return false; | |
359 } | |
360 else | |
361 { | |
362 // y is the plane y coord in the polygon geometry | |
363 // xmin and xmax are ALSO expressed in the polygon geometry | |
364 | |
365 // let's convert them to 3D world geometry... | |
366 Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + | |
367 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
368 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - | |
369 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
370 | |
371 // then to the cutting plane geometry... | |
372 slice.ProjectPoint2(x1, y1, p1); | |
373 slice.ProjectPoint2(x2, y2, p2); | |
374 return true; | |
375 } | |
376 } | |
377 else if (GeometryToolbox::IsParallelOrOpposite | |
378 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | |
379 { | |
380 // plane is constant X => Sagittal view (remember that in the | |
381 // sagittal projection, the normal must be swapped) | |
382 | |
383 | |
384 /* | |
385 Please read the comments in the section above, by taking into account | |
386 the fact that, in this case, the plane has a constant X, not Y (in | |
387 polygon geometry_ coordinates) | |
388 */ | |
389 | |
390 if (x < extent_.GetX1() || | |
391 x > extent_.GetX2()) | |
392 { | |
393 return false; | |
394 } | |
395 | |
396 bool isFirst = true; | |
397 double ymin = std::numeric_limits<double>::infinity(); | |
398 double ymax = -std::numeric_limits<double>::infinity(); | |
399 | |
400 double prevX, prevY; | |
401 geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); | |
402 | |
403 for (size_t i = 0; i < points_.size(); i++) | |
404 { | |
405 // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py | |
406 double curX, curY; | |
407 geometry_.ProjectPoint2(curX, curY, points_[i]); | |
408 | |
409 if ((prevX < x && curX > x) || | |
410 (prevX > x && curX < x)) | |
411 { | |
412 double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); | |
413 ymin = std::min(ymin, p); | |
414 ymax = std::max(ymax, p); | |
415 isFirst = false; | |
416 } | |
417 | |
418 prevX = curX; | |
419 prevY = curY; | |
420 } | |
421 | |
422 if (isFirst) | |
423 { | |
424 return false; | |
425 } | |
426 else | |
427 { | |
428 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + | |
429 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
430 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - | |
431 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
432 | |
433 slice.ProjectPoint2(x1, y1, p1); | |
434 slice.ProjectPoint2(x2, y2, p2); | |
435 | |
436 return true; | |
437 } | |
438 } | |
439 else | |
440 { | |
441 // Should not happen | |
442 return false; | |
443 } | |
444 } | |
445 | |
446 | |
447 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const | |
448 { | |
449 if (index >= structures_.size()) | |
450 { | |
451 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
452 } | |
453 | |
454 return structures_[index]; | |
455 } | |
456 | |
457 | |
458 DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) | |
459 { | |
460 if (index >= structures_.size()) | |
461 { | |
462 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
463 } | |
464 | |
465 return structures_[index]; | |
466 } | |
467 | |
468 void DicomStructureSet::Setup(const IDicomDataset& tags) | |
469 { | |
470 DicomDatasetReader reader(tags); | |
471 | |
472 size_t count, tmp; | |
473 if (!tags.GetSequenceSize(count, DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE) || | |
474 !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) || | |
475 tmp != count || | |
476 !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) || | |
477 tmp != count) | |
478 { | |
479 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
480 } | |
481 | |
482 structures_.resize(count); | |
483 for (size_t i = 0; i < count; i++) | |
484 { | |
485 structures_[i].interpretation_ = reader.GetStringValue | |
486 (DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, | |
487 DICOM_TAG_RT_ROI_INTERPRETED_TYPE), | |
488 "No interpretation"); | |
489 | |
490 structures_[i].name_ = reader.GetStringValue | |
491 (DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, | |
492 DICOM_TAG_ROI_NAME), | |
493 "No name"); | |
494 | |
495 Vector color; | |
496 if (ParseVector(color, tags, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
497 DICOM_TAG_ROI_DISPLAY_COLOR)) && | |
498 color.size() == 3) | |
499 { | |
500 structures_[i].red_ = ConvertColor(color[0]); | |
501 structures_[i].green_ = ConvertColor(color[1]); | |
502 structures_[i].blue_ = ConvertColor(color[2]); | |
503 } | |
504 else | |
505 { | |
506 structures_[i].red_ = 255; | |
507 structures_[i].green_ = 0; | |
508 structures_[i].blue_ = 0; | |
509 } | |
510 | |
511 size_t countSlices; | |
512 if (!tags.GetSequenceSize(countSlices, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
513 DICOM_TAG_CONTOUR_SEQUENCE))) | |
514 { | |
515 countSlices = 0; | |
516 } | |
517 | |
518 LOG(INFO) << "New RT structure: \"" << structures_[i].name_ | |
519 << "\" with interpretation \"" << structures_[i].interpretation_ | |
520 << "\" containing " << countSlices << " slices (color: " | |
521 << static_cast<int>(structures_[i].red_) << "," | |
522 << static_cast<int>(structures_[i].green_) << "," | |
523 << static_cast<int>(structures_[i].blue_) << ")"; | |
524 | |
525 // These temporary variables avoid allocating many vectors in the loop below | |
526 DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
527 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
528 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); | |
529 | |
530 DicomPath geometricTypePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
531 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
532 DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); | |
533 | |
534 DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
535 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
536 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); | |
537 | |
538 // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155) | |
539 DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
540 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
541 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, | |
542 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); | |
543 | |
544 DicomPath contourDataPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
545 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
546 DICOM_TAG_CONTOUR_DATA); | |
547 | |
548 for (size_t j = 0; j < countSlices; j++) | |
549 { | |
550 unsigned int countPoints; | |
551 | |
552 countPointsPath.SetPrefixIndex(1, j); | |
553 if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath)) | |
554 { | |
555 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
556 } | |
557 | |
558 //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices"; | |
559 | |
560 geometricTypePath.SetPrefixIndex(1, j); | |
561 std::string type = reader.GetMandatoryStringValue(geometricTypePath); | |
562 if (type != "CLOSED_PLANAR") | |
563 { | |
564 LOG(WARNING) << "Ignoring contour with geometry type: " << type; | |
565 continue; | |
566 } | |
567 | |
568 size_t size; | |
569 | |
570 imageSequencePath.SetPrefixIndex(1, j); | |
571 if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1) | |
572 { | |
573 LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry."; | |
574 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
575 } | |
576 | |
577 referencedInstancePath.SetPrefixIndex(1, j); | |
578 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); | |
579 | |
580 contourDataPath.SetPrefixIndex(1, j); | |
581 std::string slicesData = reader.GetMandatoryStringValue(contourDataPath); | |
582 | |
583 Vector points; | |
584 if (!LinearAlgebra::ParseVector(points, slicesData) || | |
585 points.size() != 3 * countPoints) | |
586 { | |
587 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
588 } | |
589 | |
590 // seen in real world | |
591 if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") | |
592 { | |
593 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)"; | |
594 } | |
595 | |
596 Polygon polygon(sopInstanceUid); | |
597 polygon.Reserve(countPoints); | |
598 | |
599 for (size_t k = 0; k < countPoints; k++) | |
600 { | |
601 Vector v(3); | |
602 v[0] = points[3 * k]; | |
603 v[1] = points[3 * k + 1]; | |
604 v[2] = points[3 * k + 2]; | |
605 polygon.AddPoint(v); | |
606 } | |
607 | |
608 structures_[i].polygons_.push_back(polygon); | |
609 } | |
610 } | |
611 } | |
612 | |
613 | |
614 #if ORTHANC_ENABLE_DCMTK == 1 | |
615 DicomStructureSet::DicomStructureSet(Orthanc::ParsedDicomFile& instance) | |
616 { | |
617 ParsedDicomDataset dataset(instance); | |
618 Setup(dataset); | |
619 } | |
620 #endif | |
621 | |
622 | |
623 Vector DicomStructureSet::GetStructureCenter(size_t index) const | |
624 { | |
625 const Structure& structure = GetStructure(index); | |
626 | |
627 Vector center; | |
628 LinearAlgebra::AssignVector(center, 0, 0, 0); | |
629 if (structure.polygons_.empty()) | |
630 { | |
631 return center; | |
632 } | |
633 | |
634 double n = static_cast<double>(structure.polygons_.size()); | |
635 | |
636 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
637 polygon != structure.polygons_.end(); ++polygon) | |
638 { | |
639 if (!polygon->GetPoints().empty()) | |
640 { | |
641 center += polygon->GetPoints().front() / n; | |
642 } | |
643 } | |
644 | |
645 return center; | |
646 } | |
647 | |
648 | |
649 const std::string& DicomStructureSet::GetStructureName(size_t index) const | |
650 { | |
651 return GetStructure(index).name_; | |
652 } | |
653 | |
654 | |
655 const std::string& DicomStructureSet::GetStructureInterpretation(size_t index) const | |
656 { | |
657 return GetStructure(index).interpretation_; | |
658 } | |
659 | |
660 | |
661 Color DicomStructureSet::GetStructureColor(size_t index) const | |
662 { | |
663 const Structure& s = GetStructure(index); | |
664 return Color(s.red_, s.green_, s.blue_); | |
665 } | |
666 | |
667 | |
668 void DicomStructureSet::GetStructureColor(uint8_t& red, | |
669 uint8_t& green, | |
670 uint8_t& blue, | |
671 size_t index) const | |
672 { | |
673 const Structure& s = GetStructure(index); | |
674 red = s.red_; | |
675 green = s.green_; | |
676 blue = s.blue_; | |
677 } | |
678 | |
679 | |
680 void DicomStructureSet::GetReferencedInstances(std::set<std::string>& instances) | |
681 { | |
682 for (Structures::const_iterator structure = structures_.begin(); | |
683 structure != structures_.end(); ++structure) | |
684 { | |
685 for (Polygons::const_iterator polygon = structure->polygons_.begin(); | |
686 polygon != structure->polygons_.end(); ++polygon) | |
687 { | |
688 instances.insert(polygon->GetSopInstanceUid()); | |
689 } | |
690 } | |
691 } | |
692 | |
693 | |
694 void DicomStructureSet::AddReferencedSlice(const std::string& sopInstanceUid, | |
695 const std::string& seriesInstanceUid, | |
696 const CoordinateSystem3D& geometry, | |
697 double thickness) | |
698 { | |
699 if (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()) | |
700 { | |
701 // This geometry is already known | |
702 LOG(ERROR) << "DicomStructureSet::AddReferencedSlice(): (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()). sopInstanceUid = " << sopInstanceUid; | |
703 | |
704 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
705 } | |
706 else | |
707 { | |
708 if (thickness < 0) | |
709 { | |
710 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
711 } | |
712 | |
713 if (!referencedSlices_.empty()) | |
714 { | |
715 const ReferencedSlice& reference = referencedSlices_.begin()->second; | |
716 | |
717 if (reference.seriesInstanceUid_ != seriesInstanceUid) | |
718 { | |
719 LOG(ERROR) << "This RT-STRUCT refers to several different series"; | |
720 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
721 } | |
722 | |
723 if (!GeometryToolbox::IsParallel(reference.geometry_.GetNormal(), geometry.GetNormal())) | |
724 { | |
725 LOG(ERROR) << "The slices in this RT-STRUCT are not parallel"; | |
726 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
727 } | |
728 } | |
729 | |
730 referencedSlices_[sopInstanceUid] = ReferencedSlice(seriesInstanceUid, geometry, thickness); | |
731 | |
732 for (Structures::iterator structure = structures_.begin(); | |
733 structure != structures_.end(); ++structure) | |
734 { | |
735 for (Polygons::iterator polygon = structure->polygons_.begin(); | |
736 polygon != structure->polygons_.end(); ++polygon) | |
737 { | |
738 polygon->UpdateReferencedSlice(referencedSlices_); | |
739 } | |
740 } | |
741 } | |
742 } | |
743 | |
744 | |
745 void DicomStructureSet::AddReferencedSlice(const Orthanc::DicomMap& dataset) | |
746 { | |
747 CoordinateSystem3D slice(dataset); | |
748 | |
749 double thickness = 1; // 1 mm by default | |
750 | |
751 std::string s; | |
752 Vector v; | |
753 if (dataset.LookupStringValue(s, Orthanc::DICOM_TAG_SLICE_THICKNESS, false) && | |
754 LinearAlgebra::ParseVector(v, s) && | |
755 v.size() > 0) | |
756 { | |
757 thickness = v[0]; | |
758 } | |
759 | |
760 std::string instance, series; | |
761 if (dataset.LookupStringValue(instance, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false) && | |
762 dataset.LookupStringValue(series, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false)) | |
763 { | |
764 AddReferencedSlice(instance, series, slice, thickness); | |
765 } | |
766 else | |
767 { | |
768 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
769 } | |
770 } | |
771 | |
772 | |
773 void DicomStructureSet::CheckReferencedSlices() | |
774 { | |
775 for (Structures::iterator structure = structures_.begin(); | |
776 structure != structures_.end(); ++structure) | |
777 { | |
778 for (Polygons::iterator polygon = structure->polygons_.begin(); | |
779 polygon != structure->polygons_.end(); ++polygon) | |
780 { | |
781 if (!polygon->UpdateReferencedSlice(referencedSlices_)) | |
782 { | |
783 std::string sopInstanceUid = polygon->GetSopInstanceUid(); | |
784 if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") | |
785 { | |
786 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): " | |
787 << " missing information about referenced instance " | |
788 << "(sopInstanceUid is empty!)"; | |
789 } | |
790 else | |
791 { | |
792 LOG(ERROR) << "DicomStructureSet::CheckReferencedSlices(): " | |
793 << " missing information about referenced instance " | |
794 << "(sopInstanceUid = " << sopInstanceUid << ")"; | |
795 } | |
796 //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
797 } | |
798 } | |
799 } | |
800 } | |
801 | |
802 | |
803 Vector DicomStructureSet::GetNormal() const | |
804 { | |
805 if (referencedSlices_.empty()) | |
806 { | |
807 Vector v; | |
808 LinearAlgebra::AssignVector(v, 0, 0, 1); | |
809 return v; | |
810 } | |
811 else | |
812 { | |
813 return referencedSlices_.begin()->second.geometry_.GetNormal(); | |
814 } | |
815 } | |
816 | |
817 bool DicomStructureSet::ProjectStructure( | |
818 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
819 std::vector< std::vector<Point2D> >& polygons, | |
820 #else | |
821 std::vector< std::pair<Point2D, Point2D> >& segments, | |
822 #endif | |
823 const Structure& structure, | |
824 const CoordinateSystem3D& sourceSlice) const | |
825 { | |
826 const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); | |
827 | |
828 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
829 polygons.clear(); | |
830 #else | |
831 segments.clear(); | |
832 #endif | |
833 | |
834 Vector normal = GetNormal(); | |
835 | |
836 bool isOpposite; | |
837 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) | |
838 { | |
839 // This is an axial projection | |
840 | |
841 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
842 polygon != structure.polygons_.end(); ++polygon) | |
843 { | |
844 if (polygon->IsOnSlice(slice)) | |
845 { | |
846 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
847 polygons.push_back(std::vector<Point2D>()); | |
848 | |
849 for (Points::const_iterator p = polygon->GetPoints().begin(); | |
850 p != polygon->GetPoints().end(); ++p) | |
851 { | |
852 double x, y; | |
853 slice.ProjectPoint2(x, y, *p); | |
854 polygons.back().push_back(Point2D(x, y)); | |
855 } | |
856 #else | |
857 // we need to add all the segments corresponding to this polygon | |
858 const std::vector<Vector>& points3D = polygon->GetPoints(); | |
859 if (points3D.size() >= 3) | |
860 { | |
861 Point2D prev2D; | |
862 { | |
863 Vector prev = points3D[0]; | |
864 double prevX, prevY; | |
865 slice.ProjectPoint2(prevX, prevY, prev); | |
866 prev2D = Point2D(prevX, prevY); | |
867 } | |
868 | |
869 size_t pointCount = points3D.size(); | |
870 for (size_t ipt = 1; ipt < pointCount; ++ipt) | |
871 { | |
872 Vector next = points3D[ipt]; | |
873 double nextX, nextY; | |
874 slice.ProjectPoint2(nextX, nextY, next); | |
875 Point2D next2D(nextX, nextY); | |
876 segments.push_back(std::pair<Point2D, Point2D>(prev2D, next2D)); | |
877 prev2D = next2D; | |
878 } | |
879 } | |
880 else | |
881 { | |
882 LOG(ERROR) << "Contour with less than 3 points!"; | |
883 // !!! | |
884 } | |
885 #endif | |
886 } | |
887 } | |
888 | |
889 return true; | |
890 } | |
891 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | |
892 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | |
893 { | |
894 #if 1 | |
895 // Sagittal or coronal projection | |
896 | |
897 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
898 std::vector<BoostPolygon> projected; | |
899 | |
900 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
901 polygon != structure.polygons_.end(); ++polygon) | |
902 { | |
903 double x1, y1, x2, y2; | |
904 | |
905 if (polygon->Project(x1, y1, x2, y2, slice)) | |
906 { | |
907 projected.push_back(CreateRectangle(x1, y1, x2, y2)); | |
908 } | |
909 } | |
910 #else | |
911 // this will contain the intersection of the polygon slab with | |
912 // the cutting plane, projected on the cutting plane coord system | |
913 // (that yields a rectangle) + the Z coordinate of the polygon | |
914 // (this is required to group polygons with the same Z later) | |
915 std::vector<std::pair<RtStructRectangleInSlab, double> > projected; | |
916 | |
917 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
918 polygon != structure.polygons_.end(); ++polygon) | |
919 { | |
920 double x1, y1, x2, y2; | |
921 | |
922 if (polygon->Project(x1, y1, x2, y2, slice)) | |
923 { | |
924 double curZ = polygon->GetGeometryOrigin()[2]; | |
925 | |
926 // x1,y1 and x2,y2 are in "slice" coordinates (the cutting plane | |
927 // geometry) | |
928 projected.push_back(std::make_pair(CreateRectangle( | |
929 static_cast<float>(x1), | |
930 static_cast<float>(y1), | |
931 static_cast<float>(x2), | |
932 static_cast<float>(y2)),curZ)); | |
933 } | |
934 } | |
935 #endif | |
936 | |
937 #if USE_BOOST_UNION_FOR_POLYGONS != 1 | |
938 // projected contains a set of rectangles specified by two opposite | |
939 // corners (x1,y1,x2,y2) | |
940 // we need to merge them | |
941 // each slab yields ONE polygon! | |
942 | |
943 // we need to sorted all the rectangles that originate from the same Z | |
944 // into lanes. To make sure they are grouped together in the array, we | |
945 // sort it. | |
946 std::sort(projected.begin(), projected.end(), CompareRectanglesForProjection); | |
947 | |
948 std::vector<RtStructRectanglesInSlab> rectanglesForEachSlab; | |
949 rectanglesForEachSlab.reserve(projected.size()); | |
950 | |
951 double curZ = 0; | |
952 for (size_t i = 0; i < projected.size(); ++i) | |
953 { | |
954 #if 0 | |
955 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
956 #else | |
957 if (i == 0) | |
958 { | |
959 curZ = projected[i].second; | |
960 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
961 } | |
962 else | |
963 { | |
964 // this check is needed to prevent creating a new slab if | |
965 // the new polygon is at the same Z coord than last one | |
966 if (!LinearAlgebra::IsNear(curZ, projected[i].second)) | |
967 { | |
968 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
969 curZ = projected[i].second; | |
970 } | |
971 } | |
972 #endif | |
973 | |
974 rectanglesForEachSlab.back().push_back(projected[i].first); | |
975 | |
976 // as long as they have the same y, we should put them into the same lane | |
977 // BUT in Sebastien's code, there is only one polygon per lane. | |
978 | |
979 //std::cout << "rect: xmin = " << rect.xmin << " xmax = " << rect.xmax << " ymin = " << rect.ymin << " ymax = " << rect.ymax << std::endl; | |
980 } | |
981 | |
982 // now we need to sort the slabs in increasing Y order (see ConvertListOfSlabsToSegments) | |
983 std::sort(rectanglesForEachSlab.begin(), rectanglesForEachSlab.end(), CompareSlabsY); | |
984 | |
985 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, projected.size()); | |
986 #else | |
987 BoostMultiPolygon merged; | |
988 Union(merged, projected); | |
989 | |
990 polygons.resize(merged.size()); | |
991 for (size_t i = 0; i < merged.size(); i++) | |
992 { | |
993 const std::vector<BoostPoint>& outer = merged[i].outer(); | |
994 | |
995 polygons[i].resize(outer.size()); | |
996 for (size_t j = 0; j < outer.size(); j++) | |
997 { | |
998 polygons[i][j] = Point2D(outer[j].x(), outer[j].y()); | |
999 } | |
1000 } | |
1001 #endif | |
1002 | |
1003 #else | |
1004 for (Polygons::iterator polygon = structure.polygons_.begin(); | |
1005 polygon != structure.polygons_.end(); ++polygon) | |
1006 { | |
1007 double x1, y1, x2, y2; | |
1008 if (polygon->Project(x1, y1, x2, y2, slice)) | |
1009 { | |
1010 std::vector<Point2D> p(4); | |
1011 p[0] = std::make_pair(x1, y1); | |
1012 p[1] = std::make_pair(x2, y1); | |
1013 p[2] = std::make_pair(x2, y2); | |
1014 p[3] = std::make_pair(x1, y2); | |
1015 polygons.push_back(p); | |
1016 } | |
1017 } | |
1018 #endif | |
1019 | |
1020 return true; | |
1021 } | |
1022 else | |
1023 { | |
1024 return false; | |
1025 } | |
1026 } | |
1027 | |
1028 | |
1029 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, | |
1030 const CoordinateSystem3D& plane, | |
1031 size_t structureIndex, | |
1032 const Color& color) const | |
1033 { | |
1034 #if USE_BOOST_UNION_FOR_POLYGONS == 1 | |
1035 std::vector< std::vector<Point2D> > polygons; | |
1036 if (ProjectStructure(polygons, structureIndex, plane)) | |
1037 { | |
1038 for (size_t j = 0; j < polygons.size(); j++) | |
1039 { | |
1040 std::vector<ScenePoint2D> chain; | |
1041 chain.reserve(polygons[j].size()); | |
1042 | |
1043 for (size_t k = 0; k < polygons[j].size(); k++) | |
1044 { | |
1045 chain.push_back(ScenePoint2D(polygons[j][k].x, polygons[j][k].y)); | |
1046 } | |
1047 | |
1048 layer.AddChain(chain, true, color.GetRed(), color.GetGreen(), color.GetBlue()); | |
1049 } | |
1050 } | |
1051 | |
1052 #else | |
1053 std::vector< std::pair<Point2D, Point2D> > segments; | |
1054 | |
1055 if (ProjectStructure(segments, structureIndex, plane)) | |
1056 { | |
1057 for (size_t j = 0; j < segments.size(); j++) | |
1058 { | |
1059 std::vector<ScenePoint2D> chain(2); | |
1060 chain[0] = ScenePoint2D(segments[j].first.x, segments[j].first.y); | |
1061 chain[1] = ScenePoint2D(segments[j].second.x, segments[j].second.y); | |
1062 layer.AddChain(chain, false, color.GetRed(), color.GetGreen(), color.GetBlue()); | |
1063 } | |
1064 } | |
1065 #endif | |
1066 } | |
1067 } |