与OpenCV的第十天

第一件事:离散傅里叶变换( DFT )


关于二维离散傅里叶变换的理论部分请移步我在简书的另一篇文章 《图像傅里叶变换(二维离散傅里叶变换)》, 这里只介绍在 OpenCV 中是如何计算并可视化二维离散傅里叶变换的


变换步骤如下:

第一步:扩展矩阵面积至最优面积

通过改变图片矩阵的面积可以增加 DFT 的执行效率,当图片面积为 2、3、5 的倍数时 DFT 执行效率最快,因此为了达到 DFT 的执行效率最快,我们需要对原矩阵进行扩展。

  • 计算需要扩展的行数和列数
    OpenCV 为我们提供了这样一个函数 int getOptimalDFTSize(int vecsize),这个函数传入一个原矩阵的行数或列数,返回最优的行数或列数。
    示例
    int opr = cv::getOptimalDFTSize(srcImg.rows);
    int opc = cv::getOptimalDFTSize(srcImg.cols);
  • 扩展面积至最优
    这一步,我们可以通过 OpenCV 为我们提供的 copyMakeBorder 函数进行扩展。
    示例
cv::Mat paddedImg ;
cv::copyMakeBorder(srcImg, paddedImg, 0, opr - srcImg.rows, 0, opc - srcImg.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));

函数原型

void copyMakeBorder(InputArray src, OutputArray dst,
                                 int top, int bottom, int left, int right,
                                 int borderType, const Scalar& value = Scalar() );
