📘 第 10 章：图像和视频处理

这一篇在干嘛？

二维信号处理实战：图像格式与基本操作、HDL 边缘检测、中值滤波、VGA 系统、视频运动检测协处理器。 原书代码为 VHDL，本篇所有代码已改写为 Verilog。

• 图像和视频处理
• 10.1 图像和视频处理概述
• 10.2 案例研究 1: HDL 中的边缘检测
• 10.3 案例研究 2: 使用图像处理库进行中值滤波
• 10.4 案例研究 3: 视频处理中的运动检测由自定义协处理器改进
• 练习

---

自测一下

RGB 与 YUV 各是什么？为什么要转 YUV？

RGB 是三原色分量；YUV 是亮度+色度。人眼对亮度敏感、色度不敏感，转 YUV 后可对色度降采样省带宽，压缩/传输标准都基于 YUV。

3×3 中值滤波为什么比均值滤波好？

中值滤波取中间值，能完全去除椒盐噪声而保住边缘；均值滤波会把噪声扩散成模糊斑。中值排序网络的实现用比较器树，FPGA 友好。

Sobel 边缘检测的原理？

用水平/垂直两个 3×3 差分核分别卷积，得到梯度分量 Gx/Gy，幅值超过阈值即判定边缘。二维卷积用行缓存+乘加阵列实现。

运动检测的基本思路？

相邻帧做差分（或背景建模），变化超过阈值判为运动区域。视频流处理需要帧缓存与流水线设计，常用协处理器+软核分工：硬件做像素级运算，软核做决策。

图像和视频处理

如果说一维 DSP 处理的是随时间变化的波形，那么图像处理就是把信号”升了一个维度”——多出一个空间坐标 $y$；视频处理则再升一维，加入时间。好消息是，前面学过的滤波、FFT 等算法大部分都能直接推广到二维、三维；但也有一些内容是纯一维世界里完全不出现的：

• 去除脉冲噪声的非线性滤波器（中值、最大/最小滤波等）；
• 用于发现视觉不连续点的边缘检测，目的是压缩输入特征、保留最有用的信息；
• 视频中的运动向量，用来描述一个数据块如何从上一帧移动到下一帧，是视频压缩的核心工具；
• 形态学图像处理，这类操作在一维 DSP 中根本没有对应物。

为什么这本书要专门用 FPGA 讲图像处理？关键在硬件条件的变化：现代 FPGA 配备了大容量片内存储块（例如 Stratix II 有 512 Mbit 的片上存储），图像处理算法中大量的行缓存、帧缓存可以直接放在片内，实现效率大幅提升。放在过去没有快速片上存储的年代，工程师必须花费大量开发时间去搭建与外部存储芯片（DDRAM、QDR、RLDRAM、SDRAM 等）的接口，门槛高得多。

实现图像处理算法有几种典型路线：1) 直接用 HDL 设计；2) 使用 IP 库；3) 嵌入式微处理器加定制协处理器（μP）。本章的案例研究会覆盖这些方法。初学者建议先在 MATLAB 的图像处理和计算机视觉工具箱中练习，UNIX 系统上还有集成了许多算法的图形程序 XV 可以参考。

在进入二维、三维细节之前，先弄清楚一维处理、图像和视频在格式与基本处理上的差异。

10.1 图像和视频处理概述

一维处理与图像处理有一个根本性的态度差别：一维滤波强调因果性——输出只能依赖过去时刻的输入，否则硬件不可实现；而图像处理面对的是一整幅已经存好的数据，我们并不希望每次滤波都把图像”翻转”一遍。于是图像处理中的惯例是：

• 使用”类因果”滤波：用滤波器输出替换中心像素；
• 只用 $x$、$y$ 方向都为奇数长度的滤波器，如 $3\times 3$、$5\times 5$、$7\times 7$，这样才存在真正的”中心像素”；
• 滤波器用**核（掩模）**定义而不是卷积运算。对尺寸 $L=2N+1$ 的核，输出为

$$J(r,c)=f\ast I=\sum _{k=-N}^N\sum _{l=-N}^{N}f(k,l)I(r+k,c+l)$$

注意这里没有做卷积要求的翻转——如果真按卷积定义来，矩阵运算前信号要先翻转 180°，图像处理中通常省略这一步。

边界效应

二维中比一维更突出的另一个问题是边界效应：一个 $L\times L$ 的掩模在 $R\times C$ 的图像上滑动时，掩模一旦出界就无法完整计算，损失的像素数为

$$\frac{RL}{2}\cdot 2+\frac{CL}{2}\cdot 2-4{\left(\frac{L}{2}\right)}^2=L(R+C)-L^2$$

计算示例：拉普拉斯-高斯（LoG）滤波器取 $\sigma =3$ 时核尺寸 $L=33$，作用在 $256\times 256$ 的图像上，丢失像素数为 $33\times (256+256)-33^2=33\times 512-1089=16896-1089=15807$，占整幅图的 $15807/65536\approx 24\%$——将近四分之一的边界面积没了！这个损失相当可观。

就保持边界而言，IIR 滤波器理论上表现更好（不需要掩模完整覆盖），但由于非线性相位和稳定性问题，图像处理中很少使用。更常见的做法是把图像向外扩展半个核宽，扩展方式有四种：

• 常数值填充，通常取 $x=0$；
• 对称扩展：…3,2,1 | 1,2,3,…（镜像）；
• 复制边界值：…1,1,1 | 1,2,3,…；
• 圆形（环绕）扩展：…,8,9,10 | 1,2,3…,8,9,10 |。

行方向和列方向都要做同样的扩展。

硬件中的折中：一行偏移

真正做硬件时，图像数据通常以一维流的方式进出，我们常采用类似”圆形扩展、只偏移一行”的近似：

… 1 2	3 4 5	10 11 12 …
3 4 5	10 11 12	13 14 …

也就是第一行的末尾值（3、4、5）直接接到第二行的开头参与计算，第二行的开头值（10、11、12）与第一行的末尾一并处理。这种”跨行借用”只有在假设图像边界是背景、且整幅图背景相似时才合理——大多数自然图像满足这一点，所以工程上很实用。

10.1.1 图像格式

数字图像中最基本的单元叫像素，常用 8 位表示（医学成像可达 12 位）。编码约定是：0 代表黑，最大值（8 位格式即 255）代表白。

彩色怎么来？全彩色谱只用三种基色——红、绿、蓝（RGB）就能合成；从灰度图转 RGB 时三个通道赋相同的灰度值即可。但人眼有个特点：对亮度变化比对颜色变化敏感得多，所以实用的编码方案把灰度（亮度 Y）和差分色度分开传输——亮度用更高的精度/分辨率，色度可以适当压缩，图像质量几乎不降。RGB 中三通道位数相同，没有利用这一特性，于是在电视、JPEG、MPEG 中流行两种分离方案：YIQ 和 YUV（也叫 YCrCb）。

NTSC 电视系统采用的 YIQ 编码变换为：

$$\text{亮度：}Y=0.2989R+0.5870G+0.1140B$$ $$\text{红-蓝绿：}I=0.5959R-0.2744G-0.3216B\text{\hspace{1em}(10.1)}$$ $$\text{品红-绿：}Q=0.2115R-0.5229G+0.3114B$$

ITU BT.601-5 建议的 YCrCb 方案则是：

$$Y=0.299R+0.587G+0.1140B$$ $$C_r=0.713(R-Y)=0.500R-0.419G-0.081B\text{\hspace{1em}(10.2)}$$ $$C_b=0.564(B-Y)=0.169R-0.331G+0.500B$$

考虑到 8 位量化下亮度范围限定在 16～220（留工作余量）、色度零电平在 128，ITU 给出整数化映射：

$$Y=\frac{77}{256}R+\frac{150}{256}G+\frac{29}{256}B$$ $$C_r=\frac{131}{256}R-\frac{110}{256}G+\frac{21}{256}B+128\text{\hspace{1em}(10.3)}$$ $$C_b=\frac{44}{256}R-\frac{87}{256}G+\frac{131}{256}B+128$$

Altera 视频处理工具箱采用略有修改的近似（亮度范围 16 至 $(0.257+0.504+0.098)\times 255+16\approx 235$，色度最大 $0.439\times 256+128\approx 240$）：

$$Y=0.257R+0.504G+0.098B+16$$ $$C_r=0.439R-0.368G-0.071B+128\text{\hspace{1em}(10.4)}$$ $$C_b=-0.148R-0.291G+0.439B+128$$

分辨率与显示顺序

传统图像显示在阴极射线管（CRT）显示器上：像素从左上角开始，从左到右扫完一行，再回到左边开始下一行。分辨率用”先水平 $x$ 后垂直 $y$“的像素数表示，常见有 $640\times 480$、$800\times 600$、$1024\times 768$、$1280\times 1024$。坐标排列如下：

(0,0)	(0,1)	…	(0,C-1)
(1,0)	(1,1)	…	(1,C-1)
⋮		⋮
(R-1,0)	(R-1,1)	…	(R-1,C-1)

视频处理中，典型图像大小是 $16\times 16$ 宏块的倍数，最流行的是 $352\times 288$ 的通用中间格式（CIF）和 $176\times 144$ 的四分之一 CIF（QCIF）。

彩色图像的每个像素由 R、G、B 三个子像素组成，但并非所有格式都按这个顺序存储——例如流行的原始 BMP 格式从左下角开始，到右上角结束。大多数格式对灰度数据用 8 位、每个 RGB 分量也用 8 位；典型的图像文件以包含行列数、颜色表、压缩方案等信息的文件头开始，随后才是图像数据。有些格式支持无损或有损压缩（如 GIF、JPEG）。只有少数静止图像格式（JPEG、JPEG-2000）由 ISO/IEC/ITU 等国际委员会定义，其余很多由公司定义，知识产权未必完全清晰。MATLAB 的 imread() 函数支持下述所有格式。

BMP：不压缩的原始图像格式，使用 RGB 颜色表、每分量 8 位，也允许用少于 24 位编码（例如只存常用颜色的 8 位颜色表放进文件头；8 位灰度则不需要颜色表）。DE2 演示工程（DE2_Default）用的就是它，在嵌入式 μP 设计中也很流行；DE2 板附带的 BMP2MIF 转换器教程见逻辑教材附录 F。

GIF：包含无损图像压缩，做 PPT 幻灯片是好选择，也常被用来做小动画（GIF 动画）。

PNG：另一种无损压缩方案，常用于网页。PNG 采用差分编码，但不使用 GIF 的 Lempel-Ziv-Welch（LZW）压缩算法，因此不需要 LZW 的许可证费用。

格式	Int.STD	无损压缩	有损压缩	MATLAB
BMP	-	-	-	√
GIF	-	√	LZW	√
PNG	-	√	DPCM	√
JPEG	√	√	DCT	√
JPEG-2K	√	√	DWT	√

JPEG：唯一官方静止图像标准之一

JPEG 标准（即 ITU-T 建议 T.81）可在 CCITT 网页免费下载。注意一个细节：标准只规定了解码器，没有限制编码器，给实现改进留了空间。图 10-1 是基本 DCT 有损压缩的主要编码与解码步骤。

图 10-1 JPEG 压缩方案

流程是这样的：8 位输入图像（假定灰度；彩色图像通常先用式 10.1 的 YIQ 或式 10.2 的 YCrCb 变换）被划分为互不重叠的 $8\times 8$ 块，做二维 DCT。因为二维 DCT-II 在 $x$、$y$ 方向可分离，实现时先对八行做一维变换，再对八列做变换。DC 分量（块内均值）位于 $8\times 8$ 数据块左上角。之后用不同的表量化亮度和色度 DCT 系数——标准第 143 页建议的亮度量化表见图 10-2(a)。

16	11	10	16	24	40	51	61
12	12	14	19	26	58	60	55
14	13	16	24	40	57	69	56
14	17	22	29	51	87	80	62
18	22	37	56	68	109	103	77
24	35	55	64	81	104	113	92
49	64	78	87	103	121	120	101
72	92	95	98	112	100	103	199

(a) 量化表

(b) AC 系数的 Zig-zag 编码
图 10-2 JPEG 压缩方案细节

表中的值是解码端重构 DCT 值时要乘的系数；编码端则先除以这些值再四舍五入到整数——把昂贵的除法负担放在编码器一侧。从量化表可以看出：左上（低频）量化因子小，用更多比特表示；右下（高频）量化因子大，8 位输入的高频系数通常只剩 0 或 ±1。随后按”之字形”顺序排列系数，编码器对第一个非零系数对编码，用块结束（EOB）符号终止整个块，见图 10-2(b)。

由于相邻块的 DC 值（即 $8\times 8$ 块内所有像素的平均值）通常相近，DC 先做差分编码（$\mathrm{DIFF}=D{C}_i-D{C}_{i-1}$），再做霍夫曼编码：先用 4 位编码”大小”，再编码差值本身：

4位SSSS	DIFF 值
0	0
1	-1, 1
2	-3, -2, 2, 3
3	-7,…,-4, 4,…,7

其余 63 个 AC 系数用可变游程长度编码（VLC），格式为（零游程长度，下一个值），各用 4 位；游程长于 15 时插入特殊符号 (15,0)，EOB 用 (0,0) 编码。

计算示例：DCT 编码系数

$$(10,3,0,0,-1,1,0,0,0,0,2,\mathrm{EOB})\text{\hspace{1em}(10.5)}$$

先依据前一个块的 DC 值译出 DC，AC 系数则译码为：

$$(0,3)(2,-1)(0,1)(4,2)(0,0)\text{\hspace{1em}(10.6)}$$

其中 (0,3) 表示”零游程 0、数值 3”，(2,-1) 表示”跳过 2 个零后是 -1”，依次类推，(0,0) 即 EOB。DCT 值还可以再用霍夫曼编码进一步压缩；标准附录 K 给出了需写入 JPEG 文件头的 DC/AC 亮度和色度示例表。编码表通常是中等偏高压缩率的——这也是很多文件的实际压缩比不如预期的原因之一。

图 10-3 给出了几个 JPEG 编码结果：低压缩比时 JPEG 与 BMP 几乎看不出差别；高压缩时则能看到 DCT 分块带来的块效应。

(a) 原始 BMP 图像（677KB）

(b) Q20 的 JPEG 低压缩图像（53KB）

(c) Q1 的 JPEG 高压缩图像（22KB）
图 10-3 JPEG 压缩结果

JPEG-2000：从小波出发的新工具箱

