CN105389476A - 基于梯度特征的调强放射治疗计划剂量数据的插值算法 - Google Patents

Abstract

本发明涉及一种基于梯度特征的调强放射治疗计划剂量数据的插值算法，属于调强放射治疗技术领域。该插值算法根据调强放射治疗计划剂量数据平面的梯度特征通过采用改进传统的Canny算法获取调强放射治疗计划剂量数据平面上各个梯度边缘点和非梯度边缘点；并根据所获得的梯度边缘点所对应梯度剖面的锐度以及非梯度边缘点所对应的偏离系数得到调强放射治疗计划剂量数据平面上每一点对应的双三次插值核的系数；并对每一个待插值点使用所述双三次插值核的系数进行双三次插值，从而求得每一个待插值点的调强放射治疗计划剂量数据。本发明的插值算法在减小误差的同时，克服了传统双线性插值的平滑效应，保留了原有的梯度信息。

Description

基于梯度特征的调强放射治疗计划剂量数据的插值算法

技术领域

本发明涉及一种调强放射治疗计划剂量数据的处理方法，特别是涉及一种基于梯度特征的调强放射治疗计划剂量数据的插值算法，属于调强放射治疗技术领域。。

背景技术

调强放射治疗(IntensityModulatedRadiationTherapy,IMRT)剂量验证作为放射治疗计划可靠性的重要保证，剂量验证的准确性对放射治疗效果以及患者的安全有着重要意义。剂量验证的基本原理是对于所制定的调强放射治疗计划，利用DD(DoseDifference)、DTA(DistanceToAgreement)、Gamma因子等剂量验证方法来比较治疗计划系统(TreatmentPlanSystem，TPS)对所制定的治疗计划进行计算产生的调强放射治疗计划剂量数据和按治疗计划进行实施时测量硬件得到的测量剂量数据；调强放射治疗计划剂量数据从产生到进行剂量验证的过程中有时需要经过多次插值计算。一方面，TPS使用相关模型对调强放射治疗计划进行计算得到一个较为精确的调强放射治疗计划剂量数据平面，而TPS内部有时还需要改变这个较为精确的调强放射治疗计划剂量数据平面上数据点之间的物理间距，以满足特定的输出要求，因而需要对这个较为精确的调强放射治疗计划剂量数据平面进行插值计算；另一方面，剂量验证时需要先将调强放射治疗计划剂量数据平面和测量剂量数据平面按等剂量中心进行重合之后再比较相同位置上的点，但由于其两者的物理间距不相同，因此会导致所述调强放射治疗计划剂量数据平面和测量剂量数据平面两者重合之后其剂量数据平面上的点将会错位分布，此时仍需要进行插值计算，从而使在调强放射治疗计划剂量数据平面上得到与测量剂量数据平面相同位置的点的剂量值，进而方可以使用DD、DTA、Gamma因子等剂量验证方法进行剂量验证。

在实际应用中主要通过如DoseLab，MapCkeck等软件进行剂量验证，它们对调强放射治疗计划剂量数据进行插值时主要采用的方法多为双线性插值算法，该方法是利用待插值点的四个相邻点的值作线性内插，其公式表示如下(贾永红，《数字图像处理》武汉大学出版社，2003.)：

f(i+u,j+v)＝(1-u)(1-v)f(i,j)+(1-u)vf(i,j+1)

A(1)

+u(1-v)f(i+1,j)+uvf(i+1,j+1)

其中，u和v分别为待插值点在水平方向、垂直方向上与点(i,j)的距离；但是由于调强放射治疗为了达到较好的治疗效果同时保护正常组织，一般在病灶区和正常区域之间形成较大的剂量梯度，因此调强放射治疗计划剂量数据中往往存在较多的梯度边缘点，而传统的双线性插值由于其算法本身具有平滑效应因而无法正确的处理这些梯度边缘，从而会产生较大的计算误差，并且使整个调强放射治疗计划剂量数据分布趋于平滑而失去原有的梯度信息，进而对治疗医生和物理师的判断产生误导。

发明内容

本发明的目的就是在于克服上述所述传统双线性插值算法存在的缺陷和不足，提出一种基于梯度特征的调强放射治疗计划剂量数据的插值算法。该插值算法能够在减小计算误差的同时，克服传统双线性插值的平滑效应，且保留了原有的梯度信息；并有利于治疗医生和物理师对调强放射治疗计划剂量数据的正确判断，从而避免产生误导。

为实现上述目的，本发明采用由以下技术措施构成的技术方案来实现的。

本发明所述基于梯度特征的调强放射治疗计划剂量数据的插值算法，是根据调强放射治疗计划剂量数据平面的梯度信息，使用改进传统的Canny边缘探测算法获取该剂量数据平面上各梯度边缘点和非梯度边缘点，根据获取的梯度边缘点所对应梯度剖面的锐度以及非梯度边缘点所对应的偏离系数，获得调强放射治疗计划剂量数据平面上每一点为双三次插值中心点时对应的双三次插值核的系数，对每一个待插值点使用该双三次插值核的系数进行双三次插值，从而得到每一个待插值点的调强放射治疗计划剂量数据；包括以下步骤：

步骤1：对治疗计划系统产生的调强放射治疗计划剂量数据平面I_L，计算该调强放射治疗计划剂量数据平面上每一点(i,j)的梯度▽I_L(i,j)以及该点对应的梯度模||▽I_L(i,j)||；

步骤2：对传统Canny边缘探测算法进行改进，根据步骤1中得到的调强放射治疗计划剂量数据平面上每一点的梯度信息，使用改进的Canny边缘探测算法获取调强放射治疗计划剂量数据平面上各个梯度边缘点(i′,j′)和非梯度边缘点(i″,j″)；

步骤3：对步骤2所获取的调强放射治疗计划剂量数据平面上各个梯度边缘点(i′,j′)进行梯度剖面的追踪，从而得到每一个梯度边缘点对应的梯度剖面；

