摘要和正文在实际应用当中设备捕获到频率范围是受成像传感器最大采样

成本和元件的物理特性极限制约了大多数传统的成像设备提高分辨率的成像能力。为了处理这样的问题，众多图像处理领域的研究正在致力于开发一本文首先通过对超分辨率重建适用的成像模型和CT成像设备的点扩散函数CT的超分辨率重建的数学模型。本文将基于变换性质的相位相关算法应用在T图像的图像配准详尽阐述了相位相关算法获得图像之间的平移偏移量的方法和借助极坐标系将图T图像的图像配1地通过减小传感器矩阵的尺寸大小和增加传感器的数目，这不是切实可行的办但是伴随着像素尺寸的缩小而来的是信噪比（ignaltonoisertio，S）和元。在实际应用当中，设备捕获到的频率范围是受成像设备传但是这样的代价是的成像元件实际的成本和元件的物理特性极限制约了大多数传统的成像设备提用最广泛的医学影像设备主要是磁成像（MRI）和X线计算机断层成（T处理中，与分辨率相关的需要去攻克的难题被认为是部分容积效应（prtilvolumeftPET分辨率的设备进行成像的时候，不可避免地会产生部分容积效应。通常，PI示的是与来自邻近组织的信号平均成正比的灰度值这样可能导致边界确的位置可能会有偏移，对脑部I解剖图像的分割会带来。同向分辨率TIPE加能有效克服或减少与PE相关的不良影响。医疗成像领域在拥有更高的分辨率的机器和高级图像内容处理工具方面发展很快医疗成像系统的开发致力于提高分辨率因为更高的分辨率是更准理解解剖结构的关键它可以支持的早期发现和增加对于和病理大和形态估计的准确度。20012002年期间[1~7]，超分辨率算法开始被尝试从计算机视觉领域行重现，而且在衍生的设备上如功能磁成像（functionalmagneticresonanceimaging，fMRI）和正电子发射型计算机断层显像（positronemissioncomputedtomography，PET）MRIPET设Stefan将基于MAP（最大后验概率）的方法应用到CTMRI上重建3D影CT图像进行超分辨率重建。本文共分为六个章节进行论述。第2章CT图像的超分辨率重T成像是指利用射线束对特定的检查层面进行扫描通过各个不同方向的探测器接收的组织对线的衰减值，经过模/数转换后将数据输入计算在显示器上显示出来。T图像的内容是特定的检查层面的组织密度分布。T图像具有较高的清晰度和密度分辨率，并且不会受到层面以外的组织结构的面，这时反复使用较高能量的射线进行成像会对患者的造成一定的射线损超分辨率重建技术是指利用一组低分辨率图像之间的包含的信息来重建生重建方法[8~17]。2003年，SungholPrk方法被HanochUr和DanielGross提出。空间域的方法进展很快，Irani和Peleg了一种基于空间域的反向投影迭代（iterativebackprojection，IBP）的方法，不给出一组低分辨率图像{yk，并且假设{ykyk=fk(x)hS+vk;k

是高分辨率原始图像x降质生成低分辨率图像时产生的相对基准图像的几何数（PSF。最后vk代表的是加性噪声。符号是卷积运算符S代表的是从接下来本节将要给出的是超分辨率的一般表达形式，是一种与Irani和的IBP方法有关的形式，同样也是接近于ML的形式。式（2.2）中，x(0)作为高 1 k fkk1k

(ykS

yk是K个采样结果之一，fk是对相同参考帧 时产生的几何变形，S代表式（2.3）所示的是超分辨率重建使用的图像成像模型，其中x(0)代表了一幅k次迭代反投影之后得到的高分辨率估计图像，y(n)表示由x(n)经过降质模型k处理，还有使用系统的采样频率进行下采样时产生的效果。就这样，一幅nx(n经过降质模型（忽略噪声）得到了一组n幅低分辨 y(n)

(x(n))h

k根据已知的降质模型产生的低分辨率图像集y(n)和真实到的低分辨率图像集{yk之间的差异，通过如式（2.4）所示的更新过程叠加到当前的高分辨kk 1k

k1

(((yk

)s)

p=hK2(x)Kk

yk

yyk2（2.1（2.3）和式（2.4描述了最大似然估计量在式描述的2范数准则下的一个解，并且假设在这个模型里面的加性噪声是服从正Friu等人已经推导出一般的p1范数解的鲁棒性更好。2

(x(n))KkK

y y

2

A(x)是包含了所期望的超分辨率图像x是正则化系数，每一种超分辨率算法需要根据其最接近真实成像系统特征和具体应用场景入图像参考的是同一幅低分辨率图像在医疗成像设备中典型的变换是床）和成像设备之间相对于初始位置的物理位移。函数（PointSpreadFunction，PSF。在成像系统中，理想的图像应该是由组成在医疗成像系统中通常有三种确定点扩散函数的方法一是用条带测试卡直接测量调制传递函数但在实际应用中难以制造出恰当的实体模型效果不好二是把有机玻璃块放在水模中根据其轮廓函数计算调制传递函数实际应用中临界效应会使得轮廓图像改变从而产生误差三是用金属丝模型确定通常确定T成像系统的点扩散函数就是用此方法将射线衰减能力很高的金属材（比如钨）制造成的细丝作为信号模型，竖直放置在T中运行扫描。由于金属丝很细每个断面近似一个点光源物体因此可以将细丝断面的成像近似为系统的脉冲响应，根据这个可以确定T成像系统的点扩散函数。T成像系统的点扩散函数与 函数十分接近。理想的点扩散函数的逆是式（2.4）的BP核函数php是收敛的，式中h

比（PeakSignaltoNoiseRatio，PSNR）和结构相似性（structuralsimilarityindex，I1I2，图像的行列数分别是i和jc个通道。式定义的是两幅图像之间的均方误差。式（2.8）定义的峰值信噪比。其中，MAXI字节，所以MAXI的取值一般为255。PSNR的值通常大30到50之间，数MSE II2ij MAX2

PSNR10log10 IMSE

20

MAXI10MSE 结构相似 X-X- ÷Y-+÷xySSIM定义如式（2.9）所示。其中lx

cxysxy为结构对比函数。α、βγ

lx,y

2xy

,ck

22 式（2.10）x为xy为yL表示像素值的取值范围，通常是255。k1为很小的数，通常为0.01cx,y 2xy2222

,c2k2L

式（2.11）定义的是对比度对比函数。其中，xxyy的方差。xy为式（2.12）定义的x和yk20.03。 1M

N

i,

sx,y2xy

x 式（2.13）α=1β=1γ=1和c3c2/2由此得到式（2.14）表示的SSIM2xyc12xyc2 1 2SSIMx,y 1 2

22c22c

3章图像配图像配准就是寻找两幅图像之间存在的关系[20,21]用变换模型将两幅图（3.1）定义的是刚体变换模型。其中，定义点(x0,y0)是基准图像I0(x,y)中的一点，点(x1,y1)是待配准图像I1(x,y)中与点(x0,y0)对应的一点，相对旋转偏移角度为θsxsy分别表示相对水平平移量和竖直平移量。图（3.1）表示的是参数为θ=45sx=3.6sy=4.3时图像的刚体变换结果。x1cossinx0sx

y cosy s1 0 y 3.1射变换是由平移、缩放、旋转和剪切的一种组合变换。式（3.2）变换模型，其中a00、a01、a10、a11、sxsy

a01x0

sxy a

s

1 11

0 ysx=0sy=42.24A）变换前图像B）变换后图像3.2仿射变换的示意图平行或者夹角关系。式（3.3）定义了一种二中的非线性变换模型，其中G表示任意一种函数。(x1,y1)G(x0,y0

算法得到待配准的图像配准结果基于变换域法通常指的是基于变换性质变换有许多成快速算法可以应用，这类方法是图像配准中常用的方确度较高、计算速度较快的基于变换的相位相关算法进行图像配准。和Y方向上的平移不会改变变换的幅值，但是变换之间的相位差与这两幅相似图像的空间位移有关。式（3.4）I1(x,yI2x,y

s(x0,y0)的关系。式（3.5）表示这两幅图像的变换存在关系，其中F1(u,v)和F2(u,v)分别I1(x,y)和I2(x,y)的变换I2(x,y)I1(xx0,yy0 F(u,v)F(u,v)e2j(ux0vy0

F*uvF(uv)的共轭。F(u,v)FF(u,v)F式（3.8）定义一个R（ratio）表示两幅图像变换之间的一种关系F(u,v)F*(u,v)R1 F(u,v)F*(u,

（3.7s(x0y0与R有关。F(u,v)F*(u,v)R1 F(u,v)F*(u, F(u,v)F*(u,v)e2j(ux0vy0 F(u,v)F*(u,v)e2j(ux0vy0 F(u,v)F*(u,v)e2j(ux0vy0 F(u,v)F*(u, e2j(ux0vy0

根据ReddyChatterji的研究[25]，式（3.9）定义的R有特殊的含义R函数，这个点代表了R的逆变换中绝对值最大的那个点，这个点的坐(x0,y0就是对基准图像最佳的配准结果。相位相关方法的关键步骤就是求取R的逆变换，搜索得到R的逆变换中最大值的坐标(x0,y0)用于确定s(x0y0。式（3.10）I1(x,y)I2(x,ys(x0,y0)还有旋其中F1(u,v)和F2(u,v)分别I1(x,y)和I2(x,y)的变换。I2(x,y)I1(xcos0ysin0x0,-xsin0ycos0y0

210F(u,v)F(u210

v

)e2j(ux0vy0

00式（3.12）表示的是这两幅图像的变换的幅值之间存在的关系，其中00M1(uvM2(uvF1(uvF2(uvM2(u,v)M1(ucos0vsin0,usin0vcos0

F1(uvF2uv)的幅值关系，从式（3.12）转换到对数极坐标系中，旋转偏移量就能有像平移在坐标系中类似的表示。式（3.13）表示的两幅图像的变的幅值在对数极坐标系中的形式M1(log(),)M2(log(),0

间的旋转偏移量θ。I1(x,y)I2x,ys(x0,y0和旋转偏移量θI1(x,yI2(x,yg1(m1m2g2(m1m2，g1(m1m2g2(m1m2)的FFTF1(uvF2(uv)；F1(uvF2(uv的绝对值M1(uv和M2(uv)M1(uv)M2(uv)Flp1(uvFlp2(uv根据式（3.8）求得Flp1(uvFlp2(uvR1计算R1的IFFT得到IR1IR1abs(IR1中最大值所在的坐标(log(),kI2(x,yI1(x,y的旋转偏移角度g2(m1m2反向旋转kg3(m1m2g3(m1m2的FFTF3(uv根据式（3.8）求得F1(uvF3(uvR2R2的IFFTIR2IR2abs(IR2中最大值所在的位置(x0,y0)就可以得到亚像素相对平移量s(x0y0。 3.3用于测试相位相关算法的图像集是由一幅高分辨率图像经过降质得到的一128*128所示。表（.1）所示的是输入图像之间的实际偏移量和相位相关算法得到的图（（3.1所示结果统计得到第一组图像配准结果的平均误差为0.1825，总体的方差为0.0136，标准差为0.167。3.1水平方向（像素000竖直方向（像素000旋转角度000水平方向（像素竖直方向（像素旋转角度水平方向（像素竖直方向（像素旋转角度水平方向（像素竖直方向（像素旋转角度 3.4第二组配准实验采用的是胸部CT256*256，0.1100.0849，标0.1098。3.2水平方向（像素000竖直方向（像素000旋转角度000水平方向（像素竖直方向（像素旋转角度水平方向（像素竖直方向（像素旋转角度水平方向（像素竖直方向（像素旋转角度 3.5第三组配准实验采用的是腹部CT256*256，0.09750.033，标0.1078。3.3水平方向（像素000竖直方向（像素000旋转角度000水平方向（像素竖直方向（像素旋转角度水平方向（像素0竖直方向（像素0旋转角度水平方向（像素竖直方向（像素旋转角度CT0.2个像CTCT4像插图像插值解决的是根据未给出像素值的像素点附近的已给出像素值的像素[2627]第三，计算出来的估计值与给出的采样数据属于同续的函数模型。最近邻插值法是一种简单易实现的多变量插值（multivariateinterpolation）在图（4.1）中点(uv是未知的像素点。在点(uv附近一共有四个已知像素点分别是左上方的(i,j)、右上方的(ij1、左下方的(i1j)和右下方的(i1j1，它们对应的灰度值分别是g(i,j)、g(i,j1、g(i1,j)和g(i1j1。通过计算未知像素点(uv与周围附近四个已知像素点点(i,j)(ij1(i1j)和(i1,j1点(uvD(uv),(i,

为点(uv与点(i,j)之间距离，式（4.1）示的是未知的像素点(uv与周围附近四个已知像素点的最近距离SSminD(u,v),(i,j),D(u,v),(i,j1),D(u,v),(i1,j),D(u,v),(i1,4.1欧氏距离、街区距离和棋盘距离等都可以用来表示平面空间上两个点的距(ui)2(v(ui)2(v3D实时渲染中对表面纹理的插值。在数字图像处理中双线性插值通过对未知像素点点附近的2*2区域4个邻4个邻近像素的示意图。4.24从图（4.2）中可以看到，未知像素点p(uv与a(x1y2、b(x1y1、c(x2y2)和d(x2y1)四个已知像素点相邻a(x1,y2)b(x1,y1)c(x2,y2)和d(x2,y1)四个已知像素点的像素值分别为a、b、c和d，待插值点p(uv)的像素值为P。首先，在X方向上进行线性插值，得到式（4.3）和式（4.4。其中，R和Rx2ubu

d,RR(u,yx x

Sx2uau

c,SS(u,yx x （4.5

Py2vRy2

vy1y2

（4.3（4.5Px2uy2

bux1y2vx2uvy1

aux1vy1

通常，在数字图像中a(x1y2、b(x1y1)、c(x2y2和d(x2y1是2*2区域内（i，j（ij+1（i+1，j+1（4.7P(i1u)(j1v)b(ui)(j1v)(i1u)(vj)a(ui)(vj)

的4个邻近已知的像素值通过求和的方式计算得到未知像素点的像素值412w2wswsw48w5w2w

w1w

w4.3（4.9式（4.10）和式（4.11）fi,j表示源图像中点i,j的像素值。fiu,jv

Su

S1u

S2u fi1,jfi1,j

fi,jfi,j

fi1,jfi1,j

fi2,jfi2,j1B fi1,j fi,j fi1,j fi2,jfi1,j fi,j fi1,j fi2,j

Su

S1u

S2u

图像缩放和显示技术中。4.4k幅尺寸大小为rowscols的低分辨率图像是Sk，以及它（rotate[k]创建一幅尺寸大小为10*rows10*cols像素值全为-1的索引图temp10倍出已知像素所在的低分辨率图像序号。如第2幅的低分辨率图像序号为2，其相4.5创建一幅尺寸大小为L*rowsL*cols的高分辨率图像ImgHR，L为高分辨率图像相对于低分辨率图像的放大倍数。在索引图temp中高分辨率ImgHRij所对应像素点为temp10*i/L,10*j/L。在该点附近的5*5区域内寻找是否存在已知的像素点。如果存在有tempi,j1且距离最近，即标记有SkiLjL作为高分辨率图像的像素值。如果没有找到，就采用距离这点最结合亚像素偏移量的双线性插值法为了得到高分辨率图像像素网格中未知率观测图像，并且构成一个正方形把待插值的这个位置像素点包围在正方形内部。图（4.6）4.6结合亚像素偏移量的双线性插值法的示意图k幅尺寸大小为rowscols的低分辨率图像是Sk，以及它（rotate[k]创建一幅尺寸大小为10*rows10*cols像素值全为-1temp创建一幅尺寸大小为L*rowsL*cols的高分辨率图像ImgHR，L为高分辨率图像相对于低分辨率图像的放大倍数。索引图temp中高分辨率ImgHRij所对应的像素点为temp10*iL,10*j/L。在该点附近5*5的区域内寻找是否存在有低分辨率图像像素点。寻找范围由近至远，当存在有tempi,j1且距离最近，即标记已知像素点所在的低分辨率图像的序号。根据索引图中低分辨率图像像素点tempi,j与待插值的高分辨率图像像素点temp10*iL,10*j/LSki/L,j/L所在的同所在的位置包围在其中。即有4种情况：若在待插值点的左上角，则选取Si,j

Si1,j

Si,j

Si1,jkLL

k L

kL

k

Si,j

Si1,j

Si,j

Si1,jkLL

k L

kL

k

；若在Si,j

Si,j

Si,jkLL

kL

kL 插值点的左下角，则选取

Si1,j

Si,j

Si1,jk

kLL

k L ；若在待插值点的右下角，则

Si,j Si1,jkL k 和 （4.7插值算法实验所用的输入图像是由原始的高分辨图像经过降质得到的低分辨率图像对低分辨率图像应用三种经典的插值法和本文结合亚像素偏移行比较得出实验结果。（256*256（b128*128（4.7（4.8是上述插值算法的实验结果。表（41）所示实验数据是插值算法结果与原始高分辨率图像进行比较，计算得到的PSRMSSI。 （1）（2）（3）双三次插值法结（1）（2）（3）双三次插值法结 （4）结合位移量的最近邻插值（5）4.8第二组插值实验采用的是胸部CT图像，原始图像（c）尺寸为512*512，降质后得到的低分辨率图像（d）尺寸为256*256，如图（4.9）所示。图（4.10）展示的是上述插值算法的实验结果。表（4.2）原始高分辨率图像进行比较，计算得到的PSNRMSSIM 4.9第二组实验的原始图像（c）和输入图像4.2（（1）最近邻插值法结 （2）双线性插值法结 （3）双三次插值法结 （4）结合位移量的最近邻插值4.10第三组插值实验采用的是腹部CT图像，原始图像（e）尺寸为512*512，降原始高分辨率图像进行比较，计算得到的PSNRMSSIM。 4.11第三组实验的原始图像（c）和输入图像4.3（1）（2）（3）双三次插值法结 （4）结合位移量的最近邻插值（5）（1）（2）（3）双三次插值法结 （4）结合位移量的最近邻插值（5）4.12量的插值算法要优于经典的插值算法。并且在视觉感受上，对于CT图像来5章迭代反投影方出就备受关注。Irani和Peleg首先将这种方法应用到图像的超分辨率重建上。迭迭代反投影算法相结合运用到图像的超分辨率重建中。Dong等人利用非局部均kkg(y)DHpsff(x)kk

式（5.1）gkyk幅图像，y表示低分f(x中一部分区域的像素点经过核函数为H的点扩散函数（PSF）抽取后得到，D与原始的高分辨率图像足够接近的超分辨率重建图像。式（5.2）定义的就是这f(n1)(x)

fn(x)(gkyk

(y)gn(y))H

kk式（5.2）fn(xngnykk示由经过nfn(xHBP降质模型得到的一组低分辨率模拟图像。插值插值生（反投影降（投影（1）（2）（3）（4）1将输入的KIk(m1m2gk(m1m2第一幅g0(m1m2)作为基准图像。利用相位相关算法进行图像配准，得到gk(m1m2g0(m1m2sk(xk,yk和相对旋转角度k2Ik(m1m2)Ik(m1m2反转k后得到校正过的低分辨率观测图像ik(m1m2得到。对校正旋转后的第一幅低分辨率观测图像i0(m1m2XYLF0(n1n23hk(m1m2与校正旋转后的低分辨率观测图像ik(m1m2sk(xk,yk)和相对旋转角度k和扩大倍数L，对当前的高分辨率估计图像Fn(n1,n2)进行降质得到K幅分别对应输入hk(m1m2sk(xk,yk)数，令xkxk分别xk的整数部分和小数部分，令ykyk分别为yk的p=F(n+[x],n+[y])(1-x`)(1-y` =F(n+[x], p=F(n+[x], =F(n+[x],n+[y])x` hk(m1,m2)=p1+p2+p3+

迭代反投影的过程就是将校正后的低分辨率观测图像ik(m1m2和低分辨率模拟图像hk(m1m2)的差值数据ek(m1m2)反投影到当前的高分辨率估计图像Fn(n1n2Fn1(n1n2ek(m1,m2)ik(m1,m2)hk(m1,m2

反投影过程具体描述如下：第一步，将误差数据ek(m1m2进行双线性差值，在XYLEk(n1n2)。第二ik(m1m2)Ek(n1n2)Ek(n1n2)sk(xk,ykxkykEk(n1n2的每个数据会依据距离加Fn(n1n24个像素点上。式（5.5）所示的 (n[x],n[y])f

[x],n[y])1E(n,n)(1x`)(1y`n1

1 (n[x]1,

[y])f(n[x]1,

[y])

E(n,n)(1y`n1

1

(n[x],n[y]1)f(n[x],

[y]1)

E(n,n)(1x`)n1

`fn1(n1[xk]1,n2[yk]1)fn(n1[xk]1,n2[yk]1)

Ek(n1,n2)xk

4Fn(n1n2降质得到的低分辨hk(m1m2与实际的校正旋转后的低分辨率观测图像ik(m1m2的误差error，判断更新后的高分辨率估计图像Fn1(n1n2error若小hk(m1m2与实际的低分辨率观测定义的是误差值error的计算公式。 meanmm ik(m1,m2)gk(m1,m2error

2m20m1

(i(m,m)g(m,m)4与原始图像相比较，分析实验结果。具体的实验步骤如下：第一步，将输入的4幅低分辨率观测图像Ik(m1m2)转换成灰度图像gk(m1m2g0(m1m23幅低分辨率观测图像相对这幅基准图像的亚像素平移量和旋测图像进行旋转校正得到ik(m1m2F0(n1n2降质后得到的低分辨率模拟图像hk(m1m2)。计算得到低分辨率模拟图像hk(m1m2与校正旋转后的低分辨率观测图像ik(m1m2之间的误差ek(m1m2ek(m1m2Ek(n1n2Ek(n1n2Fn(n1n2Fn1(n1n240。满足条件，迭