CN110298900A - 一种基于各向异性基函数建立spect重构的方法 - Google Patents

Abstract

本发明公开了一种基于各向异性基函数建立SPECT重构的方法，包括下述步骤：S1：基于各向异性基函数建立SPECT重构的离散模型；S2：采用最大似然估计准则模拟SPECT重构，得到离散模型参数的估计值，建立SPECT重构的优化模型；S3：在SPECT重构的优化模型中引入两类不同的正则项，分别是HDCR和HFNR，建立描述SPECT重构的最小泛函；S4：采用不动点算法建立迭代模型，求解SPECT重构的优化模型。本发明建立SPECT重构的方法能够在不损坏图像质量的前提下减少SPECT检查中使用的放射试剂的剂量，并且能够有效消除由全变差正则化所导致的阶梯效应，同时能保持图像细节以及真实边界。

Description

一种基于各向异性基函数建立SPECT重构的方法

技术领域

本发明涉及医学影像技术领域，具体涉及一种基于各向异性基函数建立SPECT重构的方法。

背景技术

SPECT技术是通过估计投影数据的断层重构中的示踪剂分布来提供诊断信息。数十年来，它对人类的健康做出了巨大的贡献。但是，SPECT在图像质量和辐射剂量之间涉及到一个基本权衡。辐射接触的增加将会导致损害患者健康的可能性。因此，在SPECT的研究中，降低辐射剂量的非常有必要的。虽然，更低的辐射剂量能够通过降低示踪迹的活性来实现，但这将会不可避免的导致投影数据中更少的γ光子以及重构图像中更大的噪音，而这种噪声图像在临床上也许不那么有用。SPECT重构方法大致分为两个阶段，第一阶段是形成描述SPECT过程的离散模型。对此，已有的方法大多都是基于分片常数的离散模型和基于blob基函数的离散模型。基于分片常数的离散模型是用分片常数函数来逼近示踪剂的分布函数，进而获得离散模型；基于blob基函数的离散模型用blob基函数来逼近示踪剂分布函数。第二阶段是对离散模型进行求解，对此，已有的方法有EM(Expectation Maximization)方法和全变差正则化方法。EM方法是一种迭代优化策略，计算方法是每一次迭代都分两步，其中一步为期望步(E步)，另一个为极大步(M步)；全变差正则化方法是通过引入对解的偏导数先验约束条件，建立优化模型，进而求解优化问题，它能有效的抑制噪声，但是会使图像中产生阶梯伪影。

由于分片常数逼近精度低，传统的离散模型有一个瓶颈问题，任何图像处理算法均无法克服其引入的模型误差，并且离散模型不适用于低剂量的投影数据，低剂量的辐射会导致严重的图像噪音。运用Blob基函数方法后，得到的全离散系统矩阵会缺乏稀疏性。这是因为Blob基函数缺乏局部性。所以在使用Blob方法做迭代重构时，它需要更多的计算量来实现投影和反投影操作，并且Blob不能精确地表示分片连续图像。同时，更重要的是分片常数基函数与Blob基函数均是各向同性的，这使得它们无法表达医学影像所具的有各向异性的几何特征，例如示踪剂在心肌和脑部的分布函数就具有各向异性的特征。此外，EM收敛速度慢，对噪声敏感，不适用于低剂量辐射。全变差正则化为了压制环境噪音，要求对所有的灰度值变化作同等的惩罚，这种惩罚方法同时也抑制、模糊了图像中不同解剖结构之间的真实边界。

发明内容

为了克服现有技术的缺陷，本发明提供一种基于各向异性基函数建立SPECT重构的方法，能够在不损坏图像质量的前提下减少SPECT检查中使用的放射试剂的剂量，并且能够有效消除由全变差正则化所导致的阶梯效应，同时能保持图像细节以及真实边界。

为了达到上述目的，本发明采用以下技术方案：

1、一种基于各向异性基函数建立SPECT重构的方法，包括下述步骤：

S1：基于各向异性基函数建立SPECT重构的离散模型：

