第五章 图像恢复和重建版权所有,1997 (c) Dale Carnegie & Associates,Inc.
CHAPTER 5
IMAGE RESTORATION
and RECONSTRUCTION
§ 1 退化的数学模型和对角化
§ 2 无约束恢复
§ 3 有约束恢复
§ 4 几何失真校正
§ 5 投影重建
§ 5.1 退化的数学模型和对角化
退化：图像质量的降低；
失真可看作是退化，校正是恢复；
投影可看作是退化（三维到二维平面），重建是恢复；
§ 5.1.1 简单的通用图像退化模型
n（ x,y）
f（ x,y） H + g（ x,y）
模型化：一个作用在 f（ x,y） 上的系统 H与一个加性噪声 n（ x,y） 的联合作用,导致产生退化图像 g（ x,y） 。
假设已知 n（ x,y） 的统计特性（或先求出），图像复原就是已知 g
（ x,y） 求 f（ x,y） 的问题 （近似过程）。
g（ x,y） = H [f（ x,y） ] + n（ x,y） ；
已知 退化 解 噪声
§ 5.1.2 退化模型的计算一,1D情况：
设 f（ x） 和 h（ x） 均匀采样后放入尺寸为 A,B的数组；
f（ x）,x=0,1,…,A-1； h（ x）,x=0,1,…,B-1；
为了避免卷积的周期重叠，取 M? A+B-1；
将 f（ x） 和 h（ x） 用零扩展补齐； (M周期 )
ge（ x） =? fe（ m） h（ x-m） ； x,m=0,1,…,M-1；
用矩阵形式表示的卷积为 g（ x） = H f,展开得
ge（ 0） = he（ 0） he（ -1） … he（ -M+1） fe（ 0）
ge（ 1） = he（ 1） he（ 0） … he（ -M+2） fe（ 1）
ge（ M-1） = he（ M-1） he（ M-2） … he（ 0） fe（ M-1）
§ 5.1.2 退化模型的计算 （续 1）
式中,he（ 0） he（ -1） … he（ 1）
H = he（ 1） he（ 0） … he（ 2）
he（ M-1） he（ M-2） … he（ 0）
H矩阵是一个循环矩阵，其特点有：行列首尾相接；环形
循环矩阵相加后仍是循环矩阵，相乘后仍是循环矩阵；
§ 5.1.2 退化模型的计算 （续 2）
二,2D情况（ 1D的直接推广,0?x? M-1,0? y? N-1 ）
he(x,y)= h(x,y),0?x? A-1,0? y? B-1,否则 he(x,y)= 0；
fe(x,y)= f (x,y),0?x? C-1,0? y? D-1,否则 fe(x,y)= 0；
ge (x,y) = fe (m,n) he (x-m,y-n) ； m=0,1,…,M-1； n=0,1,…,N-1；
用矩阵形式表示的卷积为 g = H f,展开得
ge（ 0） = H0 HM-1 … H1 fe（ 0）
ge（ 1） = H1 H0 … H2 fe（ 1）
ge（ M-1） = HM-1 HM-2 … H0 fe（ MN-1）
其中每个块 Hi 是由扩展函数 he (x,y)的第 i行而来，是块循环矩阵。
§ 5.1.3 循环矩阵对角化一,1D对角化设 循环矩阵 H的本征值和特征向量分别 为：
(k) = [ 1 exp[j2?/M *k] …… exp[j2?/M* (M-1)k]T；
(k) = he(0)+he(M-1)exp[j2?/M?k]+…… + he(1)exp[j2?/M? (M-1) k] ；
将 H的 M个特征向量组成一个 M*M的矩阵 W：
W=[?(0)?(1)?(2) ……?(M-1)]
H = WDW-1；
式中 D为一个对角矩阵，其元素正是 H本征值，
即 D (k,k)=?(k)。
§ 5.1.3 循环矩阵对角化 (续 1)
设对 4*4的循环矩阵 C对角化。
设?(k)为整数 1的 M次根,k为根的序号，即?(k)M = 1 = lne；
因?(k)=exp[j2?k / M]； 例中 M=4； 即?(k)4 = 1；?(k)=exp[j2?k/4]；
规定标量?(k)= c0+?(k) c1 +?(k)2 c2 +?(k)3 c3= ∑ ci exp[j2?k/4*i]；
根据线性代数分析,k=0,1,2,3时可得到四个本征值和相应的特征向量；
(k)?(k)=C?(k)； 组成?(k)矩阵 W=[?(0)?(1)?(2)?(3)]
各元素?(k) =[ 1 exp (j2?/4?k ) exp (j2?/4?2k ) exp (j2?/4?3k ) ]；
当循环矩阵已知本征值和特征向量后，可以进行对角化；
设对角化矩阵为 D,则 D= W-1CW；（ C等于 H ）；
D为对角阵,D (k,k)=?(k) =MH (k) = W-1HW；
§ 5.1.3 循环矩阵对角化 (续 2)
二,2D对角化（块循环矩阵的对角化）
借助循环矩阵的结果进行推广；
设矩阵 W的尺寸为 MN*MN,每个元素为 W(i,m)=exp(j2?/M*im)WN ；
WN为一个 N*N矩阵，每个元素为 WN(k,n)=exp(j2?/N*kn)；
H = WDW-1 ；且 HT = WD*W-1 ； D*为 D的复共轭；
∴ D = W-1HW；可转化为对角阵。
§ 5.1.4 退化模型对角化的效果一,1D情况：
H=WDW-1 ； 等式代入 g=H f,并两边同左乘 W-1,得 W-1 g= DW-1 f，
W-1 g= DW-1 f 都是 M维向量，令它们分别为 G(k)和 F(k),得 G(k)=DF(k)；
其中,G(k)= (1/M) ∑ ge (i) exp(j2?/M*ki)是 ge (x)的傅里叶变换；
F(k)= (1/M) ∑ fe (i) exp(j2?/M*ki)是 fe (x)的傅里叶变换；
已知 D (k,k)=?(k) = ∑ he (i) exp(j2?/M*ki)=MH (k),
H (k)是 he (i)的傅里叶变换，
所以 G(k)= M× H (k) F(k),式中 H (k) F(k)为 fe (x) 与 ge (x)的卷积，
可采用 FFT计算。
§ 5.1.4 退化模型对角化的效果 (续 1)
二,2D情况：
G(u,v)= (1/MN) ∑∑ ge (x,y) exp[-j2? (ux/M+vy/N)] ；
F(u,v)= (1/MN) ∑∑ fe (x,y) exp[-j2? (ux/M+vy/N)] ；
N(u,v)= (1/MN) ∑∑ ne (x,y) exp[-j2? (ux/M+vy/N)] ；
H(u,v)= (1/MN) ∑∑ he (x,y) exp[-j2? (ux/M+vy/N)] ；
D的 MN个对角元素可用下式表示，
D(k,i)= MN *H ([k/N],k mod N)；当 i=k； 否则为 0；
其中 [k/N]表示取最大整数,mod表示取余数；
所以 G(u,v)= H(u,v) F(u,v) + N(u,v)； u,v=0,…,M -1
§ 5.1.5 有、无约束恢复一、无约束恢复公式
f ‘ = (nTn) -1 HT g = H -1 (HT) –1 HT g = H-1 g ；
二、有约束恢复公式若选取 f ‘ 的一个线性操作符 Q（ 变换矩阵），
使得 ||Q f ‘ ||最小，采用拉格朗日乘数法解决，最小化下式的 f ‘,
L（ f ‘ ） = ||Qf ‘ ||2+ l（ ||g-H f ‘ || 2 -||n|| 2 ），
l 为拉格朗日乘数，该式为准则函数；
最小化得有约束恢复公式：
f ‘ = [HT H + s QT Q ]-1 HT g； s = 1 / l ；
§ 5.2 无约束恢复
§ 5.2.1 逆滤波一、无噪声时处理方法
F ‘ (u,v) = G(u,v) / H(u,v) ； u,v=0,1,…,M-1
上式给定的恢复方法称为逆滤波，条件 H –1 (u,v) 存在。
H(u,v)取零或很小时，计算上较困难。
二、有噪声时处理方法（ G = H F + N,两端同除 H）
F ‘ (u,v) = F(u,v) + N(u,v) / H(u,v) ；
H(u,v)取零或很小时，计算结果差距很大；一般情况下，
逆滤波器并不正好是 1 / H(u,v),记为 M(u,v),称为恢复转移函数；
M(u,v) = 1 / H(u,v) ；如 u2+v2 ≤?02,?0为去掉零点的值；
M(u,v) = 1 ； 如 u2+v2 >?02,
§ 5.2.1 逆滤波 (续 1)
图像退化和恢复模型
n(x,y)
f (x,y) H(u,v) + g (x,y) M(u,v) f ‘ (x,y)
退化 恢复
改进,M (u,v) = k ； 如 H(u,v) ≤ d,d和 k为小于 1的常数；
M(u,v) = 1/ H(u,v) ；其它,d越小越好
H (u,v) 已知，则可恢复原图像；
退化系统的转移函数可以用退化图像的傅里叶变换来近似，即
H (u,v) ≈ G (u,v)
FFT 去零点 IFFT
g (x,y) G(u,v) ≈ H(u,v) M(u,v) f ‘ (x,y)
点冲激下 F(u,v)=1
§ 5.2.2 运动模糊图像的恢复
运动模糊图像的恢复主要是消除匀速直线运动模糊，采用照相机摄影。
假设对平面匀速运动的景物采集一幅图像 f (x,y),并设 x0 (t) 和 y0 (t)
分别是景物在 x,y方向的运动分量；
g(x,y) =∫0T f [x- x0 (t),y- y0 (t)] dt ； T是采集时间长度；
g(x,y)即是由于运动而造成的模糊图像；
傅里叶变换 G(u,v)= F(u,v) H(u,v) ；
归结为无约束恢复类
H(u,v) = ∫0T exp[-j2?(u x0 (t) +v y0 (t))] dt；
如果知道了运动分量 x0 (t) 和 y0 (t),就可解析得到传递函数 H(u,v) 。
CHAPTER 5
IMAGE RESTORATION
and RECONSTRUCTION
§ 1 退化的数学模型和对角化
§ 2 无约束恢复
§ 3 有约束恢复
§ 4 几何失真校正
§ 5 投影重建
§ 5.1 退化的数学模型和对角化
退化：图像质量的降低；
失真可看作是退化，校正是恢复；
投影可看作是退化（三维到二维平面），重建是恢复；
§ 5.1.1 简单的通用图像退化模型
n（ x,y）
f（ x,y） H + g（ x,y）
模型化：一个作用在 f（ x,y） 上的系统 H与一个加性噪声 n（ x,y） 的联合作用,导致产生退化图像 g（ x,y） 。
假设已知 n（ x,y） 的统计特性（或先求出），图像复原就是已知 g
（ x,y） 求 f（ x,y） 的问题 （近似过程）。
g（ x,y） = H [f（ x,y） ] + n（ x,y） ；
已知 退化 解 噪声
§ 5.1.2 退化模型的计算一,1D情况：
设 f（ x） 和 h（ x） 均匀采样后放入尺寸为 A,B的数组；
f（ x）,x=0,1,…,A-1； h（ x）,x=0,1,…,B-1；
为了避免卷积的周期重叠，取 M? A+B-1；
将 f（ x） 和 h（ x） 用零扩展补齐； (M周期 )
ge（ x） =? fe（ m） h（ x-m） ； x,m=0,1,…,M-1；
用矩阵形式表示的卷积为 g（ x） = H f,展开得
ge（ 0） = he（ 0） he（ -1） … he（ -M+1） fe（ 0）
ge（ 1） = he（ 1） he（ 0） … he（ -M+2） fe（ 1）
ge（ M-1） = he（ M-1） he（ M-2） … he（ 0） fe（ M-1）
§ 5.1.2 退化模型的计算 （续 1）
式中,he（ 0） he（ -1） … he（ 1）
H = he（ 1） he（ 0） … he（ 2）
he（ M-1） he（ M-2） … he（ 0）
H矩阵是一个循环矩阵，其特点有：行列首尾相接；环形
循环矩阵相加后仍是循环矩阵，相乘后仍是循环矩阵；
§ 5.1.2 退化模型的计算 （续 2）
二,2D情况（ 1D的直接推广,0?x? M-1,0? y? N-1 ）
he(x,y)= h(x,y),0?x? A-1,0? y? B-1,否则 he(x,y)= 0；
fe(x,y)= f (x,y),0?x? C-1,0? y? D-1,否则 fe(x,y)= 0；
ge (x,y) = fe (m,n) he (x-m,y-n) ； m=0,1,…,M-1； n=0,1,…,N-1；
用矩阵形式表示的卷积为 g = H f,展开得
ge（ 0） = H0 HM-1 … H1 fe（ 0）
ge（ 1） = H1 H0 … H2 fe（ 1）
ge（ M-1） = HM-1 HM-2 … H0 fe（ MN-1）
其中每个块 Hi 是由扩展函数 he (x,y)的第 i行而来，是块循环矩阵。
§ 5.1.3 循环矩阵对角化一,1D对角化设 循环矩阵 H的本征值和特征向量分别 为：
(k) = [ 1 exp[j2?/M *k] …… exp[j2?/M* (M-1)k]T；
(k) = he(0)+he(M-1)exp[j2?/M?k]+…… + he(1)exp[j2?/M? (M-1) k] ；
将 H的 M个特征向量组成一个 M*M的矩阵 W：
W=[?(0)?(1)?(2) ……?(M-1)]
H = WDW-1；
式中 D为一个对角矩阵，其元素正是 H本征值，
即 D (k,k)=?(k)。
§ 5.1.3 循环矩阵对角化 (续 1)
设对 4*4的循环矩阵 C对角化。
设?(k)为整数 1的 M次根,k为根的序号，即?(k)M = 1 = lne；
因?(k)=exp[j2?k / M]； 例中 M=4； 即?(k)4 = 1；?(k)=exp[j2?k/4]；
规定标量?(k)= c0+?(k) c1 +?(k)2 c2 +?(k)3 c3= ∑ ci exp[j2?k/4*i]；
根据线性代数分析,k=0,1,2,3时可得到四个本征值和相应的特征向量；
(k)?(k)=C?(k)； 组成?(k)矩阵 W=[?(0)?(1)?(2)?(3)]
各元素?(k) =[ 1 exp (j2?/4?k ) exp (j2?/4?2k ) exp (j2?/4?3k ) ]；
当循环矩阵已知本征值和特征向量后，可以进行对角化；
设对角化矩阵为 D,则 D= W-1CW；（ C等于 H ）；
D为对角阵,D (k,k)=?(k) =MH (k) = W-1HW；
§ 5.1.3 循环矩阵对角化 (续 2)
二,2D对角化（块循环矩阵的对角化）
借助循环矩阵的结果进行推广；
设矩阵 W的尺寸为 MN*MN,每个元素为 W(i,m)=exp(j2?/M*im)WN ；
WN为一个 N*N矩阵，每个元素为 WN(k,n)=exp(j2?/N*kn)；
H = WDW-1 ；且 HT = WD*W-1 ； D*为 D的复共轭；
∴ D = W-1HW；可转化为对角阵。
§ 5.1.4 退化模型对角化的效果一,1D情况：
H=WDW-1 ； 等式代入 g=H f,并两边同左乘 W-1,得 W-1 g= DW-1 f，
W-1 g= DW-1 f 都是 M维向量，令它们分别为 G(k)和 F(k),得 G(k)=DF(k)；
其中,G(k)= (1/M) ∑ ge (i) exp(j2?/M*ki)是 ge (x)的傅里叶变换；
F(k)= (1/M) ∑ fe (i) exp(j2?/M*ki)是 fe (x)的傅里叶变换；
已知 D (k,k)=?(k) = ∑ he (i) exp(j2?/M*ki)=MH (k),
H (k)是 he (i)的傅里叶变换，
所以 G(k)= M× H (k) F(k),式中 H (k) F(k)为 fe (x) 与 ge (x)的卷积，
可采用 FFT计算。
§ 5.1.4 退化模型对角化的效果 (续 1)
二,2D情况：
G(u,v)= (1/MN) ∑∑ ge (x,y) exp[-j2? (ux/M+vy/N)] ；
F(u,v)= (1/MN) ∑∑ fe (x,y) exp[-j2? (ux/M+vy/N)] ；
N(u,v)= (1/MN) ∑∑ ne (x,y) exp[-j2? (ux/M+vy/N)] ；
H(u,v)= (1/MN) ∑∑ he (x,y) exp[-j2? (ux/M+vy/N)] ；
D的 MN个对角元素可用下式表示，
D(k,i)= MN *H ([k/N],k mod N)；当 i=k； 否则为 0；
其中 [k/N]表示取最大整数,mod表示取余数；
所以 G(u,v)= H(u,v) F(u,v) + N(u,v)； u,v=0,…,M -1
§ 5.1.5 有、无约束恢复一、无约束恢复公式
f ‘ = (nTn) -1 HT g = H -1 (HT) –1 HT g = H-1 g ；
二、有约束恢复公式若选取 f ‘ 的一个线性操作符 Q（ 变换矩阵），
使得 ||Q f ‘ ||最小，采用拉格朗日乘数法解决，最小化下式的 f ‘,
L（ f ‘ ） = ||Qf ‘ ||2+ l（ ||g-H f ‘ || 2 -||n|| 2 ），
l 为拉格朗日乘数，该式为准则函数；
最小化得有约束恢复公式：
f ‘ = [HT H + s QT Q ]-1 HT g； s = 1 / l ；
§ 5.2 无约束恢复
§ 5.2.1 逆滤波一、无噪声时处理方法
F ‘ (u,v) = G(u,v) / H(u,v) ； u,v=0,1,…,M-1
上式给定的恢复方法称为逆滤波，条件 H –1 (u,v) 存在。
H(u,v)取零或很小时，计算上较困难。
二、有噪声时处理方法（ G = H F + N,两端同除 H）
F ‘ (u,v) = F(u,v) + N(u,v) / H(u,v) ；
H(u,v)取零或很小时，计算结果差距很大；一般情况下，
逆滤波器并不正好是 1 / H(u,v),记为 M(u,v),称为恢复转移函数；
M(u,v) = 1 / H(u,v) ；如 u2+v2 ≤?02,?0为去掉零点的值；
M(u,v) = 1 ； 如 u2+v2 >?02,
§ 5.2.1 逆滤波 (续 1)
图像退化和恢复模型
n(x,y)
f (x,y) H(u,v) + g (x,y) M(u,v) f ‘ (x,y)
退化 恢复
改进,M (u,v) = k ； 如 H(u,v) ≤ d,d和 k为小于 1的常数；
M(u,v) = 1/ H(u,v) ；其它,d越小越好
H (u,v) 已知，则可恢复原图像；
退化系统的转移函数可以用退化图像的傅里叶变换来近似，即
H (u,v) ≈ G (u,v)
FFT 去零点 IFFT
g (x,y) G(u,v) ≈ H(u,v) M(u,v) f ‘ (x,y)
点冲激下 F(u,v)=1
§ 5.2.2 运动模糊图像的恢复
运动模糊图像的恢复主要是消除匀速直线运动模糊，采用照相机摄影。
假设对平面匀速运动的景物采集一幅图像 f (x,y),并设 x0 (t) 和 y0 (t)
分别是景物在 x,y方向的运动分量；
g(x,y) =∫0T f [x- x0 (t),y- y0 (t)] dt ； T是采集时间长度；
g(x,y)即是由于运动而造成的模糊图像；
傅里叶变换 G(u,v)= F(u,v) H(u,v) ；
归结为无约束恢复类
H(u,v) = ∫0T exp[-j2?(u x0 (t) +v y0 (t))] dt；
如果知道了运动分量 x0 (t) 和 y0 (t),就可解析得到传递函数 H(u,v) 。