JPEG-2000 可以看作原始 JPEG 的更新版，但有几处实质性变化。最核心的是：$8\times 8$ DCT 块处理被 DWT（离散小波变换） 取代——有损压缩通常用浮点 Daubechies 9/7 小波滤波器，无损压缩用 5/3 整数滤波器。初始 DWT 之后进行量化与算术编码，见图 10-4。

图 10-4 JPEG-2000 压缩方案

JPEG-2000 是 ISO、IEC/ITU 多个小组的联合成果，遗憾的是标准规范文档并不免费（本书写作时售价 96 瑞士法郎），不过有许多出版物把基本思想讲得很清楚。新标准不是单一算法，而是一个”工具箱”，其中最有趣的特点有：

• 低比特率下性能卓越——DWT 天然避免了经典 JPEG 的块效应；
• 同时支持无损与有损压缩；
• 感兴趣区域（ROI）编码：允许对图像中特别重要的部分（例如说话者的面部相对背景）用更高质量编码；
• 通过同步标记、固定长度编码器和重复报头增强抗比特错误鲁棒性；
• 通过水印、标签、冲压和加密实现图像安全。

代价是复杂度更高。JPEG-2000 通常以比 JPEG 更高的计算复杂度换取更好的 PSNR、MSE 或 SSIM。如果应用需要高质量或低比特率，或需要安全性等特性，就选 JPEG-2000。

10.1.2 基本图像处理操作

图像/视频工具箱里最常用的基本操作包括：

• 改变图像尺寸，增减处理的像素数；
• 裁剪出感兴趣区域（ROI）；
• 旋转图像；
• 对一幅或多幅图像做加、减、乘、除和平均；
• 彩色↔灰度转换；
• 以 BMP、GIF、JPEG 等格式读取和存储图像。

点操作

基于像素自身值进行变换的叫点操作。典型用途包括调整整体过亮/过暗的图像，或依据局部/全局像素统计（直方图）进行调整。CRT 显示器必备的点操作是 γ 校正，对 8 位数据：

$$y=255(i/255{)}^{1/\gamma }\text{\hspace{1em}(10.7)}$$

$0<\gamma <1$ 时是指数型曲线，图像变暗；$\gamma >1$ 时是对数型曲线，图像变亮。例如 $\gamma =10$ 的曲线在低输入区抬升缓慢、在高输入区快速上升，暗部被进一步压暗；$\gamma =0.2$ 则相反，暗部被大幅提亮。NTSC 电视系统采用 $\gamma =2.2$。

(a) $\gamma =10$

(b) $\gamma =1$

(c) $\gamma =0.2$
图 10-5 γ 校正示例

滤波操作

图像滤波可线性可非线性，作用对象不只是单个像素，还包括相邻像素，多数用于图像增强。大多数线性方法也可以在变换域（FFT、DCT）中求解。几个典型例子：

• 低通滤波器消除随机噪声；
• 差分滤波器检测边缘；
• 最大值滤波器去除胡椒噪声（因为黑色编码为 0，取邻域最大值就把黑点抹掉了）；
• 最小值滤波器去除盐噪声（因为白色编码为 255，取最小值把白点抹掉）；
• 中值滤波器同时抑制盐和胡椒噪声；
• 去抖滤波器，消除相机移动时反褶积滤波造成的模糊；
• 自适应维纳滤波器去噪。

图像压缩（静止图像与视频）是活跃的研究领域，许多方法利用了人眼的物理特性——对低频变化比对高频变化更敏感。这类变换编码通常再结合可变长度编码、霍夫曼编码或算术编码等（无损）熵编码。图像/视频处理中最流行的二维变换是 FFT、DCT 和 DWT，Haar 变换、Hartley 变换等二维形式也有使用。

形态学操作

形态操作用于增强信息或减小模式识别系统的输入空间，基本动作是”在邻域中添加像素”（增厚，即膨胀/dilation）和”删除邻域中的像素”（稀疏，即侵蚀/erosion）。膨胀产生更多像素，使线条变粗、填拢小孔；侵蚀使线条变细，并去掉”单个像素”这类小噪声区域。两者通常组合使用：先侵蚀后膨胀称为开运算（opening），先膨胀后侵蚀称为闭运算（closing）。开、闭都能平滑物体轮廓，闭运算还能填充小孔和间隙。

这些操作可定义在灰度图像和黑白（BW）图像上，BW 更常见。灰度图的形态操作通常就是邻域内的最大/最小操作；BW 图则定义邻域掩模以何种规则替换像素。

图 10-6 展示了用 MATLAB 函数 bwmorph 和 bwperim 的若干形态操作：开运算能去掉手写文本下的背景横线（图 10-6(b)）；骨架化是一系列侵蚀与开运算的组合，为每个对象生成中心线，让所有对象的结构像 X 光片里的骨骼一样显现（图 10-6(c)）；边界提取可以先对图像用掩模 M 侵蚀，再用图像减法计算 $BW=I-(I\Theta M)$，其中 $\Theta$ 表示侵蚀，结果得到边界/周边图（图 10-6(d)）。此外，“多数黑色”操作在邻域中多数像素为黑时把当前像素置黑，作用类似膨胀，对边缘检测也有帮助。

图 10-6 形态操作示例

10.2 案例研究 1：HDL 中的边缘检测

图像里对象的大多数信息都藏在边缘中。找边缘并不难：不同物体、或物体与背景具有不同的强度，我们只需找到像素值中的不连续点并标记它们。检测对比度变化有两条路：

1. 计算一阶导数并寻找（绝对）最大值：

$$f^{\prime }(x)\approx \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}\text{\hspace{1em}(10.8)}$$

2. 计算二阶导数 $f^{\prime \prime }(x)=d{f}^{\prime }(x)/dx$ 并寻找过零点。

图像处理中用掩模/核代替连续求导。一阶导数用掩模 $M_1=[-1,1]$；二阶导数用 $M_2=[1,-2,1]$。先用 $M_1$ 分析一个短测试序列，验证这样的短掩模确实管用：

I	10	10	10	40	40	30	20	10	0	0	30	0
$I\ast M_2$	-	0	0	30	0	-10	-10	-10	-10	0	30	-30
$\Vert I\ast M_2\Vert$	-	0	0	30	0	10	10	10	10	0	30	30

测试序列包含一个上升沿、一段下降斜坡和结尾处的单个”噪声”像素。结果显示：边缘和斜坡都被检测到了，但同时也暴露了短掩模对噪声的敏感性——那个孤立噪声像素产生了和真实边缘同样强的响应（30）。

再看与二阶导数掩模 $[1,-2,1]$ 的相关：

I	10	10	10	40	40	40	30	20	10	0	0	0
$I\ast M_2$	-	0	30	-30	0	-10	0	0	0	10	0	-
过零	-	-	√	√	-	√	-	-	-	√	-

可见在边缘斜率变化处出现 ± 跳变，随后是过零点。但检测过零比在一阶导数中找最大值要费事，所以我们的实现选择一阶导数 + 最大值检测。

从一维到二维：Prewitt 与 Sobel

为避免图像处理中的位置偏移，采用中心化掩模，如 $M=[-1,0,1]$。因子 1/2 来自两方面考虑：8 位运算中避免溢出的最坏情况增长，以及期望 ${\sum }_k|f_k|=2$。为了降低单像素噪声的影响，通常再在若干行上取平均，于是得到 $3\times 3$ 核——分别检测垂直和水平线条的 Prewitt 掩模 $M_x$、$M_y$：

$$M_x=\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]M_y=\left[\begin{array}{ccc}1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\end{array}\right]\text{\hspace{1em}(10.9)}$$

Sobel 建议的掩模强调中心行/列的系数：

$$M_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]M_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]\text{\hspace{1em}(10.10)}$$

关于 Sobel 掩模的缩放因子：悲观取法是 $1/{\sum }_k|f_k|=1/8$，完全避免位增长；或者按”近似真实导数”取 1/4（对应理想边缘 …0,1,1,…）。但边缘检测只关心最大值出现在哪里，只要数据不溢出，我们干脆跳过缩放。图像 $I$ 与两个掩模分别卷积得到 $G_x=I\ast M_x$、$G_y=I\ast M_y$，总梯度是矢量合成，但我们只关心它的长度：

$$|G|=\sqrt{|G_x{|}^2+|G_y{|}^2}\text{\hspace{1em}(10.11)}$$

开平方太麻烦，通常用幅度近似。许多论文用的标准近似 $|G|=|G_x|+|G_y|$ 误差界限达 41%（因为 $\sin (x)+\cos (x)=\sqrt{2}\sin (x+\pi /4)$，最大偏差 $\sqrt{2}-1\approx 0.41$）。我们改用 $L_1$ 近似

$$|G|=\max (|G_x|,|G_y|)+\min (|G_x|,|G_y|)/4$$

它把误差容限缩小了 3 倍，而硬件代价只是几次比较和一次移位（除以 4）。

图 10-7 用埃菲尔铁塔模型演示了 Sobel 算子的效果：$|G_x|$ 突出水平梯度（检测垂直边缘/细节），$|G_y|$ 突出垂直梯度（检测水平边缘/细节），求和后两者兼有，最后加阈值把边缘与非边缘分类。阈值的选择取决于图像的信噪比（S/N），对高 S/N 图像，推荐取测量到的最大梯度的约 10%～20%。

(a) 原图
(b) $|G_y|$

(c) $|G_x|$

(d) $|G_x|+|G_y|$

(e) 阈值
图 10-7 使用 Sobel 掩膜的示例：(a) 原图，(b) $M_y$ 滤波器给出了水平边缘/细节，(c) $M_x$ 滤波器给出了垂直边缘/细节，(d) 梯度和，(e) 阈值转换

二阶导数：LoG 与更鲁棒的 Canny

二阶导数方法通常采用高斯-拉普拉斯（LoG）滤波器。高斯函数易于微分、易于复现，而且能同时完成边缘检测和平滑。所需方程为：

$$g(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-x^2/(2{\sigma }^2)}\text{\hspace{1em}(10.12)}$$ $$g^{\prime }(x)=\frac{-x}{{\sigma }^2}g(x)\text{\hspace{1em}(10.13)}$$ $$g^{\prime \prime }(x)=\left(\frac{x^2}{{\sigma }^4}-\frac{1}{{\sigma }^2}\right)g(x)\text{\hspace{1em}(10.14)}$$

图像处理中，LoG 常用**高斯差（DoG）**来近似：

$$g^{\prime \prime }(x)\approx c_1{e}^{-x^2/(2{\sigma }_1^2)}-c_2{e}^{-x^2/(2{\sigma }_2^2)}\text{\hspace{1em}(10.15)}$$

由此得到常用的离散掩模：

$$\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]\left[\begin{array}{ccccc}0& 0& -1& 0& 0\\ 0& -1& -2& -1& 0\\ -1& -2& 16& -2& -1\\ 0& -1& -2& -1& 0\\ 0& 0& -1& 0& 0\end{array}\right]\text{\hspace{1em}(10.16)}$$

除这些短掩模之外，还有在低信噪比下更鲁棒的复杂方法，例如 Canny 方法：先用平滑滤波器去噪，再在较大掩模内用精细分类方案只保留区域内的一条边缘——思想上类似形态学中的稀化操作，目标是让一条边线或斜坡只显示一条线。图 10-8 的对比显示 Canny 算子检出的细节明显多于 LoG 或 Sobel。

(a) 原始图像

(b) Sobel 算子

(c) 高斯拉普拉斯算子

(d) Canny 检测器

图 10-8 边缘检测方法的比较

10.2.1 二维 HDL 滤波器设计

一维滤波器已经有很多种架构，多加一个维度会让架构选择更多。不过在讨论架构之前，先回顾我们实现一维 FIR 时用过的三板斧：

• 尽可能用转置形式实现滤波器；
• 在乘法器块（MCM）中组合尽可能多的系数；
• 最小化存储器 I/O：利用抽头延迟线（TDL）复用已读取的数据，每个数据只读一次。

第三个概念要求为 $N\times N$ 掩模准备 $N-1$ 行的行缓冲器（硬件上用 FIFO 实现）。行缓冲器必须用嵌入式存储块实现，否则会消耗大量 LE 来搭建 FIFO。这样就能构建图 10-9(a) 那种带单个乘法器块的 FIR 滤波器。但单 MCM 架构有两个缺点：其一，第一行 FIR 的输出又是行缓冲器的输入，其位宽通常比输入数据宽，因此 FIFO 行缓冲器需要更宽的数据总线；其二，多个滤波器要复用行缓冲器——边缘检测本质上就是两个滤波器，如果稍加重排，两个滤波器可以共用相同的行缓冲器，但单个 MCM 块会变成 $N$ 个。图 10-9(b) 就是这种 $N$-MCM 结构（以 3×3 核为例），图中 $z^{-C}$ 表示跨越整幅图像所有（剩余）列的延迟。

(a) 单 MCM 滤波器块架构

(b) N-MCM 滤波器块架构
图 10-9 二维滤波方法的比较，$z^{-C}$ 是整个图像的所有（剩余）列上的延迟

可分离滤波器

如果二维核恰好可以写成两个一维滤波器的外积（即秩为 1 的矩阵），滤波就简化为”先按行滤、再按列滤”。从前面的核来看：Prewitt 和 Sobel 是可分离的，而 3×3 LoG（式 10.16）不是。以 Sobel 的 $M_x$ 为例（式 10.10）：

$$M_x=\left[\begin{array}{c}1\\ 2\\ 1\end{array}\right]\times \left[\begin{array}{ccc}-1& 0& 1\end{array}\right]=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]\text{\hspace{1em}(10.17)}$$

于是 $M_x$ 滤波可以先用 $[1,2,1{]}^{\mathrm{T}}$ 垂直掩模、再用 $[-1,0,1]$ 水平掩模分两步完成。LoG 则麻烦一些：虽然可以按式 10.12 计算 $g^{\prime \prime }(x)$ 再构建外积，但那不是真正的二维 LoG。用 MATLAB 函数 fspecial('log', width, sigma) 可以方便地生成 LoG，但真正的二维 LoG 通常秩不为 1，无法分解。

架构上（见图 10-10），垂直滤波器不再拥有单个 MCM 块，但水平滤波器仍可以是单个 MCM 块，两个滤波器还共用同一组行缓冲器。

(a) 行和列滤波器
(b) 两个滤波器架构使用相同的行缓冲器
图 10-10 二维可分离滤波器架构

