纹理合成：Pyramid-Based Texture Analysis/Synthesis_综合

Pyramid-Based Texture Analysis/Synthesis

文章目录

Pyramid-Based Texture Analysis/Synthesis
- Steerable Pyramid
- Histogram Matching

受视觉感知理论的启发，Heeger和Bergen提出一种基于金字塔的图像纹理合成技术，着重与随机纹理的合成(与随机纹理对应的是确定性纹理)，此技术的优势是只需要给定目标纹理，不需要更多额外信息即可自动完成目标纹理的合成，关于纹理方面的更多内容可以阅读 Human Texture Perception。本文的基础是拉普拉斯金字塔，关于金字塔的更多内容可以阅读文章几种常见的金字塔变换总结，本算法的输入时给定纹理图像和噪声图像(一般是高斯白噪声)，目的是将噪声图像通过拉普拉斯金字塔转化为与给定纹理在人眼视觉上具有一致性的合成纹理图像。另外由于《Pyramid-Based Texture Analysis/Synthesis》写得还是比较简要的，很多内容都不是很清晰，关于此篇论文的更详细内容可以参考2014年发表的文章《The Heeger-Bergen Pyramid-Based Texture Synthesis Algorithm》。
纹理合成算法可以分为基于邻域的合成方式和基于统计信息的合成方式，基于邻域的合成方式主要是通过对纹理图像的重新排布来完成的，基于统计信息的合成方式主要基于纹理辨别理论( texture discrimination theories)，一般包含一下两个步骤，第一步是统计纹理图像的视觉信息，例如直方图、协方差函数等，第二步是将输入的目标图像或者噪声图像变成与纹理图像统计信息相同的合成图像，本文的纹理合成技术即属于此类方法。

Steerable Pyramid

一般将Steerable Pyramid译作方向可操作金字塔，与拉普拉斯金字塔类似，方向可操作金字塔也是图像的一种多尺度完备表示，在论文The steerable pyramid: A flflexible architecture for multi-scale derivative computation中首次提出，但不同于拉普拉斯的各向同性，方向可操作金字塔具有各向异性，是图像的线性多尺度、多方向分解，其

Algorithm 1: Steerable Pyramid Decomposition

Input：离散数字图像 $u$ ，图像尺寸为 $M \times N$ ；分解尺度数 $P$ ；方向数 $Q$ ；高通滤波器 $h_{0}$ ，低通滤波器 $l_{0}、l$ ，方向带通滤波器 $b_{q}$
Output：图像 $u$ 的方向可控制金字塔由 $P \times Q + 2$ 张不同尺寸的图像组成

1.计算图像的高频部分 $h_{0}*u$

2.计算图像的低频部分 $v←l_{0}*u$

3.计算图像的高频方向带通部分 $b_{q}*v$

4.计算图像的低频部分 $v = l ? v$

5.对图像 $v$ 进行降采样 $v=↓2vv=\downarrow _{2}v$

6.重复步骤3-5共 $P$ 次

纹理合成：Pyramid-Based Texture Analysis/Synthesis

构造Steerable Pyramid示意图：左边为图像分解部分，右边为图像重建部分

