Mercurial > hg > orthanc-stl

comparison Sources/Plugin.cpp @ 32:976da5476810

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 {