FFmpeg提取与合并命令使用小结

前言

  前面有一段时间需要对视频做一些简单的处理,所以就去查了一下一些常用的视频处理工具,因为不想再另外装个什么软件,所以就决定直接采用 FFmpeg 了,毕竟其它的一些视频处理软件也极有可能只是对 FFmpeg 进行一些图形化界面的封装而已。

前言

  前面有一段时间需要对视频做一些简单的处理,所以就去查了一下一些常用的视频处理工具,因为不想再另外装个什么软件,所以就决定直接采用 FFmpeg 了,毕竟其它的一些视频处理软件也极有可能只是对 FFmpeg 进行一些图形化界面的封装而已。

  而在查找 FFmpeg 相关资料的时候,也同时发现了 Libav , 对于 FFmpeg 和 Libav Shaun 只想说,开源界的战争总是这么莫名其妙,有各种各种一些奇怪的原因(或许大佬们都是各有各的脾性吧 ╮(╯▽╰)╭),像两个小飞机的战争也是如此,不过神仙打架,凡人享福也好 ~\(≧▽≦)/~。神仙们的想法我等凡人是无法揣测了,不过这种事还是少出点为好,不然下次变为小鬼遭殃就不好了。

  使用 FFmpeg 而不是 Libav 是因为 FFmpeg 的相关资料相对来说要多很多;而且Libav 更新的特性,FFmpeg 全都支持;更重要的是 FFmpeg 是更新特别频繁,差不多是日更了,几乎每天都会发布一个稳定版;而且从两者的下载页面来看, FFmpeg 好像更会照顾跨平台用户一点。下文主要是对 FFmpeg 的部分命令进行记录总结。

Shaun 的系统环境为 Win10 1607,测试的 FFmpeg 版本为 ffmpeg-20180429-19c3df0-win64-static

  在 Win10 下使用 FFmpeg 前首先需要配置系统环境变量,Windows 下 FFmpeg 的安装是下载完 Windows 版编译好的 FFmpeg 压缩包后直接解压即可,将解压后的 FFmpeg 文件夹中的 bin 目录添加到系统环境变量 Path 中,不然系统会无法找到 ffmpeg 命令。

提取命令

  FFmpeg 的提取命令的参数很多,因此可以使用不同参数排列组合达到不同的提取效果,具体如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# (1)快速提取video.mp4从第一分钟开始持续两分钟的视频,即到第三分钟,并将提取结果输出为cut.mp4
ffmpeg -ss 00:01:00 -i video.mp4 -to 00:02:00 -c copy cut.mp4

# (2)快速提取video.mp4从第一分钟开始持续两分钟的视频,即到第三分钟,并将提取结果输出为cut.mp4
ffmpeg -ss 00:01:00 -i video.mp4 -t 00:02:00 -c copy cut.mp4

# (3)快速提取video.mp4从第一分钟开始到第二分钟的视频,并将提取结果输出为cut.mp4
ffmpeg -ss 00:01:00 -i video.mp4 -to 00:02:00 -c copy -copyts cut.mp4

# (4)精确提取video.mp4从第一分钟开始到第二分钟的视频,并将提取结果输出为cut.mp4
ffmpeg -i video.mp4 -ss 00:01:00 -to 00:02:00 -c copy cut.mp4

# (5)精确提取video.mp4从第一分钟开始持续两分钟的视频,即到第三分钟,并将提取结果输出为cut.mp4
ffmpeg -i video.mp4 -ss 00:01:00 -t 00:02:00 -acodec copy -vcodec copy cut.mp4

# (6)快速提取video.mp4从第三分钟开始持续60秒的视频,即到第四分钟,并将提取结果输出为cut.mp4
ffmpeg -ss 00:03:00 -i video.mp4 -t 60 -c copy -avoid_negative_ts 1 cut.mp4

以下为各参数的含义:

  • -i:用法为 -i INPUT_VIDEO,代表输入的视频,该视频应为 MPEG 编码(h.264, MPEG4/divx/xvid, MPEG2; MP2, MP3, AAC) ;
  • -ss:用法为 -ss START_TIME,代表提取的开始时间,时间格式有两种写法:1、纯数字格式,以秒为单位(eg: -ss 90,代表从第90秒开始提取);2、时:分:秒 格式(eg: -ss 00:01:30,代表从 0 时 1 分 30 秒开始提取);
  • -to:用法为 -to STOP_TIME,代表提取的结束时间,时间格式同样有两种写法:1、纯数字格式,以秒为单位(eg: -to 180,代表第180秒结束提取);2、时:分:秒 格式(eg: -to 00:03:00,代表 0 时 3 分 00 秒结束提取);
  • -t:用法为 -t DURATION_TIME,代表提取的持续时间,时间格式同样有两种写法:1、纯数字格式,以秒为单位(eg: -to 180,代表提取持续180秒);2、时:分:秒 格式(eg: -to 00:03:00,代表提取持续 0 时 3 分 00 秒);
  • -c:用法为 -c CODEC,代表音视频编码格式(若 CODEC 为 copy 则代表输出视频的音视频编码格式与原视频一样),其实 -c-codec 的缩写,其中包含音频编码参数 -acodec 和 视频编码参数 -vcodec
  • -copyts:保持原有时间戳,即使当 -ss-i 之前时,仍使 -t-to 保持原有效果,不会被同化;
  • -avoid_negative_ts:当开启该参数时,提取视频会找到首尾的相邻关键帧(这样会造成提取时间不精确),补全视频,当 -ss-i 之前时该参数默认开启;
  • -accurate_seek:使提取时间更精确,在转码时该参数默认启用。
  • 更多参数可参考:Format-Options

※注: 1、-t-to 不能同时使用,若同时使用,将以 -t 为准;  2、当 -ss-i 之前时,可以实现快速提取,但提取的时间点不精确,此时-to-t 的效果一样,都表现为 -t ,此时可以加上 -copyts 参数使两者效果不一样;  3、当 -ss-i 之后时,提取的时间点比较精确,但提取速度比较慢。

  以上命令(1)、(2)是一样的提取结果,命令(3)的 -t-to 会产生不一样的结果,(6)的 -t-to 会产生一样的结果。因此若只考虑速度,可以使用命令(3),若需使提取时间精确,则需使用命令(4)和(5)

PS:当然还有一种更精确的方式,通过命令 ffmpeg -i input.mp4 -strict -2 -qscale 0 -intra output.mp4 将输入视频由原来的帧间编码转换为帧内编码,使每一帧都成为关键帧,如此转换之后再进行提取,可使提取时间十分精确,但该转换方式会造成视频文件空间成倍增大。

合并命令

  FFmpeg 的合并命令大概有 4 种,可参考: FFMpeg无损合并视频的多种方法 。这里主要介绍两种:

方法一:使用 FFmpeg concat 分离器(推荐)

  该方法使用命令为:ffmpeg -f concat -i filelist.txt -c copy output.mp4,其中 filelist.txt 文件最好和待合并的视频在同一个文件夹中,文件中内容就是待合并的视频的描述,具体内容如下:

1
2
file 'input1.mp4'
file 'input2.mp4'

1
2
file ./input1.mp4
file ./input2.mp4

※注: 待合并的视频文件名最好由英文、数字、连接符(-)或下划线(_)组成,且中间最好不要有空格,不然可能会因为文件名而在合并时出现一些奇怪的错误

方法二:利用中间格式合并

  所谓利用中间格式进行合并是因为可能有些待合并的视频编码不一致,这时采用方法一可能无法进行合并,这时需要对待合并的视频的进行统一转码,都转成同一个编码格式。若要采用这种方法,建议都转成 mpg 格式,因为 mpg 格式可以直接 cat 命令进行合并,具体如下:

1
2
3
4
5
6
7
8
# 将input1.avi转成intermediate1.mpg,输出视频质量为1
ffmpeg -i input1.avi -qscale:v 1 intermediate1.mpg
# 将input2.mp4转成intermediate2.mpg,输出视频质量为1
ffmpeg -i input2.mp4 -qscale:v 1 intermediate2.mpg
# 合并intermediate1.mpg和intermediate2.mpg,输出合并结果intermediate_all.mpg
cat intermediate1.mpg intermediate2.mpg > intermediate_all.mpg
# 将intermediate_all.mpg转成output.avi,输出视频质量为2
ffmpeg -i intermediate_all.mpg -qscale:v 2 output.avi

  其中参数 -qscale 是指使用固定的视频量化标度 ,取值范围为 0.01 ~ 255 ,越小代表质量越好,一般推荐取值为 2 ~ 5,实际使用不能超过 50, -qscale:v 代表设置视频输出质量。

PS: 方法一为无损合并,方法二为有损合并

后记

  以后有用到新的命令再继续记录吧 ↖(^ω^)↗ 。纯命令行不太方便的话,可以试试 Avidemux,全平台通用。

参考资料

[1] ffmpeg Documentationhttps://ffmpeg.org/documentation.html

[2] FFmpeg wiki: Seekinghttp://trac.ffmpeg.org/wiki

[3] FFmpeg:视频转码、剪切、合并、播放速调整

[4] 通过 ffmpeg 无损剪切/拼接视频

[5] 使用ffmpeg合并视频文件的三种方法

解决VSCode使用Cmder作为默认终端问题

前言

  Shaun 最近想换换新口味,想尝试用 Cmder 作为 VSCode 下的默认终端,不想再继续使用 git-bash 了,因为 git-bash 有时会出现一些乱码问题。但是在用 VSCode 集成 Cmder 时出现了几个问题。

前言

  Shaun 最近想换换新口味,想尝试用 Cmder 作为 VSCode 下的默认终端,不想再继续使用 git-bash 了,因为 git-bash 有时会出现一些乱码问题。但是在用 VSCode 集成 Cmder 时出现了几个问题。

问题篇

  如果在 VSCode 用户设置文件中直接添加 Cmder.exe 及其路径,那么在使用 VSCode 终端时会重新打开一个 Cmder 窗口,而不是直接显示在 VSCode 的「终端」里。Shaun 想要的是 Cmder 就是 VSCode 的终端,就在 VSCode 里,就和原有cmd 终端一样,在 VSCode 下的终端可以直接输入相关命令,而不是另外弹出一个命令行窗口。

解决方案篇

方法一

将 Cmder 放进一个文件夹中,文件夹名不带空格,比如 Shaun 所有的绿色软件全部放在 D:\ProgramFiles 下,在 VSCode 用户设置文件中添加:

1
2
3
4
"terminal.integrated.shell.windows": "C:\\Windows\\Sysnative\\cmd.exe",
"terminal.integrated.shellArgs.windows": [
"/k D:\\ProgramFiles\\Cmder\\vendor\\init.bat"
],

方法二

