矩阵求逆函数实现,继续扫盲

Devil

如果转载请注明出处
http://www.azure.com.cn/

矩阵表示

[ A11  A12  A13  A14 ]
[ A21  A22  A23  A24 ]
[ A31  A32  A33  A34 ]  /  det(A)
[ A41  A42  A43  A44 ]


函数实现（假设我们已经实现了Matrix4类)

1. bool MatrixInvert(Matrix4& inv, Matrix4& source)
   
2. {
   
3.      Matrix4 TempMat;
   
4.      
   
5.      //检验原矩阵的有效性
   
6.      if(fabs(a._44 - 1.0f) > 0.001f)
   
7.          return false;
   
8.      if(fabs(a._14)>0.001f || fabs(a._24)>0.001f || fabs(a._34)>0.001f)
   
9.          return false;
   
10. 
    
11.      //求矩阵行列式的值的倒数
    
12.      float fDetInv = 1.0f / (source._11 * (source._22 * source._33 - source._23 * source._32) -
    
13.                                     source._12 * (source._21 * source._33 - source._23 * source._31) +
    
14.                                     source._13 * (source._21 * source._32 - source._22 * source._31));
    
15. 
    
16.      //将愿矩阵与矩阵行列式的值的倒数相乘
    
17.      TempMat._11 = fDetInv * (source._22 * source._33 - source._23 * source._32);
    
18.      TempMat._12 = -fDetInv * (source._12 * source._33 - source._13 * source._32);
    
19.      TempMat._13 = fDetInv * (source._12 * source._23 - source._13 * source._22);
    
20.      TempMat._14 = 0.0f;
    
21. 
    
22.      TempMat._21 = -fDetInv * (source._21 * source._33 - source._23 * source._31);
    
23.      TempMat._22 = fDetInv * (source._11 * source._23 - source._13 * source._31);
    
24.      TempMat._23 = -fDetInv * (source._11 * source._23 - source._13 * source._21);
    
25.      TempMat._24 = 0.0f;
    
26. 
    
27.      TempMat._31 = fDetInv * (source._21 * source._32 - source._22 * source._31);
    
28.      TempMat._32 = -fDetInv * (source._11 * source._32 - source._12 * source._31);
    
29.      TempMat._33 = fDetInv * (source._11 * source._22 - source._12 * source._21);
    
30.      TempMat._34 = 0.0f;
    
31. 
    
32.      TempMat._41 = -(source._41 * TempMat._11 + source._42 * TempMat._21 + source._43 * TempMat._31);
    
33.      TempMat._42 = -(source._41 * TempMat._12 + source._42 * TempMat._22 + source._43 * TempMat._32);
    
34.      TempMat._43 = -(source._41 * TempMat._13 + source._42 * TempMat._23 + source._43 * TempMat._33);
    
35.      TempMat._44 = 1.0f;
    
36. 
    
37.      memcpy(&inv, &TempMat, sizeof(Matrix4));
    
38.      return true;
    
39. }

复制代码

bracelet

D3D下就直接用DX函数了，何必劳心费神，效率还不一定高

Devil

一句话：
知其然，知其所以然。

tarkey

给你个
1 2 3 0
4 5 6 0
7 8 9 0
1 1 1 1
你算算Inverse试试

mousesure

2楼的，基础不牢是不好嘀~

Jeminai

Devil: Re:矩阵求逆函数实现,继续扫盲

一句话：
知其然，知其所以然。

对！了解一下还是必要的。

ddrmsdos

我的这个有更多的函数........包括求逆........
//CGraphicX.h
#ifndef CGRAPHICX_H
#define CGRAPHICX_H
//定义实用结构和计算类


#define LEFT 1
#define RIGHT 2
#define BOTTOM 4
#define TOP 8

#ifndef NULL
#define NULL 0
#endif
//3*3矩阵，是2维齐次坐标系
typedef struct Matrix_2_2
{
        Matrix_2_2();
        float mat[3][3];
        const Matrix_2_2 operator*(const Matrix_2_2& opMat);
} Matrix_2_2;

typedef struct GVertex
{
        float x,y,z;
        GVertex();
        //顶点与矩阵乘法
        const GVertex operator*(const Matrix_2_2& opMat);
}GVertex;


