Mercurial > hg > orthanc-stone
comparison Framework/Volumes/VolumeReslicer.cpp @ 154:7a52c968ea1b wasm
reorganization
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Fri, 02 Feb 2018 09:39:13 +0100 |
parents | Framework/Toolbox/VolumeReslicer.cpp@ae531ab5dcd9 |
children | a053ca7fa5c6 |
comparison
equal
deleted
inserted
replaced
153:ae531ab5dcd9 | 154:7a52c968ea1b |
---|---|
1 #include "VolumeReslicer.h" | |
2 | |
3 #include <Core/Logging.h> | |
4 #include <Core/OrthancException.h> | |
5 | |
6 #include <boost/math/special_functions/round.hpp> | |
7 | |
8 #if defined(_MSC_VER) | |
9 # define ORTHANC_STONE_FORCE_INLINE __forceinline | |
10 #elif defined(__GNUC__) || defined(__clang__) || defined(__EMSCRIPTEN__) | |
11 # define ORTHANC_STONE_FORCE_INLINE inline __attribute((always_inline)) | |
12 #else | |
13 # error Please support your compiler here | |
14 #endif | |
15 | |
16 | |
17 namespace OrthancStone | |
18 { | |
19 // Anonymous namespace to avoid clashes between compilation modules | |
20 namespace | |
21 { | |
22 enum TransferFunction | |
23 { | |
24 TransferFunction_Copy, | |
25 TransferFunction_Float, | |
26 TransferFunction_Linear | |
27 }; | |
28 | |
29 | |
30 template <Orthanc::PixelFormat Format> | |
31 struct PixelTraits; | |
32 | |
33 template <> | |
34 struct PixelTraits<Orthanc::PixelFormat_Grayscale8> | |
35 { | |
36 typedef uint8_t PixelType; | |
37 | |
38 static void SetOutOfVolume(PixelType& target) | |
39 { | |
40 target = 0; | |
41 } | |
42 | |
43 static void GetVoxel(PixelType& target, | |
44 const ImageBuffer3D& image, | |
45 unsigned int x, | |
46 unsigned int y, | |
47 unsigned int z) | |
48 { | |
49 assert(x < image.GetWidth() && y < image.GetHeight() && z < image.GetDepth()); | |
50 target = image.GetVoxelGrayscale8Unchecked(x, y, z); | |
51 } | |
52 | |
53 static float GetFloatVoxel(const ImageBuffer3D& image, | |
54 unsigned int x, | |
55 unsigned int y, | |
56 unsigned int z) | |
57 { | |
58 assert(x < image.GetWidth() && y < image.GetHeight() && z < image.GetDepth()); | |
59 return static_cast<float>(image.GetVoxelGrayscale8Unchecked(x, y, z)); | |
60 } | |
61 }; | |
62 | |
63 template <> | |
64 struct PixelTraits<Orthanc::PixelFormat_Grayscale16> | |
65 { | |
66 typedef uint16_t PixelType; | |
67 | |
68 static void SetOutOfVolume(PixelType& target) | |
69 { | |
70 target = 0; | |
71 } | |
72 | |
73 static void GetVoxel(PixelType& target, | |
74 const ImageBuffer3D& image, | |
75 unsigned int x, | |
76 unsigned int y, | |
77 unsigned int z) | |
78 { | |
79 assert(x < image.GetWidth() && y < image.GetHeight() && z < image.GetDepth()); | |
80 target = image.GetVoxelGrayscale16Unchecked(x, y, z); | |
81 } | |
82 | |
83 static float GetFloatVoxel(const ImageBuffer3D& image, | |
84 unsigned int x, | |
85 unsigned int y, | |
86 unsigned int z) | |
87 { | |
88 assert(x < image.GetWidth() && y < image.GetHeight() && z < image.GetDepth()); | |
89 return static_cast<float>(image.GetVoxelGrayscale16Unchecked(x, y, z)); | |
90 } | |
91 }; | |
92 | |
93 | |
94 template <> | |
95 struct PixelTraits<Orthanc::PixelFormat_SignedGrayscale16> | |
96 { | |
97 typedef int16_t PixelType; | |
98 | |
99 static void SetOutOfVolume(PixelType& target) | |
100 { | |
101 target = std::numeric_limits<PixelType>::min(); | |
102 } | |
103 | |
104 static void GetVoxel(PixelType& target, | |
105 const ImageBuffer3D& image, | |
106 unsigned int x, | |
107 unsigned int y, | |
108 unsigned int z) | |
109 { | |
110 assert(x < image.GetWidth() && y < image.GetHeight() && z < image.GetDepth()); | |
111 target = image.GetVoxelSignedGrayscale16Unchecked(x, y, z); | |
112 } | |
113 | |
114 static float GetFloatVoxel(const ImageBuffer3D& image, | |
115 unsigned int x, | |
116 unsigned int y, | |
117 unsigned int z) | |
118 { | |
119 assert(x < image.GetWidth() && y < image.GetHeight() && z < image.GetDepth()); | |
120 return static_cast<float>(image.GetVoxelSignedGrayscale16Unchecked(x, y, z)); | |
121 } | |
122 }; | |
123 | |
124 | |
125 template <> | |
126 struct PixelTraits<Orthanc::PixelFormat_BGRA32> | |
127 { | |
128 struct PixelType | |
129 { | |
130 uint8_t blue_; | |
131 uint8_t green_; | |
132 uint8_t red_; | |
133 uint8_t alpha_; | |
134 }; | |
135 }; | |
136 | |
137 | |
138 | |
139 template <Orthanc::PixelFormat InputFormat, | |
140 Orthanc::PixelFormat OutputFormat> | |
141 class PixelWriter | |
142 { | |
143 public: | |
144 typedef typename PixelTraits<InputFormat>::PixelType InputPixelType; | |
145 typedef PixelTraits<OutputFormat> OutputPixelTraits; | |
146 typedef typename PixelTraits<OutputFormat>::PixelType OutputPixelType; | |
147 | |
148 private: | |
149 template <typename T> | |
150 static void SetValueInternal(OutputPixelType* pixel, | |
151 const T& value) | |
152 { | |
153 if (value < std::numeric_limits<OutputPixelType>::min()) | |
154 { | |
155 *pixel = std::numeric_limits<OutputPixelType>::min(); | |
156 } | |
157 else if (value > std::numeric_limits<OutputPixelType>::max()) | |
158 { | |
159 *pixel = std::numeric_limits<OutputPixelType>::max(); | |
160 } | |
161 else | |
162 { | |
163 *pixel = static_cast<OutputPixelType>(value); | |
164 } | |
165 } | |
166 | |
167 public: | |
168 ORTHANC_STONE_FORCE_INLINE | |
169 void SetFloatValue(OutputPixelType* pixel, | |
170 float value) const | |
171 { | |
172 SetValueInternal<float>(pixel, value); | |
173 } | |
174 | |
175 ORTHANC_STONE_FORCE_INLINE | |
176 void SetValue(OutputPixelType* pixel, | |
177 const InputPixelType& value) const | |
178 { | |
179 SetValueInternal<InputPixelType>(pixel, value); | |
180 } | |
181 }; | |
182 | |
183 | |
184 template <Orthanc::PixelFormat InputFormat> | |
185 class PixelWriter<InputFormat, Orthanc::PixelFormat_BGRA32> | |
186 { | |
187 public: | |
188 typedef typename PixelTraits<InputFormat>::PixelType InputPixelType; | |
189 typedef PixelTraits<Orthanc::PixelFormat_BGRA32> OutputPixelTraits; | |
190 typedef PixelTraits<Orthanc::PixelFormat_BGRA32>::PixelType OutputPixelType; | |
191 | |
192 private: | |
193 template <typename T> | |
194 static void SetValueInternal(OutputPixelType* pixel, | |
195 const T& value) | |
196 { | |
197 uint8_t v; | |
198 if (value < 0) | |
199 { | |
200 v = 0; | |
201 } | |
202 else if (value >= 255.0f) | |
203 { | |
204 v = 255; | |
205 } | |
206 else | |
207 { | |
208 v = static_cast<uint8_t>(value); | |
209 } | |
210 | |
211 pixel->blue_ = v; | |
212 pixel->green_ = v; | |
213 pixel->red_ = v; | |
214 pixel->alpha_ = 255; | |
215 } | |
216 | |
217 public: | |
218 ORTHANC_STONE_FORCE_INLINE | |
219 void SetFloatValue(OutputPixelType* pixel, | |
220 float value) const | |
221 { | |
222 SetValueInternal<float>(pixel, value); | |
223 } | |
224 | |
225 ORTHANC_STONE_FORCE_INLINE | |
226 void SetValue(OutputPixelType* pixel, | |
227 const InputPixelType& value) const | |
228 { | |
229 SetValueInternal<InputPixelType>(pixel, value); | |
230 } | |
231 }; | |
232 | |
233 | |
234 | |
235 class VoxelReaderBase | |
236 { | |
237 protected: | |
238 const ImageBuffer3D& source_; | |
239 unsigned int sourceWidth_; | |
240 unsigned int sourceHeight_; | |
241 unsigned int sourceDepth_; | |
242 | |
243 public: | |
244 VoxelReaderBase(const ImageBuffer3D& source) : | |
245 source_(source), | |
246 sourceWidth_(source.GetWidth()), | |
247 sourceHeight_(source.GetHeight()), | |
248 sourceDepth_(source.GetDepth()) | |
249 { | |
250 } | |
251 | |
252 unsigned int GetSourceWidth() const | |
253 { | |
254 return sourceWidth_; | |
255 } | |
256 | |
257 unsigned int GetSourceHeight() const | |
258 { | |
259 return sourceHeight_; | |
260 } | |
261 | |
262 unsigned int GetSourceDepth() const | |
263 { | |
264 return sourceDepth_; | |
265 } | |
266 | |
267 bool GetNearestCoordinates(unsigned int& sourceX, | |
268 unsigned int& sourceY, | |
269 unsigned int& sourceZ, | |
270 float worldX, | |
271 float worldY, | |
272 float worldZ) const | |
273 { | |
274 if (worldX >= 0 && | |
275 worldY >= 0 && | |
276 worldZ >= 0) | |
277 { | |
278 sourceX = static_cast<unsigned int>(worldX * static_cast<float>(sourceWidth_)); | |
279 sourceY = static_cast<unsigned int>(worldY * static_cast<float>(sourceHeight_)); | |
280 sourceZ = static_cast<unsigned int>(worldZ * static_cast<float>(sourceDepth_)); | |
281 | |
282 return (sourceX < sourceWidth_ && | |
283 sourceY < sourceHeight_ && | |
284 sourceZ < sourceDepth_); | |
285 } | |
286 else | |
287 { | |
288 return false; | |
289 } | |
290 } | |
291 }; | |
292 | |
293 | |
294 template <Orthanc::PixelFormat InputFormat, | |
295 ImageInterpolation Interpolation> | |
296 class VoxelReader; | |
297 | |
298 | |
299 template <Orthanc::PixelFormat InputFormat> | |
300 class VoxelReader<InputFormat, ImageInterpolation_Nearest> : | |
301 public VoxelReaderBase | |
302 { | |
303 public: | |
304 typedef typename PixelTraits<InputFormat>::PixelType InputPixelType; | |
305 | |
306 VoxelReader(const ImageBuffer3D& source) : | |
307 VoxelReaderBase(source) | |
308 { | |
309 } | |
310 | |
311 ORTHANC_STONE_FORCE_INLINE | |
312 float GetFloatValue(float worldX, | |
313 float worldY, | |
314 float worldZ) const | |
315 { | |
316 InputPixelType value; | |
317 GetValue(value, worldX, worldY, worldZ); | |
318 return static_cast<float>(value); | |
319 } | |
320 | |
321 ORTHANC_STONE_FORCE_INLINE | |
322 void GetValue(InputPixelType& target, | |
323 float worldX, | |
324 float worldY, | |
325 float worldZ) const | |
326 { | |
327 unsigned int sourceX, sourceY, sourceZ; | |
328 | |
329 if (GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ)) | |
330 { | |
331 PixelTraits<InputFormat>::GetVoxel(target, source_, sourceX, sourceY, sourceZ); | |
332 } | |
333 else | |
334 { | |
335 PixelTraits<InputFormat>::SetOutOfVolume(target); | |
336 } | |
337 } | |
338 }; | |
339 | |
340 | |
341 template <Orthanc::PixelFormat InputFormat> | |
342 class VoxelReader<InputFormat, ImageInterpolation_Bilinear> : | |
343 public VoxelReaderBase | |
344 { | |
345 private: | |
346 float outOfVolume_; | |
347 | |
348 public: | |
349 VoxelReader(const ImageBuffer3D& source) : | |
350 VoxelReaderBase(source) | |
351 { | |
352 typename PixelTraits<InputFormat>::PixelType value; | |
353 PixelTraits<InputFormat>::SetOutOfVolume(value); | |
354 outOfVolume_ = static_cast<float>(value); | |
355 } | |
356 | |
357 void SampleVoxels(float& f00, | |
358 float& f10, | |
359 float& f01, | |
360 float& f11, | |
361 unsigned int sourceX, | |
362 unsigned int sourceY, | |
363 unsigned int sourceZ) const | |
364 { | |
365 f00 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX, sourceY, sourceZ); | |
366 | |
367 if (sourceX + 1 < sourceWidth_) | |
368 { | |
369 f01 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY, sourceZ); | |
370 } | |
371 else | |
372 { | |
373 f01 = f00; | |
374 } | |
375 | |
376 if (sourceY + 1 < sourceHeight_) | |
377 { | |
378 f10 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX, sourceY + 1, sourceZ); | |
379 } | |
380 else | |
381 { | |
382 f10 = f00; | |
383 } | |
384 | |
385 if (sourceX + 1 < sourceWidth_ && | |
386 sourceY + 1 < sourceHeight_) | |
387 { | |
388 f11 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY + 1, sourceZ); | |
389 } | |
390 else | |
391 { | |
392 f11 = f00; | |
393 } | |
394 } | |
395 | |
396 float GetOutOfVolume() const | |
397 { | |
398 return outOfVolume_; | |
399 } | |
400 | |
401 float GetFloatValue(float worldX, | |
402 float worldY, | |
403 float worldZ) const | |
404 { | |
405 unsigned int sourceX, sourceY, sourceZ; | |
406 | |
407 if (GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ)) | |
408 { | |
409 float f00, f10, f01, f11; | |
410 SampleVoxels(f00, f10, f01, f11, sourceX, sourceY, sourceZ); | |
411 return GeometryToolbox::ComputeBilinearInterpolation(worldX, worldY, f00, f01, f10, f11); | |
412 } | |
413 else | |
414 { | |
415 return outOfVolume_; | |
416 } | |
417 } | |
418 }; | |
419 | |
420 | |
421 template <Orthanc::PixelFormat InputFormat> | |
422 class VoxelReader<InputFormat, ImageInterpolation_Trilinear> | |
423 { | |
424 private: | |
425 typedef VoxelReader<InputFormat, ImageInterpolation_Bilinear> Bilinear; | |
426 | |
427 Bilinear bilinear_; | |
428 unsigned int sourceDepth_; | |
429 | |
430 public: | |
431 VoxelReader(const ImageBuffer3D& source) : | |
432 bilinear_(source), | |
433 sourceDepth_(source.GetDepth()) | |
434 { | |
435 } | |
436 | |
437 float GetFloatValue(float worldX, | |
438 float worldY, | |
439 float worldZ) const | |
440 { | |
441 unsigned int sourceX, sourceY, sourceZ; | |
442 | |
443 if (bilinear_.GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ)) | |
444 { | |
445 float f000, f010, f001, f011, f100, f110, f101, f111; | |
446 | |
447 bilinear_.SampleVoxels(f000, f010, f001, f011, sourceX, sourceY, sourceZ); | |
448 | |
449 if (sourceZ + 1 < sourceDepth_) | |
450 { | |
451 bilinear_.SampleVoxels(f100, f110, f101, f111, sourceX, sourceY, sourceZ + 1); | |
452 return GeometryToolbox::ComputeTrilinearInterpolation | |
453 (worldX, worldY, worldZ, f000, f001, f010, f011, f100, f110, f101, f111); | |
454 } | |
455 else | |
456 { | |
457 return GeometryToolbox::ComputeBilinearInterpolation(worldX, worldY, f000, f001, f010, f011); | |
458 } | |
459 } | |
460 else | |
461 { | |
462 return bilinear_.GetOutOfVolume(); | |
463 } | |
464 } | |
465 }; | |
466 | |
467 | |
468 template <typename VoxelReader, | |
469 typename PixelWriter, | |
470 TransferFunction Function> | |
471 class PixelShader; | |
472 | |
473 | |
474 template <typename VoxelReader, | |
475 typename PixelWriter> | |
476 class PixelShader<VoxelReader, PixelWriter, TransferFunction_Copy> | |
477 { | |
478 private: | |
479 VoxelReader reader_; | |
480 PixelWriter writer_; | |
481 | |
482 public: | |
483 PixelShader(const ImageBuffer3D& image, | |
484 float /* scaling */, | |
485 float /* offset */) : | |
486 reader_(image) | |
487 { | |
488 } | |
489 | |
490 ORTHANC_STONE_FORCE_INLINE | |
491 void Apply(typename PixelWriter::OutputPixelType* pixel, | |
492 float worldX, | |
493 float worldY, | |
494 float worldZ) | |
495 { | |
496 typename VoxelReader::InputPixelType source; | |
497 reader_.GetValue(source, worldX, worldY, worldZ); | |
498 writer_.SetValue(pixel, source); | |
499 } | |
500 }; | |
501 | |
502 | |
503 template <typename VoxelReader, | |
504 typename PixelWriter> | |
505 class PixelShader<VoxelReader, PixelWriter, TransferFunction_Float> | |
506 { | |
507 private: | |
508 VoxelReader reader_; | |
509 PixelWriter writer_; | |
510 | |
511 public: | |
512 PixelShader(const ImageBuffer3D& image, | |
513 float /* scaling */, | |
514 float /* offset */) : | |
515 reader_(image) | |
516 { | |
517 } | |
518 | |
519 ORTHANC_STONE_FORCE_INLINE | |
520 void Apply(typename PixelWriter::OutputPixelType* pixel, | |
521 float worldX, | |
522 float worldY, | |
523 float worldZ) | |
524 { | |
525 writer_.SetFloatValue(pixel, reader_.GetFloatValue(worldX, worldY, worldZ)); | |
526 } | |
527 }; | |
528 | |
529 | |
530 template <typename VoxelReader, | |
531 typename PixelWriter> | |
532 class PixelShader<VoxelReader, PixelWriter, TransferFunction_Linear> | |
533 { | |
534 private: | |
535 VoxelReader reader_; | |
536 PixelWriter writer_; | |
537 float scaling_; | |
538 float offset_; | |
539 | |
540 public: | |
541 PixelShader(const ImageBuffer3D& image, | |
542 float scaling, | |
543 float offset) : | |
544 reader_(image), | |
545 scaling_(scaling), | |
546 offset_(offset) | |
547 { | |
548 } | |
549 | |
550 ORTHANC_STONE_FORCE_INLINE | |
551 void Apply(typename PixelWriter::OutputPixelType* pixel, | |
552 float worldX, | |
553 float worldY, | |
554 float worldZ) | |
555 { | |
556 writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(worldX, worldY, worldZ) + offset_); | |
557 } | |
558 }; | |
559 | |
560 | |
561 | |
562 class FastRowIterator : public boost::noncopyable | |
563 { | |
564 private: | |
565 float position_[3]; | |
566 float offset_[3]; | |
567 | |
568 public: | |
569 FastRowIterator(const Orthanc::ImageAccessor& slice, | |
570 const Extent2D& extent, | |
571 const CoordinateSystem3D& plane, | |
572 const OrientedBoundingBox& box, | |
573 unsigned int y) | |
574 { | |
575 const double width = static_cast<double>(slice.GetWidth()); | |
576 const double height = static_cast<double>(slice.GetHeight()); | |
577 assert(y < height); | |
578 | |
579 Vector q1 = plane.MapSliceToWorldCoordinates | |
580 (extent.GetX1() + extent.GetWidth() * static_cast<double>(0) / static_cast<double>(width + 1), | |
581 extent.GetY1() + extent.GetHeight() * static_cast<double>(y) / static_cast<double>(height + 1)); | |
582 | |
583 Vector q2 = plane.MapSliceToWorldCoordinates | |
584 (extent.GetX1() + extent.GetWidth() * static_cast<double>(width - 1) / static_cast<double>(width + 1), | |
585 extent.GetY1() + extent.GetHeight() * static_cast<double>(y) / static_cast<double>(height + 1)); | |
586 | |
587 Vector r1, r2; | |
588 box.ToInternalCoordinates(r1, q1); | |
589 box.ToInternalCoordinates(r2, q2); | |
590 | |
591 position_[0] = static_cast<float>(r1[0]); | |
592 position_[1] = static_cast<float>(r1[1]); | |
593 position_[2] = static_cast<float>(r1[2]); | |
594 | |
595 Vector tmp = (r2 - r1) / static_cast<double>(width - 1); | |
596 offset_[0] = static_cast<float>(tmp[0]); | |
597 offset_[1] = static_cast<float>(tmp[1]); | |
598 offset_[2] = static_cast<float>(tmp[2]); | |
599 } | |
600 | |
601 ORTHANC_STONE_FORCE_INLINE | |
602 void Next() | |
603 { | |
604 position_[0] += offset_[0]; | |
605 position_[1] += offset_[1]; | |
606 position_[2] += offset_[2]; | |
607 } | |
608 | |
609 ORTHANC_STONE_FORCE_INLINE | |
610 void GetNearestCoordinates(float& x, | |
611 float& y, | |
612 float& z) const | |
613 { | |
614 x = position_[0]; | |
615 y = position_[1]; | |
616 z = position_[2]; | |
617 } | |
618 }; | |
619 | |
620 | |
621 class SlowRowIterator : public boost::noncopyable | |
622 { | |
623 private: | |
624 const Orthanc::ImageAccessor& slice_; | |
625 const Extent2D& extent_; | |
626 const CoordinateSystem3D& plane_; | |
627 const OrientedBoundingBox& box_; | |
628 unsigned int x_; | |
629 unsigned int y_; | |
630 | |
631 public: | |
632 SlowRowIterator(const Orthanc::ImageAccessor& slice, | |
633 const Extent2D& extent, | |
634 const CoordinateSystem3D& plane, | |
635 const OrientedBoundingBox& box, | |
636 unsigned int y) : | |
637 slice_(slice), | |
638 extent_(extent), | |
639 plane_(plane), | |
640 box_(box), | |
641 x_(0), | |
642 y_(y) | |
643 { | |
644 assert(y_ < slice_.GetHeight()); | |
645 } | |
646 | |
647 void Next() | |
648 { | |
649 x_++; | |
650 } | |
651 | |
652 void GetNearestCoordinates(float& x, | |
653 float& y, | |
654 float& z) const | |
655 { | |
656 assert(x_ < slice_.GetWidth()); | |
657 | |
658 const double width = static_cast<double>(slice_.GetWidth()); | |
659 const double height = static_cast<double>(slice_.GetHeight()); | |
660 | |
661 Vector q = plane_.MapSliceToWorldCoordinates | |
662 (extent_.GetX1() + extent_.GetWidth() * static_cast<double>(x_) / (width + 1.0), | |
663 extent_.GetY1() + extent_.GetHeight() * static_cast<double>(y_) / (height + 1.0)); | |
664 | |
665 Vector r; | |
666 box_.ToInternalCoordinates(r, q); | |
667 | |
668 x = static_cast<float>(r[0]); | |
669 y = static_cast<float>(r[1]); | |
670 z = static_cast<float>(r[2]); | |
671 } | |
672 }; | |
673 | |
674 | |
675 template <typename RowIterator, | |
676 Orthanc::PixelFormat InputFormat, | |
677 Orthanc::PixelFormat OutputFormat, | |
678 ImageInterpolation Interpolation, | |
679 TransferFunction Function> | |
680 static void ProcessImage(Orthanc::ImageAccessor& slice, | |
681 const Extent2D& extent, | |
682 const ImageBuffer3D& source, | |
683 const CoordinateSystem3D& plane, | |
684 const OrientedBoundingBox& box, | |
685 float scaling, | |
686 float offset) | |
687 { | |
688 typedef VoxelReader<InputFormat, Interpolation> Reader; | |
689 typedef PixelWriter<InputFormat, OutputFormat> Writer; | |
690 typedef PixelShader<Reader, Writer, Function> Shader; | |
691 | |
692 const unsigned int outputWidth = slice.GetWidth(); | |
693 const unsigned int outputHeight = slice.GetHeight(); | |
694 | |
695 Shader shader(source, scaling, offset); | |
696 | |
697 for (unsigned int y = 0; y < outputHeight; y++) | |
698 { | |
699 typename Writer::OutputPixelType* p = reinterpret_cast<typename Writer::OutputPixelType*>(slice.GetRow(y)); | |
700 | |
701 RowIterator it(slice, extent, plane, box, y); | |
702 | |
703 for (unsigned int x = 0; x < outputWidth; x++, p++) | |
704 { | |
705 float worldX, worldY, worldZ; | |
706 it.GetNearestCoordinates(worldX, worldY, worldZ); | |
707 shader.Apply(p, worldX, worldY, worldZ); | |
708 it.Next(); | |
709 } | |
710 } | |
711 } | |
712 | |
713 | |
714 template <typename RowIterator, | |
715 Orthanc::PixelFormat InputFormat, | |
716 Orthanc::PixelFormat OutputFormat> | |
717 static void ProcessImage(Orthanc::ImageAccessor& slice, | |
718 const Extent2D& extent, | |
719 const ImageBuffer3D& source, | |
720 const CoordinateSystem3D& plane, | |
721 const OrientedBoundingBox& box, | |
722 ImageInterpolation interpolation, | |
723 bool hasLinearFunction, | |
724 float scaling, | |
725 float offset) | |
726 { | |
727 if (hasLinearFunction) | |
728 { | |
729 switch (interpolation) | |
730 { | |
731 case ImageInterpolation_Nearest: | |
732 ProcessImage<RowIterator, InputFormat, OutputFormat, | |
733 ImageInterpolation_Nearest, TransferFunction_Linear> | |
734 (slice, extent, source, plane, box, scaling, offset); | |
735 break; | |
736 | |
737 case ImageInterpolation_Bilinear: | |
738 ProcessImage<RowIterator, InputFormat, OutputFormat, | |
739 ImageInterpolation_Bilinear, TransferFunction_Linear> | |
740 (slice, extent, source, plane, box, scaling, offset); | |
741 break; | |
742 | |
743 case ImageInterpolation_Trilinear: | |
744 ProcessImage<RowIterator, InputFormat, OutputFormat, | |
745 ImageInterpolation_Trilinear, TransferFunction_Linear> | |
746 (slice, extent, source, plane, box, scaling, offset); | |
747 break; | |
748 | |
749 default: | |
750 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
751 } | |
752 } | |
753 else | |
754 { | |
755 switch (interpolation) | |
756 { | |
757 case ImageInterpolation_Nearest: | |
758 ProcessImage<RowIterator, InputFormat, OutputFormat, | |
759 ImageInterpolation_Nearest, TransferFunction_Copy> | |
760 (slice, extent, source, plane, box, 0, 0); | |
761 break; | |
762 | |
763 case ImageInterpolation_Bilinear: | |
764 ProcessImage<RowIterator, InputFormat, OutputFormat, | |
765 ImageInterpolation_Bilinear, TransferFunction_Float> | |
766 (slice, extent, source, plane, box, 0, 0); | |
767 break; | |
768 | |
769 case ImageInterpolation_Trilinear: | |
770 ProcessImage<RowIterator, InputFormat, OutputFormat, | |
771 ImageInterpolation_Trilinear, TransferFunction_Float> | |
772 (slice, extent, source, plane, box, 0, 0); | |
773 break; | |
774 | |
775 default: | |
776 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
777 } | |
778 } | |
779 } | |
780 | |
781 | |
782 template <typename RowIterator> | |
783 static void ProcessImage(Orthanc::ImageAccessor& slice, | |
784 const Extent2D& extent, | |
785 const ImageBuffer3D& source, | |
786 const CoordinateSystem3D& plane, | |
787 const OrientedBoundingBox& box, | |
788 ImageInterpolation interpolation, | |
789 bool hasLinearFunction, | |
790 float scaling, | |
791 float offset) | |
792 { | |
793 if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && | |
794 slice.GetFormat() == Orthanc::PixelFormat_Grayscale8) | |
795 { | |
796 ProcessImage<RowIterator, | |
797 Orthanc::PixelFormat_Grayscale16, | |
798 Orthanc::PixelFormat_Grayscale8> | |
799 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); | |
800 } | |
801 else if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && | |
802 slice.GetFormat() == Orthanc::PixelFormat_Grayscale16) | |
803 { | |
804 ProcessImage<RowIterator, | |
805 Orthanc::PixelFormat_Grayscale16, | |
806 Orthanc::PixelFormat_Grayscale16> | |
807 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); | |
808 } | |
809 else if (source.GetFormat() == Orthanc::PixelFormat_SignedGrayscale16 && | |
810 slice.GetFormat() == Orthanc::PixelFormat_BGRA32) | |
811 { | |
812 ProcessImage<RowIterator, | |
813 Orthanc::PixelFormat_SignedGrayscale16, | |
814 Orthanc::PixelFormat_BGRA32> | |
815 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); | |
816 } | |
817 else if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && | |
818 slice.GetFormat() == Orthanc::PixelFormat_BGRA32) | |
819 { | |
820 ProcessImage<RowIterator, | |
821 Orthanc::PixelFormat_Grayscale16, | |
822 Orthanc::PixelFormat_BGRA32> | |
823 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); | |
824 } | |
825 else | |
826 { | |
827 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
828 } | |
829 } | |
830 } | |
831 | |
832 | |
833 | |
834 void VolumeReslicer::CheckIterators(const ImageBuffer3D& source, | |
835 const CoordinateSystem3D& plane, | |
836 const OrientedBoundingBox& box) const | |
837 { | |
838 for (unsigned int y = 0; y < slice_->GetHeight(); y++) | |
839 { | |
840 FastRowIterator fast(*slice_, extent_, plane, box, y); | |
841 SlowRowIterator slow(*slice_, extent_, plane, box, y); | |
842 | |
843 for (unsigned int x = 0; x < slice_->GetWidth(); x++) | |
844 { | |
845 float px, py, pz; | |
846 fast.GetNearestCoordinates(px, py, pz); | |
847 | |
848 float qx, qy, qz; | |
849 slow.GetNearestCoordinates(qx, qy, qz); | |
850 | |
851 Vector d; | |
852 GeometryToolbox::AssignVector(d, px - qx, py - qy, pz - qz); | |
853 double norm = boost::numeric::ublas::norm_2(d); | |
854 if (norm > 0.0001) | |
855 { | |
856 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
857 } | |
858 | |
859 fast.Next(); | |
860 slow.Next(); | |
861 } | |
862 } | |
863 } | |
864 | |
865 | |
866 void VolumeReslicer::Reset() | |
867 { | |
868 success_ = false; | |
869 extent_.Reset(); | |
870 slice_.reset(NULL); | |
871 } | |
872 | |
873 | |
874 float VolumeReslicer::GetMinOutputValue() const | |
875 { | |
876 switch (outputFormat_) | |
877 { | |
878 case Orthanc::PixelFormat_Grayscale8: | |
879 case Orthanc::PixelFormat_Grayscale16: | |
880 case Orthanc::PixelFormat_BGRA32: | |
881 return 0.0f; | |
882 break; | |
883 | |
884 default: | |
885 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
886 } | |
887 } | |
888 | |
889 | |
890 float VolumeReslicer::GetMaxOutputValue() const | |
891 { | |
892 switch (outputFormat_) | |
893 { | |
894 case Orthanc::PixelFormat_Grayscale8: | |
895 case Orthanc::PixelFormat_BGRA32: | |
896 return static_cast<float>(std::numeric_limits<uint8_t>::max()); | |
897 break; | |
898 | |
899 case Orthanc::PixelFormat_Grayscale16: | |
900 return static_cast<float>(std::numeric_limits<uint16_t>::max()); | |
901 break; | |
902 | |
903 default: | |
904 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
905 } | |
906 } | |
907 | |
908 | |
909 VolumeReslicer::VolumeReslicer() : | |
910 outputFormat_(Orthanc::PixelFormat_Grayscale8), | |
911 interpolation_(ImageInterpolation_Nearest), | |
912 fastMode_(true), | |
913 success_(false) | |
914 { | |
915 ResetLinearFunction(); | |
916 } | |
917 | |
918 | |
919 void VolumeReslicer::GetLinearFunction(float& scaling, | |
920 float& offset) const | |
921 { | |
922 if (hasLinearFunction_) | |
923 { | |
924 scaling = scaling_; | |
925 offset = offset_; | |
926 } | |
927 else | |
928 { | |
929 scaling = 1.0f; | |
930 offset = 0.0f; | |
931 } | |
932 } | |
933 | |
934 | |
935 void VolumeReslicer::ResetLinearFunction() | |
936 { | |
937 Reset(); | |
938 hasLinearFunction_ = false; | |
939 scaling_ = 1.0f; | |
940 offset_ = 0.0f; | |
941 } | |
942 | |
943 | |
944 void VolumeReslicer::SetLinearFunction(float scaling, | |
945 float offset) | |
946 { | |
947 Reset(); | |
948 hasLinearFunction_ = true; | |
949 scaling_ = scaling; | |
950 offset_ = offset; | |
951 } | |
952 | |
953 | |
954 void VolumeReslicer::SetWindow(float low, | |
955 float high) | |
956 { | |
957 //printf("Range in pixel values: %f->%f\n", low, high); | |
958 float scaling = (GetMaxOutputValue() - GetMinOutputValue()) / (high - low); | |
959 float offset = GetMinOutputValue() - scaling * low; | |
960 | |
961 SetLinearFunction(scaling, offset); | |
962 | |
963 /*float x = scaling_ * low + offset_; | |
964 float y = scaling_ * high + offset_; | |
965 printf("%f %f (should be %f->%f)\n", x, y, GetMinOutputValue(), GetMaxOutputValue());*/ | |
966 } | |
967 | |
968 | |
969 void VolumeReslicer::FitRange(const ImageBuffer3D& image) | |
970 { | |
971 float minInputValue, maxInputValue; | |
972 | |
973 if (!image.GetRange(minInputValue, maxInputValue) || | |
974 maxInputValue < 1) | |
975 { | |
976 ResetLinearFunction(); | |
977 } | |
978 else | |
979 { | |
980 SetWindow(minInputValue, maxInputValue); | |
981 } | |
982 } | |
983 | |
984 | |
985 void VolumeReslicer::SetWindowing(ImageWindowing windowing, | |
986 const ImageBuffer3D& image, | |
987 float rescaleSlope, | |
988 float rescaleIntercept) | |
989 { | |
990 if (windowing == ImageWindowing_Custom || | |
991 windowing == ImageWindowing_Default) | |
992 { | |
993 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
994 } | |
995 | |
996 float center, width; | |
997 ComputeWindowing(center, width, windowing, 0, 0); | |
998 | |
999 float a = (center - width / 2.0f - rescaleIntercept) / rescaleSlope; | |
1000 float b = (center + width / 2.0f - rescaleIntercept) / rescaleSlope; | |
1001 SetWindow(a, b); | |
1002 } | |
1003 | |
1004 | |
1005 void VolumeReslicer::SetOutputFormat(Orthanc::PixelFormat format) | |
1006 { | |
1007 if (format != Orthanc::PixelFormat_Grayscale8 && | |
1008 format != Orthanc::PixelFormat_Grayscale16 && | |
1009 format != Orthanc::PixelFormat_BGRA32) | |
1010 { | |
1011 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
1012 } | |
1013 | |
1014 if (hasLinearFunction_) | |
1015 { | |
1016 LOG(WARNING) << "Calls to VolumeReslicer::SetOutputFormat() should be done before VolumeReslicer::FitRange()"; | |
1017 } | |
1018 | |
1019 outputFormat_ = format; | |
1020 Reset(); | |
1021 } | |
1022 | |
1023 | |
1024 void VolumeReslicer::SetInterpolation(ImageInterpolation interpolation) | |
1025 { | |
1026 if (interpolation != ImageInterpolation_Nearest && | |
1027 interpolation != ImageInterpolation_Bilinear && | |
1028 interpolation != ImageInterpolation_Trilinear) | |
1029 { | |
1030 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
1031 } | |
1032 | |
1033 interpolation_ = interpolation; | |
1034 Reset(); | |
1035 } | |
1036 | |
1037 | |
1038 const Extent2D& VolumeReslicer::GetOutputExtent() const | |
1039 { | |
1040 if (success_) | |
1041 { | |
1042 return extent_; | |
1043 } | |
1044 else | |
1045 { | |
1046 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
1047 } | |
1048 } | |
1049 | |
1050 | |
1051 const Orthanc::ImageAccessor& VolumeReslicer::GetOutputSlice() const | |
1052 { | |
1053 if (success_) | |
1054 { | |
1055 assert(slice_.get() != NULL); | |
1056 return *slice_; | |
1057 } | |
1058 else | |
1059 { | |
1060 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
1061 } | |
1062 } | |
1063 | |
1064 | |
1065 Orthanc::ImageAccessor* VolumeReslicer::ReleaseOutputSlice() | |
1066 { | |
1067 if (success_) | |
1068 { | |
1069 assert(slice_.get() != NULL); | |
1070 success_ = false; | |
1071 return slice_.release(); | |
1072 } | |
1073 else | |
1074 { | |
1075 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
1076 } | |
1077 } | |
1078 | |
1079 | |
1080 void VolumeReslicer::Apply(const ImageBuffer3D& source, | |
1081 const CoordinateSystem3D& plane) | |
1082 { | |
1083 // Choose the default voxel size as the finest voxel dimension | |
1084 // of the source volumetric image | |
1085 const OrthancStone::Vector dim = source.GetVoxelDimensions(OrthancStone::VolumeProjection_Axial); | |
1086 double voxelSize = dim[0]; | |
1087 | |
1088 if (dim[1] < voxelSize) | |
1089 { | |
1090 voxelSize = dim[1]; | |
1091 } | |
1092 | |
1093 if (dim[2] < voxelSize) | |
1094 { | |
1095 voxelSize = dim[2]; | |
1096 } | |
1097 | |
1098 if (voxelSize <= 0) | |
1099 { | |
1100 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
1101 } | |
1102 | |
1103 Apply(source, plane, voxelSize); | |
1104 } | |
1105 | |
1106 | |
1107 void VolumeReslicer::Apply(const ImageBuffer3D& source, | |
1108 const CoordinateSystem3D& plane, | |
1109 double voxelSize) | |
1110 { | |
1111 Reset(); | |
1112 | |
1113 // Firstly, compute the intersection of the source volumetric | |
1114 // image with the reslicing plane. This leads to a polygon with 3 | |
1115 // to 6 vertices. We compute the extent of the intersection | |
1116 // polygon, with respect to the coordinate system of the reslicing | |
1117 // plane. | |
1118 OrientedBoundingBox box(source); | |
1119 | |
1120 if (!box.ComputeExtent(extent_, plane)) | |
1121 { | |
1122 // The plane does not intersect with the bounding box of the volume | |
1123 slice_.reset(new Orthanc::Image(outputFormat_, 0, 0, false)); | |
1124 success_ = true; | |
1125 return; | |
1126 } | |
1127 | |
1128 // Secondly, the extent together with the voxel size gives the | |
1129 // size of the output image | |
1130 unsigned int width = boost::math::round(extent_.GetWidth() / voxelSize); | |
1131 unsigned int height = boost::math::round(extent_.GetHeight() / voxelSize); | |
1132 | |
1133 slice_.reset(new Orthanc::Image(outputFormat_, width, height, false)); | |
1134 | |
1135 //CheckIterators(source, plane, box); | |
1136 | |
1137 if (fastMode_) | |
1138 { | |
1139 ProcessImage<FastRowIterator>(*slice_, extent_, source, plane, box, | |
1140 interpolation_, hasLinearFunction_, scaling_, offset_); | |
1141 } | |
1142 else | |
1143 { | |
1144 ProcessImage<SlowRowIterator>(*slice_, extent_, source, plane, box, | |
1145 interpolation_, hasLinearFunction_, scaling_, offset_); | |
1146 } | |
1147 | |
1148 success_ = true; | |
1149 } | |
1150 } |