SVD分解 opencv实现

于 2022-12-13 11:27:50 发布
阅读量730 收藏 1
点赞数 1
文章标签： opencv
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_71792169/article/details/128299339
版权

头文件

#ifndef DEBUG_LRN_SVD_H
#define DEBUG_LRN_SVD_H
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
#include <opencv2/opencv.hpp>
#include </usr/local/include/eigen3/Eigen/Dense>

using namespace std;
using namespace cv;

int test_SVD();

int opencv_svd(vector<vector<float>>& vec);

int eigen_svd(vector<vector<float>>& vec);

int cpp_svd(vector<vector<float>>& vec);

void print_matrix(const vector<vector<float>>& vec);

#endif //DEBUG_LRN_SVD_H

源文件

#include "svd.h"

void print_matrix(const vector<vector<float>>& vec){
    if(vec.empty()){
        return;
    }
    for ( auto row: vec){
        if(row.empty()){
            return;
        }
        for ( auto elem: row){
            cout << elem << ", ";
        }
        cout << endl;
    }
    cout << endl;
}

int opencv_svd(vector<vector<float>>& vec){
    cout << "source matrix：" << endl;
    print_matrix(vec);

    cout << "opencv singular value decomposition:" << endl;
    const auto rows = (int)vec.size();
    const auto cols = (int)vec[0].size();
    cv::Mat mat(rows, cols, CV_32FC1);
    for (int y = 0; y < rows; ++y) {
        for (int x = 0; x < cols; ++x) {
            mat.at<float>(y, x) = vec.at(y).at(x);
        }
    }

    cv::Mat matD, matU, matVt;
    cv::SVD::compute(mat, matD, matU, matVt, 4);
    cout << "opencv singular values:" << endl;
    cout << matD << endl;
    cout << "left singular vectors:" << endl;
    cout << matU << endl;
    cout << "transposed matrix of right singular values:" << endl;
    cout << matVt << endl;
//    cv::Mat matUt, matV;
//    cv::transpose(matU, matUt);
//    cv::transpose(matVt, matV);
//
//    cout << "verify if the result is correct:" << endl;
//    cv::Mat reconMatD = matUt * mat * matV;
//    cout << reconMatD << endl;
    return 0;
}


int eigen_svd(vector<vector<float>>& vec){
//    cout << "source matrix：" << endl;
//    print_matrix(vec);

    cout << "opencv singular value decomposition:" << endl;
    const auto rows = (int)vec.size();
    const auto cols = (int)vec[0].size();

    std::vector<float> vec_;
    for (int i = 0; i < rows; ++i) {
        vec_.insert(vec_.begin() + i * cols, vec[i].begin(), vec[i].end());
    }

    Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> m(vec_.data(), rows, cols);

    cout << "source matrix:" << endl;
    std::cout << m << std::endl;

    Eigen::JacobiSVD<Eigen::MatrixXf> svd(m, Eigen::ComputeFullV | Eigen::ComputeFullU); // ComputeThinU | ComputeThinV
    Eigen::MatrixXf singular_values = svd.singularValues();
    Eigen::MatrixXf left_singular_vectors = svd.matrixU();
    Eigen::MatrixXf right_singular_vectors = svd.matrixV();

    cout << "eigen singular values:" << singular_values << endl;
    cout << "left singular vectors:" << left_singular_vectors << endl;
    cout << "right singular vectors:" << right_singular_vectors << endl;
    return 0;
}



