Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 191:46cb2eedc2e0 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ShearWarpProjectiveTransform
author Sebastien Jodogne <s.jodogne@gmail.com>
date Fri, 16 Mar 2018 15:01:52 +0100
parents
children 4abddd083374
comparison
equal deleted inserted replaced
190:465b294a55f0 191:46cb2eedc2e0
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-2018 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 <Core/OrthancException.h>
30 #include <Core/Logging.h>
31
32 #include <boost/numeric/ublas/matrix_proxy.hpp>
33 #include <cassert>
34
35
36 namespace OrthancStone
37 {
38 static bool IsValidShear(const Matrix& M_shear)
39 {
40 return (LinearAlgebra::IsCloseToZero(M_shear(0, 1)) &&
41 LinearAlgebra::IsCloseToZero(M_shear(1, 0)) &&
42 LinearAlgebra::IsCloseToZero(M_shear(2, 0)) &&
43 LinearAlgebra::IsCloseToZero(M_shear(2, 1)) &&
44 LinearAlgebra::IsNear(1.0, M_shear(2, 2)) &&
45 LinearAlgebra::IsCloseToZero(M_shear(2, 3)) &&
46 LinearAlgebra::IsCloseToZero(M_shear(3, 0)) &&
47 LinearAlgebra::IsCloseToZero(M_shear(3, 1)) &&
48 LinearAlgebra::IsNear(1.0, M_shear(3, 3)));
49 }
50
51
52 static void ComputeShearParameters(double& scaling,
53 double& offsetX,
54 double& offsetY,
55 const Matrix& shear,
56 double z)
57 {
58 // Check out: ../../Resources/Computations/ComputeShearParameters.py
59
60 if (!LinearAlgebra::IsShearMatrix(shear))
61 {
62 LOG(ERROR) << "Not a valid shear matrix";
63 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
64 }
65
66 scaling = 1.0 / (shear(3,2) * z + 1.0);
67 offsetX = shear(0,2) * z * scaling;
68 offsetY = shear(1,2) * z * scaling;
69 }
70
71
72 ShearWarpProjectiveTransform::
73 ShearWarpProjectiveTransform(const Matrix& M_view,
74 //const Matrix& P, // Permutation applied to the volume
75 unsigned int volumeWidth,
76 unsigned int volumeHeight,
77 unsigned int volumeDepth,
78 double pixelSpacingX,
79 double pixelSpacingY,
80 unsigned int imageWidth,
81 unsigned int imageHeight)
82 {
83 eye_o.resize(4);
84
85 {
86 // Find back the camera center given the "M_view" matrix
87 const double m11 = M_view(0, 0);
88 const double m12 = M_view(0, 1);
89 const double m13 = M_view(0, 2);
90 const double m14 = M_view(0, 3);
91 const double m21 = M_view(1, 0);
92 const double m22 = M_view(1, 1);
93 const double m23 = M_view(1, 2);
94 const double m24 = M_view(1, 3);
95 const double m41 = M_view(3, 0);
96 const double m42 = M_view(3, 1);
97 const double m43 = M_view(3, 2);
98 const double m44 = M_view(3, 3);
99
100 // Equations (A.8) to (A.11) on page 203. Also check out
101 // "Finding the camera center" in "Multiple View Geometry in
102 // Computer Vision - 2nd edition", page 163.
103 const double vx[9] = { m12, m13, m14, m22, m23, m24, m42, m43, m44 };
104 const double vy[9] = { m11, m13, m14, m21, m23, m24, m41, m43, m44 };
105 const double vz[9] = { m11, m12, m14, m21, m22, m24, m41, m42, m44 };
106 const double vw[9] = { m11, m12, m13, m21, m22, m23, m41, m42, m43 };
107
108 Matrix m;
109
110 LinearAlgebra::FillMatrix(m, 3, 3, vx);
111 eye_o[0] = -LinearAlgebra::ComputeDeterminant(m);
112
113 LinearAlgebra::FillMatrix(m, 3, 3, vy);
114 eye_o[1] = LinearAlgebra::ComputeDeterminant(m);
115
116 LinearAlgebra::FillMatrix(m, 3, 3, vz);
117 eye_o[2] = -LinearAlgebra::ComputeDeterminant(m);
118
119 LinearAlgebra::FillMatrix(m, 3, 3, vw);
120 eye_o[3] = LinearAlgebra::ComputeDeterminant(m);
121
122 if (LinearAlgebra::IsCloseToZero(eye_o[3]))
123 {
124 LOG(ERROR) << "The shear-warp projective transform is not applicable to affine cameras";
125 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
126 }
127 }
128
129 #if 0
130 // Assume "T_shift = I" (the eye does not lie on plane k = 0)
131 const Matrix T_shift = LinearAlgebra::IdentityMatrix(4);
132
133 // Equation (A.13) on page 204, given that the inverse of a
134 // permutation matrix is its transpose (TODO CHECK). If no T_shift
135 // or permutation P is applied, M'_view == M_view
136 const Matrix MM_view = LinearAlgebra::Product(
137 M_view,
138 LinearAlgebra::Transpose(P),
139 LinearAlgebra::InvertScalingTranslationMatrix(T_shift));
140 #else
141 // This is a shortcut, as we take "T_shift = I" and "P = I"
142 const Matrix MM_view = M_view;
143 #endif
144
145 // Equation (A.14) on page 207
146 Matrix MM_shear = LinearAlgebra::IdentityMatrix(4);
147 MM_shear(0, 2) = -eye_o[0] / eye_o[2];
148 MM_shear(1, 2) = -eye_o[1] / eye_o[2];
149 MM_shear(3, 2) = -eye_o[3] / eye_o[2];
150
151
152 // Compute the extent of the intermediate image
153 Extent2D extent;
154 double maxScaling = 1;
155
156 {
157 // Compute the shearing factors of the two extreme planes of the
158 // volume (z=0 and z=volumeDepth)
159 double scaling, offsetX, offsetY;
160 ComputeShearParameters(scaling, offsetX, offsetY, MM_shear, 0);
161
162 if (scaling > 0)
163 {
164 extent.AddPoint(offsetX, offsetY);
165 extent.AddPoint(offsetX + static_cast<double>(volumeWidth) * scaling,
166 offsetY + static_cast<double>(volumeHeight) * scaling);
167
168 if (scaling > maxScaling)
169 {
170 maxScaling = scaling;
171 }
172 }
173
174 ComputeShearParameters(scaling, offsetX, offsetY, MM_shear, volumeDepth);
175
176 if (scaling > 0)
177 {
178 extent.AddPoint(offsetX, offsetY);
179 extent.AddPoint(offsetX + static_cast<double>(volumeWidth) * scaling,
180 offsetY + static_cast<double>(volumeHeight) * scaling);
181
182 if (scaling > maxScaling)
183 {
184 maxScaling = scaling;
185 }
186 }
187 }
188
189 if (LinearAlgebra::IsCloseToZero(extent.GetWidth()) ||
190 LinearAlgebra::IsCloseToZero(extent.GetHeight()))
191 {
192 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
193 }
194
195 intermediateWidth_ = std::ceil(extent.GetWidth() / maxScaling);
196 intermediateHeight_ = std::ceil(extent.GetHeight() / maxScaling);
197
198 // This is the product "T * S" in Equation (A.16) on page 209
199 Matrix TS = LinearAlgebra::Product(
200 GeometryToolbox::CreateTranslationMatrix(static_cast<double>(intermediateWidth_) / 2.0,
201 static_cast<double>(intermediateHeight_) / 2.0, 0),
202 GeometryToolbox::CreateScalingMatrix(1.0 / maxScaling, 1.0 / maxScaling, 1),
203 GeometryToolbox::CreateTranslationMatrix(-extent.GetCenterX(), -extent.GetCenterY(), 0));
204
205 // This is Equation (A.16) on page 209. WARNING: There is an
206 // error in Lacroute's thesis: "inv(MM_shear)" is used instead
207 // of "MM_shear".
208 M_shear = LinearAlgebra::Product(TS, MM_shear);
209
210 if (!IsValidShear(M_shear))
211 {
212 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
213 }
214
215 // This is Equation (A.17) on page 209
216 Matrix tmp;
217 LinearAlgebra::InvertMatrix(tmp, M_shear);
218 M_warp = LinearAlgebra::Product(MM_view, tmp);
219
220 // Intrinsic parameters of the camera
221 k_ = LinearAlgebra::ZeroMatrix(3, 4);
222 k_(0, 0) = 1.0 / pixelSpacingX;
223 k_(0, 3) = static_cast<double>(imageWidth) / 2.0;
224 k_(1, 1) = 1.0 / pixelSpacingY;
225 k_(1, 3) = static_cast<double>(imageHeight) / 2.0;
226 k_(2, 3) = 1.0;
227 }
228
229
230 FiniteProjectiveCamera *ShearWarpProjectiveTransform::CreateCamera() const
231 {
232 Matrix p = LinearAlgebra::Product(k_, M_warp, M_shear);
233 return new FiniteProjectiveCamera(p);
234 }
235
236
237 void ShearWarpProjectiveTransform::ComputeShearOnSlice(double& a11,
238 double& b1,
239 double& a22,
240 double& b2,
241 double& shearedZ,
242 const double sourceZ)
243 {
244 // Check out: ../../Resources/Computations/ComputeShearOnSlice.py
245 assert(IsValidShear(M_shear));
246
247 const double s11 = M_shear(0, 0);
248 const double s13 = M_shear(0, 2);
249 const double s14 = M_shear(0, 3);
250 const double s22 = M_shear(1, 1);
251 const double s23 = M_shear(1, 2);
252 const double s24 = M_shear(1, 3);
253 const double s43 = M_shear(3, 2);
254
255 double scaling = 1.0 / (s43 * sourceZ + 1.0);
256 shearedZ = sourceZ * scaling;
257
258 a11 = s11 * scaling;
259 a22 = s22 * scaling;
260
261 b1 = (s13 * sourceZ + s14) * scaling;
262 b2 = (s23 * sourceZ + s24) * scaling;
263 }
264
265
266 Matrix ShearWarpProjectiveTransform::CalibrateView(const Vector& camera,
267 const Vector& principalPoint,
268 double angle)
269 {
270 if (camera.size() != 3 ||
271 principalPoint.size() != 3)
272 {
273 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
274 }
275
276 const double sid = boost::numeric::ublas::norm_2(camera - principalPoint);
277
278 Matrix a;
279 GeometryToolbox::AlignVectorsWithRotation(a, camera - principalPoint,
280 LinearAlgebra::CreateVector(0, 0, -1));
281
282 Matrix r = LinearAlgebra::Product(GeometryToolbox::CreateRotationMatrixAlongZ(angle), a);
283
284 a = LinearAlgebra::ZeroMatrix(4, 4);
285 boost::numeric::ublas::subrange(a, 0, 3, 0, 3) = r;
286
287 const Vector v = LinearAlgebra::Product(r, -camera);
288 a(0, 3) = v[0];
289 a(1, 3) = v[1];
290 a(2, 3) = v[2];
291 a(3, 3) = 1;
292
293 Matrix perspective = LinearAlgebra::ZeroMatrix(4, 4);
294 // https://stackoverflow.com/questions/5267866/calculation-of-a-perspective-transformation-matrix
295 perspective(0, 0) = sid;
296 perspective(1, 1) = sid;
297 perspective(2, 2) = sid;
298 perspective(3, 2) = 1;
299
300 Matrix M_view = LinearAlgebra::Product(perspective, a);
301 assert(M_view.size1() == 4 &&
302 M_view.size2() == 4);
303
304 {
305 // Sanity checks
306 Vector p1 = LinearAlgebra::CreateVector(camera[0], camera[1], camera[2], 1.0);
307 Vector p2 = LinearAlgebra::CreateVector(principalPoint[0], principalPoint[1], principalPoint[2], 1.0);
308
309 Vector v1 = LinearAlgebra::Product(M_view, p1);
310 Vector v2 = LinearAlgebra::Product(M_view, p2);
311
312 if (!LinearAlgebra::IsCloseToZero(v1[3]) || // Must be mapped to singularity (w=0)
313 LinearAlgebra::IsCloseToZero(v2[3]))
314 {
315 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
316 }
317
318 // The principal point must be mapped to (0,0,z,1)
319 v2 /= v2[3];
320 if (!LinearAlgebra::IsCloseToZero(v2[0]) ||
321 !LinearAlgebra::IsCloseToZero(v2[1]))
322 {
323 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
324 }
325 }
326
327 return M_view;
328 }
329 }