基于深度学习方法实现SPECT放射性核素定量测量（三）-图像重建及衰减校正

（1）图像重建

使用MLEM迭代重建算法进行图像重建，MLEM算法原理如图1所示。

图像重建可以借助重建软件进行，目前可以用于SPECT重建的软件有STIR、Castor、OMEGA、QSPECT等，这些软件大部分是在linux系统下使用的。可以使用QSPECT对SPECT图像进行图像重建，QSPECT是专用于SPECT图像重建的软件，该软件最为轻便，而且其3.0版本可以在windows系统下使用。
借助QSPECT可以方便的将projection data转换为sinogram data，可以很简单的进行迭代计算，该软件使用简单，一看就会。
QSPECT下载网址：http://www.bme.teiath.gr/ni/EnglishVersion/Downloads_files/Downloads.html

（2）衰减校正

由于光子在体模中存在衰减，因此会导致探测器探测到的计数降低，重建图像计数降低。因此有必要对重建图像进行衰减校正以提高图像的精度和准确度。使用chang氏衰减校正方法对重建后图像进行衰减校正，校正公式为：

c ( x , y ) = [ 1 m ∑ i = 1 m e x p ( − μ l ( x , y , θ i ) ) ] − 1 c(x,y)=[\frac{1}{m}\sum_{i=1}^{m}exp(-\mu l(x,y,\theta_i))]^{-1} c(x,y)=[m1∑i=1mexp(−μl(x,y,θi))]−1
f c ( x , y ) = f ( x , y ) c ( x , y ) f_c(x,y)=f(x,y)c(x,y) fc(x,y)=f(x,y)c(x,y)

式中 f ( x , y ) f(x,y) f(x,y)为重建图像， f c ( x , y ) f_c(x,y) fc(x,y)为衰减校正后图像， c ( x , y ) c(x,y) c(x,y)为校正因子分布， l ( x , y , θ i ) l(x,y,\theta_i) l(x,y,θi)为点 ( x , y ) (x,y) (x,y)在 θ i \theta_i θi方向到体模边界的距离， μ \mu μ为线性衰减系数。使用这种方法进行衰减校正需要先获取体模横断面的外轮廓线以求出 l l l，并且假定衰减系数均匀分布。对于衰减系数非均匀分布的体模可以采用下式计算校正因子：

c ( x , y ) = [ 1 m ∑ i = 1 m e x p ( − ∑ j = 1 n μ j l j ( x , y , θ i ) ) ] − 1 c(x,y)=[\frac{1}{m}\sum_{i=1}^m exp(-\sum_{j=1}^n\mu_j l_j(x,y,\theta_i))]^{-1} c(x,y)=[m1∑i=1mexp(−∑j=1nμjlj(x,y,θi))]−1
(该式是根据chang式衰减原理自行推导得到的，不一定正确！如有错误，敬请指出)

衰减校正部分主要是根据CT图像计算衰减校正因子c，以下给出计算衰减系数非均匀分布情况下计算衰减校正因子c的代码。由于探测过程使用的是平行孔准直器，因此可以考虑将3D图像分割成若干个2D图像简化计算过程，对每个2D的图像计算衰减校正因子。

在计算校准因子c的过程中， l ( x , y , θ i ) l(x,y,\theta_i) l(x,y,θi)的计算较为复杂，在计算中考虑借鉴迭代重建算法中投影矩阵的计算方法，在投影矩阵计算方法中有一种是使用快速网格遍历法求解权重因子，具体可以参考书籍《医学断层图像仿真实验》中第13次实验。

l ( x , y , θ i ) l(x,y,\theta_i) l(x,y,θi)求解的matlab代码如下：