S11：SPECT重构的离散模型采用连续积分方程表示，所述连续积分方程根据示踪剂分布活性函数和积分核函数构建，所述积分核函数包括衰减校正，康普顿散射校正和系统空间分辨率校正的影响；

S12：采用各向异性基函数逼近连续积分方程中的示踪剂分布活性函数，构建离散积分模型；

S13：构造分片多项式子空间X的基底函数和SPECT系统矩阵；

S14：将离散积分模型表示为线性方程；

S2：采用最大似然估计准则模拟SPECT重构，得到离散积分模型参数的估计值，建立SPECT重构的优化模型；

S3：在SPECT重构的优化模型中引入两类不同的正则项，分别是HDCR和HFNR，建立描述SPECT重构的最小泛函；

S4：采用不动点算法建立迭代模型，求解SPECT重构的优化模型。

作为优选的技术方案，所述各向异性基函数是在四组不同的网格上分别建立高阶分片多项式，然后组合得到的一组基函数。

作为优选的技术方案，步骤S11中所述的SPECT重构的离散模型采用连续积分方程表示，具体公式为：

其中，f(p)是示踪剂分布活性函数，是投影分布函数，h(u,θ,p)是积分核函数，方形图像域 表示探测器表面的闭区间，Θ＝[0,2π]为旋转角范围，u表示探测器上的点，θ表示闪烁照相机逆时针旋转的角度，p:＝(x,y)表示点源，x和y为x0y平面上的坐标；

其中：

h(u,θ,p)＝R(u′,θ,p)K(u,u′,v)；

R(u′,θ,p)：＝exp{-∫_{l(u′，θ，p)}μ(p′)ds(p′)}；

其中，u′＝xcosθ+ysinθ,v＝xsinθ+ycosθ,l(u′,θ,p)表示点源p和探测板之间沿着l(u′,θ)的衰减路径，表示将xy坐标系上的点p(x,y)投影到u′v坐标系上所做的变换；l(u′,θ)表示经过点(u′cosθ，u′sinθ)，并垂直于向量(cosθ，sinθ)的直线，刻画xy平面内的任意一条直线，σ_s是高斯型系统点扩展函数的与距离相关的标准偏差。

作为优选的技术方案，步骤S12中所述采用各向异性基函数逼近连续积分方程中的示踪剂分布活性函数，构建离散积分模型，具体步骤为：

所述各向异性基函数的表达式如下所述：

其中，表示探测板在投影角度记录的gamma光子数，为在摄影角度的连续投影分布函数；

设置用于解连续积分方程，X表示有限维度的分片多项式空间L^∞(Ω)的子空间，构建离散积分模型为：

函数值h_i(p)表示第i个探测元在点源p处射出的辐射能量的响应，m＝qr表示离散探测元的数目，q为探测器上单元的个数，r为投影视角的个数。

作为优选的技术方案，步骤S13中所述构造分片多项式子空间X的基底函数和SPECT系统矩阵，具体步骤为：

分片多项式子空间X的基底的具体获取方式为：

先取一个正方形区域Ω的网格，网格包含了等面积的q×q个正方形格子，面积都为h²，从左下角的格子开始，逆时针索引网格的矩阵，索引序列中第j个元素对应着矩阵中第(j₁,j₂)个格子，j₁,j₂＝1,2,…,q,其中j＝j₁+(j₂-1)q；

给定第j个格子,选择一组张量积多项式作为基函数：其中，(x_j,y_j)表示网格w_j的左下角的坐标，有k＝k₁+1+k₂n，k₁，k₂＝0，1,…，n-1，n为一维网格上基函数的个数；

将所得到的q×q个正方形格子及定义的基函数分别以正方形区域Ω的中心逆时针旋转30°和60°,再得到另外两组基函数和三组基函数合并后作为解空间X的基底，即

构造的SPECT系统矩阵为：

其中，A_k＝[a_iljk:i＝1,2,3；l＝1,2…,m；j＝1,2,…,q²]为一个m×q²矩阵，分块矩阵A_k是系统矩阵的第k列向量，分块矩阵A_k中的每一个元素