步骤4：对步骤3所得到的调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)所对应梯度剖面的锐度σ(P(i′,j′))进行计算；对步骤2所述调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)自定义的偏离系数ρ(i″,j″)进行计算；

步骤5：根据步骤4得到调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)对应梯度剖面的锐度σ(P(i′,j′))，确定锐度σ(P(i′,j′))与以梯度边缘点(i′,j′)为双三次插值中心点时所对应双三次插值核的系数a(i′,j′)之间的函数关系；

步骤6：根据步骤4得到调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)的偏离系数ρ(i″,j″)，确定偏离系数ρ(i″,j″)与以非梯度边缘点(i″,j″)为双三次插值中心点时所对应的双三次插值核的系数a(i″,j″)之间的函数关系；

步骤7：将步骤5和步骤6所得到的双三次插值核的系数的两函数结果合并，从而得出以调强放射治疗计划剂量数据平面上每一点(i,j)为双三次插值中心点时对应的双三次插值核的系数a(i,j)；

步骤8：对于待插值点，设置其坐标为(a,b)，若a∈[i,i+1)且b∈[j,j+1)，根据步骤7得出的所述双三次插值核的系数a(i,j)，令点(i,j)为双三次插值中心点，并以该双三次插值核的系数a(i,j)对待插值点进行双三次插值计算，从而得出每一个待插值点的调强放射治疗计划剂量数据。

上述技术方案中，所述步骤2中对传统Canny算法进行改进，其采用如下操作：

(a)略去传统Canny算法中高斯滤波的过程，获得略去高斯滤波后的结果；

(b)对略去传统Canny算法中高斯滤波的结果进行判断，假设(i′,j′)点为传统Canny算法略去高斯滤波后结果中的调强放射治疗计划剂量数据平面上的一个梯度边缘点，且以(i′,j′)点为中心的3×3邻域中含有其他梯度边缘点，若(i′,j′)点梯度模小于所述邻域中其他梯度边缘点梯度模，则重新将(i′,j′)点标记为非梯度边缘点。

上述技术方案中，所述步骤5中调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)对应梯度剖面的锐度σ(P(i′,j′))与以点(i′,j′)为双三次插值中心点时对应的双三次插值核的系数a(i′,j′)之间的函数关系，由以下公式(2)表示：

$a(i^{\′},j^{\′})=-0.5\frac{1}{1+\ln ({\mathrm{\σ}}_{\max }/\mathrm{\σ}(i^{\′},j^{\′}\left)\right)}\exp \left\{{(1-\frac{\mathrm{\σ}(i^{\′},j^{\′})}{{\mathrm{\σ}}_{\max }})}^2\right\}---\left(2\right)$

其中σ_max表示该调强放射治疗计划剂量数据平面上锐度的最大值。

上述技术方案中，所述步骤5中调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)的偏离系数ρ(i″,j″)与以点(i″,j″)为双三次插值中心点时对应的双三次插值核的系数a(i″,j″)之间的函数关系，由以下公式(3)表示：

$a(i^{\′\′},j^{\′\′})=-0.5\exp \{-{\left(\frac{\mathrm{\ρ}(i^{\′\′},j^{\′\′})-{\mathrm{\ρ}}_{\min }}{{\mathrm{\ρ}}_{\max }-{\mathrm{\ρ}}_{\min }}\right)}^2\}---\left(3\right)$

其中ρ_max为调强放射治疗计划剂量数据平面上最大偏离系数，ρ_min为调强放射治疗计划剂量数据平面上最小偏离系数。

本发明所述的插值算法与现有技术相比较所具有的有益技术效果：

1、本发明所述基于梯度特征的调强放射治疗计划剂量数据的插值算法，采用对Canny边缘探测算法进行了改进，使用改进的Canny边缘探测算法以获取调强放射治疗计划剂量数据平面上每一点的梯度边缘点和非梯度边缘点，以此减小假梯度边缘点的产生；通过对调强放射治疗计划剂量数据平面上非梯度边缘点使用自定义偏离系数，可以从一定程度反映出减小假梯度边缘过程中造成缺失的梯度边缘点所具有的剂量数据分布特点，从而使这种改进达到较好的效果。

2、本发明所述基于梯度特征的调强放射治疗计划剂量数据的插值算法，经实施例所得结果可以看出，这种对传统Canny边缘探测算法的改进提高了计算精确度，由于其减小了梯度边缘点数，从而可以减小计算的复杂度。

3、本发明所述基于梯度特征的调强放射治疗计划剂量数据的插值算法，通过该插值算法所得的结果进一步得知，与传统的双线性插值算法相比较，在减小插值计算误差的同时，不仅保持了整个调强放射治疗计划剂量数据平面的梯度特征，而且克服了双线性插值的平滑效应。

4、本发明所述基于梯度特征的调强放射治疗计划剂量数据的插值算法，由于减小了计算误差，克服了传统双线性插值的平滑效应，保留了原有的梯度信息；从而避免治疗医生和物理师在对调强放射治疗计划剂量数据进行判断时产生误导。

附图说明

图1为剂量验证时测量剂量数据平面和调强放射治疗计划剂量数据平面按等剂量中心对齐示意图；

图2为本发明所述插值算法的流程示意框图；

图3为传统梯度剖面追踪示意图；

图4为不同传统双三次插值核的系数下的双三次插值核的变化趋势以及与双线性插值对比曲线；

图5为传统梯度剖面上的插值以及其梯度模的分布；其中，图(a)展示了一个梯度剖面上使用双线性插值所造成的误差，其中X₀为调强放射治疗计划剂量数据平面上梯度边缘点，X₁和X₂分别为两个已知数值的点，P₁和P₂分别为两个梯度剖面所对应的剂量分布，图(b)展示了图(a)中的两个梯度剖面P1和梯度剖面P2的梯度模的分布曲线；