// ================================= 矩阵奇异值分解 =================================
template<typename _Tp>
static void JacobiSVD(std::vector<std::vector<_Tp>>& At,
                      std::vector<std::vector<_Tp>>& _W, std::vector<std::vector<_Tp>>& Vt)
{
    double minval = FLT_MIN;
    auto eps = (_Tp)(FLT_EPSILON * 2);
    const int m = At[0].size();
    const auto n = (int)_W.size();
    const int n1 = m; // urows
    std::vector<double> W(_W.size(), 0.);

    for (int i = 0; i < n; i++) {
        double sd{0.};
        for (int k = 0; k < m; k++) {
            _Tp t = At[i][k];
            sd += (double)t*t;
        }
        W[i] = sd;

        for (int k = 0; k < n; k++)
            Vt[i][k] = 0;
        Vt[i][i] = 1;
    }

    int max_iter = std::max(m, 30);
    for (int iter = 0; iter < max_iter; iter++) {
        bool changed = false;
        _Tp c, s;

        for (int i = 0; i < n - 1; i++) {
            for (int j = i + 1; j < n; j++) {
                _Tp *Ai = At[i].data(), *Aj = At[j].data();
                double a = W[i], p = 0, b = W[j];

                for (int k = 0; k < m; k++)
                    p += (double)Ai[k] * Aj[k];

                if (std::abs(p) <= eps * std::sqrt((double)a*b))
                    continue;

                p *= 2;
                double beta = a - b, gamma = hypot((double)p, beta);
                if (beta < 0) {
                    double delta = (gamma - beta)*0.5;
                    s = (_Tp)std::sqrt(delta / gamma);
                    c = (_Tp)(p / (gamma*s * 2));
                } else {
                    c = (_Tp)std::sqrt((gamma + beta) / (gamma * 2));
                    s = (_Tp)(p / (gamma*c * 2));
                }

                a = b = 0;
                for (int k = 0; k < m; k++) {
                    _Tp t0 = c*Ai[k] + s*Aj[k];
                    _Tp t1 = -s*Ai[k] + c*Aj[k];
                    Ai[k] = t0; Aj[k] = t1;

                    a += (double)t0*t0; b += (double)t1*t1;
                }
                W[i] = a; W[j] = b;

                changed = true;

                _Tp *Vi = Vt[i].data(), *Vj = Vt[j].data();

                for (int k = 0; k < n; k++) {
                    _Tp t0 = c*Vi[k] + s*Vj[k];
                    _Tp t1 = -s*Vi[k] + c*Vj[k];
                    Vi[k] = t0; Vj[k] = t1;
                }
            }
        }

        if (!changed)
            break;
    }

    for (int i = 0; i < n; i++) {
        double sd{ 0. };
        for (int k = 0; k < m; k++) {
            _Tp t = At[i][k];
            sd += (double)t*t;
        }
        W[i] = std::sqrt(sd);
    }

    for (int i = 0; i < n - 1; i++) {
        int j = i;
        for (int k = i + 1; k < n; k++) {
            if (W[j] < W[k])
                j = k;
        }
        if (i != j) {
            std::swap(W[i], W[j]);

            for (int k = 0; k < m; k++)
                std::swap(At[i][k], At[j][k]);

            for (int k = 0; k < n; k++)
                std::swap(Vt[i][k], Vt[j][k]);
        }
    }

    for (int i = 0; i < n; i++)
        _W[i][0] = (_Tp)W[i];

    srand(time(nullptr));

    for (int i = 0; i < n1; i++) {
        double sd = i < n ? W[i] : 0;

        for (int ii = 0; ii < 100 && sd <= minval; ii++) {
            // if we got a zero singular value, then in order to get the corresponding left singular vector
            // we generate a random vector, project it to the previously computed left singular vectors,
            // subtract the projection and normalize the difference.
            const auto val0 = (_Tp)(1. / m);
            for (int k = 0; k < m; k++) {
                unsigned long int rng = rand() % 4294967295; // 2^32 - 1
                _Tp val = (rng & 256) != 0 ? val0 : -val0;
                At[i][k] = val;
            }
            for (int iter = 0; iter < 2; iter++) {
                for (int j = 0; j < i; j++) {
                    sd = 0;
                    for (int k = 0; k < m; k++)
                        sd += At[i][k] * At[j][k];
                    _Tp asum = 0;
                    for (int k = 0; k < m; k++) {
                        auto t = (_Tp)(At[i][k] - sd*At[j][k]);
                        At[i][k] = t;
                        asum += std::abs(t);
                    }
                    asum = asum > eps * 100 ? 1 / asum : 0;
                    for (int k = 0; k < m; k++)
                        At[i][k] *= asum;
                }
            }

            sd = 0;
            for (int k = 0; k < m; k++) {
                _Tp t = At[i][k];
                sd += (double)t*t;
            }
            sd = std::sqrt(sd);
        }

        _Tp s = (_Tp)(sd > minval ? 1 / sd : 0.);
        for (int k = 0; k < m; k++)
            At[i][k] *= s;
    }
}