作为优选的技术方案，步骤S14所述将离散积分模型表示为线性方程，具体步骤如下所述：

定义两个列向量，有f的分块矩阵把将离散积分模型表示为矩阵向量乘形式：

通过矩阵向量乘，求出f_i,j,k,并且用线性组合求解离散积分模型。

作为优选的技术方案，步骤S2所述的建立SPECT重构的优化模型，具体步骤如下所述：

S21：将离散积分模型采用算子方程式表示为：根据泊松分布的性质，得到服从泊松分布的投影数据的对数似然函数；

S22：对f的似然函数和惩罚因子进行求和，然后再极大化和函数，并获得模型参数的估计值，SPECT重构的优化模型用方程式表示为：

f_★＝argmin{<Af，1>-〈ln(Af)，g>+λR(f)}.

其中，数据拟合项F：＝<Af，1>-<ln(Af)，g>是KL散度，正则项R(f)是Rⁿ上的正凸函数，λ＞0用来平衡F与R(f)之间的权重。

作为优选的技术方案，步骤S3所述的建立描述SPECT重构的最小泛函，具体步骤如下所述：

引入HDCR，惩罚高阶方向导数不连续性，

其中

将三角形不等式应用于Rⁱ _I(f)所有绝对值，并且计算u幂的积分，引入一个与R_I(f)相等的替代正则项，

R_surr(f)分解为的形式，其中Φ是范数，B是一个线性算子；

引入HFNR，惩罚在Ω域中重构图像的n阶微分，

带入的线性表达式为：

对于每个j∈{1,2,...,q²},定义索引集S_j:＝{j,j+1,j+p,j+p+1},并且X的子空间为分别以三种不同的基函数张成的子空间，{ψ_i,j,k：i＝1,2,3,j＝1,2,...,q²,k＝1,2,...,3n²}为W_i,j的一组正交基，即有其中Y_j:＝span{v_j,k:k＝1,2,...n²}为粗尺度空间，有张量积k＝(k₁-1)n+k₂。

作为优选的技术方案，步骤S4所述的采用不动点算法建立迭代模型，具体步骤如下所述：

其中，Prox表示的迫近算子，b^k是定义在正则化转换域的对偶迭代，f^k为k次迭代时f的值，B^T(2b^k+1-b^k+1)表示图像域中的噪音；Q:＝μ^-1I_q,P:＝β^-1S^-1，μ和β为正参数，I_q为q维单位矩阵，S是d×d对角正定预处理矩阵，迭代相对误差定义为：

tol＝||f^k-f^k+1||₂/||f^k+1||₂。

本发明与现有技术相比，具有如下优点和有益效果：

(1)本发明采用基于各向异性基函数建立SPECT重构的离散模型，比起传统的分片常数逼近方法，在示踪剂剂量少的情况下，能够提高重构图像的精度，改善图像的质量，减少不可约误差。

(2)本发明引入两类正则化方法，代替了全变差正则化，能够有效消除由全变差正则化所导致的阶梯效应，同时能保持图像细节以及真实边界。

(3)本发明在优化图像时运用了不动点算法，基于不动点刻画，发展出高效的数值求解算法，方便优化问题的数值求解。

附图说明

图1为本实施例的基于各向异性基函数建立SPECT重构方法流程图；

图2为本实施例SPECT的成像系统的固定坐标系和旋转坐标系示意图；

图3为本实施例构造基底选取正方形区域Ω的网格示意图；

图4为本实施例三个不同角度的正方形区域Ω示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