图6为本发明实施例所述在不同调强放射治疗计划剂量数据物理间距下，基于梯度特征的调强放射治疗计划剂量数据插值算法、双线性插值算法以及传统双三次插值算法当插值物理间距为0.4mm时的误差对比；

图7为本发明实施例所述在不同调强放射治疗计划剂量数据物理间距下，基于梯度特征的调强放射治疗计划剂量数据插值算法、双线性插值算法以及传统双三次插值算法当插值物理间距为0.48mm时的误差对比；

图8为本发明实施例所述在不同调强放射治疗计划剂量数据物理间距下，基于梯度特征的调强放射治疗计划剂量数据插值算法、双线性插值算法以及传统双三次插值算法当插值物理间距为0.5mm时的误差对比；

图9为本发明实施例所述在不同调强放射治疗计划剂量数据物理间距下，基于梯度特征的调强放射治疗计划剂量数据插值算法、双线性插值算法以及传统双三次插值算法当插值物理间距为0.6mm时的误差对比；

图10为本发明实施例所述在不同调强放射治疗计划剂量数据物理间距下，基于梯度特征的调强放射治疗计划剂量数据插值算法、双线性插值算法以及传统双三次插值算法在几种插值物理间距下平均误差对比；

图11为本发明实施例所述在不同调强放射治疗计划剂量数据物理间距下，基于梯度特征的调强放射治疗计划剂量数据插值算法与双线性插值算法结果的平均梯度模的对比。

图中，1为调强放射治疗计划剂量数据点，2为测量剂量数据点，3为等剂量中心点，4为测量剂量数据物理间距，5为调强放射治疗计划剂量数据物理间距，6为梯度剖面与调强放射治疗计划剂量数据坐标网格相交点，该点为已存在的数据点，7为追踪形成的梯度剖面，8为调强放射治疗计划剂量数据平面上梯度边缘点，9为梯度剖面路径上与调强放射治疗计划剂量数据坐标网格相交点，且该点为原来不存在的数据点，10为一维双线性插值权重随距离变化趋势，11为双三次插值核的系数为-0.5时双三次插值核的变化趋势，12为双三次插值核的系数为-1时双三次插值核的变化趋势，13为双三次插值核的系数为-1.5时双三次插值核的变化趋势，14为调强放射治疗计划剂量数据平面上梯度边缘点对应梯度剖面P1的真实剂量数据分布，15为调强放射治疗计划剂量数据平面上梯度边缘点对应梯度剖面P2的真实剂量数据分布，16为双线性插值所拟合出的剂量数据分布，17为待插值点处双线性插值的计算结果，18为梯度剖面P1在待插值点处的真实剂量数据，19为梯度剖面P2在待插值点处的真实剂量数据。

具体实施方式

下面结合附图并用具体实施例对本发明作进一步详细描述，所述的实施例只是对本发明的一个具体说明，而不应理解为是对本发明保护内容的任何限制。

图1展示了剂量验证过程中测量剂量数据平面和调强放射治疗计划剂量数据平面按等剂量中心对齐示意图，它们两者的等剂量中心点3的相对位置以及测量剂量数据物理间距4和调强放射治疗计划剂量数据物理间距5不相同，因此会导致两者重合之后其两个平面上的点即调强放射治疗计划剂量数据点1和测量剂量数据点2的错位分布；此时需要进行插值计算从而在调强放射治疗计划剂量数据平面上得到与测量剂量数据平面相同位置的点的剂量数据值。图1中所述的测量剂量数据物理间距4是指测量剂量数据平面在水平方向上两个相邻测量剂量数据点之间的物理距离，或测量剂量数据平面在垂直方向上两个相邻测量剂量数据点之间的物理距离。所述的调强放射治疗计划剂量数据物理间距5是指调强放射治疗计划剂量数据平面在水平方向上两个相邻计划剂量数据点之间的物理距离，或调强放射治疗计划剂量数据平面在垂直方向上两个相邻计划剂量数据点之间的物理距离。

实施例

图2展示了本发明所述基于梯度特征的调强放射治疗计划剂量数据的插值算法的流程图，本实施例根据此流程图进行详细描述，其具体步骤如下：

(1)调强放射治疗计划剂量平面上梯度以及梯度模的计算

读取由调强放射治疗计划系统产生的调强放射治疗计划剂量数据平面I_L，对调强放射治疗计划剂量数据平面上每一点(i,j)，由传统的梯度表达式可知其水平方向的梯度▽I_Lx和垂直方向的梯度▽I_Ly，由下面公式(4)表示(张广军，《机器视觉》科学出版社：2005.)：

▽I_Lx＝I_L(i+1，j)-I_L(i，j)

(4)

▽I_Ly＝I_L(i，j)-I_L(i，j+1)

本实施例中采用调强放射治疗计划剂量数据平面上每一点(i,j)的前后两个点并结合公式(4)得到点(i,j)的梯度(GradientProfilePriorandItsApplicationsinImageSuper-ResolutionandEnhancement[J]，IEEETransactionsOnImageProces-sing，vol.20，num.6.pp.1529-1542，2011)，由下面公式(5)表示：

$\▿I_L(i,j)=\frac{I_L(i+1,j)-I_L(i-1,j)}{2}\stackrel{\→}{e_i}+\frac{I_L(i,j+1)-I_L(i,j-1)}{2}\stackrel{\→}{e_j}---\left(5\right)$

调强放射治疗计划剂量数据平面上每一点(i,j)的梯度模||▽I_L(i,j)||，由下面公式(6)表示：

$\left|\right|\▿I_L(i,j)\left|\right|=\sqrt{{(\▿I_L(i,j)\·\stackrel{\→}{e_i})}^2+{(\▿I_L(i,j)\·\stackrel{\→}{e_j})}^2}---\left(6\right)$