// matSrc为原始矩阵，支持非方阵，matD存放奇异值，matU存放左奇异向量，matVt存放转置的右奇异向量
template<typename _Tp>
int svd(const std::vector<std::vector<_Tp>>& matSrc,
        std::vector<std::vector<_Tp>>& matD, std::vector<std::vector<_Tp>>& matU, std::vector<std::vector<_Tp>>& matVt)
{
    int m = matSrc.size();
    int n = matSrc[0].size();
    for (const auto& sz : matSrc) {
        if (n != sz.size()) {
            fprintf(stderr, "matrix dimension dismatch\n");
            return -1;
        }
    }

    bool at = false;
    if (m < n) {
        std::swap(m, n);
        at = true;
    }

    matD.resize(n);
    for (int i = 0; i < n; ++i) {
        matD[i].resize(1, (_Tp)0);
    }
    matU.resize(m);
    for (int i = 0; i < m; ++i) {
        matU[i].resize(m, (_Tp)0);
    }
    matVt.resize(n);
    for (int i = 0; i < n; ++i) {
        matVt[i].resize(n, (_Tp)0);
    }
    std::vector<std::vector<_Tp>> tmp_u = matU, tmp_v = matVt;

    std::vector<std::vector<_Tp>> tmp_a, tmp_a_;
    if (!at)
        transpose(matSrc, tmp_a);
    else
        tmp_a = matSrc;

    if (m == n) {
        tmp_a_ = tmp_a;
    } else {
        tmp_a_.resize(m);
        for (int i = 0; i < m; ++i) {
            tmp_a_[i].resize(m, (_Tp)0);
        }
        for (int i = 0; i < n; ++i) {
            tmp_a_[i].assign(tmp_a[i].begin(), tmp_a[i].end());
        }
    }
    JacobiSVD(tmp_a_, matD, tmp_v);

    if (!at) {
        transpose(tmp_a_, matU);
        matVt = tmp_v;
    } else {
//        transpose(tmp_v, matVt);
//        matU = tmp_a_;
    }
    return 0;
}

int cpp_svd(vector<vector<float>>& vec){
    std::vector<std::vector<float>> matD, matU, matVt;
    if (svd(vec, matD, matU, matVt) != 0) {
        fprintf(stderr, "C++ implement singular value decomposition fail\n");
        return -1;
    }
//    cout << "singular values:" << endl;
//    print_matrix(matD);
//    cout << "left singular vectors" << endl;
//    print_matrix(matU);
//    cout << "transposed matrix of right singular values:" << endl;
//    print_matrix(matVt);
    return 0;
}


