专利名称：：迭代间曲率流滤波加权最小二乘正电子发射断层成像方法
技术领域：
：本发明属于正电子发射计算机断层(PETPositronemissiontomography)成像
技术领域：
，本发明尤其涉及一种用于构建成像图像的迭代间曲率流滤波加权最小二乘正电子发射断层成像方法。
背景技术：
：在现代PET成像系统中，对原始泊松观测数据的随机符合预校正操作，尤其是延迟窗符合校正，已成为必需的一个重要环节。它可以有效地抑制大部分随机符合事件(尤其是背景干扰)提高观测数据的信噪比。然而由于预校正的处理，使得观测数据不在服从原有的泊松统计特性。这样基于原有泊松最大似然的估计方法不在可行。在这种情况下，我们通常考虑加权最小二乘方法。现有最小二乘方法的观测数据模型为yi=Σjaijfj+ϵi---(1)]]>其中yi表示预校正过的第i个探测器通道探测到的光子数，1≤i≤M，M为总的探测器通道个数；fj表示第j个象素出发出的光子数，fj≥0，1≤j≤N，N为总的象素数；aij表示第j个象素出发出的光子能被第i个探测器通道检测到的概率；εi表示附加在第i个探测器通道探测到的光子数目上的随机误差，其统计特性服从零均值，方差为di的正态分布。根据(1)，我们重建的目标函数可以表示为J(f)=argminf(y-Af)TΣ-1(y-Af),]]>以及非负约束条件f≥0(2)这里y由yi，1≤i≤M组成的观测数据向量；f是由fj，1≤j≤N组成的图像向量；A是由aij组成的M×N的系统概率矩阵；∑为M×M加权对角矩阵，其对角线上的第i个元素即为di。加权矩阵的选取非常重要，它的优劣直接影响到重建图像的质量。我们可以采用如下的方差估计方法di＝max(yi，1)，1≤i≤M，(3)在非负约束条件下，如果图像为目标函数(2)的加权最小二乘解，那么必须满足Kuhn-Tucker条件a)---f^jΣi1diaij(yi-Σpaipf^p)=f^jΣi1diaijyi,]]>若f^j>0,j=1,Λ,N---(4)]]>b)---f^jΣi1diaij(yi-Σpaipf^p)≤f^jΣi1diaijyi,]]>若f^j≤0,j=1,Λ,N---(5)]]>根据(4)我们可以直接得到现有的最小二乘迭代方法f^jn+1=f^jnΣiaijyi/diΣiaijΣpaipf^pn/di,j=1,Λ,N,n=0,1,2,Λ---(6)]]>此处n表示迭代次数。由于投影数据本身的不完备性以及重建方法的不适定性，仅使用加权最小二乘方法很大程度上得不到理想的重建图像。其结果容易产生图像不光滑，噪声较明显等现象。一般做法是正则化技术，即同时考虑与图像平滑有关的正则项来抑制噪声的影响。然而其实现的复杂度随正则项的不同而不同。例如，一些正则化方法可以带来良好的重建效果，但实现较为复杂，且平滑力度难以控制。另外正则项的引入并不能较好的抑制背景噪声。为了说明现有技术的缺点，我们根据图1所示的模板数据进行实验仿真。图2给出了利用仿真数据，按照现有加权最小二乘方法重建的结果。我们可以明显地看出重建图像不够理想，噪化现象较为明显。因此如何重建出高质量的PET图像是目前迫切需要解决的问题之一。
发明内容本发明提供一种能够提高重建图像质量的迭代间曲率流滤波加权最小二乘正电子发射断层成像方法。本发明采用如下技术方案一种用于构建成像图像的迭代间曲率流滤波加权最小二乘正电子发射断层成像方法，采用下列步骤1)获取观测数据，记录其规模M，并将它保存在向量y＝[y1，Λ，yM]T中；根据待重建图像的尺寸要求，确定初始图像的规模N，给定它的初始灰度值大于1，并将其保存在向量f＝[f1，Λ，fN]T中，2)根据观测数据y和初始图像f的大小，计算M×N的系统概率矩阵A，3)将观测数据y中的各个分量进行逐个与1比较如果比1大则保留，比1小则设置为1，并将比较结果保存在向量d中，d＝[d1，Λ，dM]T，作为观测数据的误差方差，4)将观测数据y中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的观测投影数据y~:y~=[y~1,Λ,y~M]T,]]>其中y~i=yi/di,]]>i＝1，Λ，M，5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g，g=ATy~,]]>g＝[g1，Λ，gN]T，6)将系统概率矩阵A和初始图像f相乘，得到前向投影p，p＝[p1，Λ，pM]T，再将该前向投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的前向投影p~:p~=[p~1,Λ,p~M]T,]]>其中p~i=pi/di,]]>i＝1，Λ，M，然后利用系统概率矩阵A的转置与该前向投影相乘得到反投影hh=ATp~,]]>h＝[h1，Λ，hN]T，7)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值cc＝[c1，Λ，cN]T，其中cj＝gj/hj，j＝1，Λ，N，8)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得到中间图像f^:f^=[f^1,Λ,f^N]T,]]>其中f^j=cj×fj,]]>j＝1，Λ，N，9)将中间图像进行至少一次曲率流滤波，可得到滤波图像10)将这个滤波图像作为初始图像，返回到第6步，重复这个过程直到重建后的图像收敛，可得到最终成像的图像。本发明结合当前PET成像数据的特点，针对现有方法存在的弊端，提出了一种新的解决方法迭代间曲率流滤波加权最小二乘成像方法。本发明在现有最小二乘方法迭代过程中加入曲率流滤波机制，以此来抑制图像噪化，提高成像图像的质量。与现有的技术相比，本发明具有的突出效果是实现较为简单，其中曲率流滤波过程直接作用在二维图像上，不必进行变换，因此处理速度较快；对所需重建图像的平滑度的控制可以调节滤波的次数得以完成。由于曲率流的滤波是一种良好的图像滤波技术，它对平滑图像的同时并不影响图像原有的基本特征，比如边缘，线，角等，因此重建的图像同样可以保持原始图像中应有的基本特征信息；此外，对抑制重建图像的背景干扰也有着良好的效果。图3给出了本发明对仿真数据进行重建的结果，其中进行3次迭代间滤波处理。将其与图2相对比，我们可以清楚的发现，按本发明所述的方法重建的图像清晰、光滑，能反映原始图像的本来特征。为了对本发明的有效性有一个定量认识，我们用重建图像信噪比准则来衡量这种成像方法的优劣。该准则的定义如下SNR(dB)=10log10(||f||2||f*-f||2)---(7)]]>其中f*表示仿真模板图像数据。重建信噪比就越高，则说明成像方法越好。通过计算，我们将图2和图3的重建信噪比分别列入表1中，以作比较。表1<tablesid="table1"num="001"><tablewidth="528">现有技术本发明重建信噪比(SNR)8.28dB12.53dB</table></tables>显然，用本发明所述方法可以得到较高的信噪比，其值比现有的加权最小二乘方法高出将近4dB。这已经明显地说明其优越于一般的乘性加权最小二乘方法。表2本发明重建信噪比随着滤波的次数增加而增加。但增加的幅度却随着次数的增加而逐渐减少，因此滤波次数可以控制一定的范围以内，比如10次。此外表明我们发明的方法最终将收敛到一致的重建结果。为了进一步验证我们发明的方法的优越性，我们用实际获取的3个临床数据加以说明。图4、图5、图6是用现有技术重建结果。很明显，图像质量较差，不便于临床研究。图7、图8和图9分别对应图4、图5和图6，它们是用本发明的方法成像的结果。从中我们可以清晰地发现，本发明的迭代间曲率流滤波加权最小二乘方法成像质量相当优越，它不仅平滑了图像，保持了边缘，而且对抑制背景干扰也相当有效。因此本发明更适于临床诊断需要。图1为用来测试成像方法的胸腔模板图像图2为用现有加权最小二乘方法重建的结果图3为用本发明成像后的结果(3次滤波)图4为用现有加权最小二乘方法对临床数据1进行成像的结果图5为用现有加权最小二乘方法对临床数据2进行成像的结果图6为用现有加权最小二乘方法对临床数据3进行成像的结果图7为用本发明对临床数据1进行成像的结果(5次滤波)图8为用本发明对临床数据2进行成像的结果(5次滤波)图9为用本发明对临床数据3进行成像的结果(5次滤波)具体实施例方式本发明的具体实施方式如下1)获取已进行随机符合校正的观测数据，确定其规模M，并保存为在向量y＝[y1，Λ，yM]T。根据待重建图像的尺寸要求，确定初始图像的规模N，给定它的初始灰度值大于1，并将其保存在向量f＝[f1，Λ，fN]T中。2)根据观测数据y和初始图像f的大小，计算M×N的系统概率矩阵A。为了计算该系统概率矩阵，我们应已知扫描仪的几何结构，然后可以按照一定的方法，例如视角法(Angleofview)进行计算。在实际情况中，观测数据采集自PET显像扫描仪。为了评估方法的优劣，我们也可以通过计算机仿真实验来生成需要的观测数据。这里我们用图1所示的一个计算机仿真的PET胸腔模板来说明。该胸腔模板大小为128×128象素矩阵，即N＝128×128。我们假定投影数据规模为M＝192×192，它表示192个投影角度，每个角度上有192个探测器通道。我们进一步假设通道与通道间的间距为一个象素，这样便于我们计算系统概率矩阵A。随机符合校正的方法很多，现代PET扫描仪均采用延迟窗符合校正。按其原理，我们可以作如下仿真实验以生成所需的观测数据。我们先用系统概率矩阵与胸腔模板图像相乘得到无噪声的观测数据，然后将所有的数据求和得到总的光子数。将总光子数除以总的探测器通道数M，估算出每个探测器通道所捕获的平均光子数。取该平均值的25％作为背景随机符号均值(25％是根据临产数据的分析所得的经验值)。我们将每一个无噪声观测数据加上该背景随机符号均值作为仿真用的泊松观测投影数据的均值。用该均值结合计算机上的泊松伪随机数生成器，产生所需的未校正的观测投影数据。用同样的背景随机符号均值，结合伪随机数生成器，我们可以模拟出另外一组泊松数据作为延迟窗随机符合数据。将未校正的观测投影数据减去该延迟窗随机符号数据得到最终的预校正的观测投影数据y。3)利用观测投影数据y，结合公式(3)所给的方法计算观测数据的误差方差d＝[d1，Λ，dM]T。4)将观测数据y中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的观测投影数据y~:y~=[y~1,Λ,y~M]T,]]>其中y~i=yi/di,]]>i＝1，Λ，M，5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g，即g=ATy~,]]>g＝[g1，Λ，gN]T，6)将系统概率矩阵A和初始图像f相乘，得到前向投影p，p＝[p1，Λ，pM]T，再将该前向投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的前向投影p~:p~=[p~1,Λ,p~M]T,]]>其中p~i=pi/di,]]>i＝1，Λ，M。然后利用系统概率矩阵A的转置与该前向投影相乘得到反投影hh=ATp~,]]>h＝[h1，Λ，hN]T，7)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值cc＝[c1，Λ，cN]T，其中cj＝gj/hj，j＝1，Λ，N，8)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得到中间图像f^:f^=[f^1,Λ,f^N]T,]]>其中f^j=cj×fj,]]>j＝1，Λ，N，9)对中间图像进行曲率流滤波曲率流滤波能有效的去除噪声，但不损害图像原有的特征信息，例如边缘、线、角等。它的滤波过程可以用如下的偏微分方程表示∂f(h,v;t)∂t=▿·(▿f(h,v;t)|▿f(h,v;t)|)---(7)]]>以及边界条件初始条件f(h，v；0)＝f0。这里表示梯度算子；为了便于区别，我们用h，v表示图像的水平和垂直两个方向的坐标；Ω为图像定义域；Ω为图像边界，为边界法向量；t表示时间。引入时间进化步长Δt(这里取Δt＝0.0005较为合适)并将方程(7)右边项展开得ft+1=ft+Δt×(fhht(fvt)2-2fhvtfhtfvt+fvvt(fht)2(fht)2+(fvt)2)---(8)]]>其中ft≡f(h，v；t)；fht≡∂f(h,v;t)/∂h,fvt≡∂f(h,v;t)/∂v]]>为t时刻函数f的水平和垂直两个方向的一阶偏导数；fhht≡∂2f(h,v;t)/∂h2,fhvt≡∂f(h,v;t)/∂h∂v]]>和fvvt≡∂2f(h,v;t)/∂v2]]>为二阶偏导数。在实际情况中，我们通常考虑离散图像。为了运用公式(8)进行滤波，我们需要引入合适的差分格式来替换公式(8)中的偏导数。为了方便叙述，我们可以假设图像为K×K(K×K＝N)的方阵且象素间距为单位长度1。记(fht)j为函数f(h，v；t)在点(hj，vj)(其中1≤hj≤K，1≤vj≤K)处的h方向的一阶偏导数，则它可以用如下一阶差分表示(fht)j=f(hj+1,vj;t)-f(hj,vj;t),j=1,Λ,N---(9)]]>其中hj＝int[j/K]+1，vj＝j-int[j/N]×N(int[x]表示取x的整数部分)以及j＝(hj-1)×K+vj。按同样的方法，我们还可以得到(fvt)j=f(hj,vj+1;t)-f(hj,vj;t),---(10)]]>(fhht)j=f(hj+1,vj;t)-2f(hj,vj;t)+f(hj-1,vj;t),---(11)]]>(fvvt)j=f(hj,vj+1;t)-2f(hj,vj;t)+f(hj,vj-1;t),---(12)]]>(fhvt)j=(f(hj+1,vj+1;t)+f(hj-1,vj-1;t))-(f(hj-1,vj+1;t)+f(hj+1,vj-1;t))4---(13)]]>同理令fjt=f(hj,vj;t),]]>j＝1，Λ，N。这样我们可以将(8)改写为fjt+1=fjt+Δt×(fhht)j(fvt)j2-2(fhvt)j(fht)j(fvt)j+(fvvt)j(fht)j2(fht)j2+(fvt)j2,j=1,Λ,N,---(14)]]>以及t＝0，Λ，tmax，这里tmax为总的时间即滤波次数。对于图像边缘，我们做边缘复制处理，即对于公式(9)至公式(13)若hj-1≤0，则置hj-1为1；若vj-1≤0，则置vj-1为1；若hj+1≥K，则置hj+1为K；若vj+1≥K，则置vj+1为K。根据上述的曲率流滤波原理，对中间图像进行曲率流滤波可以如下操作设置总的滤波次数tmax(一般1≤tmax≤10较为合理)；设置初始条件fj0=f^j,]]>j＝1，Λ，N，并置t＝0；按公式(14)进行滤波，置t＝t+1；重复公式(14)直至t＝tmax结束。将最终滤波结果保存在向量中，f~=[f~1,Λ,f~N]T.]]>10)将这个滤波图像作为初始图像f，返回到第6步，重复这个过程直到重建后的图像收敛，可得到最终成像的图像。权利要求1.一种用于构建成像图像的迭代间曲率流滤波加权最小二乘正电子发射断层成像方法，其特征在于采用下列步骤1)获取观测数据，记录其规模M，并将它保存在向量y＝[y1，Λ，yM]T中；根据待重建图像的尺寸要求，确定初始图像的规模N，给定它的初始灰度值大于1，并将其保存在向量f＝[f1，Λ，fN]T中，2)根据观测数据y和初始图像f的大小，计算M×N的系统概率矩阵A，3)将观测数据y中的各个分量进行逐个与1比较如果比1大则保留，比1小则设置为1，并将比较结果保存在向量d中，d＝[d1，Λ，dM]T，作为观测数据的误差方差，4)将观测数据y中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的观测投影数据y~=[y~1,Λ,y~M]T,]]>其中y~i=yi/di,]]>i＝1，Λ，M，5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g，g=ATy~,]]>g＝[g1，Λ，gN]T，6)将系统概率矩阵A和初始图像f相乘，得到前向投影p，p＝[p1，Λ，pM]T，再将该前向投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的前向投影p~=[p~1,Λ,p~M]T,]]>其中p~i=pi/di,]]>i＝1，Λ，M，然后利用系统概率矩阵A的转置与该前向投影相乘得到反投影hh=ATp~,]]>h＝[h1，Λ，hN]T，7)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值cc＝[c1，Λ，cN]T，其中cj＝gj/hj，j＝1，Λ，N，8)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得到中间图像f^=[f^1,Λ,f^N]T,]]>其中f^j=cj×fj,]]>j＝1，Λ，N，9)将中间图像进行至少一次曲率流滤波，可得到滤波图像10)将这个滤波图像作为初始图像，返回到第6步，重复这个过程直到重建后的图像收敛，可得到最终成像的图像。全文摘要本发明公开了一种用于构建成像图像的迭代间曲率流滤波加权最小二乘正电子发射断层成像方法，首先获取观测数据，确定初始图像的规模，将观测数据中的各个分量进行逐个与1比较，与误差方差中的各个分量逐个对应相除等后将中间图像进行曲率流滤波，可得到作为初始图像的滤波图像，重复这个过程直到重建后的图像收敛，可得到最终成像的图像；本发明结合当前PET成像数据的特点，在现有最小二乘方法迭代过程中加入曲率流滤波机制，以此来抑制图像噪化，提高成像图像的质量，实现较为简单，处理速度较快，对抑制重建图像的背景干扰也有着良好的效果。文档编号G06T5/00GK1799506SQ200510122610公开日2006年7月12日申请日期2005年11月29日优先权日2005年11月29日发明者周健,罗立民,李松毅申请人:东南大学