Implement the 3D image engine from scratch: (7) matrix function library

1. Mathematical Analysis

1) What is a matrix and what is it used?

Do not think about the matrix as complicated. To put it bluntly, he and the abacus are a class of things, and the matrix has no practical significance at all. He is just a mathematical tool, just like an abacus, the same result can be obtained when the beads are played up and down according to certain rules. Matrix is a tool invented by mathematicians. This tool also has a series of rules, these rules can also calculate some results, but unlike the abacus, it is used for addition, subtraction, multiplication, and division of numbers. It is used to solve equations. After learning the calculation rules of the matrix, we can find that all the transformations of vectors and coordinate points can be completed through matrix operations, which is why 3D image operations are closely related to the matrix.

 

For example, equations:

X + 2y = 3

4x + 5y = 6

We can use three matrices to represent them.

 

Coefficient Matrix: A =

| 1 2 |

| 4 5 |

 

Variable matrix: B =

| X |

| Y |

 

Constant Matrix: c =

| 3 |

| 6 |

 

So the above equations can also be written as: a * B = C

Therefore, the problem of solving the equations also becomes B. Next I will introduce the inverse of the matrix, and I will talk about how to use matrix A and matrix C to calculate matrix B.

 

 

2) Unit Matrix

Definition: all elements on the main diagonal line of the matrix are 1, and the rest are 0. It is represented by I.

For example:

| 1 0 |

| 0 1 |

Or:

| 1 0 0 |

| 0 1 0 |

| 0 0 1 |

The important purpose of the unit matrix is that it satisfies the following formula:

A * I = I * A =

 

 

3) zero matrix

Very simple:

| 0 0 |

| 0 0 |

All elements are zero.

 

 

4) matrix addition and subtraction

Add or subtract the corresponding elements of the two matrices.

 

 

5) transpose the Matrix

Transpose the rows and column elements of a matrix, for example, a =

| 0 1 2 |

| 3 4 5 |

Transpose to a' =

| 0 3 |

| 1 4 |

| 2 5 |

 

 

6) Multiply the matrix and scalar

Each element is multiplied by a scalar.

 

 

7) Multiply the matrix and Matrix

Not any two matrices can be multiplied. matrix A and matrix B are multiplied. It makes sense only when their inner dimensions are equal.

For example, if A is an m x nj matrix, B must be an n X R matrix.

The algorithm is as follows: a * B = C. Then, the vector of line I in CIJ = A and the dot product of the vector of column J in B.

For example, a =

| 0 1 |

| 2 3 |

 

B =

| 0 1 2 |

| 3 4 5 |

 

Then c =

| (). () |

| (). () | =

 

| (0*0 + 1*3) (0*1 + 1*4) (0*2 + 1*5) |

| (2*0 + 3*3) (2*1 + 3*4) (2*2 + 3*5) | =

 

| 3 4 5 |

| 9 14 19 |

 

8) matrix calculation Law

1. A + B = B +

2. A + (B + C) = (a + B) + C

3. A * (B * C) = (a * B) * C

4. A * (B + C) = a * B + A * C

5. K * (a + B) = K * A + K * B

6. (A + B) * C = A * C + B * C

7. A * I = I * A =

 

Note: A matrix generally does not meet the multiplication exchange law, that is, a * B is not equal to B *.

 

 

9) the determining factor of the matrix

The source of the determinant: each square matrix can correspond to a real number that becomes the determinant of the matrix. This value tells us whether the matrix is singular.

The determinant of matrix A is det ()

In practical application, we need to calculate the inverse of the matrix by finding the deciding factor, so the significance of the determining factor is very significant.

How can we use the above source to derive the formula for determining the determinant? In any linear algebraic number, I will omit it. The formula is as follows:

For N x n matrix,

Det (A) = A11 (when n = 1)

Det (A) = A11 * A11 + A12 * A12 +... + A1N * A1N (when n> 1)

Among them, a1j = (-1) 1 + J power * DETA (m1j) (where j = 1,..., n)

It is the remainder of the first line element.

 

Let's talk about the details, 2x2 matrix.

| A B |

| C d |

Det (A) = A * D-B * C

 

3x3 matrix

| A00 A01 A02 |

| A10 A11 A12 |

| A20 A21 A22 |