在本实施例中，提供一种基于各向异性基函数建立SPECT重构的方法，包含三个部分：基于各向异性基函数建立SPECT重构的离散模型，以及与之对应的正则化方法和不动点迭代格式。所述各向异性基函数是在四组不同的网格上分别建立高阶分片多项式，然后组合起来得到的一组基函数。采用向异性基函数能够有效的表示医学影像中的几何特征，同时，能得到具有更高精度的离散模型，并且能够满足在低剂量辐射的条件下提高图像的质量，减少图像噪音。基于各向异性基函数的正则化方法分为两类，第一类是惩罚高阶方向导数不连续性的正则化方法，表示为HDCR；第二类是惩罚图像高频噪声的正则化方法表示为HFNR。基于各向异性基函数建立SPECT重构的离散模型有如下几个优点：首先，由物理定律直接推导而来，因而更复合实际的物理机制；其次，采用各向异性基函数允许采用更高阶的精度对其进行离散，从而能极大的减少模型误差；最后，各向异性基函数有利于自适应逼近算法的实现，这在医学图像处理中尤为重要。所述HDCR和HFNR能够惩罚医学影像的各阶方向导数不连续性，比如：在网格边界的跳跃，或网格边界的法向导数的跳跃。HDCR是惩罚两个相邻网格在其共同边界上高阶法向导数之间的差异；HFNR能够有效消除由全变差正则化所导致的阶梯效应，同时能保持图像细节以及真实边界。在HFNR正则化方法中，通过建立多点的高阶差分格式从而惩罚f的n阶微分，从而实现对高频噪声的抑制。优化图像时运用了不动点算法，基于不动点刻画，本实施例提供一种高效率的数值求解算法，方便优化问题的数值求解。

如图1所示，在本实施例中，基于各向异性基函数建立SPECT重构的方法具体步骤如下所述：

S1：基于各向异性基函数建立SPECT重构的离散模型；

不同模态的SPECT成像过程均可刻画成如下一般形式的连续模型g＝Af，其中g表示探测到的降质投影数据，f表示待重构的示踪剂分布函数，线性算子A则描述平行射线束SPECT成像系统的物理降质过程，A包含衰减校正，康普顿散射，空间分辨率或者偶然重合等降质因素或其复合。具体地,本实施例将Radon变换、光电效应、康普顿散射以及空间分辨率等物理降质因素，描述为积分算子，从而建立起平行射线束SPECT系统的积分方程成像模型,具体地积分方程为:

其中，f(p)是示踪剂分布活性函数，是投影分布函数，h(u,θ,p)是积分核函数，它包括衰减校正，康普顿散射校正和系统空间分辨率校正的影响。方形图像域 表示探测器表面的闭区间，Θ＝[0,2π]为旋转角范围，u表示探测器上的点，θ表示闪烁照相机逆时针旋转的角度，p:＝(x,y)表示点源，x和y为x0y平面上的坐标。SPECT重构的核心就是成像模型的完全离散，由在成像系统的离散点取样得到。

如图2所示，先介绍SPECT的两个成像系统，一个固定坐标系xyz和一个旋转坐标系u′vw，闪烁照相机绕着z轴旋转，旋转坐标系坐落在闪烁照相机的表面。当旋转角为0的时候，u′坐标轴与x坐标轴平行，w与z轴相一致，且与旋转角度无关。坐标轴z和w是无关紧要的，这是因为平行孔准直器把gamma射线束的线积分均限定在xy平面内。

其中：

h(u,θ,p)＝R(u′,θ,p)K(u,u′,v)；

其中，u′＝xcosθ+ysinθ,v＝xsinθ+ycosθ,l(u′,θ,p)表示点源p和探测板之间沿着l(u′,θ)的衰减路径，表示将xy坐标系上的点p(x,y)投影到u′v坐标系上所做的变换。l(u′,θ)表示经过点(u′cosθ，u′sinθ)，并垂直于向量(cosθ，sinθ)的直线，这可以刻画xy平面内的任意一条直线，σ_s是高斯型系统点扩展函数的与距离相关的标准偏差；

采用了各向异性基函数来逼近示踪剂分布活性函数f(p),X表示有限维度的分片多项式空间L^∞(Ω)的子空间。在本实施例设置一个用于解方程(1)，各向异性基函数的表达式如下所述：

其中，表示探测板在投影角度记录的gamma光子数，为在摄影角度的连续投影分布函数；

设置用于解连续积分方程，X表示有限维度的分片多项式空间L^∞(Ω)的子空间，构建离散积分模型为：