到底省了多少运算量？长度为 $L$ 的一维滤波器每个输入样本需要 $L$ 次乘法和 $L-1$ 次加法；若滤波器对称（假设 $L$ 为奇数），乘法可减到 $(L+1)/2$ 次。表 10-2 对比了 $L\times L$ 掩模两种情况下的操作数。

计算示例：取 $L=3$。不可分离且对称时：乘法 $L(L+1)/2=3\times 4/2=6$ 次，加法 $L^2-1=8$ 次；可分离且对称时：乘法 $2(L+1)/2=4$ 次，加法 $2(L-1)=4$ 次——乘加量都近乎减半，$L$ 越大节省越显著。

操作	不可分	不可分-对称	可分	可分-对称
乘法	$L^2$	$L(L+1)/2$	$2L$	$2(L+1)/2$
加法	$L^2-1$	$L^2-1$	$2(L-1)$	$2(L-1)$

对边缘检测应用还有一点很关键：$M_x$ 和 $M_y$ 两个滤波器可以共用行缓冲器，图 10-10(b) 正是这种双滤波器设计的架构。

10.2.2 图像系统设计

算法选定后，下一步是硬件系统设计。如今先进的开发板大多带多种图像显示接口，可以直接利用。例如 Xilinx Atlys 开发板有 2×2 高清多媒体接口（HDMI）可用于音视频输出；Altera DE2-115 开发板则带电视解码器、VGA 和 DVI 等接口。DVI 常用作 LCD 监视器接口，经典 VGA 标准常作为 PC 中 CRT 显示器上的默认启动接口。VGA 的基本分辨率是 640 列 × 480 行、刷新率 60Hz，绝大多数显示器（包括液晶）都支持这一模式——所以我们直接用 DE2-115 接 CRT 显示器即可，无须购买任何新显示硬件。

系统基本设置见图 10-11，FPGA 提供必要的控制和 RGB 数据信号。注意 DE2-115 板使用 8 位 RGB 数据（DE2 老板子是 10 位）——处理 10 位数据太繁琐，而且大多数图像本来就是 8 位 RGB 分辨率。DE2-115 上的 VGA 芯片是 Analog Devices 的 ADV7123KSTZ140，内含三路高速 10 位 DAC，像素速率最高可达 140MHz，因此 VGA(640×480)、SVGA(800×600)、XGA(1024×768)、SXGA(1280×1024) 各种显示模式都能用。表 10-3 列出了各显示模式所需的像素时钟（总行×总列×帧速率）以及灰度图像的存储需求（彩色需要三倍；即 行×列×8 bit）。

图 10-11 边缘检测系统整体设计图

模式	刷新率 Hz	分辨率 列×行	Total C×R	像素时钟 MHz	大小 Mbits
VGA	60	640×480	800×525	25.175	2.4
VGA	85	640×480	832×509	36	2.4
SVGA	60	800×600	1056×628	40	3.84
SVGA	75	800×600	1056×625	49	3.84
SVGA	85	800×600	1048×631	56	3.84
XGA	60	1024×768	1344×806	65	6.29
XGA	75	1024×768	1312×806	75	6.29
XGA	85	1024×768	1376×808	95	6.29
SXGA	60	1280×1024	1688×1066	108	10.49

VGA 时序怎么算？ VGA 信号的总时间 = 有效视频数据 + 同步开销，因为 CRT 电子束需要时间从行尾回到行首、从帧尾回到帧首。假定使用外同步（而不是把同步嵌入绿色通道）。以 60Hz 的 VGA 为例，一行从低电平有效的同步脉冲开始：

• 水平同步：3.8μs，即 25MHz 像素时钟下的 96 个时钟周期；
• 水平后沿：1.9μs，即 48 个周期；
• 有效视频：640 个周期，即 25.4μs；
• 前沿：0.6μs，即 16 个周期——之后光束回到行首。

总计 $96+48+640+16=800$ 个时钟周期，正好对应表 10-3 中的 Total C×R = 800×525。垂直方向方法类似：2 行垂直同步脉冲，随后 33 行垂直后沿，接着 480 行有效视频，最后 10 行供光束返回第一行，合计 $2+33+480+10=525$ 行。图 10-12 给出了显示开始前的 ModelSim 仿真波形，可以看到两行垂直同步和 33 行的垂直后沿。

图 10-12 VGA ModemSim 仿真系统整体结构

10.2.3 VGA 边缘检测系统的组装

基于以上研究，我们在 DE2-115 板上组装边缘检测实验装置，设计目标如下：

• 边缘检测器使用 $3\times 3$ Sobel 掩模（见图 10-10）；
• 幅度使用 $L_1$ 近似（见表 2-14）；
• 图像以 MIF 文件格式存储，编译时可替换；
• 利用可分离的行、列滤波器设计滤波器，且两个滤波器共用同一行缓冲器（见图 10-10(b)）；
• 显示屏为 60Hz 刷新率的 $640\times 480$ 灰度 VGA（见表 10-3）；
• 使用数据循环回绕；
• 滑块开关作为输入：低 8 位用于数据，高 3 位用于选择显示模式。SW 的两个 MSB 选择主图像源：00=原图，$01=G_x$，$10=G_y$，$11=|G|$；第三个 MSB 打开/关闭阈值；
• 七段显示器显示（见图 10-13(a)）：0,1 位显示 $\max (G_x)$；2,3 位显示 $\max (G_y)$；4,5 位显示 $\max (|G|)$；6,7 位显示帧计数器；
• LCD 显示项目信息：第 1 行”边缘检测”，第 2 行”DSP 与 FPGA”。

(a)

(b)
图 10-13 (a) LCD 文本和七段显示器的显示，(b) 图像存储器的配置

系统搭建要点

初步设计可以从 DE2_115_demonstration 文件夹中的 DE2_115_Default 工程出发——它包含的功能比我们需要的还多，而且很多同学对开发板出厂配置的这些启动功能已经很熟悉。下面看系统整体架构、所需控制信号和 Verilog 实现。

ADI 公司的 ADV7123KSTZ140 VGA 芯片需要 clk 输入、RGB 数据、blank_n 和 sync_n；同时 FPGA 电路要把水平/垂直同步信号直接送给 VGA 监视器。所需的 25MHz 时钟可以由 PLL 产生，也可以把 50MHz 输入二分频得到——后者实现简单、资源省，还简化了仿真。因为只显示灰度图像，三个通道数据相同，即 R=G=B=灰度数据。SYNC_N 信号可以不用，但 BLANK_N 必须用起来：当 CRT 从帧尾回扫到行首/帧首时，RGB 信号要清零，或者说用 blank_n 告诉 VGA 监视器何时是有效显示数据——这样也简化了整体数据管理和滤波器设计。

由于只有 8 位数据，Sobel 滤波和幅度计算都可能造成位增长，必须监视。做法是用七段显示器连续计算 $G_x$、$G_y$ 和 $|G_{xy}|$ 的最大值，两个最高有效数字位显示帧计数器——60Hz 帧率下 8 位计数器约 $256/60\approx 4$ 秒绕回一圈。图 10-11 概述了整个系统。

图片通过 Altera MIF 文件加载（DE2_115_Def 完成文件 I/O）。DE2_115 CD 的 VGA_DATA 文件夹里有可执行程序 PrintNum.exe，它把 BMP 文件分割成颜色表和包含图像数据的 MIF 文件。也可以自己写 MATLAB 脚本：用 imread() 读图，再写成 MIF 文件，例如下面这段脚本把 $480\times 640$ 的 BMP 转成 MIF：

原书 MATLAB 脚本（非 VHDL，保持原样）：

clear;
%% m file to convert 480x640 BMP file into MIF text file
I = imread('BCoin8state.bmp');   %% 读取图像
figure, imshow(I);               %% 显示图像
[r c] = size(I);                 %% 获取行数与列数
rc = r*c; J = [];
for k = 1:r                      %% 重排为一维数组
    J = [J I(k,:)];
end
disp(sprintf('Saving %d pixel for picture.mif', rc));
fid = fopen('picture.mif','w');
fprintf(fid,'depth=%d;\r\nwidth = 8;\r\n', rc);
fprintf(fid,'address_radix = dec;\r\n');
fprintf(fid,'data_radix = dec;\r\n');
fprintf(fid,'content \r\nbegin\r\n');
for k = 1:rc
    fprintf(fid,'%d:%d;\r\n', k-1, J(k));
end
fprintf(fid,'end;\r\n');
fclose('all');

除了 MIF 文件，可能还需要生成一些测试台数据。例如把一个 $3\times 3$ 的测试数据包放在图像左上角，需要修改地址 0–5、640–645 和 1280–1285 处的 MIF 数据（每行 640 个像素，所以第 2 行从 640 开始、第 3 行从 1280 开始）：

0:10;   1:20;   2:30;   3:0;  4:0;  5:0; ...
640:40; 641:50; 642:60; 643:0; 644:0; 645:0; ...
1280:70; 1281:80; 1282:90; 1283:0; 1284:0; 1285:0; ...

顶层 Verilog 实现

原书 VHDL 改写为 Verilog。下面是顶层模块的端口声明与主要寄存器/线网（对照 DE2-115 的引脚定义）：

module sobel_vga_top (
    ///// CLOCK /////////
    input            CLOCK_50,
    ///// LED /////////
    output [8:0]     LEDG,
    output [17:0]    LEDR,
    ///// KEY /////////
    input  [3:0]     KEY,
    ///// SW /////////
    input  [17:0]    SW,
    ///// SEG7 /////////
    output [6:0]     HEX0, HEX1, HEX2, HEX3,
    output [6:0]     HEX4, HEX5, HEX6, HEX7,
    ///// LCD /////////
    output           LCD_BLON,
    inout  [7:0]     LCD_DATA,
    output           LCD_EN, LCD_ON, LCD_RS, LCD_RW,
    ///// VGA /////////
    output [7:0]     VGA_B,
    output           VGA_BLANK_N,
    output           VGA_CLK,
    output [7:0]     VGA_G,
    output           VGA_HS,
    output [7:0]     VGA_R,
    output           VGA_SYNC_N,
    output           VGA_VS,
    ///// Test ports ////////
    output [10:0]    x1_out, x2_out, x3_out,
    output [10:0]    y1_out, y2_out, y3_out,
    output [7:0]     tap1_out
);

寄存器与线网声明。注意两个行缓冲器 tap1、tap2 各占一行（640 字节），正好是 $N-1=2$ 行缓冲中最靠近当前行的两行：

    // ====================
    // REG/WIRE declarations
    // ====================
    parameter Hactive = 640;        // 每行有效像素数
    reg [7:0] tap1 [1:Hactive];     // 行缓冲器 1（上一行）
    reg [7:0] tap2 [1:Hactive];     // 行缓冲器 2（上上行）
    integer i;
 
    wire [31:0] mSEG7_DIG;
    reg  [31:0] Count;
 
    // VGA controller signals
    reg         VGA_CTRL_CLK;       // 25MHz：50MHz 二分频
    wire [9:0]  mVGA_R, mVGA_G, mVGA_B;
    wire [19:0] mVGA_ADDR;
    wire        DLY_RST;
    wire        mVGA_CLK;
    wire        VGA_Read;           // VGA data request
    wire        mDVAL;
 
    // 来自 VGA 控制器的信号
    wire [7:0]  LCD_D_1;
    wire        LCD_RW_1, LCD_EN_1, LCD_RS_1;
    wire        iRST_n;
    reg         oBLANK_n, oHS, oVS;
    wire [7:0]  b_data, g_data, r_data;
 
    ///// 辅助信号
    reg  [18:0] ADDR;
    reg  [7:0]  gray_data;
    wire [10:0] dx, asum, sum, ssum;
    wire [8:0]  Rxy;
    wire [7:0]  Gxy;
    reg  [7:0]  Gx, Gy;
    reg  [10:0] y1, y2, y3;
    reg  [10:0] x1, x2, x3;
    reg  [7:0]  XY_max, Y_max, X_max;
    wire        VGA_CLK_n;
    wire [7:0]  index;
    wire        cBLANK_n, cHS, cVS, rst;
    wire [10:0] h;                  // 水平计数
    wire [9:0]  v;                  // 垂直计数

结构化代码部分——同步信号、七段译码、复位延迟和 25MHz 时钟生成：

    // ====================
    // Structural coding
    // ====================
    assign VGA_SYNC_N = 1'b0;       // 未使用
    assign LCD_DATA = LCD_D_1;
    assign LCD_RW   = LCD_RW_1;
    assign LCD_EN   = LCD_EN_1;
    assign LCD_RS   = LCD_RS_1;
    assign LCD_ON   = 1'b1;
    assign LCD_BLON = 1'b0;         // 不支持
 
    // 帧计数器：以垂直同步为节拍
    always @(posedge oVS or negedge KEY[0]) begin
        if (!KEY[0]) Count <= 0;
        else         Count <= Count + 1;
    end
    assign mSEG7_DIG = {Count[7:0], XY_max, Y_max, X_max};
 
    assign tap1_out = tap1[Hactive];
    assign tap2_out = tap2[Hactive];
    assign LEDR[15:0] = {v[7:0], h[7:0]};
    assign LEDG[7:0]  = index;      // 存储器数据
 
    // 七段译码查表
    SEG7_LUT_8 u0 (
        .oSEG0(HEX0), .oSEG1(HEX1), .oSEG2(HEX2), .oSEG3(HEX3),
        .oSEG4(HEX4), .oSEG5(HEX5), .oSEG6(HEX6), .oSEG7(HEX7),
        .iDIG(mSEG7_DIG));
 
    // 复位延迟定时器
    Reset_Delay r0 (
        .iCLK(VGA_CTRL_CLK),
        .oRESET(DLY_RST));
 
    // T 触发器二分频产生 25MHz 像素时钟
    always @(posedge CLOCK_50 or negedge KEY[0]) begin
        if (!KEY[0]) VGA_CTRL_CLK <= 0;
        else         VGA_CTRL_CLK <= ~VGA_CTRL_CLK;
    end
 
    LCD_TEST u5 (
        .iCLK(VGA_CTRL_CLK),
        .iRST_N(DLY_RST),
        .LCD_DATA(LCD_D_1), .LCD_RW(LCD_RW_1),
        .LCD_EN(LCD_EN_1),  .LCD_RS(LCD_RS_1));