公式(5)和公式(6)中，为水平方向上的单位矢量，为垂直方向的单位矢量。调强放射治疗计划剂量数据平面上的每一点(i,j)表示以调强放射治疗计划剂量数据平面的左下角(0,0)点为坐标原点，在水平方向上与坐标原点(0,0)之间有i个调强放射治疗计划剂量数据物理间距，在垂直方向上与坐标原点(0,0)之间有j个调强放射治疗计划剂量数据物理间距；假设水平方向和垂直方向上各自的总调强放射治疗计划剂量数据点数分别为i_max和j_max，那么i的取值范围为0到i_max-1，j的取值范围为0到j_max-1。

(2)调强放射治疗计划剂量数据平面上梯度边缘点和非梯度边缘点的获取

对调强放射治疗计划剂量数据平面进行梯度边缘点探测，首先根据传统Canny边缘探测算法(Acomputationalapproachtoedgedetection[J]，IEEETransactionsonPatternAnalysisAndMachineIntelligence，vol.8，num.6.pp.679-698，1986)，对调强放射治疗计划剂量数据平面上每一点(i,j)取以该点为中心的3×3邻域，并作如下处理：

①对该3×3邻域进行高斯滤波以消除噪声所形成的假梯度边缘点；

②沿着所述点(i,j)的梯度方向作直线与3×3邻域的边界交于两点，若点(i,j)的梯度模||▽I_L(i,j)||同时大于两个交点处所对应的梯度模，则认为点(i,j)为调强放射治疗计划剂量数据平面上的梯度边缘点；

所述两个交点处所对应的梯度模可以使用线性插值得出；但是在本实施例中，一方面，计算数据是不存在噪声的；另一方面，对于传统Canny算法得出的某些梯度边缘点，虽然满足传统Canny算法的相关条件，但是其梯度模较小而无法作为梯度边缘点，为了避免这种情况，在一个3×3邻域中若存在多个梯度边缘点，则只取其中梯度模最大点作为所述3×3区域内的唯一梯度边缘点，这种评判标准虽然较为严格，但是本算法中对于非梯度边缘点采用了自定义的偏离系数，因此可以从一定程度反映出减小假梯度边缘过程中造成缺失的梯度边缘点所具有的剂量数据分布特点，从而达到较好的计算结果。因此本插值算法首先对传统Canny算法作如下改进：

(a)略去传统Canny算法中高斯滤波的过程，获得略去高斯滤波后的结果；

(b)对略去传统Canny算法中高斯滤波的结果进行判断，假设点(i′,j′)为传统Canny算法略去高斯滤波后结果中的调强放射治疗计划剂量数据平面上的一个梯度边缘点，且以点(i′,j′)为中心的3×3邻域中含有其他梯度边缘点，若点(i′,j′)梯度模小于所述邻域中其他梯度边缘点的梯度模，则重新将(i′,j′)标记为非梯度边缘点。

进行上述操作之后，调强放射治疗计划剂量数据平面上非梯度边缘点由以下两部分组成：

(i)传统Canny算法中不是梯度边缘点的点；

(ii)改进Canny算法后所排除掉的梯度边缘点。

事实上，对传统Canny梯度检测算法进行改进之后，一方面，计算结果的误差减小；另一方面，由于本插值算法对于梯度边缘点的处理较为复杂和费时，因此对传统Canny算法进行改进之后减小了梯度边缘点数，从而后续处理中会减小计算时间。

(3)调强放射治疗计划剂量数据平面上梯度边缘点上梯度剖面的追踪

对于所获取的调强放射治疗计划剂量数据平面上梯度边缘点8，如图3所示，从该梯度边缘点8开始沿着两个方向对其梯度剖面进行追踪：先从该梯度边缘点8开始，沿着该梯度边缘点8梯度的正方向作直线与坐标网格相交于一点，再从这个交点开始沿着该交点处梯度的正方向作直线得到下一个交点，重复以上过程直到梯度模不再减小，从而得到沿着梯度正方向的梯度剖面；再从该梯度边缘点8开始，沿着该梯度边缘点8梯度的负方向作直线与坐标网格相交于一点，从这个交点处开始沿着该交点梯度的负方向作直线得到下一个交点，重复以上过程直到梯度模不再减小，从而得到沿着梯度负方向的梯度剖面。

最终调强放射治疗计划剂量数据平面上梯度边缘点8追踪形成的梯度剖面7由上述沿着梯度正方向追踪得到的梯度剖面和沿着梯度负方向追踪得到的梯度剖面两部分共同组成。

上述过程中，追踪形成的梯度剖面7与调强放射治疗计划剂量数据坐标网格的各个交点中，交点6为已存在的点，交点9不为已存在的数据点，交点9的梯度由与之相邻点梯度的线性插值得到；

(4)梯度边缘点对应梯度剖面的锐度和非梯度边缘点自定义的偏离系数的计算对每一个调强放射治疗计划剂量数据平面上梯度边缘点(i′,j′)，即图3中的调强放射治疗计划剂量数据平面上梯度边缘点8对应的追踪形成的梯度剖面7按如下公式(7)计算锐度σ(i′,j′)(GradientProfilePriorandItsApplicationsinImageSuper-ResolutionandEnhancement[J]，IEEETransactionsOnImageProces-sing，vol.20，num.6.pp.1529-1542，2011)：

$\mathrm{\σ}(i^{\′},j^{\′})=\sqrt{{\mathrm{\Σ}}_{x\∈P(i^{\′},j^{\′})}\frac{m\left(x\right)}{M(i^{\′},j^{\′})}{d}_c^2(x,(i^{\′},j^{\′}\left)\right)}---\left(7\right)$

其中，x是调强放射治疗计划剂量数据平面上梯度边缘点(i′,j′)所对应梯度剖面P(i′,j′)上与调强放射治疗计划剂量数据坐标网格相交的每一点，m(x)表示x点处的梯度模，M(i′,j′)表示梯度剖面P(i′,j′)上梯度模的和，d_c(x,(i′,j′))表示点x与梯度边缘点(i′,j′)之间的曲线距离；

