一种基于CT图像边缘的有理插值缩放方法与流程

本发明涉及图像处理技术领域，尤其涉及基于ct图像边缘特征的有理插值缩放方法。

背景技术：

图像缩放是图像处理的基本问题，并在实际中得到了广泛应用，如在医学、航天以及视频应用等系统中便常常用到图像放大，其主要目的是使得到的新图像更好地反映原场景的形状、尤其是细节信息。特别地，在医疗诊断中，医生对临床疾病诊断重要辅助诊断手段-ct图像的观察，往往需要放大操作，以便适用于细节，尤其是病灶的观察，从而提高医生对疾病诊断的准确性。

图像放大的本质是根据已知的低分辨率原始图像像素点的灰度值产生未知高分辨率像素点的灰度值。事实上，图像数据可看成是从一个原表面上采样得到的，因此图像放大问题可变为由图像数据反向重建原表面，并进行重采样的问题。目前有很多构造图像曲面的方法，其中多项式函数由于其形式简单、计算容易等优点，成为图像放大中最常用的方法之一，如早期的双线性插值(bilinearinterpolatio)、双三次样条插值(bicubicsplineinterpolation)，以及之后改进的连续的数学模型如基于bézier曲面的插值、基于b样条曲面的插值等。这些方法最大的优势是它们相对简单，但同时也容易产生比较明显的马赛克，模糊以及锯齿状现象。尤其是对边缘明显的图像，一般得不到理想的结果。

图像边缘信息是影响视觉效果的关键因素，也是图像匹配、识别等图像处理中的关键因素。为使插值后图像的边缘保持良好的清晰度，一些研究人员提出了基于边缘的插值方法。这些方法往往以原始图像的边缘作为条件进行插值用以放大图像，如可根据相关性最小准则确定物体的边界，然后根据边界像素调整其周围像素点所占插值的权重。这些方法可以改进图像的主观质量，改善视觉效果。但边缘算法常常会将一些假边缘判定为边缘，从而造成错误，视觉效果比较差。

在医学辅助诊断中，医生除了应用其专业知识进行目测等常规操作外，还需要对对图像进行缩放、旋转等操作，从不同角度观察图像，甚至精确估计病灶的位置和大小。除了医学硬件设备和辐射剂量限制等硬缺陷，ct图像病变组织和正常组织的边缘本身还存在比较模糊的特点，如果采用上述各方法进行插值放大，放大后的图像很可能无法正确显示病变组织和正常组织之间的关系。

技术实现要素：

本发明的目的就是为了解决上述问题，提供基于边缘信息的ct图像的有理插值缩放方法，通过边缘检测算法获得边界曲线以及关键像素，放大时获得边界像素的拟合曲线，重采样获得新边界，以保持视觉上边界的锐化度。

为了实现上述目的，本发明采用如下技术方案：

一种基于ct图像边缘的有理插值缩放方法，包括

步骤一，通过边缘检测算法获取图像单像素边界；

步骤二，获得边界像素的拟合曲线，重采样获得放大后的新边界；

步骤三，根据不同情况对非边界像素进行重采样。

所述步骤一中，利用自适应阈值的canny算子获得单像素边界。

所述步骤二中，获得边界像素拟合曲线的方法包括，求解各边界像素点在u向和v向上的导数，选择相邻两边界像素pi,j，ps,t进行三次hermite插值，位置插值曲线记为l(t)，像素插值曲线记为p(t)，其中t为插值点的参数，分别用来确定像素位置及像素值。

三次hermite插值形式为：

l(t)＝f0(t)li，j+f1(t)ls，t+g0(t)l′i，j+g1(t)l′s，tt∈[0，1](1)

p(t)＝f0(t)pi，j+f1(t)ps，t+g0(t)p′i，j+g1(t)p′s，tt∈[0，1](2)

其中：

f0(t)＝2t³-3t²+1,f1(t)＝-2t³+3t²

