cv::Mat temp0 = R.t(); cv::Mat temp1 = R * temp0; cv::Mat temp2(960,960,CV_32FC1); temp2 = temp1.inv(); cv::Size s = temp1.size(); std::cout<<s.height<<" "<<s.width<<std::endl; std::cout<<cv::format(temp1,"numpy" ) << std::endl; std::cout<<cv::format(temp2,"numpy" ) << std::endl;
Transpose工作正常,矩阵乘法也是如此.因此,Mat temp1的大小为960×960.但是,当我执行temp2 = temp1.inv()时,我收到temp2中的所有零.我的意思是零是所有960×960细胞.此外,R仅为CV_32FC1类型.所以它可能不是数据类型问题.我无法理解这里的问题.我用Google搜索了这么多.你能帮忙吗?
编辑
我正在复制Mat :: inv()函数的gdb输出下面.我很难搞清楚这一切,但如果有人更熟悉OpenCV,也许它会有所帮助:)
Breakpoint 1,CreateShares::ConstructShares (this=0x80556d0,channel=...,k=2,n=4) at CreateShares.cpp:165 165 temp2 = temp1.inv(); (gdb) step cv::Mat::operator= (this=0xbffff294,e=...) at /usr/include/opencv2/core/mat.hpp:1373 1373 e.op->assign(e,*this); (gdb) 1374 return *this; (gdb) step 1375 } (gdb) step cv::MatExpr::~MatExpr (this=0xbfffef64,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:1167 1167 class CV_EXPORTS MatExpr (gdb) step cv::Mat::~Mat (this=0xbfffefdc,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:295 295 release(); (gdb) step cv::Mat::release (this=0xbfffefdc) at /usr/include/opencv2/core/mat.hpp:381 381 if( refcount && CV_XADD(refcount,-1) == 1 ) (gdb) step 383 data = datastart = dataend = datalimit = 0; (gdb) step 384 size.p[0] = 0; (gdb) step 385 refcount = 0; (gdb) step 386 } (gdb) step cv::Mat::~Mat (this=0xbfffefdc,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:296 296 if( step.p != step.buf ) (gdb) step 298 } (gdb) step cv::Mat::~Mat (this=0xbfffefa4,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:295 295 release(); (gdb) step cv::Mat::release (this=0xbfffefa4) at /usr/include/opencv2/core/mat.hpp:381 381 if( refcount && CV_XADD(refcount,-1) == 1 ) (gdb) step 383 data = datastart = dataend = datalimit = 0; (gdb) step 384 size.p[0] = 0; (gdb) step 385 refcount = 0; (gdb) step 386 } (gdb) step cv::Mat::~Mat (this=0xbfffefa4,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:296 296 if( step.p != step.buf ) (gdb) step 298 } (gdb) step cv::Mat::~Mat (this=0xbfffef6c,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:295 295 release(); (gdb) step cv::Mat::release (this=0xbfffef6c) at /usr/include/opencv2/core/mat.hpp:381 381 if( refcount && CV_XADD(refcount,-1) == 1 ) (gdb) step 383 data = datastart = dataend = datalimit = 0; (gdb) step 384 size.p[0] = 0; (gdb) step 385 refcount = 0; (gdb) step 386 } (gdb) step cv::Mat::~Mat (this=0xbfffef6c,__in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:296 296 if( step.p != step.buf ) (gdb) step 298 } (gdb) step CreateShares::ConstructShares (this=0x80556d0,n=4) at CreateShares.cpp:167 167 cv::Size s = temp1.size(); (gdb) step cv::Mat::MSize::operator() (this=0xbffff284) at /usr/include/opencv2/core/mat.hpp:705 705 return Size(p[1],p[0]); (gdb) step cv::Size_<int>::Size_ (this=0xbffff2f8,_width=960,_height=960) at /usr/include/opencv2/core/operations.hpp:1624 1624 : width(_width),height(_height) {} (gdb) step cv::Mat::MSize::operator() (this=0xbffff284) at /usr/include/opencv2/core/mat.hpp:706 706 } (gdb) step
解决方法
来自维基百科:
A square matrix that is not invertible is called singular or
degenerate. A square matrix is singular if and only if its determinant
is 0.
你可以像这样显示行列式……
std::cout<<"determinant(temp1)="<<cv::determinant(temp1)<<"\n";
从the documentation for Mat::inv()开始,有三种方法可供选择:
> DECOMP_LU(默认值)是LU分解.矩阵必须是非单数的.
> DECOMP_CHOLESKY只是对称正定义矩阵的Cholesky分解.这种类型比大矩阵上的LU快两倍.
> DECOMP_SVD是SVD分解.如果矩阵是奇异的或甚至是非正方形,则计算伪反转.
从the documentation for invert()开始,这可能是由Mat :: inv()在内部使用的:
In the case of DECOMP_LU method,the function returns the src
determinant ( src must be square). If it is 0,the matrix is not
inverted and dst is filled with zeros.
这与您所看到的结果一致.
关于数学的注释
我不是数学家,但我觉得反转矩阵可能是一个混乱的事情 – 如果你的矩阵非常大,更是如此.实际上,原则上存在这些反转可能是正确的,但实际上不可能以任何精度计算.在使用您的代码进行一些实验时,我发现在很多情况下,我会得到不完全为零的行列式,但非常接近于零 – 可能表明数值精度可能是限制因素.我尝试使用64位值而不是32来指定矩阵,并且得到了不同,但不一定是更好的答案.
根据计算temp1矩阵的方式,识别它总是对称的可能是有用的. DECOMP_CHOLESKY方法专门设计用于对称positive definite矩阵,因此使用它可能会提供一些优势.
在实验上,我发现归一化(由@cedrou建议)使得反函数更有可能返回非零矩阵(使用DECOMP_LU但不使用DECOMP_CHOLESKY).但是,根据我对如何初始化R矩阵的猜测,得到的矩阵似乎永远不会满足逆矩阵的定义:A * inverse(A)= Identity.但是你并不一定关心这一点 – 这也许是SVD方法计算伪逆的原因.
最后,似乎关于为什么反演失败的这个更深层次的问题可能是数学问题而不是编程问题.基于此我在数学网站上做了一些搜索,结果发现有人已经问过如何做这件事:https://math.stackexchange.com/questions/182662
关于调试的说明
根据您的调试跟踪,我倾向于认为您感兴趣的部分已编译为不可跟踪的库,并在您运行步骤时跳过.换句话说,在第一步之后的那个神秘的空白行代表它实际运行inv()函数的部分.之后,它将结果分配给temp2并破坏临时对象.因此,您的调试跟踪不会告诉我们有关inv()内部发生的事情.
165 temp2 = temp1.inv(); (gdb) step cv::Mat::operator= (this=0xbffff294,*this);
我自己运行了一个调试器,并能够通过内部调用遍历invert()并观察它根据矩阵的内部分析(确定它不可逆)决定失败 – 并因此返回填充的矩阵使用零,匹配您报告的内容.
如果您有兴趣查看源代码,则可以在cxlapack.cpp中定义invert()函数.