本实施例的插值算法中自定义的调强放射治疗计划剂量数据平面上非梯度边缘点(i″,j″)的偏离系数ρ(i″,j″)为每一点与周围3×3邻域整体的平均计划剂量值偏差，以此来度量该点的偏离程度，如以下公式(8)所示：

$\mathrm{\ρ}(i^{\′\′},j^{\′\′})=\frac{|\stackrel{\‾}{I_L(i^{\′\′},j^{\′\′})}-I_L(i^{\′\′},j^{\′\′})|}{I_L(i^{\′\′},j^{\′\′})}---\left(8\right)$

其中为以非梯度边缘点(i″,j″)为中心的3×3邻域中除非梯度边缘点(i″,j″)点以外各个点的平均计划剂量值。I_L(i″,j″)为该非梯度边缘点(i″,j″)的计划剂量值。

(5)确定双三次插值核的系数

本实施例对双三次插值核的系数对梯度剖面的调节作用进行研究和分析。

所述传统双三次插值利用待插值点周围的16个点作为基准点，使用如下公式(9)(贾永红，《数字图像处理》武汉大学出版社，2003.)得到待插值点的剂量值：

f(i+u,j+v)＝ABC(9)

式中，f(i_C+u,j_C+v)表示与双三次插值的中心点(i_C,j_C)在水平方向和垂直方向距离分别为u和v的待插值点的计划剂量值，A、B、C为三个矩阵，其具体表示如下式：

A＝[s(1+u)s(u)s(1-u)s(2-u)](10)

$B=\left[\begin{array}{cccc}f(i_C-1,j_C-1)& f(i_C-1,j_C)& f(i_C-1,j_C+1)& f(i_C-1,j_C+2)\\ f(i_C,j_C-1)& f(i_C,j_C)& f(i_C,j_C+1)& f(i_C,j_C+2)\\ f(i_C+1,j_C-1)& f(i_C+1,j_C)& f(i_C+1,j_C+1)& f(i_C+1,j_C+2)\\ f(i_C+2,j_C-1)& f(i_C+2,j_C)& f(i_C+2,j_C+1)& f(i_C+2,j_C+2)\end{array}\right]---\left(11\right)$

C＝[s(1+v)s(v)s(1-v)s(2-v)]^T(12)

其中u和v分别表示水平方向和垂直方向上待插值点与双三插值的中心点(i_C,j_C)的距离，f(i_C,j_C)表示坐标为(i_C,j_C)的点的计划剂量值，s(w)表示双三次插值核。

传统的双三次插值为了逼近最佳插值函数sin(w)/w，其对应的插值核s(w)表达式如下(13)式：

