首页 / 专利库 / 诊断设备和程序 / 正电子发射断层成像 / 迭代间曲率流滤波加权最小二乘正电子发射断层成像方法

迭代曲率流滤波加权最小二乘电子发射断层成像方法

阅读:301发布:2020-05-19

专利汇可以提供迭代曲率流滤波加权最小二乘电子发射断层成像方法专利检索,专利查询,专利分析的服务。并且本 发明 公开了一种用于构建成像图像的 迭代 间 曲率 流滤波加权最小二乘 正 电子 发射 断层 成像方法,首先获取观测数据,确定初始图像的规模,将观测数据中的各个分量进行逐个与1比较,与误差方差中的各个分量逐个对应相除等后将中间图像进行曲率流滤波,可得到作为初始图像的滤波图像,重复这个过程直到重建后的图像收敛,可得到最终成像的图像;本发明结合当前PET成像数据的特点,在现有最小二乘方法迭代过程中加入曲率流滤波机制,以此来抑制图像噪化,提高成像图像的 质量 ,实现较为简单,处理速度较快,对抑制重建图像的背景干扰也有着良好的效果。,下面是迭代曲率流滤波加权最小二乘电子发射断层成像方法专利的具体信息内容。

1.一种用于构建成像图像的迭代曲率流滤波加权最小二乘电子发射断层成像方法,其 特征在于采用下列步骤:
1)获取观测数据,记录其规模M,并将它保存在向量y=[y1,…,yM]T中;根据待重建图 像的尺寸要求,针对正方形图像,设它的边长为K,确定初始图像的规模N=K×K,给 定它的初始灰度值大于1,并将其保存在向量f=[f1,…,fN]T中,设定迭代的总次数I, I=3,
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 = y i / d i , i=1,…,M,
5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g, g = A T y ~ , g=[g1,…,gN]T,
6)设置m=0,
7)将系统概率矩阵A和初始图像f相乘,得到前向投影p,p=[p1,…,pM]T,再将该前向 投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的前向投影: p ~ = [ p ~ 1 , . . . , p ~ M ] T , 其中 p ~ i = p i / d i , i=1,…,M,然后利用系统概率矩阵A的转置与该前向 投影相乘得到反投影h: h = A T p ~ , h=[h1,…,hN]T,
8)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值c: c=[c1,…,cN]T,其中cj=gj/hj,j=1,…,N,
9)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得到中间图 像 f ^ : f ^ = [ f ^ 1 , . . . , f ^ N ] T , 其中 f ^ j = c j × f j , j=1,…,N,
10)设置t=0,以及滤波总次数tmax,1≤tmax≤10,
11)将中间图像重新记作一个新的临时图像向量 f ^ t : f ^ t = [ f ^ l t , . . . , f ^ N t ] T ,
12)用该临时图像计算如下导数向量:
f ^ h t = [ ( f h t ) 1 , . . . , ( f h t ) N ] T , f ^ v t = [ ( f v t ) 1 , . . . , ( f v t ) N ] T ,
f ^ hh t = [ ( f hh t ) 1 , . . . , ( f hh t ) N ] T , f ^ vv t = [ ( f vv t ) 1 , . . . , ( f vv t ) N ] T ,
f ^ hv t = [ ( f hv t ) 1 , . . . , ( f hv t ) N ] T ,
其中
( f h t ) j = f ^ h j × K + v j t - f ^ ( h j - 1 ) × K + v j t ,
( f v t ) j = f ^ ( h j - 1 ) × K + v j + 1 t - f ^ ( h j - 1 ) × K + v j t ,
( f hh t ) j = f ^ h j × K + v j t - 2 × f ^ ( h j - 1 ) × K + v j t + f ^ ( h j - 2 ) × K + v j t ,
( f vv t ) j = f ^ ( h j - 1 ) × K + v j + 1 t - 2 × f ^ ( h j - 1 ) × K + v j t + f ^ ( h j - 1 ) × K + v j - 1 t ,
( f hv t ) j = 1 4 × ( f ^ h j × K + v j + 1 t + f ^ ( h j - 2 ) × K + v j - 1 t - f ^ ( h j - 2 ) × K + v j + 1 t - f ^ h j × K + v j - 1 t ) ,
j=1,…,N,这里的hj和vj用下面两式计算:
hj=int[j/K]+1,vj=j-int[j/k]×K
此处的int[j/K]表示取j除以K的商,
其中ft≡f(h,v;t); f h t f ( h , v ; t ) / h , f v t f ( h , v ; t ) / v 为t时刻函数f的平和垂直两个 方向的一阶偏导数; f hh t 2 f ( h , v ; t ) / h 2 , f hv t f ( h , v ; t ) / h v f vv t 2 f ( h , v ; t ) / v 2 为二阶偏导 数,
13)计算中间滤波图像 f ~ t : f ~ t = [ f ~ 1 t , . . . , f ~ N t ] T , 其中
f ~ j t = f ^ j t + Δt × ( f hh t ) j × ( ( f v t ) j ) 2 - 2 × ( f hv t ) j × ( f h t ) j × ( f v t ) j + ( f vv t ) j × ( ( f h t ) j ) 2 ( ( f h t ) j ) 2 + ( ( f v t ) j ) 2 , j=1,…,N
14)将该得到的中间滤波图像赋值给临时图像,并将t的值增加1,
15)如果t的值大于或等于滤波总次数tmax,则进行第16步,否则返回第11步,
16)将临时图像保存为滤波图像,将这个滤波图像作为初始图像,m的值相应增加 1,返回到第7步,重复这个过程直到m达到迭代的总次数I,可得到最终成像的图像。