int test_SVD()
{
    std::vector<std::vector<float>> vec{ { 0.68f, 0.597f, 0.2752f, 0.68f, 0.597f, 0.2752f },
                                         { -0.211f, 0.823f, -0.9273f, 0.68f, 0.597f, 0.2752f },
                                         { 0.566f, -0.605f, 0.1260f, 0.68f, 0.597f, 0.2752f } };
    auto time = (double)getTickCount();
    const int repeat = 100;
    for(int i = 0;i<repeat;i++){
        opencv_svd(vec); // 对比下来这个最快
    }
    auto time2 = ((double)getTickCount() - time) / (repeat * getTickFrequency());
    for(int i = 0;i<repeat;i++){
        eigen_svd(vec);
    }
    auto time3 = ((double)getTickCount() - time2) / (repeat * getTickFrequency());
    for(int i = 0;i<repeat;i++){
        cpp_svd(vec);
    }
    auto time4 = ((double)getTickCount() - time3) / (repeat * getTickFrequency());
    time = 0;
    return 0;
}
OpenCV学习之SVD技术进行图像分解与重构
AutoXTruck专栏
04-21 2680
简述 　　上一篇博文算法学习之SVD理论推导介绍了SVD矩阵分解的完整过程，感兴趣的可以去阅读一下。本文以图像为基础，利用SVD理论对二维图像矩阵进行分解，同时选取不同的特征值个数来对图像进行重构。主要使用OpenCV第三方视觉库完成SVD应用过程，原谅我这个调包小虾米。 　　使用OpenCV库的SVD功能API函数时，简短说一下SVD的分解矩阵，我们通过对图像矩阵AAA(m×nm×nm×n，其中...
OpenCV——相当通俗易懂的SVD奇异值分解
coco的博客
08-14 9125
在说明前，我也是查了大量文档，弄清楚各个名词的意思，才写下这篇博客。。。 特征值和特征向量： 根据公式： A是n*n的方阵（必须是方阵），x是特征向量，λ是特征值。 一般情况下，会有n个特征值，和n个特征向量。（这里的一般是指方阵是满秩，各行各列都是线性无关）。 引出问题：如果不是n*n维的矩阵，怎么求？ 假设矩阵X是m*n，则可以求或矩阵的特征分解。 一般的非方阵的特征分解称为奇异...
opencv SVD分解
11-03 6920
调用opencv  cvSvd函数对矩阵分解之后，形成的矩阵还matlab分解之后得到的矩阵有所不同，matlab对矩阵A分解之后的三个矩阵分别为U   W   V，而V的最后一列即为矩阵的最小特征值对应的特征向量，也是AX=0的解。然而，opencv对矩阵A分解之后所得的矩阵不光顺序与matlab分解的有所不同，而且某些矩阵元素会相差一个负号，比如opencv分解出来的矩阵顺序是W U V,W和
OpenCV中使用SVD分解与重构
dianmocui5883的博客
11-07 1158
OpenCV中SVD分解函数compute C++: static void SVD::compute(InputArray src, OutputArray w, OutputArray u, OutputArray vt, int flags=0 ) src – Decomposed matrix w – Computed singular values u – Comp...
基于opencv的SVD分解求解变换矩阵
zhoufm260613的博客
07-23 2876
基于opencv的svd分解，求解变换矩阵，通过输入多个点在两个坐标系下的坐标，输出两个坐标系的旋转矩阵和平移向量
SVD分解结果优化
02-27
在OpenCV库中，Mat类提供了内置的SVD函数，可以对任意大小的矩阵进行SVD分解。然而，原生的OpenCV 2.4.9版本在某些情况下，如特征值为0时，可能会导致U矩阵的不完整表示，从而影响到后续的计算和应用。 在SVD分解中...
SVD的C++代码实现
02-09
SVD分解算法的C++代码实现 For an m-by-n matrix A, the singular value decomposition is an m-by-m， orthogonal matrix U, an m-by-n diagonal matrix S, and an n-by-n orthogonal ，matrix V so that A = U*S*V...
CUDA 7.5的SVD奇异值分解
04-22
CUDA 7.5的SVD奇异值分解实现，使用cuSolver实现的，工程是VS2013的工程，需要自己设置VC++库的include、library目录，里面有一点opencv的svd分解比较，不想看的就注释掉对应的部分就行。
opencv学习日记——SVD奇异值分解与reshape
weixin_51326570的博客
03-09 1113
opencv学习日记——SVD奇异值分解与reshape #include<iostream> #include<opencv2/opencv.hpp> using namespace std; using namespace cv; int main() { //定义一个3*3的矩阵a float b[] = {1, 2, 3, 4, 5, 6, 7, 8, 9}; Mat a = Mat(3, 3, CV_32FC1, b); Mat w, u
SVD学习：压缩图像（opencv）
zfjBIT的专栏
06-28 784
奇异值分解在图像处理中有着重要应用。假定一幅图像有n∗n个像素，如果将这n^2个数据一起传送，往往会显得数据量很大。因此我们希望能够改传送另外一些比较少的数据，在接收端利用这些数据重构原图像。 假定对矩阵A进行奇异值分解，便得到A=UΣV^2，其中，奇异值从小到大的顺序排列。如果从中选取k个大奇异值以及这些对应的左右奇异向量逼近原图像，便可以共使用k(2n+1)个数值取代原来的n∗n个图像数据。...
opencv3之SVD
windistance的博客
07-22 9133
cv::SVD::compute static void cv::SVD::compute ( InputArray src, OutputArray w, OutputArray u, OutputArray vt, int flags = 0 ) 参数： //A=u*w*vt src decomposed matrix w calcul...
OpenCV——svd速度优化
coco的博客
08-23 1941
在发现opencv运行SVD::compute()过程太慢，一张299*299的图片Debug要几s，release也要0.5s。 第一次尝试： 将299*299图像分成13份23*299的图像运行svd。release模式下能达到0.25s，将近快一倍。 第二次： 加上UMat，成绩感人。release模式下速度0.02s，又快了10倍。 void main() { Mat ...
opencv mat初始化_14、OpenCV主成分分析、奇异值分解和随机数生成
weixin_39941262的博客
11-26 250
最近太忙了，今天抽空更新一篇。随着OpenCV库的发展，引入新封装的对象变得越来越普遍，这些封装功能过于复杂，无法与单个功能相关联，并且如果作为一组功能实施，会导致库的整体功能空间变得太混乱。因此，新功能通常由相关的新对象类型来表示，它可以执行任何功能。大多数这些都有一个重载的operator()。他们可以设置需要的任何数据或配置，并且他们被要求通过普通成员函数来执行，或者通过重载operator...
opencv SVD算法
zhanggusheng的博客
04-13 4339
参照：https://www.cnblogs.com/pinard/p/6251584.html #include <iostream> #include <opencv2/opencv.hpp> using namespace std; using namespace cv; //SVD算法 int main() { string filename = "C:/...
OpenCV 和 Eigen求解线性方程 Ax=B
douzi949389的博客
02-14 2739
#include <iostream> #include <Eigen/Core> #include <Eigen/Eigen> #include <Eigen/Geometry> #include<opencv2/opencv.hpp> int test_opencv() { printf("\nSolve equation:AX=b\n\n"); cv::Mat A = (cv::Mat_<float>(4,.
基于openCV数字图像与机器视觉（SVD特征值提取，检测硬币、细胞、条码）
伊木子曦的博客
12-10 3329
文章目录用奇异值分解（SVD）进行特征值提取（降维）处理；检测图像中硬币、细胞的个数。图片中的条形码检测与识别 用奇异值分解（SVD）进行特征值提取（降维）处理； 代码 import numpy as np import os from PIL import Image import matplotlib.pyplot as plt import matplotlib as mpl from pprint import pprint def restore1(sigma, u, v, K): #
opencv2.4中SVD分解的几种调用方法
weixin_34342992的博客
12-02 484
原帖地址： http://blog.sina.com.cn/s/blog_6109b5d00101ag7a.html       在摄影测量和计算机视觉中，考虑最优解问题时，经常要用到SVD分解。奇异值分解 (singular value decomposition,SVD) 是一种可靠地正交矩阵分解法，但它比QR分解法要花上近十倍的计算时间。在matlab中，[U,S,V]=svd(A)...
opencv实现奇异值分解图片压缩
最新发布
10-15
在OpenCV中，奇异值分解（SVD，Singular Value Decomposition）常用于图像处理中的降维和压缩。SVD将一个矩阵分解成三个部分：U、Σ和V^T，其中U和V是正交矩阵，Σ是对角矩阵，包含原始数据的主要特征值。对于图像来...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

是原来的你吗

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值