其中，函数值h_i(p)表示第i个探测元在点源p处射出的辐射能量的响应，m＝qr表示离散探测元的数目，q为探测器上单元的个数，r为投影视角的个数。

一旦X的基底选定了，便成为了一个完整的离散线性系统。在本实施例中，采用各向异性基函数作为X的基底。因为g函数是通过采样得到的离散数据，所以要先得到离散的f，然后再通过和基函数的线性组合得到连续函数f，在本实施例中，f表示连续函数，f是离散形式。

在本实施例中，分片多项式子空间X的基底的具体获取方式为：

如图3所示，先取一个正方形区域Ω的网格，网格包含了等面积的q×q个正方形格子，面积都为h²，从左下角的格子开始，逆时针索引网格的矩阵，索引序列中第j个元素对应着矩阵中第(j₁,j₂)个格子，j₁,j₂＝1,2,…,q,其中j＝j₁+(j₂-1)q；给定第j个格子,选择一组张量积多项式作为基函数：其中，(x_j,y_j)表示网格w_j的左下角的坐标，有k＝k₁+1+k₂n，k₁、k₂＝0，1,…，n-1，n为一维网格上基函数的个数。

将所得到的q×q个正方形格子及由其定义的基函数分别以Ω的中心逆时针旋转30°和60°,如图4所示的两组虚线方形区域,再得到另外两组基函数和一起合并为三组基函数，来作为解空间X的基底，即

生成的SPECT系统矩阵为：

其中，，A_k＝[a_iljk:i＝1,2,3；l＝1,2…,m；j＝1,2,…,q²]为一个m×q²矩阵，分块矩阵A_k是系统矩阵的第k列向量，分块矩阵A_k中的每一个元素

定义两个列向量，有定义的f的分块矩阵把离散方程写为矩阵向量乘形式：通过矩阵向量乘，求出f_i,j,k,并且用线性组合来解离散方程。

S2：采用最大似然估计准则模拟SPECT重构，得到离散模型参数的估计值，建立SPECT重构的优化模型；对于特定的统计模型和给定的一组观测数据，最大似然(ML)估计准则通过极大化似然函数来获得模拟参数的估计值，本实施例用最大似然估计准则来模拟图像重构问题。

具体步骤为：

对于离散的随机变量来说，ML估计实则是在极大化观测数据在模型参数取定的前提条件下的条件概率。将ML估计准则应用到ECT图像重构问题，ML估计值即为重构图像。

由于粒子的计数过程是一个泊松过程，将离散模型采用算子方程式表示为：根据泊松分布的性质，得到服从泊松分布的投影数据的对数似然函数；

最后，ML估计准则通过对f的似然函数和惩罚因子进行求和，然后再极大化这个和函数，并获得模型参数的估计值，SPECT重构的优化模型用方程式表示为：

f_*＝argmin{<Af,1>-<ln(Af),g>+λR(f)}

其中，数据拟合项F：＝<Af,1>-<ln(Af)，g>是KL散度，正则项R(f)是Rⁿ上的正凸函数，λ＞0，用来平衡F与R(f)之间的权重。

S3：在SPECT重构的优化模型中引入两类不同的正则项，分别是HDCR和HFNR，建立描述SPECT重构的最小泛函；

引进两类不同的正则项，第一类是惩罚高阶方向导数不连续性的正则化方法(HDCR)；第二类是惩罚图像高频噪声的正则化方法(HFNR)。具体过程如下：

①HDCR

HDCR是惩罚高阶方向导数不连续性，具体是分别惩罚三组方形网格中两个相邻格子在共同边界上解f局限的差异，

其中

离散f得到的分量函数

将三角形不等式应用于Rⁱ _I(f)所有绝对值，并且计算u幂的积分，接着，引进了一个与R_I(f)相等的替代正则项：

这里，替代正则项R_surr(f)可以分解为的形式，其中Φ是范数，B是一个线性算子，以便于后期算法的计算。

②HFNR

HFNR是惩罚在Ω域中重构图像的n阶微分，这能够促进在重建系统中解的稀疏性。

