Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/VolumeReslicer.cpp @ 153:ae531ab5dcd9 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
new class: VolumeReslicer
author Sebastien Jodogne <s.jodogne@gmail.com>
date Fri, 02 Feb 2018 09:35:07 +0100
parents
children
comparison
equal deleted inserted replaced
152:2e023be0563c 153:ae531ab5dcd9
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 }