1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; -*- */
2 /* ***** BEGIN LICENSE BLOCK *****
3 * This file is part of openfx-supportext <https://github.com/devernay/openfx-supportext>,
4 * Copyright (C) 2013-2018 INRIA
5 *
6 * openfx-supportext is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * openfx-supportext is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with openfx-supportext. If not, see <http://www.gnu.org/licenses/gpl-2.0.html>
18 * ***** END LICENSE BLOCK ***** */
19
20 /*
21 * OFX Matrix utilities.
22 */
23
24 #ifndef openfx_supportext_ofxsMatrix2D_h
25 #define openfx_supportext_ofxsMatrix2D_h
26
27 #include <cmath>
28 #include <cfloat>
29
30 namespace OFX {
31 // NEVER define a variable/constant in a header (said it 100 times already)
32 // an inline function is OK
33 inline double
ofxsPi()34 ofxsPi()
35 {
36 return 3.14159265358979323846264338327950288419717;
37 }
38
39 /**
40 * @brief A simple 3 * 3 matrix class layed out as such:
41 * a b c
42 * d e f
43 * g h i
44 **/
45
46 inline double
ofxsToDegrees(double rad)47 ofxsToDegrees(double rad)
48 {
49 rad = rad * 180.0 / ofxsPi();
50
51 return rad;
52 }
53
54 inline double
ofxsToRadians(double deg)55 ofxsToRadians(double deg)
56 {
57 deg = deg * ofxsPi() / 180.0;
58
59 return deg;
60 }
61
62 struct Point3D
63 {
64 double x, y, z;
65
Point3DPoint3D66 Point3D()
67 : x(0), y(0), z(0)
68 {
69 }
70
Point3DPoint3D71 Point3D(double x,
72 double y,
73 double z)
74 : x(x), y(y), z(z)
75 {
76 }
77
Point3DPoint3D78 Point3D(const Point3D & p)
79 : x(p.x), y(p.y), z(p.z)
80 {
81 }
82
83 Point3D& operator=(const Point3D & p)
84 {
85 x = p.x;
86 y = p.y;
87 z = p.z;
88 return *this;
89 }
90
91 bool operator==(const Point3D & other)
92 {
93 return x == other.x && y == other.y && z == other.z;
94 }
95 };
96
97
98 /**
99 * \brief Compute the cross-product of two vectors
100 *
101 */
102 inline OFX::Point3D
crossprod(const OFX::Point3D & a,const OFX::Point3D & b)103 crossprod(const OFX::Point3D & a,
104 const OFX::Point3D & b)
105 {
106 OFX::Point3D c;
107
108 c.x = a.y * b.z - a.z * b.y;
109 c.y = a.z * b.x - a.x * b.z;
110 c.z = a.x * b.y - a.y * b.x;
111
112 return c;
113 }
114
115 struct Point4D
116 {
117 double x, y, z, w;
118
Point4DPoint4D119 Point4D()
120 : x(0), y(0), z(0), w(0)
121 {
122 }
123
Point4DPoint4D124 Point4D(double x,
125 double y,
126 double z,
127 double w)
128 : x(x), y(y), z(z), w(w)
129 {
130 }
131
Point4DPoint4D132 Point4D(const Point4D & o)
133 : x(o.x), y(o.y), z(o.z), w(o.w)
134 {
135 }
136
operatorPoint4D137 double & operator() (int i)
138 {
139 switch (i) {
140 case 0:
141
142 return x;
143 case 1:
144
145 return y;
146 case 2:
147
148 return z;
149 case 3:
150
151 return w;
152 default:
153 assert(false);
154
155 return x;
156 }
157 }
158
operatorPoint4D159 double operator() (int i) const
160 {
161 switch (i) {
162 case 0:
163
164 return x;
165 case 1:
166
167 return y;
168 case 2:
169
170 return z;
171 case 3:
172
173 return w;
174 default:
175 assert(false);
176
177 return x;
178 }
179 }
180
181 bool operator==(const Point4D & o)
182 {
183 return x == o.x && y == o.y && z == o.z && w == o.w;
184 }
185 };
186
187 struct Matrix3x3
188 {
189 double m[3*3];
190
Matrix3x3Matrix3x3191 Matrix3x3()
192 {
193 std::fill(m, m + 3*3, 0.);
194 }
195
Matrix3x3Matrix3x3196 Matrix3x3(double a,
197 double b,
198 double c,
199 double d,
200 double e,
201 double f,
202 double g,
203 double h,
204 double i)
205 {
206 m[0] = a; m[1] = b; m[2] = c;
207 m[3] = d; m[4] = e; m[5] = f;
208 m[6] = g; m[7] = h; m[8] = i;
209 }
210
Matrix3x3Matrix3x3211 Matrix3x3(const Matrix3x3 & mat)
212 {
213 std::copy(mat.m, mat.m + 3*3, m);
214 }
215
216 /// Contruct from columns
Matrix3x3Matrix3x3217 Matrix3x3(const OFX::Point3D &m0,
218 const OFX::Point3D &m1,
219 const OFX::Point3D &m2)
220 {
221 m[0] = m0.x; m[1] = m1.x; m[2] = m2.x;
222 m[3] = m0.y; m[4] = m1.y; m[5] = m2.y;
223 m[6] = m0.z; m[7] = m1.z; m[8] = m2.z;
224 }
225
226 Matrix3x3 & operator=(const Matrix3x3 & mat)
227 {
228 std::copy(mat.m, mat.m + 3*3, m); return *this;
229 }
230
isIdentityMatrix3x3231 bool isIdentity() const
232 {
233 return m[0] == 1 && m[1] == 0 && m[2] == 0 &&
234 m[3] == 0 && m[4] == 1 && m[5] == 0 &&
235 m[6] == 0 && m[7] == 0 && m[8] == 1;
236 }
237
operatorMatrix3x3238 double & operator()(int row,
239 int col)
240 {
241 assert(row >= 0 && row < 3 && col >= 0 && col < 3);
242
243 return m[row * 3 + col];
244 }
245
operatorMatrix3x3246 double operator()(int row,
247 int col) const
248 {
249 assert(row >= 0 && row < 3 && col >= 0 && col < 3);
250
251 return m[row * 3 + col];
252 }
253
254 Matrix3x3 operator*(const Matrix3x3 & m2) const
255 {
256 return Matrix3x3(m[0] * m2.m[0] + m[1] * m2.m[3] + m[2] * m2.m[6],
257 m[0] * m2.m[1] + m[1] * m2.m[4] + m[2] * m2.m[7],
258 m[0] * m2.m[2] + m[1] * m2.m[5] + m[2] * m2.m[8],
259 m[3] * m2.m[0] + m[4] * m2.m[3] + m[5] * m2.m[6],
260 m[3] * m2.m[1] + m[4] * m2.m[4] + m[5] * m2.m[7],
261 m[3] * m2.m[2] + m[4] * m2.m[5] + m[5] * m2.m[8],
262 m[6] * m2.m[0] + m[7] * m2.m[3] + m[8] * m2.m[6],
263 m[6] * m2.m[1] + m[7] * m2.m[4] + m[8] * m2.m[7],
264 m[6] * m2.m[2] + m[7] * m2.m[5] + m[8] * m2.m[8]);
265 }
266
267 Point3D operator*(const Point3D & p) const
268 {
269 Point3D ret;
270
271 ret.x = m[0] * p.x + m[1] * p.y + m[2] * p.z;
272 ret.y = m[3] * p.x + m[4] * p.y + m[5] * p.z;
273 ret.z = m[6] * p.x + m[7] * p.y + m[8] * p.z;
274
275 return ret;
276 }
277
278 Matrix3x3 operator *(double s)
279 {
280 return Matrix3x3(m[0] * s,
281 m[1] * s,
282 m[2] * s,
283 m[3] * s,
284 m[4] * s,
285 m[5] * s,
286 m[6] * s,
287 m[7] * s,
288 m[8] * s);
289 }
290
291 Matrix3x3& operator *=(double s)
292 {
293 for (int i=0; i < 9; ++i) {
294 m[i] *= s;
295 }
296 return *this;
297 }
298
299
determinantMatrix3x3300 double determinant() const
301 {
302 return m[0] * (m[4] * m[8] - m[7] * m[5])
303 - m[1] * (m[3] * m[8] - m[6] * m[5])
304 + m[2] * (m[3] * m[7] - m[6] * m[4]);
305 }
306
scaledAdjointMatrix3x3307 Matrix3x3 scaledAdjoint(double s) const
308 {
309 Matrix3x3 ret;
310
311 ret.m[0] = (s) * (m[4] * m[8] - m[7] * m[5]);
312 ret.m[3] = (s) * (m[5] * m[6] - m[3] * m[8]);
313 ret.m[6] = (s) * (m[3] * m[7] - m[4] * m[6]);
314
315 ret.m[1] = (s) * (m[2] * m[7] - m[1] * m[8]);
316 ret.m[4] = (s) * (m[0] * m[8] - m[2] * m[6]);
317 ret.m[7] = (s) * (m[1] * m[6] - m[0] * m[7]);
318
319 ret.m[2] = (s) * (m[1] * m[5] - m[2] * m[4]);
320 ret.m[5] = (s) * (m[2] * m[3] - m[0] * m[5]);
321 ret.m[8] = (s) * (m[0] * m[4] - m[1] * m[3]);
322
323 return ret;
324 }
325
inverseMatrix3x3326 bool inverse(Matrix3x3* invOut) const
327 {
328 double inv[9], det;
329 int i;
330
331 inv[0] = (m[4] * m[8] - m[7] * m[5]);
332 inv[3] = (m[5] * m[6] - m[3] * m[8]);
333 inv[6] = (m[3] * m[7] - m[4] * m[6]);
334
335 det = m[0] * inv[0] + m[1] * inv[3] + m[2] * inv[6];
336
337 if (det == 0) {
338 return false;
339 }
340
341 inv[1] = (m[2] * m[7] - m[1] * m[8]);
342 inv[4] = (m[0] * m[8] - m[2] * m[6]);
343 inv[7] = (m[1] * m[6] - m[0] * m[7]);
344
345 inv[2] = (m[1] * m[5] - m[2] * m[4]);
346 inv[5] = (m[2] * m[3] - m[0] * m[5]);
347 inv[8] = (m[0] * m[4] - m[1] * m[3]);
348
349 det = 1.0 / det;
350
351 for (i = 0; i < 9; ++i) {
352 invOut->m[i] = inv[i] * det;
353 }
354 return true;
355
356 /* old version
357 double det = determinant();
358 if (det == 0) {
359 return false;
360 }
361 *invOut = scaledAdjoint( 1. / det );
362 return true;
363 */
364 }
365
setIdentityMatrix3x3366 void setIdentity()
367 {
368 m[0] = 1; m[1] = 0; m[2] = 0;
369 m[3] = 0; m[4] = 1; m[5] = 0;
370 m[6] = 0; m[7] = 0; m[8] = 1;
371 }
372
373 /**
374 * \brief Compute a homography from 4 points correspondences
375 * \param p1 source point
376 * \param p2 source point
377 * \param p3 source point
378 * \param p4 source point
379 * \param q1 target point
380 * \param q2 target point
381 * \param q3 target point
382 * \param q4 target point
383 * \return the homography matrix that maps pi's to qi's
384 *
385 Using four point-correspondences pi ↔ pi^, we can set up an equation system to solve for the homography matrix H.
386 An algorithm to obtain these parameters requiring only the inversion of a 3 × 3 equation system is as follows.
387 From the four point-correspondences pi ↔ pi^ with (i ∈ {1, 2, 3, 4}),
388 compute h1 = (p1 × p2 ) × (p3 × p4 ), h2 = (p1 × p3 ) × (p2 × p4 ), h3 = (p1 × p4 ) × (p2 × p3 ).
389 Also compute h1^ , h2^ , h3^ using the same principle from the points pi^.
390 Now, the homography matrix H can be obtained easily from
391 H · [h1 h2 h3] = [h1^ h2^ h3^],
392 which only requires the inversion of the matrix [h1 h2 h3].
393
394 Algo from:
395 http://www.dirk-farin.net/publications/phd/text/AB_EfficientComputationOfHomographiesFromFourCorrespondences.pdf
396 */
setHomographyFromFourPointsMatrix3x3397 bool setHomographyFromFourPoints(const OFX::Point3D &p1,
398 const OFX::Point3D &p2,
399 const OFX::Point3D &p3,
400 const OFX::Point3D &p4,
401 const OFX::Point3D &q1,
402 const OFX::Point3D &q2,
403 const OFX::Point3D &q3,
404 const OFX::Point3D &q4)
405 {
406 OFX::Matrix3x3 invHp;
407 OFX::Matrix3x3 Hp( crossprod( crossprod(p1, p2), crossprod(p3, p4) ),
408 crossprod( crossprod(p1, p3), crossprod(p2, p4) ),
409 crossprod( crossprod(p1, p4), crossprod(p2, p3) ) );
410
411 if ( !Hp.inverse(&invHp) ) {
412 return false;
413 }
414 OFX::Matrix3x3 Hq( crossprod( crossprod(q1, q2), crossprod(q3, q4) ),
415 crossprod( crossprod(q1, q3), crossprod(q2, q4) ),
416 crossprod( crossprod(q1, q4), crossprod(q2, q3) ) );
417 double detHq = Hq.determinant();
418 if (detHq == 0.) {
419 return false;
420 }
421 *this = Hq * invHp;
422
423 return true;
424 }
425
setAffineFromThreePointsMatrix3x3426 bool setAffineFromThreePoints(const OFX::Point3D &p1,
427 const OFX::Point3D &p2,
428 const OFX::Point3D &p3,
429 const OFX::Point3D &q1,
430 const OFX::Point3D &q2,
431 const OFX::Point3D &q3)
432 {
433 OFX::Matrix3x3 invHp;
434 OFX::Matrix3x3 Hp(p1, p2, p3);
435 if ( !Hp.inverse(&invHp) ) {
436 return false;
437 }
438 OFX::Matrix3x3 Hq(q1, q2, q3);
439 double detHq = Hq.determinant();
440 if (detHq == 0.) {
441 return false;
442 }
443 *this = Hq * invHp;
444
445 return true;
446 }
447
setSimilarityFromTwoPointsMatrix3x3448 bool setSimilarityFromTwoPoints(const OFX::Point3D &p1,
449 const OFX::Point3D &p2,
450 const OFX::Point3D &q1,
451 const OFX::Point3D &q2)
452 {
453 // Generate a third point so that p1p3 is orthogonal to p1p2, and compute the affine transform
454 OFX::Point3D p3, q3;
455
456 p3.x = p1.x - (p2.y - p1.y);
457 p3.y = p1.y + (p2.x - p1.x);
458 p3.z = 1.;
459 q3.x = q1.x - (q2.y - q1.y);
460 q3.y = q1.y + (q2.x - q1.x);
461 q3.z = 1.;
462
463 return setAffineFromThreePoints(p1, p2, p3, q1, q2, q3);
464 /*
465 there is probably a better solution.
466 we have to solve for H in
467 [x1 x2]
468 [ h1 -h2 h3] [y1 y2] [x1' x2']
469 [ h2 h1 h4] [ 1 1] = [y1' y2']
470
471 which is equivalent to
472 [x1 -y1 1 0] [h1] [x1']
473 [x2 -y2 1 0] [h2] [x2']
474 [y1 x1 0 1] [h3] = [y1']
475 [y2 x2 0 1] [h4] [y2']
476 The 4x4 matrix should be easily invertible
477
478 with(linalg);
479 M := Matrix([[x1, -y1, 1, 0], [x2, -y2, 1, 0], [y1, x1, 0, 1], [y2, x2, 0, 1]]);
480 inverse(M);
481 */
482 /*
483 double det = p1.x*p1.x - 2*p2.x*p1.x + p2.x*p2.x +p1.y*p1.y -2*p1.y*p2.y +p2.y*p2.y;
484 if (det == 0.) {
485 return false;
486 }
487 double h1 = (p1.x-p2.x)*(q1.x-q2.x) + (p1.y-p2.y)*(q1.y-q2.y);
488 double h2 = (p1.x-p2.x)*(q1.y-q2.y) - (p1.y-p2.y)*(q1.x-q2.x);
489 double h3 =
490 todo...
491 */
492 }
493
setTranslationFromOnePointMatrix3x3494 bool setTranslationFromOnePoint(const OFX::Point3D &p1,
495 const OFX::Point3D &q1)
496 {
497 m[0] = 1.;
498 m[1] = 0.;
499 m[2] = q1.x - p1.x;
500 m[3] = 0.;
501 m[4] = 1.;
502 m[5] = q1.y - p1.y;
503 m[6] = 0.;
504 m[7] = 0.;
505 m[8] = 1.;
506
507 return true;
508 }
509 };
510
511 inline Matrix3x3 ofxsMatRotation(double rads);
512 inline Matrix3x3 ofxsMatRotationAroundPoint(double rads, double pointX, double pointY);
513 inline Matrix3x3 ofxsMatTranslation(double translateX, double translateY);
514 inline Matrix3x3 ofxsMatScale(double scaleX, double scaleY);
515 inline Matrix3x3 ofxsMatScale(double scale);
516 inline Matrix3x3 ofxsMatScaleAroundPoint(double scaleX, double scaleY, double pointX, double pointY);
517 inline Matrix3x3 ofxsMatSkewXY(double skewX, double skewY, bool skewOrderYX);
518
519 // matrix transform from destination to source, in canonical coordinates
520 inline Matrix3x3 ofxsMatInverseTransformCanonical(double translateX, double translateY, double scaleX, double scaleY, double skewX, double skewY, bool skewOrderYX, double rads, double centerX, double centerY);
521
522 // matrix transform from source to destination in canonical coordinates
523 inline Matrix3x3 ofxsMatTransformCanonical(double translateX, double translateY, double scaleX, double scaleY, double skewX, double skewY, bool skewOrderYX, double rads, double centerX, double centerY);
524
525 /// transform from pixel coordinates to canonical coordinates
526 inline Matrix3x3 ofxsMatPixelToCanonical(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
527 double renderscaleX, //!< 0.5 for a half-resolution image
528 double renderscaleY,
529 bool fielded); //!< true if the image property kOfxImagePropField is kOfxImageFieldLower or kOfxImageFieldUpper (apply 0.5 field scale in Y
530 /// transform from canonical coordinates to pixel coordinates
531 inline Matrix3x3 ofxsMatCanonicalToPixel(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
532 double renderscaleX, //!< 0.5 for a half-resolution image
533 double renderscaleY,
534 bool fielded); //!< true if the image property kOfxImagePropField is kOfxImageFieldLower or kOfxImageFieldUpper (apply 0.5 field scale in Y
535
536 // matrix transform from destination to source, in pixel coordinates
537 inline Matrix3x3 ofxsMatInverseTransformPixel(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
538 double renderscaleX, //!< 0.5 for a half-resolution image
539 double renderscaleY,
540 bool fielded,
541 double translateX, double translateY,
542 double scaleX, double scaleY,
543 double skewX,
544 double skewY,
545 bool skewOrderYX,
546 double rads,
547 double centerX, double centerY);
548
549 // matrix transform from source to destination in pixel coordinates
550 inline Matrix3x3 ofxsMatTransformPixel(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
551 double renderscaleX, //!< 0.5 for a half-resolution image
552 double renderscaleY,
553 bool fielded,
554 double translateX, double translateY,
555 double scaleX, double scaleY,
556 double skewX,
557 double skewY,
558 bool skewOrderYX,
559 double rads,
560 double centerX, double centerY);
561 struct Matrix4x4
562 {
563 double m[16];
564
Matrix4x4Matrix4x4565 Matrix4x4()
566 {
567 std::fill(m, m + 16, 0.);
568 }
569
Matrix4x4Matrix4x4570 Matrix4x4(const double d[16])
571 {
572 std::copy(d, d + 16, m);
573 }
574
Matrix4x4Matrix4x4575 Matrix4x4(const Matrix4x4 & o)
576 {
577 std::copy(o.m, o.m + 16, m);
578 }
579
operatorMatrix4x4580 double & operator()(int row,
581 int col)
582 {
583 assert(row >= 0 && row < 4 && col >= 0 && col < 4);
584
585 return m[row * 4 + col];
586 }
587
operatorMatrix4x4588 double operator()(int row,
589 int col) const
590 {
591 assert(row >= 0 && row < 4 && col >= 0 && col < 4);
592
593 return m[row * 4 + col];
594 }
595
determinantMatrix4x4596 double determinant() const
597 {
598 double inv0 = (m[5] * m[10] * m[15] -
599 m[5] * m[11] * m[14] -
600 m[9] * m[6] * m[15] +
601 m[9] * m[7] * m[14] +
602 m[13] * m[6] * m[11] -
603 m[13] * m[7] * m[10]);
604
605 double inv4 = (-m[4] * m[10] * m[15] +
606 m[4] * m[11] * m[14] +
607 m[8] * m[6] * m[15] -
608 m[8] * m[7] * m[14] -
609 m[12] * m[6] * m[11] +
610 m[12] * m[7] * m[10]);
611
612 double inv8 = (m[4] * m[9] * m[15] -
613 m[4] * m[11] * m[13] -
614 m[8] * m[5] * m[15] +
615 m[8] * m[7] * m[13] +
616 m[12] * m[5] * m[11] -
617 m[12] * m[7] * m[9]);
618
619 double inv12 = (-m[4] * m[9] * m[14] +
620 m[4] * m[10] * m[13] +
621 m[8] * m[5] * m[14] -
622 m[8] * m[6] * m[13] -
623 m[12] * m[5] * m[10] +
624 m[12] * m[6] * m[9]);
625
626 return m[0] * inv0 + m[1] * inv4 + m[2] * inv8 + m[3] * inv12;
627 }
628
inverseMatrix4x4629 bool inverse(Matrix4x4* invOut) const
630 {
631 double inv[16], det;
632 int i;
633
634 inv[0] = (m[5] * m[10] * m[15] -
635 m[5] * m[11] * m[14] -
636 m[9] * m[6] * m[15] +
637 m[9] * m[7] * m[14] +
638 m[13] * m[6] * m[11] -
639 m[13] * m[7] * m[10]);
640
641 inv[4] = (-m[4] * m[10] * m[15] +
642 m[4] * m[11] * m[14] +
643 m[8] * m[6] * m[15] -
644 m[8] * m[7] * m[14] -
645 m[12] * m[6] * m[11] +
646 m[12] * m[7] * m[10]);
647
648 inv[8] = (m[4] * m[9] * m[15] -
649 m[4] * m[11] * m[13] -
650 m[8] * m[5] * m[15] +
651 m[8] * m[7] * m[13] +
652 m[12] * m[5] * m[11] -
653 m[12] * m[7] * m[9]);
654
655 inv[12] = (-m[4] * m[9] * m[14] +
656 m[4] * m[10] * m[13] +
657 m[8] * m[5] * m[14] -
658 m[8] * m[6] * m[13] -
659 m[12] * m[5] * m[10] +
660 m[12] * m[6] * m[9]);
661
662 det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
663
664 if (det == 0) {
665 return false;
666 }
667
668 inv[1] = (-m[1] * m[10] * m[15] +
669 m[1] * m[11] * m[14] +
670 m[9] * m[2] * m[15] -
671 m[9] * m[3] * m[14] -
672 m[13] * m[2] * m[11] +
673 m[13] * m[3] * m[10]);
674
675 inv[5] = (m[0] * m[10] * m[15] -
676 m[0] * m[11] * m[14] -
677 m[8] * m[2] * m[15] +
678 m[8] * m[3] * m[14] +
679 m[12] * m[2] * m[11] -
680 m[12] * m[3] * m[10]);
681
682 inv[9] = (-m[0] * m[9] * m[15] +
683 m[0] * m[11] * m[13] +
684 m[8] * m[1] * m[15] -
685 m[8] * m[3] * m[13] -
686 m[12] * m[1] * m[11] +
687 m[12] * m[3] * m[9]);
688
689 inv[13] = (m[0] * m[9] * m[14] -
690 m[0] * m[10] * m[13] -
691 m[8] * m[1] * m[14] +
692 m[8] * m[2] * m[13] +
693 m[12] * m[1] * m[10] -
694 m[12] * m[2] * m[9]);
695
696 inv[2] = (m[1] * m[6] * m[15] -
697 m[1] * m[7] * m[14] -
698 m[5] * m[2] * m[15] +
699 m[5] * m[3] * m[14] +
700 m[13] * m[2] * m[7] -
701 m[13] * m[3] * m[6]);
702
703 inv[6] = (-m[0] * m[6] * m[15] +
704 m[0] * m[7] * m[14] +
705 m[4] * m[2] * m[15] -
706 m[4] * m[3] * m[14] -
707 m[12] * m[2] * m[7] +
708 m[12] * m[3] * m[6]);
709
710 inv[10] = (m[0] * m[5] * m[15] -
711 m[0] * m[7] * m[13] -
712 m[4] * m[1] * m[15] +
713 m[4] * m[3] * m[13] +
714 m[12] * m[1] * m[7] -
715 m[12] * m[3] * m[5]);
716
717 inv[14] = (-m[0] * m[5] * m[14] +
718 m[0] * m[6] * m[13] +
719 m[4] * m[1] * m[14] -
720 m[4] * m[2] * m[13] -
721 m[12] * m[1] * m[6] +
722 m[12] * m[2] * m[5]);
723
724 inv[3] = (-m[1] * m[6] * m[11] +
725 m[1] * m[7] * m[10] +
726 m[5] * m[2] * m[11] -
727 m[5] * m[3] * m[10] -
728 m[9] * m[2] * m[7] +
729 m[9] * m[3] * m[6]);
730
731 inv[7] = (m[0] * m[6] * m[11] -
732 m[0] * m[7] * m[10] -
733 m[4] * m[2] * m[11] +
734 m[4] * m[3] * m[10] +
735 m[8] * m[2] * m[7] -
736 m[8] * m[3] * m[6]);
737
738 inv[11] = (-m[0] * m[5] * m[11] +
739 m[0] * m[7] * m[9] +
740 m[4] * m[1] * m[11] -
741 m[4] * m[3] * m[9] -
742 m[8] * m[1] * m[7] +
743 m[8] * m[3] * m[5]);
744
745 inv[15] = (m[0] * m[5] * m[10] -
746 m[0] * m[6] * m[9] -
747 m[4] * m[1] * m[10] +
748 m[4] * m[2] * m[9] +
749 m[8] * m[1] * m[6] -
750 m[8] * m[2] * m[5]);
751
752 det = 1.0 / det;
753
754 for (i = 0; i < 16; i++) {
755 invOut->m[i] = inv[i] * det;
756 }
757
758 return true;
759 }
760 };
761
762 inline Matrix4x4
763 operator*(const Matrix4x4 & m1,
764 const Matrix4x4 & m2)
765 {
766 Matrix4x4 ret;
767
768 for (int i = 0; i < 4; ++i) {
769 for (int j = 0; j < 4; ++j) {
770 for (int x = 0; x < 4; ++x) {
771 ret(i, j) += m1(i, x) * m2(x, j);
772 }
773 }
774 }
775
776 return ret;
777 }
778
779 inline Point4D
780 operator*(const Matrix4x4 & m,
781 const Point4D & p)
782 {
783 Point4D ret;
784
785 for (int i = 0; i < 4; ++i) {
786 ret(i) = 0.;
787 for (int j = 0; j < 4; ++j) {
788 ret(i) += m(i, j) * p(j);
789 }
790 }
791
792 return ret;
793 }
794
795 inline Matrix4x4
matrix4x4FromMatrix3x3(const Matrix3x3 & m)796 matrix4x4FromMatrix3x3(const Matrix3x3 & m)
797 {
798 Matrix4x4 ret;
799
800 for (int i = 0; i < 3; ++i) {
801 for (int j = 0; j < 3; ++j) {
802 ret(i,j) = m(i,j);
803 }
804 }
805 ret(0, 3) = 0.;
806 ret(1, 3) = 0.;
807 ret(2, 3) = 0.;
808 ret(3, 0) = 0.; ret(3, 1) = 0.; ret(3, 2) = 0.; ret(3, 3) = 1.;
809
810 return ret;
811 }
812
813 ////////////////////
814 // IMPLEMENTATION //
815 ////////////////////
816
817
818 inline Matrix3x3
ofxsMatRotation(double rads)819 ofxsMatRotation(double rads)
820 {
821 double c = std::cos(rads);
822 double s = std::sin(rads);
823
824 return Matrix3x3(c, s, 0, -s, c, 0, 0, 0, 1);
825 }
826
827 inline Matrix3x3
ofxsMatRotationAroundPoint(double rads,double px,double py)828 ofxsMatRotationAroundPoint(double rads,
829 double px,
830 double py)
831 {
832 return ofxsMatTranslation(px, py) * ( ofxsMatRotation(rads) * ofxsMatTranslation(-px, -py) );
833 }
834
835 inline Matrix3x3
ofxsMatTranslation(double x,double y)836 ofxsMatTranslation(double x,
837 double y)
838 {
839 return Matrix3x3(1., 0., x,
840 0., 1., y,
841 0., 0., 1.);
842 }
843
844 inline Matrix3x3
ofxsMatScale(double x,double y)845 ofxsMatScale(double x,
846 double y)
847 {
848 return Matrix3x3(x, 0., 0.,
849 0., y, 0.,
850 0., 0., 1.);
851 }
852
853 inline Matrix3x3
ofxsMatScale(double s)854 ofxsMatScale(double s)
855 {
856 return ofxsMatScale(s, s);
857 }
858
859 inline Matrix3x3
ofxsMatScaleAroundPoint(double scaleX,double scaleY,double px,double py)860 ofxsMatScaleAroundPoint(double scaleX,
861 double scaleY,
862 double px,
863 double py)
864 {
865 return ofxsMatTranslation(px, py) * ( ofxsMatScale(scaleX, scaleY) * ofxsMatTranslation(-px, -py) );
866 }
867
868 inline Matrix3x3
ofxsMatSkewXY(double skewX,double skewY,bool skewOrderYX)869 ofxsMatSkewXY(double skewX,
870 double skewY,
871 bool skewOrderYX)
872 {
873 return Matrix3x3(skewOrderYX ? 1. : (1. + skewX * skewY), skewX, 0.,
874 skewY, skewOrderYX ? (1. + skewX * skewY) : 1, 0.,
875 0., 0., 1.);
876 }
877
878 // matrix transform from destination to source
879 inline Matrix3x3
ofxsMatInverseTransformCanonical(double translateX,double translateY,double scaleX,double scaleY,double skewX,double skewY,bool skewOrderYX,double rads,double centerX,double centerY)880 ofxsMatInverseTransformCanonical(double translateX,
881 double translateY,
882 double scaleX,
883 double scaleY,
884 double skewX,
885 double skewY,
886 bool skewOrderYX,
887 double rads,
888 double centerX,
889 double centerY)
890 {
891 ///1) We translate to the center of the transform.
892 ///2) We scale
893 ///3) We apply skewX and skewY in the right order
894 ///4) We rotate
895 ///5) We apply the global translation
896 ///5) We translate back to the origin
897
898 // avoid division by zero (use FLT_MIN instead of DBL_MIN in case there are further casts to float)
899 if (std::fabs(scaleX) < FLT_MIN ) {
900 scaleX = scaleX > 0 ? FLT_MIN : -FLT_MIN;
901 }
902 if (std::fabs(scaleY) < FLT_MIN ) {
903 scaleY = scaleY > 0 ? FLT_MIN : -FLT_MIN;
904 }
905 // since this is the inverse, oerations are in reverse order
906 return ofxsMatTranslation(centerX, centerY) *
907 ofxsMatScale(1. / scaleX, 1. / scaleY) *
908 ofxsMatSkewXY(-skewX, -skewY, !skewOrderYX) *
909 ofxsMatRotation(rads) *
910 ofxsMatTranslation(-translateX, -translateY) *
911 ofxsMatTranslation(-centerX, -centerY);
912 }
913
914 // matrix transform from source to destination
915 inline Matrix3x3
ofxsMatTransformCanonical(double translateX,double translateY,double scaleX,double scaleY,double skewX,double skewY,bool skewOrderYX,double rads,double centerX,double centerY)916 ofxsMatTransformCanonical(double translateX,
917 double translateY,
918 double scaleX,
919 double scaleY,
920 double skewX,
921 double skewY,
922 bool skewOrderYX,
923 double rads,
924 double centerX,
925 double centerY)
926 {
927 ///1) We translate to the center of the transform.
928 ///2) We scale
929 ///3) We apply skewX and skewY in the right order
930 ///4) We rotate
931 ///5) We apply the global translation
932 ///5) We translate back to the origin
933
934 return ofxsMatTranslation(centerX, centerY) *
935 ofxsMatTranslation(translateX, translateY) *
936 ofxsMatRotation(-rads) *
937 ofxsMatSkewXY(skewX, skewY, skewOrderYX) *
938 ofxsMatScale(scaleX, scaleY) *
939 ofxsMatTranslation(-centerX, -centerY);
940 }
941
942 // The transforms between pixel and canonical coordinated
943 // http://openfx.sourceforge.net/Documentation/1.3/ofxProgrammingReference.html#MappingCoordinates
944
945 /// transform from pixel coordinates to canonical coordinates
946 inline Matrix3x3
ofxsMatPixelToCanonical(double pixelaspectratio,double renderscaleX,double renderscaleY,bool fielded)947 ofxsMatPixelToCanonical(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
948 double renderscaleX, //!< 0.5 for a half-resolution image
949 double renderscaleY,
950 bool fielded) //!< true if the image property kOfxImagePropField is kOfxImageFieldLower or kOfxImageFieldUpper (apply 0.5 field scale in Y
951 {
952 /*
953 To map an X and Y coordinates from Pixel coordinates to Canonical coordinates, we perform the following multiplications...
954
955 X' = (X * PAR)/SX
956 Y' = Y/(SY * FS)
957 */
958
959 // FIXME: when it's the Upper field, showuldn't the first pixel start at canonical coordinate (0,0.5) ?
960 return ofxsMatScale( pixelaspectratio / renderscaleX, 1. / ( renderscaleY * (fielded ? 0.5 : 1.0) ) );
961 }
962
963 /// transform from canonical coordinates to pixel coordinates
964 inline Matrix3x3
ofxsMatCanonicalToPixel(double pixelaspectratio,double renderscaleX,double renderscaleY,bool fielded)965 ofxsMatCanonicalToPixel(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
966 double renderscaleX, //!< 0.5 for a half-resolution image
967 double renderscaleY,
968 bool fielded) //!< true if the image property kOfxImagePropField is kOfxImageFieldLower or kOfxImageFieldUpper (apply 0.5 field scale in Y
969 {
970 /*
971 To map an X and Y coordinates from Canonical coordinates to Pixel coordinates, we perform the following multiplications...
972
973 X' = (X * SX)/PAR
974 Y' = Y * SY * FS
975 */
976
977 // FIXME: when it's the Upper field, showuldn't the first pixel start at canonical coordinate (0,0.5) ?
978 return ofxsMatScale( renderscaleX / pixelaspectratio, renderscaleY * (fielded ? 0.5 : 1.0) );
979 }
980
981 // matrix transform from destination to source
982 inline Matrix3x3
ofxsMatInverseTransformPixel(double pixelaspectratio,double renderscaleX,double renderscaleY,bool fielded,double translateX,double translateY,double scaleX,double scaleY,double skewX,double skewY,bool skewOrderYX,double rads,double centerX,double centerY)983 ofxsMatInverseTransformPixel(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
984 double renderscaleX, //!< 0.5 for a half-resolution image
985 double renderscaleY,
986 bool fielded,
987 double translateX,
988 double translateY,
989 double scaleX,
990 double scaleY,
991 double skewX,
992 double skewY,
993 bool skewOrderYX,
994 double rads,
995 double centerX,
996 double centerY)
997 {
998 ///1) We go from pixel to canonical
999 ///2) we apply transform
1000 ///3) We go back to pixels
1001
1002 return ofxsMatCanonicalToPixel(pixelaspectratio, renderscaleX, renderscaleY, fielded) *
1003 ofxsMatInverseTransformCanonical(translateX, translateY, scaleX, scaleY, skewX, skewY, skewOrderYX, rads, centerX, centerY) *
1004 ofxsMatPixelToCanonical(pixelaspectratio, renderscaleX, renderscaleY, fielded);
1005 }
1006
1007 // matrix transform from source to destination
1008 inline Matrix3x3
ofxsMatTransformPixel(double pixelaspectratio,double renderscaleX,double renderscaleY,bool fielded,double translateX,double translateY,double scaleX,double scaleY,double skewX,double skewY,bool skewOrderYX,double rads,double centerX,double centerY)1009 ofxsMatTransformPixel(double pixelaspectratio, //!< 1.067 for PAL, where 720x576 pixels occupy 768x576 in canonical coords
1010 double renderscaleX, //!< 0.5 for a half-resolution image
1011 double renderscaleY,
1012 bool fielded,
1013 double translateX,
1014 double translateY,
1015 double scaleX,
1016 double scaleY,
1017 double skewX,
1018 double skewY,
1019 bool skewOrderYX,
1020 double rads,
1021 double centerX,
1022 double centerY)
1023 {
1024 ///1) We go from pixel to canonical
1025 ///2) we apply transform
1026 ///3) We go back to pixels
1027
1028 return ofxsMatCanonicalToPixel(pixelaspectratio, renderscaleX, renderscaleY, fielded) *
1029 ofxsMatTransformCanonical(translateX, translateY, scaleX, scaleY, skewX, skewY, skewOrderYX, rads, centerX, centerY) *
1030 ofxsMatPixelToCanonical(pixelaspectratio, renderscaleX, renderscaleY, fielded);
1031 }
1032 };
1033
1034 #endif // openfx_supportext_ofxsMatrix2D_h
1035