g0(t)＝t³-2t²+t,g1(t)＝t³-t²；li,j为ij点的位置，ls,t为st点的位置，l’i,j为ij点位置导数，l’s,t为st点位置导数，pi,j为ij点的像素值，ps,t为st点的像素值，p’i,j为ij点的像素导数，p’s,t为st点的像素导数。

将参数t＝1/s，t＝2/s……t＝s-1/s代入公式(1)(2)，放大s倍再四舍五入得到放大后的曲线所经过的像素位置及像素值：

各边界像素点在u向和v向上的导数由保型插值方法估计得到，求解方法为：

用相邻象素点pi-1,j，pi,j和pi+1,j三点估计pi,j处的u向导数(p′i,j)u，用象素pi,j-1，pi,j和pi,j+1三点估计pi,j处的v向导数(p′i,j)v，对(p′i,j)u，一阶导数计算方法为：(p′i,j)u＝(pi+1,j-pi-1,j)/2。

所述步骤三中，对于放大后ct图像邻近新边界像素的关键像素采用分片双变量有理插值算法进行重采样；对于处于平滑区域的普通像素采用双线性插值算法进行重采样。

关键像素点重采样的方法为，

pi,j(x,y)是由九个点决定的[xi,xi+1；yj,yj+1]的双变量有理插值函数，其x-方向上的插值曲线定义为：

其中，α为插值点在x方向上的参数，为简化公式，记：hi＝xi+1-xi,lj＝yj+1-yj,θ＝(x-xi)/hi,η＝(y-yj)/lj，

则，

则pi,j(x,y)为：

其中β为插值点在y方向上的参数，ωr,s(θ,αi；η,βj)＝ωr(θ,αi)ωs(η,βj)，

对这些像素进行放缩，就是根据放缩比例对原始图像进行重采样；放大s倍，则原图像i(x,y)在x，y方向按照1/s的间隔重采样，得到放大s倍后的图像i′(x′,y′)。

对于处于平滑区域的普通像素的求解步骤具体为：对网格内的点p的插值计算公式为：

p＝(1-dx)pi,j+(1-dx)(1-dy)pi,j+1+dx(1-dx)pi+1,j+1+dxdypi+1,j，其中dx为p点x值与pi点x值得距离，dy为p点y值与pi点y值的距离。

本发明的有益效果：

人们的视觉感知特性是非线性的，因此用非线性函数中的有理函数建立数学模型会得到更理想的效果。有理函数插值尽管比多项式复杂，但用它近似表示函数时，却比多项式灵活、更能反映函数的一些特性，尤其是一些突变特性。

通过边缘检测算法获得边界曲线以及关键像素，放大时获得边界像素的拟合曲线，重采样获得新边界，以保持视觉上边界的锐化度；非边界像素则根据不同情况分为有理插值和普通双线性插值两种。由于对图像边缘作了特殊处理，从而避免在图像缩放过程中丢失边缘信息，从而更好地展现图像的细节，能较好的显示出病变的影像信息。放大后的ct图像能清晰的显示病灶与周围正常组织器官的关系，及早发现一些不容易发现的病变，为临床医师做出正确的临床治疗决策提供重要的参考，并对手术可行性难度及风险评价也有重要指导意义。

附图说明

图1是原始图像的边界轮廓示意图，黑色连线代表求解的边界轮廓。

图2是放大2倍后形成的新的边界轮廓示意图，黑色连线代表求解的新边界像素，其中白色圆点是原始数据点放大后的位置，黑色圆点则是放大2倍后新增加密度的采样点。

图3是非边界像素的两种情况。三角形覆盖的像素点代表邻近新边界像素的关键像素，菱形覆盖的像素点代表平滑区域的普通像素点。

图4是平滑区域像素点的双线性插值示意图，其四个角点均不是边界点，且dx，dy分别为待插值点距离区域上边界及做边界的距离。

图5(a)是采用样条插值方法对ct图像放大2倍后得到的，图5(b)是由本发明的方法对ct图像放大2倍后得到的。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

