摄影测量解析相对定向及绝对定向原理(C++实现)

声明：本文非原创，原创为相对定向-绝对定向（C++实现）--摄影测量学，本文只是参照该文章代码学习了相对定向到绝对定向的实现的过程,在此对csuhan表示感谢！此外关于摄影测量解析相对定向及绝对定向原理及公式部分参见乖摸摸头_1951

原始数据为：
焦距：150mm
左片像点坐标：
0.016012 0.079963
0.08856 0.081134
0.013362 -0.07937
0.08224 -0.080027
0.051758 0.080555
0.014618 -0.000231
0.04988 -0.000782
0.08614 -0.001346
0.048035 -0.079962
右片像点坐标：
-0.07393 0.078706
-0.005252 0.078184
-0.079122 -0.078879
-0.009887 -0.080089
-0.039953 0.078463
-0.076006 0.000036
-0.042201 -0.001022
-0.007706 -0.002112
-0.044438 -0.079736
地面点坐标：
5083.205 5852.099 527.925
5780.02 5906.365 571.549
5210.879 4258.446 461.81
5909.264 4314.283 455.484
计算结果为

计算结果.png

#include "pch.h"
#include <iostream>
#include"Matrix.h"
#include"Rotation.h"
#include<fstream>
using  namespace std;
int main()
{
    int m1=9, m2=4;
    double f=0.15;
    double Left[9][2] = 
    {
        0.016012,0.079963,
        0.08856,0.081134,
        0.013362 ,-0.07937,
        0.08224 ,-0.080027,
        0.051758,0.080555,
        0.014618 ,- 0.000231,
        0.04988 ,-0.000782,
        0.08614,-0.001346,
        0.048035,-0.079962
    };
    double Right[9][2] =
    {
        -0.07393   ,0.078706,
        -0.005252  ,0.078184,
        -0.079122 ,-0.078879,
        -0.009887, -0.080089,
        -0.039953 ,0.078463,
        -0.076006  ,0.000036,
        -0.042201 ,-0.001022,
        -0.007706 ,-0.002112,
        -0.044438 ,-0.079736,
    };
    double Ground[4][3]=
    {
        5083.205 ,5852.099,527.925,
        5780.02  ,5906.365,571.549,
        5210.879 ,4258.446,461.81,
        5909.264 ,4314.283,455.484,
    };
    Model resolve(m1,Left, Right, m2, Ground, f);
    resolve.relative();//相对定位
    resolve.model();//求算模型点坐标
    resolve.absolute();//绝对定位
    resolve.transabs();
    return 0;
}

定向过程