typedef struct FRect
{
        FRect():left(0),top(0),right(0),bottom(0){}
        float left,top,right,bottom;
}FRect;

class CGraphicX
{
public:
        //单位矩阵
        static void MatrixIdentility(Matrix_2_2 *op);
        //矩阵乘法
        static void MatrixMul(Matrix_2_2 *op1,Matrix_2_2 *op2,Matrix_2_2 *result);
        //求逆
        static void MatrixOther(Matrix_2_2 *op,Matrix_2_2 *result);
        //求位移
        static void MatrixTrans(Matrix_2_2 *result,float x,float y);
        //求旋转
        static void MatrixRotate(Matrix_2_2 *result,float rad);
        //求缩放
        static void MatrixScale(Matrix_2_2 *result,float x,float y);
        //错切变换
        static void MatrixTiltedX(Matrix_2_2 *result,float c);
        static void MatrixTiltedY(Matrix_2_2 *result,float c);
        //设置视矩阵
        static void MatrixView(Matrix_2_2 *result,int width,int height,int left,int top,int right,int bottom);
        static int Clip_Line(int &x1,int &y1,int &x2, int &y2,FRect *pRect);
private:
        //求2*2矩阵行列式
        static float MatrixValue(float a1,float a2,float a3,float a4);
        //求3*3矩阵行列式
        static float MatrixValue_2_2(Matrix_2_2 *op);       
       
};


#endif



#include "CGraphicX.h"

#include <math.h>