一部分用户可能有点强迫症 ๑乛◡乛๑,硬是要把绿色软件还放入系统盘中的 Program Files 文件夹里,这样在 VSCode 里配置 Cmder 作为默认终端时就会出现问题。主要是因为 Program Files 文件夹名中有空格(这里吐槽一下带空格的文件名真鸡儿坑爹,命令行中根本无法访问,这应该是巨硬的历史遗留问题了, 特立独行的支持带空格的路径名,想显摆一下自己,但以目前的情况看,这种支持简直无力吐槽,造成了一堆问题,和 Windows 路径名中的 \ 有的一拼,都是逼死现代程序员的设计 _(´ཀ`」 ∠)_)。好了,吐槽的话也就说到这里了,如果配置路径里无法避免 Program Files 文件夹,这里有三种解决方案:

  1. Windows 在支持带空格的长文件名的同时,也会分配一个短名称,可以称为该文件夹的别名,通过这个别名就可以在命令行中访问该文件夹,获取这个别名的方法有:在命令行中输入 dir /X ,即可在文件夹名之前的一列的看到该别名,若没有别名则为空白,若要添加别名则需要加入 /N 参数。Shaun 这里显示别名如下:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    $ dir /X
    驱动器 C 中的卷是 System
    卷的序列号是 XXXX-XXXX

    C:\ 的目录

    .........

    2018/05/16 19:21 <DIR> PROGRA~1 Program Files
    2018/05/17 15:14 <DIR> PROGRA~2 Program Files (x86)

    ..........

    由上面可知 Program Files 文件夹的别名为 PROGRA~1,而 Program Files (x86) 文件夹的别名为 PROGRA~2,在配置路径时只需要用别名替换相应的文件夹名即可,如下:

    1
    2
    3
    4
    5
    6
    {
    "terminal.integrated.shell.windows": "C:\\Windows\\Sysnative\\cmd.exe",
    "terminal.integrated.shellArgs.windows": [
    "/k C:\\PROGRA~1\\Cmder\\vendor\\init.bat"
    ]
    }
  2. 通过转义符添加 " " 使 Program Files 作为一个整体,如下:

    1
    2
    "terminal.integrated.shell.windows": "cmd.exe",
    "terminal.integrated.shellArgs.windows": ["/k", "C:\\\"Program Files\"\\Cmder\\vendor\\init.bat"],
  3. 直接在系统环境变量中新建一个变量,将 Cmder 的根目录加进去,如下:

    变量名变量值
    CMDER_ROOTC:\Program Files

    再在 VSCode 用户设置文件中添加:

    1
    2
    3
    4
    "terminal.integrated.shell.windows": "C:\\Windows\\system32\\cmd.exe",
    "terminal.integrated.shellArgs.windows": [
    "/k %CMDER_ROOT%\\vendor\\init.bat"
    ]

通过以上几种方法(推荐直接使用方法一),就能成功在 VSCode 中集成 Cmder,可以直接在 VSCode 的终端里享受 Cmder 了。

后记

  没有什么好说的了,该吐槽也已经吐槽完了。方法二 Shaun 试过,只是感觉应该可行,有问题的小伙伴可以在下方留言(这是不可能的,这辈子都不会有人留言的,反正也没人会看到的 ╮(╯▽╰)╭)。以上解决方案全部来自网络,Shaun 只是做个总结记录,万一以后碰到类似问题至少有多种方案可以尝试。

参考资料

[1] How to use Cmder in Visual Studio Code?

[2] Setting Cmder.exe as integrated shell still opens in separate window

[3] IntelliJ idea webstrom Visual Studio Code vscode 设置cmder为默认终端 Terminal

图割算法之NCuts浅见

前言

  图割算法主要有两个发展方向,一个是 Shaun 上次写过的 Graph Cuts ,通过边界项和区域项(有的可能还会添加一个约束项,或者叫标签项)之和建立能量方程,并利用一定程度的交互式构建一个 S/T 图,最后采用最大流/最小割(Max-flow/Min-cut )算法寻找 S/T 图的最小“割”,从而对图像进行分割;而另一个方向就是 NCuts(Normalized Cuts),中文一般叫规范割、标准割等等,该方法主要出自:Shi J, Malik J. Normalized cuts and image segmentation[J]. IEEE Transactions on pattern analysis and machine intelligence, 2000, 22(8): 888-905.

前言

  图割算法主要有两个发展方向,一个是 Shaun 上次写过的 Graph Cuts ,通过边界项和区域项(有的可能还会添加一个约束项,或者叫标签项)之和建立能量方程,并利用一定程度的交互式构建一个 S/T 图,最后采用最大流/最小割(Max-flow/Min-cut )算法寻找 S/T 图的最小“割”,从而对图像进行分割;而另一个方向就是 NCuts(Normalized Cuts),中文一般叫规范割、标准割等等,该方法主要出自:Shi J, Malik J. Normalized cuts and image segmentation[J]. IEEE Transactions on pattern analysis and machine intelligence, 2000, 22(8): 888-905.

提到 NCuts,就不得不提「谱聚类」,因为 NCuts 可以说是谱聚类的一个应用。

谱聚类篇

谱聚类(Spectral Clustering, SC)是一种基于图论的聚类方法——将无向带权图划分为两个或两个以上的最优子图(sub-Graph),使子图内部尽量相似,而子图间距离尽量远,以达到常见的聚类的目的。

在了解谱聚类之前需要先了解两个概念:拉普拉斯矩阵(Laplace Matrix)和瑞利熵(Rayleigh quotient)。

拉普拉斯矩阵

拉普拉斯矩阵的定义如下:

\(W\) 为无向带权图 \(G=(V,E)\) 的邻接矩阵,\(D\) 为无向带权图 \(G\) 的度矩阵,(所谓的度矩阵即为图中各顶点邻接的所有边权值之和构成的矩阵,是一个对角矩阵,若设顶点 \(i\) 和顶点 \(j\) 之间的权值为 \(w(i,j)\),则度矩阵对角线上的元素 \(d_{i,i}=\sum \limits_{j=1}^n w(i,j)\) ),则拉普拉斯矩阵 \(L=D-W\)。如下图:

若图的邻接矩阵:\(\mathbf{W}=\begin{pmatrix} 0 & 0.1 & 0 & 0 & 0.2 & 0 \\ 0.1 & 0 & 0 & 0.5 & 0.1 & 0 \\ 0 & 0 & 0 & 0.3 & 0 & 0.4 \\ 0 & 0.5 & 0.3 & 0 & 0 & 0.1 \\ 0.2 & 0.1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0.4 & 0.1 & 0 & 0 \end{pmatrix}\) ,其中图中不连通的顶点之间的权值为 0;

则对应的度矩阵为:\(\mathbf{D}=\begin{pmatrix} 0.3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0.7 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0.7 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.9 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0.3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0.5 \end{pmatrix}\)

则拉普拉斯矩阵 \(L\) 为:\(\mathbf{L}=\mathbf{D}-\mathbf{W}=\begin{pmatrix} 0.3 & -0.1 & 0 & 0 & -0.2 & 0 \\ -0.1 & 0.7 & 0 & -0.5 & -0.1 & 0 \\ 0 & 0 & 0.7 & -0.3 & 0 & -0.4 \\ 0 & -0.5 & -0.3 & 0.9 & 0 & -0.1 \\ -0.2 & -0.1 & 0 & 0 & 0.3 & 0 \\ 0 & 0 & -0.4 & -0.1 & 0 & 0.5 \end{pmatrix}\)

以上三个矩阵都为实对称矩阵。对于拉普拉斯矩阵,由 $L* [1,1,1,1,1,1]^T =0* [1,1,1,1,1,1]^T= $ 可知其一定存在一个特征值为 0,对应的特征向量为 \([1,1,1,1,1,1]^T\)

拉普拉斯矩阵是半正定的,且对应的 n 个实数特征值都大于等于 0,即 \(0=λ_1≤λ_2≤\cdots≤λ_n\), 且最小的特征值为 0 。

至于拉普拉斯矩阵的更多性质,可参考:Laplacian matrix

瑞利熵

  瑞利熵的公式形如:\(R(M,x)=\frac{x^*Mx}{x^*x}\) ,其中 M 是厄米特矩阵(Hermitian matrix),满足 \(M=M^H\) ,即矩阵 \(M\) 与其共轭转置矩阵相等,\(x\) 为非 0 向量,\(x^*\) 是指 \(x\) 向量的共轭转置向量。对于实数而言,厄米特矩阵就是对称阵 \(M=M^T\)\(x^*\) 就是 \(x\) 向量的转置向量 \(x'\)\(x^T\)

  若设 \(x^*x=c\),其中 \(c\) 为一个常数,则 \(R(M,x)=\frac{x^*Mx}{c}\) ,则其极值问题可转化为:\(\min R(M,x) = \min x^*Mx \qquad s.t.\ x^*x=c\) ,可利用拉格朗日乘数法求取 \(x^*Mx\) 的极值。 具体求法如下: \[ \mathcal{L}(x) = x^T M x -\lambda \left (x^Tx - c \right) \\ \begin{align} &\frac{d\mathcal{L}(x)}{dx} = 0 \\ &\Rightarrow x^TM - \lambda x^T = 0 \\ &\Rightarrow Mx - \lambda x = 0 \\ &\Rightarrow M x = \lambda x \end{align} \\ \therefore R(M,x) = \frac{x^T M x}{x^T x} = \lambda \frac{x^Tx}{x^T x} = \lambda. \]

瑞利熵有如下几个性质:

  • R 的最大值就是 M 最大特征值,R 的最小值就是 M 最小特征值 ;
  • x 的解,就是 M 对应的特征向量。

至于瑞利熵的更多性质,可参考:Rayleigh quotient

  以上只是一般意义上的普通瑞利熵,还有一个广义瑞利熵(generalized Rayleigh quotient ),其定义为:\(R(M,D;x)=\frac{x^*Mx}{x^*Dx}\) ,即在分母上添加一个 D 矩阵相乘,其中 D 为 Hermite 正定矩阵 ,满足 \(D=D^*\),普通瑞利熵是广义瑞利熵中 D 矩阵为单位矩阵时的情况。广义瑞利熵可以通过以下变换转换为普通瑞利熵: \[ \begin{equation} \begin{aligned} R(M,D;x) &= \frac{x^*Mx}{x^*Dx} \xlongequal{x=D^{-\frac{1}{2}}y} \frac{(D^{-\frac{1}{2}}y)^*M(D^{-\frac{1}{2}}y)}{(D^{-\frac{1}{2}}y)^*D(D^{-\frac{1}{2}}y)} \\ &= \frac{y^*(D^{-\frac{1}{2}})^*M(D^{-\frac{1}{2}}y)}{y^*(D^{-\frac{1}{2}})^*D(D^{-\frac{1}{2}}y)}=\frac{y^*D^{-\frac{1}{2}}MD^{-\frac{1}{2}}y}{y^*D^{-\frac{1}{2}}DD^{-\frac{1}{2}}y} \\ &= \frac{y^*(D^{-\frac{1}{2}}MD^{-\frac{1}{2}})y}{y^*y} \end{aligned} \end{equation} \]

只需要求出矩阵 \(D^{-\frac{1}{2}}MD^{-\frac{1}{2}}\) 的特征值和特征向量即可。


好了,这两个概念介绍完了,就可以真正介绍谱聚类了。

  上面说了,谱聚类是一种基于图论的聚类方法,既然有图,则必有相应的邻接矩阵。设图 G 的顶点被聚类成两类,即将图 G 分割为子图 A 和 B,则所要断开的边的权值之和为代价函数(也叫损失函数),类似于“Graph Cuts”中的能量方程。割边 Cut(A,B) 的具体表示为:\(Cut(A,B)=\sum \limits_{i \in A,j \in B} w(i,j)\)

  设图 \(G\) 中共有 \(n\) 个顶点,则需要构建一个 \(n \times n\) 的邻接矩阵 \(W\),其相应的度矩阵为 \(D\),对应的拉普拉斯矩阵为 \(L=D-W\),令 \(p_i = \begin{cases} l_1 & i \in A \\ l_2 & i \in B\end{cases}\) ,则 \[ Cut(A,B)=\sum \limits_{i \in A,j \in B} w(i,j) = \frac{\sum \limits_{i=1}^n \sum \limits_{j=i}^n w(i,j) (p_i-p_j)^2}{(l_1-l_2)^2}= \frac{\sum \limits_{i=1}^n \sum \limits_{j=1}^n w(i,j) (p_i-p_j)^2}{2(l_1-l_2)^2} \]

当且仅当 \(i\)\(j\) 不属于同一子图时,\((p_i-p_j)^2/(l_1-l_2)^2=1\),否则 \((p_i-p_j)^2/(l_1-l_2)^2=0\),(为什么采取这种计算方式,是因为在无法直接确定割边时,用这种计算方式能确保全部割边且只有割边会被加入权重之和),至于为什么等式最后还要除以 2 ,是因为等式最后的那种写法每条边会被遍历两次,每条割边权值会被计入两次,所以还要除以2。

而: \[ \begin{equation} \begin{aligned} & \sum \limits_{i=1}^n \sum \limits_{j=1}^n w(i,j) (p_i-p_j)^2 = \sum \limits_{i=1}^n \sum \limits_{j=1}^n w(i,j) (p_i^2-2p_ip_j+p_j^2) \\ &= \sum \limits_{i=1}^n \sum \limits_{j=1}^n w(i,j) (p_i^2) + \sum \limits_{i=1}^n \sum \limits_{j=1}^n w(i,j) (p_j^2) - \sum \limits_{i=1}^n \sum \limits_{j=1}^n w(i,j) (2p_ip_j) \\ &= \sum \limits_{i=1}^n p_i^2 (\sum \limits_{j=1}^n w(i,j)) + \sum \limits_{j=1}^n p_j^2 (\sum \limits_{i=1}^n w(i,j)) - 2\sum \limits_{i=1}^n \sum \limits_{j=1}^n p_iw(i,j) p_j \\ &= p^TDp+p^TDp-2p^TWp = 2p^T(D-W)p = 2p^TLp \end{aligned} \end{equation} \]

则:\(Cut(A,B)=\frac{p^TLp}{(l_1-l_2)^2}\) ,当 \(l_1=1,l_2=-1\) 时, \(\min Cut(A,B) = \min p^TLp \qquad s.t.\ p^Tp=n\) ,n 为图的顶点个数。由上可知,该最小割问题即是一个瑞利熵问题,只需求取 \(L\) 的特征值和特征向量。

  对 \(L\) 的特征值进行从小到大排列(即取最小 k 个特征向量),若要分成 k 类,则需要取前 k 个特征值(除 0 以外)对应的特征向量,并归一化,将每一个特征向量按列排列构成一个 \(n \times k\) 的特征矩阵,对特征矩阵的行向量使用 k-means 或其它聚类算法,将 n 个行向量聚成 k 类,每一行都属于某一类,根据聚类结果为每个顶点分配相应的类标签,从而完成谱聚类。所以利用谱聚类求的解相当于是一个近似解,它将连续的 n 维问题离散化为 k 维问题,p 是一个 n 维的标签向量,标签值为 \(\{1,2,\cdots,k\}\) ,而 \(L\) 的特征向量中的元素并不是离散化的 \(\{1,2,\cdots,k\}\) ,无法直接用 \(L\) 的特征向量作为标签向量,所以需要将其离散化,对特征矩阵的行向量进行聚类,从而近似的生成标签向量;由于只取前 k 个特征向量,所以在对特征矩阵进行聚类时被聚类的 n 条数据只有 k 维,而本来 \(L\) 的特征向量至少有 n 个,即被聚类的 n 条数据至少有 n 维,所以从某种程度上,谱聚类同时也对数据进行了降维处理。

※注: 推荐使用 SVD 代替特征值和特征向量的求取。至于谱聚类的更多性质,可参考:Spectral clustering

  因为有时候简单的全局最小割,可能并不是最优割,所以需要对割做一个归一化,即 Normalized Cuts ,简称 NCuts。

图篇

  同样,Shaun 还从图的构造开始,NCuts 所需要构造的图就是最普通的无向带权图,图的顶点由图像中像素点构成,相邻的像素点之间互相连接构成图的边。边的权值计算公式为: \[ w(i,j) = e^{\dfrac{-\|\mathbf{F}(i)-\mathbf{F}(j)\|_2^2}{\sigma_I}} * \begin{cases} e^{\dfrac{-\|\mathbf{X}(i)-\mathbf{X}(j)\|_2^2}{\sigma_X}} & \text{if } \|\mathbf{X}(i)-\mathbf{X}(j)\|_2 < r \\ 0 & \text{otherwise}. \end{cases} \]   其中 \(\mathbf{X}(·)\) 是指该顶点的空间位置向量,即图像中像素点的坐标; \(\mathbf{F}(·)\) 可指像素点的强度(灰度)值,颜色向量或纹理值;\(\|\mathbf{X}(i)-\mathbf{X}(j)\|_2\) 表示向量的「2-范数」,即欧氏距离。当 \(\mathbf{F}(·)\) 表示强度(灰度)值时,NCuts 分割的是灰度图;当 \(\mathbf{F}(·)\) 表示颜色向量(一般在 HSV 颜色空间,有 \((h,s,v)\) 三个颜色分量)时,NCuts 分割的是彩色图。对于这个权值的计算说人话就是:当两个像素点之间的距离大于一个人为指定的 \(r\) 时(其实也就用这个 \(r\) 判断到底是连接 4 邻域还是连接 8 邻域,一般不会连接 8 邻域之外的像素点),就不连接,此时权值设为 0 ;否则,则计算两像素点颜色向量的 RBF 核的值和两像素点坐标向量的 RBF 核的值,并取其乘积作为权值。

好了,Ncuts 的图的构成大致就是这样。接下来就是它的割法了。

割篇

  NCuts 中定义割 \(Cut(A,B)=\sum \limits_{i \in A,j \in B} w(i,j)\) ,定义 \(assoc(A,V)=\sum\limits_{i \in A,t \in V}w(i,t)\) 为子图 A 内所有像素点连接的所有边权重之和,最终的 NCuts 目标函数如下: \[ Ncut(A,B)= \frac{cut(A,B)}{assoc(A,V)}+\frac{cut(A,B)}{assoc(B,V)} \] 即通过除以两个子图内所有像素点连接的所有边权值之和对割做了个归一化处理。

  令 \(\mathbf{W}\) 为图的邻接矩阵, \(\mathbf{D}\) 为其对应的度矩阵,\(L=D-W\) 为其对应的拉普拉斯矩阵,令 \(S_1 = assoc(A,V) = \sum\limits_{i \in A}D_{ii}\)\(S_2 = assoc(B,V) = \sum\limits_{i \in B}D_{ii}\),则 \(S = S_1+S_2=\sum\limits_{i \in A}D_{ii}+\sum\limits_{i \in B}D_{ii}=\sum D_{ii}= \sum \limits_{i=1}^n\sum \limits_{j=1}^n w(i,j)\) ,则:\(NCut(A,B)=\frac{p^TLp}{(l_1-l_2)^2}*(\frac{1}{S_1}+\frac{1}{S_2})\) ,其中向量 \(p\)\(N=|V|\) 维的标签向量, \(p_i = \begin{cases} l_1 & i \in A \\ l_2 & i \in B\end{cases}\) ,令 \(l_1 = \sqrt{\frac{S_2}{S_1S}},l_2=-\sqrt{\frac{S_1}{S_2S}}\) ,则 \(NCut(A,B)=p^TLp\) ,此时 \(p^TDp=\sum p_i^2D_{ii}=\sum\limits_{i \in A}l_1^2D_{ii}+\sum\limits_{i \in B}l_2^2D_{ii}=l_1^2S_1+l_2^2S_2=1\) ,则可化为 \(NCut(A,B)=\frac{p^TLp}{p^TDp}\) ,即 \(\min NCut(A,B) = \min p^TLp \qquad s.t.\ p^TDp=1\) ,该式即为一个广义瑞利熵,只需要求出矩阵 \(D^{-\frac{1}{2}}LD^{-\frac{1}{2}}\) 的特征值和特征向量,按照谱聚类最后的聚类方式将顶点(像素点)聚成 k 类,从而完成图像的分割。而 NCuts 有个特殊的变换,具体如下: \[ \begin{equation} \begin{aligned} NCut(A,B) &=\frac{p^TLp}{p^TDp}\xlongequal{p=D^{-\frac{1}{2}}x} \frac{x^T(D^{-\frac{1}{2}}LD^{-\frac{1}{2}})x}{x^Tx} \\ &= \frac{x^T(D^{-\frac{1}{2}}(D-W)D^{-\frac{1}{2}})x}{x^Tx} \\ &= \frac{x^T(D^{-\frac{1}{2}}DD^{-\frac{1}{2}})x}{x^Tx} - \frac{x^T(D^{-\frac{1}{2}}WD^{-\frac{1}{2}})x}{x^Tx} \\ &= I - \frac{x^T(D^{-\frac{1}{2}}WD^{-\frac{1}{2}})x}{x^Tx} \end{aligned} \end{equation} \] 使用该变换可得出:\(\min NCut(A,B) = \max x^T(D^{-\frac{1}{2}}WD^{-\frac{1}{2}})x \qquad s.t.\ x^Tx=1\) ,这样可不求出拉普拉斯矩阵,直接求出 \(D^{-\frac{1}{2}}WD^{-\frac{1}{2}}\) 的特征值和特征矩阵即可,而此时在谱聚类最后聚类的时候要将其特征值从大到小排列(即取最大 k 个特征向量),取前 k 个特征值对应的特征向量构造特征矩阵进行聚类。

附录

  Ratio Cuts(RCuts) 和 NCuts 的目标函数差不多,两者间的差异就是:RCuts 是除以子图内顶点的个数;NCuts 是除以子图内所有顶点的边权值之和;最后变换完之后 RCuts 是个普通瑞利熵,而 NCuts 是个广义瑞利熵。具体的 RCuts 割法等以后有机会用到再写吧,毕竟子图顶点个数多不代表所占权重就大,而且在割时一般是基于权重的,所以一般而言 NCuts 的结果要比 RCuts 的结果要好。

总结

  由上文中应该可以看出 NCuts 与 Graph Cuts 需要构造的图是不一样,所使用的能量方程从形式上看也是完全不一样的(虽然本质上都是求“最小割”),其构成的图主要差别在于 Graph Cuts 多了 S 和 T 两个顶点,而这两个顶点由交互式获取,NCuts 不需要交互式,因此也没有这两个顶点,Ncuts 构造的图就相当于只由 Graph Cuts 的普通顶点和边界项 n-links 构成图,因此 Ncuts 没有 Graph Cuts 所谓的区域项 t-links ,也不需要求区域项 t-links 的权值。而这两种方法边界项 n-links 的权值求法都是差不多的,都是用 (高斯)RBF 核函数进行相似性度量,即 \(w_{\{i,j\}} \propto exp\left(-\frac{(I_i-I_j)^2}{2\sigma^2}\right)\)

后记

  后面查阅资料发现图割其实还有好几个方向(井底之蛙了o(╯□╰)o),等有机会把用图割做超像素分割的那篇论文看一下吧 😜 。

  5 月份的时候自己抽空实现了一下该算法,发现效果很糟糕,分割效果很不理想,当然不排除 Shaun 代码有 BUG 的原因,或者是还有一些步骤或 trick 漏掉了 ╮(╯▽╰)╭。而且不管是用 svds 还是 eigs 其分割结果都是随机的,只能保证邻域内像素点大概在一起,但具体谁和谁在一起这就是随机的了,这导致每次运行程序都可能有不同的分割结果 _(´ཀ`」 ∠)_ ,后面发现分割结果随机应该是 k-means 的锅(与 svd 和 eig 无关),原因是因为 k-means 每次初始化 k 个种子点时都是随机的,这自然会导致每次分割结果不完全一致 o(╯□╰)o。

参考资料

[1] 【聚类算法】谱聚类(Spectral Clustering)

[2] 谱聚类算法(Spectral Clustering)

[3] 拉普拉斯矩阵(Laplace Matrix)与瑞利熵(Rayleigh quotient)http://www.cnblogs.com/xingshansi/category/958994.html

[4] Normalized cuts and image segmentation 简单实现

[5] 谱聚类(spectral clustering)原理总结http://www.cnblogs.com/pinard/category/894692.html

图割算法之Graph Cuts浅见

前言

  以后可能需要用图割算法做一些事情,所以就简单阅读了一下图割算法中最基础的一篇论文,Boykov Y Y, Jolly M P. Interactive graph cuts for optimal boundary & region segmentation of objects in ND images. Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on. IEEE, 2001, 1: 105-112.

前言

  以后可能需要用图割算法做一些事情,所以就简单阅读了一下图割算法中最基础的一篇论文,Boykov Y Y, Jolly M P. Interactive graph cuts for optimal boundary & region segmentation of objects in ND images. Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on. IEEE, 2001, 1: 105-112.

图篇

  图割算法(Graph Cuts),顾名思义,是基于图的一种图像分割方法,这里的图既是指数据结构中的图结构也是指数学中的图论,毕竟数据结构中的图结构也来自于数学中图论。既然是图,必然离不开图的构建,论文中图的构建如下,可以简称为 S/T 图:

论文中图结构

和普通的无向图 \(G=(V,E)\) 一样,该图也是由顶点集合 \(V\) 和边集合 \(E\) 构成,不一样的是,该图有两种顶点和两种边:

  1. 第一种顶点是普通顶点。由图像中各像素点组成,每个顶点对应于图像中每个像素点。每两个邻域顶点,对应于图像中每相邻两个像素点(对2维图像来说,是8邻域;对3维图像来说,是26邻域),的连接就是一条边,这种边叫做 n-links (neighborhood links)。
  2. 第二种顶点是两个终端顶点,代表目标(object)的 source terminal(简称 S)和代表背景(background)的 sink terminal(简称 T),每个终端顶点都与所有普通顶点相连,这种相连构成的边叫做 t-links (terminal links)。

  设图像 \(P\) 中每个像素点为 \(p\),则 \(p \in P\),目标终端顶点 \(S\),背景终端顶点 \(T\) ,则构建的图 \(G\) 中顶点集合可表示为:\(V = P \cup \{S,T\}\) ;每个 \(p\) 都有两种 t-links 边,分别为 \(\{p,S\}\)\(\{p,T\}\) ,每个 \(p\) 与其邻域内像素 \(q\) 构成的 n-links 边为 \(\{p,q\}\) ,设 \(N\) 代表 n-links 边集合,即邻域边集合,则 \(\{p,q\} \in N\) ,则构建的图 \(G\) 中边集合可表示为:\(E = N \bigcup\limits_{p \in P} \{\{p,S\},\{p,T\}\}\)

因为这构建的是一个无向带权图,所以还需要设定边的权值,论文中边的权值设定如下:

权值(代价,能量)区域
\(\{p,q\}\)\(B_{\{p,q\}}\)\(\{p,q\} \in N\)
\(\{p,S\}\)\(\lambda \cdot R_p("bkg")\)\(p \in P, p \notin O \cup B\)
\(\{p,S\}\)\(K\)\(p \in O\)
\(\{p,S\}\)0\(p \in B\)
\(\{p,T\}\)\(\lambda \cdot R_p("obj")\)\(p \in P, p \notin O \cup B\)
\(\{p,T\}\)0\(p \in O\)
\(\{p,T\}\)\(K\)\(p \in B\)
  • 其中集合 \(O\) 代表确定为目标的像素点的集合,即在交互时手动标注为目标的像素点;
  • 集合 \(B\) 代表确定为背景的像素点的集合,即在交互时手动标注为背景的像素点;
  • \(B_{\{p,q\}}\) 可以表示 {p,q} 不连续的惩罚因子,其正比于一个指数函数,具体为:\(B_{\{p,q\}} \propto exp\left(-\frac{(I_p-I_q)^2}{2\sigma^2}\right) \cdot \frac{1}{dist(p,q)}\) ,其中 \(I_p\)\(I_q\) 代表像素 \(p\) 和像素 \(q\) 的像素值,\(dist(p,q)\) 代表其两像素点之间的距离。从 \(B_{\{p,q\}}\) 对应的函数中可以看出当两像素点之间的差异越大时,其权值越小;两像素点之间距离越大时,权值越小(这里对于二维图像来说,一般只考虑周围 8 邻域内像素点,因此可以简单认为相邻两像素点之间距离为 1)。
  • \(R_p("bkg")\)\(R_p("obj")\) 分别是指像素 p 分配给背景和前景目标的惩罚因子,该惩罚因子与像素 p 属于前景目标的概率 \(Pr(I|O)\) 和背景的概率 \(Pr(I|B)\) 有关,一般取概率的负对数值,具体如下:\(R_p(“obj”) = −lnPr(I_p|O)\)\(R_p(“bkg”) = −lnPr(I_p|B)\) 。至于概率的计算方法可以用简单的直方图概率模型,因为 \(O\)\(B\) 为手动标注的确定的前景目标和背景的像素点集合,因此可以分别统计其灰度直方图,再对直方图频数进行归一化得到分布概率直方图,将像素 p 与 背景分布概率直方图和前景目标概率分布直方图进行对比,即可得到其属于背景的概率和属于前景目标的概率(当然更好的一种计算前景和背景概率的方法是用高斯混合模型)。
  • \(K\) 表示整个无向带权图中最大的权值,具体计算方法如下:\(K=1+\max\limits_{p \in P} \sum \limits_{q:\{p,q\} \in N} B_{\{p,q\}}\) ,即取图像中像素点邻域边权值之和的最大值再加 1 (对于二维图像而言,邻域边权值之和是指 8 邻域边权值之和,图像中每个像素点都有其邻域边权值之和,取所有像素点中的最大值),至于为什么要令 K 为整个无向带权图中权值最大值?请详见下文。
  • 至于其中的参数 \(\lambda\)\(\sigma\) 论文中没有明确指定,在初始计算时可以将其置为 1,随后再慢慢调整。

至此,整个 Graph Cuts 算法中关于图的构建已经完全确定,接下来就是割(cut)的算法了。

割篇

  割(cut)是构建的无向带权图中边集合 \(E\) 的一个子集 \(C\),即割 \(C \subset E\) ,该 cut 的 代价(cost)可表示为:\(|C|=\sum \limits_{e \in C} W_e\),即该割边集合的所有边权值之和。

如果一个割,其包含的所有边的权值之和最小,那么这个割就称为最小割,也就是图割的结果。在图分割里,是要求两个终端顶点被分离开的,即最小割把图的顶点划分为两个不相交的子集 \(S\)\(T\) ,其中 \(s∈S\)\(t∈T\)\(T=V/S\)。事实上,这两个子集对应于图像的前景像素集和背景像素集,也就相当于完成了图像分割。

  割相关的算法有很多,eg:Max-flow/Min-cut、GrabCut、One cut、 Normalized cut、Ratio cut 等等,但是其最关键的地方在于能量方程(或称 代价函数,损失函数)的优化,能量优化的目的在于最小化能量函数,而最小化能量函数时取得割就是最小割,即图割的结果。论文中能量方程的定义如下: \[ E(A)= \lambda \cdot R(A)+B(A) \] 其中: \[ R(A)= \sum_{p \in P} R_p(A_p) \\ B(A)= \sum_{\{p,q\} \in N} B_{\{p,q\}} \cdot \delta(A_p,A_q) \\ 其中 \delta(A_p,A_q) = \begin{cases} 1 & \text{if } A_p \neq A_q \\ 0 & \text{otherwise}. \end{cases} \]

  • 其中令 \(A=(A_1,\cdots,A_p,\cdots,A_{|P|})\) 为二值向量,每个 \(A_p\) 可赋值为 “obj”(可用 1 表示前景) 或 “bkg”(可用 0 表示背景);
  • \(R(A)\) 表示区域项,即上文中的两种 t-links 边相连的图像像素点权值,代表像素点分配给 “obj” 或 “bkg” 的惩罚项,对应于 \(R_p("obj")\)\(R_p("bkg")\) ,求和后即可得到 \(R(A)\) 。当像素点 p 属于前景目标的概率 \(Pr(I|O)\) 大于背景的概率 \(Pr(I|B)\) 时,由上文 \(R_p(A_p) = −lnPr(I_p|A_p)\) 可知,此时 \(R_p("obj")\) 小于 \(R_p("bkg")\) ,即当像素 p 更有可能属于目标时,将 p 归类为目标就会使能量 \(R(A)\) 小,也因此要取概率的负对数值,由此,如果所有像素点都能正确划分为前景和背景,则总能量会达到最小。当像素点 p 已经人为手动标注为确定的前景目标时,这时其与顶点 \(S\) 相连的边的权值为 K,为整个无向带权图中最大的权值,即能量也最大,此时就不可能被分割,这就是为什么要令 K 为整个无向带权图中权值最大值,而其与顶点 \(T\) 相连的权值为 0 ,能量也为 0,此时一定会被分割。
  • \(B(A)\) 表示边界项,即上文中 n-links 边的权值,代表邻域像素点 \({p,q}\) 不连续的惩罚,当其差异越大时, \(B(A)\) 越趋近于 0,即当邻域像素点差异越大时,由上文 \(B_{\{p,q\}} \propto exp\left(-\frac{(I_p-I_q)^2}{2\sigma^2}\right) \cdot \frac{1}{dist(p,q)}\) 可知其权值越小,能量 \(B(A)\) 也越小,则越应该被分割。
  • 参数 \(\lambda\) 用来表示区域项和边界项的重要性,初始计算时同样可以将其置为 1。
  • 至于 \(\delta(A_p,A_q)\) 的个人理解为:当邻域像素点 \(p,q\) 被分配的二值变量不同时,即一个是前景另一个是背景,这时在进行能量计算时需要考虑连接亮点的边 \(\{p,q\}\) 的权值,否则可以不考虑其权值。

最后,本文割的具体实施是采用基于「最大流/最小割」算法进行的,而「最大流/最小割」算法针对的是有向图,所以需要把本文 S/T 图的每条无向边都转化为一去一回的两条有向边。

附录

  该论文还有一个有意思的地方就是其证明了能量函数 E(A) 的值只与割 \(C\) 有关,最小割与最小化能量函数对应。证明过程如下:

由上文可知,对于图像中每个像素点,割 \(C\) 切断且仅切断一条 t-links 边,即要么切断 \(\{p,S\}\) 边,要么切断 \(\{p,T\}\) 边,毕竟一个像素点不可能同时属于前景和背景,也不可能既不属于前景也不属于背景;不属于同一终端顶点相连的 n-links 边也会被割 \(C\) 切断,因为如果不断开,则整个图还是相连的,而如果属于同一终端顶点相连的 n-links 边不应该被割 \(C\) 切断,否则,该割就不是最小割。

对于任意割 \(C\) ,可定义其对应的图像分割结果向量 \(A(C)\) ,如下: \[ A_p(C) = \begin{cases} "obj" & \text{if } \{p,T\} \in C \\ "bkg" & \text{if } \{p,S\} \in C \end{cases} \] 因此当 C 确定之后,分割结果 A 也就确定了。

利用上文给出的边的权值计算方法结合定义的 A 可得出割的代价(cost): \[ \begin{align} |C| &= \sum_{p \notin O \cup B} \lambda \cdot R_p(A_p(C)) + \sum_{p \in O} \lambda \cdot R_p(A_p(C)) + \sum_{p \in B} \lambda \cdot R_p(A_p(C)) + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \\ &= \sum_{p \notin O \cup B} \lambda \cdot R_p(A_p(C)) +0 + 0 + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \end{align} \] 这里是因为 \(O \cap B = \emptyset\) ,而且如果确定像素点 p 属于前景或背景,则其割边的权值必为 0,相应的代价也为 0,eg:设 p 确定为前景,则必定断开 {p,T} 边,而 {p,T} 边的权值为 0。

又因为:\(E(A) = \lambda \cdot R(A) + B(A)\) ,则: \[ \begin{align} E(A(C)) &= \lambda \cdot R(A(C)) + B(A(C)) \\ &= \sum_{p \in P} \lambda \cdot R_p(A_p(C)) + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \\ &= \sum_{p \notin O \cup B} \lambda \cdot R_p(A_p(C)) + \sum_{p \in O} \lambda \cdot R_p("obj") + \sum_{p \in B} \lambda \cdot R_p("bkg") + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \end{align} \] 所以:\(|C| = E(A(C)) - \sum \limits_{p \in O} \lambda \cdot R_p("obj") - \sum \limits_{p \in B} \lambda \cdot R_p("bkg")\)

又因为对于标注确定的前景和背景,任意割 C 的 \(\sum \limits_{p \in O} \lambda \cdot R_p("obj")\)\(\sum \limits_{p \in B} \lambda \cdot R_p("bkg")\) 都是完全确定不变的,可以令其为 const,则 \(E(A) = |C| + const\) ,当 \(C\) 取最小时,对应的 \(E(A)\) 也最小,即当割最小时,能量函数也最小。

后记

  虽说关于 Graph Cuts 网上已经有了无数篇博客对其进行描述,但其中或多或少感觉还是有些令人困惑的地方,于是写下本文将那些个人感觉有些困惑的地方重新梳理一下,再改变一下原论文的行文思路,个人感觉更好理解一点,不然一上来就是神马能量方程,简直一脸懵逼 (・ัω・ั) 。

参考资料

[1] Interactive Graph Cuts for Optimal Boundary and Region Segmentation of Objects in N-D ImageS

[2] Graph Cuts

[3] 图像分割之(二)Graph Cut(图割)

[4] 图像分割算法——Graph Cuts

[5] Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D Images 阅读笔记

C++中static用法小结

前言

  static 是 C++ 中很常用的一个关键字,它的用法也很多,时常会将其弄混,索性做个小结,以免以后忘记了或者继续弄混 (。・ω・。)。

前言

  static 是 C++ 中很常用的一个关键字,它的用法也很多,时常会将其弄混,索性做个小结,以免以后忘记了或者继续弄混 (。・ω・。)。

预备篇

  首先要了解程序中数据的存储形式,一般而言数据的存储形式有三种:

  • 栈区(stack)—— 由编译器自动分配释放,一般用来存放函数的参数值,局部变量的值等;
  • 堆区(heap)—— 由程序员分配释放,对应于对象的 new 或 malloc 和 delete 或 free,若程序员忘记释放,则在程序完全退出之后由操作系统回收;
  • 静态存储区(static)—— 在编译时由编译器分配,在程序完全退出时由操作系统回收,一般用来存放全局变量和 static 变量。

  一般1声明的变量默认(如果变量类型,eg: int, double, ... 等,前不加 static 或其它关键字)都是 auto 2的,其一般存放在栈区,生存周期就只在包围其的 { } 内,在包围其的 { } 外就无法使用该变量。而 static 存放在静态存储区,其生存周期是全局的,它要等整个程序完全退出时才会销毁,在程序运行过程中,每次调用 static 变量都保持上一次调用结束后的值

类中篇

静态成员变量

  类中的静态成员变量被该类的所有实例共享,也可以不通过类的实例使用,在使用时首先需要对其初始化,也必须对其进行初始化,因为类中的静态成员变量只是声明,而且,类中的静态成员变量和普通静态变量一样是在程序初始化的时候分配的,在程序完全退出时由操作系统回收。具体用法如下:

1
2
3
4
5
6
7
8
class Test
{
private:
static int s_value; // 注意,这里不能初始化!因为其不属于类对象,只属于类作用域,独立于该类的任何实例
};

// 在cpp中或类定义体外必须对它进行定义和初始化,因为在程序编译时首先执行的就是对其初始化并分配内存:
int Test::s_value = 0; // 注意,这里没有static的修饰!

总而言之就是:类中的静态成员变量可以简单理解为一个名为 Test::s_value 的全局变量,被所有该类的实例共用,但独立于该类的任何实例,只属于该类作用域,在类的定义中能且只能被声明,不能在类定义体中进行初始化,必须要在类定义体外被定义和初始化

静态成员函数

  类中的静态成员函数和类中的静态成员变量有点类似,其在实现时不需要再加 static 修饰,同样能被该类的所有实例复用,同样只属于类作用域中的全局函数,同样不需要类的实例即可调用。类中的静态成员函数不能访问类的普通成员变量,只能访问类的静态成员变量(可以参考 C++静态成员函数访问非静态成员的几种方法 中的小 trick 访问普通成员变量,但非特殊情况不建议这么做)。具体用法如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
class Test
{
private:
static void func(int i);

// 静态成员函数调用非静态成员变量方法
static void staticTest(Test *t)
{
t->value += 1;
}
private:
int value;
};

// 在cpp中可以不通过类的实例进行调用:
void Test::func(int);

总而言之就是:类中的静态成员函数可以简单理解为一个名为 Test::func(int) 的全局函数,能被该类的所有实例复用,但独立于该类的任何实例,只属于该类作用域,可以不通过类的实例进行调用,也可以像普通成员函数一样通过类的实例进行调用

特定范围篇

  为了使全局变量或函数只在特定 cpp 文件中起作用,需要在 cpp 文件中相应变量或函数前添加static 修饰,如下表:

类型.h 文件中.cpp 文件中
全局变量不使用 static 修饰,使用 extern 修饰使用 static 修饰
全局函数不使用 static 修饰使用 static 修饰
  • 如果在头文件中声明 static 全局变量,则在包含该头文件的每个 .cpp 文件中都会生成一个独立的同名变量,而这种写法没有任何意义;如果在 .cpp 文件中不使用 static 声明全局变量,则该全局变量可能会被其它 .cpp 文件共享,也可能不会,造成该变量的不确定性;所以如果该全局变量要被所有 .cpp 文件共享,则需要在头文件中声明 extern 全局变量(eg:extern int g_value; // 注意,不要初始化值!),再在某个 .cpp 文件中单独进行定义和初始化(仅一次)(eg:int g_value = 0; // 不要extern修饰!),如此即可在每个 .cpp 文件中共享该全局变量;而若只想在单个 .cpp 文件中使用全局变量,则需要在该 .cpp 文件中全局范围类声明和定义 static int g_value = 0;,如此可保证该变量能且只能被该 .cpp 文件使用。
  • 如果在 .cpp 文件中不使用 static 声明全局函数,则该全局函数可能会被其它 .cpp 文件共享,也可能不会,这样在别的 .cpp 文件调用同名函数时可能会出现问题;而在头文件中使用 static 声明全局函数同样没有任何意义;所以如果要被多个 .cpp 文件复用,就将其声明移到头文件中,且不需要 static 修饰,而若只想在特定 .cpp 文件中使用该全局函数,则需要在声明时添加 static 修饰。

最后,若是在 .hpp 文件中,则需要去除全局对象,将全局函数封装为类的静态方法。

  PS:若在函数中使用 static 修饰变量,则该函数无法做到线程安全,在程序运行过程中,每次调用该函数,函数内的 static 变量都将保持上一次调用结束后的值,所以在函数中慎用 static 变量,除非需要这个特性

后记

  写这篇文章的初衷在于时常需要 static 时老是忘记或弄混它的用法,不得不去网上查找,虽说网上的相关资料也有很多,但在找的时候还是有点麻烦,毕竟有很多不是自己需要的,而且自己总结一下对其理解又更深一些,下次要用时也能马上找到自己所需。

参考资料

[1] c/c++ static 用法总结(三版本合一)

[2] C++中static的用法总结

[3] C++ 类中的static成员的初始化和特点

[4] C++静态成员函数访问非静态成员的几种方法


  1. 这里的一般是指局部变量,若为全局变量则默认为 extern ,局部变量没有默认初值,其初值不确定,一般需要人为明确的赋初值,而全局变量默认初值为 0 ,一个比较好的编程习惯是声明一个变量就对其进行初始化(赋初值),尽量少用全局变量,全局变量显示声明 extern。↩︎

  2. ※注:这里的 auto 与 C++11 中的意义不同,这里的 auto 指的是变量的存储形式,而不是 C++11 那种可以当做任意的变量类型,eg: int, double, std::vector<std::vector<double>>, ...... ,与其对应的还有 extern 和 register 关键字,其中 register 关键字基本不用 。↩︎

斐波那契数列的三种写法

前言

  本文预示着 Shaun 开始着手准备找工作的事了,初步计划是先把『剑指Offer』上的题先做一遍,对照着 牛客网 上的题进行测试,尽量争取先把书上的题都能 AC 。一般定义的斐波那契数列数列为:0,1,1,2,3,5,8……(对应 \(F(0)=0, F(1)=1, F(2)=1, \cdots \cdots\)),用数学公式表示即为:\(F(n)=F(n-1)+F(n-2)\)。以下代码均用 C++ 实现,且均通过牛客的测试。

前言

  本文预示着 Shaun 开始着手准备找工作的事了,初步计划是先把『剑指Offer』上的题先做一遍,对照着 牛客网 上的题进行测试,尽量争取先把书上的题都能 AC 。一般定义的斐波那契数列数列为:0,1,1,2,3,5,8……(对应 \(F(0)=0, F(1)=1, F(2)=1, \cdots \cdots\)),用数学公式表示即为:\(F(n)=F(n-1)+F(n-2)\)。以下代码均用 C++ 实现,且均通过牛客的测试。

循环写法

1
2
3
4
5
6
7
8
9
10
11
int fibonacci_loop(const unsigned int &n)
{
int fn = 0, f1 = 0, f2 = 1;
for (int i = 0; i < n; i++)
{
fn = f1 + f2;
f2 = f1;
f1 = fn;
}
return fn;
}

递归写法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int fibonacci_recursive(const unsigned int &n)
{
if (n == 0)
{
return 0;
}
else if (n == 1 || n == 2)
{
return 1;
}
else
{
return fibonacci_recursive(n - 1) + fibonacci_recursive(n - 2);
}
}

  ※注:这里 Shaun 在牛客上进行测试的时候,如果把 || n == 2 去掉的话,就没法通过,可见多递归一次花费的时间并不是线性增长的。

尾递归写法

  说来惭愧,这个概念还是在一个小学弟那里得知的,后面才逐渐了解并学会使用。

  尾递归,简而言之就是最后会且仅会调用函数本身,递归调用函数之后没有其它的语句需要执行。就像上面的递归,它在递归调用之后还会执行加法运算,而尾递归在执行递归调用之后就没有其它的运算了。

1
2
3
4
5
6
7
8
9
10
11
int fibonacci_tailRecursive(unsigned int n, unsigned int f1 = 1, unsigned int fn = 0)
{
if (n == 0)
{
return fn;
}
else
{
return fibonacci_tailRecursive(n - 1, fn, fn + f1);
}
}

总结

  循环和尾递归花费的时间和空间都差不多,都要比普通的递归要小,普通的递归优势在于便于理解,代码好写,在不强调性能的前提下,用递归写法的代码可读性可能要好些。

参考资料

[1] 递归与尾递归总结http://www.cnblogs.com/Anker/category/436371.html

LaTex语法指北

前言

  本文就简单记录一下使用 LaTex 可能会用到的一些语法,LaTex 的语法公式需要用符号 $ 括起来,其中 $...$ 代表行内公式,$$...$$ 代表行间公式。

前言

  本文就简单记录一下使用 LaTex 可能会用到的一些语法,LaTex 的语法公式需要用符号 $ 括起来,其中 $...$ 代表行内公式,$$...$$ 代表行间公式。

常用篇

§ 数学符号

名称示例预览名称示例预览
属于\in\(\in\)不属于\notin\(\notin\)
并集\cup\bigcup\(\cup\)\(\bigcup\)交集\cap\bigcap\(\cap\)\(\bigcap\)
包含于\subset\(\subset\)真包含于\subsetneqq\(\subsetneqq\)
包含\supset\(\supset\)真包含\supsetneqq\(\supsetneqq\)
不包含于\not\subset\(\not\subset\)空集\emptyset\(\emptyset\)
左大括号\{\(\{\)右大括号\}\(\}\)
正比于\propto\(\propto\)省略号\cdots\ddots\vdots\(\cdots\)\(\ddots\)\(\vdots\)
不等于\neq\(\neq\)等价于\Leftrightarrow\Longleftrightarrow\(\Leftrightarrow\)\(\Longleftrightarrow\)
左推出\Leftarrow\Longleftarrow\(\Leftarrow\)\(\Longleftarrow\)右推出\Rightarrow\Longrightarrow\(\Rightarrow\)\(\Longrightarrow\)
大于等于\ge\geq\(\ge\)\(\geq\)小于等于\le\leq\(\le\)\(\leq\)

§ 希腊字母

名称示例(小写和大写)预览(小写和大写)名称示例(小写和大写)预览(小写和大写)
alpha\alpha\Alpha\(\alpha\)\(\Alpha\)beta\beta\Beta\(\beta\)\(\Beta\)
gamma\gamma\Gamma\(\gamma\)\(\Gamma\)delta\delta\Delta\(\delta\)\(\Delta\)
theta\theta\vartheta\Theta\(\theta\)\(\vartheta\)\(\Theta\)lambda\lambda\Lambda\(\lambda\)\(\Lambda\)
mu\mu\Mu\(\mu\)\(\Mu\)pi\pi\Pi\(\pi\)\(\Pi\)
sigma\sigma\Sigma\(\sigma\)\(\Sigma\)omega\omega\Omega\(\omega\)\(\Omega\)
chi(卡方分布)\chi\Chi\(\chi\)\(\Chi\)phi\phi\Phi\(\phi\)\(\Phi\)

§ 运算符

名称示例预览名称示例预览
点乘\cdot\(\cdot\)叉乘\times\(\times\)
除号\div\(\div\)分数\frac{c}{a+b}\(\frac{c}{a+b}\)
累加\sum\(\sum\)累乘\prod\(\prod\)
积分\int\(\int\)

特殊篇

§ 公式书写

名称示例预览
大括号多行公式\begin{cases} x=v , \\ y=w \end{cases}\(\begin{cases} x=v , \\ y=w \end{cases}\)

§ 矩阵书写

示例语法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
\begin{matrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9
\end{matrix}
\quad
\begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix}
\quad
\begin{Bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{Bmatrix}
\quad
\begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}
\quad
\begin{vmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{vmatrix}
\qquad
\begin{Vmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{Vmatrix}

显示效果: \[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \quad \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix} \quad \begin{Bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{Bmatrix} \quad \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} \quad \begin{vmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{vmatrix} \quad \begin{Vmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{Vmatrix} \]

扩展篇

技巧篇

将下标放在正下方:

将普通下标 \max_{x<y} \(\max_{x<y}\) 改为 \max\limits_{x<y} \(\max\limits_{x<y}\),即在下标符号 _ 前添加 \limits

至于分块矩阵的虚线可参考: 数式: 行列、ベクトル matrix, vector TeX

字符加粗: 通过 \mathbf{X}, eg. \(\mathbf{X}\)\(X\) ,

后记

  写 Wiki 毕竟不是一蹴而就的事,等以后碰到一些常用的或有时间在继续更新吧 ↖(^ ω ^)↗。

参考资料

[1] 常用数学符号的 LaTeX 表示方法

[2] latex中的希腊字母

[3] Latex中将下标放在正下方

搜索技巧

前言

  Shaun 一直以为自己的搜索能力还可以,基本上自己想要的东西都能搜到,但是自从接触到这个世界,才知道自己大概只是个入门的水平(或者说连入门都说不上 /つ∇T)),和网上的一些大神相比还有比较大的差距,此文只是网上一些资料的整理,方便 Shaun 熟练使用,提高驾驶技巧,毕竟信息检索能力还是很重要的。

前言

  Shaun 一直以为自己的搜索能力还可以,基本上自己想要的东西都能搜到,但是自从接触到这个世界,才知道自己大概只是个入门的水平(或者说连入门都说不上 /つ∇T)),和网上的一些大神相比还有比较大的差距,此文只是网上一些资料的整理,方便 Shaun 熟练使用,提高驾驶技巧,毕竟信息检索能力还是很重要的。

搜索篇

Google 可以通过添加一些字符优化搜索结果,如:

搜索需求对应字符
搜索社交媒体在用于搜索社交媒体的字词前加上 @。例如:@twitter
从搜索结果中排除特定字词在您要排除的字词前加上 -。例如:jaguar speed -car
搜索完全匹配的结果为字词或短语加上引号。例如:"tallest building"
搜索通配符或未知字词在字词或短语中您要放置占位符的地方加上 *。例如:"largest * in the world"
在某个数字范围内执行搜索在两个数字之间加上 ..。例如:camera $50..$100
组合搜索在各个搜索查询之间加上“OR”。例如:marathon OR race
搜索特定网站在相应网站或域名前加上“site:”。例如:site:youtube.comsite:.gov
搜索相关网站在已知网址前加上“related:”。例如:related:time.com
查找链接到某个特定网页的网页在已知网址前加上“link:”。例如:link:chongbuluo.com,就能查到哪些网页中包含链接chongbuluo.com
查找在URL地址里有搜索关键词的页面在已知网址前加上“inurl:”。例如:inurl:chongbuluo,就能查到哪些网页 url 中包含链接chongbuluo
查找在网页标题里有搜索关键词的页面在已知网址前加上“intitle:”。例如:intitle:chongbuluo,就能查到哪些网页标题中包含链接chongbuluo
查找在网页正文里有搜索关键词的页面在已知网址前加上“intext:”。例如:intext:chongbuluo,就能查到哪些网页正文中包含链接chongbuluo
查找pdf,xml,xls,txt,doc,csv等特定格式的结果在特定文件格式前加上filetype:。例如 filetype:pdf
查找关键词的定义在关键词前加上define:。例如 define:搜索

  经 Shaun 实测,以上大部分字符对百度搜索引擎同样适用。当然,最好的搜索方式是使用高级搜索,高级搜索可以进行一系列设置,比如时间范围,从而使搜索结果更精确,更容易得到想要的结果,不过在不十分确定的时候,不要做太多限制,不然可能会过滤掉关键信息。

  Google 和百度对英文字符大小写都不敏感,搜索 QQqq 所得到的结果是一样的。搜索是否成功最关键的地方还是在于关键词的选取,关键词的选取这没什么好说的,只能提高搜索熟练度及对问题的把握程度了。

  哦,还有一点忘记说了,就是在 Google 中如果要搜索那个的话,必须要在搜索设置里面关闭安全搜索功能,如果简体中文不能关闭安全搜索功能的话,就在搜索设置里将语言更换到繁体中文或英语应该就能关闭安全搜索功能。

技巧篇

  有时候第一个页面没有想要的结果,于是需要看第二个、第三个页面的结果,是不是觉得翻页很麻烦?(如果不觉得麻烦可以直接跳过🙄),但是设置增加搜索结果条目又会提高显示延迟,这时可以使用 Super_preloaderPlus_one 脚本(Chrome 中可以使用 AutoPagerize 插件),只需鼠标继续往下滚轮就会自动加载下页搜索结果,无需点击翻页。

  当点击网页链接却出现404错误或无法显示页面的错误时,这个时候可以使用搜索引擎的「网页快照」功能,Google的这个功能点击搜索结果页面标题下的 绿色小三角 即可看到,百度的这个叫 百度快照,在搜索结果 url 地址的末端,即结果最后面。快照功能最强大的是 互联网档案馆Internet Archive),又叫『网站时光倒流机器』(Wayback Machine),在知道链接的情况下,将链接输入到 Internet Archive 点击搜索就能查看历史快照了。

  如果要查找某个网页出现的关键词,可以利用 Chrome 和 Firefox 的网页搜索功能,一般可以通过 Ctrl + F 去检索。也可以通过 F12 或「鼠标右键」(Firefox 点击「查看元素」,Chrome 点击「检查」,都是右键弹框最后一个选项)进入浏览器控制台开发工具,Chrome 中还要通过 Ctrl + F 才能检索元素,而 Firefox 可以直接搜索 HTML,输入要查找的关键词,回车即可。

专项篇

  图像搜索:这个一般用 Google 和百度的以图搜图就可以了,不过也有些特殊的图像搜索引擎,如 TinEye 等,当然,图像搜索最好还是安装相应的插件,在 Chrome 中可以安装一个叫 二箱 的插件,在 Firefox 中可以安装一个叫 Search by Image 的插件。

  音乐搜索:这个如果能听出歌词的话就直接去搜索引擎上搜听到的歌词,如果无法听出歌词的话,可以尝试各大音乐软件的「听歌识曲」功能。

  当然更精确的专项搜索一般只存在于特定的网站,这就要看对整个互联网的了解程度了,这里强推「虫部落-快搜」,不管是日常搜索需求还是特殊搜索需求基本都能满足。

后记

  尽量优先使用 Google,百度或许能找到自己想要的,但太浪费时间了,浪费时间就是浪费生命,不过能看到本文的看官,应该也是能熟练使用 Google 的了。

参考资料

[1] 如何在 Google 中进行搜索(https://support.google.com/websearch/#topic=3081620)

[2] 优化网页搜索

[3] 提高搜索能力的关键技巧(如何查找可靠出处)

[4] 谷歌搜索技巧:搜索语法+隐藏彩蛋+高级设置

矩阵的应用之图像仿射变换

前言

  好像很久没写新的东西了,主要是最近期末有一大堆事情要做,忙着写各种结课论文和复习数学,又加上最近忙着把《奥日与黑暗森林:终极版》剧情通关,这游戏不管是画面还是音乐都特别棒,但是对键盘用户太不友好了,不能改键位,需要一只手控制键盘,一只手控制鼠标,如果只是普通的键位就算了,它还要万恶的 shift 键配合,手残党完全吃不消 (´・ω・`)。最终死亡 658 次好歹剧情通关了,最后悔的是没把三段跳点出来 /つ∇T)。这些事情一搞完,Shaun 这不就又开始写了嘛( 说什么也摆脱不了你拖延症晚期的事实 ( ̄ε(# ̄)☆╰╮( ̄▽ ̄///) )。 选修的《数字图像处理》是也早结课了,虽然教的东西大都事前已经了解了,但是好像还没写过图像处理相关的 blog,所以谨以此篇最基础的 blog 来表示一下。

前言

  好像很久没写新的东西了,主要是最近期末有一大堆事情要做,忙着写各种结课论文和复习数学,又加上最近忙着把《奥日与黑暗森林:终极版》剧情通关,这游戏不管是画面还是音乐都特别棒,但是对键盘用户太不友好了,不能改键位,需要一只手控制键盘,一只手控制鼠标,如果只是普通的键位就算了,它还要万恶的 shift 键配合,手残党完全吃不消 (´・ω・`)。最终死亡 658 次好歹剧情通关了,最后悔的是没把三段跳点出来 /つ∇T)。这些事情一搞完,Shaun 这不就又开始写了嘛( 说什么也摆脱不了你拖延症晚期的事实 ( ̄ε(# ̄)☆╰╮( ̄▽ ̄///) )。 选修的《数字图像处理》是也早结课了,虽然教的东西大都事前已经了解了,但是好像还没写过图像处理相关的 blog,所以谨以此篇最基础的 blog 来表示一下。

预备篇

  首先需要了解的是图像中坐标系和数学课本中常用的坐标系略有不同,图像中坐标系是以左上角为原点,水平向右为 X 轴,垂直向下为 Y 轴;而数学课本中常见的坐标系是以图像中心为原点,水平向右为 X 轴,垂直向上为 Y 轴。所以由图像中坐标 \((x,y)\) 转数学课本中常见的坐标 \((x',y')\) 的公式为 $x' = x - C/2; y' = -y + R/2; $ 其中 \(C\) 表示原图像总列数,即原图像宽度,\(R\) 表示原图像总行数,即原图像高度。

X 轴Y 轴X 轴Y 轴数学坐标系图像坐标系

  ※注:值得注意的是因为 MATLAB 和 OpenCV 的像素索引坐标形式为 (行坐标,列坐标) ,所以若以本文这样定的图像坐标,则图像中坐标 \((x,y)\) 对应的像素值为 \(f(y,x)\)

变换篇

据 OpenCV 文档所说:

  1. 一个任意的仿射变换都能表示为 乘以一个矩阵 (线性变换) 接着再 加上一个向量 (平移).

  2. 综上所述, 我们能够用仿射变换来表示:

    1. 旋转 (线性变换)
    2. 平移 (向量加)
    3. 缩放操作 (线性变换)

    你现在可以知道, 事实上, 仿射变换代表的是两幅图之间的 关系 .

  3. 我们通常使用 2 矩阵来表示仿射变换.

    A = \begin{bmatrix}     a_{00} & a_{01} \\     a_{10} & a_{11}     \end{bmatrix}_{2 \times 2} B = \begin{bmatrix}     b_{00} \\     b_{10}     \end{bmatrix}_{2 \times 1} M = \begin{bmatrix}     A & B     \end{bmatrix} =\begin{bmatrix}     a_{00} & a_{01} & b_{00} \\     a_{10} & a_{11} & b_{10}\end{bmatrix}_{2 \times 3}  

    考虑到我们要使用矩阵 AB 对二维向量 X = \begin-{bmatrix}x \ y\end-{bmatrix} 做变换, 所以也能表示为下列形式:

    T = A \begin-{bmatrix}x \ y\end-{bmatrix} + B or T = M ^{T}

    T=\begin-{bmatrix}a_{00}x + a_{01}y + b_{00} \ a_{10}x + a_{11}y + b_{10}\end-{bmatrix}  

而在冈萨雷斯的《数字图像处理_第三版》里有:

最常用的空间坐标变换之一是仿射变换,其一般形式如下: \[ \begin{bmatrix} x & y & 1 \end{bmatrix} = \begin{bmatrix} v & w & 1 \end{bmatrix} \textbf{T} = \begin{bmatrix} v & w & 1 \end{bmatrix} \begin{bmatrix} t_{11} & t_{12} & 0 \\ t_{21} & t_{22} & 0 \\ t_{31} & t_{32} & 1 \end{bmatrix} \] 其中 \((v,w)\) 为原坐标,\((x,y)\) 为变换后的坐标,可根据变换矩阵 \(\textbf{T}\) 中的元素选择的值,对一组坐标点做尺度、旋转、平移或偏移变换。一些常见的变换矩阵及作用如下表:

变换名称仿射变换矩阵\(\textbf{T}\)坐标公式
恒等变换\(\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\)\(\begin{cases} x=v , \\ y=w \end{cases}\)
尺度变换\(\begin{bmatrix} c_x & 0 & 0 \\ 0 & c_y & 0 \\ 0 & 0 & 1 \end{bmatrix}\)\(\begin{cases} x=vc_x , \\ y=wc_y \end{cases}\)
旋转变换(以逆时针为正)\(\begin{bmatrix} cos(\theta) & sin(\theta) & 0 \\ -sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}\)\(\begin{cases} x=vcos(\theta)-wsin(\theta) , \\ y=vsin(\theta)+wcos(\theta) \end{cases}\)
旋转变换(以顺时针为正)\(\begin{bmatrix} cos(\theta) & -sin(\theta) & 0 \\ sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}\)\(\begin{cases} x=vcos(\theta)+wsin(\theta) , \\ y=-vsin(\theta)+wcos(\theta) \end{cases}\)
平移变换\(\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ t_x & t_y & 1 \end{bmatrix}\)\(\begin{cases} x=v+t_x , \\ y=w+t_y \end{cases}\)
偏移变换(水平)\(\begin{bmatrix} 1 & 0 & 0 \\ s_h & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\)\(\begin{cases} x=v+ws_h , \\ y=w \end{cases}\)
偏移变换(垂直)\(\begin{bmatrix} 1 & s_v & 0 \\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\)\(\begin{cases} x=v , \\ y=vs_v+w \end{cases}\)

仿射变换的实现由两种方式:一种是 前向映射(Forward Mapping):直接采用利用原图像坐标 \((v,w)\) 通过 \(\begin{bmatrix} x & y & 1\end{bmatrix}=\begin{bmatrix} v & w & 1\end{bmatrix} \textbf{T}\) 得到变换后的坐标 \((x,y)\),使用前向映射会导致一些问题:可能会有多个像素坐标映射到输出图像的同一位置,也可能输出图像的某些位置完全没有相应的输入图像像素与它匹配,也就是没有被映射到,造成有规律的空洞(黑色的花纹状);更好的一种方式是采用 反向映射(Inverse Mapping):扫描输出图像的位置 \((x,y)\),通过 \(\begin{bmatrix} v & w & 1\end{bmatrix}= \begin{bmatrix} x & y & 1\end{bmatrix}\textbf{T}^{-1}\)(其中 \(\textbf{T}^{-1}\)\(\textbf{T}\) 的逆矩阵)计算输入图像对应的位置 \((v,w)\),通过插值方法决定输出图像该位置的灰度值。

  本文这里采取冈萨雷斯的《数字图像处理_第三版》的变换矩阵方式,毕竟所学的矩阵论也是将变换矩阵放在后面作为第二个因子。虽然仿射变换都有现成的 API 可以调用,而且速度一般要比自己写的要快,但是知其然终究也要知其所以然。

  下面就以旋转变换为例了,因为尺度变换和平移变换只需要相应的缩放图像即可,而旋转变换不仅需要更改图像大小,还要确定旋转中心,而旋转中心一般以图像中心为标准。

旋转变换

  旋转变换首先需要确定旋转中心,若以图像左上角为旋转中心,则只需要像做尺度变换和平移变换那样通过 \(\begin{bmatrix} v & w & 1\end{bmatrix}= \begin{bmatrix} x & y & 1\end{bmatrix}\textbf{T}^{-1}\) 做普通的变换即可。而以图像中心为旋转中心,首先需要做坐标变换,将以左上角为原点,水平向右为 X 轴,垂直向下为 Y 轴的图像坐标系转换为以图像中心为原点,水平向右为 X 轴,垂直向下为 Y 轴的数学坐标系;再做正常的旋转变换;随后再将数学坐标系转换为图像坐标系,所以图像中心为旋转中心的旋转变换总共需要做三次变换。这里就以图像中心为旋转中心为例,由于有三次变换,所以应该有三个变换矩阵相乘,设 图像坐标系==》数学坐标系的变换矩阵为 T1,旋转变换矩阵为 T2,数学坐标系==》图像坐标系的变换矩阵为 T3,设顺时针旋转角度为 \(\theta\) ,原图像宽度为 \(C\),高度为 \(R\),旋转后图像宽度为 \(W\),高度为 \(H\),则: \[ T1=\begin{bmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ -0.5C & 0.5R & 1 \end{bmatrix}   T2=\begin{bmatrix} cos(\theta) & -sin(\theta) & 0 \\ sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}   T3=\begin{bmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0.5W & 0.5H & 1 \end{bmatrix} \] 则旋转变换最终形式为:\(\begin{bmatrix} x & y & 1 \end{bmatrix}=\begin{bmatrix} v & w & 1\end{bmatrix} \textbf{T}=\begin{bmatrix} v & w & 1 \end{bmatrix}T1*T2*T3\)

※BTW:旋转变换中,旋转后图像宽度 \(W\),高度 \(H\) 与 原图像宽度 \(C\),高度 \(R\) 的关系为: \[ \begin{cases} H = |R*cos(\theta)| + |C*sin(\theta)| , \\ W = |C*cos(\theta)| + |R*sin(\theta)| \end{cases} \]

插值篇

  因为经过采用反向映射1方式的仿射变换之后,得到的原图像坐标 \((v,w)\) 往往不是整数值,所以无法知道其对应的像素值 \(f(w,v)\),这时需要采取插值的方式近似估计该坐标位置的像素值。

  常用的插值方法有最近邻插值(nearest neighbor interpolation)、双线性插值(bilinear interpolation)和双三次插值(bicubic interpolation),其中双三次插值在保持图像细节方面最好,但花费时间也最多,PS 中的消除锯齿和羽化效果好像就采用了双三次插值。

  最近邻插值只考虑相邻最近的像素,双线性插值考虑相邻的 4 个像素点,双三次插值则考虑相邻的 16 个像素点。

最近邻插值

  最近邻插值最简单,将变换后图像的坐标 \((x,y)\) 通过反向映射得到原图像坐标 \((v,w)\) ,直接对 \((v,w)\) 进行四舍五入得到相应的整数坐标 \((⌊v+0.5⌋,⌊w+0.5⌋)\)2,用该整数坐标的像素值近似估计 \((v,w)\) 的像素值,令 \(f(y,x)=f(⌊w+0.5⌋,⌊v+0.5⌋)\) ,从而得到变换后图像每个像素点的像素值。

双线性插值

  双线性插值是线性插值方法的一种扩展,它是 X 和 Y 两个方向上线性插值的组合。

X 轴Y 轴P1P2P3P4Z1Z2P(v,w)

如上图,设变换后图像的坐标 \((x,y)\) 通过反向映射得到原图像坐标 \((v,w)\) ,即点 \(P\) 正好处于四个像素点 \(P1(v_0, w_0)\)\(P2(v_0+1, w_0)\)\(P3(v_0+1, w_0+1)\)\(P4(v_0, w_0+1)\) 的中间,其中 \(v_0=⌊v⌋\)\(w_0=⌊w⌋\) ,点 \(P\) 对应的像素值为 \(f(P)\) 因为双线性插值即在 \(X\)\(Y\) 两个方向进行线性插值,首先计算 \(X\) 方向的插值: \[ \begin{cases} \frac{f(P2)-f(P1)}{P2.x-P1.x}=\frac{f(Z1)-f(P1)}{Z1.x-P1.x} , \\ \frac{f(P3)-f(P4)}{P3.x-P4.x}=\frac{f(Z2)-f(P4)}{Z2.x-P4.x} \end{cases} \] 即: \[ \begin{cases} f(Z1)=\frac{Z1.x-P1.x}{P2.x-P1.x}f(P2)+\frac{P2.x-Z1.x}{P2.x-P1.x}f(P1) =(v-v_0)f(P2)+(v_0+1-v)f(P1), \\ f(Z2)=\frac{Z2.x-P4.x}{P3.x-P4.x}f(P3)+\frac{P3.x-Z2.x}{P3.x-P4.x}f(P4) =(v-v_0)f(P3)+(v_0+1-v)f(P4) \end{cases}\tag{1} \] 然后计算 \(Y\) 方向的插值: \[ \begin{equation} \frac{f(Z2)-f(Z1)}{Z2.y-Z1.y}=\frac{f(P)-f(Z1)}{P.y-Z1.y} \end{equation} \] 即: \[ \begin{equation} f(P)=\frac{P.y-Z1.y}{Z2.y-Z1.y}f(Z2)+\frac{Z2.y-P.y}{Z2.y-Z1.y}f(Z1) =(w-w_0)f(Z2)+(w_0+1-w)f(Z1) \end{equation}\tag{2} \] 结合式(1)和式(2)可得: \[ \begin{equation} f(P)=(v_0+1-v)(w_0+1-w)f(P1)+(v-v_0)(w_0+1-w)f(P2)+(v-v_0)(w-w_0)f(P3)+(v_0+1-v)(w-w_0)f(P4) \end{equation} \] 用矩阵形式可表示为: \[ f(P)=\begin{bmatrix} v_0+1-v & v-v_0 \end{bmatrix} \begin{bmatrix} f(P1) & f(P4) \\ f(P2) & f(P3) \end{bmatrix} \begin{bmatrix} (w_0+1-w) \\ (w-w_0) \end{bmatrix} \] 即: \[ f(y,x)=f(w,v)=\begin{bmatrix} ⌊v⌋+1-v & v-⌊v⌋ \end{bmatrix} \begin{bmatrix} f(⌊w⌋,⌊v⌋) & f(⌊w⌋+1,⌊v⌋) \\ f(⌊w⌋,⌊v⌋+1) & f(⌊w⌋+1,⌊v⌋+1) \end{bmatrix} \begin{bmatrix} (⌊w⌋+1-w) \\ (w-⌊w⌋) \end{bmatrix} \]

双三次插值

X 轴Y 轴P11P(v,w)

  如上图,双三次插值需要考虑相邻16个像素(4×4),用双三次插值重采样的图像更平滑并且更能保留图像细节,在这三种插值算法中,双三次插值效果最好,但处理速度最慢。同样设变换后图像的坐标 \((x,y)​\) 通过反向映射得到原图像坐标 \((v,w)​\) ,与其左上角相邻最近的 点P11 坐标则为 \((⌊v⌋,⌊w⌋)​\) ,该插值方法需要选取一个合适的插值基函数,参照维基百科 Bicubic interpolation 的一般为: \[ W(x) = \begin{cases} (a+2)|x|^3-(a+3)|x|^2+1 & \text{for } |x| \leq 1, \\ a|x|^3-5a|x|^2+8a|x|-4a & \text{for } 1 < |x| < 2, \\ 0 & \text{otherwise}, \end{cases} \] 其中 \(a\) 一般取 -0.5 、-0.75 或 -1;则:\(f(y,x)=f(w,v)=A*B*C\) ,其中: \[ A=\begin{bmatrix} W( v-(⌊v⌋-1) ) & W(v-⌊v⌋) & W( (⌊v⌋+1)-v ) & W( (⌊v⌋+2)-v ) \end{bmatrix} \\ B=\begin{bmatrix} f(⌊w⌋-1,⌊v⌋-1) & f(⌊w⌋,⌊v⌋-1) & f(⌊w⌋+1,⌊v⌋-1) & f(⌊w⌋+2,⌊v⌋-1) \\ f(⌊w⌋-1,⌊v⌋) & f(⌊w⌋,⌊v⌋) & f(⌊w⌋+1,⌊v⌋) & f(⌊w⌋+2,⌊v⌋) \\ f(⌊w⌋-1,⌊v⌋+1) & f(⌊w⌋,⌊v⌋+1) & f(⌊w⌋+1,⌊v⌋+1) & f(⌊w⌋+2,⌊v⌋+1) \\ f(⌊w⌋-1,⌊v⌋+2) & f(⌊w⌋,⌊v⌋+2) & f(⌊w⌋+1,⌊v⌋+2) & f(⌊w⌋+2,⌊v⌋+2) \end{bmatrix} \\ C=\begin{bmatrix} W( w-(⌊w⌋-1) ) \\ W(w-⌊w⌋) \\ W( (⌊w⌋+1)-w ) \\ W( (⌊w⌋+2)-w ) \end{bmatrix} \]

即: \[ f(y,x)=f(w,v)= \sum\limits_{row=-1}^2\sum\limits_{col=-1}^2f(⌊w⌋+row,⌊v⌋+col)W(row-(w-⌊w⌋))W(col-(v-⌊v⌋)) \] 另附:网上也有人中间那个矩阵 \(B\) 是本文中间矩阵 \(B\) 的转置,经过下文实践,感觉效果差不多,但从理论上来说,应该本文这样写才是对的吧🤔。

Lanczos 插值

  Lanczos 插值和双三次插值本质上差不多,不同的是插值权重基函数不一样,Lanczos 插值算法的基函数为: \[ L(x) = \begin{cases} 1 & \text{for } x = 0, \\ a*sin(\pi x)sin(\pi x/a)/(\pi x)^2 & \text{for } 0 < |x| \leq a, \\ 0 & \text{otherwise}, \end{cases} \] 其中 a 代表插值窗口半径,若 a = 2,则代表取周围 4×4 邻域像素点进行插值,和双三次插值除权重不一样外其它都一样,此时适用于对图像下采样(缩小),若 a = 3,则代表取周围 9×9 邻域像素点进行插值,此时适用于对图像上采样(放大),OpenCV 的 INTER_LANCZOS4 插值参数就代表取周围 8×8 邻域的像素点进行插值。

实践篇

  本次实践采用 Matlab R2016b,具体 matlab 实现代码为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
clc;clear;close all;

img = imread('lena_gray.jpg'); % 读取图像
[R, C] = size(img); % 获取图像大小

theta = 45 * pi / 180.0; % 旋转角度
H = ceil(abs(R*cos(theta)) + abs(C*sin(theta))); % 变换后图像的高度
W = ceil(abs(C*cos(theta)) + abs(R*sin(theta))); % 变换后图像的宽度
res = zeros(H, W); % 构造结果矩阵。每个像素点默认初始化为0(黑色)

T1 = [1 0 0; 0 -1 0; -0.5*C 0.5*R 1]; % 将原图像坐标映射到数学笛卡尔坐标
T2 = [cos(theta) -sin(theta) 0; sin(theta) cos(theta) 0; 0 0 1]; % 数学笛卡尔坐标下顺时针旋转的变换矩阵
T3 = [1 0 0; 0 -1 0; 0.5*W 0.5*H 1]; % 将数学笛卡尔坐标映射到旋转后的图像坐标
T = T1*T2*T3;
inv_T = inv(T); % 求逆矩阵
% inv_T = [cos(theta) -sin(theta) 0; sin(theta) cos(theta) 0; -0.5*W*cos(theta)-0.5*H*sin(theta)+0.5*C 0.5*W*sin(theta)-0.5*H*cos(theta)+0.5*R 1];

for y = 1 : H % 变换后图像的纵坐标,行,高
for x = 1 : W % 变换后图像的横坐标,列,宽
original_coordinate = [x y 1] * inv_T; % 矩阵乘法
v = original_coordinate(1); % 原图像的横坐标,列,宽
w = original_coordinate(2); % 原图像的纵坐标,行,高
% 变换后的位置判断是否越界
if v>=1 && w>=1 && v<=C && w<=R
res(y, x) = img(round(w), round(v)); % 用原图像对应坐标的像素值填充变换后的图像(最邻近插值)

% ------------- 双线性插值(bilinear interpolation)-----------------
left = floor(v); right = ceil(v); top = floor(w); bottom = ceil(w);
dC = v-left; % 列偏差
dR = w-top; % 行偏差
res(y, x) = (1-dC)*(1-dR)*img(top, left) + dC*(1-dR)*img(top,right) + (1-dC)*dR*img(bottom, left) + dC*dR*img(bottom, right);

% ------------- 双三次插值(bicubic interpolation) -------------------------
if left>=2 && top>=2 && left<=(C-2) && top<=(R-2)
img = double(img);
MA = [bicubic(1+dC) bicubic(dC) bicubic(1-dC) bicubic(2-dC)];
MB = [img(top-1,left-1) img(top,left-1) img(top+1,left-1) img(top+2,left-1);
img(top-1,left) img(top,left) img(top+1,left) img(top+2,left);
img(top-1,left+1) img(top,left+1) img(top+1,left+1) img(top+2,left+1);
img(top-1,left+2) img(top,left+2) img(top+1,left+2) img(top+2,left+2)];
% MB = MB'; % 求转置矩阵
MC = [bicubic(1+dR); bicubic(dR); bicubic(1-dR); bicubic(2-dR)];
res(y, x) = MA*MB*MC;
end
end
end
end;

figure, imshow(uint8(res)); % 显示图像

BiCubic 基函数 Matlab 代码为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function W = bicubic(x)
%bicubic 双三次插值基函数
a = -1; % 默认取a为-1
x1 = abs(x);
x2 = x1*x1;
x3 = x1*x2;

if x1 <= 1
W = 1 - (a+3)*x2 + (a+2)*x3;
elseif x1>1 && x1<=2
W = -4*a + 8*a*x1 - 5*a*x2 + a*x3;
else
W = 0;
end

  旋转变换中感觉插值的作用没体现出来,以肉眼来看感觉三种插值方法的效果差不多,可能是 Shaun 选取的示例不好,为了体现插值效果,应该采用尺度变换(缩放变换)的。以上代码改为尺度变换也简单,自定义图像缩放后的宽高,以两倍为例,H = R * 2; W = C * 2;,再将旋转变换矩阵改为尺度变换矩阵,尺度变换矩阵中 \(c_x=W/C;c_y=H/R\)

为了便于理解,Shaun 对代码就不进行优化了(其实是你懒吧 _(:з」∠)_)。

后记

  本文算是数字图像处理中最基础的知识了,但 Shaun 在写时还是查阅了大量相关的资料,有些地方理解的还不是很透彻,行文思路有点混乱 ╮(╯▽╰)╭。本来是不想使用图片的,但本文不用图片很难理解清楚,又为了不使用外部图,最后只得参考 SVG 教程如何创建SVG箭头和polymarker——marker元素 采用 SVG 绘制相应图片了。等有时间再把用 OpenCV 实现的 C++ 代码也贴上吧。最后再感叹一下 Matlab 确实是做科研的好工具(°Д°)Ъ,吐槽一下 MathJax 排版好痛苦啊,太多需要转义符\的地方了吧。,搞错了 Σ(゚д゚;),这主要和 markdown 渲染有关,hexo 默认的 markdown 渲染插件 hexo-renderer-marked 太普通了,有些东西根本没办法渲染或者渲染有问题 (╯‵□′)╯︵┴─┴,Shaun 最后决定使用 hexo-renderer-pandoc 插件渲染 markdown,这样就完美了 (๑•̀ㅂ•́)و✧ 。

  至于具体怎么使用 hexo-renderer-pandoc 替换默认的渲染器可参考:如何禁止 hexo 在 html 代码里插入<br>标签?。具体如下:

1
2
3
4
5
# 1、安装 Pandoc,可以不顺带安装 MiKTex
# 2、卸载默认渲染器
npm uninstall hexo-renderer-marked --save
# 3、安装 hexo-renderer-pandoc
npm install hexo-renderer-pandoc --save

附录

  挖的坑总是要填的,呐,这就是用 OpenCV 实现的旋转变换,实现语言为 C++ :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <opencv2/opencv.hpp>

#define M_PI 3.14159265358979323846

double bicubic(double x)
{
// bicubic 双三次插值基函数
int a = -1; // 默认取a为 - 1
double x1 = fabs(x);
double x2 = x1*x1;
double x3 = x1*x2;

if (x1 <= 1)
{
return 1 - (a + 3)*x2 + (a + 2)*x3;
}
else if (x1 > 1 && x1 <= 2)
{
return -4 * a + 8 * a*x1 - 5 * a*x2 + a*x3;
}
else
{
return 0;
}
}

int main(int argc, char *argv[])
{
cv::Mat img = cv::imread("../Data/lena_gray.jpg", 0); // 以灰度模式读取图片

int R = img.rows; // 获取原图像高度
int C = img.cols; // 获取原图像宽度

double theta = 45 * M_PI / 180.0; // 旋转角度

int H = ceil(fabs(R*cos(theta)) + fabs(C*sin(theta))); // 变换后图像的高度
int W = ceil(fabs(C*cos(theta)) + fabs(R*sin(theta))); // 变换后图像的宽度

cv::Mat res = cv::Mat::zeros(H, W, CV_8UC1); // 构造结果矩阵。每个像素点默认初始化为0(黑色)

cv::Mat T1 = (cv::Mat_<double>(3, 3) << 1, 0, 0, 0, -1, 0, -0.5*C, 0.5*R, 1); // 将原图像坐标映射到数学笛卡尔坐标
cv::Mat T2 = (cv::Mat_<double>(3, 3) << cos(theta), -sin(theta), 0, sin(theta), cos(theta), 0, 0, 0, 1); // 数学笛卡尔坐标下顺时针旋转的变换矩阵
double t3[3][3] = { { 1, 0, 0 }, { 0, -1, 0 }, { 0.5*W, 0.5*H, 1 } }; // 将数学笛卡尔坐标映射到旋转后的图像坐标
cv::Mat T3 = cv::Mat(3, 3, CV_64FC1, t3);
cv::Mat T = T1*T2*T3;
cv::Mat inv_T = T.inv(); // 求逆矩阵
//cv::Mat inv_T = (cv::Mat_<double>(3, 3) << cos(theta), -sin(theta), 0, sin(theta), cos(theta), 0, -0.5*W*cos(theta) - 0.5*H*sin(theta) + 0.5*C, 0.5*W*sin(theta) - 0.5*H*cos(theta) + 0.5*R, 1);

for (int y = 0; y < H; y++)
{
for (int x = 0; x < W; x++)
{
cv::Mat point = (cv::Mat_<double>(1, 3) << x, y, 1);
cv::Mat original_coordinate = point * inv_T; // 矩阵乘法
double v = original_coordinate.at<double>(0, 0); // 原图像的横坐标,列,宽
double w = original_coordinate.at<double>(0, 1); // 原图像的纵坐标,行,高

// 变换后的位置判断是否越界
if (v >= 0 && w >= 0 && v <= C - 1 && w <= R - 1)
{
res.at<uchar>(y, x) = img.at<uchar>(round(w), round(v)); // 用原图像对应坐标的像素值填充变换后的图像(最邻近插值)

// ------------ - 双线性插值(bilinear interpolation)---------------- -
int left = floor(v), right = ceil(v), top = floor(w), bottom = ceil(w);
double dC = v - left; // 列偏差
double dR = w - top; // 行偏差
res.at<uchar>(y, x) = (1 - dC)*(1 - dR)*img.at<uchar>(top, left) + dC*(1 - dR)*img.at<uchar>(top, right) + (1 - dC)*dR*img.at<uchar>(bottom, left) + dC*dR*img.at<uchar>(bottom, right);

// ------------ - 双三次插值(bicubic interpolation)------------------------ -
if (left >= 1 && top >= 1 && left <= (C - 3) && top <= (R - 3))
{
cv::Mat MA = (cv::Mat_<double>(1, 4) << bicubic(1 + dC), bicubic(dC), bicubic(1 - dC), bicubic(2 - dC));
cv::Mat MB = img(cv::Rect(left - 1, top - 1, 4, 4)); // 提取当前相邻区域16个像素点做插值
MB.convertTo(MB, CV_64FC1); // 变换为浮点型数据
MB = MB.t(); // 求转置矩阵
cv::Mat MC = (cv::Mat_<double>(4, 1) << bicubic(1 + dR), bicubic(dR), bicubic(1 - dR), bicubic(2 - dR));
cv::Mat result = MA*MB*MC;
res.at<uchar>(y, x) = static_cast<uchar>(result.at<double>(0, 0));
}
}
}
}

cv::imshow("result", res); // 显示变换后图像
cv::waitKey(0);

return 0;
}

  以上 C++ 代码在 VS2013 下能完美运行,不管是用 OpenCV-2.4.11 还是 OpenCV-3.2.0。其实完全理解的话,不管用什么工具都能实现,只是看哪个工具方便一点而已,就这个而言,感觉 Matlab 要方便很多,Shaun 就不继续挖 Python 的坑了,毕竟如果要用 Python 实现其实还是用 OpenCV,只是用 OpenCV Python 版的接口而已。

参考资料

[1] 第4章 图像几何变换

[2] 图像旋转原理及实现

[3] 仿射变换http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/tutorials.html

[4] 图像处理常用插值方法总结

[5] Wikipedia Bilinear interpolation

[6] 双线性插值算法的详细总结

[7] Wikipedia Bicubic interpolation

[8] 双三次插值(bicubic interpolation)原理及MATLAB源码实现

[9] 图像缩放】双立方(三次)卷积插值https://dailc.github.io/blog/tags.html#%E6%8F%92%E5%80%BC%E7%AE%97%E6%B3%95


  1. 为啥不说前向映射呢?这是因为若原图像坐标 \((v,w)\) 通过前向映射方式得到变换后图像的坐标 \((x,y)\) ,而且这个坐标为小数的话,一般采用四舍五入的方式得到变换后图像对应的整数坐标 \((⌊x+0.5⌋,⌊y+0.5⌋)\),令 \(f(⌊y+0.5⌋,⌊x+0.5⌋)=f(w,v)\)↩︎

  2. \(⌊x⌋\) 表示向下取整,称为 Floor,指的是小于或等于 \(x\) 的最大整数;\(⌈x⌉\) 表示向上取整,称为 Ceil,指的是大于或等于 \(x\) 的最小整数,\(eg:⌊5.6⌋ = 5,⌊-5.6⌋ = -6;⌈5.6⌉ = 6,⌈-5.6⌉ = -5。\)↩︎