Is det () = a00 * A11 * A22 + A01 * A12 * A20 + A02 * A10 * A21-A02 * A11 * A20-A01 * A10 * A22-a00 * A12 * A21

This formula is obtained by using the Klein law. Why is the remainder formula described above, because 9 multiplication and 5 addition are required for the determination of the 3 × 3 matrix using the Klein law. However, if we use the remainder formula to calculate the determinant, we can reduce the multiplication by two times:

Det (A) = a00 * (A11 * A22-A21 * A12)-A01 * (A10 * A22-A20 * A12)-A02 * (A10 * A21-A20 * A11)

 

 

10) inverse of the matrix

Definition: For matrix A, there is a matrix A-1 that makes a * A-1 = I, then the A-1 is the inverse of matrix. In fact, this A-1 is a multiplication inverse element.

There is also a saying in mathematics: If a n × n matrix does not have a multiplication inverse element, it is called a strange matrix.

His role is very important, such as the equations mentioned above:

A * B = C, matrix B

By multiplying the two sides of the equation by the A-1, then:

(A-1 * A) * B = A-1 * C

B = A-1 * C, so we can obtain B.

 

 

2. Code Implementation

1) I provide a matrix of 1x3, 1X4, 3x3, and 4x4, as for 2D points or vectors, homogeneous coordinates can be stored in the 3x3 matrix. Likewise, for 3D points or vectors, homogeneous coordinates can be stored in the 4x4 matrix, therefore, you only need to write the calculation functions of the above four matrices to meet all the 3D matrix transformations. This time, I have modified the structure of previously written vectors. Because vectors and vertices can usually be expressed as a matrix of 1x3 or 1X4, so I changed their memory structure to the same as the matrix, so that the conversion can be forced directly.

The following is the definition of a matrix structure:

Typedef struct matrix1x3_type // 1 × 3 matrix <br/>{< br/> Union <br/>{< br/> double M [3]; <br/> struct <br/>{< br/> double m00, M01, M02; <br/>}; <br/>} matrix1x3, * matrix1x3_ptr; </P> <p> typedef struct matrix1x4_type // 1X4 matrix <br/> {<br/> Union <br/>{< br/> double M [4]; <br/> struct <br/>{< br/> double m00, M01, M02, M03; <br/> }; <br/>} matrix1x4, * matrix1x4_ptr; </P> <p> typedef struct matrix3_type // 3 × 3 matrix <br/> {<br/> Union <br/> {<br/> double M [3] [3]; <br/> struct <br/>{< br/> double m00, M01, M02; <br/> double M10, M11, M12; <br/> double M20, m21, m22; <br/>}; <br/>} matrix3x3, * matrix3x3_ptr; </P> <p> typedef struct matrix4x4_type // 4x4 matrix <br/> {<br/> Union <br/>{< br/> double M [4] [4]; <br/> struct <br/> {<br/> double m00, M01, M02, M03; <br/> double M10, M11, M12, M13; <br/> double M20, m21, m22, M23; <br/> double m30, m31, M32, M33; <br/>}; <br/> }; <br/>} matrix4x4, * matrix4x4_ptr;

 

The following is a change to the structure of the original vector and vertex:

 Typedef struct point2d_type // 2D Cartesian coordinate, 2D vector <br/>{< br/> Union <br/>{< br/> double M [2]; // array mode <br/> struct // variable mode <br/>{< br/> double X; <br/> double y; <br/> }; <br/>}; <br/>} point2d, * point2d_ptr, vector2d, * vector2d_ptr; </P> <p> typedef struct point3d_type/3D Cartesian coordinate, 3D vectors <br/>{< br/> Union <br/>{< br/> double M [3]; <br/> struct <br/>{< br/> double X; <br/> double y; <br/> double Z; <br/> }; <br/>}; <br/>} point3d, * point3d_ptr, vector3d, * vector3d_ptr; </P> <p> typedef struct vector4d_type // 4D vector <br/>{< br/> Union <br/>{< br/> double M [4]; <br/> struct <br/>{< br/> double X; <br/> double y; <br/> double Z; <br/> double W; <br/>}; <br/>} vector4d, * vector4d_ptr;

 

2) matrix operation function implementation