$s\left(w\right)=\left\{\begin{array}{cc}{\left|w\right|}^3-2{\left|w\right|}^2+1,& 0\&le;\left|w\right|<1\\ -{\left|w\right|}^3+5{\left|w\right|}^2-8\left|w\right|+4,& 1\&le;\left|w\right|<2\\ 0,& \left|w\right|\&GreaterEqual;2\end{array}\right.---\left(13\right)$

其中w项的实际意义为待插值点到每一个基准点之间的距离，因此s(w)表示16个基准点按照与待插值点之间的距离对应的加权系数。

而RobertGK(ConvolutionInterpolationforDigitalImageProcessing[J]，IEEETransactionsOnAcoustics，Speech，AndSignalProcessing，vol.29，num.6.pp.1153-1160，1981)提出双三次插值核s(w)可以有如下公式(14)，当a<0时：

$s\left(w\right)=\left\{\begin{array}{cc}(a+2){\left|w\right|}^3-(a+3){\left|w\right|}^2+1\&GreaterEqual;0,& 0\&le;\left|w\right|<1\\ a{\left|w\right|}^3-5a{\left|w\right|}^2+8a\left|w\right|-4a\&le;0,& 1\&le;\left|w\right|<2\\ 0,& \left|w\right|\&GreaterEqual;2\end{array}\right.---\left(14\right)$

可以看出，当a＝-1时就是传统的双三次插值公式对应的插值核。

以下为本实施例对双三次插值核的系数调节作用的研究分析：

假设有两个双三次插值核的系数a₁和a₂且满足a₁<a₂<0，则对应的双三次插值核之差为下式：

$s_{a1}\left(w\right)-s_{a2}\left(w\right)=\left\{\begin{array}{cc}(a_1-a_2)w^2(w-1)\&GreaterEqual;0,& 0\&le;\left|w\right|<1\\ (a_1-a_2)(w-1){(w-2)}^2\&le;0,& 1\&le;\left|w\right|<2\\ 0,& \left|w\right|\&GreaterEqual;2\end{array}\right.---\left(15\right)$

从公式(15)可以看出0≤|w|＜1时，双三次插值核的系数a绝对值越小，双三次插值核值越小，此时双三次插值的中心点(i_C,j_C)所对应的权重越小，从而待插值点与双三次插值的中心点(i_C,j_C)的差异越大，因而中心点(i_C,j_C)处容易形成较大的计划剂量的梯度模；从公式(15)可以看出，1≤|w|＜2时，双三次插值核的系数a绝对值越小，双三次插值核的值越大，而且由公式(14)可知此时双三次插值核的值小于0，故此时双三次插值核的绝对值越接近于0，此时双三插值的中心点(i_C,j_C)之外的其他基础点对待插值点的影响越小，从而双三次插值核的系数a对于待插值点与双三次插值的中心点(i_C,j_C)之间在0≤|w|＜1时的调节作用越明显，从而可由双三次插值核的系数a在中心点(i_C,j_C)处调节出较大幅度的计划剂量的梯度模。

从图4中可以看出，双三次插值核的系数为-0.5时双三次插值核的变化趋势11，双三次插值核的系数为-1时双三次插值核的变化趋势12和双三次插值核的系数为-1.5时双三次插值核的变化趋势13所表示的双三次插值核变化趋势中，双三次插值核的系数绝对值越小所对应双三次插值核变化越快，且双三次插值核s(w)的函数值越小，因此梯度边缘点所对应的权重越小，从而待插值点与梯度边缘点的差异也越大，；因而在中心点(i_C,j_C)处容易形成较大幅度的计划剂量的梯度模。

在实际过程中，如图5(a)所示，对于调强放射治疗计划剂量数据平面上每一个梯度边缘点8所对应的梯度剖面P1或梯度剖面P2，用双线性插值所拟合出的剂量数据分布16由于无法真实反映调强放射治疗计划剂量数据平面上梯度边缘点对应梯度剖面P1的真实剂量数据分布14，或调强放射治疗计划剂量数据平面上梯度边缘点对应梯度剖面P2的真实剂量数据分布15，从而使梯度边缘附近的插值计算产生较大的误差。结合图5(b)可以看出，梯度剖面P2的梯度模变化比梯度剖面P1快，而待插值点处双线性插值的计算结果17与梯度剖面P2在待插值点处的真实剂量数据19之间的差异大于其与梯度剖面P1在待插值点处的真实剂量数据18之间的差异，从而可以看出梯度模变化越快的梯度剖面对应的误差越大。

而从调强放射治疗计划剂量数据平面上梯度边缘点(i′,j′)对应梯度剖面的锐度σ(i′,j′)的定义可以看出，锐度的值越小，对应剖面上的梯度模变化越快，从图5可以看出，待插值点的值会偏离梯度边缘点(i′,j′)而接近于两旁的点，结合对图4以及公式(14)和公式(15)的分析可以看出，若以该梯度边缘点(i′,j′)为双三次插值中心点(i_C,j_C)，则需要双三次插值核的系数的绝对值越小；通过这样的处理，可以使调强放射治疗计划剂量数据平面上梯度边缘点对应梯度剖面的剂量数据分布更加接近于真实的剂量数据分布；从自定义的非梯度边缘点(i″,j″)的偏离系数来看，偏离系数ρ(i″,j″)越大则待插值点与非梯度边缘点(i″,j″)之间的差异越大，结合对图4以及公式(14)和公式(15)的分析可以看出，若以该非梯度边缘点(i″,j″)为双三次插值中心点(i_C,j_C)，则需要双三次插值核的系数的绝对值越小。

通过以上的研究分析结果拟合出锐度或偏离程度与每一点双三次插值核的系数之间的函数关系。这就是本插值算法使用双三次插值核的系数进行调节的思路。

(6)根据前面所述对双三次插值核的系数调节作用的分析，本插值算法对调强放射治疗计划剂量数据平面上每一个梯度边缘点和非梯度边缘点分别构造函数以确定以每一点为双三次插值的中心点时对应的双三次插值核的系数。

对于调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)，以点(i′,j′)为双三次插值中心点时，对双三次插值核的系数a(i′,j′)与该梯度边缘点对应梯度剖面的锐度σ(i′,j′)之间构造函数关系如下式所示：

$a(i^{\&prime;},j^{\&prime;})=-0.5\frac{1}{1+ln({\mathrm{\&sigma;}}_{max}/\mathrm{\&sigma;}(i^{\&prime;},j^{\&prime;}\left)\right)}\exp \left\{{(1-\frac{\mathrm{\&sigma;}(i^{\&prime;},j^{\&prime;})}{{\mathrm{\&sigma;}}_{max}})}^2\right\}---\left(16\right)$

其中，σ_max表示该调强放射治疗计划剂量数据平面上锐度的最大值；可以看出，式(16)中双三次插值核的系数a(i′,j′)的绝对值随锐度的减小而减小，满足前面得出的规律，即双三次插值核的值随双三次插值核的系数的绝对值减小而减小。

对于调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)，以点(i″,j″)为双三次插值中心点时，对双三次插值核的系数a(i″,j″)与该非梯度边缘点的偏离系数ρ(i″,j″)之间构造函数关系如下式所示：

$a(i^{\&prime;\&prime;},j^{\&prime;\&prime;})=-0.5\exp \{-{\left(\frac{\mathrm{\&rho;}(i^{\&prime;\&prime;},j^{\&prime;\&prime;})-{\mathrm{\&rho;}}_{min}}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{min}}\right)}^2\}---\left(17\right)$

其中ρ_max为调强放射治疗计划剂量数据平面上最大偏离系数，ρ_min为调强放射治疗计划剂量数据平面上最小偏离系数。

将(16)式和(17)式两者结果合并，从而得出以调强放射治疗计划剂量数据平面上每一点(i,j)为双三次插值中心点时所对应的双三次插值核的系数a(i,j)；

(7)待插值点调强放射治疗计划剂量数据的获得

从以上过程可知，对于调强放射治疗计划剂量数据平面上梯度边缘点和非梯度边缘点可以分别使用公式(16)和公式(17)所对应的函数关系确定以调强放射治疗计划剂量数据平面上每一点(i,j)为双三次插值中心点时对应的双三次插值核的系数。因此，对于调强放射治疗计划剂量数据平面上每个点(i,j)，都可以得出以点(i,j)为双三次插值中心点时对应的双三次插值核的系数a(i,j)。对于每一个待插值点，设其坐标为(a,b)，若a∈[i,i+1)且b∈[j,j+1)，则以点(i,j)为双三次插值中心点，以a(i,j)为双三次插值核的系数，即可得到该点所对应的双三次插值核如下式所示：