带入的线性表达式为：

对于每个j∈{1,2,...,q²},定义索引集S_j:＝{j,j+1,j+p,j+p+1},并且X的子空间为分别以三种不同的基函数张成的子空间，{ψ_i,j,k：i＝1,2,3,j＝1,2,...,q²,k＝1,2,...,3n²}为W_i,j的一组正交基，即有其中Y_j:＝span{v_j,k:k＝1,2,...n²}为粗尺度空间，有张量积k＝(k₁-1)n+k₂。

S4：采用不动点算法建立迭代模型，求解SPECT重构的优化模型。

用HDCR和HFNR建立描述重构图像最小泛函后，优化问题的求解需要有效的迭代方法，本实施例可以用不动点迫近算法解决不光滑的优化问题，迭代算法如下：

Prox表示φ^*迫近算子，b^k是定义在正则化转换域的对偶迭代，f^k为k次迭代时f的值，B^T(2b^k+1-b^k+1)表示图像域中的噪音；Q:＝μ^-1I_q,P:＝β^-1S^-1，μ和β为正参数，I_q为q维单位矩阵，S是d×d对角正定预处理矩阵，迭代相对误差定义为：

tol＝||f^k-f^k+1||₂/||f^k+1||₂。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (9)

1.一种基于各向异性基函数建立SPECT重构的方法，其特征在于，包括下述步骤：

S1：基于各向异性基函数建立SPECT重构的离散模型：

S11：SPECT重构的离散模型采用连续积分方程表示，所述连续积分方程根据示踪剂分布活性函数和积分核函数构建，所述积分核函数包括衰减校正，康普顿散射校正和系统空间分辨率校正的影响；

S12：采用各向异性基函数逼近连续积分方程中的示踪剂分布活性函数，构建离散积分模型；

S13：构造分片多项式子空间X的基底函数和SPECT系统矩阵；

S14：将离散积分模型表示为线性方程；

S2：采用最大似然估计准则模拟SPECT重构，得到离散积分模型参数的估计值，建立SPECT重构的优化模型；

S3：在SPECT重构的优化模型中引入两类不同的正则项，分别是HDCR和HFNR，建立描述SPECT重构的最小泛函；

S4：采用不动点算法建立迭代模型，求解SPECT重构的优化模型。

2.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，所述各向异性基函数是在四组不同的网格上分别建立高阶分片多项式，然后组合得到的一组基函数。

3.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S11中所述的SPECT重构的离散模型采用连续积分方程表示，具体公式为：

其中，f(p)是示踪剂分布活性函数，是投影分布函数，h(u,θ,p)是积分核函数，方形图像域 表示探测器表面的闭区间，Θ＝[0,2π]为旋转角范围，u表示探测器上的点，θ表示闪烁照相机逆时针旋转的角度，p:＝(x,y)表示点源，x和y为x0y平面上的坐标；

其中：

h(u,θ,p)＝R(u′,θ,p)K(u,u′,v)；

R(u′,θ,p)∶＝exp{-∫_{l(u′，θ，p)}μ(p′)ds(p′)}；

其中，u′＝xcosθ+ysinθ,v＝xsinθ+ycosθ,l(u′,θ,p)表示点源p和探测板之间沿着l(u′,θ)的衰减路径，表示将xy坐标系上的点p(x,y)投影到u′v坐标系上所做的变换；l(u′,θ)表示经过点(u′cosθ，u′sinθ)，并垂直于向量(cosθ，sinθ)的直线，刻画xy平面内的任意一条直线，σ_s是高斯型系统点扩展函数的与距离相关的标准偏差。

4.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S12中所述采用各向异性基函数逼近连续积分方程中的示踪剂分布活性函数，构建离散积分模型，具体步骤为：

所述各向异性基函数的表达式如下所述：

其中，表示探测板在投影角度记录的gamma光子数，为在摄影角度的连续投影分布函数；

设置用于解连续积分方程，X表示有限维度的分片多项式空间L^∞(Ω)的子空间，构建离散积分模型为：