说明书全文

技术领域

发明属于电子发射计算机断层(PET:Positron emission tomography)成像 技术领域,本发明尤其涉及一种用于构建成像图像的迭代曲率流滤波加权最小 二乘正电子发射断层成像方法。

背景技术

在现代PET成像系统中,对原始泊松观测数据的随机符合预校正操作,尤 其是延迟窗符合校正,已成为必需的一个重要环节。它可以有效地抑制大部分随 机符合事件(尤其是背景干扰)提高观测数据的信噪比。然而由于预校正的处理, 使得观测数据不在服从原有的泊松统计特性。这样基于原有泊松最大似然的估计 方法不在可行。在这种情况下,我们通常考虑加权最小二乘方法。
现有最小二乘方法的观测数据模型为:
y i = Σ j a ij f i + ϵ 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 ) = arg min f ( 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 Σ i 1 d i a ij ( y i - Σ p a ip f ^ p ) = f ^ j Σ i 1 d i a ij y i , f ^ j > 0 , j = 1 , . . . , N - - - ( 4 )
b) f ^ j Σ i 1 d i a ij ( y i - Σ p a ip f ^ p ) f ^ j Σ i 1 d i a ij y i , f ^ j 0 , j = 1 , . . . , N - - - ( 5 )
根据(4)我们可以直接得到现有的最小二乘迭代方法:
f ^ j n + 1 = f ^ j n Σ i a ij y i / d i Σ i a ij Σ p a ip f ^ p n / d i , 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 ~ 1 , . . . , y ~ M ] T , 其中 y ~ i = y i / d i , i=1,…,M,
5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g, g = A T y ~ , g=[g1,…,gN]T,
6)将系统概率矩阵A和初始图像f相乘,得到前向投影p,p=[p1,…pM]T,再 将该前向投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到 修正的前向投影 p ~ = [ p ~ 1 , . . . , p ~ M ] T , 其中 p ~ i = p i / d i , i=1,…,M,然后利用 系统概率矩阵A的转置与该前向投影相乘得到反投影h: h = A T p ~ , h=[h1,…,hN]T,
7)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值 c:c=[c1,…,cN]T,其中cj=gj/hj,j=1,…,N,
8)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得 到中间图像 f ^ = [ f ^ 1 , . . . , f ^ N ] T , 其中 f ^ j = c j × f j , j=1,…,N,
9)将中间图像进行至少一次曲率流滤波,可得到滤波图像
10)将这个滤波图像作为初始图像,返回到第6步,重复这个过程直到重建后 的图像收敛,可得到最终成像的图像。
本发明结合当前PET成像数据的特点,针对现有方法存在的弊端,提出了 一种新的解决方法:迭代间曲率流滤波加权最小二乘成像方法。本发明在现有最 小二乘方法迭代过程中加入曲率流滤波机制,以此来抑制图像噪化,提高成像图 像的质量。
与现有的技术相比,本发明具有的突出效果是:实现较为简单,其中曲率流 滤波过程直接作用在二维图像上,不必进行变换,因此处理速度较快;对所需重 建图像的平滑度的控制可以调节滤波的次数得以完成。由于曲率流的滤波是一种 良好的图像滤波技术,它对平滑图像的同时并不影响图像原有的基本特征,比如 边缘,线,角等,因此重建的图像同样可以保持原始图像中应有的基本特征信息; 此外,对抑制重建图像的背景干扰也有着良好的效果。
图3给出了本发明对仿真数据进行重建的结果,其中进行3次迭代间滤波处 理。将其与图2相对比,我们可以清楚的发现,按本发明所述的方法重建的图像 清晰、光滑,能反映原始图像的本来特征。
为了对本发明的有效性有一个定量认识,我们用重建图像信噪比准则来衡量 这种成像方法的优劣。该准则的定义如下:
SNR ( dB ) = 10 log 10 ( || f || 2 || f * - f || 2 ) - - - ( 7 )
其中f*表示仿真模板图像数据。重建信噪比就越高,则说明成像方法越好。通 过计算,我们将图2和图3的重建信噪比分别列入表1中,以作比较。
表1

