ClangFormat: apply to source, most of intern
[blender.git] / intern / cycles / util / util_transform.cpp
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_projection.h"
50 #include "util/util_transform.h"
51
52 #include "util/util_boundbox.h"
53 #include "util/util_math.h"
54
55 CCL_NAMESPACE_BEGIN
56
57 /* Transform Inverse */
58
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
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
158 Transform transform_inverse(const Transform &tfm)
159 {
160   ProjectionTransform projection(tfm);
161   return projection_to_transform(projection_inverse(projection));
162 }
163
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
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
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   Transform R = M;
231   float norm;
232   int iteration = 0;
233
234   do {
235     Transform Rnext;
236     Transform Rit = transform_transposed_inverse(R);
237
238     for (int i = 0; i < 3; i++)
239       for (int j = 0; j < 4; j++)
240         Rnext[i][j] = 0.5f * (R[i][j] + Rit[i][j]);
241
242     norm = 0.0f;
243     for (int i = 0; i < 3; i++) {
244       norm = max(norm,
245                  fabsf(R[i][0] - Rnext[i][0]) + fabsf(R[i][1] - Rnext[i][1]) +
246                      fabsf(R[i][2] - Rnext[i][2]));
247     }
248
249     R = Rnext;
250     iteration++;
251   } while (iteration < 100 && norm > 1e-4f);
252
253   if (transform_negative_scale(R))
254     R = R * transform_scale(-1.0f, -1.0f, -1.0f);
255
256   decomp->x = transform_to_quat(R);
257
258   /* extract scale and pack it */
259   Transform scale = transform_inverse(R) * M;
260   decomp->y.w = scale.x.x;
261   decomp->z = make_float4(scale.x.y, scale.x.z, scale.y.x, scale.y.y);
262   decomp->w = make_float4(scale.y.z, scale.z.x, scale.z.y, scale.z.z);
263 }
264
265 void transform_motion_decompose(DecomposedTransform *decomp, const Transform *motion, size_t size)
266 {
267   for (size_t i = 0; i < size; i++) {
268     transform_decompose(decomp + i, motion + i);
269
270     if (i > 0) {
271       /* Ensure rotation around shortest angle, negated quaternions are the same
272        * but this means we don't have to do the check in quat_interpolate */
273       if (dot(decomp[i - 1].x, decomp[i].x) < 0.0f)
274         decomp[i - 1].x = -decomp[i - 1].x;
275     }
276   }
277 }
278
279 Transform transform_from_viewplane(BoundBox2D &viewplane)
280 {
281   return transform_scale(1.0f / (viewplane.right - viewplane.left),
282                          1.0f / (viewplane.top - viewplane.bottom),
283                          1.0f) *
284          transform_translate(-viewplane.left, -viewplane.bottom, 0.0f);
285 }
286
287 CCL_NAMESPACE_END