Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/ShearWarpProjectiveTransform.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/Toolbox/ShearWarpProjectiveTransform.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 "ShearWarpProjectiveTransform.h"
23
24 #include "ImageGeometry.h"
25 #include "Extent2D.h"
26 #include "FiniteProjectiveCamera.h"
27 #include "GeometryToolbox.h"
28
29 #include <Images/PixelTraits.h>
30 #include <Images/ImageProcessing.h>
31 #include <OrthancException.h>
32 #include <Logging.h>
33
34 #include <boost/numeric/ublas/matrix_proxy.hpp>
35 #include <boost/math/special_functions/round.hpp>
36 #include <cassert>
37
38
39 namespace OrthancStone
40 {
41 static bool IsValidShear(const Matrix& M_shear)
42 {
43 return (LinearAlgebra::IsCloseToZero(M_shear(0, 1)) &&
44 LinearAlgebra::IsCloseToZero(M_shear(1, 0)) &&
45 LinearAlgebra::IsCloseToZero(M_shear(2, 0)) &&
46 LinearAlgebra::IsCloseToZero(M_shear(2, 1)) &&
47 LinearAlgebra::IsNear(1.0, M_shear(2, 2)) &&
48 LinearAlgebra::IsCloseToZero(M_shear(2, 3)) &&
49 LinearAlgebra::IsCloseToZero(M_shear(3, 0)) &&
50 LinearAlgebra::IsCloseToZero(M_shear(3, 1)) &&
51 LinearAlgebra::IsNear(1.0, M_shear(3, 3)));
52 }
53
54
55 static void ComputeShearParameters(double& scaling,
56 double& offsetX,
57 double& offsetY,
58 const Matrix& shear,
59 double z)
60 {
61 // Check out: ../../Resources/Computations/ComputeShearParameters.py
62
63 if (!LinearAlgebra::IsShearMatrix(shear))
64 {
65 LOG(ERROR) << "Not a valid shear matrix";
66 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
67 }
68
69 scaling = 1.0 / (shear(3,2) * z + 1.0);
70 offsetX = shear(0,2) * z * scaling;
71 offsetY = shear(1,2) * z * scaling;
72 }
73
74
75 ShearWarpProjectiveTransform::
76 ShearWarpProjectiveTransform(const Matrix& M_view,
77 //const Matrix& P, // Permutation applied to the volume
78 unsigned int volumeWidth,
79 unsigned int volumeHeight,
80 unsigned int volumeDepth,
81 double pixelSpacingX,
82 double pixelSpacingY,
83 unsigned int imageWidth,
84 unsigned int imageHeight)
85 {
86 eye_o.resize(4);
87
88 {
89 // Find back the camera center given the "M_view" matrix
90 const double m11 = M_view(0, 0);
91 const double m12 = M_view(0, 1);
92 const double m13 = M_view(0, 2);
93 const double m14 = M_view(0, 3);
94 const double m21 = M_view(1, 0);
95 const double m22 = M_view(1, 1);
96 const double m23 = M_view(1, 2);
97 const double m24 = M_view(1, 3);
98 const double m41 = M_view(3, 0);
99 const double m42 = M_view(3, 1);
100 const double m43 = M_view(3, 2);
101 const double m44 = M_view(3, 3);
102
103 // Equations (A.8) to (A.11) on page 203. Also check out
104 // "Finding the camera center" in "Multiple View Geometry in
105 // Computer Vision - 2nd edition", page 163.
106 const double vx[9] = { m12, m13, m14, m22, m23, m24, m42, m43, m44 };
107 const double vy[9] = { m11, m13, m14, m21, m23, m24, m41, m43, m44 };
108 const double vz[9] = { m11, m12, m14, m21, m22, m24, m41, m42, m44 };
109 const double vw[9] = { m11, m12, m13, m21, m22, m23, m41, m42, m43 };
110
111 Matrix m;
112
113 LinearAlgebra::FillMatrix(m, 3, 3, vx);
114 eye_o[0] = -LinearAlgebra::ComputeDeterminant(m);
115
116 LinearAlgebra::FillMatrix(m, 3, 3, vy);
117 eye_o[1] = LinearAlgebra::ComputeDeterminant(m);
118
119 LinearAlgebra::FillMatrix(m, 3, 3, vz);
120 eye_o[2] = -LinearAlgebra::ComputeDeterminant(m);
121
122 LinearAlgebra::FillMatrix(m, 3, 3, vw);
123 eye_o[3] = LinearAlgebra::ComputeDeterminant(m);
124
125 if (LinearAlgebra::IsCloseToZero(eye_o[3]))
126 {
127 LOG(ERROR) << "The shear-warp projective transform is not applicable to affine cameras";
128 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
129 }
130 }
131
132 #if 0
133 // Assume "T_shift = I" (the eye does not lie on plane k = 0)
134 const Matrix T_shift = LinearAlgebra::IdentityMatrix(4);
135
136 // Equation (A.13) on page 204, given that the inverse of a
137 // permutation matrix is its transpose (TODO CHECK). If no T_shift
138 // or permutation P is applied, M'_view == M_view
139 const Matrix MM_view = LinearAlgebra::Product(
140 M_view,
141 LinearAlgebra::Transpose(P),
142 LinearAlgebra::InvertScalingTranslationMatrix(T_shift));
143 #else
144 // This is a shortcut, as we take "T_shift = I" and "P = I"
145 const Matrix MM_view = M_view;
146 #endif
147
148 // Equation (A.14) on page 207
149 Matrix MM_shear = LinearAlgebra::IdentityMatrix(4);
150 MM_shear(0, 2) = -eye_o[0] / eye_o[2];
151 MM_shear(1, 2) = -eye_o[1] / eye_o[2];
152 MM_shear(3, 2) = -eye_o[3] / eye_o[2];
153
154
155 // Compute the extent of the intermediate image
156 Extent2D extent;
157 double maxScaling = 1;
158
159 {
160 // Compute the shearing factors of the two extreme planes of the
161 // volume (z=0 and z=volumeDepth)
162 double scaling, offsetX, offsetY;
163 ComputeShearParameters(scaling, offsetX, offsetY, MM_shear, 0);
164
165 if (scaling > 0)
166 {
167 extent.AddPoint(offsetX, offsetY);
168 extent.AddPoint(offsetX + static_cast<double>(volumeWidth) * scaling,
169 offsetY + static_cast<double>(volumeHeight) * scaling);
170
171 if (scaling > maxScaling)
172 {
173 maxScaling = scaling;
174 }
175 }
176
177 ComputeShearParameters(scaling, offsetX, offsetY, MM_shear, volumeDepth);
178
179 if (scaling > 0)
180 {
181 extent.AddPoint(offsetX, offsetY);
182 extent.AddPoint(offsetX + static_cast<double>(volumeWidth) * scaling,
183 offsetY + static_cast<double>(volumeHeight) * scaling);
184
185 if (scaling > maxScaling)
186 {
187 maxScaling = scaling;
188 }
189 }
190 }
191
192 if (LinearAlgebra::IsCloseToZero(extent.GetWidth()) ||
193 LinearAlgebra::IsCloseToZero(extent.GetHeight()))
194 {
195 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
196 }
197
198 intermediateWidth_ =
199 static_cast<unsigned int>(std::ceil(extent.GetWidth() / maxScaling));
200 intermediateHeight_ =
201 static_cast<unsigned int>(std::ceil(extent.GetHeight() / maxScaling));
202
203 // This is the product "T * S" in Equation (A.16) on page 209
204 Matrix TS = LinearAlgebra::Product(
205 GeometryToolbox::CreateTranslationMatrix(
206 static_cast<double>(intermediateWidth_) / 2.0,
207 static_cast<double>(intermediateHeight_) / 2.0, 0),
208 GeometryToolbox::CreateScalingMatrix(
209 1.0 / maxScaling, 1.0 / maxScaling, 1),
210 GeometryToolbox::CreateTranslationMatrix(
211 -extent.GetCenterX(), -extent.GetCenterY(), 0));
212
213 // This is Equation (A.16) on page 209. WARNING: There is an
214 // error in Lacroute's thesis: "inv(MM_shear)" is used instead
215 // of "MM_shear".
216 M_shear = LinearAlgebra::Product(TS, MM_shear);
217
218 if (!IsValidShear(M_shear))
219 {
220 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
221 }
222
223 // This is Equation (A.17) on page 209
224 Matrix tmp;
225 LinearAlgebra::InvertMatrix(tmp, M_shear);
226 M_warp = LinearAlgebra::Product(MM_view, tmp);
227
228 // Intrinsic parameters of the camera
229 k_ = LinearAlgebra::ZeroMatrix(3, 4);
230 k_(0, 0) = 1.0 / pixelSpacingX;
231 k_(0, 3) = static_cast<double>(imageWidth) / 2.0;
232 k_(1, 1) = 1.0 / pixelSpacingY;
233 k_(1, 3) = static_cast<double>(imageHeight) / 2.0;
234 k_(2, 3) = 1.0;
235 }
236
237
238 FiniteProjectiveCamera *ShearWarpProjectiveTransform::CreateCamera() const
239 {
240 Matrix p = LinearAlgebra::Product(k_, M_warp, M_shear);
241 return new FiniteProjectiveCamera(p);
242 }
243
244
245 void ShearWarpProjectiveTransform::ComputeShearOnSlice(double& a11,
246 double& b1,
247 double& a22,
248 double& b2,
249 double& shearedZ,
250 const double sourceZ)
251 {
252 // Check out: ../../Resources/Computations/ComputeShearOnSlice.py
253 assert(IsValidShear(M_shear));
254
255 const double s11 = M_shear(0, 0);
256 const double s13 = M_shear(0, 2);
257 const double s14 = M_shear(0, 3);
258 const double s22 = M_shear(1, 1);
259 const double s23 = M_shear(1, 2);
260 const double s24 = M_shear(1, 3);
261 const double s43 = M_shear(3, 2);
262
263 double scaling = 1.0 / (s43 * sourceZ + 1.0);
264 shearedZ = sourceZ * scaling;
265
266 a11 = s11 * scaling;
267 a22 = s22 * scaling;
268
269 b1 = (s13 * sourceZ + s14) * scaling;
270 b2 = (s23 * sourceZ + s24) * scaling;
271 }
272
273
274 Matrix ShearWarpProjectiveTransform::CalibrateView(const Vector& camera,
275 const Vector& principalPoint,
276 double angle)
277 {
278 if (camera.size() != 3 ||
279 principalPoint.size() != 3)
280 {
281 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
282 }
283
284 const double sid = boost::numeric::ublas::norm_2(camera - principalPoint);
285
286 Matrix a;
287 GeometryToolbox::AlignVectorsWithRotation(a, camera - principalPoint,
288 LinearAlgebra::CreateVector(0, 0, -1));
289
290 Matrix r = LinearAlgebra::Product(GeometryToolbox::CreateRotationMatrixAlongZ(angle), a);
291
292 a = LinearAlgebra::ZeroMatrix(4, 4);
293 boost::numeric::ublas::subrange(a, 0, 3, 0, 3) = r;
294
295 const Vector v = LinearAlgebra::Product(r, -camera);
296 a(0, 3) = v[0];
297 a(1, 3) = v[1];
298 a(2, 3) = v[2];
299 a(3, 3) = 1;
300
301 Matrix perspective = LinearAlgebra::ZeroMatrix(4, 4);
302 // https://stackoverflow.com/questions/5267866/calculation-of-a-perspective-transformation-matrix
303 perspective(0, 0) = sid;
304 perspective(1, 1) = sid;
305 perspective(2, 2) = sid;
306 perspective(3, 2) = 1;
307
308 Matrix M_view = LinearAlgebra::Product(perspective, a);
309 assert(M_view.size1() == 4 &&
310 M_view.size2() == 4);
311
312 {
313 // Sanity checks
314 Vector p1 = LinearAlgebra::CreateVector(camera[0], camera[1], camera[2], 1.0);
315 Vector p2 = LinearAlgebra::CreateVector(principalPoint[0], principalPoint[1], principalPoint[2], 1.0);
316
317 Vector v1 = LinearAlgebra::Product(M_view, p1);
318 Vector v2 = LinearAlgebra::Product(M_view, p2);
319
320 if (!LinearAlgebra::IsCloseToZero(v1[3]) || // Must be mapped to singularity (w=0)
321 LinearAlgebra::IsCloseToZero(v2[3]))
322 {
323 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
324 }
325
326 // The principal point must be mapped to (0,0,z,1)
327 v2 /= v2[3];
328 if (!LinearAlgebra::IsCloseToZero(v2[0]) ||
329 !LinearAlgebra::IsCloseToZero(v2[1]))
330 {
331 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
332 }
333 }
334
335 return M_view;
336 }
337
338
339 template <Orthanc::PixelFormat SourceFormat,
340 Orthanc::PixelFormat TargetFormat,
341 bool MIP>
342 static void ApplyAxialInternal(Orthanc::ImageAccessor& target,
343 float& maxValue,
344 const Matrix& M_view,
345 const ImageBuffer3D& source,
346 const VolumeImageGeometry& geometry,
347 double pixelSpacing,
348 unsigned int countSlices,
349 ImageInterpolation shearInterpolation,
350 ImageInterpolation warpInterpolation)
351 {
352 typedef Orthanc::PixelTraits<SourceFormat> SourceTraits;
353 typedef Orthanc::PixelTraits<TargetFormat> TargetTraits;
354
355 /**
356 * Step 1: Precompute some information.
357 **/
358
359 if (target.GetFormat() != TargetFormat ||
360 source.GetFormat() != SourceFormat ||
361 !std::numeric_limits<float>::is_iec559 ||
362 sizeof(float) != 4)
363 {
364 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
365 }
366
367 if (countSlices > source.GetDepth())
368 {
369 countSlices = source.GetDepth();
370 }
371
372 if (countSlices == 0)
373 {
374 maxValue = 0;
375 Orthanc::ImageProcessing::Set(target, 0);
376 return;
377 }
378
379 LOG(INFO) << "Number of rendered slices: " << countSlices;
380
381
382 /**
383 * Step 2: Extract the shear-warp transform corresponding to
384 * M_view.
385 **/
386
387 // Compute the "world" matrix that maps the source volume to the
388 // (0,0,0)->(1,1,1) unit cube
389 Vector origin = geometry.GetCoordinates(0, 0, 0);
390 Vector ps = geometry.GetVoxelDimensions(VolumeProjection_Axial);
391 Matrix world = LinearAlgebra::Product(
392 GeometryToolbox::CreateScalingMatrix(1.0 / ps[0], 1.0 / ps[1], 1.0 / ps[2]),
393 GeometryToolbox::CreateTranslationMatrix(-origin[0], -origin[1], -origin[2]));
394
395 Matrix worldInv;
396 LinearAlgebra::InvertMatrix(worldInv, world);
397
398 ShearWarpProjectiveTransform shearWarp(LinearAlgebra::Product(M_view, worldInv),
399 /*LinearAlgebra::IdentityMatrix(4),*/
400 source.GetWidth(),
401 source.GetHeight(),
402 source.GetDepth(),
403 pixelSpacing, pixelSpacing,
404 target.GetWidth(), target.GetHeight());
405
406 const unsigned int intermediateWidth = shearWarp.GetIntermediateWidth();
407 const unsigned int intermediateHeight = shearWarp.GetIntermediateHeight();
408
409
410 /**
411 * Step 3: Apply the "shear" part of the transform to form the
412 * intermediate image. The sheared images are accumulated into the
413 * Float32 image "accumulator". The number of samples available
414 * for each pixel is stored in the "counter" image.
415 **/
416
417 std::unique_ptr<Orthanc::ImageAccessor> accumulator, counter, intermediate;
418
419 accumulator.reset(new Orthanc::Image(Orthanc::PixelFormat_Float32,
420 intermediateWidth, intermediateHeight, false));
421 counter.reset(new Orthanc::Image(Orthanc::PixelFormat_Grayscale16,
422 intermediateWidth, intermediateHeight, false));
423 intermediate.reset(new Orthanc::Image(SourceFormat, intermediateWidth, intermediateHeight, false));
424
425 Orthanc::ImageProcessing::Set(*accumulator, 0);
426 Orthanc::ImageProcessing::Set(*counter, 0);
427
428 // Loop around the slices of the volume
429 for (unsigned int i = 0; i <= countSlices; i++)
430 {
431 // (3.a) Compute the shear for this specific slice
432 unsigned int z = static_cast<unsigned int>(
433 boost::math::iround(static_cast<double>(i) /
434 static_cast<double>(countSlices) *
435 static_cast<double>(source.GetDepth() - 1)));
436
437 double a11, b1, a22, b2, vz;
438 shearWarp.ComputeShearOnSlice(a11, b1, a22, b2, vz, static_cast<double>(z) + 0.5);
439
440
441 {
442 // (3.b) Detect the "useful" portion of the intermediate image
443 // for this slice (i.e. the bounding box where the source
444 // slice is mapped to by the shear), so as to update "counter"
445 Matrix a = LinearAlgebra::ZeroMatrix(3, 3);
446 a(0,0) = a11;
447 a(0,2) = b1;
448 a(1,1) = a22;
449 a(1,2) = b2;
450 a(2,2) = 1;
451
452 unsigned int x1, y1, x2, y2;
453 if (GetProjectiveTransformExtent(x1, y1, x2, y2, a,
454 source.GetWidth(), source.GetHeight(),
455 intermediateWidth, intermediateHeight))
456 {
457 for (unsigned int y = y1; y <= y2; y++)
458 {
459 uint16_t* p = reinterpret_cast<uint16_t*>(counter->GetRow(y)) + x1;
460 for (unsigned int x = x1; x <= x2; x++, p++)
461 {
462 if (MIP)
463 {
464 // TODO - In the case of MIP, "counter" could be
465 // reduced to "PixelFormat_Grayscale8" to reduce
466 // memory usage
467 *p = 1;
468 }
469 else
470 {
471 *p += 1;
472 }
473 }
474 }
475 }
476 }
477
478
479 {
480 // (3.c) Shear the source slice into a temporary image
481 ImageBuffer3D::SliceReader reader(source, VolumeProjection_Axial, z);
482 ApplyAffineTransform(*intermediate, reader.GetAccessor(),
483 a11, 0, b1,
484 0, a22, b2,
485 shearInterpolation, true);
486 }
487
488
489 for (unsigned int y = 0; y < intermediateHeight; y++)
490 {
491 // (3.d) Accumulate the pixels of the sheared image into "accumulator"
492 const typename SourceTraits::PixelType* p =
493 reinterpret_cast<const typename SourceTraits::PixelType*>(intermediate->GetConstRow(y));
494
495 float* q = reinterpret_cast<float*>(accumulator->GetRow(y));
496
497 for (unsigned int x = 0; x < intermediateWidth; x++)
498 {
499 float pixel = SourceTraits::PixelToFloat(*p);
500
501 if (MIP)
502 {
503 // Get maximum for MIP
504 if (*q < pixel)
505 {
506 *q = pixel;
507 }
508 }
509 else
510 {
511 *q += pixel;
512 }
513
514 p++;
515 q++;
516 }
517 }
518 }
519
520
521 /**
522 * Step 4: The intermediate image (that will be transformed by the
523 * "warp") is now available as an accumulator image together with
524 * a counter image. "Flatten" these two images into one.
525 **/
526
527 intermediate.reset(new Orthanc::Image
528 (TargetFormat, intermediateWidth, intermediateHeight, false));
529
530 maxValue = 0;
531
532 for (unsigned int y = 0; y < intermediateHeight; y++)
533 {
534 const float *qacc = reinterpret_cast<const float*>(accumulator->GetConstRow(y));
535 const uint16_t *qcount = reinterpret_cast<const uint16_t*>(counter->GetConstRow(y));
536 typename TargetTraits::PixelType *p =
537 reinterpret_cast<typename TargetTraits::PixelType*>(intermediate->GetRow(y));
538
539 for (unsigned int x = 0; x < intermediateWidth; x++)
540 {
541 if (*qcount == 0)
542 {
543 TargetTraits::SetZero(*p);
544 }
545 else
546 {
547 *p = static_cast<typename TargetTraits::PixelType>
548 (*qacc / static_cast<float>(*qcount));
549
550 if (*p > maxValue)
551 {
552 maxValue = *p;
553 }
554 }
555
556 p++;
557 qacc++;
558 qcount++;
559 }
560 }
561
562 // We don't need the accumulator images anymore
563 accumulator.reset(NULL);
564 counter.reset(NULL);
565
566
567 /**
568 * Step 6: Apply the "warp" part of the transform to map the
569 * intermediate image to the final image.
570 **/
571
572 Matrix warp;
573
574 {
575 // (5.a) Compute the "warp" matrix by removing the 3rd row and
576 // 3rd column from the GetWarp() matrix
577 // Check out: ../../Resources/Computations/ComputeWarp.py
578
579 Matrix fullWarp = LinearAlgebra::Product
580 (shearWarp.GetIntrinsicParameters(), shearWarp.GetWarp());
581
582 const double v[] = {
583 fullWarp(0,0), fullWarp(0,1), fullWarp(0,3),
584 fullWarp(1,0), fullWarp(1,1), fullWarp(1,3),
585 fullWarp(2,0), fullWarp(2,1), fullWarp(2,3)
586 };
587
588 LinearAlgebra::FillMatrix(warp, 3, 3, v);
589 }
590
591 // (5.b) Apply the projective transform to the image
592 ApplyProjectiveTransform(target, *intermediate, warp, warpInterpolation, true);
593 }
594
595
596 template <Orthanc::PixelFormat SourceFormat,
597 Orthanc::PixelFormat TargetFormat>
598 static void ApplyAxialInternal2(Orthanc::ImageAccessor& target,
599 float& maxValue,
600 const Matrix& M_view,
601 const ImageBuffer3D& source,
602 const VolumeImageGeometry& geometry,
603 bool mip,
604 double pixelSpacing,
605 unsigned int countSlices,
606 ImageInterpolation shearInterpolation,
607 ImageInterpolation warpInterpolation)
608 {
609 if (mip)
610 {
611 ApplyAxialInternal<SourceFormat, TargetFormat, true>
612 (target, maxValue, M_view, source, geometry, pixelSpacing,
613 countSlices, shearInterpolation, warpInterpolation);
614 }
615 else
616 {
617 ApplyAxialInternal<SourceFormat, TargetFormat, false>
618 (target, maxValue, M_view, source, geometry, pixelSpacing,
619 countSlices, shearInterpolation, warpInterpolation);
620 }
621 }
622
623
624 Orthanc::ImageAccessor*
625 ShearWarpProjectiveTransform::ApplyAxial(float& maxValue,
626 const Matrix& M_view,
627 const ImageBuffer3D& source,
628 const VolumeImageGeometry& geometry,
629 Orthanc::PixelFormat targetFormat,
630 unsigned int targetWidth,
631 unsigned int targetHeight,
632 bool mip,
633 double pixelSpacing,
634 unsigned int countSlices,
635 ImageInterpolation shearInterpolation,
636 ImageInterpolation warpInterpolation)
637 {
638 std::unique_ptr<Orthanc::ImageAccessor> target
639 (new Orthanc::Image(targetFormat, targetWidth, targetHeight, false));
640
641 if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 &&
642 targetFormat == Orthanc::PixelFormat_Grayscale16)
643 {
644 ApplyAxialInternal2<Orthanc::PixelFormat_Grayscale16,
645 Orthanc::PixelFormat_Grayscale16>
646 (*target, maxValue, M_view, source, geometry, mip, pixelSpacing,
647 countSlices, shearInterpolation, warpInterpolation);
648 }
649 else if (source.GetFormat() == Orthanc::PixelFormat_SignedGrayscale16 &&
650 targetFormat == Orthanc::PixelFormat_SignedGrayscale16)
651 {
652 ApplyAxialInternal2<Orthanc::PixelFormat_SignedGrayscale16,
653 Orthanc::PixelFormat_SignedGrayscale16>
654 (*target, maxValue, M_view, source, geometry, mip, pixelSpacing,
655 countSlices, shearInterpolation, warpInterpolation);
656 }
657 else
658 {
659 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
660 }
661
662 return target.release();
663 }
664 }