Void _ cppyin_math: matrixcreate (matrix3x3_ptr pm, <br/> double m00, double M01, double M02, <br/> double M10, double M11, double M12, <br/> double M20, double m21, double m22) <br/>{< br/> PM-> m00 = m00; <br/> PM-> M01 = M01; <br/> PM-> M02 = M02; <br/> PM-> M10 = M10; <br/> PM-> M11 = M11; <br/> PM-> M12 = m12; <br/> PM-> M20 = M20; <br/> PM-> m21 = m21; <br/> PM-> m22 = m22; <br/>}</P> <p> void _ cppyin_math: matrixcreate (matrix4x4_ptr pm, <br/> double m00, double M01, double M02, double M03, <br/> double M10, double M11, double M12, double M13, <br/> double M20, double m21, double m22, double M23, <br/> double m30, double m31, double M32, double M33) <br/>{< br/> PM-> m00 = m00; <br/> PM-> M01 = M01; <br/> PM-> M02 = M02; <br/> PM-> M03 = M03; <br/> PM-> M10 = M10; <br/> PM-> M11 = M11; <br/> PM-> M12 = m12; <br/> PM-> M13 = M13; <br/> PM-> M20 = M20; <br/> PM-> m21 = m21; <br/> PM-> m22 = m22; <br/> PM-> M23 = M23; <br/> PM-> m30 = m30; <br/> PM-> m31 = m31; <br/> PM-> M32 = M32; <br/> PM-> M33 = M33; <br/>}</P> <p> void _ cppyin_math: matrixadd (matrix3x3_ptr Ma, matrix3x3_ptr MB, matrix3x3_ptr msum) <br/> {<br/> for (int row = 0; row <3; ++ row) <br/> {<br/> for (INT Col = 0; col <3; ++ col) <br/> {<br/> msum-> M [row] [col] = ma-> M [row] [col] + MB-> M [row] [col ]; <br/>}</P> <p> void _ cppyin_math: matrixadd (matrix4x4_ptr Ma, matrix4x4_ptr MB, matrix4x4_ptr msum) <br/> {<br/> for (int row = 0; row <4; ++ row) <br/> {<br/> for (INT Col = 0; col <4; ++ col) <br/> {<br/> msum-> M [row] [col] = ma-> M [row] [col] + MB-> M [row] [col ]; <br/>}</P> <p> void _ cppyin_math: matrixmul (matrix1x3_ptr Ma, matrix3x3_ptr MB, matrix1x3_ptr mprod) <br/>{< br/> for (INT Col = 0; Col <3; ++ col) <br/>{< br/> double sum = 0; <br/> for (INT Index = 0; index <3; ++ index) <br/>{< br/> sum + = (MA-> M [Index] * MB-> M [Index] [col]); <br/>}< br/> mprod-> M [col] = sum; <br/>}</P> <p> void _ cppyin_math:: matrixmul (matrix1x4_ptr Ma, matrix4x4_ptr MB, matrix1x4_ptr mprod) <br/>{< br/> for (INT Col = 0; Col <4; ++ col) <br/> {<br/> double sum = 0; <br/> for (INT Index = 0; index <4; ++ index) <br/>{< br/> sum + = (MA-> M [Index] * MB-> M [Index] [col]); <br/>}< br/> mprod-> M [col] = sum; <br/>}</P> <p> void _ cppyin_math:: matrixmul (matrix3x3_ptr Ma, matrix3x3_ptr MB, matrix3x3_ptr mprod) <br/>{< br/> for (int row = 0; row <3; ++ row) <br/>{< br/> for (INT Col = 0; Col <3; ++ col) <br/>{< br/> double sum = 0.0; <br/> for (INT Index = 0; index <3; ++ index) <br/> {<br/> sum + = ma-> M [row] [Index] * MB-> M [Index] [col]; <br/>}< br/> mprod-> M [row] [col] = sum; <br/>}</P> <p> void _ cppyin_math: matrixmul (matrix4_ptr Ma, matrix4x4_ptr MB, matrix4x4_ptr mprod) <br/> {<br/> for (int row = 0; row <4; ++ row) <br/> {<br/> for (INT Col = 0; col <4; ++ col) <br/>{< br/> double sum = 0.0; <br/> for (INT Index = 0; index <4; ++ index) <br/>{< br/> sum + = ma-> M [row] [Index] * MB-> M [Index] [col]; <br/>}< br/> mprod-> M [row] [col] = sum; <br/>}</P> <p> double _ cppyin_math: matrixdet (matrix3x3_ptr m) <br/> {<br/> return (<br/> M-> m00 * (m-> M11 * m-> m22-m-> m21 * m-> M12) -M-> M01 * (m-> M10 * m-> m22-m-> M20 * m-> M12) + M-> M02 * (m-> M10 * m-> m21-m-> M20 * m-> m11) <br/> ); <br/>}</P> <p> int _ cppyin_math: matrixinverse (matrix3x3_ptr M, matrix3x3_ptr mi) <br/> {<br/> double det = m-> m00 * (m-> M11 * m-> m22-m-> m21 * m-> M12) -<br/> M-> M01 * (m-> M10 * m-> m22-m-> M20 * m-> M12) + <br/> M-> M02 * (m-> M10 * m-> m21-m-> M20 * m-> m11 ); </P> <p> If (ABS (DET) <epsilon) <br/> return 0; </P> <p> double det_inv = 1.0/det; </P> <p> mi-> m00 = det_inv * (m-> M11 * m-> m22-m-> m21 * m-> M12 ); <br/> mi-> M10 =-det_inv * (m-> M10 * m-> m22-m-> M20 * m-> M12 ); <br/> mi-> M20 = det_inv * (m-> M10 * m-> m21-m-> M20 * m-> m11 ); </P> <p> mi-> M01 =-det_inv * (m-> M01 * m-> m22-m-> m21 * m-> M02 ); <br/> mi-> M11 = det_inv * (m-> m00 * m-> m22-m-> M20 * m-> M02 ); <br/> mi-> m21 =-det_inv * (m-> m00 * m-> m21-m-> M20 * m-> M01 ); </P> <p> mi-> M02 = det_inv * (m-> M01 * m-> M12-m-> M11 * m-> M02 ); <br/> mi-> M12 =-det_inv * (m-> m00 * m-> M12-m-> M10 * m-> M02 ); <br/> mi-> m22 = det_inv * (m-> m00 * m-> M11-m-> M10 * m-> M01 ); </P> <p> return 1; <br/>}</P> <p> int _ cppyin_math: matrixinverse (matrix4x4_ptr M, matrix4x4_ptr mi) <br/>{< br/> double det = (m-> m00 * (m-> M11 * m-> m22-m-> M12 * m-> m21) -<br/> M-> M01 * (m-> M10 * m-> m22-m-> M12 * m-> M20) + <br/> M-> M02 * (m-> M10 * m-> m21-m-> M11 * m-> M20 )); </P> <p> If (ABS (DET) <epsilon) <br/> return 0; </P> <p> double det_inv = 1.0/det; </P> <p> mi-> m00 = det_inv * (m-> M11 * m-> m22-m-> M12 * m-> m21 ); <br/> mi-> M01 =-det_inv * (m-> M01 * m-> m22-m-> M02 * m-> m21 ); <br/> mi-> M02 = det_inv * (m-> M01 * m-> M12-m-> M02 * m-> m11 ); <br/> mi-> M03 = 0.0; </P> <p> mi-> M10 =-det_inv * (m-> M10 * m-> m22-m-> M12 * m-> M20 ); <br/> mi-> M11 = det_inv * (m-> m00 * m-> m22-m-> M02 * m-> M20 ); <br/> mi-> M12 =-det_inv * (m-> m00 * m-> M12-m-> M02 * m-> M10 ); <br/> mi-> M13 = 0.0; </P> <p> mi-> M20 = det_inv * (m-> M10 * m-> m21-m-> M11 * m-> M20 ); <br/> mi-> m21 =-det_inv * (m-> m00 * m-> m21-m-> M01 * m-> M20 ); <br/> mi-> m22 = det_inv * (m-> m00 * m-> M11-m-> M01 * m-> M10 ); <br/> mi-> M23 = 0.0; </P> <p> mi-> m30 =-(m-> m30 * mi-> m00 + M-> m31 * mi-> M10 + M-> M32 * Mi -> M20 ); <br/> mi-> m31 =-(m-> m30 * mi-> M01 + M-> m31 * mi-> M11 + M-> M32 * mi-> m21 ); <br/> mi-> M32 =-(m-> m30 * mi-> M02 + M-> m31 * mi-> m12 + M-> M32 * mi-> m22 ); <br/> mi-> M33 = 1.0; </P> <p> return 1; <br/>}

 

 

 

3. Download Code

Download complete project code:> click to enter the download page <