显然,用本发明所述方法可以得到较高的信噪比,其值比现有的加权最小二 乘方法高出将近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。为了计 算该系统概率矩阵,我们应已知扫描仪的几何结构,然后可以按照一定的方 法,例如视角法(Angle of view)进行计算。
在实际情况中,观测数据采集自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 ~ 1 , . . . , y ~ M ] T ; 其中 y ~ i = y i / d i , i=1,…,M,
5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g,即: g = A T y ~ , g=[g1,…,gN]T,
6)将系统概率矩阵A和初始图像f相乘,得到前向投影p,p=[p1,…,pM]T,再 将该前向投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到 修正的前向投影 p ~ = [ p ~ 1 , . . . , p ~ M ] T , 其中 p ~ i = p i / d i , i=1,…,M。然后利用 系统概率矩阵A的转置与该前向投影相乘得到反投影h: h = A T p ~ , h=[h1,…,hN]T,
7)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值 c:c=[c1,…,cN]T,其中cj=gj/hj,j=1,…,N,
8)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得 到中间图像 f ^ = [ f ^ 1 , . . . , f ^ N ] T , 其中 f ^ j = c j × f j , j=1,…,N,
9)对中间图像进行曲率流滤波:曲率流滤波能有效的去除噪声,但不损害图像 原有的特征信息,例如边缘、线、角等。它的滤波过程可以用如下的偏微分 方程表示:
f ( h , v ; t ) t =▿· ( f ( h , v ; t ) | f ( h , v ; t ) | ) - - - ( 7 )
以及边界条件: f ( h , v ; t ) N | f = Ω = 0 , 初始条件:f(h,v;0)=f0。这里表示梯度 算子;为了便于区别,我们用h,v表示图像的平和垂直两个方向的坐标; Ω为图像定义域;为图像边界,为边界法向量;t表示时间。引入时间 进化步长Δt(这里取Δt=0.0005较为合适)并将方程(7)右边项展开得:
f t + 1 = f t + Δt × ( f hh t ( f v t ) 2 - 2 f hv t f h t f v t + f vv t ( f h t ) 2 ( f h t ) 2 + ( f v t ) 2 ) - - - ( 8 )
其中ft≡f(h,v;t); f h t = f ( h , v ; t ) / h , f v t f ( h , v ; t ) / v 为t时刻函数f的水平和垂直 两个方向的一阶偏导数; f hh t 2 f ( h , v ; t ) / h 2 , f hv t f ( h , v ; t ) / h v f vv t 2 f ( h , v ; t ) / v 2 为二阶偏导数。
在实际情况中,我们通常考虑离散图像。为了运用公式(8)进行滤波,我们 需要引入合适的差分格式来替换公式(8)中的偏导数。为了方便叙述,我们可以 假设图像为K×K(K×K=N)的方阵且象素间距为单位长度1。记(fh t)j为函数 f(h,v;t)在点(hj,vj)(其中1≤hj≤K,1≤vj≤K)处的h方向的一阶偏导数,则它 可以用如下一阶差分表示:
( f h t ) j = f ( h j + 1 , v j ; t ) - f ( h j , v j ; t ) , j = 1 , . . . , N - - - ( 9 )
其中:hj=int[j/K]+1,vj=j-int[j/N]×N(int[x]表示取x的整数部分)以及 j=(hj-1)×K+vj。按同样的方法,我们还可以得到:
( f v t ) j = f ( h j , v j + 1 ; t ) - f ( h j , v j ; t ) , - - - ( 10 )
( f hh t ) j = f ( h j + 1 , v j ; t ) - 2 f ( h j , v j ; t ) + f ( h j - 1 , v j ; t ) , - - - ( 11 )
( f vv t ) j = f ( h j , v j + 1 ; t ) - 2 f ( h j , v j ; t ) + f ( h j , v j - 1 ; t ) , - - - ( 12 )
( f hv t ) j = ( f ( h j + 1 , v j + 1 ; t ) + f ( h j - 1 , v j - 1 ; t ) ) - ( f ( h j - 1 , v j + 1 ; t ) + f ( h j + 1 , v j - 1 ; t ) ) 4 - - - ( 13 )
同理令 f j t = f ( h j , v j ; t ) , j=1,…,N。这样我们可以将(8)改写为:
f j t + 1 = f j i + Δt × ( f hh t ) j ( f v t ) j 2 - 2 ( f hv t ) j ( f h t ) j ( f v t ) j + ( f vv t ) j ( f h t ) j 2 ( f h t ) j 2 + ( f v t ) j 2 , 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较为合理);设置初始条件: f j 0 = f ^ j , j=1,…,N,并置t=0;按公式(14)进行滤波,置t=t+1;重复公式(14)直至t=tmax 结束。将最终滤波结果保存在向量中, f ~ = [ f ~ 1 , . . . , f ~ N ] T .
10)将这个滤波图像作为初始图像f,返回到第6步,重复这个过程直到重建 后的图像收敛,可得到最终成像的图像。
高效检索全球专利

专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台,是提供专利分析,专利查询,专利检索等数据服务功能的知识产权数据服务商。

我们的产品包含105个国家的1.26亿组数据,免费查、免费专利分析。

申请试用

分析报告

专利汇分析报告产品可以对行业情报数据进行梳理分析,涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度,系统通过AI智能系统对图表进行解读,只需1分钟,一键生成行业专利分析报告。

申请试用

QQ群二维码
意见反馈