参数 意义
InputArray src 输入矩阵
OutputArray dst 扩展之后的矩阵
int top 顶部扩展行数
int bottom 下部扩展行数
int left 左侧扩展列数
int right 右侧扩展列数
int borderType 扩展的边界类型
BORDER_CONSTANT ( iiiiii [ abcdefgh ] iiiiiii 以指定的 i 值进行扩展 )
BORDER_REPLICATE( aaaaaa [ abcdefgh ] hhhhhhh
BORDER_REFLECT( fedcba [ abcdefgh ] hgfedcb
BORDER_WRAP( cdefgh [ abcdefgh ] abcdefg
BORDER_REFLECT_101( gfedcb [ abcdefgh ] gfedcba
BORDER_TRANSPARENT ( uvwxyz [ absdefgh ] ijklmno
const Scalar& value 指定的填充值
第二步:将单通道扩展至双通道,以接收 DFT 的复数结果

通过 与OpenCV的第九天 中所学到的 merge 函数进行通道扩展。
示例

    cv::Mat complexI;
    cv::Mat planes[] = {cv::Mat_<float>(paddedImg), cv::Mat::zeros(paddedImg.size(), CV_32F)};
    cv::merge(planes, 2, complexI);
第三步:就地(输入作为输出)傅里叶变换

示例

cv::dft(complexI,complexI);

dft 函数详解
该函数的作用是对一维或二维数组进行正向或逆向傅里叶变换。
函数原型:

void dft(InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0);
  • 第一个参数:InputArray src,输入的矩阵
  • 第二个参数:OutputArray dst,最后运算所得结果
  • 第三个参数:int flags 转换标识符,默认值为 0 ,源码及释义:
enum DftFlags {
    /** performs an inverse 1D or 2D transform instead of the default forward
        transform. */
/** 译:
执行一维或二维的逆变换,
而不是默认的正向变换。 */
    DFT_INVERSE        = 1,
    /** scales the result: divide it by the number of array elements. Normally, it is
        combined with DFT_INVERSE. */
/** 译:
通过除以数组元素的数量对输出结果进行缩放,
通常会结合 DFT_INVERSE 使用 */
    DFT_SCALE          = 2,
    /** performs a forward or inverse transform of every individual row of the input
        matrix; this flag enables you to transform multiple vectors simultaneously and can be used to
        decrease the overhead (which is sometimes several times larger than the processing itself) to
        perform 3D and higher-dimensional transformations and so forth.*/
/** 译:
 对输入矩阵的每一行进行正向或逆向变换,
设定这个标识符将允许你同时对多个向量进行变换,
还可以减少三维和多维变换等操作的开销。*/
    DFT_ROWS           = 4,
    /** performs a forward transformation of 1D or 2D real array; the result,
        though being a complex array, has complex-conjugate symmetry (*CCS*, see the function
        description below for details), and such an array can be packed into a real array of the same
        size as input, which is the fastest option and which is what the function does by default;
        however, you may wish to get a full complex array (for simpler spectrum analysis, and so on) -
        pass the flag to enable the function to produce a full-size complex output array. */
/**这里需要详细解释一下,
首先我们知道的是傅里叶变换的结果矩阵是关于原点共轭对称的,
那么两个关于原点对称的点的值,实部相等,虚部相反,
所以我们可以在一个象限中只存储实部,在另一个对称象限中只存储虚部。
这样便可以将原来需要两个通道的存储的复数值压缩进只有一个通道的矩阵进行存储,
这也是该参数的意义所在。
dft 函数的默认执行方式是:
当输入阵列为单通道时,将使用单通道阵列压缩存储输出结果;
当输入阵列为双通道时,将使用双通道阵列存储输出结果;
当指定该参数之后,将使用双通道阵列存储输出结果。
*/
    DFT_COMPLEX_OUTPUT = 16,
    /** performs an inverse transformation of a 1D or 2D complex array; the
        result is normally a complex array of the same size, however, if the input array has
        conjugate-complex symmetry (for example, it is a result of forward transformation with
        DFT_COMPLEX_OUTPUT flag), the output is a real array; while the function itself does not
        check whether the input is symmetrical or not, you can pass the flag and then the function
        will assume the symmetry and produce the real output array (note that when the input is packed
        into a real array and inverse transformation is executed, the function treats the input as a
        packed complex-conjugate symmetrical array, and the output will also be a real array). */
/** 译:
对一维或二维复数阵列进行傅里叶变换的结果通常是一个复数阵列,
对一维或二维的共轭对称阵列(如正向变换生成的矩阵)进行变换时结果是一个实数阵列。
但是 dft 函数并不会判断输入矩阵是否是共轭对称的,
我们可以通过指定这个参数来让 dft 函数假定输入阵列是具有对称性的,来生成实数阵列。
(注意:当输入阵列是一个压缩过的一通道阵列,如 dft 函数的默认正向变换结果,
那么 dft 函数会将其视为具有共轭对称性,从而输出实数阵列)
*/
    DFT_REAL_OUTPUT    = 32,
    /** performs an inverse 1D or 2D transform instead of the default forward transform. */
    DCT_INVERSE        = DFT_INVERSE,
    /** performs a forward or inverse transform of every individual row of the input
        matrix. This flag enables you to transform multiple vectors simultaneously and can be used to
        decrease the overhead (which is sometimes several times larger than the processing itself) to
        perform 3D and higher-dimensional transforms and so forth.*/
    DCT_ROWS           = DFT_ROWS
};
  • 第四个参数:int nonzeroRows,默认值为 0 ,当此参数设置为非零值时,dft 函数假定输入阵列只有第一个非零行有非零元素(正向变换)或输出阵列只有第一个非零行有非零元素(逆向变换)。
    这个方法在计算数组的互关联和卷积时非常有用。
第四步:将双通道复数矩阵转换为单通道幅度矩阵
  • 先分解为单通道实部矩阵和虚部矩阵,通过在《与OpenCV的第九天》 所学到的 split 函数进行通道分离。
    示例
//这里直接使用了第二步中的 planes
cv::split(complexI, planes);
  • 使用 magnitude 函数计算幅值矩阵
    示例
cv::magnitude(planes[0], planes[1], planes[0]);

函数原型

void magnitude(InputArray x, InputArray y, OutputArray magnitude);
  • 第一个参数,矢量的浮点型 x 坐标矩阵
  • 第二个参数,矢量的浮点型 y 坐标矩阵
  • 第三个参数,输出的结果矩阵

以下步骤都是加强视觉效果


第五步:转换幅值至对数级

将计算得到的幅值范围大到难以在屏幕上显示,为了可以通过图像观察到幅值变化,我们需要将高值变为亮点,而低值变为暗点。所以我们需要将幅值转至对数级。
示例

magI += cv::Scalar::all(1);
cv::log(magI, magI);
第六步:重分布矩阵
  • 取行列为偶数
magI = magI(cv::Rect(0, 0,magI.cols & -2,magI.rows & -2));//取小于行列的最大偶数
  • 矩阵重排列,将对角矩阵进行交换。
    int cx = magI.cols/2;
    int cy = magI.rows/2;
    
    cv::Mat q0(magI,cv::Rect(0,0,cx,cy));
    cv::Mat q1(magI,cv::Rect(cx,0,cx,cy));
    cv::Mat q2(magI,cv::Rect(0,cy,cx,cy));
    cv::Mat q3(magI,cv::Rect(cx,cy,cx,cy));
    
    //重排列
    cv::Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
  • 虽然进行了指数变换,但是仍然超出了 0-1 的范围。所以需要用 normalize 函数 进行规范化
    cv::normalize(magI, magI, 0, 1, CV_MINMAX);
第六步:显示幅值图
    cv::imshow("dstImg", magI);
    cv::waitKey();

OpenCV 计算并可视化二维 DFT 结束,以下是完整代码和运行效果图。


完整代码

#include <iostream>
#include <opencv2/opencv.hpp>
int main(int argc, const char * argv[]) {
    //傅里叶变换
    //读取为灰度矩阵
    cv::Mat srcImg = cv::imread("1.jpg",CV_LOAD_IMAGE_GRAYSCALE);
    cv::imshow("srcImg", srcImg);
    cv::Mat paddedImg ;
    //计算最优大小
    int opr = cv::getOptimalDFTSize(srcImg.rows);
    int opc = cv::getOptimalDFTSize(srcImg.cols);
    
    //扩展矩阵面积至最优
    cv::copyMakeBorder(srcImg, paddedImg, 0, opr - srcImg.rows, 0, opc - srcImg.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));
    
    //dft 函数计算结果为复数矩阵,需要将原单通道灰度矩阵扩展至双通道以存储复数,用之前学过的 merge 函数将单通道矩阵与同等大小单通道零矩阵合并产生双通道矩阵。
    cv::Mat complexI;
    cv::Mat planes[] = {cv::Mat_<float>(paddedImg), cv::Mat::zeros(paddedImg.size(), CV_32F)};
    cv::merge(planes, 2, complexI);
    
    //就地傅里叶变换
    cv::dft(complexI,complexI);
    
    //结果是复数,为了方便分析,我们将其转换成幅度矩阵
    //先分解为实部和虚部
    cv::split(complexI, planes);
    //通过magnitude 函数进行幅度值计算
    cv::magnitude(planes[0], planes[1], planes[0]);
    cv::Mat magI = planes[0];
    
    //将计算得到的幅值转换到对数级,转换得到的幅值范围大到难以在屏幕上显示。为了可以通过图像观察幅值变化,我们需要将高值变为亮点,而低值变为暗点。进行对数变换:
    magI += cv::Scalar::all(1);
    cv::log(magI, magI);
    
    //裁剪重分布矩阵,将原点移到图像中心。因为之前学到的傅里叶变换的周期性,所以我们只需要将图像分为四个 1/4 子矩阵,然后将对角矩阵交换位置即可。
    //裁剪,将第一步扩展的矩阵进行裁剪
    magI = magI(cv::Rect(0, 0,magI.cols & -2,magI.rows & -2));
    
    int cx = magI.cols/2;
    int cy = magI.rows/2;
    
    cv::Mat q0(magI,cv::Rect(0,0,cx,cy));
    cv::Mat q1(magI,cv::Rect(cx,0,cx,cy));
    cv::Mat q2(magI,cv::Rect(0,cy,cx,cy));
    cv::Mat q3(magI,cv::Rect(cx,cy,cx,cy));
    
    //重排列
    cv::Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
    
    //虽然进行了指数变换,但是仍然超出了 0-1 的范围。所以需要用 normalize 函数 进行规范化
    cv::normalize(magI, magI, 0, 1, CV_MINMAX);
    
    cv::imshow("dstImg", magI);
    
    cv::waitKey();
    
    return 0;
}
运行效果图【原始灰度图】
运行效果图【可视化之后的幅度图】
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 201,924评论 5 474
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,781评论 2 378
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 148,813评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,264评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,273评论 5 363
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,383评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,800评论 3 393
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,482评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,673评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,497评论 2 318
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,545评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,240评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,802评论 3 304
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,866评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,101评论 1 258
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,673评论 2 348
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,245评论 2 341

推荐阅读更多精彩内容