§ 5.3 有约束恢复
§ 5.3.1 维纳滤波 (Winner 滤波 )
Winner 滤波器是一种最小均匀误差滤波器；
维纳滤波是基于图像和噪声的相关矩阵的（统计方法），属于有约束恢复。
设 Rf和 Rn分别是原图像 f (x,y)和噪声图像 n(x,y)的相关矩阵，
Rf和 Rn均为实对称矩阵。
考虑一般图像像素间的相关不超过 20-30像素，故相关矩阵在主对角线方向有不为 0的条带,Rf和 Rn都可用块循环矩阵表达。
Rf = WAW-1 ； Rn = WBW-1 ； A,B类似于 H → D 对角化。
§ 5.3.1 维纳滤波 ()续 1)
定义 QTQ= Rf-1Rn,并代入一般公式 f ‘ = [HTH +s QTQ]-1 HT g；
得 f ‘ = [HTH +s Rf-1Rn ]-1 HT g
= [WD*DW-1 + sWA-1BW-1 ]-1 WD*W-1 g；
整理得
F‘ (u,v)={(1/ H(u,v) )*|H(u,v)|2/ [| H(u,v)|2 + s [Sf (u,v)/Sn(u,v) ] }*
G(u,v)；
讨论（ 1） s=1,{·}中的项就是维纳滤波器；
（ 2） s是变量，称为参数维纳滤波器；
（ 3）无噪声时,Sn(u,v)=0； 原式 =1/ H(u,v)，
即维纳滤波器成为了理想逆滤波器。
若 Sf (u,v)和 Sn(u,v) 未知时，可近似计算：
F‘ (u,v) ≈ { (1/ H(u,v) )*|H(u,v)|2/ [| H(u,v)|2 + K ] }* G(u,v)；
式中,K为预先设定的常数；
比较：维纳滤波需要事先知道统计量 Rf和 Rn，在图像受噪声影响时，恢复效果比逆滤波要好，而且噪声越大，优势越明显。
§ 5.3.2 有约束最小平方恢复有约束最小平方恢复也称为最小乘方滤波恢复；
有约束最小平方恢复方法只需要有关噪声均值和方差即可。
将 pe(x,y)扩展，得
pe(x,y)= p(x,y)； 当 0≤x≤2和 0≤y≤2时
pe(x,y)=0； 当 3≤x≤M-1和 3≤y≤N-1时用矩阵形式表示，构造分块循环矩阵 C；
用 W进行对角化，有 E = W-1CW； E 是一个对角矩阵。
对角矩阵元素 E(k,i) = P([k/N],k mod N) ； 如 i=k,否则为 0；
其中,P(u,v)是 pe(x,y)的 2D傅里叶变换。
§ 5.3.2 有约束最小平方恢复（续 1）
求满足约束条件 ||g-H f ‘ ||2 = || n ||2 的最优解，得
f ‘ = [HTH +s CT C ]-1 HT g
= [WD*DW-1 + sWE*EW-1]-1 WD*W-1 g；
Wf ‘ = [D*D + sE*E]-1 D*W-1 g；（ 简化计算，转移到频域）
F‘ (u,v) = H*(u,v) / [| H(u,v)|2 + s P(u,v)2] G(u,v)；
可见与维纳滤波器相似，但仅与噪声均值和方差有关。
与维纳滤波的比较：
有模糊又有噪声时，有约束最小平方滤波的效果较好；
仅有模糊没有噪声时，两者基本一致。
§ 5.4 几何失真校正
几何失真在广义上也是一种图像退化，校正就是一种图像恢复的过程。
几何失真校正主要包括二个步骤：
空间变换，恢复像素原空间的位置；
灰度插值，对空间变换后像素赋予原位置的灰度值；
§ 5.4.1 空间变换
按照一幅标准图像或一组基准点来校正一幅几何失真图像；
设基准（原）图像为 f (x,y),几何形变图像为 g (x’,y’) ；
x’,y’为失真图像的坐标，即 x’= s(x,y),y’= t (x,y) ；
s (x,y)和 t (x,y) 代表空间变换；若已知 s 和 t,就可以通过反变换来恢复图像。
§ 5.4.1 空间变换（续 1）
一、线性失真（三角形线性法）
s (x,y) = k1x+ k2y + k3 或 ax + by + c；
t (x,y) = k4x+ k5y + k6 或 dx + ey + f；
将小三角形区域三个顶点作为控制点,3组（ 6个）已知点建立方程组，
xi’ = k1xi+ k2yi + k3,i=1,2,3
yi’ = k4xi+ k5yi + k6,i=1,2,3
由待定系数法，求出 k1,k2,…,k6，即可实现三角形内各像素的校正。
关键,找出控制点的位置；
二、双线性等式（四边形区域）
s (x,y) = k1x+ k2y + k3xy + k4 ；
t (x,y) = k5x+ k6y + k7xy + k8 ；
§ 5.4.1 空间变换（续 2）
将四边形区域的顶点作为控制点,4组（ 8个）已知点建立方程组，
xi’ = k1xi+ k2yi + k3 xiyi + k4,i=1,2,3,4
yi’ = k5xi+ k6yi + k7 xiyi + k8,i=1,2,3,4
由待定系数法，求出 k1,k2,…,k8 8个系数，这些系数可映射到四边形内的所有点。
三、二元多项式法将原图像的空间坐标和被校正图像的空间坐标之间的关系用下式描述：
xi’ = ∑∑aij xi yi = a00+a01y+a02y2+a10x+a11xy+a20x2
yi’ = ∑∑bij xi yi = b00+b01y+b02y2+b10x+b11xy+b20x2
二元多项式中,aij 和 bij 为待定系数,i=0,…,n; j=0,…,n-i；
同线性法类似，采用已知的控制点对，用曲面拟合方法求解系数，补充拟合误差平方和的一阶导数为 0的条件，构成线性方程组，求出 aij 和 bij,变换。
§ 5.4.2 灰度插值
灰度插值确定校正点处的亮度，以内插法为主。
一、最近邻插值（零阶插值）
将像素邻点中距离最近的邻点灰度作为该点的灰度。
方法简单，但不够准确。
二、双线性插值
计算公式 g(i+u,j+v)= (1-u) (1-v) g(i,j)+ (1-u) v g(i,j+1)+
u (1-v) g(i+1,j) + uv g(i+1,j+1)；
§ 5.5 投影重建
§ 5.5.1 概述一、投影重建定义投影重建指从一个物体的多个（轴向）投影图重建目标图像的过程。输入是一系列图像，输出是一幅重建图。
投影重建包括：
透射断层成像 TCT,简称 CT,transmission computer
tomography； 对大脑不同层面进行多次扫描，合成图像。
发射断层成像 ECT,Emission CT；
反射断层成像 RCT,Reflection CT；
核磁共振成像 MRI,Magnetic Resonance Imaging
§ 5.5.1 概述（续 1）
二、投影重建图像模型设重建图像 f (x,y) 在一个以原点为圆心的单位圆 Q外为零。
考虑有一条由发射源到接受器的直线在平面上与 f (x,y)在 Q内相交，这条直线可用 2条参数来确定：
1）它与原点的距离 s,2） 它与 X轴的夹角?；
用 g(s,?)表示沿直线 (s,?) 对 f (x,y) 积分，经坐标变换有
g(s,?)=∫f (x,y) dt ；
= ∫f (s * cos? - t * sin?,s * sin? - t * cos?) dt
这个积分结果含义就是 f (x,y) 沿 t方向的投影，积分限取决于 s,?
和 Q,若?转一个 360度,s从 0到边缘，整个图形全部被投影了。
§ 5.5.1 概述（续 2）
当 Q是一个单位圆时,s2 + t 2=1；
因此，用 f (x,y)表示要被重建的目标，由 (s,?)确定的积分路线对应一条从发射器到接受器的射线，接受器所获得的积分测量值是 g(s,?)。
对给定的 g(s,?)，要确定 f (x,y)，数学上讲，就是解 g(s,?)方程。
三、重建方法分类：
1,变换重建方法：先在连续域解析处理，在离散化计算；
包括傅里叶反变换重建法，卷积逆投影重建法
2,级数展开重建法：离散化分析，直接得到数字解；
3,综合重建法：
§ 5.5.2 傅里叶反变换重建
§ 5.5.2.1 基本步骤基本步骤
1）建立数学模型；已知量、未知量均为连续实数的函数；
2）利用反变换公式解未知量 f (x,y) ；
3）调节反变换公式，以适应离散、有噪声应用的要求；
测量的数据：许多个离散点 (s,?)上的估计值 g(s,?) ；
重建的图像：离散的 2D数组；
受测量离散的近似和等价公式的不同，计算结果会不同；
§ 5.5.2.2 傅里叶变换投影定理
傅里叶变换投影定理是变换方法的基础。
设图像区被一个直角网格覆盖,X方向的点数从 K-到 K+,Y方向的点数从 L-到 L+,重建算法就是要通过 M*N个测量值估计出在 K*L个采样点的 f (k△ x,l△ y) ；
设 G(R,?)是 g(s,?)对应第一个变量 s的一维傅里叶变换，即
G(R,?)= ∫g(s,?) exp[-2j?Rs] ds； 沿?角的直线上的投影的 FT
设 F(X,Y)是 f (x,y)的二维傅里叶变换，即
F(X,Y) = ∫∫f (x,y) exp[-2j?(xX+yY) ] dxdy；
则 G(R,?)= F(Rcos?,Rsin?) ；证明过程略。
含义,f (x,y)在与 X轴成? 角的直线上投影的傅里叶变换 G(R,?)是 f
(x,y)的傅里叶变换 F(X,Y)在朝向角?上的一个截面。
§ 5.5.2.3 傅里叶反变换重建公式
对投影定理公式两端在直角坐标系中取傅里叶反变换，就得到了 f (x,y)
的重建公式；
f (x,y)=∫∫G((X2+Y2)1/2,arcty(Y/X)) exp[-2j?(xX+yY) ] dXdY；
即由 g(s,?)进行傅里叶变换，到 G(R,?)，反变换到 f (x,y)；
G(·)的计算：
G(R,?)用一系列采样点 (m△ s,?n)对 g(·)求和得到，
G∑(R,?)= △ s∑g(m△ s,?n) exp[-2j?R(m△ s)] ；
令 R=k△ R； △ R为采样间距,k为整数，△ R=1/ (m△ s)；
则 G∑(△ R,?)= △ s∑g(m△ s,?n) exp[-2j?km/M)] ；
§ 5.5.2.3 傅里叶反变换重建公式（续 1）
反变换得
f W(k△ x,k△ y) = △ x△ y ∑∑FW(u△ x,v△ y) exp[-2j?(ku/U+lv/V) ] ；
综合起来，基于傅里叶反变换的重建技术主要步骤：
（ 1）对角度?n ； n=0,1,…,N的投影进行一维傅里叶变换，获得极坐标下的采样值；
（ 2）通过内插估计直角坐标下的采样值；
（ 3）进行二维傅里叶反变换求得重建图像；
优点：计算量小；
缺点：需要二维插值实现不易，重建图像质量较差；
§ 5.5.2.4 卷积逆投影重建
滤波逆投影重建
一、公式
由傅里叶变换投影定理，在极坐标系中取 G(R,?)= F(Rcos?,Rsin?)
的傅里叶反变换，得
f (x,y)= ∫∫ G(R,?) exp[-2j?R(xcos? + ysin?) ] |R| dR d?；
引入加窗对 G(R,?)进行估计，得
f W(x,y)= ∫ g’ (xcos? + ysin?,?) d?；
这就是卷积逆投影重建公式，表明将各个?角方向的投影先滤波，然后将不同方向投影的滤波结果积分，从而得到重建图像。
§ 5.5.2.4 卷积逆投影重建（续 1）
公式一部分是 f (x,y)在?角方向的投影与 h(s) 的卷积 ； h(s) 称为卷积函数；
公式另一部分是一条以?角通过 (x,y)点的射线；
二、计算的实现
f W(k△ x,l△ y) ≈△?∑g’ (k△ xcos?n + l△ ysin?n,?n) ； 近似计算公式
离散化简化计算,g’ (m△ s,?n),M-≤m≤M+,根据 M个 g’ 值按插值获得 K*L个 g’ 值，再求和得到 f W ；
即 gc’(m’△ s,?n) ≈ △ s∑g(m△ s,?n) h[(m’- m )]△ s ； 离散卷积
gI’(s’,?n) ≈ △ s∑ gc’(m△ s,?n) I(s’ - m△ s ) ； 一次插值
§ 5.3.1 维纳滤波 (Winner 滤波 )
Winner 滤波器是一种最小均匀误差滤波器；
维纳滤波是基于图像和噪声的相关矩阵的（统计方法），属于有约束恢复。
设 Rf和 Rn分别是原图像 f (x,y)和噪声图像 n(x,y)的相关矩阵，
Rf和 Rn均为实对称矩阵。
考虑一般图像像素间的相关不超过 20-30像素，故相关矩阵在主对角线方向有不为 0的条带,Rf和 Rn都可用块循环矩阵表达。
Rf = WAW-1 ； Rn = WBW-1 ； A,B类似于 H → D 对角化。
§ 5.3.1 维纳滤波 ()续 1)
定义 QTQ= Rf-1Rn,并代入一般公式 f ‘ = [HTH +s QTQ]-1 HT g；
得 f ‘ = [HTH +s Rf-1Rn ]-1 HT g
= [WD*DW-1 + sWA-1BW-1 ]-1 WD*W-1 g；
整理得
F‘ (u,v)={(1/ H(u,v) )*|H(u,v)|2/ [| H(u,v)|2 + s [Sf (u,v)/Sn(u,v) ] }*
G(u,v)；
讨论（ 1） s=1,{·}中的项就是维纳滤波器；
（ 2） s是变量，称为参数维纳滤波器；
（ 3）无噪声时,Sn(u,v)=0； 原式 =1/ H(u,v)，
即维纳滤波器成为了理想逆滤波器。
若 Sf (u,v)和 Sn(u,v) 未知时，可近似计算：
F‘ (u,v) ≈ { (1/ H(u,v) )*|H(u,v)|2/ [| H(u,v)|2 + K ] }* G(u,v)；
式中,K为预先设定的常数；
比较：维纳滤波需要事先知道统计量 Rf和 Rn，在图像受噪声影响时，恢复效果比逆滤波要好，而且噪声越大，优势越明显。
§ 5.3.2 有约束最小平方恢复有约束最小平方恢复也称为最小乘方滤波恢复；
有约束最小平方恢复方法只需要有关噪声均值和方差即可。
将 pe(x,y)扩展，得
pe(x,y)= p(x,y)； 当 0≤x≤2和 0≤y≤2时
pe(x,y)=0； 当 3≤x≤M-1和 3≤y≤N-1时用矩阵形式表示，构造分块循环矩阵 C；
用 W进行对角化，有 E = W-1CW； E 是一个对角矩阵。
对角矩阵元素 E(k,i) = P([k/N],k mod N) ； 如 i=k,否则为 0；
其中,P(u,v)是 pe(x,y)的 2D傅里叶变换。
§ 5.3.2 有约束最小平方恢复（续 1）
求满足约束条件 ||g-H f ‘ ||2 = || n ||2 的最优解，得
f ‘ = [HTH +s CT C ]-1 HT g
= [WD*DW-1 + sWE*EW-1]-1 WD*W-1 g；
Wf ‘ = [D*D + sE*E]-1 D*W-1 g；（ 简化计算，转移到频域）
F‘ (u,v) = H*(u,v) / [| H(u,v)|2 + s P(u,v)2] G(u,v)；
可见与维纳滤波器相似，但仅与噪声均值和方差有关。
与维纳滤波的比较：
有模糊又有噪声时，有约束最小平方滤波的效果较好；
仅有模糊没有噪声时，两者基本一致。
§ 5.4 几何失真校正
几何失真在广义上也是一种图像退化，校正就是一种图像恢复的过程。
几何失真校正主要包括二个步骤：
空间变换，恢复像素原空间的位置；
灰度插值，对空间变换后像素赋予原位置的灰度值；
§ 5.4.1 空间变换
按照一幅标准图像或一组基准点来校正一幅几何失真图像；
设基准（原）图像为 f (x,y),几何形变图像为 g (x’,y’) ；
x’,y’为失真图像的坐标，即 x’= s(x,y),y’= t (x,y) ；
s (x,y)和 t (x,y) 代表空间变换；若已知 s 和 t,就可以通过反变换来恢复图像。
§ 5.4.1 空间变换（续 1）
一、线性失真（三角形线性法）
s (x,y) = k1x+ k2y + k3 或 ax + by + c；
t (x,y) = k4x+ k5y + k6 或 dx + ey + f；
将小三角形区域三个顶点作为控制点,3组（ 6个）已知点建立方程组，
xi’ = k1xi+ k2yi + k3,i=1,2,3
yi’ = k4xi+ k5yi + k6,i=1,2,3
由待定系数法，求出 k1,k2,…,k6，即可实现三角形内各像素的校正。
关键,找出控制点的位置；
二、双线性等式（四边形区域）
s (x,y) = k1x+ k2y + k3xy + k4 ；
t (x,y) = k5x+ k6y + k7xy + k8 ；
§ 5.4.1 空间变换（续 2）
将四边形区域的顶点作为控制点,4组（ 8个）已知点建立方程组，
xi’ = k1xi+ k2yi + k3 xiyi + k4,i=1,2,3,4
yi’ = k5xi+ k6yi + k7 xiyi + k8,i=1,2,3,4
由待定系数法，求出 k1,k2,…,k8 8个系数，这些系数可映射到四边形内的所有点。
三、二元多项式法将原图像的空间坐标和被校正图像的空间坐标之间的关系用下式描述：
xi’ = ∑∑aij xi yi = a00+a01y+a02y2+a10x+a11xy+a20x2
yi’ = ∑∑bij xi yi = b00+b01y+b02y2+b10x+b11xy+b20x2
二元多项式中,aij 和 bij 为待定系数,i=0,…,n; j=0,…,n-i；
同线性法类似，采用已知的控制点对，用曲面拟合方法求解系数，补充拟合误差平方和的一阶导数为 0的条件，构成线性方程组，求出 aij 和 bij,变换。
§ 5.4.2 灰度插值
灰度插值确定校正点处的亮度，以内插法为主。
一、最近邻插值（零阶插值）
将像素邻点中距离最近的邻点灰度作为该点的灰度。
方法简单，但不够准确。
二、双线性插值
计算公式 g(i+u,j+v)= (1-u) (1-v) g(i,j)+ (1-u) v g(i,j+1)+
u (1-v) g(i+1,j) + uv g(i+1,j+1)；
§ 5.5 投影重建
§ 5.5.1 概述一、投影重建定义投影重建指从一个物体的多个（轴向）投影图重建目标图像的过程。输入是一系列图像，输出是一幅重建图。
投影重建包括：
透射断层成像 TCT,简称 CT,transmission computer
tomography； 对大脑不同层面进行多次扫描，合成图像。
发射断层成像 ECT,Emission CT；
反射断层成像 RCT,Reflection CT；
核磁共振成像 MRI,Magnetic Resonance Imaging
§ 5.5.1 概述（续 1）
二、投影重建图像模型设重建图像 f (x,y) 在一个以原点为圆心的单位圆 Q外为零。
考虑有一条由发射源到接受器的直线在平面上与 f (x,y)在 Q内相交，这条直线可用 2条参数来确定：
1）它与原点的距离 s,2） 它与 X轴的夹角?；
用 g(s,?)表示沿直线 (s,?) 对 f (x,y) 积分，经坐标变换有
g(s,?)=∫f (x,y) dt ；
= ∫f (s * cos? - t * sin?,s * sin? - t * cos?) dt
这个积分结果含义就是 f (x,y) 沿 t方向的投影，积分限取决于 s,?
和 Q,若?转一个 360度,s从 0到边缘，整个图形全部被投影了。
§ 5.5.1 概述（续 2）
当 Q是一个单位圆时,s2 + t 2=1；
因此，用 f (x,y)表示要被重建的目标，由 (s,?)确定的积分路线对应一条从发射器到接受器的射线，接受器所获得的积分测量值是 g(s,?)。
对给定的 g(s,?)，要确定 f (x,y)，数学上讲，就是解 g(s,?)方程。
三、重建方法分类：
1,变换重建方法：先在连续域解析处理，在离散化计算；
包括傅里叶反变换重建法，卷积逆投影重建法
2,级数展开重建法：离散化分析，直接得到数字解；
3,综合重建法：
§ 5.5.2 傅里叶反变换重建
§ 5.5.2.1 基本步骤基本步骤
1）建立数学模型；已知量、未知量均为连续实数的函数；
2）利用反变换公式解未知量 f (x,y) ；
3）调节反变换公式，以适应离散、有噪声应用的要求；
测量的数据：许多个离散点 (s,?)上的估计值 g(s,?) ；
重建的图像：离散的 2D数组；
受测量离散的近似和等价公式的不同，计算结果会不同；
§ 5.5.2.2 傅里叶变换投影定理
傅里叶变换投影定理是变换方法的基础。
设图像区被一个直角网格覆盖,X方向的点数从 K-到 K+,Y方向的点数从 L-到 L+,重建算法就是要通过 M*N个测量值估计出在 K*L个采样点的 f (k△ x,l△ y) ；
设 G(R,?)是 g(s,?)对应第一个变量 s的一维傅里叶变换，即
G(R,?)= ∫g(s,?) exp[-2j?Rs] ds； 沿?角的直线上的投影的 FT
设 F(X,Y)是 f (x,y)的二维傅里叶变换，即
F(X,Y) = ∫∫f (x,y) exp[-2j?(xX+yY) ] dxdy；
则 G(R,?)= F(Rcos?,Rsin?) ；证明过程略。
含义,f (x,y)在与 X轴成? 角的直线上投影的傅里叶变换 G(R,?)是 f
(x,y)的傅里叶变换 F(X,Y)在朝向角?上的一个截面。
§ 5.5.2.3 傅里叶反变换重建公式
对投影定理公式两端在直角坐标系中取傅里叶反变换，就得到了 f (x,y)
的重建公式；
f (x,y)=∫∫G((X2+Y2)1/2,arcty(Y/X)) exp[-2j?(xX+yY) ] dXdY；
即由 g(s,?)进行傅里叶变换，到 G(R,?)，反变换到 f (x,y)；
G(·)的计算：
G(R,?)用一系列采样点 (m△ s,?n)对 g(·)求和得到，
G∑(R,?)= △ s∑g(m△ s,?n) exp[-2j?R(m△ s)] ；
令 R=k△ R； △ R为采样间距,k为整数，△ R=1/ (m△ s)；
则 G∑(△ R,?)= △ s∑g(m△ s,?n) exp[-2j?km/M)] ；
§ 5.5.2.3 傅里叶反变换重建公式（续 1）
反变换得
f W(k△ x,k△ y) = △ x△ y ∑∑FW(u△ x,v△ y) exp[-2j?(ku/U+lv/V) ] ；
综合起来，基于傅里叶反变换的重建技术主要步骤：
（ 1）对角度?n ； n=0,1,…,N的投影进行一维傅里叶变换，获得极坐标下的采样值；
（ 2）通过内插估计直角坐标下的采样值；
（ 3）进行二维傅里叶反变换求得重建图像；
优点：计算量小；
缺点：需要二维插值实现不易，重建图像质量较差；
§ 5.5.2.4 卷积逆投影重建
滤波逆投影重建
一、公式
由傅里叶变换投影定理，在极坐标系中取 G(R,?)= F(Rcos?,Rsin?)
的傅里叶反变换，得
f (x,y)= ∫∫ G(R,?) exp[-2j?R(xcos? + ysin?) ] |R| dR d?；
引入加窗对 G(R,?)进行估计，得
f W(x,y)= ∫ g’ (xcos? + ysin?,?) d?；
这就是卷积逆投影重建公式，表明将各个?角方向的投影先滤波，然后将不同方向投影的滤波结果积分，从而得到重建图像。
§ 5.5.2.4 卷积逆投影重建（续 1）
公式一部分是 f (x,y)在?角方向的投影与 h(s) 的卷积 ； h(s) 称为卷积函数；
公式另一部分是一条以?角通过 (x,y)点的射线；
二、计算的实现
f W(k△ x,l△ y) ≈△?∑g’ (k△ xcos?n + l△ ysin?n,?n) ； 近似计算公式
离散化简化计算,g’ (m△ s,?n),M-≤m≤M+,根据 M个 g’ 值按插值获得 K*L个 g’ 值，再求和得到 f W ；
即 gc’(m’△ s,?n) ≈ △ s∑g(m△ s,?n) h[(m’- m )]△ s ； 离散卷积
gI’(s’,?n) ≈ △ s∑ gc’(m△ s,?n) I(s’ - m△ s ) ； 一次插值
