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

摘要

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

主权项

1.一种用于构建成像图像的迭代间曲率流滤波加权最小二乘正电子发射断层成像方法，其特征在于采用下列步骤：1)获取观测数据，记录其规模M，并将它保存在向量y＝[y₁，…，y_M]^T中；根据待重建图像的尺寸要求，针对正方形图像，设它的边长为K，确定初始图像的规模N＝K×K，给定它的初始灰度值大于1，并将其保存在向量f＝[f₁，…，f_N]^T中，设定迭代的总次数I，I＝3，2)根据观测数据y和初始图像f的大小，计算M×N的系统概率矩阵A，3)将观测数据y中的各个分量进行逐个与1比较：如果比1大则保留，比1小则设置为1，并将比较结果保存在向量d中，d＝[d₁，…，d_M]^T，作为观测数据的误差方差，4)将观测数据y中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的观测投影数据$\tilde{y}:\tilde{y}={[{\tilde{y}}_1,...,{\tilde{y}}_M]}^T,$其中${\tilde{y}}_i=y_i/d_i,$i＝1，…，M，5)将系统概率矩阵A的转置与修正的观测数据相乘得到向量g，$g=A^T\tilde{y},$g＝[g₁，…，g_N]^T，6)设置m＝0，7)将系统概率矩阵A和初始图像f相乘，得到前向投影p，p＝[p₁，…，p_M]^T，再将该前向投影p中的各个分量与误差方差d中的各个分量逐个对应相除得到修正的前向投影：$\tilde{p}={[{\tilde{p}}_1,...,{\tilde{p}}_M]}^T,$其中${\tilde{p}}_i=p_i/d_i,$i＝1，…，M，然后利用系统概率矩阵A的转置与该前向投影相乘得到反投影h：$h=A^T\tilde{p},$h＝[h₁，…，h_N]^T，8)将向量g的各个分量与反投影h中的各个分量逐个对应相除得到图像修正值c：c＝[c₁，…，c_N]^T，其中c_j＝g_j/h_j，j＝1，…，N，9)将该图像修正值c中的各个分量与初始图像f中的各个分量逐个对应相乘得到中间图像$\hat{f}:\hat{f}={[{\hat{f}}_1,...,{\hat{f}}_N]}^T,$其中${\hat{f}}_j=c_j\times f_j,$j＝1，…，N，10)设置t＝0，以及滤波总次数t_max，1≤t_max≤10，11)将中间图像重新记作一个新的临时图像向量${\hat{f}}^t:{\hat{f}}^t={[{\hat{f}}_l^t,...,{\hat{f}}_N^t]}^T,$12)用该临时图像计算如下导数向量：${\hat{f}}_h^t={[{\left(f_h^t\right)}_1,...,{\left(f_h^t\right)}_N]}^T,$${\hat{f}}_v^t={[{\left(f_v^t\right)}_1,...,{\left(f_v^t\right)}_N]}^T,$${\hat{f}}_{\mathrm{hh}}^t={[{\left(f_{\mathrm{hh}}^t\right)}_1,...,{\left(f_{\mathrm{hh}}^t\right)}_N]}^T,$${\hat{f}}_{\mathrm{vv}}^t={[{\left(f_{\mathrm{vv}}^t\right)}_1,...,{\left(f_{\mathrm{vv}}^t\right)}_N]}^T,$${\hat{f}}_{\mathrm{hv}}^t={[{\left(f_{\mathrm{hv}}^t\right)}_1,...,{\left(f_{\mathrm{hv}}^t\right)}_N]}^T,$其中${\left(f_h^t\right)}_j={\hat{f}}_{h_j\times K+v_j}^t-{\hat{f}}_{(h_j-1)\times K+v_j}^t,$${\left(f_v^t\right)}_j={\hat{f}}_{(h_j-1)\times K+v_j+1}^t-{\hat{f}}_{(h_j-1)\times K+v_j}^t,$${\left(f_{\mathrm{hh}}^t\right)}_j={\hat{f}}_{h_j\times K+v_j}^t-2\times {\hat{f}}_{(h_j-1)\times K+v_j}^t+{\hat{f}}_{(h_j-2)\times K+v_j}^t,$${\left(f_{\mathrm{vv}}^t\right)}_j={\hat{f}}_{(h_j-1)\times K+v_j+1}^t-2\times {\hat{f}}_{(h_j-1)\times K+v_j}^t+{\hat{f}}_{(h_j-1)\times K+v_j-1}^t,$${\left(f_{\mathrm{hv}}^t\right)}_j=\frac{1}{4}\times ({\hat{f}}_{h_j\times K+v_j+1}^t+{\hat{f}}_{(h_j-2)\times K+v_j-1}^t-{\hat{f}}_{(h_j-2)\times K+v_j+1}^t-{\hat{f}}_{h_j\times K+v_j-1}^t),$j＝1，…，N，这里的h_j和v_j用下面两式计算：h_j＝int[j/K]+1，v_j＝j-int[j/k]×K此处的int[j/K]表示取j除以K的商，其中f^t≡f(h，v；t)；$f_h^t\equiv \partial f(h,v;t)/\partial h,$$f_v^t\equiv \partial f(h,v;t)/\partial v$为t时刻函数f的水平和垂直两个方向的一阶偏导数；$f_{\mathrm{hh}}^t\equiv {\partial }^{2}f(h,v;t)/{\partial h}^2,$$f_{\mathrm{hv}}^t\equiv \partial f(h,v;t)/\partial h\partial v$和$f_{\mathrm{vv}}^t\equiv {\partial }^{2}f(h,v;t)/{\partial v}^2$为二阶偏导数，13)计算中间滤波图像${\tilde{f}}^t:{\tilde{f}}^t={[{\tilde{f}}_1^t,...,{\tilde{f}}_N^t]}^T,$其中${\tilde{f}}_j^t={\hat{f}}_j^t+\mathrm{\Delta t}\times \frac{{\left(f_{\mathrm{hh}}^t\right)}_j\times {\left((f_v^t{)}_j\right)}^2-2\times {\left(f_{\mathrm{hv}}^t\right)}_j\times {\left(f_h^t\right)}_j\times {\left(f_v^t\right)}_j+{\left(f_{\mathrm{vv}}^t\right)}_j\times {\left({\left(f_h^t\right)}_j\right)}^2}{\sqrt{{\left({\left(f_h^t\right)}_j\right)}^2+{\left({\left(f_v^t\right)}_j\right)}^2}},$j＝1，…，N14)将该得到的中间滤波图像赋值给临时图像，并将t的值增加1，15)如果t的值大于或等于滤波总次数t_max，则进行第16步，否则返回第11步，16)将临时图像保存为滤波图像，将这个滤波图像作为初始图像，m的值相应增加1，返回到第7步，重复这个过程直到m达到迭代的总次数I，可得到最终成像的图像。