Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/ShearWarpProjectiveTransform.cpp @ 1512:244ad1e4e76a
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 } |