Mercurial > hg > orthanc-stl
comparison Sources/Plugin.cpp @ 32:976da5476810
reorganization
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Thu, 04 Apr 2024 18:35:54 +0200 |
parents | ab231760799d |
children | 2460b376d3f7 |
comparison
equal
deleted
inserted
replaced
31:ab231760799d | 32:976da5476810 |
---|---|
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 "VTKToolbox.h" | |
26 #include "Vector3D.h" | |
27 #include "Toolbox.h" | |
28 #include "Extent2D.h" | |
29 | |
25 #include "../Resources/Orthanc/Plugins/OrthancPluginCppWrapper.h" | 30 #include "../Resources/Orthanc/Plugins/OrthancPluginCppWrapper.h" |
26 | 31 |
27 #include <EmbeddedResources.h> | 32 #include <EmbeddedResources.h> |
28 | 33 |
29 #include <DicomFormat/DicomInstanceHasher.h> | 34 #include <DicomFormat/DicomInstanceHasher.h> |
30 #include <Compression/GzipCompressor.h> | |
31 #include <ChunkedBuffer.h> | |
32 #include <DicomParsing/FromDcmtkBridge.h> | 35 #include <DicomParsing/FromDcmtkBridge.h> |
33 #include <DicomParsing/ParsedDicomFile.h> | 36 #include <DicomParsing/ParsedDicomFile.h> |
34 #include <Images/ImageProcessing.h> | 37 #include <Images/ImageProcessing.h> |
35 #include <Logging.h> | 38 #include <Logging.h> |
36 #include <OrthancFramework.h> | 39 #include <OrthancFramework.h> |
37 #include <SerializationToolbox.h> | 40 #include <SerializationToolbox.h> |
38 #include <SystemToolbox.h> | 41 #include <SystemToolbox.h> |
39 | 42 |
40 #include <vtkImageConstantPad.h> | |
41 #include <vtkImageData.h> | |
42 #include <vtkImageResize.h> | |
43 #include <vtkMarchingCubes.h> | |
44 #include <vtkNew.h> | |
45 #include <vtkPolyData.h> | |
46 #include <vtkPolyDataNormals.h> | |
47 #include <vtkSmoothPolyDataFilter.h> | |
48 #include <vtkTriangle.h> | |
49 | |
50 #include <boost/thread/shared_mutex.hpp> | 43 #include <boost/thread/shared_mutex.hpp> |
51 | |
52 #include <nifti1_io.h> | |
53 | 44 |
54 #define ORTHANC_PLUGIN_NAME "stl" | 45 #define ORTHANC_PLUGIN_NAME "stl" |
55 | 46 |
56 | 47 |
57 // Forward declaration | 48 // Forward declaration |
160 #include <dcmtk/dcmdata/dcfilefo.h> | 151 #include <dcmtk/dcmdata/dcfilefo.h> |
161 #include <dcmtk/dcmdata/dcitem.h> | 152 #include <dcmtk/dcmdata/dcitem.h> |
162 #include <dcmtk/dcmdata/dcsequen.h> | 153 #include <dcmtk/dcmdata/dcsequen.h> |
163 #include <dcmtk/dcmdata/dcuid.h> | 154 #include <dcmtk/dcmdata/dcuid.h> |
164 | 155 |
165 class Extent2D : public boost::noncopyable | |
166 { | |
167 private: | |
168 bool isEmpty_; | |
169 double x1_; | |
170 double y1_; | |
171 double x2_; | |
172 double y2_; | |
173 | |
174 void CheckNotEmpty() const | |
175 { | |
176 if (isEmpty_) | |
177 { | |
178 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
179 } | |
180 } | |
181 | |
182 public: | |
183 Extent2D() : | |
184 isEmpty_(true), | |
185 x1_(0), | |
186 y1_(0), | |
187 x2_(0), | |
188 y2_(0) | |
189 { | |
190 } | |
191 | |
192 bool IsEmpty() const | |
193 { | |
194 return isEmpty_; | |
195 } | |
196 | |
197 double GetMinX() const | |
198 { | |
199 CheckNotEmpty(); | |
200 return x1_; | |
201 } | |
202 | |
203 double GetMaxX() const | |
204 { | |
205 CheckNotEmpty(); | |
206 return x2_; | |
207 } | |
208 | |
209 double GetMinY() const | |
210 { | |
211 CheckNotEmpty(); | |
212 return y1_; | |
213 } | |
214 | |
215 double GetMaxY() const | |
216 { | |
217 CheckNotEmpty(); | |
218 return y2_; | |
219 } | |
220 | |
221 double GetWidth() const | |
222 { | |
223 CheckNotEmpty(); | |
224 return x2_ - x1_; | |
225 } | |
226 | |
227 double GetHeight() const | |
228 { | |
229 CheckNotEmpty(); | |
230 return y2_ - y1_; | |
231 } | |
232 | |
233 void Add(double x, | |
234 double y) | |
235 { | |
236 if (isEmpty_) | |
237 { | |
238 x1_ = x2_ = x; | |
239 y1_ = y2_ = y; | |
240 isEmpty_ = false; | |
241 } | |
242 else | |
243 { | |
244 x1_ = std::min(x1_, x); | |
245 x2_ = std::max(x2_, x); | |
246 y1_ = std::min(y1_, y); | |
247 y2_ = std::max(y2_, y); | |
248 } | |
249 } | |
250 }; | |
251 | 156 |
252 static std::string GetStringValue(DcmItem& item, | 157 static std::string GetStringValue(DcmItem& item, |
253 const DcmTagKey& key) | 158 const DcmTagKey& key) |
254 { | 159 { |
255 const char* s = NULL; | 160 const char* s = NULL; |
287 else | 192 else |
288 { | 193 { |
289 target.insert(GetStringValue(*item, DCM_ROIName)); | 194 target.insert(GetStringValue(*item, DCM_ROIName)); |
290 } | 195 } |
291 } | 196 } |
292 } | |
293 | |
294 | |
295 static bool IsNear(double a, | |
296 double b) | |
297 { | |
298 return std::abs(a - b) < 10.0 * std::numeric_limits<double>::epsilon(); | |
299 } | |
300 | |
301 | |
302 class Vector3D | |
303 { | |
304 private: | |
305 double x_; | |
306 double y_; | |
307 double z_; | |
308 | |
309 public: | |
310 Vector3D() : | |
311 x_(0), | |
312 y_(0), | |
313 z_(0) | |
314 { | |
315 } | |
316 | |
317 Vector3D(double x, | |
318 double y, | |
319 double z) : | |
320 x_(x), | |
321 y_(y), | |
322 z_(z) | |
323 { | |
324 } | |
325 | |
326 Vector3D(const Vector3D& from, | |
327 const Vector3D& to) : | |
328 x_(to.x_ - from.x_), | |
329 y_(to.y_ - from.y_), | |
330 z_(to.z_ - from.z_) | |
331 { | |
332 } | |
333 | |
334 double GetX() const | |
335 { | |
336 return x_; | |
337 } | |
338 | |
339 double GetY() const | |
340 { | |
341 return y_; | |
342 } | |
343 | |
344 double GetZ() const | |
345 { | |
346 return z_; | |
347 } | |
348 | |
349 double ComputeSquaredNorm() const | |
350 { | |
351 return x_ * x_ + y_ * y_ + z_ * z_; | |
352 } | |
353 | |
354 double ComputeNorm() const | |
355 { | |
356 return sqrt(ComputeSquaredNorm()); | |
357 } | |
358 | |
359 void Normalize() | |
360 { | |
361 double norm = ComputeNorm(); | |
362 if (!IsNear(norm, 0)) | |
363 { | |
364 x_ /= norm; | |
365 y_ /= norm; | |
366 z_ /= norm; | |
367 } | |
368 } | |
369 | |
370 static Vector3D CrossProduct(const Vector3D& u, | |
371 const Vector3D& v) | |
372 { | |
373 return Vector3D(u.GetY() * v.GetZ() - u.GetZ() * v.GetY(), | |
374 u.GetZ() * v.GetX() - u.GetX() * v.GetZ(), | |
375 u.GetX() * v.GetY() - u.GetY() * v.GetX()); | |
376 } | |
377 | |
378 static double DotProduct(const Vector3D& a, | |
379 const Vector3D& b) | |
380 { | |
381 return a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ(); | |
382 } | |
383 | |
384 static bool AreParallel(const Vector3D& a, | |
385 const Vector3D& b) | |
386 { | |
387 if (IsNear(a.ComputeSquaredNorm(), 1) && | |
388 IsNear(b.ComputeSquaredNorm(), 1)) | |
389 { | |
390 return IsNear(std::abs(Vector3D::DotProduct(a, b)), 1); | |
391 } | |
392 else | |
393 { | |
394 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange, | |
395 "Only applicable to normalized vectors"); | |
396 } | |
397 } | |
398 }; | |
399 | |
400 | |
401 | |
402 static bool MyParseDouble(double& value, | |
403 const std::string& s) | |
404 { | |
405 #if 1 | |
406 // This version is much faster, as "ParseDouble()" internally uses | |
407 // "boost::lexical_cast<double>()" | |
408 char* end = NULL; | |
409 value = strtod(s.c_str(), &end); | |
410 return (end == s.c_str() + s.size()); | |
411 #else | |
412 return Orthanc::SerializationToolbox::ParseDouble(value, s); | |
413 #endif | |
414 } | 197 } |
415 | 198 |
416 | 199 |
417 class StructurePolygon : public boost::noncopyable | 200 class StructurePolygon : public boost::noncopyable |
418 { | 201 { |
494 points_.reserve(countPoints); | 277 points_.reserve(countPoints); |
495 | 278 |
496 for (size_t i = 0; i < tokens.size(); i += 3) | 279 for (size_t i = 0; i < tokens.size(); i += 3) |
497 { | 280 { |
498 double x, y, z; | 281 double x, y, z; |
499 if (!MyParseDouble(x, tokens[i]) || | 282 if (!Toolbox::MyParseDouble(x, tokens[i]) || |
500 !MyParseDouble(y, tokens[i + 1]) || | 283 !Toolbox::MyParseDouble(y, tokens[i + 1]) || |
501 !MyParseDouble(z, tokens[i + 2])) | 284 !Toolbox::MyParseDouble(z, tokens[i + 2])) |
502 { | 285 { |
503 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 286 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
504 } | 287 } |
505 | 288 |
506 points_.push_back(Vector3D(x, y, z)); | 289 points_.push_back(Vector3D(x, y, z)); |
548 | 331 |
549 for (size_t i = 0; i < points_.size(); i++) | 332 for (size_t i = 0; i < points_.size(); i++) |
550 { | 333 { |
551 normal = Vector3D::CrossProduct(Vector3D(points_[1], points_[0]), | 334 normal = Vector3D::CrossProduct(Vector3D(points_[1], points_[0]), |
552 Vector3D(points_[2], points_[0])); | 335 Vector3D(points_[2], points_[0])); |
553 if (!IsNear(normal.ComputeNorm(), 0)) | 336 if (!Toolbox::IsNear(normal.ComputeNorm(), 0)) |
554 { | 337 { |
555 normal.Normalize(); | 338 normal.Normalize(); |
556 hasNormal = true; | 339 hasNormal = true; |
557 } | 340 } |
558 } | 341 } |
565 double a = Vector3D::DotProduct(points_[0], normal); | 348 double a = Vector3D::DotProduct(points_[0], normal); |
566 | 349 |
567 for (size_t i = 1; i < points_.size(); i++) | 350 for (size_t i = 1; i < points_.size(); i++) |
568 { | 351 { |
569 double b = Vector3D::DotProduct(points_[i], normal); | 352 double b = Vector3D::DotProduct(points_[i], normal); |
570 if (!IsNear(a, b)) | 353 if (!Toolbox::IsNear(a, b)) |
571 { | 354 { |
572 return false; | 355 return false; |
573 } | 356 } |
574 } | 357 } |
575 | 358 |
578 | 361 |
579 void Add(Extent2D& extent, | 362 void Add(Extent2D& extent, |
580 const Vector3D& axisX, | 363 const Vector3D& axisX, |
581 const Vector3D& axisY) const | 364 const Vector3D& axisY) const |
582 { | 365 { |
583 assert(IsNear(1, axisX.ComputeNorm())); | 366 assert(Toolbox::IsNear(1, axisX.ComputeNorm())); |
584 assert(IsNear(1, axisY.ComputeNorm())); | 367 assert(Toolbox::IsNear(1, axisY.ComputeNorm())); |
585 | 368 |
586 for (size_t i = 0; i < points_.size(); i++) | 369 for (size_t i = 0; i < points_.size(); i++) |
587 { | 370 { |
588 extent.Add(Vector3D::DotProduct(axisX, points_[i]), | 371 extent.Add(Vector3D::DotProduct(axisX, points_[i]), |
589 Vector3D::DotProduct(axisY, points_[i])); | 372 Vector3D::DotProduct(axisY, points_[i])); |
590 } | 373 } |
591 } | 374 } |
592 }; | 375 }; |
593 | 376 |
594 | |
595 | |
596 struct IsNearPredicate | |
597 { | |
598 bool operator() (const double& a, | |
599 const double& b) | |
600 { | |
601 return IsNear(a, b); | |
602 } | |
603 }; | |
604 | |
605 static void RemoveDuplicateValues(std::vector<double>& v) | |
606 { | |
607 IsNearPredicate predicate; | |
608 std::vector<double>::iterator last = std::unique(v.begin(), v.end(), predicate); | |
609 v.erase(last, v.end()); | |
610 } | |
611 | 377 |
612 | 378 |
613 class StructureSet : public boost::noncopyable | 379 class StructureSet : public boost::noncopyable |
614 { | 380 { |
615 private: | 381 private: |
787 assert(slicesSpacing_ > 0 && | 553 assert(slicesSpacing_ > 0 && |
788 minProjectionAlongNormal_ < maxProjectionAlongNormal_); | 554 minProjectionAlongNormal_ < maxProjectionAlongNormal_); |
789 | 555 |
790 double d = (z - minProjectionAlongNormal_) / slicesSpacing_; | 556 double d = (z - minProjectionAlongNormal_) / slicesSpacing_; |
791 | 557 |
792 if (IsNear(d, round(d))) | 558 if (Toolbox::IsNear(d, round(d))) |
793 { | 559 { |
794 if (d < 0.0 || | 560 if (d < 0.0 || |
795 d > static_cast<double>(slicesCount_) - 1.0) | 561 d > static_cast<double>(slicesCount_) - 1.0) |
796 { | 562 { |
797 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | 563 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); |
855 } | 621 } |
856 | 622 |
857 // Only keep unique projections | 623 // Only keep unique projections |
858 | 624 |
859 std::sort(projections.begin(), projections.end()); | 625 std::sort(projections.begin(), projections.end()); |
860 RemoveDuplicateValues(projections); | 626 Toolbox::RemoveDuplicateValues(projections); |
861 assert(!projections.empty()); | 627 assert(!projections.empty()); |
862 | 628 |
863 if (projections.size() == 1) | 629 if (projections.size() == 1) |
864 { | 630 { |
865 // Volume with one single slice | 631 // Volume with one single slice |
882 spacings[i] = projections[i + 1] - projections[i]; | 648 spacings[i] = projections[i + 1] - projections[i]; |
883 assert(spacings[i] > 0); | 649 assert(spacings[i] > 0); |
884 } | 650 } |
885 | 651 |
886 std::sort(spacings.begin(), spacings.end()); | 652 std::sort(spacings.begin(), spacings.end()); |
887 RemoveDuplicateValues(spacings); | 653 Toolbox::RemoveDuplicateValues(spacings); |
888 | 654 |
889 if (spacings.empty()) | 655 if (spacings.empty()) |
890 { | 656 { |
891 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | 657 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); |
892 } | 658 } |
919 it++; | 685 it++; |
920 | 686 |
921 while (it != candidates.end()) | 687 while (it != candidates.end()) |
922 { | 688 { |
923 double d = (projections[*it] - projections[reference]) / slicesSpacing_; | 689 double d = (projections[*it] - projections[reference]) / slicesSpacing_; |
924 if (IsNear(d, round(d))) | 690 if (Toolbox::IsNear(d, round(d))) |
925 { | 691 { |
926 countSupport ++; | 692 countSupport ++; |
927 } | 693 } |
928 else | 694 else |
929 { | 695 { |
958 maxProjectionAlongNormal_ = bestProjection; | 724 maxProjectionAlongNormal_ = bestProjection; |
959 | 725 |
960 for (size_t i = 0; i < projections.size(); i++) | 726 for (size_t i = 0; i < projections.size(); i++) |
961 { | 727 { |
962 double d = (projections[i] - bestProjection) / slicesSpacing_; | 728 double d = (projections[i] - bestProjection) / slicesSpacing_; |
963 if (IsNear(d, round(d))) | 729 if (Toolbox::IsNear(d, round(d))) |
964 { | 730 { |
965 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]); | 731 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]); |
966 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]); | 732 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]); |
967 } | 733 } |
968 } | 734 } |
969 | 735 |
970 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; | 736 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; |
971 if (IsNear(d, round(d))) | 737 if (Toolbox::IsNear(d, round(d))) |
972 { | 738 { |
973 slicesCount_ = static_cast<size_t>(round(d)) + 1; | 739 slicesCount_ = static_cast<size_t>(round(d)) + 1; |
974 } | 740 } |
975 else | 741 else |
976 { | 742 { |
1075 } | 841 } |
1076 } | 842 } |
1077 }; | 843 }; |
1078 | 844 |
1079 | 845 |
1080 static void WriteFloat(Orthanc::ChunkedBuffer& buffer, | |
1081 float value, | |
1082 Orthanc::Endianness endianness) | |
1083 { | |
1084 switch (endianness) | |
1085 { | |
1086 case Orthanc::Endianness_Little: | |
1087 buffer.AddChunk(&value, sizeof(float)); | |
1088 break; | |
1089 | |
1090 case Orthanc::Endianness_Big: | |
1091 { | |
1092 uint8_t tmp[4]; | |
1093 tmp[0] = reinterpret_cast<uint8_t*>(&value) [3]; | |
1094 tmp[1] = reinterpret_cast<uint8_t*>(&value) [2]; | |
1095 tmp[2] = reinterpret_cast<uint8_t*>(&value) [1]; | |
1096 tmp[3] = reinterpret_cast<uint8_t*>(&value) [0]; | |
1097 buffer.AddChunk(&tmp, sizeof(float)); | |
1098 break; | |
1099 } | |
1100 | |
1101 default: | |
1102 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
1103 } | |
1104 } | |
1105 | |
1106 | |
1107 static void WriteInteger(Orthanc::ChunkedBuffer& buffer, | |
1108 uint32_t value, | |
1109 Orthanc::Endianness endianness) | |
1110 { | |
1111 switch (endianness) | |
1112 { | |
1113 case Orthanc::Endianness_Little: | |
1114 buffer.AddChunk(&value, sizeof(uint32_t)); | |
1115 break; | |
1116 | |
1117 case Orthanc::Endianness_Big: | |
1118 { | |
1119 uint8_t tmp[4]; | |
1120 tmp[0] = reinterpret_cast<uint8_t*>(&value) [3]; | |
1121 tmp[1] = reinterpret_cast<uint8_t*>(&value) [2]; | |
1122 tmp[2] = reinterpret_cast<uint8_t*>(&value) [1]; | |
1123 tmp[3] = reinterpret_cast<uint8_t*>(&value) [0]; | |
1124 buffer.AddChunk(&tmp, sizeof(uint32_t)); | |
1125 break; | |
1126 } | |
1127 | |
1128 default: | |
1129 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
1130 } | |
1131 } | |
1132 | |
1133 | |
1134 static void EncodeSTL(std::string& target /* out */, | |
1135 vtkPolyData& mesh /* in */) | |
1136 { | |
1137 const Orthanc::Endianness endianness = Orthanc::Toolbox::DetectEndianness(); | |
1138 | |
1139 Orthanc::ChunkedBuffer buffer; | |
1140 | |
1141 uint8_t header[80]; | |
1142 memset(header, 0, sizeof(header)); | |
1143 buffer.AddChunk(header, sizeof(header)); // This doesn't depend on endianness | |
1144 | |
1145 WriteInteger(buffer, mesh.GetNumberOfCells(), endianness); | |
1146 | |
1147 for (vtkIdType i = 0; i < mesh.GetNumberOfCells(); i++) | |
1148 { | |
1149 vtkCell* cell = mesh.GetCell(i); | |
1150 vtkTriangle* triangle = dynamic_cast<vtkTriangle*>(cell); | |
1151 | |
1152 double p0[3]; | |
1153 double p1[3]; | |
1154 double p2[3]; | |
1155 triangle->GetPoints()->GetPoint(0, p0); | |
1156 triangle->GetPoints()->GetPoint(1, p1); | |
1157 triangle->GetPoints()->GetPoint(2, p2); | |
1158 | |
1159 double normal[3]; | |
1160 vtkTriangle::ComputeNormal(p0, p1, p2, normal); | |
1161 | |
1162 WriteFloat(buffer, normal[0], endianness); | |
1163 WriteFloat(buffer, normal[1], endianness); | |
1164 WriteFloat(buffer, normal[2], endianness); | |
1165 WriteFloat(buffer, p0[0], endianness); | |
1166 WriteFloat(buffer, p0[1], endianness); | |
1167 WriteFloat(buffer, p0[2], endianness); | |
1168 WriteFloat(buffer, p1[0], endianness); | |
1169 WriteFloat(buffer, p1[1], endianness); | |
1170 WriteFloat(buffer, p1[2], endianness); | |
1171 WriteFloat(buffer, p2[0], endianness); | |
1172 WriteFloat(buffer, p2[1], endianness); | |
1173 WriteFloat(buffer, p2[2], endianness); | |
1174 | |
1175 uint16_t a = 0; | |
1176 buffer.AddChunk(&a, sizeof(a)); // This doesn't depend on endianness | |
1177 } | |
1178 | |
1179 buffer.Flatten(target); | |
1180 } | |
1181 | |
1182 | |
1183 bool EncodeVolume(std::string& stl, | |
1184 vtkImageData* volume, | |
1185 unsigned int resolution, | |
1186 bool smooth) | |
1187 { | |
1188 if (volume == NULL) | |
1189 { | |
1190 throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer); | |
1191 } | |
1192 | |
1193 vtkNew<vtkImageResize> resize; | |
1194 resize->SetOutputDimensions(resolution, resolution, resolution); | |
1195 resize->SetInputData(volume); | |
1196 resize->Update(); | |
1197 | |
1198 vtkNew<vtkImageConstantPad> padding; | |
1199 padding->SetConstant(0); | |
1200 padding->SetOutputNumberOfScalarComponents(1); | |
1201 padding->SetOutputWholeExtent(-1, resolution, -1, resolution, -1, resolution); | |
1202 padding->SetInputData(resize->GetOutput()); | |
1203 padding->Update(); | |
1204 | |
1205 double range[2]; | |
1206 padding->GetOutput()->GetScalarRange(range); | |
1207 | |
1208 const double isoValue = (range[0] + range[1]) / 2.0; | |
1209 | |
1210 vtkNew<vtkMarchingCubes> surface; | |
1211 surface->SetInputData(padding->GetOutput()); | |
1212 surface->ComputeNormalsOn(); | |
1213 surface->SetValue(0, isoValue); | |
1214 surface->Update(); | |
1215 | |
1216 if (smooth) | |
1217 { | |
1218 vtkNew<vtkSmoothPolyDataFilter> smoothFilter; | |
1219 // Apply volume smoothing | |
1220 // https://examples.vtk.org/site/Cxx/PolyData/SmoothPolyDataFilter/ | |
1221 smoothFilter->SetInputConnection(surface->GetOutputPort()); | |
1222 smoothFilter->SetNumberOfIterations(15); | |
1223 smoothFilter->SetRelaxationFactor(0.1); | |
1224 smoothFilter->FeatureEdgeSmoothingOff(); | |
1225 smoothFilter->BoundarySmoothingOn(); | |
1226 smoothFilter->Update(); | |
1227 | |
1228 vtkNew<vtkPolyDataNormals> normalGenerator; | |
1229 normalGenerator->SetInputConnection(smoothFilter->GetOutputPort()); | |
1230 normalGenerator->ComputePointNormalsOn(); | |
1231 normalGenerator->ComputeCellNormalsOn(); | |
1232 normalGenerator->Update(); | |
1233 | |
1234 EncodeSTL(stl, *normalGenerator->GetOutput()); | |
1235 } | |
1236 else | |
1237 { | |
1238 EncodeSTL(stl, *surface->GetOutput()); | |
1239 } | |
1240 | |
1241 return true; | |
1242 } | |
1243 | |
1244 | |
1245 static Orthanc::ParsedDicomFile* LoadInstance(const std::string& instanceId) | 846 static Orthanc::ParsedDicomFile* LoadInstance(const std::string& instanceId) |
1246 { | 847 { |
1247 std::string dicom; | 848 std::string dicom; |
1248 | 849 |
1249 if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false)) | 850 if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false)) |
1299 Orthanc::Toolbox::TokenizeString(items, imageOrientation, '\\'); | 900 Orthanc::Toolbox::TokenizeString(items, imageOrientation, '\\'); |
1300 | 901 |
1301 double x1, x2, x3, y1, y2, y3; | 902 double x1, x2, x3, y1, y2, y3; |
1302 | 903 |
1303 if (items.size() == 6 && | 904 if (items.size() == 6 && |
1304 MyParseDouble(x1, items[0]) && | 905 Toolbox::MyParseDouble(x1, items[0]) && |
1305 MyParseDouble(x2, items[1]) && | 906 Toolbox::MyParseDouble(x2, items[1]) && |
1306 MyParseDouble(x3, items[2]) && | 907 Toolbox::MyParseDouble(x3, items[2]) && |
1307 MyParseDouble(y1, items[3]) && | 908 Toolbox::MyParseDouble(y1, items[3]) && |
1308 MyParseDouble(y2, items[4]) && | 909 Toolbox::MyParseDouble(y2, items[4]) && |
1309 MyParseDouble(y3, items[5])) | 910 Toolbox::MyParseDouble(y3, items[5])) |
1310 { | 911 { |
1311 axisX = Vector3D(x1, x2, x3); | 912 axisX = Vector3D(x1, x2, x3); |
1312 axisY = Vector3D(y1, y2, y3); | 913 axisY = Vector3D(y1, y2, y3); |
1313 } | 914 } |
1314 } | 915 } |
1335 structureSet.GetPolygonsCount() < 1) | 936 structureSet.GetPolygonsCount() < 1) |
1336 { | 937 { |
1337 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | 938 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); |
1338 } | 939 } |
1339 | 940 |
1340 if (!IsNear(1, geometry.GetSlicesNormal().ComputeNorm())) | 941 if (!Toolbox::IsNear(1, geometry.GetSlicesNormal().ComputeNorm())) |
1341 { | 942 { |
1342 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | 943 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); |
1343 } | 944 } |
1344 | 945 |
1345 Vector3D axisX, axisY; | 946 Vector3D axisX, axisY; |
1346 GetReferencedVolumeAxes(axisX, axisY, structureSet); | 947 GetReferencedVolumeAxes(axisX, axisY, structureSet); |
1347 | 948 |
1348 Vector3D axisZ = Vector3D::CrossProduct(axisX, axisY); | 949 Vector3D axisZ = Vector3D::CrossProduct(axisX, axisY); |
1349 | 950 |
1350 if (!IsNear(1, axisX.ComputeNorm()) || | 951 if (!Toolbox::IsNear(1, axisX.ComputeNorm()) || |
1351 !IsNear(1, axisY.ComputeNorm()) || | 952 !Toolbox::IsNear(1, axisY.ComputeNorm()) || |
1352 !Vector3D::AreParallel(axisZ, geometry.GetSlicesNormal())) | 953 !Vector3D::AreParallel(axisZ, geometry.GetSlicesNormal())) |
1353 { | 954 { |
1354 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | 955 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); |
1355 } | 956 } |
1356 | 957 |
1415 } | 1016 } |
1416 } | 1017 } |
1417 | 1018 |
1418 volume->SetSpacing(pixelSpacingX, pixelSpacingY, pixelSpacingZ); | 1019 volume->SetSpacing(pixelSpacingX, pixelSpacingY, pixelSpacingZ); |
1419 | 1020 |
1420 return EncodeVolume(stl, volume.Get(), resolution, smooth); | 1021 return VTKToolbox::EncodeVolume(stl, volume.Get(), resolution, smooth); |
1421 } | 1022 } |
1422 | 1023 |
1423 | 1024 |
1424 void ListStructures(OrthancPluginRestOutput* output, | 1025 void ListStructures(OrthancPluginRestOutput* output, |
1425 const char* url, | 1026 const char* url, |
1672 } | 1273 } |
1673 } | 1274 } |
1674 | 1275 |
1675 | 1276 |
1676 | 1277 |
1677 class NiftiHeader : public boost::noncopyable | |
1678 { | |
1679 private: | |
1680 nifti_image* image_; | |
1681 | |
1682 public: | |
1683 NiftiHeader(const std::string& nifti) | |
1684 { | |
1685 nifti_1_header header; | |
1686 if (nifti.size() < sizeof(header)) | |
1687 { | |
1688 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
1689 } | |
1690 | |
1691 memcpy(&header, nifti.c_str(), sizeof(header)); | |
1692 if (!nifti_hdr_looks_good(&header)) | |
1693 { | |
1694 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
1695 } | |
1696 | |
1697 image_ = nifti_convert_nhdr2nim(header, "dummy_filename"); | |
1698 if (image_ == NULL) | |
1699 { | |
1700 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
1701 } | |
1702 } | |
1703 | |
1704 ~NiftiHeader() | |
1705 { | |
1706 nifti_image_free(image_); | |
1707 } | |
1708 | |
1709 const nifti_image& GetInfo() const | |
1710 { | |
1711 assert(image_ != NULL); | |
1712 return *image_; | |
1713 } | |
1714 }; | |
1715 | |
1716 | |
1717 static void LoadNifti(vtkImageData* volume, | |
1718 std::string& nifti) | |
1719 { | |
1720 if (volume == NULL) | |
1721 { | |
1722 throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer); | |
1723 } | |
1724 | |
1725 const uint8_t* p = reinterpret_cast<const uint8_t*>(nifti.c_str()); | |
1726 | |
1727 if (nifti.size() >= 2 && | |
1728 p[0] == 0x1f && | |
1729 p[1] == 0x8b) | |
1730 { | |
1731 Orthanc::GzipCompressor compressor; | |
1732 std::string uncompressed; | |
1733 Orthanc::IBufferCompressor::Uncompress(uncompressed, compressor, nifti); | |
1734 nifti.swap(uncompressed); | |
1735 } | |
1736 | |
1737 NiftiHeader header(nifti); | |
1738 | |
1739 if (header.GetInfo().ndim != 3) | |
1740 { | |
1741 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, | |
1742 "Only 3D NIfTI volumes are allowed"); | |
1743 } | |
1744 | |
1745 size_t itemSize; | |
1746 int vtkType; | |
1747 | |
1748 switch (header.GetInfo().datatype) | |
1749 { | |
1750 case DT_UNSIGNED_CHAR: | |
1751 itemSize = 1; | |
1752 vtkType = VTK_UNSIGNED_CHAR; | |
1753 break; | |
1754 | |
1755 case DT_FLOAT: | |
1756 itemSize = sizeof(float); | |
1757 vtkType = VTK_FLOAT; | |
1758 break; | |
1759 | |
1760 case DT_DOUBLE: | |
1761 itemSize = sizeof(double); | |
1762 vtkType = VTK_DOUBLE; | |
1763 break; | |
1764 | |
1765 default: | |
1766 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
1767 } | |
1768 | |
1769 assert(static_cast<int>(header.GetInfo().nvox) == header.GetInfo().nx * header.GetInfo().ny * header.GetInfo().nz); | |
1770 | |
1771 const size_t pixelDataOffset = sizeof(nifti_1_header) + 4 /* extension */; | |
1772 | |
1773 if (nifti.size() != pixelDataOffset + header.GetInfo().nvox * itemSize) | |
1774 { | |
1775 throw Orthanc::OrthancException(Orthanc::ErrorCode_CorruptedFile); | |
1776 } | |
1777 | |
1778 volume->SetDimensions(header.GetInfo().nx, header.GetInfo().ny, header.GetInfo().nz); | |
1779 volume->AllocateScalars(vtkType, 1); | |
1780 volume->SetSpacing(header.GetInfo().dx, header.GetInfo().dy, header.GetInfo().dz); | |
1781 memcpy(volume->GetScalarPointer(), &nifti[pixelDataOffset], header.GetInfo().nvox * itemSize); | |
1782 } | |
1783 | |
1784 | |
1785 void EncodeNifti(OrthancPluginRestOutput* output, | 1278 void EncodeNifti(OrthancPluginRestOutput* output, |
1786 const char* url, | 1279 const char* url, |
1787 const OrthancPluginHttpRequest* request) | 1280 const OrthancPluginHttpRequest* request) |
1788 { | 1281 { |
1789 static const char* const KEY_NIFTI = "Nifti"; | 1282 static const char* const KEY_NIFTI = "Nifti"; |
1816 const unsigned int resolution = (body.isMember(KEY_RESOLUTION) ? | 1309 const unsigned int resolution = (body.isMember(KEY_RESOLUTION) ? |
1817 Orthanc::SerializationToolbox::ReadUnsignedInteger(body, KEY_RESOLUTION) : | 1310 Orthanc::SerializationToolbox::ReadUnsignedInteger(body, KEY_RESOLUTION) : |
1818 256 /* default value */); | 1311 256 /* default value */); |
1819 | 1312 |
1820 vtkNew<vtkImageData> volume; | 1313 vtkNew<vtkImageData> volume; |
1821 LoadNifti(volume.Get(), nifti); | 1314 VTKToolbox::LoadNifti(volume.Get(), nifti); |
1822 | 1315 |
1823 std::string stl; | 1316 std::string stl; |
1824 if (!EncodeVolume(stl, volume.Get(), resolution, smooth)) | 1317 if (!VTKToolbox::EncodeVolume(stl, volume.Get(), resolution, smooth)) |
1825 { | 1318 { |
1826 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, "Cannot encode STL from NIfTI"); | 1319 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, "Cannot encode STL from NIfTI"); |
1827 } | 1320 } |
1828 else | 1321 else |
1829 { | 1322 { |