$s\left(w\right)=\left\{\begin{array}{cc}\left(a\right(i,j)+2){\left|w\right|}^3-\left(a\right(i,j)+3){\left|w\right|}^2+1\&GreaterEqual;0,& 0\&le;\left|w\right|<1\\ a(i,j){\left|w\right|}^3-5a(i,j){\left|w\right|}^2+8a(i,j)\left|w\right|-4a\&le;0,& 1\&le;\left|w\right|<2\\ 0,& \left|w\right|\&GreaterEqual;2\end{array}\right.---\left(18\right)$

其中w项的实际意义为待插值点到每一个基准点之间的距离，将公式(18)带入公式(10)、公式(11)和公式(12)中，以点(i,j)为双三次插值中心点对该待插值点进行双三次插值计算，从而得出每一个待插值点的调强放射治疗计划剂量数据。

(7)所得结果验证

为了验证本发明所述插值算法的正确性，最直接的办法是将各个待插值点处的插值所得计划剂量数据与该点处的计划剂量数据的真实值进行比.较，然而治疗计划系统(TPS)对同一个计划只能产生一个原始调强放射治疗计划剂量数据平面，各个待插值点的计划剂量计算值是经过线性插值得到的，因此TPS系统无法得到待插值点的计划剂量数据的真实值。但是蒙特卡罗程序(BEAMnrc/Dosxyznrc)作为剂量计算的“金标准”可以用来模拟真实的剂量分布，因此本实施例中使用蒙特卡罗程序BEAMnrc/Dosxyznrc替代TPS产生原始计划剂量数据和待插值点的计划剂量数据的真实值。分别使用基于梯度特征的调强放射治疗计划剂量数据插值算法以及双线性插值算法对原始调强放射治疗计划剂量数据进行计算，得到各个待插值点处的插值所得计划剂量数据，再与蒙特卡罗程序BEAMnrc/Dosxyznrc输出的待插值点处计划剂量数据的真实值进行对比。

对于一个总大小为96mm*96mm调强放射治疗照射区域，以0.4mm到4.0mm范围内可以被96mm整除的物理间距有0.4mm、0.48mm、0.5mm、0.6mm、0.64mm、0.75mm、0.8mm、0.96mm、1.0mm、1.2mm、1.5mm、1.6mm、2.0mm、2.4mm、3.0mm、3.2mm、4.0mm，共17组不同的物理间距作为原始调强放射治疗计划剂量数据平面的物理间距，模拟出该照射区域中17组原始调强放射治疗计划剂量数据平面的计划剂量数据分布；在相同环境和条件下，模拟出17组不同物理间距的待插值点平面上各个待插值点处的计划剂量数据的真实值。在每一个原始调强放射治疗计划剂量数据平面中分别运用基于梯度特征的调强放射治疗计划剂量数据的插值算法、双线性插值算法和传统双三次插值算法，得出17种插值间距下的插值结果，总共形成17×17组结果，并与模拟得出的待插值平面上各点的计划剂量数据的真实值进行比较，得出三种算法各自的误差。

图6展示了用17组原始调强放射治疗计划剂量数据平面以0.4mm为插值间距作插值计算时三种方法的误差对比；图7展示了用17组原始调强放射治疗计划剂量数据平面以0.48mm为插值间距作插值计算时三种方法的误差对比；图8展示了用17组原始调强放射治疗计划剂量数据平面以0.5mm为插值间距作插值计算时三种方法的误差对比；图9展示了用17组原始调强放射治疗计划剂量数据平面以0.6mm为插值间距作插值计算时三种方法的误差对比。在每一个原始调强放射治疗计划剂量数据平面下分别将基于梯度特征的调强放射治疗计划剂量数据的插值算法误差以及双三次插值的误差减去双线性插值误差，从结果可以看到，基于梯度特征的调强放射治疗计划剂量数据插值算法的误差减去双线性插值误差的结果基本都小于0，也就是说，基于梯度特征的调强放射治疗计划剂量数据插值算法的误差小于双线性插值方法，而且可以看出而传统的双三次插值则在很多情况下误差还大于双线性插值。

对17组原始调强放射治疗计划剂量数据平面每组都使用上述三种插值方法分别进行17种间距的插值计算，求出每一个原始调强放射治疗计划剂量数据平面在不同插值间距下的误差平均值，将其误差平均值进行对比，结果展示在图10中，从图10中可以看出，除了原始调强放射治疗计划剂量数据物理间距为1.2mm和1.5mm时，基于梯度特征的调强放射治疗计划剂量数据的插值算法误差大于双线性插值，其余情况下基于梯度特征的调强放射治疗计划剂量数据的插值算法误差都小于双线性插值，而且每种情况下基于梯度特征的调强放射治疗计划剂量数据插值算法的误差都小于传统双三次插值结果；而传统的双三次插值结果误差基本都大于双线性插值。因而可以看到基于梯度特征的调强放射治疗计划剂量数据插值算法的结果在误差控制方面较双线性插值和传统的双三次插值均有一定优势。

再比较基于梯度特征的调强放射治疗计划剂量数据插值算法和双线性插值算法得出平面剂量的平均梯度模的值。图11展示了对17组原始计划剂量数据每组分别进行17种间距的插值计算后，每种原始调强放射治疗计划剂量数据间距下平均梯度模的对比，可以发现各个原始调强放射治疗计划剂量数据物理间距下基于梯度特征的调强放射治疗计划剂量数据插值算法得出结果的平均梯度模的值均大于双线性算法，同时可以发现，随着原始调强放射治疗计划剂量数据物理间距的增加，双线性算法对于梯度的平滑效应越发严重，而基于梯度特征的调强放射治疗计划剂量数据插值算法则良好的维持了数据平面的梯度模，这种梯度上的改善一方面可以为物理师作判断时提供良好的梯度信息。另一方面，有些剂量验证软件如MapCkeck在将测量剂量数据与TPS输出的调强放射治疗计划剂量数据对齐时会采用双线性插值，也就是说在TPS中对调强放射治疗计划剂量数据进行了插值之后，剂量验证软件还会对计划剂量数据再一次进行插值计算，此时若使用双线性插值，则会产生多次平滑效应；由于基于梯度特征的调强放射治疗计划剂量数据插值算法对计划剂量数据在梯度保持方面的效果较好，因此可以有效的减小这种多次插值所带来的多次平滑效应。

对于基于梯度特征的调强放射治疗计划剂量数据插值算法、双线性插值算法和传统双三次插值算法这三种方法在所有的情况，即17×17组实验情况下的平均误差以及平均梯度模进行计算，其结果展示在表1中。可以看出，基于梯度特征的调强放射治疗计划剂量数据插值算法的整体平均误差小于双线性插值和双三次插值；基于梯度特征的调强放疗计划剂量数据插值算法整体平均梯度模大于双线性插值，略小于传统双三次插值。同时可以看出，改进Canny边缘探测算法后进一步减小了误差，这样可以在提高准确度的同时减小梯度边缘点的数量从而减小计算复杂度。实际结果表明，从总体上来看，基于梯度特征的调强放射治疗计划剂量数据插值算法兼具误差控制方面和梯度保持方面的优势。

表1所述几种插值法对所有组别(17×17组)结果的平均值

Claims (4)

1.一种基于梯度特征的调强放射治疗计划剂量数据的插值算法，其特征在于根据计划剂量数据平面梯度信息，使用改进传统Canny边缘探测算法获取其上各梯度边缘点和非梯度边缘点，以获取的梯度边缘点对应的梯度剖面锐度及非梯度边缘点对应的偏离系数，得到以计划剂量数据平面上每一点为双三次插值中心点时对应的双三次插值核系数，对每一待插值点进行双三次插值，即得每一待插值点的调强放射治疗计划剂量数据；包括以下步骤：

步骤1：对治疗计划系统产生的调强放射治疗计划剂量数据平面I_L，计算该调强放射治疗计划剂量数据平面上每一点(i,j)的梯度以及该点对应的梯度模

步骤2：对传统Canny边缘探测算法进行改进，根据步骤1中得到的调强放射治疗计划剂量数据平面上每一点的梯度信息，使用改进的Canny边缘探测算法获取调强放射治疗计划剂量数据平面上各个梯度边缘点(i′,j′)和非梯度边缘点(i″,j″)；

步骤3：对步骤2所获取的调强放射治疗计划剂量数据平面上各个梯度边缘点(i′,j′)进行梯度剖面的追踪，从而得到每一个梯度边缘点对应的梯度剖面；

步骤4：对步骤3所得到的调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)所对应梯度剖面的锐度σ(P(i′,j′))进行计算；对步骤2所述调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)自定义的偏离系数ρ(i″,j″)进行计算；