函数值h_i(p)表示第i个探测元在点源p处射出的辐射能量的响应，m＝qr表示离散探测元的数目，q为探测器上单元的个数，r为投影视角的个数。

5.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S13中所述构造分片多项式子空间X的基底函数和SPECT系统矩阵，具体步骤为：

分片多项式子空间X的基底的具体获取方式为：

先取一个正方形区域Ω的网格，网格包含了等面积的q×q个正方形格子，面积都为h²，从左下角的格子开始，逆时针索引网格的矩阵，索引序列中第j个元素对应着矩阵中第(j₁,j₂)个格子，j₁,j₂＝1,2,…,q,其中j＝j₁+(j₂-1)q；

给定第j个格子,选择一组张量积多项式作为基函数：其中，(x_j,y_j)表示网格w_j的左下角的坐标，有k＝k₁+1+k₂n，k₁，k₂＝0，1,…，n-1，n为一维网格上基函数的个数；

将所得到的q×q个正方形格子及定义的基函数分别以正方形区域Ω的中心逆时针旋转30°和60°,再得到另外两组基函数和三组基函数合并后作为解空间X的基底，即

构造的SPECT系统矩阵为：

其中，A_k＝[a_iljk:i＝1,2,3；l＝1,2…,m；j＝1,2,…,q²]为一个m×q²矩阵，分块矩阵A_k是系统矩阵的第k列向量，分块矩阵A_k中的每一个元素

6.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S14所述将离散积分模型表示为线性方程，具体步骤如下所述：

定义两个列向量，有f的分块矩阵把将离散积分模型表示为矩阵向量乘形式：

通过矩阵向量乘，求出f_i,j,k,并且用线性组合求解离散积分模型。

7.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S2所述的建立SPECT重构的优化模型，具体步骤如下所述：

S21：将离散积分模型采用算子方程式表示为：根据泊松分布的性质，得到服从泊松分布的投影数据的对数似然函数；

S22：对f的似然函数和惩罚因子进行求和，然后再极大化和函数，并获得模型参数的估计值，SPECT重构的优化模型用方程式表示为：

f_★＝arg min{<Af，1>-<ln(Af)，g>+λR(f)}.

其中，数据拟合项F∶＝<Af,1>-<ln(Af)，g>是KL散度，正则项R(f)是Rⁿ上的正凸函数，λ>0用来平衡F与R(f)之间的权重。

8.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S3所述的建立描述SPECT重构的最小泛函，具体步骤如下所述：

引入HDCR，惩罚高阶方向导数不连续性，

其中

将三角形不等式应用于Rⁱ _I(f)所有绝对值，并且计算u幂的积分，引入一个与R_I(f)相等的替代正则项，

R_surr(f)分解为的形式，其中Φ是范数，B是一个线性算子；

引入HFNR，惩罚在Ω域中重构图像的n阶微分，

带入的线性表达式为：

对于每个j∈{1,2,...,q²},定义索引集S_j:＝{j,j+1,j+p,j+p+1},并且X的子空间为分别以三种不同的基函数张成的子空间，{ψ_i,j,k：i＝1,2,3,j＝1,2,...,q²,k＝1,2,...,3n²}为W_i,j的一组正交基，即有其中Y_j:＝span{v_j,k:k＝1,2,...n²}为粗尺度空间，有张量积k＝(k₁-1)n+k₂。

9.根据权利要求1所述的基于各向异性基函数建立SPECT重构的方法，其特征在于，步骤S4所述的采用不动点算法建立迭代模型，具体步骤如下所述：

其中，Prox表示的迫近算子，b^k是定义在正则化转换域的对偶迭代，f^k为k次迭代时f的值，B^T(2b^k+1-b^k+1)表示图像域中的噪音；Q:＝μ^-1I_q,P:＝β^-1S^-1，μ和β为正参数，I_q为q维单位矩阵，S是d×d对角正定预处理矩阵，迭代相对误差定义为：

tol＝||f^k-f^k+1||₂/||f^k+1||₂。