第 5 章图像复原

第5章图像复原

图像复原－又称为图像恢复 • 与图像增强相似－都要得到在某种意义上改进的图像，或者说，希望要改进输入图像的视觉质量 • 不同之处－图像增强技术一般要借助人的视觉系统的特性，以取得看起来好的视觉结果，而图像复原则认为图像是在某种情况下退化或恶化了(图像品质下降了)，现在需要根据相应的退化模型和知识重建或恢复原始的图像 • 图像恢复技术是要将图像退化的过程模型化，并据此采取相反的过程以得到原始的图像

图像恢复的内容 • 退化模型和循环矩阵对角化 • 复原的代数方法 • 逆滤波 • 最小二乘方滤波 • 交互式恢复 • 空间复原技术

退化模型和循环矩阵对角化 传感器非线性畸变 • 退化模型光学系统中的衍射光学系统的像差产生原因几何畸变摄影胶片的非线性图像运动造成的模糊大气流的扰动效应

定义： • f[x,y]: 原始图像 • g[x,y]: 退化图像 • n[x,y]: 加性噪声 • g[x,y] = H{f[x,y]} + n[x,y] • H{ }: 系统或操作图像恢复就是在给定g(x,y)和代表退化的H的基础上，得到对f(x,y)的某个近似的过程

简单的通用退化模型 n(x,y) g(x,y) f(x,y) H{.} +

常数 2幅图像 H多具有的性质 • 线性： • H{k1f1 + k2f2} = k1H{f1} + k2H{f2} • 相加性：令 k1 = k2 = 1，则 • H{f1 + f2} = H{f1} + H{f2} • 一致性：令 f2 = 0，则 • H{k1f1} = k1H{f1} • 位置（空间）不变性： • H{f[x-a, y-b] } = g[x-a, y-b] 图像任意位置的响应只与在该位置的输入值有关，而与位置本身无关

常见具体退化模型示例 线性空间不变非线性摄影胶片的冲洗过程模糊退化光学成像系统，由于孔径衍射产生的退化目标运动造成的模糊退化随机噪声迭加，随机性的退化

退化模型的计算 假设对2个函数f(x)和h(x)进行均匀采样，其结果放到尺寸为A和B地2个数组。对f(x)，x的取值范围是0,1,2…A-1；对h(x)，x的取值范围是0,1,2,….B-1。利用卷积计算g(x)。为了避免卷积的各个周期重叠，取M≥A+B-1，并将函数用0扩展补齐 fe(x)和he(x)表示扩展函数，卷积为 ge(x)=∑fe(m)he(x-m)x=0,1,…M-1 矩阵表示 g和f是M维列矢量： fT = [ f[0], f[1], …, f[M-1] ] gT = [ g[0], g[1], …, g[M-1] ] g=Hf

H称为M×M循环矩阵 H= 考虑噪声 g=Hf+n (1)

循环矩阵对角化 如果直接对式(1)进行计算求解f，计算量达，如M=N=512 ,则H的尺寸为262144×262144，可以通过对角化H来简化当k=0,1…M-1时，循环矩阵H(设为M×M)的特征矢量和特征值分别为

将H的M个特征矢量组成1个M×M的矩阵W W = [w(0) w(1) w(2) … w(M-2) w(M-1)] H = WDW-1D = W-1HW 其中：D(k,k) = λ(k)，WW-1 =W-1W=I

复原的代数方法 • 图像复原的主要目的是当给定退化的图像g以及H和n的某种假设，估计出原始图像f • 代数复原方法的中心是寻找一个估计的f^，它使事先确定的某种优度准则为最小

无约束复原方法 • 由退化模型可知，其噪声项为： n= g-Hf 在并不知道n的情况下，希望找到一个f^，使得Hf^在最小二乘方意义上来说近似于g，也就是说，寻找一个f^，使得 ||n||2 = ||g – Hf^||2  J(f^)

或 nTn = (g – Hf^)T(g – Hf^) 实际上是求J(f^)的极小值问题，除了要求J(f^)为最小外，不受任何其它条件约束，因此称为无约束复原 dJ(f^ )/df^ = 0 = -2HT(g – Hf^) 即 f^ = (HTH)-1 HTg (2) M=N时，则有 f^ = H-1(HT)-1 HTg = H-1 g

