Blender  V3.3
transform.cpp
Go to the documentation of this file.
1 /* SPDX-License-Identifier: Apache-2.0
2  * Copyright 2011-2022 Blender Foundation */
3 
4 #include "util/transform.h"
5 #include "util/projection.h"
6 
7 #include "util/boundbox.h"
8 #include "util/math.h"
9 
11 
12 /* Transform Inverse */
13 
14 static bool projection_matrix4_inverse(float R[][4], float M[][4])
15 {
16  /* SPDX-License-Identifier: BSD-3-Clause
17  * Adapted from code:
18  * Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
19  * Digital Ltd. LLC. All rights reserved. */
20 
21  /* forward elimination */
22  for (int i = 0; i < 4; i++) {
23  int pivot = i;
24  float pivotsize = M[i][i];
25 
26  if (pivotsize < 0)
27  pivotsize = -pivotsize;
28 
29  for (int j = i + 1; j < 4; j++) {
30  float tmp = M[j][i];
31 
32  if (tmp < 0)
33  tmp = -tmp;
34 
35  if (tmp > pivotsize) {
36  pivot = j;
37  pivotsize = tmp;
38  }
39  }
40 
41  if (UNLIKELY(pivotsize == 0.0f))
42  return false;
43 
44  if (pivot != i) {
45  for (int j = 0; j < 4; j++) {
46  float tmp;
47 
48  tmp = M[i][j];
49  M[i][j] = M[pivot][j];
50  M[pivot][j] = tmp;
51 
52  tmp = R[i][j];
53  R[i][j] = R[pivot][j];
54  R[pivot][j] = tmp;
55  }
56  }
57 
58  for (int j = i + 1; j < 4; j++) {
59  float f = M[j][i] / M[i][i];
60 
61  for (int k = 0; k < 4; k++) {
62  M[j][k] -= f * M[i][k];
63  R[j][k] -= f * R[i][k];
64  }
65  }
66  }
67 
68  /* backward substitution */
69  for (int i = 3; i >= 0; --i) {
70  float f;
71 
72  if (UNLIKELY((f = M[i][i]) == 0.0f))
73  return false;
74 
75  for (int j = 0; j < 4; j++) {
76  M[i][j] /= f;
77  R[i][j] /= f;
78  }
79 
80  for (int j = 0; j < i; j++) {
81  f = M[j][i];
82 
83  for (int k = 0; k < 4; k++) {
84  M[j][k] -= f * M[i][k];
85  R[j][k] -= f * R[i][k];
86  }
87  }
88  }
89 
90  return true;
91 }
92 
94 {
96  float M[4][4], R[4][4];
97 
98  memcpy(R, &tfmR, sizeof(R));
99  memcpy(M, &tfm, sizeof(M));
100 
102  return projection_identity();
103  }
104 
105  memcpy(&tfmR, R, sizeof(R));
106 
107  return tfmR;
108 }
109 
111 {
112  ProjectionTransform iprojection(transform_inverse(tfm));
113  return projection_to_transform(projection_transpose(iprojection));
114 }
115 
116 /* Motion Transform */
117 
119 {
120  double trace = (double)(tfm[0][0] + tfm[1][1] + tfm[2][2]);
121  float4 qt;
122 
123  if (trace > 0.0) {
124  double s = sqrt(trace + 1.0);
125 
126  qt.w = (float)(s / 2.0);
127  s = 0.5 / s;
128 
129  qt.x = (float)((double)(tfm[2][1] - tfm[1][2]) * s);
130  qt.y = (float)((double)(tfm[0][2] - tfm[2][0]) * s);
131  qt.z = (float)((double)(tfm[1][0] - tfm[0][1]) * s);
132  }
133  else {
134  int i = 0;
135 
136  if (tfm[1][1] > tfm[i][i])
137  i = 1;
138  if (tfm[2][2] > tfm[i][i])
139  i = 2;
140 
141  int j = (i + 1) % 3;
142  int k = (j + 1) % 3;
143 
144  double s = sqrt((double)(tfm[i][i] - (tfm[j][j] + tfm[k][k])) + 1.0);
145 
146  double q[3];
147  q[i] = s * 0.5;
148  if (s != 0.0)
149  s = 0.5 / s;
150 
151  double w = (double)(tfm[k][j] - tfm[j][k]) * s;
152  q[j] = (double)(tfm[j][i] + tfm[i][j]) * s;
153  q[k] = (double)(tfm[k][i] + tfm[i][k]) * s;
154 
155  qt.x = (float)q[0];
156  qt.y = (float)q[1];
157  qt.z = (float)q[2];
158  qt.w = (float)w;
159  }
160 
161  return qt;
162 }
163 
164 static void transform_decompose(DecomposedTransform *decomp, const Transform *tfm)
165 {
166  /* extract translation */
167  decomp->y = make_float4(tfm->x.w, tfm->y.w, tfm->z.w, 0.0f);
168 
169  /* extract rotation */
170  Transform M = *tfm;
171  M.x.w = 0.0f;
172  M.y.w = 0.0f;
173  M.z.w = 0.0f;
174 
175 #if 0
176  Transform R = M;
177  float norm;
178  int iteration = 0;
179 
180  do {
181  Transform Rnext;
183 
184  for (int i = 0; i < 3; i++)
185  for (int j = 0; j < 4; j++)
186  Rnext[i][j] = 0.5f * (R[i][j] + Rit[i][j]);
187 
188  norm = 0.0f;
189  for (int i = 0; i < 3; i++) {
190  norm = max(norm,
191  fabsf(R[i][0] - Rnext[i][0]) + fabsf(R[i][1] - Rnext[i][1]) +
192  fabsf(R[i][2] - Rnext[i][2]));
193  }
194 
195  R = Rnext;
196  iteration++;
197  } while (iteration < 100 && norm > 1e-4f);
198 
200  R = R * transform_scale(-1.0f, -1.0f, -1.0f);
201 
202  decomp->x = transform_to_quat(R);
203 
204  /* extract scale and pack it */
205  Transform scale = transform_inverse(R) * M;
206  decomp->y.w = scale.x.x;
207  decomp->z = make_float4(scale.x.y, scale.x.z, scale.y.x, scale.y.y);
208  decomp->w = make_float4(scale.y.z, scale.z.x, scale.z.y, scale.z.z);
209 #else
210  float3 colx = transform_get_column(&M, 0);
211  float3 coly = transform_get_column(&M, 1);
212  float3 colz = transform_get_column(&M, 2);
213 
214  /* extract scale and shear first */
215  float3 scale, shear;
216  scale.x = len(colx);
217  colx = safe_divide(colx, scale.x);
218  shear.z = dot(colx, coly);
219  coly -= shear.z * colx;
220  scale.y = len(coly);
221  coly = safe_divide(coly, scale.y);
222  shear.y = dot(colx, colz);
223  colz -= shear.y * colx;
224  shear.x = dot(coly, colz);
225  colz -= shear.x * coly;
226  scale.z = len(colz);
227  colz = safe_divide(colz, scale.z);
228 
229  transform_set_column(&M, 0, colx);
230  transform_set_column(&M, 1, coly);
231  transform_set_column(&M, 2, colz);
232 
234  scale *= -1.0f;
235  M = M * transform_scale(-1.0f, -1.0f, -1.0f);
236  }
237 
238  decomp->x = transform_to_quat(M);
239 
240  decomp->y.w = scale.x;
241  decomp->z = make_float4(shear.z, shear.y, 0.0f, scale.y);
242  decomp->w = make_float4(shear.x, 0.0f, 0.0f, scale.z);
243 #endif
244 }
245 
246 void transform_motion_decompose(DecomposedTransform *decomp, const Transform *motion, size_t size)
247 {
248  /* Decompose and correct rotation. */
249  for (size_t i = 0; i < size; i++) {
250  transform_decompose(decomp + i, motion + i);
251 
252  if (i > 0) {
253  /* Ensure rotation around shortest angle, negated quaternions are the same
254  * but this means we don't have to do the check in quat_interpolate */
255  if (dot(decomp[i - 1].x, decomp[i].x) < 0.0f)
256  decomp[i].x = -decomp[i].x;
257  }
258  }
259 
260  /* Copy rotation to decomposed transform where scale is degenerate. This avoids weird object
261  * rotation interpolation when the scale goes to 0 for a time step.
262  *
263  * Note that this is very simple and naive implementation, which only deals with degenerated
264  * scale happening only on one frame. It is possible to improve it further by interpolating
265  * rotation into s degenerated range using rotation from time-steps from adjacent non-degenerated
266  * time steps. */
267  for (size_t i = 0; i < size; i++) {
268  const float3 scale = make_float3(decomp[i].y.w, decomp[i].z.w, decomp[i].w.w);
269  if (!is_zero(scale)) {
270  continue;
271  }
272 
273  if (i > 0) {
274  decomp[i].x = decomp[i - 1].x;
275  }
276  else if (i < size - 1) {
277  decomp[i].x = decomp[i + 1].x;
278  }
279  }
280 }
281 
283 {
284  return transform_scale(1.0f / (viewplane.right - viewplane.left),
285  1.0f / (viewplane.top - viewplane.bottom),
286  1.0f) *
287  transform_translate(-viewplane.left, -viewplane.bottom, 0.0f);
288 }
289 
typedef float(TangentPoint)[2]
sqrt(x)+1/max(0
#define UNLIKELY(x)
typedef double(DMatrix)[4][4]
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint y
float float4[4]
static DBVT_INLINE btScalar size(const btDbvtVolume &a)
Definition: btDbvt.cpp:52
SIMD_FORCE_INLINE const btScalar & w() const
Return the w value.
Definition: btQuadWord.h:119
SIMD_FORCE_INLINE btScalar norm() const
Return the norm (length) of the vector.
Definition: btVector3.h:263
float bottom
Definition: boundbox.h:195
float top
Definition: boundbox.h:196
float right
Definition: boundbox.h:194
float left
Definition: boundbox.h:193
#define CCL_NAMESPACE_END
Definition: cuda/compat.h:9
ccl_device_inline Transform projection_to_transform(const ProjectionTransform &a)
ccl_device_inline ProjectionTransform projection_identity()
ccl_device_inline ProjectionTransform projection_transpose(const ProjectionTransform &a)
int len
Definition: draw_manager.c:108
ccl_device_inline void transform_set_column(Transform *t, int column, float3 value)
ccl_device_inline float3 transform_get_column(const Transform *t, int column)
ccl_device_inline bool transform_negative_scale(const Transform &tfm)
ccl_device_inline Transform transform_translate(float3 t)
ccl_device_inline Transform transform_inverse(const Transform tfm)
ccl_device_inline Transform transform_scale(float3 s)
#define M
#define R
#define make_float4(x, y, z, w)
Definition: metal/compat.h:205
#define fabsf(x)
Definition: metal/compat.h:219
#define make_float3(x, y, z)
Definition: metal/compat.h:204
T dot(const vec_base< T, Size > &a, const vec_base< T, Size > &b)
bool is_zero(const T &a)
T safe_divide(const T &a, const T &b)
float z
float y
float x
Transform transform_transposed_inverse(const Transform &tfm)
Definition: transform.cpp:110
static void transform_decompose(DecomposedTransform *decomp, const Transform *tfm)
Definition: transform.cpp:164
float4 transform_to_quat(const Transform &tfm)
Definition: transform.cpp:118
void transform_motion_decompose(DecomposedTransform *decomp, const Transform *motion, size_t size)
Definition: transform.cpp:246
Transform transform_from_viewplane(BoundBox2D &viewplane)
Definition: transform.cpp:282
static CCL_NAMESPACE_BEGIN bool projection_matrix4_inverse(float R[][4], float M[][4])
Definition: transform.cpp:14
ProjectionTransform projection_inverse(const ProjectionTransform &tfm)
Definition: transform.cpp:93
float max