Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet2.cpp @ 998:38b6bb0bdd72
added a new set of classes that correctly handle non-convex polygons (not
used yet because of limitations in coordinates computing): DicomStructure2,
DicomStructureSet2, DicomStructurePolygon2, DicomStructureSetSlicer2. Too
many shortcuts have been taken when computing the actual position.
author | Benjamin Golinvaux <bgo@osimis.io> |
---|---|
date | Fri, 20 Sep 2019 11:58:00 +0200 |
parents | 1f74bc3459ba |
children | 4f28d9459e31 |
comparison
equal
deleted
inserted
replaced
995:9893fa8cd7a6 | 998:38b6bb0bdd72 |
---|---|
16 * | 16 * |
17 * You should have received a copy of the GNU Affero General Public License | 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/>. | 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. |
19 **/ | 19 **/ |
20 | 20 |
21 | |
22 #include "DicomStructureSet2.h" | 21 #include "DicomStructureSet2.h" |
23 | 22 |
24 #include "../Toolbox/GeometryToolbox.h" | 23 #include "../Toolbox/LinearAlgebra.h" |
24 #include "../StoneException.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 <Core/Toolbox.h> |
29 #include <Core/DicomFormat/DicomTag.h> | |
30 | |
29 #include <Plugins/Samples/Common/FullOrthancDataset.h> | 31 #include <Plugins/Samples/Common/FullOrthancDataset.h> |
30 #include <Plugins/Samples/Common/DicomDatasetReader.h> | 32 #include <Plugins/Samples/Common/DicomDatasetReader.h> |
31 | |
32 #include <limits> | |
33 #include <stdio.h> | |
34 #include <boost/geometry.hpp> | |
35 #include <boost/geometry/geometries/point_xy.hpp> | |
36 #include <boost/geometry/geometries/polygon.hpp> | |
37 #include <boost/geometry/multi/geometries/multi_polygon.hpp> | |
38 | |
39 typedef boost::geometry::model::d2::point_xy<double> BoostPoint; | |
40 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; | |
41 typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; | |
42 | |
43 | |
44 static void Union(BoostMultiPolygon& output, | |
45 std::vector<BoostPolygon>& input) | |
46 { | |
47 for (size_t i = 0; i < input.size(); i++) | |
48 { | |
49 boost::geometry::correct(input[i]); | |
50 } | |
51 | |
52 if (input.size() == 0) | |
53 { | |
54 output.clear(); | |
55 } | |
56 else if (input.size() == 1) | |
57 { | |
58 output.resize(1); | |
59 output[0] = input[0]; | |
60 } | |
61 else | |
62 { | |
63 boost::geometry::union_(input[0], input[1], output); | |
64 | |
65 for (size_t i = 0; i < input.size(); i++) | |
66 { | |
67 BoostMultiPolygon tmp; | |
68 boost::geometry::union_(output, input[i], tmp); | |
69 output = tmp; | |
70 } | |
71 } | |
72 } | |
73 | |
74 | |
75 static BoostPolygon CreateRectangle(float x1, float y1, | |
76 float x2, float y2) | |
77 { | |
78 BoostPolygon r; | |
79 boost::geometry::append(r, BoostPoint(x1, y1)); | |
80 boost::geometry::append(r, BoostPoint(x1, y2)); | |
81 boost::geometry::append(r, BoostPoint(x2, y2)); | |
82 boost::geometry::append(r, BoostPoint(x2, y1)); | |
83 return r; | |
84 } | |
85 | |
86 | |
87 | 33 |
88 namespace OrthancStone | 34 namespace OrthancStone |
89 { | 35 { |
90 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); | 36 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); |
91 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); | 37 static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); |
98 static const OrthancPlugins::DicomTag DICOM_TAG_ROI_NAME(0x3006, 0x0026); | 44 static const OrthancPlugins::DicomTag DICOM_TAG_ROI_NAME(0x3006, 0x0026); |
99 static const OrthancPlugins::DicomTag DICOM_TAG_RT_ROI_INTERPRETED_TYPE(0x3006, 0x00a4); | 45 static const OrthancPlugins::DicomTag DICOM_TAG_RT_ROI_INTERPRETED_TYPE(0x3006, 0x00a4); |
100 static const OrthancPlugins::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080); | 46 static const OrthancPlugins::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080); |
101 static const OrthancPlugins::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020); | 47 static const OrthancPlugins::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020); |
102 | 48 |
103 | 49 static inline uint8_t ConvertAndClipToByte(double v) |
104 static uint8_t ConvertColor(double v) | |
105 { | 50 { |
106 if (v < 0) | 51 if (v < 0) |
107 { | 52 { |
108 return 0; | 53 return 0; |
109 } | 54 } |
115 { | 60 { |
116 return static_cast<uint8_t>(v); | 61 return static_cast<uint8_t>(v); |
117 } | 62 } |
118 } | 63 } |
119 | 64 |
120 | 65 static bool ReadDicomToVector(Vector& target, |
121 static bool ParseVector(Vector& target, | 66 const OrthancPlugins::IDicomDataset& dataset, |
122 const OrthancPlugins::IDicomDataset& dataset, | 67 const OrthancPlugins::DicomPath& tag) |
123 const OrthancPlugins::DicomPath& tag) | |
124 { | 68 { |
125 std::string value; | 69 std::string value; |
126 return (dataset.GetStringValue(value, tag) && | 70 return (dataset.GetStringValue(value, tag) && |
127 LinearAlgebra::ParseVector(target, value)); | 71 LinearAlgebra::ParseVector(target, value)); |
128 } | 72 } |
129 | 73 |
130 void DicomStructureSet2::Polygon::CheckPointIsOnSlice(const Vector& v) const | 74 std::ostream& operator<<(std::ostream& s, const OrthancPlugins::DicomPath& dicomPath) |
131 { | 75 { |
132 if (hasSlice_) | 76 std::stringstream tmp; |
133 { | 77 for (size_t i = 0; i < dicomPath.GetPrefixLength(); ++i) |
134 double magnitude = | 78 { |
135 GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()); | 79 OrthancPlugins::DicomTag tag = dicomPath.GetPrefixTag(i); |
136 if(!LinearAlgebra::IsNear( | |
137 magnitude, | |
138 projectionAlongNormal_, | |
139 sliceThickness_ / 2.0 /* in mm */ )) | |
140 { | |
141 LOG(ERROR) << "This RT-STRUCT contains a point that is off the " | |
142 << "slice of its instance | " | |
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 DicomStructureSet2::Polygon::IsPointOnSliceIfAny(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 } | |
177 | |
178 void DicomStructureSet2::Polygon::AddPoint(const Vector& v) | |
179 { | |
180 #if 1 | |
181 // BGO 2019-09-03 | |
182 if (IsPointOnSliceIfAny(v)) | |
183 { | |
184 points_.push_back(v); | |
185 } | |
186 #else | |
187 CheckPoint(v); | |
188 points_.push_back(v); | |
189 #endif | |
190 } | |
191 | |
192 | |
193 bool DicomStructureSet2::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices) | |
194 { | |
195 if (hasSlice_) | |
196 { | |
197 return true; | |
198 } | |
199 else | |
200 { | |
201 ReferencedSlices::const_iterator it = slices.find(sopInstanceUid_); | |
202 | 80 |
203 if (it == slices.end()) | 81 // We use this other object to be able to use GetMainTagsName |
204 { | 82 // and Format |
205 return false; | 83 Orthanc::DicomTag tag2(tag.GetGroup(), tag.GetElement()); |
206 } | 84 size_t index = dicomPath.GetPrefixIndex(i); |
207 else | 85 tmp << tag2.GetMainTagsName() << " (" << tag2.Format() << ") [" << index << "] / "; |
208 { | 86 } |
209 const CoordinateSystem3D& geometry = it->second.geometry_; | 87 const OrthancPlugins::DicomTag& tag = dicomPath.GetFinalTag(); |
210 | 88 Orthanc::DicomTag tag2(tag.GetGroup(), tag.GetElement()); |
211 hasSlice_ = true; | 89 tmp << tag2.GetMainTagsName() << " (" << tag2.Format() << ")"; |
212 geometry_ = geometry; | 90 s << tmp.str(); |
213 projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal()); | 91 return s; |
214 sliceThickness_ = it->second.thickness_; | 92 } |
215 | 93 |
216 extent_.Reset(); | 94 |
217 | 95 void DicomStructureSet2::SetContents(const OrthancPlugins::FullOrthancDataset& tags) |
218 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) | 96 { |
219 { | 97 FillStructuresFromDataset(tags); |
220 if (IsPointOnSliceIfAny(*it)) | 98 ComputeDependentProperties(); |
221 { | 99 } |
222 double x, y; | 100 |
223 geometry.ProjectPoint(x, y, *it); | 101 void DicomStructureSet2::ComputeDependentProperties() |
224 extent_.AddPoint(x, y); | 102 { |
225 } | 103 for (size_t i = 0; i < structures_.size(); ++i) |
226 } | 104 { |
227 return true; | 105 structures_[i].ComputeDependentProperties(); |
228 } | 106 } |
229 } | 107 } |
230 } | 108 |
231 | 109 void DicomStructureSet2::FillStructuresFromDataset(const OrthancPlugins::FullOrthancDataset& tags) |
232 bool DicomStructureSet2::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const | |
233 { | |
234 bool isOpposite = false; | |
235 | |
236 if (points_.empty() || | |
237 !hasSlice_ || | |
238 !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) | |
239 { | |
240 return false; | |
241 } | |
242 | |
243 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); | |
244 | |
245 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, | |
246 sliceThickness_ / 2.0)); | |
247 } | |
248 | |
249 bool DicomStructureSet2::Polygon::Project(double& x1, | |
250 double& y1, | |
251 double& x2, | |
252 double& y2, | |
253 const CoordinateSystem3D& slice) const | |
254 { | |
255 // TODO: optimize this method using a sweep-line algorithm for polygons | |
256 | |
257 if (!hasSlice_ || | |
258 points_.size() <= 1) | |
259 { | |
260 return false; | |
261 } | |
262 | |
263 double x, y; | |
264 geometry_.ProjectPoint(x, y, slice.GetOrigin()); | |
265 | |
266 bool isOpposite; | |
267 if (GeometryToolbox::IsParallelOrOpposite | |
268 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | |
269 { | |
270 if (y < extent_.GetY1() || | |
271 y > extent_.GetY2()) | |
272 { | |
273 // The polygon does not intersect the input slice | |
274 return false; | |
275 } | |
276 | |
277 bool isFirst = true; | |
278 double xmin = std::numeric_limits<double>::infinity(); | |
279 double xmax = -std::numeric_limits<double>::infinity(); | |
280 | |
281 double prevX, prevY; | |
282 geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); | |
283 | |
284 for (size_t i = 0; i < points_.size(); i++) | |
285 { | |
286 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py | |
287 double curX, curY; | |
288 geometry_.ProjectPoint(curX, curY, points_[i]); | |
289 | |
290 if ((prevY < y && curY > y) || | |
291 (prevY > y && curY < y)) | |
292 { | |
293 double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); | |
294 xmin = std::min(xmin, p); | |
295 xmax = std::max(xmax, p); | |
296 isFirst = false; | |
297 } | |
298 | |
299 prevX = curX; | |
300 prevY = curY; | |
301 } | |
302 | |
303 if (isFirst) | |
304 { | |
305 return false; | |
306 } | |
307 else | |
308 { | |
309 Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + | |
310 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
311 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - | |
312 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
313 | |
314 slice.ProjectPoint(x1, y1, p1); | |
315 slice.ProjectPoint(x2, y2, p2); | |
316 return true; | |
317 } | |
318 } | |
319 else if (GeometryToolbox::IsParallelOrOpposite | |
320 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | |
321 { | |
322 if (x < extent_.GetX1() || | |
323 x > extent_.GetX2()) | |
324 { | |
325 return false; | |
326 } | |
327 | |
328 bool isFirst = true; | |
329 double ymin = std::numeric_limits<double>::infinity(); | |
330 double ymax = -std::numeric_limits<double>::infinity(); | |
331 | |
332 double prevX, prevY; | |
333 geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); | |
334 | |
335 for (size_t i = 0; i < points_.size(); i++) | |
336 { | |
337 // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py | |
338 double curX, curY; | |
339 geometry_.ProjectPoint(curX, curY, points_[i]); | |
340 | |
341 if ((prevX < x && curX > x) || | |
342 (prevX > x && curX < x)) | |
343 { | |
344 double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); | |
345 ymin = std::min(ymin, p); | |
346 ymax = std::max(ymax, p); | |
347 isFirst = false; | |
348 } | |
349 | |
350 prevX = curX; | |
351 prevY = curY; | |
352 } | |
353 | |
354 if (isFirst) | |
355 { | |
356 return false; | |
357 } | |
358 else | |
359 { | |
360 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + | |
361 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
362 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - | |
363 sliceThickness_ / 2.0 * geometry_.GetNormal()); | |
364 | |
365 slice.ProjectPoint(x1, y1, p1); | |
366 slice.ProjectPoint(x2, y2, p2); | |
367 | |
368 // TODO WHY THIS??? | |
369 y1 = -y1; | |
370 y2 = -y2; | |
371 | |
372 return true; | |
373 } | |
374 } | |
375 else | |
376 { | |
377 // Should not happen | |
378 return false; | |
379 } | |
380 } | |
381 | |
382 | |
383 const DicomStructureSet2::Structure& DicomStructureSet2::GetStructure(size_t index) const | |
384 { | |
385 if (index >= structures_.size()) | |
386 { | |
387 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
388 } | |
389 | |
390 return structures_[index]; | |
391 } | |
392 | |
393 | |
394 DicomStructureSet2::Structure& DicomStructureSet2::GetStructure(size_t index) | |
395 { | |
396 if (index >= structures_.size()) | |
397 { | |
398 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
399 } | |
400 | |
401 return structures_[index]; | |
402 } | |
403 | |
404 DicomStructureSet2::DicomStructureSet2(const OrthancPlugins::FullOrthancDataset& tags) | |
405 { | 110 { |
406 OrthancPlugins::DicomDatasetReader reader(tags); | 111 OrthancPlugins::DicomDatasetReader reader(tags); |
407 | 112 |
408 size_t count, tmp; | 113 // a few sanity checks |
114 size_t count = 0, tmp = 0; | |
115 | |
116 // DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE (0x3006, 0x0080); | |
117 // DICOM_TAG_ROI_CONTOUR_SEQUENCE (0x3006, 0x0039); | |
118 // DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE (0x3006, 0x0020); | |
409 if (!tags.GetSequenceSize(count, DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE) || | 119 if (!tags.GetSequenceSize(count, DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE) || |
410 !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) || | 120 !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) || |
411 tmp != count || | 121 tmp != count || |
412 !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) || | 122 !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) || |
413 tmp != count) | 123 tmp != count) |
414 { | 124 { |
415 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 125 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
416 } | 126 } |
417 | 127 |
128 // let's now parse the structures stored in the dicom file | |
129 // DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE (0x3006, 0x0080) | |
130 // DICOM_TAG_RT_ROI_INTERPRETED_TYPE (0x3006, 0x00a4) | |
131 // DICOM_TAG_ROI_DISPLAY_COLOR (0x3006, 0x002a) | |
132 // DICOM_TAG_ROI_NAME (0x3006, 0x0026) | |
418 structures_.resize(count); | 133 structures_.resize(count); |
419 for (size_t i = 0; i < count; i++) | 134 for (size_t i = 0; i < count; i++) |
420 { | 135 { |
136 // (0x3006, 0x0080)[i]/(0x3006, 0x00a4) | |
421 structures_[i].interpretation_ = reader.GetStringValue | 137 structures_[i].interpretation_ = reader.GetStringValue |
422 (OrthancPlugins::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, | 138 (OrthancPlugins::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, |
423 DICOM_TAG_RT_ROI_INTERPRETED_TYPE), | 139 DICOM_TAG_RT_ROI_INTERPRETED_TYPE), |
424 "No interpretation"); | 140 "No interpretation"); |
425 | 141 |
142 // (0x3006, 0x0020)[i]/(0x3006, 0x0026) | |
426 structures_[i].name_ = reader.GetStringValue | 143 structures_[i].name_ = reader.GetStringValue |
427 (OrthancPlugins::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, | 144 (OrthancPlugins::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, |
428 DICOM_TAG_ROI_NAME), | 145 DICOM_TAG_ROI_NAME), |
429 "No name"); | 146 "No name"); |
430 | 147 |
431 Vector color; | 148 Vector color; |
432 if (ParseVector(color, tags, OrthancPlugins::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 149 // (0x3006, 0x0039)[i]/(0x3006, 0x002a) |
433 DICOM_TAG_ROI_DISPLAY_COLOR)) && | 150 if (ReadDicomToVector(color, tags, OrthancPlugins::DicomPath( |
434 color.size() == 3) | 151 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, DICOM_TAG_ROI_DISPLAY_COLOR)) |
435 { | 152 && color.size() == 3) |
436 structures_[i].red_ = ConvertColor(color[0]); | 153 { |
437 structures_[i].green_ = ConvertColor(color[1]); | 154 structures_[i].red_ = ConvertAndClipToByte(color[0]); |
438 structures_[i].blue_ = ConvertColor(color[2]); | 155 structures_[i].green_ = ConvertAndClipToByte(color[1]); |
156 structures_[i].blue_ = ConvertAndClipToByte(color[2]); | |
439 } | 157 } |
440 else | 158 else |
441 { | 159 { |
442 structures_[i].red_ = 255; | 160 structures_[i].red_ = 255; |
443 structures_[i].green_ = 0; | 161 structures_[i].green_ = 0; |
444 structures_[i].blue_ = 0; | 162 structures_[i].blue_ = 0; |
445 } | 163 } |
446 | 164 |
447 size_t countSlices; | 165 size_t countSlices; |
448 if (!tags.GetSequenceSize(countSlices, OrthancPlugins::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 166 // DICOM_TAG_ROI_CONTOUR_SEQUENCE (0x3006, 0x0039); |
449 DICOM_TAG_CONTOUR_SEQUENCE))) | 167 // DICOM_TAG_CONTOUR_SEQUENCE (0x3006, 0x0040); |
450 { | 168 if (!tags.GetSequenceSize(countSlices, OrthancPlugins::DicomPath( |
169 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, DICOM_TAG_CONTOUR_SEQUENCE))) | |
170 { | |
171 LOG(WARNING) << "DicomStructureSet2::SetContents | structure \"" << structures_[i].name_ << "\" has no slices!"; | |
451 countSlices = 0; | 172 countSlices = 0; |
452 } | 173 } |
453 | 174 |
454 LOG(INFO) << "New RT structure: \"" << structures_[i].name_ | 175 LOG(INFO) << "New RT structure: \"" << structures_[i].name_ |
455 << "\" with interpretation \"" << structures_[i].interpretation_ | 176 << "\" with interpretation \"" << structures_[i].interpretation_ |
456 << "\" containing " << countSlices << " slices (color: " | 177 << "\" containing " << countSlices << " slices (color: " |
457 << static_cast<int>(structures_[i].red_) << "," | 178 << static_cast<int>(structures_[i].red_) << "," |
458 << static_cast<int>(structures_[i].green_) << "," | 179 << static_cast<int>(structures_[i].green_) << "," |
459 << static_cast<int>(structures_[i].blue_) << ")"; | 180 << static_cast<int>(structures_[i].blue_) << ")"; |
460 | 181 |
461 // These temporary variables avoid allocating many vectors in the loop below | 182 // These temporary variables avoid allocating many vectors in the loop below |
462 OrthancPlugins::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
463 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
464 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); | |
465 | |
466 OrthancPlugins::DicomPath geometricTypePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
467 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
468 DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); | |
469 | 183 |
470 OrthancPlugins::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 184 // (0x3006, 0x0039)[i]/(0x3006, 0x0040)[0]/(0x3006, 0x0046) |
471 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 185 OrthancPlugins::DicomPath countPointsPath( |
472 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); | 186 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, |
187 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
188 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); | |
189 | |
190 OrthancPlugins::DicomPath geometricTypePath( | |
191 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
192 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
193 DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); | |
194 | |
195 OrthancPlugins::DicomPath imageSequencePath( | |
196 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
197 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
198 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); | |
473 | 199 |
474 // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155) | 200 // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155) |
475 OrthancPlugins::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 201 OrthancPlugins::DicomPath referencedInstancePath( |
476 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 202 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, |
477 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, | 203 DICOM_TAG_CONTOUR_SEQUENCE, 0, |
478 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); | 204 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, |
479 | 205 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); |
480 OrthancPlugins::DicomPath contourDataPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 206 |
481 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 207 OrthancPlugins::DicomPath contourDataPath( |
482 DICOM_TAG_CONTOUR_DATA); | 208 DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, |
209 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
210 DICOM_TAG_CONTOUR_DATA); | |
483 | 211 |
484 for (size_t j = 0; j < countSlices; j++) | 212 for (size_t j = 0; j < countSlices; j++) |
485 { | 213 { |
486 unsigned int countPoints; | 214 unsigned int countPoints = 0; |
487 | 215 |
488 countPointsPath.SetPrefixIndex(1, j); | 216 countPointsPath.SetPrefixIndex(1, j); |
489 if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath)) | 217 if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath)) |
490 { | 218 { |
219 LOG(ERROR) << "Dicom path " << countPointsPath << " is not valid (should contain an unsigned integer)"; | |
491 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 220 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
492 } | 221 } |
493 | 222 |
494 //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices"; | 223 //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices"; |
495 | 224 |
496 geometricTypePath.SetPrefixIndex(1, j); | 225 geometricTypePath.SetPrefixIndex(1, j); |
497 std::string type = reader.GetMandatoryStringValue(geometricTypePath); | 226 std::string type = reader.GetMandatoryStringValue(geometricTypePath); |
498 if (type != "CLOSED_PLANAR") | 227 if (type != "CLOSED_PLANAR") |
499 { | 228 { |
229 // TODO: support points!! | |
500 LOG(WARNING) << "Ignoring contour with geometry type: " << type; | 230 LOG(WARNING) << "Ignoring contour with geometry type: " << type; |
501 continue; | 231 continue; |
502 } | 232 } |
503 | 233 |
504 size_t size; | 234 size_t size = 0; |
505 | 235 |
506 imageSequencePath.SetPrefixIndex(1, j); | 236 imageSequencePath.SetPrefixIndex(1, j); |
507 if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1) | 237 if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1) |
508 { | 238 { |
509 LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry."; | 239 LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry."; |
510 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | 240 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); |
511 } | 241 } |
512 | 242 |
513 referencedInstancePath.SetPrefixIndex(1, j); | 243 referencedInstancePath.SetPrefixIndex(1, j); |
514 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); | 244 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); |
515 | 245 |
516 contourDataPath.SetPrefixIndex(1, j); | 246 contourDataPath.SetPrefixIndex(1, j); |
517 std::string slicesData = reader.GetMandatoryStringValue(contourDataPath); | 247 std::string slicesData = reader.GetMandatoryStringValue(contourDataPath); |
518 | 248 |
519 Vector points; | 249 Vector points; |
520 if (!LinearAlgebra::ParseVector(points, slicesData) || | 250 if (!LinearAlgebra::ParseVector(points, slicesData) || |
521 points.size() != 3 * countPoints) | 251 points.size() != 3 * countPoints) |
522 { | 252 { |
523 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 253 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
524 } | 254 } |
525 | 255 |
526 // seen in real world | 256 // seen in real world |
527 if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") | 257 if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") |
528 { | 258 { |
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)"; | 259 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)"; |
530 } | 260 } |
531 | 261 |
532 Polygon polygon(sopInstanceUid); | 262 DicomStructurePolygon2 polygon(sopInstanceUid,type); |
533 polygon.Reserve(countPoints); | 263 polygon.Reserve(countPoints); |
534 | 264 |
535 for (size_t k = 0; k < countPoints; k++) | 265 for (size_t k = 0; k < countPoints; k++) |
536 { | 266 { |
537 Vector v(3); | 267 Vector v(3); |
538 v[0] = points[3 * k]; | 268 v[0] = points[3 * k]; |
539 v[1] = points[3 * k + 1]; | 269 v[1] = points[3 * k + 1]; |
540 v[2] = points[3 * k + 2]; | 270 v[2] = points[3 * k + 2]; |
541 polygon.AddPoint(v); | 271 polygon.AddPoint(v); |
542 } | 272 } |
543 | 273 structures_[i].AddPolygon(polygon); |
544 structures_[i].polygons_.push_back(polygon); | 274 } |
545 } | 275 } |
546 } | 276 } |
547 } | 277 |
548 | 278 |
549 | 279 void DicomStructureSet2::Clear() |
550 Vector DicomStructureSet2::GetStructureCenter(size_t index) const | 280 { |
551 { | 281 structures_.clear(); |
552 const Structure& structure = GetStructure(index); | 282 } |
553 | 283 |
554 Vector center; | |
555 LinearAlgebra::AssignVector(center, 0, 0, 0); | |
556 if (structure.polygons_.empty()) | |
557 { | |
558 return center; | |
559 } | |
560 | |
561 double n = static_cast<double>(structure.polygons_.size()); | |
562 | |
563 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
564 polygon != structure.polygons_.end(); ++polygon) | |
565 { | |
566 if (!polygon->GetPoints().empty()) | |
567 { | |
568 center += polygon->GetPoints().front() / n; | |
569 } | |
570 } | |
571 | |
572 return center; | |
573 } | |
574 | |
575 | |
576 const std::string& DicomStructureSet2::GetStructureName(size_t index) const | |
577 { | |
578 return GetStructure(index).name_; | |
579 } | |
580 | |
581 | |
582 const std::string& DicomStructureSet2::GetStructureInterpretation(size_t index) const | |
583 { | |
584 return GetStructure(index).interpretation_; | |
585 } | |
586 | |
587 | |
588 Color DicomStructureSet2::GetStructureColor(size_t index) const | |
589 { | |
590 const Structure& s = GetStructure(index); | |
591 return Color(s.red_, s.green_, s.blue_); | |
592 } | |
593 | |
594 | |
595 void DicomStructureSet2::GetStructureColor(uint8_t& red, | |
596 uint8_t& green, | |
597 uint8_t& blue, | |
598 size_t index) const | |
599 { | |
600 const Structure& s = GetStructure(index); | |
601 red = s.red_; | |
602 green = s.green_; | |
603 blue = s.blue_; | |
604 } | |
605 | |
606 | |
607 void DicomStructureSet2::GetReferencedInstances(std::set<std::string>& instances) | |
608 { | |
609 for (Structures::const_iterator structure = structures_.begin(); | |
610 structure != structures_.end(); ++structure) | |
611 { | |
612 for (Polygons::const_iterator polygon = structure->polygons_.begin(); | |
613 polygon != structure->polygons_.end(); ++polygon) | |
614 { | |
615 instances.insert(polygon->GetSopInstanceUid()); | |
616 } | |
617 } | |
618 } | |
619 | |
620 | |
621 void DicomStructureSet2::AddReferencedSlice(const std::string& sopInstanceUid, | |
622 const std::string& seriesInstanceUid, | |
623 const CoordinateSystem3D& geometry, | |
624 double thickness) | |
625 { | |
626 if (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()) | |
627 { | |
628 // This geometry is already known | |
629 LOG(ERROR) << "DicomStructureSet2::AddReferencedSlice(): (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()). sopInstanceUid = " << sopInstanceUid; | |
630 | |
631 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
632 } | |
633 else | |
634 { | |
635 if (thickness < 0) | |
636 { | |
637 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
638 } | |
639 | |
640 if (!referencedSlices_.empty()) | |
641 { | |
642 const ReferencedSlice& reference = referencedSlices_.begin()->second; | |
643 | |
644 if (reference.seriesInstanceUid_ != seriesInstanceUid) | |
645 { | |
646 LOG(ERROR) << "This RT-STRUCT refers to several different series"; | |
647 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
648 } | |
649 | |
650 if (!GeometryToolbox::IsParallel(reference.geometry_.GetNormal(), geometry.GetNormal())) | |
651 { | |
652 LOG(ERROR) << "The slices in this RT-STRUCT are not parallel"; | |
653 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
654 } | |
655 } | |
656 | |
657 referencedSlices_[sopInstanceUid] = ReferencedSlice(seriesInstanceUid, geometry, thickness); | |
658 | |
659 for (Structures::iterator structure = structures_.begin(); | |
660 structure != structures_.end(); ++structure) | |
661 { | |
662 for (Polygons::iterator polygon = structure->polygons_.begin(); | |
663 polygon != structure->polygons_.end(); ++polygon) | |
664 { | |
665 polygon->UpdateReferencedSlice(referencedSlices_); | |
666 } | |
667 } | |
668 } | |
669 } | |
670 | |
671 | |
672 void DicomStructureSet2::AddReferencedSlice(const Orthanc::DicomMap& dataset) | |
673 { | |
674 CoordinateSystem3D slice(dataset); | |
675 | |
676 double thickness = 1; // 1 mm by default | |
677 | |
678 std::string s; | |
679 Vector v; | |
680 if (dataset.LookupStringValue(s, Orthanc::DICOM_TAG_SLICE_THICKNESS, false) && | |
681 LinearAlgebra::ParseVector(v, s) && | |
682 v.size() > 0) | |
683 { | |
684 thickness = v[0]; | |
685 } | |
686 | |
687 std::string instance, series; | |
688 if (dataset.LookupStringValue(instance, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false) && | |
689 dataset.LookupStringValue(series, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false)) | |
690 { | |
691 AddReferencedSlice(instance, series, slice, thickness); | |
692 } | |
693 else | |
694 { | |
695 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
696 } | |
697 } | |
698 | |
699 | |
700 void DicomStructureSet2::CheckReferencedSlices() | |
701 { | |
702 for (Structures::iterator structure = structures_.begin(); | |
703 structure != structures_.end(); ++structure) | |
704 { | |
705 for (Polygons::iterator polygon = structure->polygons_.begin(); | |
706 polygon != structure->polygons_.end(); ++polygon) | |
707 { | |
708 if (!polygon->UpdateReferencedSlice(referencedSlices_)) | |
709 { | |
710 std::string sopInstanceUid = polygon->GetSopInstanceUid(); | |
711 if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") | |
712 { | |
713 LOG(ERROR) << "DicomStructureSet2::CheckReferencedSlices(): " | |
714 << " missing information about referenced instance " | |
715 << "(sopInstanceUid is empty!)"; | |
716 } | |
717 else | |
718 { | |
719 LOG(ERROR) << "DicomStructureSet2::CheckReferencedSlices(): " | |
720 << " missing information about referenced instance " | |
721 << "(sopInstanceUid = " << sopInstanceUid << ")"; | |
722 } | |
723 //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
724 } | |
725 } | |
726 } | |
727 } | |
728 | |
729 | |
730 Vector DicomStructureSet2::GetNormal() const | |
731 { | |
732 if (referencedSlices_.empty()) | |
733 { | |
734 Vector v; | |
735 LinearAlgebra::AssignVector(v, 0, 0, 1); | |
736 return v; | |
737 } | |
738 else | |
739 { | |
740 return referencedSlices_.begin()->second.geometry_.GetNormal(); | |
741 } | |
742 } | |
743 | |
744 | |
745 bool DicomStructureSet2::ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons, | |
746 const Structure& structure, | |
747 const CoordinateSystem3D& slice) const | |
748 { | |
749 polygons.clear(); | |
750 | |
751 Vector normal = GetNormal(); | |
752 | |
753 bool isOpposite; | |
754 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) | |
755 { | |
756 // This is an axial projection | |
757 | |
758 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
759 polygon != structure.polygons_.end(); ++polygon) | |
760 { | |
761 if (polygon->IsOnSlice(slice)) | |
762 { | |
763 polygons.push_back(std::vector<PolygonPoint2D>()); | |
764 | |
765 for (Points::const_iterator p = polygon->GetPoints().begin(); | |
766 p != polygon->GetPoints().end(); ++p) | |
767 { | |
768 double x, y; | |
769 slice.ProjectPoint(x, y, *p); | |
770 polygons.back().push_back(std::make_pair(x, y)); | |
771 } | |
772 } | |
773 } | |
774 | |
775 return true; | |
776 } | |
777 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | |
778 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | |
779 { | |
780 #if 1 | |
781 // Sagittal or coronal projection | |
782 std::vector<BoostPolygon> projected; | |
783 | |
784 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | |
785 polygon != structure.polygons_.end(); ++polygon) | |
786 { | |
787 double x1, y1, x2, y2; | |
788 if (polygon->Project(x1, y1, x2, y2, slice)) | |
789 { | |
790 projected.push_back(CreateRectangle( | |
791 static_cast<float>(x1), | |
792 static_cast<float>(y1), | |
793 static_cast<float>(x2), | |
794 static_cast<float>(y2))); | |
795 } | |
796 } | |
797 | |
798 BoostMultiPolygon merged; | |
799 Union(merged, projected); | |
800 | |
801 polygons.resize(merged.size()); | |
802 for (size_t i = 0; i < merged.size(); i++) | |
803 { | |
804 const std::vector<BoostPoint>& outer = merged[i].outer(); | |
805 | |
806 polygons[i].resize(outer.size()); | |
807 for (size_t j = 0; j < outer.size(); j++) | |
808 { | |
809 polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); | |
810 } | |
811 } | |
812 #else | |
813 for (Polygons::iterator polygon = structure.polygons_.begin(); | |
814 polygon != structure.polygons_.end(); ++polygon) | |
815 { | |
816 double x1, y1, x2, y2; | |
817 if (polygon->Project(x1, y1, x2, y2, slice)) | |
818 { | |
819 std::vector<PolygonPoint2D> p(4); | |
820 p[0] = std::make_pair(x1, y1); | |
821 p[1] = std::make_pair(x2, y1); | |
822 p[2] = std::make_pair(x2, y2); | |
823 p[3] = std::make_pair(x1, y2); | |
824 polygons.push_back(p); | |
825 } | |
826 } | |
827 #endif | |
828 | |
829 return true; | |
830 } | |
831 else | |
832 { | |
833 return false; | |
834 } | |
835 } | |
836 } | 284 } |