1 /*
2 * Copyright 2011-2013 Blender Foundation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 /*
18 * Adapted from code with license:
19 *
20 * Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
21 * Digital Ltd. LLC. All rights reserved.
22 *
23 * Redistribution and use in source and binary forms, with or without
24 * modification, are permitted provided that the following conditions are
25 * met:
26 * * Redistributions of source code must retain the above copyright
27 * notice, this list of conditions and the following disclaimer.
28 * * Redistributions in binary form must reproduce the above copyright
29 * notice, this list of conditions and the following disclaimer in
30 * the documentation and/or other materials provided with the
31 * distribution.
32 * * Neither the name of Industrial Light & Magic nor the names of its
33 * contributors may be used to endorse or promote products derived
34 * from this software without specific prior written permission.
35 *
36 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
37 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
38 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
39 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
40 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
41 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
42 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
43 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
45 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
46 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
47 */
48
49 #include "util/util_transform.h"
50 #include "util/util_projection.h"
51
52 #include "util/util_boundbox.h"
53 #include "util/util_math.h"
54
55 CCL_NAMESPACE_BEGIN
56
57 /* Transform Inverse */
58
transform_matrix4_gj_inverse(float R[][4],float M[][4])59 static bool transform_matrix4_gj_inverse(float R[][4], float M[][4])
60 {
61 /* forward elimination */
62 for (int i = 0; i < 4; i++) {
63 int pivot = i;
64 float pivotsize = M[i][i];
65
66 if (pivotsize < 0)
67 pivotsize = -pivotsize;
68
69 for (int j = i + 1; j < 4; j++) {
70 float tmp = M[j][i];
71
72 if (tmp < 0)
73 tmp = -tmp;
74
75 if (tmp > pivotsize) {
76 pivot = j;
77 pivotsize = tmp;
78 }
79 }
80
81 if (UNLIKELY(pivotsize == 0.0f))
82 return false;
83
84 if (pivot != i) {
85 for (int j = 0; j < 4; j++) {
86 float tmp;
87
88 tmp = M[i][j];
89 M[i][j] = M[pivot][j];
90 M[pivot][j] = tmp;
91
92 tmp = R[i][j];
93 R[i][j] = R[pivot][j];
94 R[pivot][j] = tmp;
95 }
96 }
97
98 for (int j = i + 1; j < 4; j++) {
99 float f = M[j][i] / M[i][i];
100
101 for (int k = 0; k < 4; k++) {
102 M[j][k] -= f * M[i][k];
103 R[j][k] -= f * R[i][k];
104 }
105 }
106 }
107
108 /* backward substitution */
109 for (int i = 3; i >= 0; --i) {
110 float f;
111
112 if (UNLIKELY((f = M[i][i]) == 0.0f))
113 return false;
114
115 for (int j = 0; j < 4; j++) {
116 M[i][j] /= f;
117 R[i][j] /= f;
118 }
119
120 for (int j = 0; j < i; j++) {
121 f = M[j][i];
122
123 for (int k = 0; k < 4; k++) {
124 M[j][k] -= f * M[i][k];
125 R[j][k] -= f * R[i][k];
126 }
127 }
128 }
129
130 return true;
131 }
132
projection_inverse(const ProjectionTransform & tfm)133 ProjectionTransform projection_inverse(const ProjectionTransform &tfm)
134 {
135 ProjectionTransform tfmR = projection_identity();
136 float M[4][4], R[4][4];
137
138 memcpy(R, &tfmR, sizeof(R));
139 memcpy(M, &tfm, sizeof(M));
140
141 if (UNLIKELY(!transform_matrix4_gj_inverse(R, M))) {
142 /* matrix is degenerate (e.g. 0 scale on some axis), ideally we should
143 * never be in this situation, but try to invert it anyway with tweak */
144 M[0][0] += 1e-8f;
145 M[1][1] += 1e-8f;
146 M[2][2] += 1e-8f;
147
148 if (UNLIKELY(!transform_matrix4_gj_inverse(R, M))) {
149 return projection_identity();
150 }
151 }
152
153 memcpy(&tfmR, R, sizeof(R));
154
155 return tfmR;
156 }
157
transform_inverse(const Transform & tfm)158 Transform transform_inverse(const Transform &tfm)
159 {
160 ProjectionTransform projection(tfm);
161 return projection_to_transform(projection_inverse(projection));
162 }
163
transform_transposed_inverse(const Transform & tfm)164 Transform transform_transposed_inverse(const Transform &tfm)
165 {
166 ProjectionTransform projection(tfm);
167 ProjectionTransform iprojection = projection_inverse(projection);
168 return projection_to_transform(projection_transpose(iprojection));
169 }
170
171 /* Motion Transform */
172
transform_to_quat(const Transform & tfm)173 float4 transform_to_quat(const Transform &tfm)
174 {
175 double trace = (double)(tfm[0][0] + tfm[1][1] + tfm[2][2]);
176 float4 qt;
177
178 if (trace > 0.0) {
179 double s = sqrt(trace + 1.0);
180
181 qt.w = (float)(s / 2.0);
182 s = 0.5 / s;
183
184 qt.x = (float)((double)(tfm[2][1] - tfm[1][2]) * s);
185 qt.y = (float)((double)(tfm[0][2] - tfm[2][0]) * s);
186 qt.z = (float)((double)(tfm[1][0] - tfm[0][1]) * s);
187 }
188 else {
189 int i = 0;
190
191 if (tfm[1][1] > tfm[i][i])
192 i = 1;
193 if (tfm[2][2] > tfm[i][i])
194 i = 2;
195
196 int j = (i + 1) % 3;
197 int k = (j + 1) % 3;
198
199 double s = sqrt((double)(tfm[i][i] - (tfm[j][j] + tfm[k][k])) + 1.0);
200
201 double q[3];
202 q[i] = s * 0.5;
203 if (s != 0.0)
204 s = 0.5 / s;
205
206 double w = (double)(tfm[k][j] - tfm[j][k]) * s;
207 q[j] = (double)(tfm[j][i] + tfm[i][j]) * s;
208 q[k] = (double)(tfm[k][i] + tfm[i][k]) * s;
209
210 qt.x = (float)q[0];
211 qt.y = (float)q[1];
212 qt.z = (float)q[2];
213 qt.w = (float)w;
214 }
215
216 return qt;
217 }
218
transform_decompose(DecomposedTransform * decomp,const Transform * tfm)219 static void transform_decompose(DecomposedTransform *decomp, const Transform *tfm)
220 {
221 /* extract translation */
222 decomp->y = make_float4(tfm->x.w, tfm->y.w, tfm->z.w, 0.0f);
223
224 /* extract rotation */
225 Transform M = *tfm;
226 M.x.w = 0.0f;
227 M.y.w = 0.0f;
228 M.z.w = 0.0f;
229
230 #if 0
231 Transform R = M;
232 float norm;
233 int iteration = 0;
234
235 do {
236 Transform Rnext;
237 Transform Rit = transform_transposed_inverse(R);
238
239 for (int i = 0; i < 3; i++)
240 for (int j = 0; j < 4; j++)
241 Rnext[i][j] = 0.5f * (R[i][j] + Rit[i][j]);
242
243 norm = 0.0f;
244 for (int i = 0; i < 3; i++) {
245 norm = max(norm,
246 fabsf(R[i][0] - Rnext[i][0]) + fabsf(R[i][1] - Rnext[i][1]) +
247 fabsf(R[i][2] - Rnext[i][2]));
248 }
249
250 R = Rnext;
251 iteration++;
252 } while (iteration < 100 && norm > 1e-4f);
253
254 if (transform_negative_scale(R))
255 R = R * transform_scale(-1.0f, -1.0f, -1.0f);
256
257 decomp->x = transform_to_quat(R);
258
259 /* extract scale and pack it */
260 Transform scale = transform_inverse(R) * M;
261 decomp->y.w = scale.x.x;
262 decomp->z = make_float4(scale.x.y, scale.x.z, scale.y.x, scale.y.y);
263 decomp->w = make_float4(scale.y.z, scale.z.x, scale.z.y, scale.z.z);
264 #else
265 float3 colx = transform_get_column(&M, 0);
266 float3 coly = transform_get_column(&M, 1);
267 float3 colz = transform_get_column(&M, 2);
268
269 /* extract scale and shear first */
270 float3 scale, shear;
271 scale.x = len(colx);
272 colx = safe_divide_float3_float(colx, scale.x);
273 shear.z = dot(colx, coly);
274 coly -= shear.z * colx;
275 scale.y = len(coly);
276 coly = safe_divide_float3_float(coly, scale.y);
277 shear.y = dot(colx, colz);
278 colz -= shear.y * colx;
279 shear.x = dot(coly, colz);
280 colz -= shear.x * coly;
281 scale.z = len(colz);
282 colz = safe_divide_float3_float(colz, scale.z);
283
284 transform_set_column(&M, 0, colx);
285 transform_set_column(&M, 1, coly);
286 transform_set_column(&M, 2, colz);
287
288 if (transform_negative_scale(M)) {
289 scale *= -1.0f;
290 M = M * transform_scale(-1.0f, -1.0f, -1.0f);
291 }
292
293 decomp->x = transform_to_quat(M);
294
295 decomp->y.w = scale.x;
296 decomp->z = make_float4(shear.z, shear.y, 0.0f, scale.y);
297 decomp->w = make_float4(shear.x, 0.0f, 0.0f, scale.z);
298 #endif
299 }
300
transform_motion_decompose(DecomposedTransform * decomp,const Transform * motion,size_t size)301 void transform_motion_decompose(DecomposedTransform *decomp, const Transform *motion, size_t size)
302 {
303 /* Decompose and correct rotation. */
304 for (size_t i = 0; i < size; i++) {
305 transform_decompose(decomp + i, motion + i);
306
307 if (i > 0) {
308 /* Ensure rotation around shortest angle, negated quaternions are the same
309 * but this means we don't have to do the check in quat_interpolate */
310 if (dot(decomp[i - 1].x, decomp[i].x) < 0.0f)
311 decomp[i].x = -decomp[i].x;
312 }
313 }
314
315 /* Copy rotation to decomposed transform where scale is degenerate. This avoids weird object
316 * rotation interpolation when the scale goes to 0 for a time step.
317 *
318 * Note that this is very simple and naive implementation, which only deals with degenerated
319 * scale happening only on one frame. It is possible to improve it further by interpolating
320 * rotation into s degenerated range using rotation from time-steps from adjacent non-degenerated
321 * time steps. */
322 for (size_t i = 0; i < size; i++) {
323 const float3 scale = make_float3(decomp[i].y.w, decomp[i].z.w, decomp[i].w.w);
324 if (!is_zero(scale)) {
325 continue;
326 }
327
328 if (i > 0) {
329 decomp[i].x = decomp[i - 1].x;
330 }
331 else if (i < size - 1) {
332 decomp[i].x = decomp[i + 1].x;
333 }
334 }
335 }
336
transform_from_viewplane(BoundBox2D & viewplane)337 Transform transform_from_viewplane(BoundBox2D &viewplane)
338 {
339 return transform_scale(1.0f / (viewplane.right - viewplane.left),
340 1.0f / (viewplane.top - viewplane.bottom),
341 1.0f) *
342 transform_translate(-viewplane.left, -viewplane.bottom, 0.0f);
343 }
344
345 CCL_NAMESPACE_END
346