function Par_exp = distance(x0, y0, theta, n, delta, u)
%非均匀衰减
%计算(x0,y0)沿theta方向每一段u衰减的距离并与l相乘
%返回par_exp=u_j*l_j
%x0, y0：待计算点坐标
%theta：投影角度，角度制
%n：像素大小
%delta：每次步进距离，default:1
%l：射线长度
%u：衰减系数(矩阵)
Theta = theta*pi/180;     %角度转换成弧度
m = tan(Theta);
b = y0 - m*x0;            %射线方程：y=mx+b
X = x0;
Y = y0;
d1=0;
Par_exp=0;
%分类：
% (1)0<m<1( 0~45度；180~225度)
% (2)1<m(45~90度；225~270度)
% (3)-1<m<0(135~180度；315~360度)
% (4)m<-1(90~135度；270~315)
% 当前网格(X,Y)function judge=judgement(x, y, N)if x>N || y>N || x<1 || y<1judge=0;elsejudge=1;endend%%%%%%%%%%%%%%% 0度 %%%%%%%%%%%%%%%
if Theta==0for i = 1:nif X+1<=128Par_exp=Par_exp+delta*u(X,Y);X=X+1;elsebreak;endend
end
%%%%%%%%%%%%%%% 90度 %%%%%%%%%%%%%%%
if Theta==pi/2for i = 1:nif Y+1<=128Par_exp=Par_exp+delta*u(X,Y);Y=Y+1;elsebreak;endend
end
%%%%%%%%%%%%%%% 180度 %%%%%%%%%%%%%%%
if Theta==pifor i = 1:nif X-1>=1Par_exp=Par_exp+delta*u(X,Y);X=X-1;elsebreak;endend
end
%%%%%%%%%%%%%%% 270度 %%%%%%%%%%%%%%%
if Theta==3*pi/2for i = 1:nif Y-1>=1Par_exp=Par_exp+delta*u(X,Y);Y=Y-1;elsebreak;endend
end %%%%%%%%%%%%%%%0<m<1( 0~45度) (0, pi/4]%%%%%%%%%%%%%%%
if Theta>0 && Theta<= pi/3for i = 1:n*nd2 = d1+m;if d1>=0 && d2>deltaif judgement(X+1,Y+1,n)==0break;elseBD = delta*sqrt(m^2+1);Par_exp=Par_exp+BD*u(X,Y);d1=d2-delta;X=X+1;Y=Y+1;endendif d1>=0 && d2==deltaif judgement(X+1,Y+1, n)==0break;elseBE = delta*sqrt(m^2+1);Par_exp=Par_exp+BE*u(X,Y);d1 = 0;X=X+1;Y=Y+1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X+1, Y, n)==0break;elseBF = delta*sqrt(m^2+1);Par_exp=Par_exp+BF*u(X,Y);d1 = d2;X=X+1;Y=Y;endendif d1<0if judgement(X+1,Y, n)==0break;elseMT = d2*sqrt(m^2+1)/m;Par_exp=Par_exp+MT*u(X,Y);d1 = d2;X=X+1;Y=Y;endendend
end
%%%%%%%%%%%%%%%0<m<1(180~225度) (pi, 5pi/4]%%%%%%%%%%%%%%%
if Theta>pi && Theta<= 5*pi/4for i = 1:n*nd2 = d1+m;if d1>=0 && d2>deltaif judgement(X-1,Y-1,n)==0break;elseBD = delta*sqrt(m^2+1);Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X-1;Y=Y-1;endendif d1>=0 && d2==deltaif judgement(X-1,Y-1, n)==0break;elseBE = delta*sqrt(m^2+1);Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X-1;Y=Y-1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X-1, Y, n)==0break;elseBF = delta*sqrt(m^2+1);Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X-1;Y=Y;endendif d1<0if judgement(X-1,Y, n)==0break;elseMT = d2*sqrt(m^2+1)/m;Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X-1;Y=Y;endendend
end            %%%%%%%%%%%%%%%1<m(45~90度) [pi/4, pi/2)%%%%%%%%%%%%%%%
if Theta>=pi/4 && Theta< pi/2for i = 1:n*nd2 = d1+delta/m;if d1>=0 && d2>deltaif judgement(X+1,Y+1, n)==0break;elseBD = delta*sqrt(m^2+1)/m;Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X+1;Y=Y+1;endendif d1>=0 && d2==deltaif judgement(X+1,Y+1, n)==0break;elseBE = delta*sqrt(m^2+1)/m;Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X+1;Y=Y+1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X, Y+1, n)==0break;elseBF = delta*sqrt(m^2+1)/m;Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X;Y=Y+1;endendif d1<0if judgement(X,Y+1, n)==0break;elseMT = d2/sqrt(m^2+1);Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X;Y=Y+1;endendend
end
%%%%%%%%%%%%%%%1<m(225~270度) [5pi/4, 3pi/2)%%%%%%%%%%%%%%%
if Theta>=5*pi/4 && Theta<3*pi/2for i = 1:n*nd2 = d1+delta/m;if d1>=0 && d2>deltaif judgement(X-1,Y-1, n)==0break;elseBD = delta*sqrt(m^2+1)/m;Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X-1;Y=Y-1;endendif d1>=0 && d2==deltaif judgement(X-1,Y-1, n)==0break;elseBE = delta*sqrt(m^2+1)/m;Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X-1;Y=Y-1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X, Y-1, n)==0break;elseBF = delta*sqrt(m^2+1)/m;Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X;Y=Y-1;endendif d1<0if judgement(X,Y-1, n)==0break;elseMT = d2/sqrt(m^2+1);Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X;Y=Y-1;endendend
end        %%%%%%%%%%%%%%%-1<m<0(135~180度) [3pi/4, pi)%%%%%%%%%%%%%%%
if Theta>=3*pi/4 && Theta<pifor i = 1:n*nd2 = d1-delta*m;if d1>=0 && d2>deltaif judgement(X-1,Y+1, n)==0break;elseBD = delta*sqrt(m^2+1);Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X-1;Y=Y+1;endendif d1>=0 && d2==deltaif judgement(X-1,Y+1, n)==0break;elseBE = delta*sqrt(m^2+1);Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X-1;Y=Y+1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X-1, Y, n)==0break;elseBF = delta*sqrt(m^2+1);Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X-1;Y=Y;endendif d1<0if judgement(X-1,Y, n)==0break;elseMT = d2*sqrt(m^2+1)/(-m);Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X-1;Y=Y;endendend
end
%%%%%%%%%%%%%%%-1<m<0(315~360度) [7pi/4, 2pi)%%%%%%%%%%%%%%%
if Theta>=7*pi/4 && Theta<2*pifor i = 1:n*nd2 = d1-delta*m;if d1>=0 && d2>deltaif judgement(X+1,Y-1, n)==0break;elseBD = delta*sqrt(m^2+1);Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X+1;Y=Y-1;endendif d1>=0 && d2==deltaif judgement(X+1,Y-1, n)==0break;elseBE = delta*sqrt(m^2+1);Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X+1;Y=Y-1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X+1, Y, n)==0break;elseBF = delta*sqrt(m^2+1);Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X+1;Y=Y;endendif d1<0if judgement(X+1,Y, n)==0break;elseMT = d2*sqrt(m^2+1)/(-m);Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X+1;Y=Y;endendend
end%%%%%%%%%%%%%%%m<-1(90~135度) (pi/2, 3pi/4]%%%%%%%%%%%%%%%
if Theta>pi/2 && Theta<=3*pi/4for i = 1:n*nd2 = d1-delta/m;if d1>=0 && d2>deltaif judgement(X-1,Y+1, n)==0break;elseBD = delta*sqrt(m^2+1)/(-m);Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X-1;Y=Y+1;endendif d1>=0 && d2==deltaif judgement(X-1,Y+1, n)==0break;elseBE = delta*sqrt(m^2+1)/(-m);Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X-1;Y=Y+1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X, Y+1, n)==0break;elseBF = delta*sqrt(m^2+1)/(-m);Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X;Y=Y+1;endendif d1<0if judgement(X,Y+1, n)==0break;elseMT = d2*sqrt(m^2+1);Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X;Y=Y+1;endendend
end
%%%%%%%%%%%%%%%m<-1(270~315度) (3pi/2, 7pi/4]%%%%%%%%%%%%%%%
if Theta>3*pi/2 && Theta<=7*pi/4for i = 1:n*nd2 = d1-delta/m;if d1>=0 && d2>deltaif judgement(X+1,Y-1, n)==0break;elseBD = delta*sqrt(m^2+1)/(-m);Par_exp=Par_exp+BD*u(X,Y); d1=d2-delta;X=X+1;Y=Y-1;endendif d1>=0 && d2==deltaif judgement(X+1,Y-1, n)==0break;elseBE = delta*sqrt(m^2+1)/(-m);Par_exp=Par_exp+BE*u(X,Y); d1 = 0;X=X+1;Y=Y-1;endendif d1>=0 && d2>0 && d2<deltaif judgement(X, Y-1, n)==0break;elseBF = delta*sqrt(m^2+1)/(-m);Par_exp=Par_exp+BF*u(X,Y); d1 = d2;X=X;Y=Y-1;endendif d1<0if judgement(X,Y-1, n)==0break;elseMT = d2*sqrt(m^2+1);Par_exp=Par_exp+MT*u(X,Y); d1 = d2;X=X;Y=Y-1;endendend
end
end