步骤5：根据步骤4得到调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)对应梯度剖面的锐度σ(P(i′,j′))，确定锐度σ(P(i′,j′))与以梯度边缘点(i′,j′)为双三次插值中心点时所对应双三次插值核的系数a(i′,j′)之间的函数关系；

步骤6：根据步骤4得到调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)的偏离系数ρ(i″,j″)，确定偏离系数ρ(i″,j″)与以非梯度边缘点(i″,j″)为双三次插值中心点时所对应的双三次插值核的系数a(i″,j″)之间的函数关系；

步骤7：将步骤5和步骤6所得到的双三次插值核的系数的两函数结果合并，从而得出以调强放射治疗计划剂量数据平面上每一点(i,j)为双三次插值中心点时对应的双三次插值核的系数a(i,j)；

步骤8：对于待插值点，设置其坐标为(a,b)，若a∈[i,i+1)且b∈[j,j+1)，根据步骤7得出的所述双三次插值核的系数a(i,j)，令点(i,j)为双三次插值中心点，并以该双三次插值核的系数a(i,j)对待插值点进行双三次插值计算，从而得出每一个待插值点的调强放射治疗计划剂量数据。

2.根据权利要求1所述的插值算法，其特征在于所述步骤2中对传统Canny算法进行改进，其采用如下操作：

(a)略去传统Canny算法中高斯滤波的过程，获得略去高斯滤波后的结果；

(b)对略去传统Canny算法中高斯滤波的结果进行判断，假设(i′,j′)点为传统Canny算法略去高斯滤波后结果中的调强放射治疗计划剂量数据平面上的一个梯度边缘点，且以(i′,j′)点为中心的3×3邻域中含有其他梯度边缘点，若(i′,j′)点梯度模小于所述邻域中其他梯度边缘点梯度模，则重新将(i′,j′)点标记为非梯度边缘点。

3.根据权利要求1或2所述的插值算法，其特征在于所述步骤5中调强放射治疗计划剂量数据平面上每一个梯度边缘点(i′,j′)对应梯度剖面的锐度σ(P(i′,j′))与以点(i′,j′)为双三次插值中心点时对应的双三次插值核的系数a(i′,j′)之间的函数关系，由以下公式(1)表示：

$a\left(i^{\&prime;},j^{\&prime;}\right)=-0.5\frac{1}{1+\ln \left({\mathrm{\&sigma;}}_{\max }/\mathrm{\&sigma;}\left(i^{\&prime;},j^{\&prime;}\right)\right)}\exp \left\{{\left(1-\frac{\mathrm{\&sigma;}\left(i^{\&prime;},j^{\&prime;}\right)}{{\mathrm{\&sigma;}}_{\max }}\right)}^2\right\}---\left(1\right)$

其中σ_max表示该调强放射治疗计划剂量数据平面上锐度的最大值。

4.根据权利要求1或3所述的插值算法，其特征在于所述步骤5中调强放射治疗计划剂量数据平面上每一个非梯度边缘点(i″,j″)的偏离系数ρ(i″,j″)与以点(i″,j″)为双三次插值中心点时对应的双三次插值核的系数a(i″,j″)之间的函数关系，由以下公式(2)表示：

$a(i^{\&prime;\&prime;},j^{\&prime;\&prime;})=-0.5\exp \{-{\left(\frac{\mathrm{\&rho;}(i^{\&prime;\&prime;},j^{\&prime;\&prime;})-{\mathrm{\&rho;}}_{min}}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{min}}\right)}^2\}---\left(2\right)$

其中ρ_max为调强放射治疗计划剂量数据平面上最大偏离系数，ρ_min为调强放射治疗计划剂量数据平面上最小偏离系数。