约束复原方法 • 在最小二乘方复原处理中，为了在数学上更容易处理，常常附加某种约束条件。如令Q为f的线性算子，最小二乘复原问题可看成是使形式为||Qf^||2函数，服从约束条件||g – Hf^||2 = ||n||2 的最小问题，这种带有附件条件的极值问题可用拉各朗日乘数法处理

处理过程 寻找一个f^，使下述准则函数为最小 J(f^) = ||Qf^||2 + α{||g – Hf^||2 - ||n||2} 拉各朗日系数 dJ(f^ )/df^ = 0 f^ = (HTH + γQTQ)–1 HTg α=1/λ

逆滤波 • 基本原理在不考虑噪声的情况下，假设M=N，则根据前面的公式，有反向滤波的作用 f^ = H–1 g = WD–1 W–1 g 或 W–1 f^ = D–1 W–1 g F^(u,v) = G(u,v)/H(u,v) 傅立叶变换经过傅立叶反变换，可求得原始图像f(x,y)

F^(u,v) = F(u,v) + N(u,v)/H(u,v) • 在有噪声的情况下从上面两式可以看出，在进行复原处理时可能会发生下列情况： • H(u,v)=0或H(u,v)非常小，在这种情况下，即使无噪声，也无法精确恢复f(x,y) • 在有噪声存在时，在H(u,v)的邻域内，H(u,v)的值可能比N(u,v)的值小的多，由上式得到的噪声项可能会非常大，不能使f(x,y)正确恢复

一般说，逆滤波不能正确估计H(u,v)的零点 实际中，不用1/H(u,v)，而用另外一个关于u,v的函数M(u,v) 处理框图为： F^(u,v) G(u,v) + f(x,y) M(u,v) H(u,v) N(u,v) 1/H(u,v) u2+v2≤w20 M(u,v)= 1 u2+v2＞ w20

例子散焦模糊利用大的邻域进行恢复原始图像利用原始图像的一个邻域光谱面恢复

消除匀速直线运动引起的模糊 • 在获取图像时，由于景物和摄像机之间的相对运动，造成图像的模糊。假设对平面匀速运动的景物采集1幅图像f(x,y)，并设x0(t)和y0(t)分别是景物在x和y方向的运动分量，T是采集时间长度。实际采集到的由运动而造成的模糊图像g(x,y)为： g(x,y)=∫0Tf[x-x0(t),y-y0(t)]dt G(u,v) = F(u,v) 0Texp{-j2p[ux0(t) + vy0(t)]}dt = F(u,v)H(u,v) H(u,v) = 0Texp{-j2p[ux0(t) + vy0(t)]}dt 如果知道运动分量x0(t)和y0(t)，从上式直接得到H(u,v)

只考虑x方向的运动情况，即x0(t)=ct/T，y0(t)＝0，则只考虑x方向的运动情况，即x0(t)=ct/T，y0(t)＝0，则 H(u,v) = 0Texp{-j2πux0(t)}dt = (T/πuc)sin(πuc)e-jπuc 当n为整数时，H在u=n/c处为0 根据模糊图像g(x,y)，不考虑y的变化 g(t) = 0Tf[x-x0(t)]dt = 0Tf[x-ct/T]dt = x-cxf[t]dt, 0  x  L 其中：t = x – ct/T. 因此， g’(x) = f(x) – f(x-c) or f(x) = g’(x) + f(x-c)

假设 L = Kc, 其中 K使整数. 变量x为： x = z + mc 其中 z 在[0,c]范围， m 是 (x/a)的整数部分 (m 的取值为 0, 1, …, K-1). 则, f(z+mc) = g’(z+mc) + f[z+(m-1)c] 这个等式能通过递归 f(z)解决， f(z) = f(z-a) 定义为运动范围的一部分 0  z < c. f(x) = Sk=0m g’(x-kc) + f(x-mc) 由于g(x)为已知，问题就转化为估计f(x)