int CGraphicX::Clip_Line(int &x1,int &y1,int &x2,int &y2,FRect *pRect)
{
        // this function clips the sent line using the globally defined clipping
        // region

        // internal clipping codes
#define CLIP_CODE_C  0x0000
#define CLIP_CODE_N  0x0008
#define CLIP_CODE_S  0x0004
#define CLIP_CODE_E  0x0002
#define CLIP_CODE_W  0x0001

#define CLIP_CODE_NE 0x000a
#define CLIP_CODE_SE 0x0006
#define CLIP_CODE_NW 0x0009
#define CLIP_CODE_SW 0x0005

        int xc1=x1,
                        yc1=y1,
                xc2=x2,
                yc2=y2;

        int p1_code=0,
                p2_code=0;

        int min_clip_y;
        int max_clip_y;
        int min_clip_x;
        int max_clip_x;
        if(pRect->top<pRect->bottom)
        {
                min_clip_y=(int)(pRect->top);
                max_clip_y=(int)(pRect->bottom);
        }
        else
        {
                max_clip_y=(int)(pRect->top);
                min_clip_y=(int)(pRect->bottom);
        }
        if(pRect->left<pRect->right)
        {
                min_clip_x=(int)(pRect->left);
                max_clip_x=(int)(pRect->right);
        }
        else
        {
                max_clip_x=(int)(pRect->left);
                min_clip_x=(int)(pRect->right);
        }
       

        // determine codes for p1 and p2
        if (y1 < min_clip_y)
                p1_code|=CLIP_CODE_N;
        else
                if (y1 > max_clip_y)
                        p1_code|=CLIP_CODE_S;

        if (x1 < min_clip_x)
                p1_code|=CLIP_CODE_W;
        else
                if (x1 > max_clip_x)
                        p1_code|=CLIP_CODE_E;

        if (y2 < min_clip_y)
                p2_code|=CLIP_CODE_N;
        else
                if (y2 > max_clip_y)
                        p2_code|=CLIP_CODE_S;

        if (x2 < min_clip_x)
                p2_code|=CLIP_CODE_W;
        else
                if (x2 > max_clip_x)
                        p2_code|=CLIP_CODE_E;

        // try and trivially reject
        if ((p1_code & p2_code))
                return(0);

        // test for totally visible, if so leave points untouched
        if (p1_code==0 && p2_code==0)
                return(1);

        // determine end clip point for p1
        switch(p1_code)
        {
        case CLIP_CODE_C: break;

        case CLIP_CODE_N:
                {
                        yc1 = min_clip_y;
                        xc1 = x1 + 0.5+(min_clip_y-y1)*(x2-x1)/(y2-y1);
                } break;
        case CLIP_CODE_S:
                {
                        yc1 = max_clip_y;
                        xc1 = x1 + 0.5+(max_clip_y-y1)*(x2-x1)/(y2-y1);
                } break;

        case CLIP_CODE_W:
                {
                        xc1 = min_clip_x;
                        yc1 = y1 + 0.5+(min_clip_x-x1)*(y2-y1)/(x2-x1);
                } break;

        case CLIP_CODE_E:
                {
                        xc1 = max_clip_x;
                        yc1 = y1 + 0.5+(max_clip_x-x1)*(y2-y1)/(x2-x1);
                } break;

                // these cases are more complex, must compute 2 intersections
        case CLIP_CODE_NE:
                {
                        // north hline intersection
                        yc1 = min_clip_y;
                        xc1 = x1 + 0.5+(min_clip_y-y1)*(x2-x1)/(y2-y1);

                        // test if intersection is valid, of so then done, else compute next
                        if (xc1 < min_clip_x || xc1 > max_clip_x)
                        {
                                // east vline intersection
                                xc1 = max_clip_x;
                                yc1 = y1 + 0.5+(max_clip_x-x1)*(y2-y1)/(x2-x1);
                        } // end if

                } break;

        case CLIP_CODE_SE:
                {
                        // south hline intersection
                        yc1 = max_clip_y;
                        xc1 = x1 + 0.5+(max_clip_y-y1)*(x2-x1)/(y2-y1);       

                        // test if intersection is valid, of so then done, else compute next
                        if (xc1 < min_clip_x || xc1 > max_clip_x)
                        {
                                // east vline intersection
                                xc1 = max_clip_x;
                                yc1 = y1 + 0.5+(max_clip_x-x1)*(y2-y1)/(x2-x1);
                        } // end if

                } break;

        case CLIP_CODE_NW:
                {
                        // north hline intersection
                        yc1 = min_clip_y;
                        xc1 = x1 + 0.5+(min_clip_y-y1)*(x2-x1)/(y2-y1);

                        // test if intersection is valid, of so then done, else compute next
                        if (xc1 < min_clip_x || xc1 > max_clip_x)
                        {
                                xc1 = min_clip_x;
                                yc1 = y1 + 0.5+(min_clip_x-x1)*(y2-y1)/(x2-x1);       
                        } // end if

                } break;

        case CLIP_CODE_SW:
                {
                        // south hline intersection
                        yc1 = max_clip_y;
                        xc1 = x1 + 0.5+(max_clip_y-y1)*(x2-x1)/(y2-y1);       

                        // test if intersection is valid, of so then done, else compute next
                        if (xc1 < min_clip_x || xc1 > max_clip_x)
                        {
                                xc1 = min_clip_x;
                                yc1 = y1 + 0.5+(min_clip_x-x1)*(y2-y1)/(x2-x1);       
                        } // end if

                } break;

        default:break;

        } // end switch

        // determine clip point for p2
        switch(p2_code)
        {
        case CLIP_CODE_C: break;

        case CLIP_CODE_N:
                {
                        yc2 = min_clip_y;
                        xc2 = x2 + (min_clip_y-y2)*(x1-x2)/(y1-y2);
                } break;

        case CLIP_CODE_S:
                {
                        yc2 = max_clip_y;
                        xc2 = x2 + (max_clip_y-y2)*(x1-x2)/(y1-y2);
                } break;

        case CLIP_CODE_W:
                {
                        xc2 = min_clip_x;
                        yc2 = y2 + (min_clip_x-x2)*(y1-y2)/(x1-x2);
                } break;

        case CLIP_CODE_E:
                {
                        xc2 = max_clip_x;
                        yc2 = y2 + (max_clip_x-x2)*(y1-y2)/(x1-x2);
                } break;

                // these cases are more complex, must compute 2 intersections
        case CLIP_CODE_NE:
                {
                        // north hline intersection
                        yc2 = min_clip_y;
                        xc2 = x2 + 0.5+(min_clip_y-y2)*(x1-x2)/(y1-y2);

                        // test if intersection is valid, of so then done, else compute next
                        if (xc2 < min_clip_x || xc2 > max_clip_x)
                        {
                                // east vline intersection
                                xc2 = max_clip_x;
                                yc2 = y2 + 0.5+(max_clip_x-x2)*(y1-y2)/(x1-x2);
                        } // end if

                } break;

        case CLIP_CODE_SE:
                {
                        // south hline intersection
                        yc2 = max_clip_y;
                        xc2 = x2 + 0.5+(max_clip_y-y2)*(x1-x2)/(y1-y2);       

                        // test if intersection is valid, of so then done, else compute next
                        if (xc2 < min_clip_x || xc2 > max_clip_x)
                        {
                                // east vline intersection
                                xc2 = max_clip_x;
                                yc2 = y2 + 0.5+(max_clip_x-x2)*(y1-y2)/(x1-x2);
                        } // end if

                } break;

        case CLIP_CODE_NW:
                {
                        // north hline intersection
                        yc2 = min_clip_y;
                        xc2 = x2 + 0.5+(min_clip_y-y2)*(x1-x2)/(y1-y2);

                        // test if intersection is valid, of so then done, else compute next
                        if (xc2 < min_clip_x || xc2 > max_clip_x)
                        {
                                xc2 = min_clip_x;
                                yc2 = y2 + 0.5+(min_clip_x-x2)*(y1-y2)/(x1-x2);       
                        } // end if

                } break;

        case CLIP_CODE_SW:
                {
                        // south hline intersection
                        yc2 = max_clip_y;
                        xc2 = x2 + 0.5+(max_clip_y-y2)*(x1-x2)/(y1-y2);       

                        // test if intersection is valid, of so then done, else compute next
                        if (xc2 < min_clip_x || xc2 > max_clip_x)
                        {
                                xc2 = min_clip_x;
                                yc2 = y2 + 0.5+(min_clip_x-x2)*(y1-y2)/(x1-x2);       
                        } // end if

                } break;

        default:break;

        } // end switch

        // do bounds check
        if ((xc1 < min_clip_x) || (xc1 > max_clip_x) ||
                (yc1 < min_clip_y) || (yc1 > max_clip_y) ||
                (xc2 < min_clip_x) || (xc2 > max_clip_x) ||
                (yc2 < min_clip_y) || (yc2 > max_clip_y) )
        {
                return(0);
        } // end if

        // store vars back
        x1 = xc1;
        y1 = yc1;
        x2 = xc2;
        y2 = yc2;

        return(1);

} // end Clip_Line


