Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Volumes/VolumeReslicer.cpp @ 1512:244ad1e4e76a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization of folders
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 07 Jul 2020 16:21:02 +0200
parents Framework/Volumes/VolumeReslicer.cpp@30deba7bc8e2
children 4fb8fdf03314
comparison
equal deleted inserted replaced
1511:9dfeee74c1e6 1512:244ad1e4e76a
1 /**
2 * Stone of Orthanc
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
4 * Department, University Hospital of Liege, Belgium
5 * Copyright (C) 2017-2020 Osimis S.A., Belgium
6 *
7 * This program is free software: you can redistribute it and/or
8 * modify it under the terms of the GNU Affero General Public License
9 * as published by the Free Software Foundation, either version 3 of
10 * the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Affero General Public License for more details.
16 *
17 * You should have received a copy of the GNU Affero General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 **/
20
21
22 #include "VolumeReslicer.h"
23
24 #include "../Toolbox/GeometryToolbox.h"
25 #include "../Toolbox/SubvoxelReader.h"
26
27 #include <Images/ImageTraits.h>
28 #include <Logging.h>
29 #include <OrthancException.h>
30
31 #include <boost/math/special_functions/round.hpp>
32
33 namespace OrthancStone
34 {
35 // Anonymous namespace to avoid clashes between compilation modules
36 namespace
37 {
38 enum TransferFunction
39 {
40 TransferFunction_Copy,
41 TransferFunction_Float,
42 TransferFunction_Linear
43 };
44
45
46 template <Orthanc::PixelFormat InputFormat,
47 Orthanc::PixelFormat OutputFormat,
48 ImageInterpolation Interpolation,
49 TransferFunction Function>
50 class PixelShader;
51
52
53 template <Orthanc::PixelFormat Format>
54 class PixelShader<Format,
55 Format,
56 ImageInterpolation_Nearest,
57 TransferFunction_Copy>
58 {
59 private:
60 typedef SubvoxelReader<Format, ImageInterpolation_Nearest> VoxelReader;
61 typedef Orthanc::PixelTraits<Format> PixelWriter;
62
63 VoxelReader reader_;
64
65 public:
66 PixelShader(const ImageBuffer3D& image,
67 float /* scaling */,
68 float /* offset */) :
69 reader_(image)
70 {
71 }
72
73 ORTHANC_FORCE_INLINE
74 void Apply(typename PixelWriter::PixelType* pixel,
75 float volumeX,
76 float volumeY,
77 float volumeZ)
78 {
79 typename VoxelReader::PixelType value;
80
81 if (!reader_.GetValue(value, volumeX, volumeY, volumeZ))
82 {
83 VoxelReader::Traits::SetMinValue(value);
84 }
85
86 *pixel = value;
87 }
88 };
89
90
91 template <Orthanc::PixelFormat InputFormat,
92 Orthanc::PixelFormat OutputFormat>
93 class PixelShader<InputFormat,
94 OutputFormat,
95 ImageInterpolation_Nearest,
96 TransferFunction_Copy>
97 {
98 private:
99 typedef SubvoxelReader<InputFormat, ImageInterpolation_Nearest> VoxelReader;
100 typedef Orthanc::PixelTraits<OutputFormat> PixelWriter;
101
102 VoxelReader reader_;
103
104 public:
105 PixelShader(const ImageBuffer3D& image,
106 float /* scaling */,
107 float /* offset */) :
108 reader_(image)
109 {
110 }
111
112 ORTHANC_FORCE_INLINE
113 void Apply(typename PixelWriter::PixelType* pixel,
114 float volumeX,
115 float volumeY,
116 float volumeZ)
117 {
118 typename VoxelReader::PixelType value;
119
120 if (!reader_.GetValue(value, volumeX, volumeY, volumeZ))
121 {
122 VoxelReader::Traits::SetMinValue(value);
123 }
124
125 PixelWriter::FloatToPixel(*pixel, VoxelReader::Traits::PixelToFloat(value));
126 }
127 };
128
129
130 template <Orthanc::PixelFormat InputFormat,
131 Orthanc::PixelFormat OutputFormat,
132 ImageInterpolation Interpolation>
133 class PixelShader<InputFormat,
134 OutputFormat,
135 Interpolation,
136 TransferFunction_Float>
137 {
138 private:
139 typedef SubvoxelReader<InputFormat, Interpolation> VoxelReader;
140 typedef Orthanc::PixelTraits<OutputFormat> PixelWriter;
141
142 VoxelReader reader_;
143 float outOfVolume_;
144
145 public:
146 PixelShader(const ImageBuffer3D& image,
147 float /* scaling */,
148 float /* offset */) :
149 reader_(image),
150 outOfVolume_(static_cast<float>(std::numeric_limits<typename VoxelReader::PixelType>::min()))
151 {
152 }
153
154 ORTHANC_FORCE_INLINE
155 void Apply(typename PixelWriter::PixelType* pixel,
156 float volumeX,
157 float volumeY,
158 float volumeZ)
159 {
160 float value;
161
162 if (!reader_.GetFloatValue(value, volumeX, volumeY, volumeZ))
163 {
164 value = outOfVolume_;
165 }
166
167 PixelWriter::FloatToPixel(*pixel, value);
168 }
169 };
170
171
172 template <Orthanc::PixelFormat InputFormat,
173 Orthanc::PixelFormat OutputFormat,
174 ImageInterpolation Interpolation>
175 class PixelShader<InputFormat,
176 OutputFormat,
177 Interpolation,
178 TransferFunction_Linear>
179 {
180 private:
181 typedef SubvoxelReader<InputFormat, Interpolation> VoxelReader;
182 typedef Orthanc::PixelTraits<OutputFormat> PixelWriter;
183
184 VoxelReader reader_;
185 float scaling_;
186 float offset_;
187 float outOfVolume_;
188
189 public:
190 PixelShader(const ImageBuffer3D& image,
191 float scaling,
192 float offset) :
193 reader_(image),
194 scaling_(scaling),
195 offset_(offset),
196 outOfVolume_(static_cast<float>(std::numeric_limits<typename VoxelReader::PixelType>::min()))
197 {
198 }
199
200 ORTHANC_FORCE_INLINE
201 void Apply(typename PixelWriter::PixelType* pixel,
202 float volumeX,
203 float volumeY,
204 float volumeZ)
205 {
206 float value;
207
208 if (reader_.GetFloatValue(value, volumeX, volumeY, volumeZ))
209 {
210 value = scaling_ * value + offset_;
211 }
212 else
213 {
214 value = outOfVolume_;
215 }
216
217 PixelWriter::FloatToPixel(*pixel, value);
218 }
219 };
220
221
222
223 class FastRowIterator : public boost::noncopyable
224 {
225 private:
226 float position_[3];
227 float offset_[3];
228
229 public:
230 FastRowIterator(const Orthanc::ImageAccessor& slice,
231 const Extent2D& extent,
232 const CoordinateSystem3D& plane,
233 const OrientedVolumeBoundingBox& box,
234 unsigned int y)
235 {
236 const double width = static_cast<double>(slice.GetWidth());
237 const double height = static_cast<double>(slice.GetHeight());
238 assert(y < height);
239
240 Vector q1 = plane.MapSliceToWorldCoordinates
241 (extent.GetX1() + extent.GetWidth() * static_cast<double>(0) / static_cast<double>(width + 1),
242 extent.GetY1() + extent.GetHeight() * static_cast<double>(y) / static_cast<double>(height + 1));
243
244 Vector q2 = plane.MapSliceToWorldCoordinates
245 (extent.GetX1() + extent.GetWidth() * static_cast<double>(width - 1) / static_cast<double>(width + 1),
246 extent.GetY1() + extent.GetHeight() * static_cast<double>(y) / static_cast<double>(height + 1));
247
248 Vector r1, r2;
249 box.ToInternalCoordinates(r1, q1);
250 box.ToInternalCoordinates(r2, q2);
251
252 position_[0] = static_cast<float>(r1[0]);
253 position_[1] = static_cast<float>(r1[1]);
254 position_[2] = static_cast<float>(r1[2]);
255
256 Vector tmp = (r2 - r1) / static_cast<double>(width - 1);
257 offset_[0] = static_cast<float>(tmp[0]);
258 offset_[1] = static_cast<float>(tmp[1]);
259 offset_[2] = static_cast<float>(tmp[2]);
260 }
261
262 ORTHANC_FORCE_INLINE
263 void Next()
264 {
265 position_[0] += offset_[0];
266 position_[1] += offset_[1];
267 position_[2] += offset_[2];
268 }
269
270 ORTHANC_FORCE_INLINE
271 void GetVolumeCoordinates(float& x,
272 float& y,
273 float& z) const
274 {
275 x = position_[0];
276 y = position_[1];
277 z = position_[2];
278 }
279 };
280
281
282 class SlowRowIterator : public boost::noncopyable
283 {
284 private:
285 const Orthanc::ImageAccessor& slice_;
286 const Extent2D& extent_;
287 const CoordinateSystem3D& plane_;
288 const OrientedVolumeBoundingBox& box_;
289 unsigned int x_;
290 unsigned int y_;
291
292 public:
293 SlowRowIterator(const Orthanc::ImageAccessor& slice,
294 const Extent2D& extent,
295 const CoordinateSystem3D& plane,
296 const OrientedVolumeBoundingBox& box,
297 unsigned int y) :
298 slice_(slice),
299 extent_(extent),
300 plane_(plane),
301 box_(box),
302 x_(0),
303 y_(y)
304 {
305 assert(y_ < slice_.GetHeight());
306 }
307
308 void Next()
309 {
310 x_++;
311 }
312
313 void GetVolumeCoordinates(float& x,
314 float& y,
315 float& z) const
316 {
317 assert(x_ < slice_.GetWidth());
318
319 const double width = static_cast<double>(slice_.GetWidth());
320 const double height = static_cast<double>(slice_.GetHeight());
321
322 Vector q = plane_.MapSliceToWorldCoordinates
323 (extent_.GetX1() + extent_.GetWidth() * static_cast<double>(x_) / (width + 1.0),
324 extent_.GetY1() + extent_.GetHeight() * static_cast<double>(y_) / (height + 1.0));
325
326 Vector r;
327 box_.ToInternalCoordinates(r, q);
328
329 x = static_cast<float>(r[0]);
330 y = static_cast<float>(r[1]);
331 z = static_cast<float>(r[2]);
332 }
333 };
334
335
336 template <typename RowIterator,
337 Orthanc::PixelFormat InputFormat,
338 Orthanc::PixelFormat OutputFormat,
339 ImageInterpolation Interpolation,
340 TransferFunction Function>
341 static void ProcessImage(Orthanc::ImageAccessor& slice,
342 const Extent2D& extent,
343 const ImageBuffer3D& source,
344 const CoordinateSystem3D& plane,
345 const OrientedVolumeBoundingBox& box,
346 float scaling,
347 float offset)
348 {
349 typedef PixelShader<InputFormat, OutputFormat, Interpolation, Function> Shader;
350
351 const unsigned int outputWidth = slice.GetWidth();
352 const unsigned int outputHeight = slice.GetHeight();
353
354 const float sourceWidth = static_cast<float>(source.GetWidth());
355 const float sourceHeight = static_cast<float>(source.GetHeight());
356 const float sourceDepth = static_cast<float>(source.GetDepth());
357
358 Shader shader(source, scaling, offset);
359
360 for (unsigned int y = 0; y < outputHeight; y++)
361 {
362 typedef typename Orthanc::ImageTraits<OutputFormat>::PixelType PixelType;
363 PixelType* p = reinterpret_cast<PixelType*>(slice.GetRow(y));
364
365 RowIterator it(slice, extent, plane, box, y);
366
367 for (unsigned int x = 0; x < outputWidth; x++, p++)
368 {
369 float volumeX, volumeY, volumeZ;
370 it.GetVolumeCoordinates(volumeX, volumeY, volumeZ);
371
372 shader.Apply(p,
373 volumeX * sourceWidth,
374 volumeY * sourceHeight,
375 volumeZ * sourceDepth);
376 it.Next();
377 }
378 }
379 }
380
381
382 template <typename RowIterator,
383 Orthanc::PixelFormat InputFormat,
384 Orthanc::PixelFormat OutputFormat>
385 static void ProcessImage(Orthanc::ImageAccessor& slice,
386 const Extent2D& extent,
387 const ImageBuffer3D& source,
388 const CoordinateSystem3D& plane,
389 const OrientedVolumeBoundingBox& box,
390 ImageInterpolation interpolation,
391 bool hasLinearFunction,
392 float scaling,
393 float offset)
394 {
395 if (hasLinearFunction)
396 {
397 switch (interpolation)
398 {
399 case ImageInterpolation_Nearest:
400 ProcessImage<RowIterator, InputFormat, OutputFormat,
401 ImageInterpolation_Nearest, TransferFunction_Linear>
402 (slice, extent, source, plane, box, scaling, offset);
403 break;
404
405 case ImageInterpolation_Bilinear:
406 ProcessImage<RowIterator, InputFormat, OutputFormat,
407 ImageInterpolation_Bilinear, TransferFunction_Linear>
408 (slice, extent, source, plane, box, scaling, offset);
409 break;
410
411 case ImageInterpolation_Trilinear:
412 ProcessImage<RowIterator, InputFormat, OutputFormat,
413 ImageInterpolation_Trilinear, TransferFunction_Linear>
414 (slice, extent, source, plane, box, scaling, offset);
415 break;
416
417 default:
418 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
419 }
420 }
421 else
422 {
423 switch (interpolation)
424 {
425 case ImageInterpolation_Nearest:
426 ProcessImage<RowIterator, InputFormat, OutputFormat,
427 ImageInterpolation_Nearest, TransferFunction_Copy>
428 (slice, extent, source, plane, box, 0, 0);
429 break;
430
431 case ImageInterpolation_Bilinear:
432 ProcessImage<RowIterator, InputFormat, OutputFormat,
433 ImageInterpolation_Bilinear, TransferFunction_Float>
434 (slice, extent, source, plane, box, 0, 0);
435 break;
436
437 case ImageInterpolation_Trilinear:
438 ProcessImage<RowIterator, InputFormat, OutputFormat,
439 ImageInterpolation_Trilinear, TransferFunction_Float>
440 (slice, extent, source, plane, box, 0, 0);
441 break;
442
443 default:
444 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
445 }
446 }
447 }
448
449
450 template <typename RowIterator>
451 static void ProcessImage(Orthanc::ImageAccessor& slice,
452 const Extent2D& extent,
453 const ImageBuffer3D& source,
454 const CoordinateSystem3D& plane,
455 const OrientedVolumeBoundingBox& box,
456 ImageInterpolation interpolation,
457 bool hasLinearFunction,
458 float scaling,
459 float offset)
460 {
461 if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 &&
462 slice.GetFormat() == Orthanc::PixelFormat_Grayscale8)
463 {
464 ProcessImage<RowIterator,
465 Orthanc::PixelFormat_Grayscale16,
466 Orthanc::PixelFormat_Grayscale8>
467 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset);
468 }
469 else if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 &&
470 slice.GetFormat() == Orthanc::PixelFormat_Grayscale16)
471 {
472 ProcessImage<RowIterator,
473 Orthanc::PixelFormat_Grayscale16,
474 Orthanc::PixelFormat_Grayscale16>
475 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset);
476 }
477 else if (source.GetFormat() == Orthanc::PixelFormat_SignedGrayscale16 &&
478 slice.GetFormat() == Orthanc::PixelFormat_BGRA32)
479 {
480 ProcessImage<RowIterator,
481 Orthanc::PixelFormat_SignedGrayscale16,
482 Orthanc::PixelFormat_BGRA32>
483 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset);
484 }
485 else if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 &&
486 slice.GetFormat() == Orthanc::PixelFormat_BGRA32)
487 {
488 ProcessImage<RowIterator,
489 Orthanc::PixelFormat_Grayscale16,
490 Orthanc::PixelFormat_BGRA32>
491 (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset);
492 }
493 else
494 {
495 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
496 }
497 }
498 }
499
500
501
502 void VolumeReslicer::CheckIterators(const ImageBuffer3D& source,
503 const CoordinateSystem3D& plane,
504 const OrientedVolumeBoundingBox& box) const
505 {
506 for (unsigned int y = 0; y < slice_->GetHeight(); y++)
507 {
508 FastRowIterator fast(*slice_, extent_, plane, box, y);
509 SlowRowIterator slow(*slice_, extent_, plane, box, y);
510
511 for (unsigned int x = 0; x < slice_->GetWidth(); x++)
512 {
513 float px, py, pz;
514 fast.GetVolumeCoordinates(px, py, pz);
515
516 float qx, qy, qz;
517 slow.GetVolumeCoordinates(qx, qy, qz);
518
519 Vector d;
520 LinearAlgebra::AssignVector(d, px - qx, py - qy, pz - qz);
521 double norm = boost::numeric::ublas::norm_2(d);
522 if (norm > 0.0001)
523 {
524 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
525 }
526
527 fast.Next();
528 slow.Next();
529 }
530 }
531 }
532
533
534 void VolumeReslicer::Reset()
535 {
536 success_ = false;
537 extent_.Reset();
538 slice_.reset(NULL);
539 }
540
541
542 float VolumeReslicer::GetMinOutputValue() const
543 {
544 switch (outputFormat_)
545 {
546 case Orthanc::PixelFormat_Grayscale8:
547 case Orthanc::PixelFormat_Grayscale16:
548 case Orthanc::PixelFormat_BGRA32:
549 return 0.0f;
550 break;
551
552 default:
553 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
554 }
555 }
556
557
558 float VolumeReslicer::GetMaxOutputValue() const
559 {
560 switch (outputFormat_)
561 {
562 case Orthanc::PixelFormat_Grayscale8:
563 case Orthanc::PixelFormat_BGRA32:
564 return static_cast<float>(std::numeric_limits<uint8_t>::max());
565 break;
566
567 case Orthanc::PixelFormat_Grayscale16:
568 return static_cast<float>(std::numeric_limits<uint16_t>::max());
569 break;
570
571 default:
572 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
573 }
574 }
575
576
577 VolumeReslicer::VolumeReslicer() :
578 outputFormat_(Orthanc::PixelFormat_Grayscale8),
579 interpolation_(ImageInterpolation_Nearest),
580 fastMode_(true),
581 success_(false)
582 {
583 ResetLinearFunction();
584 }
585
586
587 void VolumeReslicer::GetLinearFunction(float& scaling,
588 float& offset) const
589 {
590 if (hasLinearFunction_)
591 {
592 scaling = scaling_;
593 offset = offset_;
594 }
595 else
596 {
597 scaling = 1.0f;
598 offset = 0.0f;
599 }
600 }
601
602
603 void VolumeReslicer::ResetLinearFunction()
604 {
605 Reset();
606 hasLinearFunction_ = false;
607 scaling_ = 1.0f;
608 offset_ = 0.0f;
609 }
610
611
612 void VolumeReslicer::SetLinearFunction(float scaling,
613 float offset)
614 {
615 Reset();
616 hasLinearFunction_ = true;
617 scaling_ = scaling;
618 offset_ = offset;
619 }
620
621
622 void VolumeReslicer::SetWindow(float low,
623 float high)
624 {
625 //printf("Range in pixel values: %f->%f\n", low, high);
626 float scaling = (GetMaxOutputValue() - GetMinOutputValue()) / (high - low);
627 float offset = GetMinOutputValue() - scaling * low;
628
629 SetLinearFunction(scaling, offset);
630
631 /*float x = scaling_ * low + offset_;
632 float y = scaling_ * high + offset_;
633 printf("%f %f (should be %f->%f)\n", x, y, GetMinOutputValue(), GetMaxOutputValue());*/
634 }
635
636
637 void VolumeReslicer::FitRange(const ImageBuffer3D& image)
638 {
639 float minInputValue, maxInputValue;
640
641 if (!image.GetRange(minInputValue, maxInputValue) ||
642 maxInputValue < 1)
643 {
644 ResetLinearFunction();
645 }
646 else
647 {
648 SetWindow(minInputValue, maxInputValue);
649 }
650 }
651
652
653 void VolumeReslicer::SetWindowing(ImageWindowing windowing,
654 const ImageBuffer3D& image,
655 float rescaleSlope,
656 float rescaleIntercept)
657 {
658 if (windowing == ImageWindowing_Custom)
659 {
660 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
661 }
662
663 float center, width;
664 ComputeWindowing(center, width, windowing, 0, 0);
665
666 float a = (center - width / 2.0f - rescaleIntercept) / rescaleSlope;
667 float b = (center + width / 2.0f - rescaleIntercept) / rescaleSlope;
668 SetWindow(a, b);
669 }
670
671
672 void VolumeReslicer::SetOutputFormat(Orthanc::PixelFormat format)
673 {
674 if (format != Orthanc::PixelFormat_Grayscale8 &&
675 format != Orthanc::PixelFormat_Grayscale16 &&
676 format != Orthanc::PixelFormat_BGRA32)
677 {
678 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
679 }
680
681 if (hasLinearFunction_)
682 {
683 LOG(WARNING) << "Calls to VolumeReslicer::SetOutputFormat() should be done before VolumeReslicer::FitRange()";
684 }
685
686 outputFormat_ = format;
687 Reset();
688 }
689
690
691 void VolumeReslicer::SetInterpolation(ImageInterpolation interpolation)
692 {
693 if (interpolation != ImageInterpolation_Nearest &&
694 interpolation != ImageInterpolation_Bilinear &&
695 interpolation != ImageInterpolation_Trilinear)
696 {
697 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
698 }
699
700 interpolation_ = interpolation;
701 Reset();
702 }
703
704
705 const Extent2D& VolumeReslicer::GetOutputExtent() const
706 {
707 if (success_)
708 {
709 return extent_;
710 }
711 else
712 {
713 LOG(ERROR) << "VolumeReslicer::GetOutputExtent(): (!success_)";
714
715 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
716 }
717 }
718
719
720 const Orthanc::ImageAccessor& VolumeReslicer::GetOutputSlice() const
721 {
722 if (success_)
723 {
724 assert(slice_.get() != NULL);
725 return *slice_;
726 }
727 else
728 {
729 LOG(ERROR) << "VolumeReslicer::GetOutputSlice(): (!success_)";
730 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
731 }
732 }
733
734
735 Orthanc::ImageAccessor* VolumeReslicer::ReleaseOutputSlice()
736 {
737 if (success_)
738 {
739 assert(slice_.get() != NULL);
740 success_ = false;
741 return slice_.release();
742 }
743 else
744 {
745 LOG(ERROR) << "VolumeReslicer::ReleaseOutputSlice(): (!success_)";
746 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
747 }
748 }
749
750
751 void VolumeReslicer::Apply(const ImageBuffer3D& source,
752 const VolumeImageGeometry& geometry,
753 const CoordinateSystem3D& plane)
754 {
755 // Choose the default voxel size as the finest voxel dimension
756 // of the source volumetric image
757 const OrthancStone::Vector dim =
758 geometry.GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
759 double voxelSize = dim[0];
760
761 if (dim[1] < voxelSize)
762 {
763 voxelSize = dim[1];
764 }
765
766 if (dim[2] < voxelSize)
767 {
768 voxelSize = dim[2];
769 }
770
771 if (voxelSize <= 0)
772 {
773 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
774 }
775
776 Apply(source, geometry, plane, voxelSize);
777 }
778
779
780 void VolumeReslicer::Apply(const ImageBuffer3D& source,
781 const VolumeImageGeometry& geometry,
782 const CoordinateSystem3D& plane,
783 double voxelSize)
784 {
785 Reset();
786 pixelSpacing_ = voxelSize;
787
788 // Firstly, compute the intersection of the source volumetric
789 // image with the reslicing plane. This leads to a polygon with 3
790 // to 6 vertices. We compute the extent of the intersection
791 // polygon, with respect to the coordinate system of the reslicing
792 // plane.
793 OrientedVolumeBoundingBox box(geometry);
794
795 if (!box.ComputeExtent(extent_, plane))
796 {
797 // The plane does not intersect with the bounding box of the volume
798 slice_.reset(new Orthanc::Image(outputFormat_, 0, 0, false));
799 success_ = true;
800 return;
801 }
802
803 // Secondly, the extent together with the voxel size gives the
804 // size of the output image
805 unsigned int width = boost::math::iround(extent_.GetWidth() / voxelSize);
806 unsigned int height = boost::math::iround(extent_.GetHeight() / voxelSize);
807
808 slice_.reset(new Orthanc::Image(outputFormat_, width, height, false));
809
810 //CheckIterators(source, plane, box);
811
812 if (fastMode_)
813 {
814 ProcessImage<FastRowIterator>(*slice_, extent_, source, plane, box,
815 interpolation_, hasLinearFunction_, scaling_, offset_);
816 }
817 else
818 {
819 ProcessImage<SlowRowIterator>(*slice_, extent_, source, plane, box,
820 interpolation_, hasLinearFunction_, scaling_, offset_);
821 }
822
823 success_ = true;
824 }
825
826
827 double VolumeReslicer::GetPixelSpacing() const
828 {
829 if (success_)
830 {
831 return pixelSpacing_;
832 }
833 else
834 {
835 LOG(ERROR) << "VolumeReslicer::GetPixelSpacing(): (!success_)";
836 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
837 }
838 }
839 }