针对利用现有方法进行图像放大时边界容易模糊、失真等现象，本发明讨论怎样避免在图像缩放过程中丢失人体视觉的重要因素—图像边缘信息，保证边界轮廓在放大后的锐化程度，保持清晰的图像细节，获得较好的主观视觉效果。对离散图像做放大操作，最有效的方法之一还是构造图像的原始曲面，并进行高密度重采样。本发明提出一种具有良好视觉效果的保持边缘特征的ct图像缩放方法。本发明的基本思想可简述如下：一般一幅图像由突变区域和平滑区域组成。平滑区域的像素值都非常的接近，因此，用一般的插值算法，效果都已足够好，突变区，一般为边界，其像素值变化很大，如果单纯用常规的放大方法会出现边缘模糊、图像细节丢失的现象。好的插值方法应当对图像的边缘特殊考虑，才能避免或减轻边缘模糊和锯齿，很好地保持边缘的光滑性和锐利性。具体可通过以下步骤完成放缩操作：通过边缘检测算法获取ct图像的初始单像素边界，并获得边界曲线以及关键像素，放大时获得边界像素的拟合曲线，重采样获得新边界，以保持视觉上边界的锐化度；非边界像素则根据不同情况分为有理插值和普通双线性插值两种。

本发明可分为三大部分：a、获取图像单像素边界；b、边界像素重采样；c、非边界像素重采样。

则利用一种基于ct图像边缘的有理插值缩放方法进行ct图像缩放的步骤描述如下：

step1：求取图像单像素边界；

step2：求边界数据点在u向和v向上的导数；

step3：进行边界数据点三次hermite插值，并按照缩放因子进行重采样；

step4：非边界像素中关键像素的有理插值重采样。

step5：非边界像素中其他像素的双线性插值重采样。

a、获取图像单像素边界

首先获取单像素边界。通过边缘检测算法获得边界曲线以及关键像素，为了使最终的图像具有清晰的边缘，采用一种有效的边界检测算子提取图像的初始单像素边界。canny算子具有很好的边缘检测性能，它提取的边缘为单像素边缘，位置非常准确，而且边缘完整、连续性很好，因此，选择自适应阈值的canny算子完成该步骤。

b、边界像素重采样

获得边界像素的拟合曲线，重采样获得放大后的新边界，以保持视觉上边界的锐化度。

(1)求解各边界像素点在u向和v向上的导数。

分别基于位置及像素值对获得的边界像素进行曲线拟合。为了后续计算方便及获得较为精确的值，选择相邻两边界像素pi，j，ps，t进行三次hermite插值。位置插值曲线记为l(t)，像素插值曲线记为p(t)，分别用来确定像素位置及像素值。其三次hermite插值形式如下：