衰减校正因子求解的matlab代码如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 非均匀衰减校正,使用chang氏修正法                   %
% F(x,y) = f(x,y)*c(x,y)                             %
% c(x,y) = [1/m * sigma[exp(sigma(u*l))]]^-1         %
%2021/1/29                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
n = 128;    % 图像大小
s = 35;     % number of slice
c = zeros(n,n,s);          % 校正因子分布
delta=1;                   % delta：移动步长
num_t = 64;                % 投影角度个数 number of angle
theta = [0:5.625:359];     % theta：存储角度
scale = 0.2;               % pixel scale=0.2cm
u_w = 0.19;                % 能量为73keV的X射线在水中的线性衰减系数19m-1, cm-1;
CT = load('CT.mat').Q;
u = CT*u_w/1000+u_w;       % 将HU转换为衰减系数
for z = 1:sfor x = 1:nfor y = 1:nPar_exp = zeros(1,64);     % 定义存储distance的数组L = zeros(1,64);C = 0;  for i = 1:num_tPar_exp(i) = distance(x, y, theta(i), n, delta, u(:,:,z));L(i) = Par_exp(i)*scale;    % 转换为实际长度, cmC = C+exp(-L(i));endc(x,y,z) = (C/num_t)^-1;endend
end
save('attenuation_map.mat', 'c');
figure(1);
imshow3Dfull(CT);
figure(2);
imshow3Dfull(c);