Matrix_2_2::Matrix_2_2()
{
        int i,j;
        for(i=0;i<3;++i)
                for(j=0;j<3;++j)
                        mat[j]=0.0f;
}

const Matrix_2_2 Matrix_2_2:perator *(const Matrix_2_2 &opMat)
{
        Matrix_2_2 tmpMat;
        int i,j,k;
        for(i=0;i<3;++i)
                for(j=0;j<3;++j)
                        for(k=0;k<3;++k)
                        {
                                tmpMat.mat[j]+=(this->mat[k])*(opMat.mat[k][j]);
                        }
                        return tmpMat;
}

const GVertex GVertex::operator*(const Matrix_2_2& opMat)
{
        GVertex tmpVertex;
        tmpVertex.x=x*opMat.mat[0][0]+y*opMat.mat[1][0]+z*opMat.mat[2][0];
        tmpVertex.y=x*opMat.mat[0][1]+y*opMat.mat[1][1]+z*opMat.mat[2][1];
        tmpVertex.z=x*opMat.mat[0][2]+y*opMat.mat[1][2]+z*opMat.mat[2][2];
        return tmpVertex;
}

GVertex::GVertex()
{
        x=y=z=0.0f;
}


float CGraphicX::MatrixValue(float a1, float a2, float a3, float a4)
{
        return a1*a4-a2*a3;
}

float CGraphicX::MatrixValue_2_2(Matrix_2_2 *op)
{
        return (op->mat[0][0])*(op->mat[1][1])*(op->mat[2][2])+
                (op->mat[0][1])*(op->mat[1][2])*(op->mat[2][0])+
                (op->mat[0][2])*(op->mat[1][0])*(op->mat[2][1])-
                (op->mat[0][0])*(op->mat[1][2])*(op->mat[2][1])-
                (op->mat[0][1])*(op->mat[1][0])*(op->mat[2][2])-
                (op->mat[0][2])*(op->mat[1][1])*(op->mat[2][0]);
}