#pragma once
#include "pch.h"
#include"Matrix.h"
class Model
{
public:
    Model(int m1, double Left[][2], double Right[][2], int m2, double Ground[][3],double f);
    ~Model();
    Matrix romatrix(double o, double p, double k);/*计算旋转矩阵*/
    void relative();//相对定位
    void model();//求算模型点坐标
    void absolute();//绝对定位
    void transabs();//模型点转换为地面坐标值

private:
                   /*左右相片量测数据\地面已知数据模型点数据*/
    int m;
    int n;
    double **L=NULL;
    double **R=NULL;
    double **G=NULL;
    double **MO=NULL;
    double **AB=NULL;
                              /*相机内方位元素*/
    double x0;
    double y0;
    double f0;
    /************************相对定向参数**********************************/
    //相对定向线元素
    double bx;
    double by;
    double bz;
    //相对定向角元素
    double o;
    double p;
    double k;
    double u;
    double v;
    /************************绝对定向参数********************************/
    //绝对定向线元素
    double X;
    double Y;
    double Z;
    //绝对定向角元素
    double O;
    double P;
    double K;
    //比例因子
    double M;
    double MA;
    //投影系数
    double N1;
    double N2;
    //精度因子
    double c1;
    double c2;
    //迭代次数
    int count1;
    int count2;
};
#include "pch.h"
#include"Rotation.h"
#include <math.h>
using namespace std;
/*构造函数赋初值*/
Model::Model(int m1, double Left[9][2], double Right[9][2], int m2, double Ground[9][3], double f)
{
    //焦距赋值
    this->f0 = f;
    //相对定向赋初值
    bx = Left[0][0] - Right[0][0];
    o = 0, p = 0, k = 0, u = 0, v = 0, count1 = 0;
    //绝对定向初始值
    X = 0, Y = 0, Z = 0, O = 0, P = 0, K = 0, MA = 1;
    //左右相片测量数据赋值
    m = m1;
    n = m2;
    L = new double*[m];
    R = new double*[m];
    G = new double*[n];
    MO = new double*[m];
    AB = new double*[m];
    for (int i = 0;i < m;i++)
    { 
        L[i] = new double[2];
        R[i] = new double[2];
        MO[i] = new double[3];
        AB[i] = new double[3];
    }
    for (int i = 0; i < n; i++)
    {
        G[i] = new double[3];
    }
    for (int i = 0; i < m; i++)
    {
        L[i][0] = Left[i][0];
        L[i][1] = Left[i][1];
        R[i][0] = Right[i][0];
        R[i][1] = Right[i][1];
    }
    //地面已知值
    for (int j = 0; j < n; j++)
    {
        G[j][0] = Ground[j][0];
        G[j][1] = Ground[j][1];
        G[j][2] = Ground[j][2];
    }
    //定义模型初始值
    for (int i = 0; i < m; i++)
    {
        MO[i][0] = 0;
        MO[i][1] = 0;
        MO[i][2] = 0;
        AB[i][0] = 0;
        AB[i][1] = 0;
        AB[i][2] = 0;

    }
}
/*析构函数*/
Model::~Model()
{
    for (int i = 0; i < m; i++)
    {
        delete[]L[i];
        delete[]R[i];
        delete[]MO[i];
        delete[]AB[i];
    }
    delete[]L;
    delete[]R;
    delete[]MO;
    delete[]AB;
    for (int i = 0; i < n; i++)
    {
        delete[]G[i];
    }
    delete[]G;
}
/*进行相对定向*/
void Model::relative()
{
    //定义法方程中的矩阵
    Matrix A(6, 5),Lc(6, 1),V(6, 1),Xa(5, 1), xfu(3, 1);;
    while (true)
    {
        by = bx * u, bz = bx * v;
        for (int i = 0; i < 6; i++)
        {
            xfu(0, 0) = R[i][0];
            xfu(1, 0) = R[i][1];
            xfu(2, 0) = -f0;
            xfu = romatrix(o, p, k)*xfu;
            //左片
            double X1 = L[i][0];
            double Y1 = L[i][1];
            double Z1 = -f0;
            //右片
            double X2 = xfu(0, 0);
            double Y2 = xfu(1, 0);
            double Z2 = xfu(2, 0);
            //计算投影系数
            N1 = (bx* Z2 - bz * X2) / (X1 * Z2 - X2*Z1);
            N2 = (bx* Z1 - bz *X1) / (X1 * Z2 - X2 * Z1);
            //计算法方程系数
            A(i, 0) = -X2 * Y2*N2 / Z2;
            A(i, 1) = -(Z2+Y2*Y2/Z2)*N2;
            A(i, 2) =X2*N2;
            A(i, 3) = bx;
            A(i, 4) = -(Y2 *bx) /Z2;
            //计算常数项
            Lc(i, 0) = N1 *Y1 - N2 * Y2 - by;
        }
        Xa = (~A*A).Inv()*(~A)*Lc;
        o += Xa(0, 0), p += Xa(1, 0), k += Xa(2, 0), u += Xa(3, 0), v += Xa(4, 0);
        count1++;
        if (fabs(Xa(0, 0))<0.00003&&fabs(Xa(1, 0))<0.00003 && fabs(Xa(2, 0))< 0.00003&&fabs(Xa(3, 0)) < 0.00003&&fabs(Xa(4, 0))< 0.00003)
        {
          cout << "Xa:\n" << Xa(0, 0) << "\t" << Xa(1, 0) << "\t" << Xa(2, 0) << "\t" << Xa(3, 0) << "\t" << Xa(4, 0) << endl;
          V = A * Xa - Lc;
          c1=sqrt((~V*V)(0,0)/6);
          cout << "相对定向的结果为：\n" << "o=" << o << "\t" << "p=" << p << "\t" << "k=" << k << "\t" << "u:" << u<< "\t" << "v:" <<v<< "\n";
          cout << "迭代的精度以及次数为:\n" << "精度为：" << c1 << "迭代次数为：" << count1 << "\n";
          break;
        }
    }
}
/*确定模型点的坐标*/
void Model::model()
{
    //确定模型比例尺M 
    M = sqrt(((G[0][0] - G[1][0])* (G[0][0] - G[1][0]) + (G[0][1] - G[1][1])* (G[0][1] - G[1][1])) / ((L[0][0] - L[1][0])* (L[0][0] - L[1][0]) + (L[0][1] - L[1][1])* (L[0][1] - L[1][1])));
    //计算模型的坐标
    cout << "模型点坐标为：\n";
    Matrix temp(3, 1);
    for (int i = 0; i < m; i++)
    {
        temp(0, 0) = R[i][0];
        temp(1, 0) = R[i][1];
        temp(2, 0) = -f0;
        temp = romatrix(o, p, k)*temp;
        //左片
        double X1 = L[i][0];
        double Y1 = L[i][1];
        double Z1 = -f0;
        //右片
        double X2 = temp(0, 0);
        double Y2 = temp(1, 0);
        double Z2 = temp(2, 0);
        //计算投影系数
        N1 = (bx* Z2 - bz * X2) / (X1 * Z2 - X2 * Z1);
        N2 = (bx* Z1 - bz * X1) / (X1 * Z2 - X2 * Z1);
        MO[i][0] = M * N1*X1;
        MO[i][1] = 0.5*M*(N1*Y1 + N2*Y2+by);
        MO[i][2] = M * (f0 + N1 * (-f0));
        cout <<  MO[i][0] << "\t" << MO[i][1] << "\t" << MO[i][2] << endl;
    }
}
//绝对定向
void Model::absolute()
{
    Matrix A(3*n,7), Xa(7,1), V(3*n,1), La(3*n,1);
    
    while(true)
    {
        int i = 0;
     for (int j=0; j < n; j++)
      { 
         Matrix temp(3,1);
         temp(0, 0) = MO[j][0], temp(1, 0) = MO[j][1], temp(2, 0) = MO[j][2];
        A(i, 0) = 1,A(i, 3) =MO[j][0], A(i, 4) = -MO[j][2], A(i, 6) = -MO[j][1];
        A(i + 1, 1) = 1, A(i + 1, 3) = MO[j][1], A(i+1, 5) = -MO[j][2], A(i+1, 6) = MO[j][0];
        A(i + 2, 2) = 1, A(i + 2, 3) = MO[j][2], A(i+2, 4) = MO[j][0], A(i + 2, 5) = MO[j][1];
        temp=romatrix(O, P, K)*temp;
        La(i, 0) = G[j][0] - MA * temp(0, 0) - X;
        La(i + 1, 0) = G[j][1] - MA * temp(1, 0) - Y;
        La(i + 2, 0) = G[j][2] - MA * temp(2, 0) - Z;
        i += 3;
      }
      Xa= (~A*A).Inv()*(~A)*La;
      X += Xa(0, 0), Y += Xa(1, 0), Z += Xa(2, 0), O += Xa(4, 0), P += Xa(5, 0), K += Xa(6, 0), MA += Xa(3, 0);
      count2++;
     if (fabs(Xa(0, 0)) < 0.001&&fabs(Xa(1, 0)) < 0.001 && fabs(Xa(2, 0)) < 0.001&&fabs(Xa(3,0)) < 0.001&&fabs(Xa(4,0)) < 0.001&&fabs(Xa(5, 0)) < 0.001&&fabs(Xa(6, 0)) < 0.001)
      {
          V = A * Xa - La;
          c2 = sqrt((~V*V)(0, 0) / (3*n-7));
          cout << "绝对定向的结果：\n" << "X:" << X << "\t" << "Y:" << Y << "\t" << "Z:" << Z << "\t" << "MA:" << MA << "\t" << "O:" << O << "\t" << "P:" << P << "\t" << "K:" << K << "\n";
          cout << "绝对定向精度与迭代的次数为：\n" << "精度：" << c2 << "\t" << "迭代次数：" << count2 << "\n";
          break;
      }
    }
}
//计算模型点在地面的坐标
void Model::transabs()
{
    cout << "转换为的地面点的坐标为：\n";
    for(int i = 0; i < m; i++)
    {
        Matrix temp(3, 1);
        temp(0, 0) = MO[i][0];
        temp(1, 0) = MO[i][1];
        temp(2, 0) = MO[i][2];
        temp=romatrix(O, P, K)*temp;
        AB[i][0] = MA * temp(0, 0) + X;
        AB[i][1] = MA * temp(1, 0) + Y;
        AB[i][2] = MA * temp(2, 0) + Z;
        cout << AB[i][0] << "\t" << AB[i][1] << "\t" << AB[i][2] << "\n";
    }
}
//计算旋转矩阵
Matrix Model::romatrix(double o, double p, double k)
{
    Matrix rota(3, 3);
    rota(0, 0) = cos(o)*cos(k) - sin(o)*sin(p)*sin(k);
    rota(0, 1) = -cos(o)*sin(k) - sin(o)*sin(p)*cos(k);
    rota(0, 2) = -sin(o)*cos(p);
    rota(1, 0) = cos(p)*sin(k);
    rota(1, 1) = cos(p)*cos(k);
    rota(1, 2) = -sin(p);
    rota(2, 0) = sin(o)*cos(k) + cos(o)*sin(p)*sin(k);
    rota(2, 1) = -sin(o)*sin(k) + cos(o)*sin(p)*cos(k);
    rota(2, 2) = cos(o)*cos(p);
    return rota;
}