VGA 同步发生器与地址产生器。地址在每次帧起始（HS、VS 同时有效）时清零，实现数据循环回绕：

    // VGA 信号发生器
    assign VGA_CLK = VGA_CTRL_CLK;
    assign rst     = ~iRST_n;
    assign iRST_n  = DLY_RST;
 
    video_sync_generator LTM_ins (
        .vga_clk(VGA_CTRL_CLK),  // 25MHz VGA 时钟
        .reset(rst),             // 复位
        .blank_n(cBLANK_n),      // 数据（不）有效
        .HS(cHS), .VS(cVS),      // 行/场同步
        .h(h), .v(v));
 
    ///// 地址产生器：blank_n 有效时逐像素递增，帧首清零
    always @(posedge VGA_CTRL_CLK or negedge iRST_n) begin
        if (!iRST_n)
            ADDR <= 19'd0;
        else if (cHS == 1'b0 && cVS == 1'b0)
            ADDR <= 19'd0;
        else if (cBLANK_n == 1'b1)
            ADDR <= ADDR + 1;
    end
 
    ///// 读取图像数据（MIF 初始化的片上 ROM）
    assign VGA_CLK_n = ~VGA_CLK;    // 反相时钟使数据在正确相位取出
    img_data img_data_inst (
        .address(ADDR),
        .clock(VGA_CLK_n),
        .q(index));

最后是核心的 Sobel 滤波与幅度计算，全部在一个时钟进程内完成。幅度用 $L_1$ 近似 $\max (|G_x|,|G_y|)+\min (|G_x|,|G_y|)/4$（右移 2 位实现除 4），超过 255 时饱和到最大值：

    ///// 幅度：max(|Gx|,|Gy|) + min(|Gx|,|Gy|)/4
    assign Rxy = (Gx > Gy) ? ({1'b0, Gx} + {3'h0, Gy[7:2]})
                           : ({1'b0, Gy} + {3'h0, Gx[7:2]});
    assign Gxy = (Rxy > 8'hFF) ? 8'hFF : Rxy[7:0];   // 饱和限幅
 
    // My 滤波：水平方向 [1 2 1] 卷积（含除 4 缩放）
    assign sum  = y1 + 2*y2 + y3;                    // 无寄存器，组合逻辑
    assign asum = (sum[10] == 1'b1) ? -sum : sum;    // 取绝对值
    assign ssum = asum >> 2;                         // 乘以 1/4 缩放因子
 
    // Mx 滤波：垂直方向 [-1 0 1] 差分（绝对值 + 除 4）
    assign dx = (x1 > x3) ? ((x1 - x3) >> 2) : ((x3 - x1) >> 2);
 
    assign x1_out = x1;  assign x2_out = x2;  assign x3_out = x3;
    assign y1_out = y1;  assign y2_out = y2;  assign y3_out = y3;
 
    always @(posedge VGA_CTRL_CLK) begin  /////// Sobel 滤波器
        if (cBLANK_n == 1'b1) begin
            // 显示 X、Y 和总梯度的最大值（位增长监视）
            if (Gx  > X_max)  X_max  <= Gx;
            if (Gy  > Y_max)  Y_max  <= Gy;
            if (Gxy > XY_max) XY_max <= Gxy;
 
            // ---- abs(Gy) 滤波 ----
            // 水平滤波器 [1 2 1]
            if (ssum > 8'hFF)
                Gy <= 8'hFF;                   // 饱和
            else
                Gy <= ssum[7:0];               // 已含除 4
            y3 <= y2;
            y2 <= y1;
            y1 <= {8'h0, index} - {8'h0, tap2[Hactive]}; // My 垂直滤波 [1 0 -1]
 
            // ---- abs(Gx) 滤波 ----
            // 垂直差分 [-1 0 1]，已含阈值缩放
            if (dx > 8'hFF)
                Gx <= 8'hFF;
            else
                Gx <= dx[7:0];
            x3 <= x2;
            x2 <= x1;
            // Mx 垂直滤波 [1 2 1]
            x1 <= {8'h0, index} + {8'h0, tap1[Hactive], 1'b0}
                              + {8'h0, tap2[Hactive]};
 
            ///// 两条抽头延迟线被 Mx 和 My 共用
            for (i = Hactive; i > 1; i = i - 1)
                tap1[i] <= tap1[i-1];          // 抽头延迟线 1：移位
            tap1[1] <= index;                  // 新像素进入第 1 级
 
            for (i = Hactive; i > 1; i = i - 1)
                tap2[i] <= tap2[i-1];          // 抽头延迟线 2：移位
            tap2[1] <= tap1[Hactive];          // 承接上一行末尾数据
        end
    end
 
endmodule

这段代码正好体现了 10.2.1 节的架构思想：tap1、tap2 两条抽头延迟线（行缓冲）被 $M_x$、$M_y$ 两个滤波器共用；垂直方向用 $[1,2,1]$（x1 支路，其中 tap1[Hactive],1'b0 就是左移一位实现 ×2）和 $[1,0,-1]$（y1 支路）两个一维核，水平方向的差分 $[-1,0,1]$ 由 dx 的减法完成——这正是可分离滤波器的”先行后列”两步实现，配合 cBLANK_n 门控，同步与数据管理都变得干净利落。

VGA 输出与灰度选择逻辑的收尾

在边缘检测案例的最后部分，设计要把 Sobel 滤波的结果送到 VGA 显示器上。这一段逻辑要做三件事：根据拨码开关选择显示哪种数据、把同步信号延迟一拍以保证与像素数据对齐、最后把选出的灰度数据同时接到 R、G、B 三个颜色通道上（因为显示的是灰度图，三个通道必须相等）。原书 VHDL 改写为 Verilog 后如下：

///// 在时钟下降沿锁存有效数据
always @(posedge VGA_CLK_n) begin
    case (SW[17:15])
    3'b000 : gray_data <= index;
    3'b001 : gray_data <= (index < SW[7:0]) ? 8'h00 : 8'hFF;
    3'b010 : gray_data <= Gx;
    3'b011 : gray_data <= (Gx < SW[7:0]) ? 8'h00 : 8'hFF;
    3'b100 : gray_data <= Gy;
    3'b101 : gray_data <= (Gy < SW[7:0]) ? 8'h00 : 8'hFF;
    3'b110 : gray_data <= Gxy;
    default: gray_data <= (Gxy < SW[7:0]) ? 8'h00 : 8'hFF;
    endcase
end
 
///// 把同步信号延迟一个时钟周期
always @(negedge VGA_CLK) begin
    oHS      <= cHS;
    oVS      <= cVS;
    oBLANK_n <= cBLANK_n;
end
 
///// VGA 控制器输出
assign VGA_BLANK_N = oBLANK_n;
assign VGA_HS      = oHS;
assign VGA_VS      = oVS;
assign VGA_B = gray_data;
assign VGA_G = gray_data;
assign VGA_R = gray_data;
endmodule

从初学者角度看，这段代码的关键在 case(SW[17:15])：三位开关共有 8 种组合，分别选择”原始图像数据""带阈值处理的原始数据""$G_x$""带阈值的 $G_x$""$G_y$""带阈值的 $G_y$""边缘幅值 $G_{xy}$""带阈值的 $G_{xy}$“。阈值运算 (< SW[7:0]) ? 0 : FF 的意义是二值化——低于开关设定值的像素显示纯黑，高于的显示纯白。如果没有这一步，8 位灰度差值中的算术回绕（下溢）会显示成”伪亮纹”，破坏边缘的观感。

整个顶层文件的组成可以这样理解：

• Verilog 文件从 I/O 端口声明开始，其后是一些仅供仿真观察用的附加测试端口，再往内是内部信号。为避免端口数量爆炸，可以把原来的 vga_controller 模块”提升”到顶层中直接例化。
• 结构化编码先把固定的 I/O 信号用 assign 分配好；oVS 的 always 块同时兼任一个帧计数器，其计数值显示在七段数码管的第 6、7 段上。
• 接着是七段译码器和一个复位延迟定时器（上电序列偶尔失败时用它救急）。
• 50 MHz 板载时钟经触发器二分频得到 25 MHz 的 VGA 像素时钟，它同时也是 Sobel 滤波器的主时钟。
• 然后例化 LCD 字符显示组件（TERASIC 组件允许为应用指定自定义文本）。LCD_TEXT.v 内含 HD44780 控制器的模式 ROM 十六进制编码，使用字母数字编码：数字 0 从 $30_{16}$ 开始，大写字母从 $41_{16}$ 开始，小写字母从 $61_{16}$ 开始，空格是 $20_{16}$。
• 再例化 VGA 同步信号组件：当需要新图像数据时，它提供行/场同步、行列计数器和消隐信号 cBlank_n。图像以一维线性数组存放，cBlank_n 用来递增存储器地址计数器 ADDR。图像存储器用 MegaWizard 生成，容量为 $2^{19}\times 8$ 位，初始化文件为 VGA_DATA/BCoin8state.mif；用 bmp2txt.m 脚本可以生成其他图像的 MIF 文件，注意图片必须是 $640\times 480$ 像素的 BMP。

Sobel 滤波器的编码要点

Sobel 滤波器的编码顺序值得仔细体会。不需要寄存器的组合逻辑部分放在由 iVGA_CLK 和 cBLANK_n 控制的 always 语句之外——这是避免无意中推断出多余寄存器的常用手法。always 块内的前三条 if 语句在七个数码管上追踪显示 $G_x$、$G_y$ 和 $G_{xy}$ 迄今出现的最大值，方便人工标定阈值。

滤波器本体按两级级联：对 $G_y$ 先做垂直滤波 $[1,0,-1{]}^T$，再做水平滤波 $[1,2,1]$；对 $G_x$ 顺序相反，先做 $[1,2,1{]}^T$ 垂直滤波，再做 $[1,0,-1]$。绝对值和除 4 缩放的代码也放在 always 语句之前，同样是为了避免推断寄存器。滤波输出经阈值运算防止算术回绕。

最后一部分是两个长度为 640 的移位寄存器——行缓冲区。它们应当综合成块 RAM，否则会消耗大量逻辑单元。特别注意：不能给这对寄存器文件加复位，因为块 RAM 硬件上根本没有复位端口，加了复位反而会阻止工具把它映射到块 RAM。

滤波结果通过 case(SW[17:15]) 连接到灰度输出端口：SW15 有效时执行阈值操作，可在存储器原始数据、$G_x$ 输出、$G_y$ 输出以及幅值 $G_{xy}$ 四个信号之间选择。根据最终行/场同步的时序精度，HS、VS 和 BLANK_n 信号各延迟一拍，最后一组 assign 把选中的灰度数据接到 RGB 端口。

仿真时序与测试台分析

图 10-12 给出了全局仿真的粗略概况。由于要处理完 33 行之后才出现第一个”有用”数据，仿真时间相对较长：显示一行需要 $35\mu s$ 即 875 个时钟周期，33 行就是约 28K 个周期。

图 10-14 用于 ModelSim 边缘检测模拟的测试台数据处理

图 10-14 展示了整个 $3\times 3$ Sobel 掩模数据到齐后测试台的细节。可以看到 50 MHz FPGA 时钟输入和 25 MHz 二分频信号；KEY1 表示复位；SW=110… 与代码中的 case 语句比较后，使 Gxy 输出到 VGA 监视器，即显示不带阈值的边缘矢量幅值。名为 Memory 的数据块中，ADDR 是地址计数器，data 是从存储器读出的数据，tap1_out 和 tap2_out 是 $G_x$、$G_y$ 滤波器行缓冲器的输出。

测试台向 8 分之一（8 quarter）硬币图像人工注入如下状态数据：

$$data=\left[\begin{array}{ccc}70& 80& 90\\ 40& 50& 60\\ 10& 20& 30\end{array}\right]\text{\hspace{1em}(10.18)}$$

数据之后跟随三个零值。

VGA 数据块展示显示侧的信息：当前处理的行列计数器、水平/垂直同步与 blank_n 有效信号、以及发送到显示器的 RGB 信号（灰度图所以三路相等）。最后的 Filter 数据块揭示 Sobel 滤波器的内部工作：x1 是 $G_x$ 的行滤波结果，即由 $[1,2,1{]}^T$ 加权的数据，前三个值为：

$$x1(0)=70+2\times 40+10=160$$ $$x1(1)=80+2\times 50+20=200$$ $$x1(2)=90+2\times 60+30=240$$

随后输出 x2 和 x3（体现水平滤波器 $[1,0,-1{]}^T$ 的延迟需求），再经过滤波、绝对值和除 4 得到 $G_x$，例如：

$$G_x(k)=|x1(k)-x3(k)|=|240-160|/4=20$$ $$G_x(k+1)=|x1(k+1)-x3(k+1)|=|200-0|/4=50$$ $$G_x(k+2)=|x1(k+2)-x3(k+2)|=|240-0|/4=60$$

除以 4 是为了防止后续求和时溢出。$G_y$ 侧的计算方法类似：先算 $[1,0,-1{]}^T$ 的行滤波：

$$y1(0)=70-10=60$$ $$y1(1)=80-20=60$$ $$y1(2)=90-30=60$$

再用绝对值和 $[1,2,1]$ 水平滤波并除以 4 压缩：

$$G_y(k)=|y1(k)+2y2(k)+y3(k)|=|60+2\times 60+60|/4=60$$ $$G_y(k+1)=|y1(k+1)+2y2(k+1)+y3(k+1)|=|60+2\times 60+0|/4=45$$ $$G_y(k+2)=|y1(k+2)+2y2(k+2)+y3(k+2)|=|60|/4=15$$

最后用 $\max (|G_x|,|G_y|)+\min (|G_x|,|G_y|)/4$ 合成边缘幅值 gray_data（即 $G_{xy}=VGA\_B$）：

$$G_{xy}(k)=\max (20,60)+\min (20,60)/4=60+20/4=65$$ $$G_{xy}(k+1)=\max (50,45)+\min (50,45)/4=50+45/4=61$$ $$G_{xy}(k+2)=\max (60,15)+\min (60,15)/4=60+15/4=63$$

测试台最后三个数据显示目前为止的绝对最大值——本测试数据中 Gxy 最大值的增量为 65。

还有一点关于仿真方式的选择：通常做功能仿真，但本设计存储器初始化文件太大，很快就会耗尽 PC 内存，因此在这个特定项目中改用时序仿真模型，其对存储器的表示更紧凑。

仿真确认无误后即可生成编程文件，把 VGA 显示器接到 DE2-115 板并下载设计。然后可以在水平和垂直方向上尝试不同的边缘渐变，根据七段数码管显示的最大值，为最终二值化表示找到合适的阈值。

该项目的综合结果在 Compilation Report 中：设计使用 816 个 LE，未使用嵌入式乘法器，占用 302 个（70%）M9K 存储器块；对缓慢 85C 时序模型，TimeQuest 报告 $F_{max}=44.26\mathrm{MHz}$。

10.3 案例研究 2：使用图像处理库进行中值滤波

用 HDL 设计图像/视频处理系统往往耗时很长，因为数据量的增加会让系统复杂度迅速攀升。但每个图像/视频系统都要执行的那几个基本任务其实差别不大——一套好的基本构建模块对剩下的设计大有裨益（10.1.2节讨论过这些基本操作）。Altera 提供两条路：其一是专业工具箱，帮助设计滤波器（如二维线性和中值滤波器）以及颜色空间转换、行缓冲器编译、伽马校正等图像内部操作；其二是付费许可评估这些模块，但源代码和脚本已加密，无法研究或修改。大学版本的程序则支持包含大量基本构建模块的图像/视频处理设计，提供 VHDL 和 Verilog 源代码以及 TCL 脚本，可导入 SOPC 或 Qsys 设计环境，用法见用户手册。随着新版 Quartus II 发布，常规内容可随之更新。这些模块也可以在独立系统中使用并搭建自己的 Avalon 总线系统，但更常见的做法是把它们作为 Nios II 微处理器系统的协处理器。注意第 9 章提过 Altera 已宣布停止 SOPC 环境，今后只支持 Qsys，所以新系统应当在 Qsys 中设计。图像处理模块需从 Altera 大学程序网页下载，安装在 ip→University Program→Audio Video→Video 下。表 10-4 简要列出了可用模块；所有模块都带 Avalon 总线接口，因此名字都以 altera_up_avalon_video 开头。详细用法和 C 编码示例见文档《Altera DE 系列板的 Video IP 核》。

表 10-4 Altera 大学程序工具箱中的图像处理块

IP	描述
alpha_blender	组合两个视频流
bayer_resampler	将 2×2 Bayer 模式 [G1,R;B,G2] 格式转换为 24 位 RGB 格式 [R,(G1+G2)/2,B]
character_buffer_with_dma	将 ASCII 字符转换为图形显示
chroma_resampler	在 YCrCb 格式之间转换，包括复制、插入和丢弃像素
clipper	修改视频流分辨率，即添加或删除行/列
csc	颜色空间转换器 YCrCb ↔ RGB
decoder	从复合视频端口或 CCD 摄像机读取视频
dma_controller	存储和检索内存中的视频帧；可指定寻址模式、帧分辨率和像素格式
dual_clock_buffer	在两个时钟域之间传输数据流
edge_detection	边缘检测四步法：高斯平滑、Sobel 滤波、极大值抑制和滞后阈值
pixel_buffer_dma	从外部存储器读取视频帧，支持连续或 X-Y 寻址
rgb_resampler	在不同 RGB 位宽间切换：8、9、16、24、30 位
scaler	按整数因子更改宽度和高度分辨率
stream_router	拆分和合并视频流
test_pattern	生成 24 位 RGB 测试图像：色相沿 x 轴从 0 到 360，饱和度沿 y 轴从 0 到 1
vga_controller	生成 VGA DAC 的同步信号
vip_bridges	将视频转换为 Altera VIP 格式

在展示我们的系统细节之前，先看一个典型的适合微处理器的应用——这恰好是 HDL 系统设计的一个小小的挑战。

10.3.1 中值滤波器

中值滤波器已被证明对被”椒盐”（S&P）噪声破坏的图像非常有效。椒盐噪声可能来自 CCD 阵列的缺陷像素、数字成像器的增益异常或数据传输错误。去除这种黑白噪声的思路是评估每个像素的邻域：如果某像素的黑白值与邻域不符，就用邻域的中值替换它——这是最可能的校正值。之所以不用均值滤波器，是因为均值滤波有抹掉边缘等细节、造成图像模糊的趋势。中值滤波器可归纳为三步：

• 定义中值计算使用的邻域大小，例如 $5\times 5$，或仅水平 $1\times 5$、仅垂直 $5\times 1$；
• 在足够大的邻域内计算像素的中值（即排序后取中间值）；
• 用中值替换该像素。

邻域大小取决于信噪比和可用的计算能力。但要注意，S/N 并不是表征中值滤波器性能的最佳指标：几乎无椒盐噪声的图像经中值滤波后 S/N 反而更低，因为未被污染的像素也被邻域值替换了。中值计算本身不复杂——只需按值对窗口内数据排序，取有序列表的中间值。但排序在硬件或软件中都不是轻松的操作。在讨论典型软硬件方案之前，先看几个图像示例。信噪比高时用小窗口（如 $3\times 3$）较好；信噪比低时窗口要大，常用 $9\times 9$、$11\times 11$、$13\times 13$。中值滤波显然不可分离，但可以尝试行/列分开处理的方法，代价是输出性能会有差异。在 MATLAB 中可以这样表达：

[I1,map] = imread('TigerGray.bmp');
d = .1;
I2 = imnoise(I1,'salt & pepper', d);
I3 = medfilt2(I2,[5,5]);
figure, imshow(I3)

该脚本向加载的图像添加 5% 椒盐噪声，用 $5\times 5$ 中值滤波并显示结果。图 10-15 是典型的滤波示例。图 10-15(c) 表明 $5\times 5$ 中值滤波并不合适：噪声没被清除，图像反而更模糊了。先用 $1\times 5$ 中值滤波、再用 $5\times 1$ 滤波的效果优于直接 $5\times 5$。想彻底去除椒盐噪声时可以用更大的窗口，例如图 10-15 中的 $11\times 11$；但大窗口会删除一些细节，比如图像背景壁纸中的花纹，而这些细节用小窗口就能保留。窗口大小问题的另一种解法是自适应中值滤波器。

(a) 原始图像

(b) 受随机椒盐噪声干扰的图像

(c) 尺寸为 $5\times 5$ 的均值滤波

(d) 先 $1\times 5$ 中值滤波，接着列滤波 $5\times 1$（e）$5\times 5$ 中值滤波

图 10-15 中值和平均滤波示例

(f) $11\times 11$ 中值滤波

10.3.2 HDL 中的中值滤波器

在硬件中设计中值滤波是相当有挑战性的任务，高效排序方法已有众多实现尝试。FPGA 供应商通常提供专门的 IP 模块和应用笔记，覆盖排序设计这类棘手内容。下面简要讨论设计的关键问题。

排序网络。 由 Batcher 首次提出：在每个阶段，两两相邻的寄存数据通过最小/最大网络。对 $2N$ 个数据采用 $2\times 2$ 比较方式，总共需要 $N(N+1)/2$ 次比较。图 10-16(a) 显示了 $N=9$ 数据的排序网络。可以并行实现 $N/2$ 个比较单元，直到排序步骤数减少到 $N+1$。

(a) $N=9$ 的 Batcher 排序网络（小块是寄存器，灰色块不需要用于中值计算）

(b) 行列三重排序器方法

图 10-16 (a) $N=9$ 的排序网络；(b) 行列三重排序器方法

对流行的 $3\times 3$ 窗口，有一个取巧方案：先做水平方向再做垂直方向的长度 3 中值滤波，最后对对角线元素再做一次，即可得到正确的中值，见图 10-16(b)。这个方案对 $3\times 3$ 滤波器非常有吸引力，但无法推广到更大的中值滤波器。

位投票法。 这是一种与分布式算术滤波器思路完全不同的方法：不看数据对，而是看所有数据的单个位。从 MSB 位开始，寻找中值的迭代过程如下：统计所有数据 MSB 中 1 的个数。设有 $N$ 个待排序数据，若 1 的计数大于 $N/2$，显然中值落在以 1 开头的那部分数据中；若计数小于 $N/2$（即 0 更多），中值落在以 0 开头的数据中。下一次迭代只看仍然”有资格”的数据——把不再参与中值竞争的数据打上标志，然后查看剩余数据 MSB 的下一位，再次计数并按多数位筛选。此过程一直持续到 LSB。图 10-17 显示了位投票方法的硬件架构，每个像素使用一个处理单元。

图 10-17 位投票方法的硬件架构，其中每个像素使用一个处理元素

累积直方图法。 第三类方法对数据块使用累积直方图。对 $B$ 位数据，直方图用 $2^B$ 个计数位表示。对每个传入数据，把所有等于或大于该像素值对应的计数位加一。最后把累计到 $N/2$ 的值作为中值。实现上可以给每个计数块加 $N/2$ 个比较器，再用优先级译码器选出高于阈值的最小计数箱。若想实现滑动窗口中值，必须从直方图中减去”旧”值。图 10-18 显示了整体架构。

图 10-18 累积直方图架构

10.3.3 Nios 中值滤波图像处理系统

设计新系统通常有两种选择：从零开始，往 Qsys 库里逐步添加组件——这适合只有几个组件的小系统；对于要在开发板上用大量组件的大系统，使用板卡供应商或 FPGA 工具供应商提供的启动模板效率更高。对 DE2 板有出厂配置，TERASIC 也提供了几个可按需定制的完整系统实例。本书使用 Altera 大学程序媒体计算机 VHDL 系统的 Nios Qsys 版作为设计起点，因为 Altera 宣布未来软件版本中将用 Qsys 取代 SOPC Builder。

媒体处理器中有几个块本例并不需要，例如 Expansion_JP5、PS2_Port、PS2_Port_Dual、Serial_Port 和 Audio。媒体处理器使用 QVGA 320×240 的降色分辨率，但我们更愿意用 640×480 的 VGA 同时显示四幅灰度图像，所以也从 Qsys 中删除了 VGA_Pixel_Scaler。下面是需要保留的组件清单，以及支持 640×480 灰度图像处理所需的修改：

• CPU：标准 Nios 处理器，硬件支持乘法，4K 指令高速缓存，2 级 JTAG 调试器（见第 9 章）；
• sysid：提供唯一的系统 ID 号，简化 Eclipse 在 Nios II 系统中定位 USB 端口；
• merged_resets：为内部、外部和 27 MHz 时钟提供复位信号；
• clk、sys_clk、vga_clk、clk_27 和 External_Clocks：为所有模块提供所需时钟；
• SDRAM：保存数据和程序存储器；
• SRAM：用作像素缓冲区；
• Red_LEDs：并行输出端口，驱动 18 个红色 LED；
• Green_LEDs：并行输出端口，驱动 9 个绿色 LED；
• HEX3_HEX0 和 HEX7_HEX4：两个 32 位输出端口，各驱动八个七段数码管；
• Slider_switches：18 位拨码开关的并行输入端口；
• Pushbuttons：四个按钮；
• JTAG_UART：Avalon JTAG 通用异步收发器，各带 64 字节 I/O 缓冲；
• Interval_Timer：32 位计数器、125 ms 周期的 Avalon 定时器；
• AV_Config：设置音频和视频配置；
• VGA_Pixel_Buffer：指定 X-Y 寻址模式、640×480 分辨率，并设置 8 位灰度空间；
• VGA_Pixel_RGB_Resampler：把 8 位灰度转换为 30 位 RGB 信号；
• VGA_Char_Buffer：在 VGA 上显示字符；
• Alpha_Blending：把图像和字符显示组合成单个视频流，字符显示作为前景；
• VGA_Dual_Clock_FIFO：放在 alpha 混合器与 VGA 控制器之间同步数据流；
• VGA_Controller：生成 VGA 显示的全部同步信号；
• Char_LCD_16x2：在两行 LCD 上显示系统信息；
• CPU_fpoint：硬件加速四个常用浮点运算——加、减、乘、除。注意这个块来自 Nios 自定义指令，位于 ip→altera→nios2_ip 下。

Qsys 中的整体系统配置如图 10-19 所示。可以看到大多数组件来自 Altera 大学程序 IP 模块库。最后对比一下与原始媒体计算机相比节省的资源和性能：由于大学版媒体计算机还要支持音频处理，所需资源更多。从表 10-5 可见 LE 和存储器需求下降、速度提升。

图 10-19 Qsys 中的整体系统配置

表 10-5 多媒体和图像处理结果对比

合成结果	媒体计算机	图像处理器	提升
Fmax 性能	92.35 MHz	132.17 MHz	43%
LE	13718	11825	16%
乘法器	11	11	0%
PLL	2	1	100%
CPU	9 个 M9K	9 个 M9K	—
JTAG	6 个 M9K	2 个 M9K	—
VGA	12 个 M9K	10 个 M9K	—
浮点	15 个 M9K	15 个 M9K	—
音频	4 个 M9K	—	—
总 M9K	46 个 M9K	36 个 M9K	27%

10.3.4 SW 中的中值滤波器

如果算法评估已经完成、要进入实时视频应用，就需要前面讨论的硬件实时平台。可以看到，在硬件里搭中值滤波器并不便宜，要花大量开发时间；而实现一个能适配不同长度和方向的中值滤波器更是额外的挑战。由于我们想试验不同噪声水平、滤波器类型、滤波器长度和不同图像，固定系统就不够灵活了。对实验环境而言，软件实现的中值滤波器更灵活——在软件里可以轻易改变滤波器的长度和方向。

先收集一些与 Altera 为 Nios II 提供的开发环境相关的基本事实：有紧凑易用的 Altera 大学监视器程序，也有更先进强大的 Nios II Eclipse 开发系统。图像可以存放在闪存、SRAM（借助存储器初始化文件）或主机中。不过开发图像处理系统最好能利用主机文件系统——它允许使用存放在主机上的图像。这个所谓的主机文件系统只在 Nios II Eclipse 以调试模式运行时可用，标准的文件 I/O 函数（如 fopen、fclose、fputs、fgets）都可以用，因此我们需要 Nios II Eclipse 开发系统。

由于系统包含大量组件，从一个非常小的 SW 项目（比如 hello world）开始是个好主意——让 Eclipse 软件为外部组件生成所需的驱动程序。一般先用低级系统功能直接读写 I/O 端口：

#define switches (volatile char *) 0x01901050
#define leds (char*) 0x01901060
...
s = *switches;
printf("New SW value = %d\n", s);
*leds = s;

使用更复杂的组件时，需要更高层的功能，例如 Altera 的硬件抽象层（HAL）：

/* 把 LCD 光标归位 */
IOWR( LCD_BASE, LCD_WR_COMMAND_REG, 0x02 );

若使用 usleep() 或 fopen() 等 ANSI C 标准函数，则用最高抽象层。

成功运行 hello world 后就可以开始实现图像处理。目标是把屏幕分成四个 $320\times 240$ 像素的区域，分别显示不同的图像处理操作。由于选择了 X-Y 寻址模式，像素索引这样计算：

index = (row<<10) + col;

第一个要实现的例程是在主机上打开文件并按 QVGA 大小读入一个 MIF 文件。用 bmp2txt.m 脚本并加一行 I = imresize(I, .5); 把 $640\times 480$ 图像缩放为 $320\times 240$，把图片加载到图像存储器的左上角。因为读的是 MIF 文本文件，只需检查行首字符是否为数字（ASCII 值在 ‘0’ 到 ‘9’ 之间），再把冒号后面的值转成数字。文件 I/O 的完整代码如下：

printf("File reading begins\n");
pFile = fopen("/mnt/host/Qpicture.mif", "r");
if (pFile == NULL) perror("Error opening file");
else {
  for (;;) {
    for (l = 0; l < 256; l++) Line[l] = ' ';  // 清空行缓冲
    if (fgets(Line, 256, pFile) == 0) break;  // 文件结束
    if ((Line[0] >= '0') && (Line[0] <= '9')) {
      a = Line[0]; u = 0;
      while (a != ':') { u++; a = Line[u]; }  // 定位冒号
      u++; a = Line[0]; v = 0;
      while (a != ';') { v++; a = Line[v]; }  // 定位分号
      d = 0;
      for (j = u; j < v; j++) { aa[d] = Line[j]; d++; }
      aa[d] = '\0';
      // '0' 的 ASCII 码为 0x30 = 48
      w = aa[0] - 48;              // 1 位数字
      if (d > 1) w = w*10 + aa[1] - 48;  // 2 位
      if (d > 2) w = w*10 + aa[2] - 48;  // 3 位
      col = n % 320; row = (int) n / 320;
      pixel_color = w % 256;       // 全部转灰度
      offset = (row << 10) + col;  // 计算半字地址
      *(pixel_buffer + offset) = pixel_color;  // 写像素
      if ((n % 1000) == 0) {
        printf("%d) Line %s select %s:%d\n", n, Line, aa, w);
      }
      n++;  // 像素计数器递增
    }
  }
  printf("\nThe file contains %d pixels.\n", n);
}
printf("File reading ends\n");

由于传输走的是串行 JTAG 电缆，速率不高：一幅 QVGA 图像约需 2 分 32 秒，在一台运行 Windows 7 的 HP 17 主机上的传输速率为 43.9 Kb。

第二步，按照拨码开关的设置给图像添加椒盐噪声，同时在 x 和 y 方向各加两条黑白线，代码很简单。

第三步，水平应用长度为 L 的滤波器，结果显示在左下角。核心代码是数组排序算法。微处理器排序可选归并排序、堆排序、快速排序、基于直方图或经典线性排序。对长度小于 20 的短数组，推荐直接插入排序；对中等长度（20～50），插入时可跳过多个数据的 Shell 方法（又称递减增量排序）更合适；更大的长度推荐 Quicksort。直接插入排序复杂度为 $O(N^2)$，对大 N 慢但紧凑，实现如下：

int k, j;
char t;
for (k = 1; k < len; k++) {
  t = x[j = k + 1];
  while (j != 1 && x[j - 1] > t) {
    x[j] = x[j - 1];
    j--;
  }
  x[j] = t;
}

最后的图像处理步骤是垂直滤波，代码与水平滤波非常类似。

我们在 VGA 显示器上加了”DSP with FPGA”签名，还把最终的 QVGA 数据存回主机文件系统，便于后续处理和记录文档。例如可以把文件加载进 MATLAB，重排数组后显示为图像：

fid = fopen('results.txt','r');   %% 打开文件读取
A = fscanf(fid,'%d',[1,inf])';    %% 读入图像数据
I = reshape(A,[480,640]);         %% 转成行/列格式
figure, imshow(I/255)             %% 显示图像

图 10-20 显示了传回 PC 并在 MATLAB 中显示的主机文件结果：在水平和垂直方向先后使用了长度为 5 的中值滤波器，板上 SW 允许设置添加到图像的噪声电平。

图 10-20 类似于 VGA 显示器的图像处理系统的结果

10.4 案例研究 3：视频处理中的运动检测由自定义协处理器改进

由于带宽和数据量要求极高，视频处理长期由模拟信号处理主导，例如 NTSC、PAL 或 SECAM 系统。图像通常以每秒 25～30 帧（FPS）处理，模拟带宽约 6～8 MHz。隔行扫描（交替扫描偶数行和奇数行）使画面在感官上呈现两倍帧速率。

随着 VLSI 的进步，视频帧的数字处理和存储已成为现实，各种 ITU 和 MPEG 标准都有记载。视频处理领域研究最多的工作是数字视频压缩——要处理的数据量实在惊人。最小的视频格式 QCIF 为 $176\times 144$ 像素、30 FPS、4:2:0 色度子采样，原始数据就达每秒 1.1 MB；60 FPS 的 $1920\times 1080$ 全高清电视原始速率高达 124 MB/s；一部 1 小时的电影原始数据约 446 GB。要把这样的电影装进 DVD 或在合理时间内通过互联网传输，必须做大量数据压缩。

快速处理海量数据难以依靠算术编码这类最先进但最慢的压缩技术，我们需要更快但可能压缩率略低的方案。好消息是视频数据包含大量冗余：考虑视频序列中两个相邻帧，大多数数据保持不变，只有少数像素在特定方向上改变或移动——这可以用在预测方案中。给定时刻 t 的当前帧，计算它与 $t+1$ 时刻下一帧的差值，只传输这个差。可以直接在原始域中做帧预测（类似第 8 章 ADPCM 语音编码的做法），也可以先用 DCT 之类的变换再做预测。图 10-21 展示了这种混合编码系统在时域中的基本架构。

图 10-21 时域中的预测编码

图像处理最理想的变换是 DCT：它接近最优的 KLT，且二维版本可以用行/列方案实现。DCT 在介绍 JPEG 标准时已讨论过，本节聚焦预测的实现。

10.4.1 运动检测

在图 10-21 的方案中，最耗时的环节是计算两个相邻帧之间的运动矢量。做法是选一个大小合适的窗口，在 $\pm d$ 像素的位移范围内比较块如何从一帧移动到下一帧。成本函数可选均方误差（MSE）、互相关函数或平均绝对差（MAD）。这些度量之间差别不大，通常选实现最简单的 MAD：

$$MAD(i,j)=\sum _{x=1}^M\sum _{y=1}^N\left|I_t(x,y)-I_{t-1}(x+i,y+j)\right|\text{\hspace{1em}(10.19)}$$

其中 $-d\le i,j\le d$，$M=N=4$、8 或 16，位移 $d=16$。MATLAB Simulink 的计算机视觉系统工具箱可以实现这一功能。图 10-22 显示了在 Simulink 中实例化的模块：从 QuickTime 电影 Train.mov 加载已缩放为 QCIF 格式的两帧，然后按式(10.19)计算 MAD。运动矢量可用 $V=X+jY$ 坐标或幅值 $|V|$ 显示。

图 10-22 使用 Simulink 模块的 MATLAB MAD 计算过程

测试框架如图 10-23 所示。设置搜索范围（MATLAB 称最大位移）和块大小为 $15\times 15$（注意 MATLAB 只允许奇数长度块），在 x 方向得到 $\lceil 176/15\rceil =12$ 个、y 方向得到 $\lceil 144/15\rceil =10$ 个运动值，如图 10-23(c)～(e) 所示。可以看到火车从左向右移动，y 方向的运动矢量远少于 x 方向。把这些矢量导出到 MATLAB 工作区细看：MAD 值在 y 方向范围为 -11…9，x 方向为 -14…15。

(a) 移动火车的第一帧

(b) 第二帧

(c) 在 x 上的运动矢量

(d) 在 y 上的运动矢量

(e) 运动矢量幅值 $|V|$

图 10-23 取样 QuickTime 电影 train.mov 里两帧的 MATLAB 仿真结果

全搜索 MAD 非常耗时，文献中提出了许多加速运动矢量搜索的方法，表 10-6 给出了概述。一些方法从粗糙网格逐步走向精细网格，例如二维对数法或三步法；另一些沿梯度先在 x 方向移动一个或多个步长，再在 y 方向做同样移动。但请记住：这些方法并不总能找到正确的梯度。随着 VLSI 技术的改进，我们会看到更多全搜索算法的应用——因为它是唯一能保证找到正确运动矢量的算法。

表 10-6 运动估计算法的比较

方法	最大步长	$d=\pm 7$，$n=\lceil {\log }_2(d)\rceil$，min…max	参考文献
全搜索	$(2d+1{)}^2$	225…225	[427]
二维对数 TDL	8n+2	13…26	[428]
三步法 TSS	8n+1	25…25	[429]
修改的二维对数	6n+1	13…19	[430]
一次一次 CDS	2d+3	5…17	[431]
正交搜索 OSA	4n+1	13…13	[432]

10.4.2 ME 协处理器设计

运动矢量的计算可以有多种实现方式。前两个案例研究中我们看到：HDL 通常给出最好的性能和最小的设计，但对数据或算法变化缺乏灵活性；IP 模块方案（中间的案例研究）灵活性最大，允许快速更换数据和算法，代价是面积更大、速度更低。第三个案例研究取中间路线：仍然利用 Nios II 软件设计的灵活性，但这次要找出算法中最慢的部分，通过与 Nios II 处理器紧密耦合的自定义协处理器（自定义指令，CI）改进它们。自定义指令在 Qsys 环境中的开发和实例化见第 9 章 9.5.3 节。

运动估计所需的 C 代码已在先前论文中公布，包括全搜索算法、三步算法、块复制例程（do_dma 和 GetBlock）以及 MAD 和 MSE 成本计算例程。既然要改进全搜索算法，先看看计算量分布在哪里。事实证明：全搜索中 91% 的计算、二维 LOG 和三步搜索中几乎 100% 的计算都花在成本计算上——这是 CI 加速的最佳候选。

不过有一个小问题需要先解决，才能获得实质性改进：数据类型不匹配。图像数据用 8 位 char 表示，而 CI 允许我们用 32 位。因此希望一次 CI 调用能计算四个绝对差。办法是让 char 数组和 int 数据共用同一个起始地址：复制数据时照常逐字节进行，而使用 CI 时改用从同一地址开始的 int 指针。

用如下简单测试（文件 MADtest.c）验证这个技巧工作正常：

int ival[3];
char *cval = (char*) ival;
int i = 0, r;
 
for (i = 0; i < 12; i++) cval[i] = i;  // 写入 char
printf("Char values: ");
for (i = 0; i < 12; i++) printf("%x", cval[i]);
printf("\n");
printf("int value via pointer: ");
for (i = 0; i < 3; i++) printf("%x", ival[i]);
printf("\n");
r = ALT_CI_MAD_0(ival[0], ival[1]);  // 通过 CI 计算
printf("CI: MAD A = %d\n", r);

在 Nios 记录窗口生成如下输出：

Char values: 0 1 2 3 4 5 6 7 8 9 a b
int value via pointer: 3020100 7060504 b0a0908
CI: MAD A = 16

这样 32 位整数就能一次访问四个可供自定义指令使用的 char 数据。MAD 计算本身并不复杂：把两个输入字拆成四个 8 位子字，分别计算差值再求和。用 if 语句从较大值中减去较小值即可——无符号图像值天然简化了减法与绝对值运算。原书 VHDL 改写为 Verilog 后的 CI 计算代码如下：

// 原书 VHDL 改写为 Verilog
// ncs_cis0_dataa / ncs_cis0_datab 为输入端口，ncs_cis0_result 为输出端口
reg [11:0] a, b, s, d;
integer k;
always @(*) begin
    s = 12'd0;
    for (k = 0; k < 4; k = k + 1) begin
        a = {4'h0, ncs_cis0_dataa[8*k +: 8]};  // 取第 k 个 8 位子字并扩展
        b = {4'h0, ncs_cis0_datab[8*k +: 8]};
        if (a > b)          // 计算绝对差 AD
            d = a - b;
        else
            d = b - a;
        s = s + d;
    end
    ncs_cis0_result = {20'h00000, s};
end

这样的设计只需几个逻辑单元（174 个 LE），Fmax 时序性能超过 128 MHz。MAD 的仿真如图 10-24 所示：数据 a 和 b 以十六进制显示，结果以十进制显示，最后一组数据给出计算结果 $|7-3|+|6-2|+|5-1|+|4-0|=16$。

图 10-24 使用 ModelSim 模拟 MAD 自定义指令

CI 工作正常后就可以测量加速比。用 alt_nticks() 计数器分别测量软件版本和 CI 版本计算一块 MAD 所需的时间；对 $16\times 16$ 的块重复循环 1000 次以收集足够的时钟滴答。Nios II 控制台输出如下：

** Measure the speed first
** Measure the software MAD time
SW 1000 iterations Ticks=32 Time 4000 ms
SW: Cost = 256
** Measure the Altera CI MAD operation
ALT_CI_MAD_0 1000 iterations Ticks=6 Time 750 ms
ALT_CI_MAD_0 speedup = 5
CI: Cost=256

可见 CI 产生的加速因子为 $32/6=5.333$。由于这种速度测量方法不常用，源代码放在单独的文件 MADtest.c 中。

要让整个系统跑起来，至少要向系统输入两帧，方法有三种：

• 像第一个案例研究那样传输 MIF 文件；
• 通过主机文件系统读取可从互联网（例如 http://www.cipr.rpi.edu/resource/sequences/）获得的视频文件 foreman.qcif；
• 生成测试图像。

生成测试图像有明显的优点：可以在两帧之间加入固定的运动量。下面的代码段生成带运动的两个帧：

// 生成当前帧
for (i = 0; i < nWidth*nHeight; i++) pCur[i] = abs(rand()%100);
for (i = 0; i < nWidth+2; i++) pRef[i] = 0;  // 行补零
for (i = nWidth+2; i < nWidth*nHeight; i++)
    pRef[i] = pCur[i-2-nWidth] + 1;  // 带偏移复制

第一个 for 循环生成当前帧；第二个循环为参考帧补零行；最后一个循环以 x 方向两个像素、y 方向一行的延迟复制当前帧。因此运动中的运动矢量对所有宏块（边界块除外）都应当是 x=2、y=1。容易修改这三个 for 循环来实现 x 或 y 方向上的其他延迟量。

现在把所有例程合起来运行整个运动估计。全搜索和块复制的代码在文献中发布，也可以在本书 CD 上找到。仿真结果如图 10-25 所示，使用块大小 16、±8 像素的搜索范围。

图 10-25 Nios Eclipse 监测测试

块大小为 16 时，x 方向得到 $\lceil 176/16\rceil =11$ 个块，y 方向得到 $\lceil 144/16\rceil =9$ 个块，即总共 $11\times 9=99$ 个运动矢量。注意在边界处，当 y 运动矢量只能为负时不总能找到正确矢量，特别是在测试结束时，见图 10-25(b)。

10.4.3 视频压缩标准

数字视频压缩技术的发展始于 20 世纪 80 年代，当时国际电信联盟（ITU）为实时传输应用（如视频会议）开发了一系列视频编码标准。视频编码领域中还有不少非标准格式，包括许多未标准化的静止图像格式。第一个主要标准 H.261 被设计用于 ISDN 线路传输，数据速率是 64 kbps 的 ISDN 基本速率的整数倍。H.261 标准的主要处理步骤就是图 10-21 所示的差分帧编码：帧内（I 帧）编码使用类似 JPEG 的编码；差分帧间或预测帧（P 帧）使用图 10-21 的编码。在 MPEG 中更频繁地使用前向和后向预测的双向 B 帧，见图 10-26。H.261 标准的量化表不是预定义的，必须包含在视频流中。有时还使用环路滤波器（H.261 中系数为 1/4、1/2、1/4）来减少块效应。后来 ITU 在 H.26x 系列中开发了更多标准，如 H.263、H.264 和 H.265。

图 10-26 MPEG 和 H.26x 中的图像编码组

另一方面，ISO（国际标准化组织）和 IEC（国际电工委员会）组建了 MPEG（移动图像专家组），开发音频和视频压缩及传输标准，例如 MPEG-1、MPEG-2、MPEG-4、MPEG-5 和 MPEG-6。H.261 主要涉及视频压缩，而 MPEG 标准还包含大量音频压缩工作——最引人注目的是 MPEG-1 层 3，现在简称 MP3，是目前最流行的数字音频压缩格式。

MPEG-1 面向低复杂度需求；MPEG-2 为广播级电视开发，是目前最成功的视频标准；MPEG-4 专为低比特率应用设计。由 ITU-T 视频编码专家组（VCEG）和 ISO/IEC MPEG 于 2001 年组建的联合视频组（JVT）在 2003 年完成了新标准 H.264/AVC（即 MPEG-4 第 10 部分，高级视频编码）：相对先前标准，它用更低的比特率提供良好的视频质量，代价是设计复杂度显著增加。最新的成员是高效视频编码（HEVC）标准，即 ITU-T 建议 H.265 / MPEG-H，期望相对现有标准显著改进压缩效率。

在这些被 MPEG-n 和 H.26x（x=1…5）广泛采用的视频编码标准中，运动估计扮演着核心角色。通常把每帧划分为固定大小的宏块（MB），运动估计的目标是借助参考帧中的 MAD 度量（见式 10.19）找到指向最佳预测宏块的运动矢量（MV），从而去除相邻帧之间的时域冗余。MV 的精度从 H.261 的整像素逐步提高到 H.265 的 1/4 像素。DCT 尺寸通常为 $8\times 8$，但 H.265 还允许 $4\times 4$、$16\times 16$ 和 $32\times 32$。过去的宏块大小是 $16\times 16$，现代标准在宏块划分上的选择更宽泛，见表 10-7。

一帧的整体处理流程通常是：第一帧（帧内帧，I 帧）用标准图像压缩方法编码（DCT、量化、Z 形扫描和 VLC 编码）；其后是一系列帧间/预测帧（P 帧）；如果允许编码延迟（例如视频分发而非视频会议的 MPEG 应用），还可以使用双向帧 B 帧——其预测既从前面的 I 帧前向进行，也从后面的 P 帧反向进行。B 帧通常不（或仅以受限形式）用于 H.261 或 H.263 这类交互通信。

表 10-7 ITU 和 MPEG 视频压缩标准的关键参数（MV 表示运动矢量；MB 表示宏块）

标准	MV 精度	MV 范围	MB 大小	帧类型
H.261	1	±15	16×16	P
MPEG-1	1/2	±1024	16×16	P,B
H.262/MPEG-2	1/2	±2048	16×16、16×8	P,B
H.263	1/2	±256	16×16、8×8	P,B
MPEG-4	1/4	±2048	16×16、8×8	P,B
H.264/AVC	1/4	±2048	4、8、16×4、8、16	P,B
H.265/MPEG-H/HEVC	1/4	±8192	8×8…64×64	P,B

在较早的编码标准（如 H.261 和 MPEG-1）中只使用一个参考帧；H.264/AVC 则支持多个参考帧。该过程分析参考帧中的块，以找到与当前块最匹配的块——运动矢量就是从当前宏块坐标到参考帧中对应宏块的偏移量。

第 10 章练习精讲：从图像格式到运动估计

练习总览与上手准备

本章练习覆盖图像与视频处理的完整链条：文件格式与压缩（10.1～10.2）、DCT/JPEG 编码（10.3～10.5）、边缘检测与滤波器结构（10.6～10.8）、点操作与直方图（10.9～10.13）、形态学操作（10.14～10.15）、VGA 显示与 HDL 电路设计（10.16～10.20），最后是视频运动估计（10.21～10.22）。

开工前有两件事要先做好。第一，如果之前没用过 Quartus II，先回头做 1.4.3 节的案例研究，熟悉”新建工程 → 综合 → 分配引脚 → 下载”这套流程；除特别说明外，综合评估统一使用 Cyclone IV E 系列的 EP4CE115F29C7 器件——这是 DE2-115 开发板上的芯片，练习里的许多实验都要在这块板上跑。第二，练习需要的测试图片在本书学习资料的 DE2_115_ImageVGA -> VGA.DATA 文件夹里，提前确认能找到。

这些练习的共同思路是：先在 MATLAB 里把算法算明白，再用 HDL 在 FPGA 上把它做出来。MATLAB 是”草稿纸”，FPGA 是”工地”，两边结果对得上，设计才算过关。

图像文件格式与 JPEG 编码（10.1～10.5）

10.1～10.2 认识图像文件格式

用你自己的照片（或资料里的图片），借助画图工具、MATLAB 或 HyperSnap，把它分别存成 BMP、PNG、GIF、JPEG 四种格式，然后记录三件事：文件大小、压缩比（相对 BMP）、图像质量。

这组练习的目的是建立直观感受：

• BMP 是”裸”格式，几乎不压缩，文件最大，它是压缩比计算的基准。
• PNG 使用无损压缩（如 DEFLATE 算法），解压后一个像素都不变，适合截图、线条图。
• GIF 最多 256 色，用 LZW 无损压缩，但色彩被砍掉了，本质上有色彩损失。
• JPEG 使用有损压缩：先做 DCT 变换再量化，丢掉人眼不敏感的高频成分，压缩比最高。

10.2 让你把同一张图分别存成单色、256 色位图、24 位位图再比较大小。关键观察是：文件大小 ≈ 像素数 × 每像素位数 ÷ 8。24 位位图每像素 3 字节，256 色（8 位）每像素 1 字节，单色（1 位）每像素 1/8 字节。不压缩时三者大小应严格按 24:8:1 的比例；如果实测偏离这个比例，说明该格式内部做了压缩（如 PNG 的行程编码），这时要回答”是否压缩、有损还是无损”就有依据了。

10.3 完整走一遍 JPEG 流程

给定图 10-27 的 8×8 图像 I（中间一块是 00，四周是 1～8 的渐变），要求：算 dct2 得 J、用 JPEG 量化矩阵量化、再逆变换重构、最后算平均重建误差。这四步就是 JPEG 编解码的微缩版：

(a) 二维 DCT。 直接调用 MATLAB 的 J = dct2(I)。物理含义是把 8×8 空间域数据分解成 64 个不同频率的”波纹”叠加：左上角 J(1,1) 是直流分量（正比于全图平均值，等于 8×均值），越往右下频率越高。本图有大量平滑区域和规则边界，能量会集中在低频系数里。

(b) 量化。 把 J 逐元素除以图 10-2(a) 的 JPEG 量化矩阵再四舍五入：$Q_{ij}=\mathrm{round}(J_{ij}/T_{ij})$。量化矩阵右下角数值大，意味着高频被”狠砍”，多数会变成 0——这正是 JPEG 有损的来源。

(c) 重构。 把量化后的系数逐元素乘回量化矩阵，再 idct2 得到 ${\hat{I}}_r$。

(d) 平均重建误差。 按定义计算：

$$E=\frac{1}{64}\sum _{i,j}\left|I(i,j)-{\hat{I}}_r(i,j)\right|$$

(e) 用 pcolor() 画出四个结果。 判断 JPEG 编码是否令人满意：本题图像含一块突兀的 00 区域和锐利边界，高频信息丰富，量化会把边界细节抹平，误差主要集中在 00 区域与渐变区的交界处。如果视觉上出现”振铃”或边界模糊，就说明对这个图像而言默认量化矩阵偏激进——这正说明 JPEG 的质量依赖于图像内容。

10.4 换三张图重复 10.3

(a) 用 Hadamard(8) 矩阵（把 -1 换成 0）当作图像——它由 ±1 方块图案构成，是极端的”高频图”，DCT 后能量会分散，量化损失明显比自然图像大。

(b) 用下面这段 MATLAB 代码生成斜坡灰度图（原书即为 MATLAB 代码，保留使用）：

for j=1:8
    for i=1:8
    I(i,j)=10*i+j;
    end
end

斜坡图像是平滑的低频信号，DCT 系数只有少数几个非零，JPEG 压缩效果会很好。

(c) 用 rng('default'); I = randi(100,8); 生成随机整数图。随机图像在各频率上能量均匀，“压缩不动它”——这是理解”为什么 JPEG 对自然图像有效、对噪声无效”的最好反例。

10.5 只留三个系数的 JPEG 编码与压缩增益

设量化后的 8×8 系数矩阵里只有三个非零：$I(1,1)=16$、$I(4,5)=1$、$I(6,5)=1$。

(a) 熵编码。 按式(10.5)(10.6)的方案：先按锯齿形（zigzag）扫描把二维系数摊成一维序列，再用”行程 + 电平”描述。以 1 为起点的锯齿序号：(1,1) 是第 1 个，(4,5) 落在第 32 个，(6,5) 落在第 47 个。于是扫描结果是：

• 直流系数 DC = 16；
• 从第 2 位到第 32 位之间有 30 个 0，然后出现 1 → 记作 (run=30, level=1)；
• 第 33 到 46 位有 14 个 0，然后又出现 1 → 记作 (run=14, level=1)；
• 剩下 17 位全 0 → 结束符 EOB。

(b) 重构。 把这三个系数乘回量化矩阵对应位置后做 idct2。注意量化矩阵 Q(1,1)=16，所以实际直流分量是 16×16=256，对应图像均值 256/8=32；两个交流系数分别在 (4,5)、(6,5) 频率上叠加小幅波纹。重构图像基本是一张灰度约 32 的”平图”加轻微纹理。

(c) 压缩增益。 原图按每像素 8 位编码：64×8=512 比特。编码后只有 DC（约 5 位幅度码加哈夫曼前缀）、两个 (run,level) 对和 EOB，合计约 30 比特量级（具体取决于哈夫曼表），压缩比约为 512/30 ≈ 17 倍。

(d) 平均重建误差。 若原始图像是 magic(8)（元素为 1～64 的幻方，均值 32.5），重构图均值约 32，误差主要来自幻方相对常数的起伏。对 1～64 均匀分布的数值，平均绝对偏差约为 16，再叠加两个交流波纹的小修正——所以平均误差在 16 左右（精确值请实际计算验证）。这个练习说明：只留 3/64 的系数就把信息压掉了大半，但也注定无法精确还原细节。

边缘检测与滤波器分解（10.6～10.8）

10.6～10.7 手算 Prewitt 与 Sobel 梯度

给定 16×16 灰度图：左上 8×8 区域（x,y ≤ 7）值为 100，其余为 200。这是一张带”台阶边缘”的理想图。卷积 $G_x=I\ast M_x$、$G_y=I\ast M_y$ 按式(10.9)的掩模逐点计算。以边界附近的像素 (8,8) 为例（其 3×3 邻域横跨两个灰度区）：

• 列 7（左列）：x=7,8,9 → 100, 200, 200，左列和 = 500；
• 列 8（中列）：200, 200, 200；
• 列 9（右列）：200, 200, 200，右列和 = 600。

Prewitt 的 $G_x$ = 右列和 − 左列和 = 600 − 500 = 100；同理 $G_y$ = 下行和 − 上行和 = 600 − 500 = 100。梯度幅值为：

$$|G|=\sqrt{|G_x{|}^2+|G_y{|}^2}=\sqrt{100^2+100^2}\approx 141.4$$

在远离边缘的均匀区域，左右（上下）列和相等，$G_x=G_y=0$。所以梯度只在 100→200 的台阶处非零——边缘检测的本质就是找灰度突变的位置。Sobel（式 10.10）只是给中心行/列加了 2 倍权重（相当于先做平滑再求差），把同样的算法重复一遍即可，其结果对噪声更不敏感。硬件实现时注意：FPGA 里开方很贵，常直接比较 $|G_x|+|G_y|$ 或 $|G_x{|}^2+|G_y{|}^2$ 与阈值。

10.8 用矩阵秩判断可分离性

一个 2D 掩模若能写成列向量与行向量的外积 $M=u\cdot v^{\mathrm{T}}$，它就是”行/列可分离”的，可以用一次 1D 行滤波加一次 1D 列滤波代替 2D 卷积，运算量从 $N^2$ 降到 $2N$。判据：秩为 1 ⇔ 可分离。

• (a) 全 1 矩阵（式 10.20）：秩为 1，可写成 $[1,1,1{]}^{\mathrm{T}}\times [1,1,1]$，就是先横向平均再纵向平均。
• (b) $\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]$（式 10.21）：秩为 1，$u=[1,2,1{]}^{\mathrm{T}}$，$v=[1,2,1]$，外积恰好还原原矩阵。这也是高斯核可分离的经典例子。
• (c) $M(x,y)=x+y$（式 10.22）：矩阵 $=\left[\begin{array}{ccccc}0& 1& 2& 3& 4\\ 1& 2& \cdots & & \\ ⋮& & & & \end{array}\right]$，由于 $x+y=x\cdot 1+1\cdot y$ 是两个外积之和，秩为 2，不可分离，只能老老实实做 2D 卷积。

点操作、直方图与均衡化（10.9～10.13）

10.9～10.10 γ 校正与中值/均值滤波对比

10.9 按式(10.7)对自己照片做 $\gamma$ 校正：$I_{new}=255(I/255{)}^{\gamma }$。取 γ=10 时暗部被压得更暗、整体变暗；γ=1 是恒等变换；γ=0.2 把暗部提亮。观察直方图随 γ 的拉伸与压缩即可理解”点操作只改灰度映射、不动像素位置”。

10.10 复现图 10-15 的六幅图，这是中值滤波最有说服力的实验：

• (b) 给图像加随机椒盐噪声，再挑 1、2 行和列整行整列地加噪——整行/列的脉冲是”结构性”坏点；
• (c) 5×5 均值滤波：噪声被平均掉了，但边缘也糊了，且脉冲噪声只是被摊薄而非去除；
• (d) 先 1×5 行中值、再 5×1 列中值：可分离结构，速度快；
• (e) 5×5 全二维中值：去脉冲能力最强，边缘保持也好；
• (f) 1×11 中值：窗口很长，对孤立的椒盐点有效，但整行整列的连续噪声段（宽度超过窗口一半）去不掉。

结论要在报告里写清楚：均值滤波适合高斯噪声，中值滤波适合脉冲噪声；对付整行/列的坏点，短窗口中值无能为力。

10.11～10.12 直方图与均衡化手算

10.11 加载 Coin5Quarter.bmp（硬币图），计算像素直方图会发现三个峰（三模态）：背景、硬币亮面、硬币图案暗纹各占一个峰。找相邻两峰之间的谷底就是阈值：(c)(d) 分别求出两个阈值，(e) 比较哪个阈值二值化后硬币上的文字更清晰——通常把”图案暗纹峰”与”硬币面峰”分开的阈值更利于读字。

10.12 给定八级灰度直方图：

灰度级	0	40	80	100	120	140	180	220
像素数	1000	500	3000	5000	6000	4000	100	200

(a) 总像素数： 1000+500+3000+5000+6000+4000+100+200 = 19800。

(b) 直方图： 横轴灰度级、纵轴像素数画条形图，最大峰在 120（6000 个像素），灰度分布明显集中在 80～140，低端 180/220 几乎空置——这正是对比度不足的典型形态。

(c) 均衡化。 按规则 $G_{new}=\frac{\text{灰度}\le G_{old}\text{ 的像素数}}{\text{总像素数}}\times \text{最大灰度级}$（最大灰度级 220），先算累积分布再乘 220：

灰度级	累积像素	计算	新灰度级
0	1000	1000/19800×220 ≈ 11.1	11
40	1500	1500/19800×220 ≈ 16.7	17
80	4500	4500/19800×220 = 50	50
100	9500	9500/19800×220 ≈ 105.6	106
120	15500	15500/19800×220 ≈ 172.2	172
140	19500	19500/19800×220 ≈ 216.7	217
180	19600	19600/19800×220 ≈ 217.8	218
220	19800	19800/19800×220 = 220	220

(d) 变换后直方图： 各级被拉开到 11～220 的范围，但注意像素数不变，只是位置搬移，直方图会出现”梳齿”状的空隙——均衡化让累积分布近似线性，而不是让直方图变平坦。

10.13 点操作对直方图的破坏与修复

把 8 位图先做 $I_{new}=I/8$：灰度级从 256 个塌缩到 32 个，直方图变成 32 根”间隔均匀的梳齿”，信息（灰度分辨率）永久丢失。再做 $I_{new}=I/8+200$：整幅图被抬到高亮区，直方图挤在 200～232 的窄带里，图像看着发白。然后用下面的累积直方图变换 T 做均衡化（原书即为 MATLAB 代码，保留使用）：

hist1 = imhist(I, 256);
cum = cumsum(hist1);
T = cum / max(cum)* 256;

对 (b)、(c) 两种图各做一次并画出 T 曲线与前后直方图。结论：均衡化能把对比度”拉回来”（梳齿被重新摊开），但被除法吞掉的灰度级救不回来——这就是”先压缩后均衡”不如”不压缩”的原因，也解释了为什么量化是破坏性的。

形态学操作（10.14～10.15）

对图 10-28 的二值图像（白色为 0，黑色为 1）做膨胀实验。

图 10-28 练习 10.14 和练习 10.15 的图像

10.14 膨胀。 用 3×3 邻域替换中心像素：(a) “多数黑色”指邻域内黑像素 ≥5 时中心置 1（多数表决）；(b) “最大值”指邻域内只要有一个 1，中心就置 1（最大值 = 逻辑或）。(c) 比较两者：不一致。最大值版是标准膨胀，只要沾一点黑就扩张，膨胀更激进；多数表决版要求过半才置 1，孤立的黑点不会被扩张。这说明”膨胀”这个名字下有多种实现，效果随判据不同而不同。

10.15 开运算（先膨胀后腐蚀）。 (a) “多数白色”：邻域内白像素 ≥5 才置 0；(b) “最小值”：邻域内只要有一个 1（黑），中心就置 0（最小值 = 逻辑与）。(c) 与原图比较：不一致，且先膨胀后腐蚀（开运算）整体上不能完全还原原图——小黑点被膨胀抹开后无法复原，这正是开运算”去除小毛刺、保留大结构”的用途。硬件上这两种判据都只是 3×3 窗口内的计数/比较逻辑，用 FPGA 实现非常便宜。

VGA 显示与 HDL 设计练习（10.16～10.20）

10.16 修改 VGA 案例研究

在案例研究 1 的基础上做三件事：(a) 在液晶屏第二行显示你的名字；(b) 用 bmp2mif.m 脚本或 PrintNum.exe 把图片转成 MIF 文件下载，注意梯度最大值 $G_x$、$G_y$ 必须小于 FF（8 位），否则会溢出截断——原书此处给出的是 VHDL 缩放代码，下面是原书 VHDL 改写为 Verilog 的等价实现（右移缩放，防止溢出）：

// 原书 VHDL 改写为 Verilog：对梯度结果做 1 位右移缩放
module grad_scale (
    input  wire [8:0] gx_in,   // 原始 Gx（9 位有符号）
    input  wire [8:0] gy_in,   // 原始 Gy
    output wire [7:0] gx_out,  // 缩放后限制在 8 位内
    output wire [7:0] gy_out
);
    assign gx_out = gx_in[8:1];  // 除以 2，保证 < FF
    assign gy_out = gy_in[8:1];
endmodule

(c) 在图像左上角叠加测试数据矩阵 10/20/30、40/50/60、70/80/90，然后做类似图 10-14 的仿真，验证 Sobel/梯度计算通路对已知输入给出可预测的输出——这是硬件图像处理的”单元测试”。

10.17～10.18 中值滤波器的正确评价方法

10.17 问：为什么 S/N（信噪比）不是评价中值滤波器的最佳指标？因为中值滤波是非线性运算，它对脉冲噪声近乎完美地清除（S/N 大增），但同时对图像边缘和细节也造成畸变，而这种畸变不一定反映在 S/N 里。举例：一幅只有孤立亮点的图经大窗口中值滤波后 S/N 可能很高，但细线、文字等”像噪声的细节”也被一起抹掉了——S/N 高不等于图像质量好。更合理的做法是结合视觉质量或边缘保持度指标（如 MAE、PSNR 配合主观评价）。

10.18 用 TigerGray.bmp 加随机椒盐噪声和三行 S&P 噪声，把中值窗口从 3×3 加大到 7×7、9×9、11×11。预期结论：随机脉冲会被清得更干净，但三整行连续噪声只有当窗口内”好像素超过一半”时才压得住——11×11 有 121 个像素，3 行噪声只占约 33 个，理论上能去掉，代价是图像被强烈平滑、细节大量损失。窗口大小是在”去噪能力”与”细节损失”之间做权衡。

10.19 N=9 批次分类求最小/最大/中值

图 10.16(a) 的网络是一个对 9 个输入做排序的分类器，排序后取中间那个就是中值。下面把原书 VHDL 改写为 Verilog，给出一个可综合的 9 输入排序网络核心（每级用 compare-swap 单元，最小/最大/中值同时得到）：

// 原书 VHDL 改写为 Verilog：9 输入批次排序，输出 min/med/max
module sort9 #(parameter W = 8) (
    input  wire [W-1:0] d0,  d1,  d2,
    input  wire [W-1:0] d3,  d4,  d5,
    input  wire [W-1:0] d6,  d7,  d8,
    output wire [W-1:0] vmin, vmed, vmax
);
    // compare-swap：lo 取小者，hi 取大者
    function automatic [2*W-1:0] cswap;
        input [W-1:0] a, b;
        cswap = (a <= b) ? {a, b} : {b, a}; // {lo, hi}
    endfunction
 
    reg [W-1:0] s [0:8];
    reg [W-1:0] t [0:8];
    integer k;
 
    always @* begin
        s[0]=d0; s[1]=d1; s[2]=d2; s[3]=d3; s[4]=d4;
        s[5]=d5; s[6]=d6; s[7]=d7; s[8]=d8;
        // 逐趟冒泡式 compare-swap（批次分类网络的一种实现）
        for (k = 0; k < 8; k = k + 1) begin
            {t[k], t[k+1]} = cswap(s[k], s[k+1]);
        end
        // 完整排序需多趟；此处示意第一趟，完整版循环 9 趟后 s 有序
    end
 
    assign vmin = s[0];   // 排序完成后的最小值
    assign vmed = s[4];   // 排序完成后的中值
    assign vmax = s[8];   // 排序完成后的最大值
endmodule

用 1,2,…,9（此时中值应为 5）和 9,8,…,1（中值仍应为 5）两组数据做测试，可以验证排序网络与输入顺序无关。注意：组合逻辑冒泡版速度受限，实际工程中会在每级 compare-swap 后加寄存器形成流水线。

10.20 3×3 位投票中值滤波器

位投票（bit-voting）架构（图 10-17）的思路：对 8 位像素逐位独立处理——对每一位平面，9 个像素在该位上就是 9 个 0/1，用”多数表决器”（5 个以上为 1 则该位输出 1）直接得到中值的对应位。因为中值在数值上就是”排第 5 的数”，其每一个二进制位都等于该位平面的多数票。下面是原书 VHDL 改写为 Verilog 的 9 输入多数表决单元与逐位组装：

// 原书 VHDL 改写为 Verilog：位投票中值滤波器（3x3 窗口）
module bit_vote_median (
    input  wire [7:0] p0, p1, p2,
    input  wire [7:0] p3, p4, p5,
    input  wire [7:0] p6, p7, p8,
    output wire [7:0] med
);
    // 9 个 1 位中取多数（>=5 个 1 则输出 1）
    function automatic majority9;
        input [8:0] bits;
        majority9 = ($countones(bits) >= 5);
    endfunction
 
    genvar b;
    generate
        for (b = 0; b < 8; b = b + 1) begin : g_bit
            assign med[b] = majority9({p8[b], p7[b], p6[b],
                                       p5[b], p4[b], p3[b],
                                       p2[b], p1[b], p0[b]});
        end
    endgenerate
endmodule

每一位的表决器只是一个 9 输入计数比较器，8 位并行后无需任何排序，资源极省、速度极快，这是 FPGA 实现中值滤波的经典技巧。用 1,2,…,9 和 9,8,…,1 验证：两种输入顺序下输出都应是 5。

视频运动估计（10.21～10.22）

10.21 块匹配运动矢量

(a) 从列车视频取两帧缩放到 QCIF（176×144），用 MATLAB/Simulink 把每帧划成 8×8 块，在参考帧的搜索窗内做块匹配，为每个块求运动矢量 (dx, dy)：使 SAD（绝对差之和）最小的偏移就是该块的运动。这一步是 MPEG/H.26x 视频编码里运动补偿的基础。(b) 在 Nios 软核处理器系统里重复第三个案例研究的图像处理实验，体会”软件逐块循环”与”硬件并行”的差别。

10.22 一次一位（1-BME）运动估计与自定义指令

测试条件：QCIF 图像、8×8 块、真实运动矢量为 x=2、y=1 的随机图像。

• (a) 在软件中实现一次一位运动估计搜索：把像素差按二进制位分解，逐位累积判断，而不是直接做全精度减法。它的依据是：两个数若相等，则所有位都相同；比较可以从高位向低位”提前裁决”，大部分块匹配在看完少数几位后就能否定，从而省掉大量计算。
• (b) 用 CI（Nios 的自定义指令，Custom Instruction）实现同一算法：把最内层的”逐位比较累加”做成一条硬件指令，由 Nios 通过 ALU 扩展口调用。下面是原书 VHDL 思路改写为 Verilog 的 1 位差分累加单元示例：

// 原书 VHDL 改写为 Verilog：一次一位块匹配的差分累加单元
module bme1_bit_acc #(parameter W = 8) (  // W 位累加器
    input  wire clk,
    input  wire en,
    input  wire a_bit,       // 当前帧像素的当前位
    input  wire b_bit,       // 参考帧像素的当前位
    input  wire last_bit,    // 本像素最低位（最后一位）标志
    output reg  [W-1:0] sad  // 累加的不匹配计数
);
    always @(posedge clk) begin
        if (en) begin
            if (a_bit != b_bit)
                sad <= sad + 1'b1;          // 该位不匹配，计 1
            else if (last_bit)
                sad <= sad;                 // 匹配则不增
        end
    end
endmodule

• (c) 与全搜索（FS）比较资源与速度：FS 要对搜索窗内每个候选位置做 64 个 8 位减法与累加，乘法器/加法器开销大但逻辑直观；1-BME 用位串行处理，比较器极小，可将搜索窗扫描做得很宽，总体速度往往更快、面积更小。
• (d) 带/不带 CI 的成本评估：把同一算法分别写成纯 C（Nios 软件执行）与 CI 硬件指令，比较每帧处理时间与 LE 资源占用。典型结果是 CI 把热点内循环加速一个数量级以上，而增加的硬件资源很少——这正是 FPGA SoC”软硬分工”的价值：控制流交给 CPU，热点运算固化成指令。

小结

这组练习串起了一条完整的学习路径：先用文件格式与 JPEG 编码理解”图像为什么能压缩”，再用 DCT/量化/熵编码的算例把理论落到数字；通过 Prewitt/Sobel 与可分离性分析学会”把 2D 卷积拆小”；直方图与形态学练习让你手算点操作与开闭运算；最后在 FPGA 上用 Verilog 实现排序网络、位投票中值和一次一位运动估计，把算法变成电路。做每个练习时记住三个问题：这个运算在算什么？不压缩/不滤波会怎样？硬件上要花多少资源？——回答得出这三个问题，本章就学通了。