基于深度学习方法实现SPECT放射性核素定量测量（三）相关推荐

基于深度学习方法实现SPECT放射性核素定量测量（一）
基于深度学习方法实现SPECT放射性核素定量测量(一)--课题介绍放射治疗是目前癌症三大治疗手段之一,当前,世界上约有70%的癌症患者使用放射治疗.放射治疗技术包括普通X射线治疗.硼中子俘获治疗.质 ...
基于深度学习方法实现SPECT放射性核素定量测量（二）
基于深度学习方法实现SPECT放射性核素定量测量(二)-GATE仿真 (1)GATE介绍 GATE(GEANT4 Application for Emission Tomography)由OpenGA ...
基于深度学习方法实现SPECT放射性核素定量测量（四）
基于深度学习实现SPECT放射性核素定量测量(四)-pix2pix网络结构搭建 (1)pix2pix网络结构在本次工作中使用GAN网络的一种pix2pix网络模型,实现重建图像计数分布到放射性核素活 ...
基于深度学习方法的3D数据合成
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达 3D 数据简介人们普遍认为,从单一角度合成 3D 数据是人类视觉的基本功能.但这对计算机视觉算法来说 ...
基于深度学习方法的推荐系统（转载）
原文:Deep Learning based Recommender System: A Survey and New Perspectives 作者:张帅, 新南威尔士大学翻译:沈春旭,清华大学 ...
nature | 基于深度学习方法的虚拟组织染色
研究背景组织病理学可以追溯到19世纪,它一直是病理学中使用的黄金标准诊断方法之一.如果在医学检查之后或在外科手术期间需要活组织检查,则需要从患者身上取出组织样本,然后将其切成微米薄片.这些病理学切片 ...
WACV2020:开源基于深度学习方法DeOccNet用来去除透视光场中的前景遮挡
作者信息最近,国防科技大学的一个研究小组提出了一种利用阵列相机去除前景遮挡成像的新方法作为国内外第一个基于深度学习的去遮挡成像工作,作者提出了掩模嵌入的方法来解决训练数据不足的问题,并建立了仿真和 ...
【SR汇总】基于深度学习方法
1.SRCNN.FSRCNN (Learning a Deep Convolutional Network for Image Super-Resolution, ECCV2014) (Acceler ...
3D目标检测深度学习方法中voxel-represetnation内容综述（三）
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达前言前两篇文章:3D目标检测深度学习方法中voxel-represetnation内容综述(一).3 ...

基于深度学习方法实现SPECT放射性核素定量测量（三）

基于深度学习方法实现SPECT放射性核素定量测量（三）-图像重建及衰减校正

（1）图像重建

（2）衰减校正

基于深度学习方法实现SPECT放射性核素定量测量（三）相关推荐

最新文章

热门文章