矩阵类

#include "pch.h"
#include"Rotation.h"
#include <math.h>
using namespace std;
/*构造函数赋初值*/
Model::Model(int m1, double Left[9][2], double Right[9][2], int m2, double Ground[9][3], double f)
{
    //焦距赋值
    this->f0 = f;
    //相对定向赋初值
    bx = Left[0][0] - Right[0][0];
    o = 0, p = 0, k = 0, u = 0, v = 0, count1 = 0;
    //绝对定向初始值
    X = 0, Y = 0, Z = 0, O = 0, P = 0, K = 0, MA = 1;
    //左右相片测量数据赋值
    m = m1;
    n = m2;
    L = new double*[m];
    R = new double*[m];
    G = new double*[n];
    MO = new double*[m];
    AB = new double*[m];
    for (int i = 0;i < m;i++)
    { 
        L[i] = new double[2];
        R[i] = new double[2];
        MO[i] = new double[3];
        AB[i] = new double[3];
    }
    for (int i = 0; i < n; i++)
    {
        G[i] = new double[3];
    }
    for (int i = 0; i < m; i++)
    {
        L[i][0] = Left[i][0];
        L[i][1] = Left[i][1];
        R[i][0] = Right[i][0];
        R[i][1] = Right[i][1];
    }
    //地面已知值
    for (int j = 0; j < n; j++)
    {
        G[j][0] = Ground[j][0];
        G[j][1] = Ground[j][1];
        G[j][2] = Ground[j][2];
    }
    //定义模型初始值
    for (int i = 0; i < m; i++)
    {
        MO[i][0] = 0;
        MO[i][1] = 0;
        MO[i][2] = 0;
        AB[i][0] = 0;
        AB[i][1] = 0;
        AB[i][2] = 0;

    }
}
/*析构函数*/
Model::~Model()
{
    for (int i = 0; i < m; i++)
    {
        delete[]L[i];
        delete[]R[i];
        delete[]MO[i];
        delete[]AB[i];
    }
    delete[]L;
    delete[]R;
    delete[]MO;
    delete[]AB;
    for (int i = 0; i < n; i++)
    {
        delete[]G[i];
    }
    delete[]G;
}
/*进行相对定向*/
void Model::relative()
{
    //定义法方程中的矩阵
    Matrix A(6, 5),Lc(6, 1),V(6, 1),Xa(5, 1), xfu(3, 1);;
    while (true)
    {
        by = bx * u, bz = bx * v;
        for (int i = 0; i < 6; i++)
        {
            xfu(0, 0) = R[i][0];
            xfu(1, 0) = R[i][1];
            xfu(2, 0) = -f0;
            xfu = romatrix(o, p, k)*xfu;
            //左片
            double X1 = L[i][0];
            double Y1 = L[i][1];
            double Z1 = -f0;
            //右片
            double X2 = xfu(0, 0);
            double Y2 = xfu(1, 0);
            double Z2 = xfu(2, 0);
            //计算投影系数
            N1 = (bx* Z2 - bz * X2) / (X1 * Z2 - X2*Z1);
            N2 = (bx* Z1 - bz *X1) / (X1 * Z2 - X2 * Z1);
            //计算法方程系数
            A(i, 0) = -X2 * Y2*N2 / Z2;
            A(i, 1) = -(Z2+Y2*Y2/Z2)*N2;
            A(i, 2) =X2*N2;
            A(i, 3) = bx;
            A(i, 4) = -(Y2 *bx) /Z2;
            //计算常数项
            Lc(i, 0) = N1 *Y1 - N2 * Y2 - by;
        }
        Xa = (~A*A).Inv()*(~A)*Lc;
        o += Xa(0, 0), p += Xa(1, 0), k += Xa(2, 0), u += Xa(3, 0), v += Xa(4, 0);
        count1++;
        if (fabs(Xa(0, 0))<0.00003&&fabs(Xa(1, 0))<0.00003 && fabs(Xa(2, 0))< 0.00003&&fabs(Xa(3, 0)) < 0.00003&&fabs(Xa(4, 0))< 0.00003)
        {
          cout << "Xa:\n" << Xa(0, 0) << "\t" << Xa(1, 0) << "\t" << Xa(2, 0) << "\t" << Xa(3, 0) << "\t" << Xa(4, 0) << endl;
          V = A * Xa - Lc;
          c1=sqrt((~V*V)(0,0)/6);
          cout << "相对定向的结果为：\n" << "o=" << o << "\t" << "p=" << p << "\t" << "k=" << k << "\t" << "u:" << u<< "\t" << "v:" <<v<< "\n";
          cout << "迭代的精度以及次数为:\n" << "精度为：" << c1 << "迭代次数为：" << count1 << "\n";
          break;
        }
    }
}
/*确定模型点的坐标*/
void Model::model()
{
    //确定模型比例尺M 
    M = sqrt(((G[0][0] - G[1][0])* (G[0][0] - G[1][0]) + (G[0][1] - G[1][1])* (G[0][1] - G[1][1])) / ((L[0][0] - L[1][0])* (L[0][0] - L[1][0]) + (L[0][1] - L[1][1])* (L[0][1] - L[1][1])));
    //计算模型的坐标
    cout << "模型点坐标为：\n";
    Matrix temp(3, 1);
    for (int i = 0; i < m; i++)
    {
        temp(0, 0) = R[i][0];
        temp(1, 0) = R[i][1];
        temp(2, 0) = -f0;
        temp = romatrix(o, p, k)*temp;
        //左片
        double X1 = L[i][0];
        double Y1 = L[i][1];
        double Z1 = -f0;
        //右片
        double X2 = temp(0, 0);
        double Y2 = temp(1, 0);
        double Z2 = temp(2, 0);
        //计算投影系数
        N1 = (bx* Z2 - bz * X2) / (X1 * Z2 - X2 * Z1);
        N2 = (bx* Z1 - bz * X1) / (X1 * Z2 - X2 * Z1);
        MO[i][0] = M * N1*X1;
        MO[i][1] = 0.5*M*(N1*Y1 + N2*Y2+by);
        MO[i][2] = M * (f0 + N1 * (-f0));
        cout <<  MO[i][0] << "\t" << MO[i][1] << "\t" << MO[i][2] << endl;
    }
}
//绝对定向
void Model::absolute()
{
    Matrix A(3*n,7), Xa(7,1), V(3*n,1), La(3*n,1);
    
    while(true)
    {
        int i = 0;
     for (int j=0; j < n; j++)
      { 
         Matrix temp(3,1);
         temp(0, 0) = MO[j][0], temp(1, 0) = MO[j][1], temp(2, 0) = MO[j][2];
        A(i, 0) = 1,A(i, 3) =MO[j][0], A(i, 4) = -MO[j][2], A(i, 6) = -MO[j][1];
        A(i + 1, 1) = 1, A(i + 1, 3) = MO[j][1], A(i+1, 5) = -MO[j][2], A(i+1, 6) = MO[j][0];
        A(i + 2, 2) = 1, A(i + 2, 3) = MO[j][2], A(i+2, 4) = MO[j][0], A(i + 2, 5) = MO[j][1];
        temp=romatrix(O, P, K)*temp;
        La(i, 0) = G[j][0] - MA * temp(0, 0) - X;
        La(i + 1, 0) = G[j][1] - MA * temp(1, 0) - Y;
        La(i + 2, 0) = G[j][2] - MA * temp(2, 0) - Z;
        i += 3;
      }
      Xa= (~A*A).Inv()*(~A)*La;
      X += Xa(0, 0), Y += Xa(1, 0), Z += Xa(2, 0), O += Xa(4, 0), P += Xa(5, 0), K += Xa(6, 0), MA += Xa(3, 0);
      count2++;
     if (fabs(Xa(0, 0)) < 0.001&&fabs(Xa(1, 0)) < 0.001 && fabs(Xa(2, 0)) < 0.001&&fabs(Xa(3,0)) < 0.001&&fabs(Xa(4,0)) < 0.001&&fabs(Xa(5, 0)) < 0.001&&fabs(Xa(6, 0)) < 0.001)
      {
          V = A * Xa - La;
          c2 = sqrt((~V*V)(0, 0) / (3*n-7));
          cout << "绝对定向的结果：\n" << "X:" << X << "\t" << "Y:" << Y << "\t" << "Z:" << Z << "\t" << "MA:" << MA << "\t" << "O:" << O << "\t" << "P:" << P << "\t" << "K:" << K << "\n";
          cout << "绝对定向精度与迭代的次数为：\n" << "精度：" << c2 << "\t" << "迭代次数：" << count2 << "\n";
          break;
      }
    }
}
//计算模型点在地面的坐标
void Model::transabs()
{
    cout << "转换为的地面点的坐标为：\n";
    for(int i = 0; i < m; i++)
    {
        Matrix temp(3, 1);
        temp(0, 0) = MO[i][0];
        temp(1, 0) = MO[i][1];
        temp(2, 0) = MO[i][2];
        temp=romatrix(O, P, K)*temp;
        AB[i][0] = MA * temp(0, 0) + X;
        AB[i][1] = MA * temp(1, 0) + Y;
        AB[i][2] = MA * temp(2, 0) + Z;
        cout << AB[i][0] << "\t" << AB[i][1] << "\t" << AB[i][2] << "\n";
    }
}
//计算旋转矩阵
Matrix Model::romatrix(double o, double p, double k)
{
    Matrix rota(3, 3);
    rota(0, 0) = cos(o)*cos(k) - sin(o)*sin(p)*sin(k);
    rota(0, 1) = -cos(o)*sin(k) - sin(o)*sin(p)*cos(k);
    rota(0, 2) = -sin(o)*cos(p);
    rota(1, 0) = cos(p)*sin(k);
    rota(1, 1) = cos(p)*cos(k);
    rota(1, 2) = -sin(p);
    rota(2, 0) = sin(o)*cos(k) + cos(o)*sin(p)*sin(k);
    rota(2, 1) = -sin(o)*sin(k) + cos(o)*sin(p)*cos(k);
    rota(2, 2) = cos(o)*cos(p);
    return rota;
}
#include "pch.h"
#include<iostream>
#include"Matrix.h"
using namespace std;
//默认构造函数
Matrix::Matrix(int r, int c)
{
    Row = r;
    Col = c;
    Data = new double*[Row];
    for (int i = 0; i < Row; i++)
    {
        Data[i] = new double[Col];
        for (int j = 0; j < Col; j++)
        {
            Data[i][j] = 0;
        }
    }
}
/*复制构造函数*/
Matrix::Matrix(const Matrix&ma)
{
    Row = ma.R();
    Col = ma.C();
    Data = new double *[Row];
    for (int i = 0; i < Row; i++)
    {
        Data[i] = new double[Col];
        memcpy(Data[i], ma.Data[i], sizeof(double)*Col);
    }
}
/*析构函数*/
Matrix::~Matrix()
{
    for (int i = 0; i < Row; i++)
    {
        delete[]Data[i];
    }
    delete[]Data;
}
/*重载运算符（）*/
double &Matrix::operator()(int a, int b)
{
    if (a >= Row || b >= Col)
        throw("Matrix::operator:Index out of range!");
    return Data[a][b];
}
/*常量对象重载运算符（）*/
double &Matrix::operator()(int a, int b)const
{
    if (a >= Row || b >= Col)
        throw("Matrix::operator:Index out of range!");
    return Data[a][b];
}
/*重载=运算符*/
Matrix Matrix ::operator =(const Matrix &ma1)
{
    if (Row != ma1.R() || Col != ma1.C())
    {
        throw"the matrix are different!";
    }
    for (int i = 0; i < Row; i++)
    {
        for (int j = 0; j < Col; j++)
        {
            Data[i][j] = ma1(i, j);
        }
    }
    return *this;
}
/*重载+运算符*/
Matrix operator+(const Matrix &ma1, const Matrix &ma2)
{
    if (ma1.R() != ma2.R() || ma1.C() != ma2.C())
        throw("the size of two matrix is different");
    Matrix ma(ma1.R(), ma1.C());
    for (int i = 0; i < ma1.R(); i++)
    {
        for (int j = 0; j < ma1.C(); j++)
        {
            ma(i, j) = ma1(i, j) + ma2(i, j);
        }
    }
    return ma;
}
/*重载-运算符*/
Matrix operator-(const Matrix &ma1, const Matrix &ma2)
{
    if (ma1.R() != ma2.R() || ma1.C() != ma2.C())
        throw("the size of two matrix is different");
    Matrix ma(ma1.R(), ma1.C());
    for (int i = 0; i < ma1.R(); i++)
    {
        for (int j = 0; j < ma1.C(); j++)
        {
            ma(i, j) = ma1(i, j) - ma2(i, j);
        }
    }
    return ma;
}
/*重载*运算符--矩阵相乘*/
Matrix operator * (const Matrix &ma1, const Matrix &ma2)
{
    if (ma1.C() != ma2.R())
        throw("There two matrix can not *");
    Matrix ma(ma1.R(), ma2.C());
    for (int i = 0; i < ma1.R(); i++)
    {
        for (int j = 0; j < ma2.C(); j++)
        {
            for (int k = 0; k < ma1.C(); k++)
            {
                ma(i, j) += ma1(i, k)*ma2(k, j);
            }
        }
    }
    return ma;
}
/*重载*运算符--数乘矩阵左乘*/
Matrix operator *(int a, const Matrix &ma1)
{
    Matrix ma(ma1.R(), ma1.C());
    for (int i = 0; i < ma1.R(); i++)
    {
        for (int j = 0; j < ma1.C(); j++)
        {
            ma(i, j) = a * ma1(i, j);
        }
    }
    return ma;
}
/*重载*运算符--数乘矩阵右乘*/
Matrix operator *(const Matrix &ma1, int a)
{
    Matrix ma(ma1.R(), ma1.C());
    for (int i = 0; i < ma1.R(); i++)
    {
        for (int j = 0; j < ma1.C(); j++)
        {
            ma(i, j) = a * ma1(i, j);
        }
    }
    return ma;
}
/*重载~运算符--矩阵转置*/
Matrix operator ~(const Matrix &ma1)
{
    Matrix ma(ma1.C(), ma1.R());
    for (int i = 0; i < ma1.R(); i++)
    {
        for (int j = 0; j < ma1.C(); j++)
        {
            ma(j, i) = ma1(i, j);
        }
    }
    return ma;
}
/*矩阵求逆----全选主元法*/
Matrix Matrix::Inv()
{
    if (Row != Col)
    {
        throw("The matrix can not be inversion!");
    }
    Matrix cmat(Row, Row);
    cmat = *this;//复制对象
    int i, j, k, v;
    int *mainrow = new int[Row];//存放主元的行号
    int *maincol = new int[Row];//存放主元的列号
    double mid;                 //临时变量
    double maincount;          //主元值  
    /*选取主元*/
    for (k = 0; k < Row; k++)
    {
        maincount = 0;
        for (i = k; i < Row; i++)
        {
            for (j = k; j < Col; j++)
            {
                if (fabs(cmat(i, j)) >= maincount)
                {
                    maincount = fabs(cmat(i, j));
                    mainrow[k] = i;
                    maincol[k] = j;
                }

            }
        }
        if (fabs(maincount) < 0.0000000000001)
        {
             throw("The matrix can not be inversion!");
        }
        /*行交换*/
        if (mainrow[k] != k)
        {
            for (i = 0; i < Row; i++)
            {
                v = mainrow[k];
                mid = cmat(k, i);
                cmat(k, i) = cmat(v, i);
                cmat(v, i) = mid;
            }
        }
        /*列交换*/
        if (maincol[k] != k)
        {
            for (i = 0; i < Row; i++)
            {
                v = maincol[k];
                mid = cmat(i, k);
                cmat(i, k) = cmat(i, v);
                cmat(i, v) = mid;
            }
        }
        /*----计算cmat(k, k)*/
        cmat(k, k) = 1 / cmat(k, k);
        /*----计算cmat(k,j)j!=k*/
        for (j = 0; j < Row; j++)
        {
            if (j != k)
            {
                cmat(k, j) *= cmat(k, k);
            }
        }
        /*----计算cmat(i,j)i,j!=k*/
        for (i = 0; i < Row; i++)
        {
            if (i != k)
            {
                for (j = 0; j < Row; j++)
                {
                    if (j != k)
                        cmat(i, j) -= cmat(i, k)*cmat(k, j);
                }
            }
        }
        /*计算cmat(i,k)i!=k*/
        for (i = 0; i < Row; i++)
        {
            if (i != k)
                cmat(i, k) = -cmat(i, k)*cmat(k, k);
        }
    }
    /*交换行列*/
    for (k = Row - 1; k >= 0; k--)
    {
        for (i = 0; i < Row; i++)
        {
            v = maincol[k];
            mid = cmat(k, i);
            cmat(k, i) = cmat(v, i);
            cmat(v, i) = mid;

        }
        for (i = 0; i < Row; i++)
        {
            v = mainrow[k];
            mid = cmat(i, k);
            cmat(i, k) = cmat(i, v);
            cmat(i, v) = mid;
        }
    }
    delete[]mainrow;
    delete[]maincol;
    return cmat;
}
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 204,732评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 87,496评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 151,264评论 0赞 338
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,807评论 1赞 277
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,806评论 5赞 368
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,675评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,029评论 3赞 399
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,683评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 41,704评论 1赞 299
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,666评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,773评论 1赞 332
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,413评论 4赞 321
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,016评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,978评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,204评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,083评论 2赞 350
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,503评论 2赞 343

推荐阅读更多精彩内容