l(t)＝f0(t)li，j+f1(t)ls，t+(g0(t)li，j+g1(t)ls，tt∈[0，1](1)

p(t)＝f0(t)pi，j+f1(t)ps，t+g0(t)p′i，j+g1(t)p′s，tt∈[0，1](2)

其中

f0(t)＝2t³-3t²+1，f1(t)＝-2t³+3t²

g0(t)＝t³-2t²+t，g1(t)＝t³-t²(3)

l(t)则可以分解为x(t)及y(t)。公式(1)(2)中的各像素的位置偏导及像素偏导值均未知，需要估计。可选择bessel方法(抛物线方法)来估计上述各偏导值。由于求解方法类似，仅以像素pi，j的一阶偏导数(p′i，j)x为例介绍具体估计方法。

由pi-1，j，pi，j和pi+1，j三点构造二次多项式插值曲线，并对曲线求导，所求得的一阶偏导作为(p′i，j)x的近似值，即(p′i，j)x由下式定义

其中

由于ct图像一般是规则数据场，pi-1,j，pi,j和pi+1,j三点间为等步长，所以

对边界数据点p1,j和pn,j，由

得：(p′1,j)x＝2δp1,j-(p′2,j)x，同理可得：(p′n,j)x＝2δpn-1,j-(p′n-1,j)x

同理，获得y向一阶偏导(p′i,j)y，采用pi,j-1，pi,j和pi,j+1三点进行估计，则

同理可得p′s,t。对于x(t)、y(t)中的偏导则分别采用pa,b(边界像素，pi,j的前一个点)，pi,jps,t，或pi,j，ps,t，pc,d(边界像素，ps,t的后一个点)的位置x值和y值进行估计。将所求偏导代入式(1)(2)，获得位置插值曲线l(t)，像素插值曲线p(t)。由于数据点偏导值在两像素处是相同的，相邻两段曲线是c¹连续的。

(2)边界重采样

根据l(t)求解放大后的曲线，如图1中的黑色连线代表求解的边界像素。将参数t＝1/s，t＝2/s……t＝s-1/s代入公式(1)，放大s倍再四舍五入：

即为放大后的曲线所经过的像素位置，如图2中的黑色连线代表求解的新边界像素，而其中白色圆点是原始数据点放大后的位置，黑色圆点则是放大s倍后新扩充的点。

c、非边界像素重采样

上面描述了边缘的提取以及边缘像素的插值放大，处于非边界轮廓区域的像素采用如下方法进行求解。非边界像素重采样则根据不同情况分为有理插值和普通双线性插值两种。根据上述分析，边界像素附近的像素变化比较剧烈，而有理函数插值作为典型的非线性逼近方法之一，用它近似表示函数时，更能反映函数的一些特性，尤其是一些突变特性。因此，更好地反映ct图像相邻像素之间的数值关系。对于图像中平滑区域，由于各像素值会非常接近，用一般的插值算法，如双线性或双三次插值，视觉效果不会很差。因此本发明分两种情况对非边界轮廓像素进行求解，如图3，一种为邻近新边界像素的关键像素，一种为处于平滑区域的普通像素。

(1)邻近新边界像素的关键像素的求解。

对于图3中的关键像素，为了保留这些数据的突变性，本发明采用分片双变量有理插值，每片都含有参数，非常便于局部修改。具体形式如下：

pi,j(x,y)是由九个点{(xr,ys,pr,s,r＝i,i+1,i+2,s＝j,j+1,j+2)}决定的[xi,xi+1；yj,yj+1]的双变量有理插值函数，满足如下条件：

pi,j(xr,ys)＝p(xr,ys),r＝i,i+1,s＝j,j+1

则其x-方向上的插值曲线定义如下：

其中：hi＝xi+1-xi,lj＝yj+1-yj,θ＝(x-xi)/hi,η＝(y-yj)/lj.

则pi,j(x,y)定义如下：

其中ωr,s(θ,αi；η,βj)＝ωr(θ,αi)ωs(η,βj)。

对这些像素进行放缩，就是根据放缩比例对原始图像进行重采样。放大s倍，则原图像i(x,y)在x，y方向按照1/s的间隔重采样，得到放大s倍后的图像i′(x′,y′)。

(2)普通像素的求解。

如图4，网格的四个顶点均不是边缘点，为计算简单，对网格内的点p的插值可采用普通的双线性插值算法，计算公式为：

p＝(1-dx)pi,j+(1-dx)(1-dy)pi,j+1+dx(1-dx)pi+1,j+1+dxdypi+1,j(7)

我们将本发明应用在ct图像处理系统中。对ct图像进行放大，能为医生诊断带来很大的方便。其实现可通过软件编程实现。首先给出本发明编程的伪代码。

在ct图像处理系统中对ct图像进行放大，将需要放大的图像image以及需要放大的倍数s传递给zoombasedonriep方法。

作为实施例，选择医生诊断过程中需要用到的ct图像，分别对放大2倍和4倍。效果图见附图5(a)-图5(b)。从图示可知本发明可得到比较好的效果。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。