T1和T2图像及其计算-《医学图像处理》小作业一-Python代码/Matlab代码

天津中医药大学-20级医学信息工程   教师：王翌   学生：邓集亲

声明：本文章中所涉及的代码并非我个人独立完成的成果。

        在撰写的过程中，我广泛地吸取了前一辈人，尤其是学长学姐们的宝贵学习经验。通过深入研究他们的学习轨迹，以及查阅和分析了众多相关的理论文献与资料，并在王老师的悉心指导下，经过反复的实验、调试与优化，最终得以总结完成本文所展现的代码。

        建议各位学弟学妹们，先根据王老师的授课内容，独立思考一下应该如何完成。如果实在是有理解上的困难，不知道从何下手，再来学习和参考本文。

        学长我是用Python写的，如果你使用MATLAB，同样可以参考此代码进行翻译。我在代码中加入了许多注释和测试环节，以便于理解。由于学长能力有限，代码中或许存在一些疏漏或错误，请批判性地参考。

实验一：和图像

作业要求：

1.  参考“医学成像技术”课的内容，基于Bloch方程，从若干原始图像中计算出和图像。

例子：

实验图片：

见附件中相关内容

任务目标：在data中有可供生成T1分布图和T2分布图的两组文件。尝试用它们分别计算出一幅大脑切片图的T1和T2分布图出来，并提供脑部区域任意一个点的原始数据与拟合曲线的对比图（类似于如下的图片，用十字叉标出实际测量的点，再用另一条线表示拟合得到的曲线）。源代码也应一并提交。

数据及模型说明：

对用于计算分布图的各扫描图像而言，每个像素的信号满足。其中b是处于-1到1之间的未知数。对不同扫描图像中处于同一位置的像素，可认为三个未知数、和具有相同的取值。这些扫描图像使用了不同的值（每幅图所在文件夹的最后一个数即是以毫秒为单位的值），从而使得对应的像素具备了不同的数值。实际上这种单指数衰减模型仍是较为粗糙的，可参考所附的文章用双指数模型进行拟合。

对用于计算分布图的各扫描图像而言，每个像素的信号满足。对不同扫描图像中处于同一位置的像素，可认为两个未知数和具有相同的取值（这与上一段的不是一个量）。编号1~~8的各图像中，第i个对应于使用了毫秒的TE值。（更具体的TE值也存储在了各个文件中，用dicominfo读入图像的后台信息可翻找到TE值。）建议只使用偶数图像甚至是4，6，8三幅图像进行拟合，可以与使用所有图像进行拟合得到的结果做对比。

函数及程序说明：

用dicomread可以读出图像本身。dicominfo可以看到图像的扫描参数等信息（是的你眼尖的话能看到被扫描的就是本人）。

为了拟合参数可以用lsqcurvefit，这个函数实在不知道如何使用尽管来找我。

为了显示图片可以用imshow，为了画曲线和标点可以用plot。

显然图片中有些区域不是脑子只是空气，这部分像素中信号都很小且对其拟合或是值毫无意义。请考虑如何在结果中去除它们（例如直接赋值为0）以避免图像看着太乱。

Windows图片文件目录：

ImageSet文件目录下储存其他待处理图片

Python代码：

import numpy as np
import os
import pydicom
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit


images_num=20
all_images=[[]] * images_num

file_path=r"D:\deng\smalljob\job\data\T1map"
i=str(8)
j=0

#循环读入图片并转换
for n in range(8,1751):
    if(os.path.exists(file_path+'\\'+'CEST_wangyi_7_829_1_5_'+i)):
        for f in os.listdir(file_path+'\\'+'CEST_wangyi_7_829_1_5_'+i):
            if f.endswith('.dcm'):   #遍历读入

                all_images[j]=pydicom.dcmread(file_path+'\\'+'CEST_wangyi_7_829_1_5_'+i+'\\'+f)   #这里将图片读入是为了后面进行比较
                
                j=j+1
    temp=int(i)+1
    i=str(temp)

#新建一个三维数组，用来进行后续操作
images_shape=all_images[0].pixel_array.shape
images_array=np.zeros((len(all_images),images_shape[0],images_shape[1]))

#将原来读入的图像存储进新的数组
for n in range(0,20):
    images_array[n]=all_images[n].pixel_array
    #print(images_array[n][49][49])
    #plt.imshow(images_array[n])
    #print(images_array[n])



