【OpenCV3.3】特征值、奇异值分解与图像矩阵重构

      在图像处理方面,矩阵分解被广泛用于降维(压缩)、去噪、特征提取、数字水印等,是十分重要的数学工具,其中特征分解(谱分解)和奇异值分解是两种常用方法,本文简单介绍如何在OpenCV中使用它们对图像进行分解,然后重新构造图像。

      本文不会阐述两种分解的数学背景知识,但是为了方便读者唤醒记忆,会先贴出(部分)数学定义,详细的介绍和证明建议阅读矩阵理论相关书籍或者参考资料。

特征值分解

矩阵对角化定理(Matrix diagonalization theorem):对于 N×N 方阵 A ,如果它有 N 个线性无关的特征向量,那么存在一个特征分解

A=QΛQ−1

其中, Q N×N 方阵,且其第 i 列为A 的特征向量qiΛ是对角矩阵,其对角线上的元素为对应的特征值,即 Λiii

更进一步,如果方阵 A 是对称方阵,可得 Q 的每一列都是 A 的互相正交且归一化(单位长度)的特征向量,即 Q−1=QT,此即 对称对角化定理(Symmetric diagonalization theorem)。

奇异值分解

定义    设的特征值为

则称 A 的奇异值;当 A为零矩阵时它的奇异值都是0。

定理    设,则存在m阶酉矩阵 U n 阶酉矩阵V,使得

其中,而为矩阵A 的全部非零奇异值。改写上式为

称该式为矩阵 A 奇异值分解

      由上述定义(结合矩阵理论)可以看到,特征分解相比奇异值分解其应用条件要苛刻得多,特征分解要求矩阵与 对角矩阵(可对角化) 或 Jordan标准形矩阵 相似,这在图像上一般是不满足的,而奇异值分解适用于任何矩阵,应用范围也就更加广泛,但是奇异值分解本身时间复杂度相当高,随着矩阵规模的增长,计算复杂度将呈立方增长,所以很多情况下直接套用而不考虑优化是达不到性能指标的。

      落实到OpenCV中,提供了两个直接计算的方法,其中cv::eigen用于计算对称矩阵的特征值和特征向量,cv::SVD::compute用于进行奇异值分解,两者在目前版本使用的都是雅可比(Jacobi)法,下面是具体的C++代码例子,完成了两种分解与矩阵重构:

#include "stdafx.h"
#define CV_SHOW(a) cv::namedWindow(#a, cv::WINDOW_NORMAL); cv::imshow(#a, a); cv::resizeWindow(#a, a.cols * 3, a.rows * 3)

//-------------------------------------------------------------------------

INT WINAPI WinMain(_In_ HINSTANCE, _In_opt_ HINSTANCE, _In_ LPSTR, _In_ INT)
{
	cv::ipp::setUseIPP(true);

	// test.jpg是任意一副图片
	cv::Mat image = cv::imread("test.jpg", cv::IMREAD_REDUCED_GRAYSCALE_2);
	CV_SHOW(image);

	// 转换为OpenCV接受的计算类型,CV_32F或CV_64F
	cv::Mat srcimage;
	image.convertTo(srcimage, CV_32FC1);

	// 因为非方阵没有特征值,故取一个正方形, 注意未必是对称矩阵,并且 cv::eigen 计算时只会取上三角矩阵,故还原出的图片只是一半+对称
	cv::Mat square_word(srcimage, cv::Rect(0, 0, MIN(srcimage.cols, srcimage.rows), MIN(srcimage.cols, srcimage.rows)));
	cv::Mat vals, vecs;
	cv::eigen(square_word, vals, vecs);
	// 上述函数返回的特征值是列向量,转为对角矩阵,方便后面直接应用矩阵乘法
	cv::Mat diagonal_vals(vals.rows, vals.rows, vals.type(), cv::Scalar(0));
	float *__restrict vals_data = vals.ptr<float>(0);
	for (int i = 0; i < vals.rows; ++i) { // 尝试减少i的迭代次数,可以看到恢复后的图像与原图差距越来越大
		diagonal_vals.ptr<float>(i)[i] = vals_data[i];
	}
	// 根据公式,我们需要计算特征向量构成的矩阵的逆矩阵;由于属于不同特征值的特征向量线性无关,故而vecs必定可逆;
	// 注意每个特征向量以行形式存在vecs里,需要先转为列形式
	cv::transpose(vecs, vecs);
	cv::Mat vecs_inverted;
	cv::invert(vecs, vecs_inverted);
	// 恢复原图(的一半)
	cv::Mat e = vecs * diagonal_vals * vecs_inverted;
	// cv::imshow时如果输入图像是32位浮点类型,每个像素值会乘以255,这不是我们想要的
	e.convertTo(e, CV_8UC1);
	CV_SHOW(e);

	// 奇异值分解适用于任何矩阵,直接对srcimage进行计算
	cv::Mat w, U, Vt;
	cv::SVD::compute(srcimage, w, U, Vt);
	// 同样的,存放奇异值的W是列向量,转为对象矩阵
	cv::Mat W(w.rows, w.rows, w.type(), cv::Scalar(0));
	float *__restrict w_data = w.ptr<float>(0);
	for (int i = 0; i < w.rows; ++i) { // 奇异值已经按递减排序,一般取前面几个大的值就足以还原出基本原图,即有损压缩,这一过程也会将一些噪声压缩掉
		W.ptr<float>(i)[i] = w_data[i];
	}
	// 恢复原图
	cv::Mat svd = U * W * Vt;
	// 显示
	svd.convertTo(svd, CV_8UC1);
	CV_SHOW(svd);

	return cv::waitKey();
}

        参考资料

  1. 黄廷祝,钟守铭,李正良. 矩阵理论[M].北京:高等教育出版社,2003.32-37.
  2. 程云鹏,张凯院,徐仲. 矩阵论[M].西北工业大学出版社,2006.9.1
  3. 特征值分解和主成份分析
阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页