Mercurial > hg > orthanc-stl
comparison Sources/VTKToolbox.cpp @ 32:976da5476810
reorganization
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Thu, 04 Apr 2024 18:35:54 +0200 |
parents | |
children | 2460b376d3f7 |
comparison
equal
deleted
inserted
replaced
31:ab231760799d | 32:976da5476810 |
---|---|
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 } |