|
32
|
1 /**
|
|
|
2 * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium
|
|
|
3 * SPDX-License-Identifier: GPL-3.0-or-later
|
|
|
4 */
|
|
|
5
|
|
|
6 /**
|
|
|
7 * STL plugin for Orthanc
|
|
|
8 * Copyright (C) 2023-2024 Sebastien Jodogne, UCLouvain, Belgium
|
|
|
9 *
|
|
|
10 * This program is free software: you can redistribute it and/or
|
|
|
11 * modify it under the terms of the GNU General Public License as
|
|
|
12 * published by the Free Software Foundation, either version 3 of the
|
|
|
13 * License, or (at your option) any later version.
|
|
|
14 *
|
|
|
15 * This program is distributed in the hope that it will be useful, but
|
|
|
16 * WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
18 * General Public License for more details.
|
|
|
19 *
|
|
|
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/>.
|
|
|
22 **/
|
|
|
23
|
|
|
24
|
|
|
25 #include "VTKToolbox.h"
|
|
|
26
|
|
|
27 #include <ChunkedBuffer.h>
|
|
|
28 #include <Compression/GzipCompressor.h>
|
|
|
29 #include <Endianness.h>
|
|
|
30 #include <OrthancException.h>
|
|
|
31 #include <Toolbox.h>
|
|
|
32
|
|
|
33 #include <vtkTriangle.h>
|
|
|
34 #include <vtkImageResize.h>
|
|
|
35 #include <vtkImageConstantPad.h>
|
|
|
36 #include <vtkMarchingCubes.h>
|
|
|
37 #include <vtkSmoothPolyDataFilter.h>
|
|
|
38 #include <vtkPolyDataNormals.h>
|
|
|
39
|
|
|
40 #include <nifti1_io.h>
|
|
|
41
|
|
|
42
|
|
|
43 namespace VTKToolbox
|
|
|
44 {
|
|
|
45 static void WriteFloat(Orthanc::ChunkedBuffer& buffer,
|
|
|
46 float value,
|
|
|
47 Orthanc::Endianness endianness)
|
|
|
48 {
|
|
|
49 switch (endianness)
|
|
|
50 {
|
|
|
51 case Orthanc::Endianness_Little:
|
|
|
52 buffer.AddChunk(&value, sizeof(float));
|
|
|
53 break;
|
|
|
54
|
|
|
55 case Orthanc::Endianness_Big:
|
|
|
56 {
|
|
|
57 uint8_t tmp[4];
|
|
|
58 tmp[0] = reinterpret_cast<uint8_t*>(&value) [3];
|
|
|
59 tmp[1] = reinterpret_cast<uint8_t*>(&value) [2];
|
|
|
60 tmp[2] = reinterpret_cast<uint8_t*>(&value) [1];
|
|
|
61 tmp[3] = reinterpret_cast<uint8_t*>(&value) [0];
|
|
|
62 buffer.AddChunk(&tmp, sizeof(float));
|
|
|
63 break;
|
|
|
64 }
|
|
|
65
|
|
|
66 default:
|
|
|
67 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
|
|
|
68 }
|
|
|
69 }
|
|
|
70
|
|
|
71
|
|
|
72 static void WriteInteger(Orthanc::ChunkedBuffer& buffer,
|
|
|
73 uint32_t value,
|
|
|
74 Orthanc::Endianness endianness)
|
|
|
75 {
|
|
|
76 switch (endianness)
|
|
|
77 {
|
|
|
78 case Orthanc::Endianness_Little:
|
|
|
79 buffer.AddChunk(&value, sizeof(uint32_t));
|
|
|
80 break;
|
|
|
81
|
|
|
82 case Orthanc::Endianness_Big:
|
|
|
83 {
|
|
|
84 uint8_t tmp[4];
|
|
|
85 tmp[0] = reinterpret_cast<uint8_t*>(&value) [3];
|
|
|
86 tmp[1] = reinterpret_cast<uint8_t*>(&value) [2];
|
|
|
87 tmp[2] = reinterpret_cast<uint8_t*>(&value) [1];
|
|
|
88 tmp[3] = reinterpret_cast<uint8_t*>(&value) [0];
|
|
|
89 buffer.AddChunk(&tmp, sizeof(uint32_t));
|
|
|
90 break;
|
|
|
91 }
|
|
|
92
|
|
|
93 default:
|
|
|
94 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
|
|
|
95 }
|
|
|
96 }
|
|
|
97
|
|
|
98
|
|
|
99 void EncodeSTL(std::string& target /* out */,
|
|
|
100 vtkPolyData& mesh /* in */)
|
|
|
101 {
|
|
|
102 const Orthanc::Endianness endianness = Orthanc::Toolbox::DetectEndianness();
|
|
|
103
|
|
|
104 Orthanc::ChunkedBuffer buffer;
|
|
|
105
|
|
|
106 uint8_t header[80];
|
|
|
107 memset(header, 0, sizeof(header));
|
|
|
108 buffer.AddChunk(header, sizeof(header)); // This doesn't depend on endianness
|
|
|
109
|
|
|
110 WriteInteger(buffer, mesh.GetNumberOfCells(), endianness);
|
|
|
111
|
|
|
112 for (vtkIdType i = 0; i < mesh.GetNumberOfCells(); i++)
|
|
|
113 {
|
|
|
114 vtkCell* cell = mesh.GetCell(i);
|
|
|
115 vtkTriangle* triangle = dynamic_cast<vtkTriangle*>(cell);
|
|
|
116
|
|
|
117 double p0[3];
|
|
|
118 double p1[3];
|
|
|
119 double p2[3];
|
|
|
120 triangle->GetPoints()->GetPoint(0, p0);
|
|
|
121 triangle->GetPoints()->GetPoint(1, p1);
|
|
|
122 triangle->GetPoints()->GetPoint(2, p2);
|
|
|
123
|
|
|
124 double normal[3];
|
|
|
125 vtkTriangle::ComputeNormal(p0, p1, p2, normal);
|
|
|
126
|
|
|
127 WriteFloat(buffer, normal[0], endianness);
|
|
|
128 WriteFloat(buffer, normal[1], endianness);
|
|
|
129 WriteFloat(buffer, normal[2], endianness);
|
|
|
130 WriteFloat(buffer, p0[0], endianness);
|
|
|
131 WriteFloat(buffer, p0[1], endianness);
|
|
|
132 WriteFloat(buffer, p0[2], endianness);
|
|
|
133 WriteFloat(buffer, p1[0], endianness);
|
|
|
134 WriteFloat(buffer, p1[1], endianness);
|
|
|
135 WriteFloat(buffer, p1[2], endianness);
|
|
|
136 WriteFloat(buffer, p2[0], endianness);
|
|
|
137 WriteFloat(buffer, p2[1], endianness);
|
|
|
138 WriteFloat(buffer, p2[2], endianness);
|
|
|
139
|
|
|
140 uint16_t a = 0;
|
|
|
141 buffer.AddChunk(&a, sizeof(a)); // This doesn't depend on endianness
|
|
|
142 }
|
|
|
143
|
|
|
144 buffer.Flatten(target);
|
|
|
145 }
|
|
|
146
|
|
|
147
|
|
|
148 bool EncodeVolume(std::string& stl,
|
|
|
149 vtkImageData* volume,
|
|
|
150 unsigned int resolution,
|
|
|
151 bool smooth)
|
|
|
152 {
|
|
|
153 if (volume == NULL)
|
|
|
154 {
|
|
|
155 throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer);
|
|
|
156 }
|
|
|
157
|
|
|
158 vtkNew<vtkImageResize> resize;
|
|
|
159 resize->SetOutputDimensions(resolution, resolution, resolution);
|
|
|
160 resize->SetInputData(volume);
|
|
|
161 resize->Update();
|
|
|
162
|
|
|
163 vtkNew<vtkImageConstantPad> padding;
|
|
|
164 padding->SetConstant(0);
|
|
|
165 padding->SetOutputNumberOfScalarComponents(1);
|
|
|
166 padding->SetOutputWholeExtent(-1, resolution, -1, resolution, -1, resolution);
|
|
|
167 padding->SetInputData(resize->GetOutput());
|
|
|
168 padding->Update();
|
|
|
169
|
|
|
170 double range[2];
|
|
|
171 padding->GetOutput()->GetScalarRange(range);
|
|
|
172
|
|
|
173 const double isoValue = (range[0] + range[1]) / 2.0;
|
|
|
174
|
|
|
175 vtkNew<vtkMarchingCubes> surface;
|
|
|
176 surface->SetInputData(padding->GetOutput());
|
|
|
177 surface->ComputeNormalsOn();
|
|
|
178 surface->SetValue(0, isoValue);
|
|
|
179 surface->Update();
|
|
|
180
|
|
|
181 if (smooth)
|
|
|
182 {
|
|
|
183 vtkNew<vtkSmoothPolyDataFilter> smoothFilter;
|
|
|
184 // Apply volume smoothing
|
|
|
185 // https://examples.vtk.org/site/Cxx/PolyData/SmoothPolyDataFilter/
|
|
|
186 smoothFilter->SetInputConnection(surface->GetOutputPort());
|
|
|
187 smoothFilter->SetNumberOfIterations(15);
|
|
|
188 smoothFilter->SetRelaxationFactor(0.1);
|
|
|
189 smoothFilter->FeatureEdgeSmoothingOff();
|
|
|
190 smoothFilter->BoundarySmoothingOn();
|
|
|
191 smoothFilter->Update();
|
|
|
192
|
|
|
193 vtkNew<vtkPolyDataNormals> normalGenerator;
|
|
|
194 normalGenerator->SetInputConnection(smoothFilter->GetOutputPort());
|
|
|
195 normalGenerator->ComputePointNormalsOn();
|
|
|
196 normalGenerator->ComputeCellNormalsOn();
|
|
|
197 normalGenerator->Update();
|
|
|
198
|
|
|
199 EncodeSTL(stl, *normalGenerator->GetOutput());
|
|
|
200 }
|
|
|
201 else
|
|
|
202 {
|
|
|
203 EncodeSTL(stl, *surface->GetOutput());
|
|
|
204 }
|
|
|
205
|
|
|
206 return true;
|
|
|
207 }
|
|
|
208
|
|
|
209
|
|
|
210 namespace
|
|
|
211 {
|
|
|
212 class NiftiHeader : public boost::noncopyable
|
|
|
213 {
|
|
|
214 private:
|
|
|
215 nifti_image* image_;
|
|
|
216
|
|
|
217 public:
|
|
|
218 NiftiHeader(const std::string& nifti)
|
|
|
219 {
|
|
|
220 nifti_1_header header;
|
|
|
221 if (nifti.size() < sizeof(header))
|
|
|
222 {
|
|
|
223 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
|
|
|
224 }
|
|
|
225
|
|
|
226 memcpy(&header, nifti.c_str(), sizeof(header));
|
|
|
227 if (!nifti_hdr_looks_good(&header))
|
|
|
228 {
|
|
|
229 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
|
|
|
230 }
|
|
|
231
|
|
|
232 image_ = nifti_convert_nhdr2nim(header, "dummy_filename");
|
|
|
233 if (image_ == NULL)
|
|
|
234 {
|
|
|
235 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
|
|
|
236 }
|
|
|
237 }
|
|
|
238
|
|
|
239 ~NiftiHeader()
|
|
|
240 {
|
|
|
241 nifti_image_free(image_);
|
|
|
242 }
|
|
|
243
|
|
|
244 const nifti_image& GetInfo() const
|
|
|
245 {
|
|
|
246 assert(image_ != NULL);
|
|
|
247 return *image_;
|
|
|
248 }
|
|
|
249 };
|
|
|
250 }
|
|
|
251
|
|
|
252
|
|
|
253 void LoadNifti(vtkImageData* volume /* out */,
|
|
|
254 std::string& nifti /* in */)
|
|
|
255 {
|
|
|
256 if (volume == NULL)
|
|
|
257 {
|
|
|
258 throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer);
|
|
|
259 }
|
|
|
260
|
|
|
261 const uint8_t* p = reinterpret_cast<const uint8_t*>(nifti.c_str());
|
|
|
262
|
|
|
263 if (nifti.size() >= 2 &&
|
|
|
264 p[0] == 0x1f &&
|
|
|
265 p[1] == 0x8b)
|
|
|
266 {
|
|
|
267 Orthanc::GzipCompressor compressor;
|
|
|
268 std::string uncompressed;
|
|
|
269 Orthanc::IBufferCompressor::Uncompress(uncompressed, compressor, nifti);
|
|
|
270 nifti.swap(uncompressed);
|
|
|
271 }
|
|
|
272
|
|
|
273 NiftiHeader header(nifti);
|
|
|
274
|
|
|
275 if (header.GetInfo().ndim != 3)
|
|
|
276 {
|
|
|
277 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat,
|
|
|
278 "Only 3D NIfTI volumes are allowed");
|
|
|
279 }
|
|
|
280
|
|
|
281 size_t itemSize;
|
|
|
282 int vtkType;
|
|
|
283
|
|
|
284 switch (header.GetInfo().datatype)
|
|
|
285 {
|
|
|
286 case DT_UNSIGNED_CHAR:
|
|
|
287 itemSize = 1;
|
|
|
288 vtkType = VTK_UNSIGNED_CHAR;
|
|
|
289 break;
|
|
|
290
|
|
|
291 case DT_FLOAT:
|
|
|
292 itemSize = sizeof(float);
|
|
|
293 vtkType = VTK_FLOAT;
|
|
|
294 break;
|
|
|
295
|
|
|
296 case DT_DOUBLE:
|
|
|
297 itemSize = sizeof(double);
|
|
|
298 vtkType = VTK_DOUBLE;
|
|
|
299 break;
|
|
|
300
|
|
|
301 default:
|
|
|
302 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
|
|
|
303 }
|
|
|
304
|
|
|
305 assert(static_cast<int>(header.GetInfo().nvox) == header.GetInfo().nx * header.GetInfo().ny * header.GetInfo().nz);
|
|
|
306
|
|
|
307 const size_t pixelDataOffset = sizeof(nifti_1_header) + 4 /* extension */;
|
|
|
308
|
|
|
309 if (nifti.size() != pixelDataOffset + header.GetInfo().nvox * itemSize)
|
|
|
310 {
|
|
|
311 throw Orthanc::OrthancException(Orthanc::ErrorCode_CorruptedFile);
|
|
|
312 }
|
|
|
313
|
|
|
314 volume->SetDimensions(header.GetInfo().nx, header.GetInfo().ny, header.GetInfo().nz);
|
|
|
315 volume->AllocateScalars(vtkType, 1);
|
|
|
316 volume->SetSpacing(header.GetInfo().dx, header.GetInfo().dy, header.GetInfo().dz);
|
|
|
317 memcpy(volume->GetScalarPointer(), &nifti[pixelDataOffset], header.GetInfo().nvox * itemSize);
|
|
|
318 }
|
|
|
319 }
|