方向可操作金字塔可通过FFT来加速计算，记原图像 $u$ 的离散傅里叶变换为：
$(FM,N(u))m,n=u^m,n=∑k=0M?1∑l=0N?1uk,le?2iπ(kmM+lnN),(m,n)∈Ω^M,N\left(\mathcal{F}_{M, N}(u)\right)_{m, n}=\hat{u}_{m, n}=\sum_{k=0}^{M-1} \sum_{l=0}^{N-1} u_{k, l} e^{-2 i \pi\left(k \frac{m}{M}+l \frac{n}{N}\right)}, \quad(m, n) \in \hat{\Omega}_{M, N}$
其对应的离散傅里叶逆变换为：
$(FM,N?1(v))k,l=1MN∑m=?M2∑n=?N2M2?1vm,ne2πi(mkM+nlN),(k,l)∈ΩM,N\left(\mathcal{F}_{M, N}^{-1}(v)\right)_{k, l}=\frac{1}{M N} \sum_{m=-\frac{M}{2}} \sum_{n=-\frac{N}{2}}^{\frac{M}{2}-1} v_{m, n} e^{2 \pi i\left(m \frac{k}{M}+n \frac{l}{N}\right)}, \quad(k, l) \in \Omega_{M, N}$
在极坐标系下定义高通滤波器及低通滤波器：
$r?π2\begin{array}{l} L(r, \theta)=L(r)=\left\{\begin{array}{ll} 1 & \text { if } r \leqslant \frac{\pi}{4}, \\ \cos \left(\frac{\pi}{2} \log _{2}\left(\frac{4 r}{\pi}\right)\right) & \text { if } \frac{\pi}{4} \leqslant r \leqslant \frac{\pi}{2} \\ 0 & \text { if } r \geqslant \frac{\pi}{2} \end{array}\right. \\ H(r, \theta)=H(r)=\left\{\begin{array}{ll} 0 & \text { if } r \leqslant \frac{\pi}{4} \\ \cos \left(\frac{\pi}{2} \log _{2}\left(\frac{2 r}{\pi}\right)\right) & \text { if } \frac{\pi}{4} \leqslant r \leqslant \frac{\pi}{2} \\ 1 & \text { if } r \geqslant \frac{\pi}{2} \end{array}\right. \end{array}$
方向滤波器定义为：
$Gq(r,θ)=Gq(θ)=αQ(cos?(θ?πqQ)Q?11∣θ?πqQ∣?π2+cos?(θ?π(q?Q)Q)Q?11∣θ?π(q?Q)Q∣?π2)；q=1,2.....QG_{q}(r, \theta)= G_{q}(\theta)=\alpha_{Q}\left(\cos \left(\theta-\frac{\pi q}{Q}\right)^{Q-1} \mathbb{1}_{\left|\theta-\frac{\pi q}{Q}\right| \leqslant \frac{\pi}{2}}+\cos \left(\theta-\frac{\pi(q-Q)}{Q}\right)^{Q-1} \mathbb{1}_{\left|\theta-\frac{\pi(q-Q)}{Q}\right| \leqslant \frac{\pi}{2}}\right)；q=1,2.....Q$
其中 $αQ=2Q?1(Q?1)!Q(2(Q?1))!\alpha_{Q}=2^{Q-1} \frac{(Q-1) !}{\sqrt{Q(2(Q-1)) !}}$ 为归一化常数。
最终的高通、低通、方向滤波器如下：
$L0(r,θ)=L0(r)=L(r2)H0(r,θ)=H0(r)=H(r2)Bq(r,θ)=H(r)Gq(θ)\begin{array}{l} L_{0}(r, \theta)=L_{0}(r)=L\left(\frac{r}{2}\right) \\ H_{0}(r, \theta)=H_{0}(r)=H\left(\frac{r}{2}\right) \\ B_{q}(r, \theta)=H(r) G_{q}(\theta) \end{array}$
记 $\in\left\{H_{0}, L_{0}, L, B_{q}\right\}$ 表示各滤波器的连续形式， $\in\left\{h_{0}, l_{0}, l, b_{q}\right\}$ 为各滤波器的离散形式， $\in\left\{h_{0}, l_{0}, l, b_{q}\right\}$ 的离散傅里叶形式表示为 $f^∈{h^0,l^0,l^,b^q}\hat f \in\left\{\hat h_{0}, \hat l_{0}, \hat l, \hat b_{q}\right\}$ ，可由解析函数 $F$ 在频域内采样得到：
$f^m,n=F?ρ(2πmM,2πnN),(m,n)∈Ω^M,N\hat{f}_{m, n}=F \circ \rho\left(\frac{2 \pi m}{M}, \frac{2 \pi n}{N}\right), \quad(m, n) \in \hat{\Omega}_{M, N}$
$\circ \rho$ 表示沿 $x$ 轴采样 $M$ 次，沿 $y$ 轴采样 $N$ 次，其采样过程如下图所示：
纹理合成：Pyramid-Based Texture Analysis/Synthesis

离散空域的采样过程

此时可以将Algorithm 1中 $f ? u$ 可以通过 $f . ? u$ 来实现，下采样可用通过提取离散傅里叶形式的中心部分来实现。

接下来我们需要明确的是如何通过上述过程中分解得到的金字塔重建原图，方向可控金字塔相比拉普拉斯之类的金字塔最大的不同是该金字塔被设计为自转置(Self-Inverting)的，即是如果用 $A$ 表示金字塔中某高频层的分解算法，那么此高频层的重建算子为 $A^{T}$ ，即是 $A^{T}Au=u$ ，具体内容可以参考论文《The steerable pyramid: a flexible architecture for multi-scale derivative computation》，也可以参考小波逐渐分解构建图像金字塔，其重建过程如下：

Algorithm 2：Steerable Pyramid Reconstruction

Input：金字塔尺度 $P$ ，方向数 $Q$ ，原图像尺寸 $M \times N$

Output：原图像 $u$

1.将Algorithm 1中的最后剩下的低通部分 $v$ 作为初始化的 $u$

2.将 $u$ 上采样， $u=↑2uu=\uparrow _{2}u$

3. $u = u ? l$

4.当前scale下所有方向子带与对应方向的高通滤波器 $b_{q}$ 卷子并叠加到 $u$ 上

5.重复2-4直至叠加回全部带方向的高频信息，得到 $u$

6. $u=u*l_{0}+h_{0}$

Histogram Matching

关于直方图均衡化及直方图匹配的详细推导可以参考Rafael C. Gonzalez的《Digital Image Processing》一书，简单来说，对于一幅图像而言当其概率密度函数 $p_{r}(r)$ (即直方图)为均匀分布时，图像的信息熵最高，认为图像包含的信息量最大，而一幅自然图像的概率密度函数往往是非均匀的，假设存在某种映射关系 $s = T (r)$ 使得变换后的概率密度函数 $p_{s}(s)$ 服从均匀分布，此映射 $T$ 即为直方图均衡化，如果 $p_{r}、T_{r}$ 已知，且在规定值域上 $T (r)$ 连续可微，则变换后的变量 $s$ 的概率密度函数可以表示为：
$ps(s)=pr(r)∣drds∣p_{s}(s)=p_{r}(r)|\frac{dr}{ds}|$
又变量 $s 、 r$ 均单调递增，即有：
$ps(s)=pr(r)drdsp_{s}(s)=p_{r}(r)\frac{dr}{ds}$
又我们的目标是变量 $s$ 的概率密度函数为均匀分布，即是：
$ps(s)=pr(r)drds=1Np_{s}(s)=p_{r}(r)\frac{dr}{ds}=\frac{1}{N}$
其中 $N$ 为图像灰度级数目，即有：
$N p_{r}(r)dr = ds$
两边同时积分，得到映射函数 $T (r)$ ：
$T(r)=N∫prdrT(r)=N\int p_{r}dr$
即直方图均衡化的映射函数由原图像的累计分布函数决定。直方图均衡化的目的是将图像的信息熵变成和均匀分布的随机变量的信息熵相同，如果我们并不想让变换后的图像服从均匀分布，而是服从其他特定分布时，同样的只需要给出变换后想要的概率密度函数，即可，此过程称为直方图规定化或者直方图匹配，具体算法如下：

Algorithm：Histogram Matching

Input：输入图像 $u$ ，参考图像 $v$ ，图像尺寸均为 $M \times N$

Output：与图像 $v$ 直方图相同的图像 $u$

1. $L = M \times N$

2.统计图像 $v$ 的直方图，并计算其累计分布直方图，以灰度级为0~255为例，此时的累计分布直方图为一个长度为256的数组，将此数组乘以255之后取整

3.用图像 $u$ 的像素值为索引值从累计分布直方图中查询得到图像 $u$ 的新的像素值

纹理合成：Pyramid-Based Texture Analysis/Synthesis

直方图匹配示意图

本文中重要的概念是将直方图匹配用于方向可控金字塔内的各个频带，以此达到纹理合成的目的，如下图所示：
纹理合成：Pyramid-Based Texture Analysis/Synthesis

金字塔图像的直方图匹配的示意图

可以看出基于直方图匹配的纹理合成算法主要分为一下四步来完成：

给定目标纹理图像 $u$ 和用高斯白噪声初始化的合成图像 $v$
分别将纹理图像和合成图像分解为方向可控金字塔
统计图像 $v$ 的各层各方向带通图像的累计分布直方图，以从直方图对图像 $v$ 所对应的图像做直方图匹配，并将匹配得到的新图像存在图像 $v$ 的金字塔的对应位置
用直方图图匹配得到的新的方向可控金字塔重建图像，即得合成的纹理图像

另外还可以通过迭代以上步骤来优化合成的纹理图像，即用前面合成的纹理图像替换第一步中用高斯白噪声初始化的合成图像 $v$ ，因此合成纹理图像的最终结果将由迭代次数 $N$ 、金字塔层数 $P$ 、方向数 $Q$ 共同决定的，一般来说迭代次数、金字塔层数、方向数都是越大越好，但在实际应用过程中我们往往还需要考虑到效率问题，因此这三个变量的组合方式也是值得注意的。
纹理合成：Pyramid-Based Texture Analysis/Synthesis

迭代次数对合成结果的影响：原图；迭代1次；迭代5次；迭代10次

纹理合成：Pyramid-Based Texture Analysis/Synthesis

分解方向数对合成结果的影响：原图；方向数为2；方向数为4；方向数为8

纹理合成：Pyramid-Based Texture Analysis/Synthesis

分解层数对合成结果的影响：原图；1层；4层；8层