Mercurial > hg > orthanc-stl
comparison Sources/Plugin.cpp @ 1:0f03a8a0bd6f
encoding of RT-STRUCT as STL
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Mon, 17 Jul 2023 18:54:31 +0200 |
parents | 4e889a8e8be2 |
children | 2bdb9acb7dcf |
comparison
equal
deleted
inserted
replaced
0:4e889a8e8be2 | 1:0f03a8a0bd6f |
---|---|
22 **/ | 22 **/ |
23 | 23 |
24 | 24 |
25 #include "../Resources/Orthanc/Plugins/OrthancPluginCppWrapper.h" | 25 #include "../Resources/Orthanc/Plugins/OrthancPluginCppWrapper.h" |
26 | 26 |
27 #include <ChunkedBuffer.h> | |
28 #include <DicomParsing/ParsedDicomFile.h> | |
29 #include <Images/ImageProcessing.h> | |
27 #include <Logging.h> | 30 #include <Logging.h> |
31 #include <OrthancFramework.h> | |
32 #include <SerializationToolbox.h> | |
28 #include <SystemToolbox.h> | 33 #include <SystemToolbox.h> |
29 | 34 |
30 #include <EmbeddedResources.h> | 35 #include <EmbeddedResources.h> |
36 | |
37 #include <vtkImageData.h> | |
38 #include <vtkImageResize.h> | |
39 #include <vtkMarchingCubes.h> | |
40 #include <vtkNew.h> | |
41 #include <vtkPolyData.h> | |
42 #include <vtkTriangle.h> | |
43 #include <vtkSmoothPolyDataFilter.h> | |
44 #include <vtkPolyDataNormals.h> | |
45 #include <vtkImageConstantPad.h> | |
31 | 46 |
32 #include <boost/thread/shared_mutex.hpp> | 47 #include <boost/thread/shared_mutex.hpp> |
33 | 48 |
34 // Forward declaration | 49 // Forward declaration |
35 void ReadStaticAsset(std::string& target, | 50 void ReadStaticAsset(std::string& target, |
119 OrthancPluginAnswerBuffer(OrthancPlugins::GetGlobalContext(), output, s.c_str(), s.size(), Orthanc::EnumerationToString(Orthanc::MimeType_JavaScript)); | 134 OrthancPluginAnswerBuffer(OrthancPlugins::GetGlobalContext(), output, s.c_str(), s.size(), Orthanc::EnumerationToString(Orthanc::MimeType_JavaScript)); |
120 } | 135 } |
121 else | 136 else |
122 { | 137 { |
123 cache_.Answer(output, file); | 138 cache_.Answer(output, file); |
139 } | |
140 } | |
141 | |
142 | |
143 | |
144 | |
145 #include <dcmtk/dcmdata/dcdeftag.h> | |
146 #include <dcmtk/dcmdata/dcfilefo.h> | |
147 #include <dcmtk/dcmdata/dcitem.h> | |
148 #include <dcmtk/dcmdata/dcsequen.h> | |
149 | |
150 class Extent2D : public boost::noncopyable | |
151 { | |
152 private: | |
153 bool isEmpty_; | |
154 double x1_; | |
155 double y1_; | |
156 double x2_; | |
157 double y2_; | |
158 | |
159 void CheckNotEmpty() const | |
160 { | |
161 if (isEmpty_) | |
162 { | |
163 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
164 } | |
165 } | |
166 | |
167 public: | |
168 Extent2D() : | |
169 isEmpty_(true), | |
170 x1_(0), | |
171 y1_(0), | |
172 x2_(0), | |
173 y2_(0) | |
174 { | |
175 } | |
176 | |
177 bool IsEmpty() const | |
178 { | |
179 return isEmpty_; | |
180 } | |
181 | |
182 double GetMinX() const | |
183 { | |
184 CheckNotEmpty(); | |
185 return x1_; | |
186 } | |
187 | |
188 double GetMaxX() const | |
189 { | |
190 CheckNotEmpty(); | |
191 return x2_; | |
192 } | |
193 | |
194 double GetMinY() const | |
195 { | |
196 CheckNotEmpty(); | |
197 return y1_; | |
198 } | |
199 | |
200 double GetMaxY() const | |
201 { | |
202 CheckNotEmpty(); | |
203 return y2_; | |
204 } | |
205 | |
206 double GetWidth() const | |
207 { | |
208 CheckNotEmpty(); | |
209 return x2_ - x1_; | |
210 } | |
211 | |
212 double GetHeight() const | |
213 { | |
214 CheckNotEmpty(); | |
215 return y2_ - y1_; | |
216 } | |
217 | |
218 void Add(double x, | |
219 double y) | |
220 { | |
221 if (isEmpty_) | |
222 { | |
223 x1_ = x2_ = x; | |
224 y1_ = y2_ = y; | |
225 isEmpty_ = false; | |
226 } | |
227 else | |
228 { | |
229 x1_ = std::min(x1_, x); | |
230 x2_ = std::max(x2_, x); | |
231 y1_ = std::min(y1_, y); | |
232 y2_ = std::max(y2_, y); | |
233 } | |
234 } | |
235 }; | |
236 | |
237 static void GetStringValue(std::string& value, | |
238 DcmItem& item, | |
239 const DcmTagKey& key) | |
240 { | |
241 const char* s = NULL; | |
242 if (!item.findAndGetString(key, s).good() || | |
243 s == NULL) | |
244 { | |
245 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
246 } | |
247 else | |
248 { | |
249 value.assign(s); | |
250 } | |
251 } | |
252 | |
253 | |
254 static void ListStructuresNames(std::set<std::string>& target, | |
255 Orthanc::ParsedDicomFile& source) | |
256 { | |
257 target.clear(); | |
258 | |
259 DcmSequenceOfItems* sequence = NULL; | |
260 if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() || | |
261 sequence == NULL) | |
262 { | |
263 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
264 } | |
265 | |
266 for (unsigned long i = 0; i < sequence->card(); i++) | |
267 { | |
268 DcmItem* item = sequence->getItem(i); | |
269 if (item == NULL) | |
270 { | |
271 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
272 } | |
273 else | |
274 { | |
275 std::string value; | |
276 GetStringValue(value, *item, DCM_ROIName); | |
277 target.insert(value); | |
278 } | |
279 } | |
280 } | |
281 | |
282 | |
283 static bool IsNear(double a, | |
284 double b) | |
285 { | |
286 return std::abs(a - b) < 10.0 * std::numeric_limits<double>::epsilon(); | |
287 } | |
288 | |
289 | |
290 class Vector3D | |
291 { | |
292 private: | |
293 double x_; | |
294 double y_; | |
295 double z_; | |
296 | |
297 public: | |
298 Vector3D() : | |
299 x_(0), | |
300 y_(0), | |
301 z_(0) | |
302 { | |
303 } | |
304 | |
305 Vector3D(double x, | |
306 double y, | |
307 double z) : | |
308 x_(x), | |
309 y_(y), | |
310 z_(z) | |
311 { | |
312 } | |
313 | |
314 Vector3D(const Vector3D& from, | |
315 const Vector3D& to) : | |
316 x_(to.x_ - from.x_), | |
317 y_(to.y_ - from.y_), | |
318 z_(to.z_ - from.z_) | |
319 { | |
320 } | |
321 | |
322 double GetX() const | |
323 { | |
324 return x_; | |
325 } | |
326 | |
327 double GetY() const | |
328 { | |
329 return y_; | |
330 } | |
331 | |
332 double GetZ() const | |
333 { | |
334 return z_; | |
335 } | |
336 | |
337 double ComputeNorm() const | |
338 { | |
339 return sqrt(x_ * x_ + y_ * y_ + z_ * z_); | |
340 } | |
341 | |
342 void Normalize() | |
343 { | |
344 double norm = ComputeNorm(); | |
345 if (!IsNear(norm, 0)) | |
346 { | |
347 x_ /= norm; | |
348 y_ /= norm; | |
349 z_ /= norm; | |
350 } | |
351 } | |
352 | |
353 static Vector3D CrossProduct(const Vector3D& u, | |
354 const Vector3D& v) | |
355 { | |
356 return Vector3D(u.GetY() * v.GetZ() - u.GetZ() * v.GetY(), | |
357 u.GetZ() * v.GetX() - u.GetX() * v.GetZ(), | |
358 u.GetX() * v.GetY() - u.GetY() * v.GetX()); | |
359 } | |
360 | |
361 static double DotProduct(const Vector3D& a, | |
362 const Vector3D& b) | |
363 { | |
364 return a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ(); | |
365 } | |
366 }; | |
367 | |
368 | |
369 | |
370 static bool MyParseDouble(double& value, | |
371 const std::string& s) | |
372 { | |
373 #if 1 | |
374 char* end = NULL; | |
375 value = strtod(s.c_str(), &end); | |
376 return (end == s.c_str() + s.size()); | |
377 #else | |
378 return Orthanc::SerializationToolbox::ParseDouble(value, s); | |
379 #endif | |
380 } | |
381 | |
382 | |
383 class StructurePolygon : public boost::noncopyable | |
384 { | |
385 private: | |
386 std::string roiName_; | |
387 std::string referencedSopInstanceUid_; | |
388 uint8_t red_; | |
389 uint8_t green_; | |
390 uint8_t blue_; | |
391 std::vector<Vector3D> points_; | |
392 | |
393 public: | |
394 StructurePolygon(Orthanc::ParsedDicomFile& dicom, | |
395 unsigned long roiIndex, | |
396 unsigned long contourIndex) | |
397 { | |
398 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); | |
399 | |
400 DcmItem* structure = NULL; | |
401 DcmItem* roi = NULL; | |
402 DcmItem* contour = NULL; | |
403 DcmSequenceOfItems* referenced = NULL; | |
404 | |
405 if (!dataset.findAndGetSequenceItem(DCM_StructureSetROISequence, structure, roiIndex).good() || | |
406 structure == NULL || | |
407 !dataset.findAndGetSequenceItem(DCM_ROIContourSequence, roi, roiIndex).good() || | |
408 roi == NULL || | |
409 !roi->findAndGetSequenceItem(DCM_ContourSequence, contour, contourIndex).good() || | |
410 contour == NULL || | |
411 !contour->findAndGetSequence(DCM_ContourImageSequence, referenced).good() || | |
412 referenced == NULL || | |
413 referenced->card() != 1) | |
414 { | |
415 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
416 } | |
417 | |
418 GetStringValue(roiName_, *structure, DCM_ROIName); | |
419 GetStringValue(referencedSopInstanceUid_, *referenced->getItem(0), DCM_ReferencedSOPInstanceUID); | |
420 | |
421 std::string s; | |
422 GetStringValue(s, *contour, DCM_ContourGeometricType); | |
423 if (s != "CLOSED_PLANAR") | |
424 { | |
425 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
426 } | |
427 | |
428 { | |
429 std::string color; | |
430 GetStringValue(color, *roi, DCM_ROIDisplayColor); | |
431 | |
432 std::vector<std::string> tokens; | |
433 Orthanc::Toolbox::TokenizeString(tokens, color, '\\'); | |
434 | |
435 uint32_t r, g, b; | |
436 if (tokens.size() != 3 || | |
437 !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(r, tokens[0]) || | |
438 !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(g, tokens[1]) || | |
439 !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(b, tokens[2]) || | |
440 r > 255 || | |
441 g > 255 || | |
442 b > 255) | |
443 { | |
444 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
445 } | |
446 | |
447 red_ = r; | |
448 green_ = g; | |
449 blue_ = b; | |
450 } | |
451 | |
452 { | |
453 GetStringValue(s, *contour, DCM_ContourData); | |
454 | |
455 std::vector<std::string> tokens; | |
456 Orthanc::Toolbox::TokenizeString(tokens, s, '\\'); | |
457 | |
458 GetStringValue(s, *contour, DCM_NumberOfContourPoints); | |
459 | |
460 uint32_t countPoints; | |
461 if (!Orthanc::SerializationToolbox::ParseUnsignedInteger32(countPoints, s) || | |
462 tokens.size() != 3 * countPoints) | |
463 { | |
464 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
465 } | |
466 | |
467 points_.reserve(countPoints); | |
468 | |
469 for (size_t i = 0; i < tokens.size(); i += 3) | |
470 { | |
471 double x, y, z; | |
472 if (!MyParseDouble(x, tokens[i]) || | |
473 !MyParseDouble(y, tokens[i + 1]) || | |
474 !MyParseDouble(z, tokens[i + 2])) | |
475 { | |
476 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
477 } | |
478 | |
479 points_.push_back(Vector3D(x, y, z)); | |
480 } | |
481 | |
482 assert(points_.size() == countPoints); | |
483 } | |
484 } | |
485 | |
486 const std::string& GetRoiName() const | |
487 { | |
488 return roiName_; | |
489 } | |
490 | |
491 const std::string& GetReferencedSopInstanceUid() const | |
492 { | |
493 return referencedSopInstanceUid_; | |
494 } | |
495 | |
496 size_t GetPointsCount() const | |
497 { | |
498 return points_.size(); | |
499 } | |
500 | |
501 const Vector3D& GetPoint(size_t i) const | |
502 { | |
503 if (i >= points_.size()) | |
504 { | |
505 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
506 } | |
507 else | |
508 { | |
509 return points_[i]; | |
510 } | |
511 } | |
512 | |
513 bool IsCoplanar(Vector3D& normal) const | |
514 { | |
515 if (points_.size() < 3) | |
516 { | |
517 return false; | |
518 } | |
519 | |
520 bool hasNormal = false; | |
521 | |
522 for (size_t i = 0; i < points_.size(); i++) | |
523 { | |
524 normal = Vector3D::CrossProduct(Vector3D(points_[1], points_[0]), | |
525 Vector3D(points_[2], points_[0])); | |
526 if (!IsNear(normal.ComputeNorm(), 0)) | |
527 { | |
528 normal.Normalize(); | |
529 hasNormal = true; | |
530 } | |
531 } | |
532 | |
533 if (!hasNormal) | |
534 { | |
535 return false; | |
536 } | |
537 | |
538 double a = Vector3D::DotProduct(points_[0], normal); | |
539 | |
540 for (size_t i = 1; i < points_.size(); i++) | |
541 { | |
542 double b = Vector3D::DotProduct(points_[i], normal); | |
543 if (!IsNear(a, b)) | |
544 { | |
545 return false; | |
546 } | |
547 } | |
548 | |
549 return true; | |
550 } | |
551 | |
552 void Add(Extent2D& extent, | |
553 const Vector3D& axisX, | |
554 const Vector3D& axisY) const | |
555 { | |
556 assert(IsNear(1, axisX.ComputeNorm())); | |
557 assert(IsNear(1, axisY.ComputeNorm())); | |
558 | |
559 for (size_t i = 0; i < points_.size(); i++) | |
560 { | |
561 extent.Add(Vector3D::DotProduct(axisX, points_[i]), | |
562 Vector3D::DotProduct(axisY, points_[i])); | |
563 } | |
564 } | |
565 }; | |
566 | |
567 | |
568 | |
569 struct IsNearPredicate | |
570 { | |
571 bool operator() (const double& a, | |
572 const double& b) | |
573 { | |
574 return IsNear(a, b); | |
575 } | |
576 }; | |
577 | |
578 static void RemoveDuplicateValues(std::vector<double>& v) | |
579 { | |
580 IsNearPredicate predicate; | |
581 std::vector<double>::iterator last = std::unique(v.begin(), v.end(), predicate); | |
582 v.erase(last, v.end()); | |
583 } | |
584 | |
585 | |
586 class StructureSet : public boost::noncopyable | |
587 { | |
588 private: | |
589 std::vector<StructurePolygon*> polygons_; | |
590 bool hasGeometry_; | |
591 Vector3D slicesNormal_; | |
592 double slicesSpacing_; | |
593 double minProjectionAlongNormal_; | |
594 double maxProjectionAlongNormal_; | |
595 std::string patientId_; | |
596 std::string studyInstanceUid_; | |
597 std::string seriesInstanceUid_; | |
598 std::string sopInstanceUid_; | |
599 | |
600 void ComputeGeometry() | |
601 { | |
602 std::list<double> positionsList; | |
603 hasGeometry_ = false; | |
604 | |
605 for (size_t i = 0; i < polygons_.size(); i++) | |
606 { | |
607 assert(polygons_[i] != NULL); | |
608 | |
609 Vector3D n; | |
610 if (polygons_[i]->IsCoplanar(n)) | |
611 { | |
612 const Vector3D& point = polygons_[i]->GetPoint(0); | |
613 double z = Vector3D::DotProduct(point, n); | |
614 | |
615 if (!hasGeometry_) | |
616 { | |
617 hasGeometry_ = true; | |
618 slicesNormal_ = n; | |
619 minProjectionAlongNormal_ = z; | |
620 maxProjectionAlongNormal_ = z; | |
621 } | |
622 else if (!IsNear(std::abs(Vector3D::DotProduct(n, slicesNormal_)), 1)) | |
623 { | |
624 hasGeometry_ = false; | |
625 | |
626 // RT-STRUCT with non-parallel slices | |
627 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
628 } | |
629 else | |
630 { | |
631 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, z); | |
632 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, z); | |
633 } | |
634 | |
635 positionsList.push_back(Vector3D::DotProduct(n, point)); | |
636 } | |
637 } | |
638 | |
639 if (hasGeometry_) | |
640 { | |
641 std::vector<double> positions(positionsList.begin(), positionsList.end()); | |
642 assert(!positions.empty()); | |
643 | |
644 std::sort(positions.begin(), positions.end()); | |
645 RemoveDuplicateValues(positions); | |
646 assert(!positions.empty()); | |
647 | |
648 if (positions.size() == 1) | |
649 { | |
650 hasGeometry_ = false; | |
651 return; | |
652 } | |
653 | |
654 std::vector<double> offsets; | |
655 offsets.resize(positions.size() - 1); | |
656 | |
657 for (size_t i = 0; i < offsets.size(); i++) | |
658 { | |
659 offsets[i] = positions[i + 1] - positions[i]; | |
660 assert(offsets[i] > 0); | |
661 } | |
662 | |
663 std::sort(offsets.begin(), offsets.end()); | |
664 RemoveDuplicateValues(offsets); | |
665 | |
666 slicesSpacing_ = offsets[0]; | |
667 | |
668 for (size_t i = 1; i < offsets.size(); i++) | |
669 { | |
670 double d = offsets[i] / slicesSpacing_; | |
671 if (!IsNear(d, round(d))) | |
672 { | |
673 // Irregular spacing between the slices | |
674 hasGeometry_ = false; | |
675 break; | |
676 } | |
677 } | |
678 } | |
679 } | |
680 | |
681 void CheckHasGeometry() const | |
682 { | |
683 if (!hasGeometry_) | |
684 { | |
685 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
686 } | |
687 } | |
688 | |
689 public: | |
690 explicit StructureSet(Orthanc::ParsedDicomFile& dicom) : | |
691 hasGeometry_(false), | |
692 slicesSpacing_(0), | |
693 minProjectionAlongNormal_(0), | |
694 maxProjectionAlongNormal_(0) | |
695 { | |
696 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); | |
697 GetStringValue(patientId_, dataset, DCM_PatientID); | |
698 GetStringValue(studyInstanceUid_, dataset, DCM_StudyInstanceUID); | |
699 GetStringValue(seriesInstanceUid_, dataset, DCM_SeriesInstanceUID); | |
700 GetStringValue(sopInstanceUid_, dataset, DCM_SOPInstanceUID); | |
701 | |
702 DcmSequenceOfItems* rois = NULL; | |
703 if (!dataset.findAndGetSequence(DCM_ROIContourSequence, rois).good() || | |
704 rois == NULL) | |
705 { | |
706 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
707 } | |
708 | |
709 std::vector<DcmSequenceOfItems*> contours(rois->card()); | |
710 size_t countPolygons = 0; | |
711 | |
712 for (unsigned long i = 0; i < rois->card(); i++) | |
713 { | |
714 DcmSequenceOfItems* contour = NULL; | |
715 if (!rois->getItem(i)->findAndGetSequence(DCM_ContourSequence, contour).good() || | |
716 contour == NULL) | |
717 { | |
718 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
719 } | |
720 else | |
721 { | |
722 contours[i] = contour; | |
723 countPolygons += contour->card(); | |
724 } | |
725 } | |
726 | |
727 polygons_.resize(countPolygons); | |
728 | |
729 size_t pos = 0; | |
730 for (unsigned long i = 0; i < contours.size(); i++) | |
731 { | |
732 for (unsigned long j = 0; j < contours[i]->card(); j++, pos++) | |
733 { | |
734 polygons_[pos] = new StructurePolygon(dicom, i, j); | |
735 } | |
736 } | |
737 | |
738 assert(pos == countPolygons); | |
739 | |
740 ComputeGeometry(); | |
741 } | |
742 | |
743 ~StructureSet() | |
744 { | |
745 for (size_t i = 0; i < polygons_.size(); i++) | |
746 { | |
747 assert(polygons_[i] != NULL); | |
748 delete polygons_[i]; | |
749 } | |
750 } | |
751 | |
752 const std::string& GetPatient() const | |
753 { | |
754 return patientId_; | |
755 } | |
756 | |
757 const std::string& GetStudyInstanceUid() const | |
758 { | |
759 return studyInstanceUid_; | |
760 } | |
761 | |
762 const std::string& GetSeriesInstanceUid() const | |
763 { | |
764 return seriesInstanceUid_; | |
765 } | |
766 | |
767 const std::string& GetSopInstanceUid() const | |
768 { | |
769 return sopInstanceUid_; | |
770 } | |
771 | |
772 std::string HashStudy() const | |
773 { | |
774 std::string s; | |
775 Orthanc::Toolbox::ComputeSHA1(s, patientId_ + "|" + studyInstanceUid_); | |
776 return s; | |
777 } | |
778 | |
779 size_t GetPolygonsCount() const | |
780 { | |
781 return polygons_.size(); | |
782 } | |
783 | |
784 const StructurePolygon& GetPolygon(size_t i) const | |
785 { | |
786 if (i >= polygons_.size()) | |
787 { | |
788 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
789 } | |
790 else | |
791 { | |
792 assert(polygons_[i] != NULL); | |
793 return *polygons_[i]; | |
794 } | |
795 } | |
796 | |
797 bool HasGeometry() const | |
798 { | |
799 return hasGeometry_; | |
800 } | |
801 | |
802 const Vector3D& GetSlicesNormal() const | |
803 { | |
804 CheckHasGeometry(); | |
805 return slicesNormal_; | |
806 } | |
807 | |
808 double GetSlicesSpacing() const | |
809 { | |
810 CheckHasGeometry(); | |
811 return slicesSpacing_; | |
812 } | |
813 | |
814 double GetMinProjectionAlongNormal() const | |
815 { | |
816 CheckHasGeometry(); | |
817 return minProjectionAlongNormal_; | |
818 } | |
819 | |
820 double GetMaxProjectionAlongNormal() const | |
821 { | |
822 CheckHasGeometry(); | |
823 return maxProjectionAlongNormal_; | |
824 } | |
825 | |
826 double ProjectAlongNormal(const StructurePolygon& polygon) const | |
827 { | |
828 CheckHasGeometry(); | |
829 return Vector3D::DotProduct(slicesNormal_, polygon.GetPoint(0)); | |
830 } | |
831 | |
832 size_t GetSlicesCount() const | |
833 { | |
834 CheckHasGeometry(); | |
835 double c = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; | |
836 assert(c >= 0); | |
837 | |
838 if (!IsNear(c, round(c))) | |
839 { | |
840 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
841 } | |
842 else | |
843 { | |
844 return static_cast<size_t>(round(c)) + 1; | |
845 } | |
846 } | |
847 | |
848 bool LookupSliceIndex(size_t& slice, | |
849 const StructurePolygon& polygon) const | |
850 { | |
851 CheckHasGeometry(); | |
852 | |
853 double z = ProjectAlongNormal(polygon); | |
854 | |
855 if (z < minProjectionAlongNormal_ || | |
856 z > maxProjectionAlongNormal_) | |
857 { | |
858 return false; | |
859 } | |
860 else | |
861 { | |
862 double c = (z - minProjectionAlongNormal_) / slicesSpacing_; | |
863 | |
864 if (IsNear(c, round(c))) | |
865 { | |
866 slice = static_cast<size_t>(round(c)); | |
867 return true; | |
868 } | |
869 else | |
870 { | |
871 return false; | |
872 } | |
873 } | |
874 } | |
875 | |
876 bool LookupReferencedSopInstanceUid(std::string& sopInstanceUid) const | |
877 { | |
878 if (HasGeometry()) | |
879 { | |
880 for (size_t i = 0; i < polygons_.size(); i++) | |
881 { | |
882 assert(polygons_[i] != NULL); | |
883 | |
884 Vector3D n; | |
885 if (polygons_[i]->IsCoplanar(n) && | |
886 Vector3D::DotProduct(n, slicesNormal_)) | |
887 { | |
888 sopInstanceUid = polygons_[i]->GetReferencedSopInstanceUid(); | |
889 return true; | |
890 } | |
891 } | |
892 } | |
893 | |
894 return false; | |
895 } | |
896 }; | |
897 | |
898 | |
899 | |
900 | |
901 class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller | |
902 { | |
903 private: | |
904 Orthanc::ImageAccessor& target_; | |
905 | |
906 public: | |
907 XorFiller(Orthanc::ImageAccessor& target) : | |
908 target_(target) | |
909 { | |
910 } | |
911 | |
912 virtual void Fill(int y, | |
913 int x1, | |
914 int x2) ORTHANC_OVERRIDE | |
915 { | |
916 assert(x1 <= x2); | |
917 | |
918 if (y >= 0 && | |
919 y < static_cast<int>(target_.GetHeight())) | |
920 { | |
921 x1 = std::max(x1, 0); | |
922 x2 = std::min(x2, static_cast<int>(target_.GetWidth()) - 1); | |
923 | |
924 uint8_t* p = reinterpret_cast<uint8_t*>(target_.GetRow(y)) + x1; | |
925 | |
926 for (int i = x1; i <= x2; i++, p++) | |
927 { | |
928 *p = (*p ^ 0xff); | |
929 } | |
930 } | |
931 } | |
932 }; | |
933 | |
934 | |
935 static void EncodeSTL(std::string& target /* out */, | |
936 vtkPolyData& mesh /* in */) | |
937 { | |
938 // TODO - Conversion to little endian on big endian | |
939 | |
940 Orthanc::ChunkedBuffer buffer; | |
941 | |
942 uint8_t header[80]; | |
943 memset(header, 0, sizeof(header)); | |
944 buffer.AddChunk(header, sizeof(header)); | |
945 | |
946 uint32_t n = mesh.GetNumberOfCells(); | |
947 buffer.AddChunk(&n, sizeof(n)); | |
948 | |
949 for (vtkIdType i = 0; i < mesh.GetNumberOfCells(); i++) | |
950 { | |
951 vtkCell* cell = mesh.GetCell(i); | |
952 vtkTriangle* triangle = dynamic_cast<vtkTriangle*>(cell); | |
953 | |
954 double p0[3]; | |
955 double p1[3]; | |
956 double p2[3]; | |
957 triangle->GetPoints()->GetPoint(0, p0); | |
958 triangle->GetPoints()->GetPoint(1, p1); | |
959 triangle->GetPoints()->GetPoint(2, p2); | |
960 | |
961 double normal[3]; | |
962 vtkTriangle::ComputeNormal(p0, p1, p2, normal); | |
963 | |
964 float d[4 * 3] = { | |
965 static_cast<float>(normal[0]), static_cast<float>(normal[1]), static_cast<float>(normal[2]), | |
966 static_cast<float>(p0[0]), static_cast<float>(p0[1]), static_cast<float>(p0[2]), | |
967 static_cast<float>(p1[0]), static_cast<float>(p1[1]), static_cast<float>(p1[2]), | |
968 static_cast<float>(p2[0]), static_cast<float>(p2[1]), static_cast<float>(p2[2]) }; | |
969 buffer.AddChunk(d, sizeof(d)); | |
970 | |
971 uint16_t a = 0; | |
972 buffer.AddChunk(&a, sizeof(a)); | |
973 } | |
974 | |
975 buffer.Flatten(target); | |
976 } | |
977 | |
978 | |
979 bool EncodeStructureSetMesh(std::string& stl, | |
980 const StructureSet& structureSet, | |
981 const std::set<std::string>& roiNames, | |
982 unsigned int resolution, | |
983 bool smooth) | |
984 { | |
985 if (!structureSet.HasGeometry()) | |
986 { | |
987 return false; | |
988 } | |
989 | |
990 if (resolution < 1) | |
991 { | |
992 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
993 } | |
994 | |
995 if (!IsNear(1, structureSet.GetSlicesNormal().ComputeNorm())) | |
996 { | |
997 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
998 } | |
999 | |
1000 // TODO - Axes could be retrieved from the referenced CT volume | |
1001 Vector3D axisX(1, 0, 0); | |
1002 Vector3D axisY = Vector3D::CrossProduct(structureSet.GetSlicesNormal(), axisX); | |
1003 | |
1004 if (!IsNear(1, axisX.ComputeNorm()) || | |
1005 !IsNear(1, axisY.ComputeNorm())) | |
1006 { | |
1007 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
1008 } | |
1009 | |
1010 Extent2D extent; | |
1011 for (size_t i = 0; i < structureSet.GetPolygonsCount(); i++) | |
1012 { | |
1013 structureSet.GetPolygon(i).Add(extent, axisX, axisY); | |
1014 } | |
1015 | |
1016 const int depth = structureSet.GetSlicesCount(); | |
1017 | |
1018 vtkNew<vtkImageData> volume; | |
1019 volume->SetDimensions(resolution, resolution, depth); | |
1020 volume->AllocateScalars(VTK_UNSIGNED_CHAR, 1); | |
1021 | |
1022 assert(sizeof(unsigned char) == 1); | |
1023 memset(volume->GetScalarPointer(), 0, resolution * resolution * depth); | |
1024 | |
1025 for (size_t i = 0; i < structureSet.GetPolygonsCount(); i++) | |
1026 { | |
1027 const StructurePolygon& polygon = structureSet.GetPolygon(i); | |
1028 if (roiNames.find(polygon.GetRoiName()) == roiNames.end()) | |
1029 { | |
1030 // This polygon doesn't correspond to a ROI of interest | |
1031 continue; | |
1032 } | |
1033 | |
1034 size_t j; | |
1035 if (!structureSet.LookupSliceIndex(j, polygon)) | |
1036 { | |
1037 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
1038 } | |
1039 | |
1040 std::vector<Orthanc::ImageProcessing::ImagePoint> points; | |
1041 points.reserve(polygon.GetPointsCount()); | |
1042 for (size_t j = 0; j < polygon.GetPointsCount(); j++) | |
1043 { | |
1044 const Vector3D& point = polygon.GetPoint(j); | |
1045 double x = (Vector3D::DotProduct(point, axisX) - extent.GetMinX()) / extent.GetWidth() * static_cast<double>(resolution); | |
1046 double y = (Vector3D::DotProduct(point, axisY) - extent.GetMinY()) / extent.GetHeight() * static_cast<double>(resolution); | |
1047 points.push_back(Orthanc::ImageProcessing::ImagePoint(static_cast<int32_t>(std::floor(x)), | |
1048 static_cast<int32_t>(std::floor(y)))); | |
1049 } | |
1050 | |
1051 Orthanc::ImageAccessor slice; | |
1052 slice.AssignWritable(Orthanc::PixelFormat_Grayscale8, resolution, resolution, resolution /* pitch */, | |
1053 reinterpret_cast<uint8_t*>(volume->GetScalarPointer()) + j * resolution * resolution); | |
1054 | |
1055 XorFiller filler(slice); | |
1056 Orthanc::ImageProcessing::FillPolygon(filler, points); | |
1057 } | |
1058 | |
1059 vtkNew<vtkImageResize> resize; | |
1060 resize->SetOutputDimensions(resolution, resolution, resolution); | |
1061 resize->SetInputData(volume.Get()); | |
1062 resize->Update(); | |
1063 | |
1064 resize->GetOutput()->SetSpacing( | |
1065 extent.GetWidth() / static_cast<double>(resolution), | |
1066 extent.GetHeight() / static_cast<double>(resolution), | |
1067 (structureSet.GetMaxProjectionAlongNormal() - structureSet.GetMinProjectionAlongNormal()) / static_cast<double>(resolution)); | |
1068 | |
1069 // TODO | |
1070 // resize->GetOutput()->SetOrigin() | |
1071 | |
1072 vtkNew<vtkImageConstantPad> padding; | |
1073 padding->SetConstant(0); | |
1074 padding->SetOutputNumberOfScalarComponents(1); | |
1075 padding->SetOutputWholeExtent(-1, resolution, -1, resolution, -1, resolution); | |
1076 padding->SetInputData(resize->GetOutput()); | |
1077 padding->Update(); | |
1078 | |
1079 vtkNew<vtkMarchingCubes> surface; | |
1080 surface->SetInputData(padding->GetOutput()); | |
1081 surface->ComputeNormalsOn(); | |
1082 surface->SetValue(0, 128 /*isoValue*/); | |
1083 surface->Update(); | |
1084 | |
1085 if (smooth) | |
1086 { | |
1087 vtkNew<vtkSmoothPolyDataFilter> smoothFilter; | |
1088 // Apply volume smoothing | |
1089 // https://examples.vtk.org/site/Cxx/PolyData/SmoothPolyDataFilter/ | |
1090 smoothFilter->SetInputConnection(surface->GetOutputPort()); | |
1091 smoothFilter->SetNumberOfIterations(15); | |
1092 smoothFilter->SetRelaxationFactor(0.1); | |
1093 smoothFilter->FeatureEdgeSmoothingOff(); | |
1094 smoothFilter->BoundarySmoothingOn(); | |
1095 smoothFilter->Update(); | |
1096 | |
1097 vtkNew<vtkPolyDataNormals> normalGenerator; | |
1098 normalGenerator->SetInputConnection(smoothFilter->GetOutputPort()); | |
1099 normalGenerator->ComputePointNormalsOn(); | |
1100 normalGenerator->ComputeCellNormalsOn(); | |
1101 normalGenerator->Update(); | |
1102 | |
1103 EncodeSTL(stl, *normalGenerator->GetOutput()); | |
1104 } | |
1105 else | |
1106 { | |
1107 EncodeSTL(stl, *surface->GetOutput()); | |
1108 } | |
1109 | |
1110 return true; | |
1111 } | |
1112 | |
1113 | |
1114 void ListStructures(OrthancPluginRestOutput* output, | |
1115 const char* url, | |
1116 const OrthancPluginHttpRequest* request) | |
1117 { | |
1118 const std::string instanceId(request->groups[0]); | |
1119 | |
1120 if (request->method != OrthancPluginHttpMethod_Get) | |
1121 { | |
1122 OrthancPluginSendMethodNotAllowed(OrthancPlugins::GetGlobalContext(), output, "GET"); | |
1123 return; | |
1124 } | |
1125 | |
1126 std::string dicom; | |
1127 if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false)) | |
1128 { | |
1129 throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource); | |
1130 } | |
1131 else | |
1132 { | |
1133 Orthanc::ParsedDicomFile parsed(dicom); | |
1134 | |
1135 std::set<std::string> names; | |
1136 ListStructuresNames(names, parsed); | |
1137 | |
1138 Json::Value answer = Json::arrayValue; | |
1139 | |
1140 for (std::set<std::string>::const_iterator it = names.begin(); it != names.end(); ++it) | |
1141 { | |
1142 answer.append(*it); | |
1143 } | |
1144 | |
1145 std::string s = answer.toStyledString(); | |
1146 OrthancPluginAnswerBuffer(OrthancPlugins::GetGlobalContext(), output, s.c_str(), s.size(), Orthanc::MIME_JSON); | |
1147 } | |
1148 } | |
1149 | |
1150 | |
1151 void Encode(OrthancPluginRestOutput* output, | |
1152 const char* url, | |
1153 const OrthancPluginHttpRequest* request) | |
1154 { | |
1155 static const char* const KEY_INSTANCE = "Instance"; | |
1156 static const char* const KEY_RESOLUTION = "Resolution"; | |
1157 static const char* const KEY_ROI_NAMES = "RoiNames"; | |
1158 static const char* const KEY_SMOOTH = "Smooth"; | |
1159 static const char* const KEY_TAGS = "Tags"; | |
1160 | |
1161 if (request->method != OrthancPluginHttpMethod_Post) | |
1162 { | |
1163 OrthancPluginSendMethodNotAllowed(OrthancPlugins::GetGlobalContext(), output, "POST"); | |
1164 return; | |
1165 } | |
1166 | |
1167 Json::Value body; | |
1168 if (!Orthanc::Toolbox::ReadJson(body, request->body, request->bodySize)) | |
1169 { | |
1170 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest); | |
1171 } | |
1172 | |
1173 const std::string instanceId = Orthanc::SerializationToolbox::ReadString(body, KEY_INSTANCE); | |
1174 const bool smooth = (body.isMember(KEY_SMOOTH) ? | |
1175 Orthanc::SerializationToolbox::ReadBoolean(body, KEY_SMOOTH) : | |
1176 true /* smooth by default */); | |
1177 const unsigned int resolution = (body.isMember(KEY_RESOLUTION) ? | |
1178 Orthanc::SerializationToolbox::ReadUnsignedInteger(body, KEY_RESOLUTION) : | |
1179 256 /* default value */); | |
1180 | |
1181 std::set<std::string> roiNames; | |
1182 Orthanc::SerializationToolbox::ReadSetOfStrings(roiNames, body, KEY_ROI_NAMES); | |
1183 | |
1184 std::string dicom; | |
1185 if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false)) | |
1186 { | |
1187 throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource); | |
1188 } | |
1189 else | |
1190 { | |
1191 Orthanc::ParsedDicomFile parsed(dicom); | |
1192 StructureSet structureSet(parsed); | |
1193 | |
1194 std::string stl; | |
1195 if (!EncodeStructureSetMesh(stl, structureSet, roiNames, resolution, smooth)) | |
1196 { | |
1197 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, "Cannot encode STL"); | |
1198 } | |
1199 else | |
1200 { | |
1201 std::string content; | |
1202 Orthanc::Toolbox::EncodeDataUriScheme(content, "model/stl", stl); | |
1203 | |
1204 Json::Value create; | |
1205 create["Content"] = content; | |
1206 create["Parent"] = structureSet.HashStudy(); | |
1207 | |
1208 if (body.isMember(KEY_TAGS)) | |
1209 { | |
1210 create[KEY_TAGS] = body[KEY_TAGS]; | |
1211 } | |
1212 else | |
1213 { | |
1214 std::string description; | |
1215 for (std::set<std::string>::const_iterator it = roiNames.begin(); it != roiNames.end(); ++it) | |
1216 { | |
1217 if (!description.empty()) | |
1218 { | |
1219 description += ", "; | |
1220 } | |
1221 | |
1222 description += *it; | |
1223 } | |
1224 | |
1225 create[KEY_TAGS] = Json::objectValue; | |
1226 create[KEY_TAGS]["SeriesDescription"] = description; | |
1227 } | |
1228 | |
1229 Json::Value answer; | |
1230 if (OrthancPlugins::RestApiPost(answer, "/tools/create-dicom", create.toStyledString(), false)) | |
1231 { | |
1232 std::string s = answer.toStyledString(); | |
1233 OrthancPluginAnswerBuffer(OrthancPlugins::GetGlobalContext(), output, s.c_str(), s.size(), Orthanc::MIME_JSON); | |
1234 } | |
1235 else | |
1236 { | |
1237 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest, "Cannot create DICOM from STL"); | |
1238 } | |
1239 } | |
124 } | 1240 } |
125 } | 1241 } |
126 | 1242 |
127 | 1243 |
128 extern "C" | 1244 extern "C" |
148 Orthanc::Logging::InitializePluginContext(context); | 1264 Orthanc::Logging::InitializePluginContext(context); |
149 #else | 1265 #else |
150 Orthanc::Logging::Initialize(context); | 1266 Orthanc::Logging::Initialize(context); |
151 #endif | 1267 #endif |
152 | 1268 |
1269 Orthanc::InitializeFramework("", false); | |
1270 | |
153 hasCreateDicomStl_ = OrthancPlugins::CheckMinimalOrthancVersion(1, 12, 1); | 1271 hasCreateDicomStl_ = OrthancPlugins::CheckMinimalOrthancVersion(1, 12, 1); |
154 | 1272 |
155 if (!hasCreateDicomStl_) | 1273 if (!hasCreateDicomStl_) |
156 { | 1274 { |
157 LOG(WARNING) << "Your version of Orthanc (" << std::string(context->orthancVersion) | 1275 LOG(WARNING) << "Your version of Orthanc (" << std::string(context->orthancVersion) |
158 << ") is insufficient to create DICOM STL, it should be above 1.12.1"; | 1276 << ") is insufficient to create DICOM STL, it should be above 1.12.1"; |
159 } | 1277 } |
160 | 1278 |
161 OrthancPluginSetDescription(context, "STL plugin for Orthanc."); | 1279 OrthancPluginSetDescription(context, "STL plugin for Orthanc."); |
162 | 1280 |
163 OrthancPlugins::RegisterRestCallback<ServeFile>("/stl/(.*)", true); | 1281 OrthancPlugins::RegisterRestCallback<ServeFile>("/stl/app/(.*)", true); |
164 | 1282 |
165 // Extend the default Orthanc Explorer with custom JavaScript for STL | 1283 // Extend the default Orthanc Explorer with custom JavaScript for STL |
166 std::string explorer; | 1284 std::string explorer; |
167 Orthanc::EmbeddedResources::GetFileResource(explorer, Orthanc::EmbeddedResources::ORTHANC_EXPLORER); | 1285 Orthanc::EmbeddedResources::GetFileResource(explorer, Orthanc::EmbeddedResources::ORTHANC_EXPLORER); |
168 OrthancPluginExtendOrthancExplorer(OrthancPlugins::GetGlobalContext(), explorer.c_str()); | 1286 OrthancPluginExtendOrthancExplorer(OrthancPlugins::GetGlobalContext(), explorer.c_str()); |
169 | 1287 |
1288 OrthancPlugins::RegisterRestCallback<ListStructures>("/stl/rt-struct/([0-9a-f-]+)", true); | |
1289 | |
1290 if (hasCreateDicomStl_) | |
1291 { | |
1292 OrthancPlugins::RegisterRestCallback<Encode>("/stl/encode", true); | |
1293 } | |
1294 | |
170 return 0; | 1295 return 0; |
171 } | 1296 } |
172 | 1297 |
173 | 1298 |
174 ORTHANC_PLUGINS_API void OrthancPluginFinalize() | 1299 ORTHANC_PLUGINS_API void OrthancPluginFinalize() |
175 { | 1300 { |
1301 Orthanc::FinalizeFramework(); | |
176 } | 1302 } |
177 | 1303 |
178 | 1304 |
179 ORTHANC_PLUGINS_API const char* OrthancPluginGetName() | 1305 ORTHANC_PLUGINS_API const char* OrthancPluginGetName() |
180 { | 1306 { |