Mercurial > hg > orthanc-stl

comparison Sources/VTKToolbox.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
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 }