void CGraphicX::MatrixOther(Matrix_2_2 *op, Matrix_2_2 *result)
{
        Matrix_2_2 a;
        float sa=MatrixValue_2_2(op);
        if(op==0)
        {
                result=0;
                return;
        }
        sa=1.0f/sa;
        a.mat[0][0]=sa*MatrixValue(op->mat[1][1],op->mat[1][2],op->mat[2][1],op->mat[2][2]);
        a.mat[1][0]=(-1.0f)*sa*MatrixValue(op->mat[1][0],op->mat[1][2],op->mat[2][0],op->mat[2][2]);
        a.mat[2][0]=sa*MatrixValue(op->mat[1][0],op->mat[1][1],op->mat[2][0],op->mat[2][1]);
        a.mat[0][1]=(-1.0f)*sa*MatrixValue(op->mat[0][1],op->mat[0][2],op->mat[2][1],op->mat[2][2]);
        a.mat[1][1]=sa*MatrixValue(op->mat[0][0],op->mat[0][2],op->mat[2][0],op->mat[2][2]);
        a.mat[2][1]=(-1.0f)*sa*MatrixValue(op->mat[0][0],op->mat[0][1],op->mat[2][0],op->mat[2][1]);
        a.mat[0][2]=sa*MatrixValue(op->mat[0][1],op->mat[0][2],op->mat[1][1],op->mat[1][2]);
        a.mat[1][2]=(-1.0f)*sa*MatrixValue(op->mat[0][0],op->mat[0][2],op->mat[1][0],op->mat[1][2]);
        a.mat[2][2]=sa*MatrixValue(op->mat[0][0],op->mat[0][1],op->mat[1][0],op->mat[1][1]);
        *result=a;
}

void CGraphicX::MatrixIdentility(Matrix_2_2 *op)
{
        for(int i=0;i<3;++i)
        {
                op->mat[0]=0.0f;
                op->mat[1]=0.0f;
                op->mat[2]=0.0f;
                op->mat=1.0f;
        }
}

void CGraphicX::MatrixMul(Matrix_2_2 *op1, Matrix_2_2 *op2, Matrix_2_2 *result)
{
        Matrix_2_2 tmpMat;
        int i,j,k;
        for(i=0;i<3;++i)
                for(j=0;j<3;++j)
                        for(k=0;k<3;++k)
                        {
                                tmpMat.mat[j]+=(op1->mat[k])*(op2->mat[k][j]);
                        }
                        *result=tmpMat;
}

void CGraphicX::MatrixTrans(Matrix_2_2 *result, float x, float y)
{
        MatrixIdentility(result);
        result->mat[2][0]=x;
        result->mat[2][1]=y;
}

void CGraphicX::MatrixRotate(Matrix_2_2 *result,float rad)
{
        MatrixIdentility(result);
        result->mat[0][0]=cosf(rad);
        result->mat[0][1]=sinf(rad);
        result->mat[1][0]=(-1)*sinf(rad);
        result->mat[1][1]=cosf(rad);
}

void CGraphicX::MatrixScale(Matrix_2_2 *result, float x, float y)
{
        MatrixIdentility(result);
        result->mat[0][0]=x;
        result->mat[1][1]=y;
}

void CGraphicX::MatrixTiltedX(Matrix_2_2 *result, float c)
{
        MatrixIdentility(result);
        result->mat[1][0]=c;
}

void CGraphicX::MatrixTiltedY(Matrix_2_2 *result, float c)
{
        MatrixIdentility(result);
        result->mat[0][1]=c;
}

void CGraphicX::MatrixView(Matrix_2_2 *result,int width, int height, int left, int top, int right, int bottom)
{
        if(left<0)
                left=0;
        if(top<0)
                top=0;
        if(right<0)
                right=0;
        if(bottom<0)
                bottom=0;
        if(left>=width)
                left=width-1;
        if(top>=height)
                top=height-1;
        if(right>=width)
                right=width-1;
        if(bottom>=height)
                bottom=height-1;
        Matrix_2_2 matScale;
        Matrix_2_2 matTrans;
        float xScale=((float)right-(float)left)/(float)width;
        float yScale=((float)bottom-(float)top)/(float)height;
        MatrixScale(&matScale,xScale,yScale);
        MatrixTrans(&matTrans,(float)left,(float)top);
        *result=matScale*matTrans;
}

corpus

不可逆咋办？

qqj1228

LS说的对，我记得线性代数里说过不是每个矩阵都是可逆的。

reedlxm

不可逆的行列式的值为0