Mercurial > hg > orthanc-stl
comparison Sources/Plugin.cpp @ 35:ee3bc8f7df5b
reorganization
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Thu, 04 Apr 2024 20:37:08 +0200 |
parents | bee2017f3088 |
children | 13698d34e059 |
comparison
equal
deleted
inserted
replaced
34:bee2017f3088 | 35:ee3bc8f7df5b |
---|---|
20 * You should have received a copy of the GNU General Public License | 20 * You should have received a copy of the GNU General Public License |
21 * along with this program. If not, see <http://www.gnu.org/licenses/>. | 21 * along with this program. If not, see <http://www.gnu.org/licenses/>. |
22 **/ | 22 **/ |
23 | 23 |
24 | 24 |
25 #include "StructureSetGeometry.h" | |
26 #include "StructureSet.h" | |
25 #include "STLToolbox.h" | 27 #include "STLToolbox.h" |
26 #include "StructurePolygon.h" | 28 #include "StructurePolygon.h" |
27 #include "VTKToolbox.h" | 29 #include "VTKToolbox.h" |
28 #include "Vector3D.h" | 30 #include "Vector3D.h" |
29 #include "Toolbox.h" | 31 #include "Toolbox.h" |
152 | 154 |
153 | 155 |
154 #include <dcmtk/dcmdata/dcfilefo.h> | 156 #include <dcmtk/dcmdata/dcfilefo.h> |
155 #include <dcmtk/dcmdata/dcsequen.h> | 157 #include <dcmtk/dcmdata/dcsequen.h> |
156 #include <dcmtk/dcmdata/dcuid.h> | 158 #include <dcmtk/dcmdata/dcuid.h> |
157 | |
158 | |
159 class StructureSet : public boost::noncopyable | |
160 { | |
161 private: | |
162 std::vector<StructurePolygon*> polygons_; | |
163 std::string patientId_; | |
164 std::string studyInstanceUid_; | |
165 std::string seriesInstanceUid_; | |
166 std::string sopInstanceUid_; | |
167 bool hasFrameOfReferenceUid_; | |
168 std::string frameOfReferenceUid_; | |
169 | |
170 public: | |
171 explicit StructureSet(Orthanc::ParsedDicomFile& dicom) : | |
172 hasFrameOfReferenceUid_(false) | |
173 { | |
174 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); | |
175 patientId_ = STLToolbox::GetStringValue(dataset, DCM_PatientID); | |
176 studyInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_StudyInstanceUID); | |
177 seriesInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SeriesInstanceUID); | |
178 sopInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SOPInstanceUID); | |
179 | |
180 DcmSequenceOfItems* frame = NULL; | |
181 if (!dataset.findAndGetSequence(DCM_ReferencedFrameOfReferenceSequence, frame).good() || | |
182 frame == NULL) | |
183 { | |
184 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
185 } | |
186 | |
187 if (frame->card() == 1) | |
188 { | |
189 const char* v = NULL; | |
190 if (frame->getItem(0)->findAndGetString(DCM_FrameOfReferenceUID, v).good() && | |
191 v != NULL) | |
192 { | |
193 hasFrameOfReferenceUid_ = true; | |
194 frameOfReferenceUid_.assign(v); | |
195 } | |
196 } | |
197 | |
198 DcmSequenceOfItems* rois = NULL; | |
199 if (!dataset.findAndGetSequence(DCM_ROIContourSequence, rois).good() || | |
200 rois == NULL) | |
201 { | |
202 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
203 } | |
204 | |
205 std::vector<DcmSequenceOfItems*> contours(rois->card()); | |
206 size_t countPolygons = 0; | |
207 | |
208 for (unsigned long i = 0; i < rois->card(); i++) | |
209 { | |
210 DcmSequenceOfItems* contour = NULL; | |
211 if (!rois->getItem(i)->findAndGetSequence(DCM_ContourSequence, contour).good() || | |
212 contour == NULL) | |
213 { | |
214 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
215 } | |
216 else | |
217 { | |
218 contours[i] = contour; | |
219 countPolygons += contour->card(); | |
220 } | |
221 } | |
222 | |
223 polygons_.resize(countPolygons); | |
224 | |
225 size_t pos = 0; | |
226 for (unsigned long i = 0; i < contours.size(); i++) | |
227 { | |
228 for (unsigned long j = 0; j < contours[i]->card(); j++, pos++) | |
229 { | |
230 polygons_[pos] = new StructurePolygon(dicom, i, j); | |
231 } | |
232 } | |
233 | |
234 assert(pos == countPolygons); | |
235 } | |
236 | |
237 ~StructureSet() | |
238 { | |
239 for (size_t i = 0; i < polygons_.size(); i++) | |
240 { | |
241 assert(polygons_[i] != NULL); | |
242 delete polygons_[i]; | |
243 } | |
244 } | |
245 | |
246 const std::string& GetPatientId() const | |
247 { | |
248 return patientId_; | |
249 } | |
250 | |
251 const std::string& GetStudyInstanceUid() const | |
252 { | |
253 return studyInstanceUid_; | |
254 } | |
255 | |
256 const std::string& GetSeriesInstanceUid() const | |
257 { | |
258 return seriesInstanceUid_; | |
259 } | |
260 | |
261 const std::string& GetSopInstanceUid() const | |
262 { | |
263 return sopInstanceUid_; | |
264 } | |
265 | |
266 std::string HashStudy() const | |
267 { | |
268 Orthanc::DicomInstanceHasher hasher(patientId_, studyInstanceUid_, seriesInstanceUid_, sopInstanceUid_); | |
269 return hasher.HashStudy(); | |
270 } | |
271 | |
272 size_t GetPolygonsCount() const | |
273 { | |
274 return polygons_.size(); | |
275 } | |
276 | |
277 const StructurePolygon& GetPolygon(size_t i) const | |
278 { | |
279 if (i >= polygons_.size()) | |
280 { | |
281 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
282 } | |
283 else | |
284 { | |
285 assert(polygons_[i] != NULL); | |
286 return *polygons_[i]; | |
287 } | |
288 } | |
289 | |
290 bool HasFrameOfReferenceUid() const | |
291 { | |
292 return hasFrameOfReferenceUid_; | |
293 } | |
294 | |
295 const std::string& GetFrameOfReferenceUid() const | |
296 { | |
297 if (hasFrameOfReferenceUid_) | |
298 { | |
299 return frameOfReferenceUid_; | |
300 } | |
301 else | |
302 { | |
303 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
304 } | |
305 } | |
306 | |
307 // This static method is faster than constructing the full "StructureSet" object | |
308 static void ListStructuresNames(std::set<std::string>& target, | |
309 Orthanc::ParsedDicomFile& source) | |
310 { | |
311 target.clear(); | |
312 | |
313 DcmSequenceOfItems* sequence = NULL; | |
314 if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() || | |
315 sequence == NULL) | |
316 { | |
317 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
318 } | |
319 | |
320 for (unsigned long i = 0; i < sequence->card(); i++) | |
321 { | |
322 DcmItem* item = sequence->getItem(i); | |
323 if (item == NULL) | |
324 { | |
325 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
326 } | |
327 else | |
328 { | |
329 target.insert(STLToolbox::GetStringValue(*item, DCM_ROIName)); | |
330 } | |
331 } | |
332 } | |
333 }; | |
334 | |
335 | |
336 class StructureSetGeometry : public boost::noncopyable | |
337 { | |
338 private: | |
339 bool strict_; | |
340 Vector3D slicesNormal_; | |
341 double slicesSpacing_; | |
342 double minProjectionAlongNormal_; | |
343 double maxProjectionAlongNormal_; | |
344 size_t slicesCount_; | |
345 | |
346 bool LookupProjectionIndex(size_t& index, | |
347 double z) const | |
348 { | |
349 if (slicesCount_ == 0) | |
350 { | |
351 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
352 } | |
353 | |
354 if (z < minProjectionAlongNormal_ || | |
355 z > maxProjectionAlongNormal_) | |
356 { | |
357 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
358 } | |
359 | |
360 assert(slicesSpacing_ > 0 && | |
361 minProjectionAlongNormal_ < maxProjectionAlongNormal_); | |
362 | |
363 double d = (z - minProjectionAlongNormal_) / slicesSpacing_; | |
364 | |
365 if (STLToolbox::IsNear(d, round(d))) | |
366 { | |
367 if (d < 0.0 || | |
368 d > static_cast<double>(slicesCount_) - 1.0) | |
369 { | |
370 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
371 } | |
372 else | |
373 { | |
374 index = static_cast<size_t>(round(d)); | |
375 return true; | |
376 } | |
377 } | |
378 else | |
379 { | |
380 return false; | |
381 } | |
382 } | |
383 | |
384 public: | |
385 StructureSetGeometry(const StructureSet& structures, | |
386 bool strict) : | |
387 strict_(strict) | |
388 { | |
389 bool isValid = false; | |
390 | |
391 std::vector<double> projections; | |
392 projections.reserve(structures.GetPolygonsCount()); | |
393 | |
394 for (size_t i = 0; i < structures.GetPolygonsCount(); i++) | |
395 { | |
396 Vector3D normal; | |
397 if (structures.GetPolygon(i).IsCoplanar(normal)) | |
398 { | |
399 // Initialize the normal of the whole volume, if need be | |
400 if (!isValid) | |
401 { | |
402 isValid = true; | |
403 slicesNormal_ = normal; | |
404 } | |
405 | |
406 if (Vector3D::AreParallel(normal, slicesNormal_)) | |
407 { | |
408 // This is a valid slice (it is parallel to the normal) | |
409 const Vector3D& point = structures.GetPolygon(i).GetPoint(0); | |
410 projections.push_back(Vector3D::DotProduct(point, slicesNormal_)); | |
411 } | |
412 else | |
413 { | |
414 // RT-STRUCT with non-parallel slices | |
415 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
416 } | |
417 } | |
418 else | |
419 { | |
420 // Ignore slices that are not coplanar | |
421 } | |
422 } | |
423 | |
424 if (projections.empty()) | |
425 { | |
426 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented, | |
427 "Structure set without a valid geometry"); | |
428 } | |
429 | |
430 // Only keep unique projections | |
431 | |
432 std::sort(projections.begin(), projections.end()); | |
433 STLToolbox::RemoveDuplicateValues(projections); | |
434 assert(!projections.empty()); | |
435 | |
436 if (projections.size() == 1) | |
437 { | |
438 // Volume with one single slice | |
439 minProjectionAlongNormal_ = projections[0]; | |
440 maxProjectionAlongNormal_ = projections[0]; | |
441 slicesSpacing_ = 1; // Arbitrary value | |
442 slicesCount_ = 1; | |
443 return; | |
444 } | |
445 | |
446 | |
447 // Compute the most probable spacing between the slices | |
448 | |
449 { | |
450 std::vector<double> spacings; | |
451 spacings.resize(projections.size() - 1); | |
452 | |
453 for (size_t i = 0; i < spacings.size(); i++) | |
454 { | |
455 spacings[i] = projections[i + 1] - projections[i]; | |
456 assert(spacings[i] > 0); | |
457 } | |
458 | |
459 std::sort(spacings.begin(), spacings.end()); | |
460 STLToolbox::RemoveDuplicateValues(spacings); | |
461 | |
462 if (spacings.empty()) | |
463 { | |
464 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
465 } | |
466 | |
467 slicesSpacing_ = spacings[spacings.size() / 10]; // Take the 90% percentile of smallest spacings | |
468 assert(slicesSpacing_ > 0); | |
469 } | |
470 | |
471 | |
472 // Find the projection along the normal with the largest support | |
473 | |
474 bool first = true; | |
475 size_t bestSupport; | |
476 double bestProjection; | |
477 | |
478 std::list<size_t> candidates; | |
479 for (size_t i = 0; i < projections.size(); i++) | |
480 { | |
481 candidates.push_back(i); | |
482 } | |
483 | |
484 while (!candidates.empty()) | |
485 { | |
486 std::list<size_t> next; | |
487 | |
488 size_t countSupport = 0; | |
489 | |
490 std::list<size_t>::const_iterator it = candidates.begin(); | |
491 size_t reference = *it; | |
492 it++; | |
493 | |
494 while (it != candidates.end()) | |
495 { | |
496 double d = (projections[*it] - projections[reference]) / slicesSpacing_; | |
497 if (STLToolbox::IsNear(d, round(d))) | |
498 { | |
499 countSupport ++; | |
500 } | |
501 else | |
502 { | |
503 next.push_back(*it); | |
504 } | |
505 | |
506 it++; | |
507 } | |
508 | |
509 if (first || | |
510 countSupport > bestSupport) | |
511 { | |
512 first = false; | |
513 bestSupport = countSupport; | |
514 bestProjection = projections[reference]; | |
515 } | |
516 | |
517 if (strict && | |
518 !next.empty()) | |
519 { | |
520 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, | |
521 "Structure set with multiple support, which is not allowed in Strict mode"); | |
522 } | |
523 | |
524 candidates.swap(next); | |
525 } | |
526 | |
527 | |
528 // Compute the range of the projections | |
529 | |
530 minProjectionAlongNormal_ = bestProjection; | |
531 maxProjectionAlongNormal_ = bestProjection; | |
532 | |
533 for (size_t i = 0; i < projections.size(); i++) | |
534 { | |
535 double d = (projections[i] - bestProjection) / slicesSpacing_; | |
536 if (STLToolbox::IsNear(d, round(d))) | |
537 { | |
538 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]); | |
539 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]); | |
540 } | |
541 } | |
542 | |
543 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; | |
544 if (STLToolbox::IsNear(d, round(d))) | |
545 { | |
546 slicesCount_ = static_cast<size_t>(round(d)) + 1; | |
547 } | |
548 else | |
549 { | |
550 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
551 } | |
552 | |
553 | |
554 // Sanity check | |
555 | |
556 size_t a, b; | |
557 if (!LookupProjectionIndex(a, minProjectionAlongNormal_) || | |
558 !LookupProjectionIndex(b, maxProjectionAlongNormal_) || | |
559 a != 0 || | |
560 b + 1 != slicesCount_) | |
561 { | |
562 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
563 } | |
564 } | |
565 | |
566 const Vector3D& GetSlicesNormal() const | |
567 { | |
568 return slicesNormal_; | |
569 } | |
570 | |
571 double GetSlicesSpacing() const | |
572 { | |
573 return slicesSpacing_; | |
574 } | |
575 | |
576 double GetMinProjectionAlongNormal() const | |
577 { | |
578 return minProjectionAlongNormal_; | |
579 } | |
580 | |
581 double GetMaxProjectionAlongNormal() const | |
582 { | |
583 return maxProjectionAlongNormal_; | |
584 } | |
585 | |
586 bool ProjectAlongNormal(double& z, | |
587 const StructurePolygon& polygon) const | |
588 { | |
589 Vector3D normal; | |
590 if (polygon.IsCoplanar(normal) && | |
591 Vector3D::AreParallel(normal, slicesNormal_)) | |
592 { | |
593 z = Vector3D::DotProduct(polygon.GetPoint(0), slicesNormal_); | |
594 return true; | |
595 } | |
596 else | |
597 { | |
598 return false; | |
599 } | |
600 } | |
601 | |
602 size_t GetSlicesCount() const | |
603 { | |
604 return slicesCount_; | |
605 } | |
606 | |
607 bool LookupSliceIndex(size_t& slice, | |
608 const StructurePolygon& polygon) const | |
609 { | |
610 double z; | |
611 return (ProjectAlongNormal(z, polygon) && | |
612 LookupProjectionIndex(slice, z)); | |
613 } | |
614 }; | |
615 | |
616 | |
617 | 159 |
618 | 160 |
619 class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller | 161 class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller |
620 { | 162 { |
621 private: | 163 private: |