为了估计f(x): f(x-mc) = f(x) – f^(x). • 对每个 kc x < (k+1)c计算，并将 k = 0, 1, …, K-1 的结果加起来，得到 • Kf(x) = Sk=0K-1f(x+kc) – Sk=0K-1f^(x+kc), 0  x < c • f(x)  A – (1/K) Sk=0K-1f^(x+kc), 0  x < c • 其中 (1/K) Sk=0K-1f(x+kc) 是未知的，但是当 K很大时，其接近平均值，将其设为常数A。 • 因此, • f(x-mc)  A – (1/K) Sk=0K-1f^(x+kc-mc), 0  x  L •  A – (1/K)Sk=0K-1Sj=0kg’[x-mc+(k-j)c] • 或对 0  x,y  L • f(x,y)  A – (1/K)Sk=0K-1Sj=0kg’[x-mc+(k-j)c,y] + Sj=0mg’[x-jc,y]

最小二乘方滤波 • 最小二乘滤波也就是维纳滤波，它是使原始图像f(x,y)及其恢复图像f^(x,y)之间的均方误差最小的复原方法具体的数学公式推导过程忽略，直接给出公式 • Sf(u,v):为f[x,y]的功率普，Sh(u,v)为 n[x,y]的功率普

讨论一下上式的几种情况 （1）如果s=1，方括号中的项就是维纳滤波器（2）如果s是变量，就称为参数维纳滤波器（3）当没有噪声时，Sn(u,v)=0，维纳滤波器就退化为理想的逆滤波器（4）当Sn(u,v)和Sf(u,v)未知时，用常数K可代替因此必须调节s以满足 f^ = (HTH + sQTQ)–1 HTg

逆滤波和维纳滤波恢复比较 SNR 1 10 100 退化图像傅立叶功率普逆滤波恢复维纳滤波恢复光谱图

模糊和增加噪声 原始图像约束的最小二乘滤波逆滤波恢复

交互式恢复 • 前面讨论都是自动解析的恢复方法，在具体恢复工作中，常常需要人机结合，由人来控制恢复过程，以达到一些特殊的效果实际中，有时图像会被1种2-D的正弦干扰模式(也叫相关噪声)覆盖。令η(x,y)代表幅度为A，频率分量为(u0,v0)的正弦干扰模式，即： η(x,y)=Asin(u0x+v0y) 傅立叶变换为： N(u,v)=-jA[δ(u-u0/(2π),v-v0/v(2π))- δ(u+u0/ (2π),v+v0/ (2π)) ]

这里退化仅由噪声造成，所以有： G(u,v)=F(u,v)+N(u,v) 利用前面讲的带阻滤波器消除，以去掉正弦干扰模式影响 n(x,y) = A sin(u0x+v0y) 傅立叶光谱正弦干扰恢复图像

例子周期干扰退化图像傅立叶谱

干涉模式 处理后图像

空间复原技术 • 空间变换 • 灰度插值

空间变换 • 在图像的获取或显示过程中，产生几何失真，如成像系统有一定的几何非线性，因此会造成如图所示的枕形失真或桶形失真，另外，由于地球表面呈球形，摄取的平面图像也将会有较大的几何失真。对这些图像必须加以校正，以免影响分析精度原始图像桶形失真枕形失真

* * * * * * * * * * * * * * * * • 校正过程 + + 已校正图像观测图像实际空间畸变空间变形校正理想图像观测图像和校正图像之间对应点

设原图为f(x,y)，受到几何形变得影响变成g(x’,y’)，这里(x’,y’)表示失真图像的坐标设原图为f(x,y)，受到几何形变得影响变成g(x’,y’)，这里(x’,y’)表示失真图像的坐标 x’=s(x,y) y’=t(x,y) 线性失真 s(x,y)=k1x+k2y+k3 t(x,y)=k4x+k5y+k6 非线性失真 s(x,y)=k1+k2x+k3y+k4x2+k5xy+k6y2 t(x,y)=k7+k8x+k9y+k10x2+k11xy+k12y2 如果已知s(x,y)和t(x,y)的解析表达，可以通过反变换来恢复图像

在实际中，通常不知道解析表达，需要在恢复过程中的输入图象(失真图)和输出图(校正图)上找一些其位置确切知道的点(称为约束对应点)，然后利用这些点建立2幅图像间其它象素空间位置的对应关系在实际中，通常不知道解析表达，需要在恢复过程中的输入图象(失真图)和输出图(校正图)上找一些其位置确切知道的点(称为约束对应点)，然后利用这些点建立2幅图像间其它象素空间位置的对应关系

灰度插值 (x^,y^) (x,y) f(x,y) g(x^,y^) 最近邻内插

实例变形图像几何校正后的图像

第 5 章 图像复原