opencv实现正交匹配追踪算法OMP

//dic: 字典矩阵；
//signal :待重构信号（一次只能重构一个信号，即一个向量）
//min_residual: 最小残差
//sparsity:稀疏度
//coe:重构系数
//atom_index:字典原子选择序号
//返回最后的残差
float OMP(    Mat& dic,Mat& signal,float min_residual,int sparsity,Mat& coe,vector<int>& atom_index)

{
    if(signal.cols>1)
    {
        cout<<"wrong signal"<<endl;
        return -1;
    }
    signal=signal/norm(signal);  //信号单位化
    Mat temp(1,dic.cols,5);
    for(int i=0;i<dic.cols;i++)
    {
        temp.col(i)=norm(dic.col(i));  //每个原子的模长
    }
    divide(dic,repeat(temp,dic.rows,1),dic); //字典原子单位化
    Mat residual = signal.clone();  //初始化残差
    coe.create(0, 1, CV_32FC1);  //初始化系数
    Mat phi;    //保存已选出的原子向量
    float max_coefficient;
    unsigned int atom_id;  //每次所选择的原子的序号

    for(;;)
    {
        max_coefficient = 0;
        //取出内积最大列
        for (int i = 0; i <dic.cols; i++)
        {
            float coefficient = (float)dic.col(i).dot(residual); 

            if (abs(coefficient) > abs(max_coefficient))
            {
            max_coefficient = coefficient;
            atom_id = i;
            }
        }
        atom_index.push_back(atom_id); //添加选出的原子序号        
        Mat& temp_atom= dic.col(atom_id); //取出该原子
        if (phi.cols == 0)
            phi = temp_atom;
        else
            hconcat(phi,temp_atom,phi); //将新原子合并到原子集合中（都是列向量）
        
        coe.push_back(0.0f);    //对系数矩阵新加一项
        solve(phi, signal,coe, DECOMP_SVD);    //求解最小二乘问题

        residual = signal - phi*coe;  //更新残差
        float res_norm = (float)norm(residual);

        if (coe.rows >= sparsity || res_norm <= min_residual) //如果残差小于阈值或达到要求的稀疏度，就返回
        {
            return res_norm;
        }
    }
}