Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 191:46cb2eedc2e0 wasm
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 } |