# 构造模型函数
def model_fun(x, *p):
    return p[0] * np.abs(1 - (1 - p[1]) * np.exp(-x / p[2]))

# 定义ti数组
ti = np.array([8, 10, 12, 14, 16, 18, 20, 25, 30, 40, 60, 80, 100, 300, 500, 750, 1000, 1250, 1500, 1750])

#设置初始值
p0=np.array([10000,-1,1500])

#定义上限与下限
lb=np.array([0,-1,100])
ub=np.array([20000,1,4000])

#新建一维数组，用来存储每张图像相同像素点位的像素值
y=[0.]* images_num
y=np.array(y)

#初始化数组


fit_result = np.zeros((64, 64, 3))


#file=open('D:\deng\\test.txt','w')
for i in range(0,64):
    for j in range(0,64):
        for n in range(0,20):
            
            temp=images_array[n][i][j]
            
            y[n]=temp
        if sum(y)/20<30:#过滤平均像素值偏低的像素点位，可自行调整，此处过滤条件为低于30
            continue

        popt,pcov=curve_fit(model_fun,ti,y,p0=p0,bounds=(lb,ub))
        #file.write(str(popt)+'\n')
        #print(popt)
        #print(pcov)
        fit_result[i,j,:] = popt 
        
fit_T1 = fit_result[:,:,2]

#x = ti
#for n in range(0,20):

#    temp=images_array[n][50][50]
            
#    y[n]=temp

#popt,pcov=curve_fit(model_fun,ti,y,p0=p0,bounds=(lb,ub))
#print(y)


#plt.figure()
#plt.plot(x, model_fun(x,*popt))
#plt.ylabel('Signal Intensity')
#plt.xlabel('TI (ms)')
#plt.legend(['Original Signal', 'Fitted Curve'])
#plt.plot(x, y)
#plt.grid(True)


plt.imshow(fit_T1,cmap='gray')
plt.colorbar()
plt.title('fit_T1')
plt.axis('off')

plt.show()

运行结果：（灰度图，可以看见，除了图中心的脑图外，周围出现了其他的干扰信息）

运行结果：（无过滤条件，代码中有些注释是学长为了测试所以才写的，根据需要自行测试或删除）

运行结果：（过滤条件：过滤平均像素值低于30的点，条件可自行调整，比如设置50或更大的值，效果可能会更好，因为我们目的是计算脑图部分，而不是噪声） 

 以下是计算T2图像

import numpy as np
import os
import pydicom
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

images_num=8
all_images=[[]] * images_num

file_path=r"D:\deng\smalljob\job\data\T2map"
i=str(1)
j=0

#循环读入图片并转换
for n in range(0,8):
    for f in os.listdir(file_path):
        if f=='000'+i+'.dcm':   #遍历读入

            all_images[j]=pydicom.dcmread(file_path+'\\'+f)   #这里将图片读入是为了后面进行比较

            j=j+1
    temp=int(i)+1
    i=str(temp)

images_shape=all_images[0].pixel_array.shape
images_array=np.zeros((3,images_shape[0],images_shape[1]))
all_images_temp=[[]]*3

m=0
for n in range(2,5):
    images_array[m]=all_images[n*2-1].pixel_array
    all_images_temp[m]=all_images[n*2-1]
    m=m+1

# 构造模型函数
def model_fun(x,*p):
    return p[0]*np.exp(-x/p[1])


TE =[0.]* 3
for n in range(0,3):
    TE[n] = (all_images_temp[n]).EchoTime

# 初始化拟合参数的初值
p0=np.array([5000,100])

# 定义下限和上限

lb=np.array([1000,1])
ub=np.array([20000,3000])

y=[0.]* 3
y=np.array(y)


fit_result = np.zeros((256, 256, 2))

for i in range(0,256):
    for j in range(0,256):
        for n in range(0,3):
            #if images_array[n][i][j] != 0:
            temp=images_array[n][i][j]
            
            y[n]=temp
        if sum(y)/20<30:
            continue
        popt,pcov=curve_fit(model_fun,TE,y,p0=p0,maxfev=80000,bounds=(lb,ub))
        
        fit_result[i,j,:] = popt 
        
fit_T2 = fit_result[:,:,1]


plt.imshow(fit_T2,cmap='gray')
plt.colorbar()
plt.title('fit_T2')
plt.axis('off')

plt.show()

运行结果：（灰度图，运行时间较长，不急慢慢等）

运行结果：（彩色图像，注意：不设置过滤条件，会报错！拟合次数超过限制！）