介紹

扇形束CT重建系統(tǒng)矩陣的生成通常采用網(wǎng)格法，即根據(jù)射線與像素相交的長(zhǎng)度來(lái)確定該像素對(duì)射線投影的貢獻(xiàn)。網(wǎng)格法是一種精確但耗時(shí)的方法。最近Wang等提出根據(jù)體素中心在探測(cè)器上的投影位置與其相鄰射線在探測(cè)器投影位置之差來(lái)表示體素對(duì)射線投影的貢獻(xiàn)。作者表示此種方法生成的系統(tǒng)矩陣精度高且并行性好。文獻(xiàn)參考：

[1] Wang, X. , Li, S. , & Jiang, C. . (2021). Fast system matrix calculation based on voxel projection in cone-beam ct. Optik - International Journal for Light and Electron Optics, 231(2), 166422.

此博文僅實(shí)現(xiàn)二維扇形束以供讀者參考，三維錐束矩陣可在此基礎(chǔ)上進(jìn)行擴(kuò)展。

Python代碼實(shí)現(xiàn)

代碼1（非并行）：
將系統(tǒng)矩陣生成過(guò)程分解為每個(gè)投影角度中的每個(gè)像素對(duì)射線貢獻(xiàn)的求解，方便理解。由于每個(gè)投影角度都對(duì)圖像的所有像素進(jìn)行了遍歷，耗時(shí)較長(zhǎng)。

import numpy as np import math import pydicom from matplotlib import pyplot as pltdef readdat(file,shape):fid = open(file,'rb')nefid = fid.read()fid.close()neimage = np.reshape(np.frombuffer(nefid,dtype=np.float32),shape,'F')return neimage SDD = 1250#962.9 SOD = 980#570 ReconCentre = [0,0] channel_spacing = 0.5 channel_num = int(512) image_size = 256#512 pixel_spacing = channel_spacing*SOD/SDD step = 1 # 旋轉(zhuǎn)間隔（單位：deg） view_num = 360#720 fan_beam_angle = math.atan((channel_num/2*pixel_spacing-pixel_spacing/2)/SDD) * 2'''每個(gè)探測(cè)器單元的中心（每條Ray在探測(cè)器上的landpoint）到探測(cè)器中心的距離''' Ray2DetCentre = [] for i in range(1,channel_num+1):if i<=channel_num/2: # 探測(cè)器左側(cè)temp = -channel_spacing/2 + channel_spacing*(i-channel_num/2)else: # 探測(cè)器右側(cè)temp = channel_spacing/2 + channel_spacing*(i-1-channel_num/2)Ray2DetCentre.append(temp)'''計(jì)算被掃描物體每個(gè)體素的坐標(biāo)''' VoxelCentreCoordinate = dict() temp_dict = {'x':[], 'y':[]} VoxelCentreCoordinate.update(temp_dict) row_count = 0 for i in range(1,image_size**2+1):col_count = i%image_sizeif (i-1)%image_size==0:row_count += 1if i%image_size==0:col_count = image_size# 坐標(biāo)原點(diǎn)在左上角VoxelCentreCoordinate['x'].append((col_count-1)*pixel_spacing+pixel_spacing/2)VoxelCentreCoordinate['y'].append((row_count-1)*pixel_spacing+pixel_spacing/2) # 將坐標(biāo)原點(diǎn)移動(dòng)到中心(HFS) VoxelCentreCoordinate['x'] = list(np.array(VoxelCentreCoordinate['x']) - (image_size/2 * pixel_spacing)) VoxelCentreCoordinate['y'] = list(np.array(VoxelCentreCoordinate['y']) - (image_size/2 * pixel_spacing)) # 將坐標(biāo)原點(diǎn)移動(dòng)到中心(FFS) VoxelCentreCoordinate['x'] = list(-1*np.array(VoxelCentreCoordinate['x']))'''讀入待投影的圖像''' image = readdat(r'D:\RadonTrans\papers\code\phantom.dat',[256,256]) image = image.T image = image.flatten() proj = np.zeros((channel_num,int(view_num/step)))SourceCoordinate = dict() temp_dict = {'x':[],'y':[]} SourceCoordinate.update(temp_dict)'''生成系統(tǒng)矩陣''' for i in range(1,view_num,step):print(i)sys_m = dict()temp_dict = {'weight':dict(), 'index':dict()}sys_m.update(temp_dict)theta = (i-1)/view_num*360*math.pi/180 SourceCoordinate['x'].append(-SOD*math.cos(theta))SourceCoordinate['y'].append(-SOD*math.sin(theta))voxel_count = 0h_landpoint_temp = []for j in range(image_size**2):h = VoxelCentreCoordinate['y'][j]*math.cos(theta) - VoxelCentreCoordinate['x'][j]*math.sin(theta)# 體素j的中心在探測(cè)器上的landpoint距離探測(cè)器中心的距離（以探測(cè)器中心為原點(diǎn)，探測(cè)器右側(cè)距離為正，探測(cè)器左側(cè)距離為負(fù)）h_landpoint = h * SDD / math.sqrt((SourceCoordinate['x'][int((i-1)/step)]-VoxelCentreCoordinate['x'][j])**2+(SourceCoordinate['y'][int((i-1)/step)]-VoxelCentreCoordinate['y'][j])**2-h**2)h_landpoint_temp.append(h_landpoint)# 判斷當(dāng)前體素j是否在扇形束范圍之外 或者 當(dāng)前體素j的landpoint是否超出探測(cè)器邊界if abs(h_landpoint) > channel_spacing*channel_num/2:continuevoxel_count += 1# 將voxel權(quán)重及索引分別存入距離該voxel最近的兩根射線idx = int(np.floor((h_landpoint-(Ray2DetCentre[0]-channel_spacing/2)) / channel_spacing))try:sys_m['weight'][idx].append(abs(h_landpoint-Ray2DetCentre[idx])/channel_spacing)temp = abs(h_landpoint-Ray2DetCentre[idx])/channel_spacingsys_m['index'][idx].append(j)except:sys_m['weight'][idx] = []temp = abs(h_landpoint-Ray2DetCentre[idx-1])/channel_spacingsys_m['weight'][idx].append(abs(h_landpoint-Ray2DetCentre[idx-1])/channel_spacing)sys_m['index'][idx] = []sys_m['index'][idx].append(j)try:sys_m['weight'][idx+1].append(1-temp)sys_m['index'][idx+1].append(j)except:sys_m['weight'][idx+1] = []sys_m['weight'][idx+1].append(1-temp) sys_m['index'][idx+1] = []sys_m['index'][idx+1].append(j) channel_index = np.array(sorted(sys_m['index'].keys()))for m in channel_index:proj[m, int((i-1)/step)] = np.sum(np.array(sys_m['weight'][m]) * image[np.array(sys_m['index'][m])])

代碼2：（所有像素并行計(jì)算）
通過(guò)矩陣運(yùn)算同時(shí)計(jì)算二維圖像的所有像素對(duì)射線的貢獻(xiàn)（三維圖像有更多體素效果將更加顯著），所有像素并行計(jì)算，效率大幅提升。

import numpy as np import math import pydicom from matplotlib import pyplot as pltdef readdat(file,shape):fid = open(file,'rb')nefid = fid.read()fid.close()neimage = np.reshape(np.frombuffer(nefid,dtype=np.float32),shape,'F')return neimage SDD = 1250#962.9 SOD = 980#570 ReconCentre = [0,0] channel_spacing = 0.5 channel_num = int(512) image_size = 256#512 pixel_spacing = channel_spacing*SOD/SDD step = 1 view_num = 360#720 fan_beam_angle = math.atan((channel_num/2*pixel_spacing-pixel_spacing/2)/SDD) * 2'''每個(gè)探測(cè)器單元的中心（每條Ray在探測(cè)器上的landpoint）到探測(cè)器中心的距離''' Ray2DetCentre = [] for i in range(1,channel_num+1):if i<=channel_num/2: # 探測(cè)器左側(cè)temp = -channel_spacing/2 + channel_spacing*(i-channel_num/2)else: # 探測(cè)器右側(cè)temp = channel_spacing/2 + channel_spacing*(i-1-channel_num/2)Ray2DetCentre.append(temp) Ray2DetCentre = np.array(Ray2DetCentre)'''計(jì)算被掃描物體每個(gè)體素的坐標(biāo)''' VoxelCentreCoordinate = dict() temp_dict = {'x':[], 'y':[]} VoxelCentreCoordinate.update(temp_dict) row_count = 0 for i in range(1,image_size**2+1):col_count = i%image_sizeif (i-1)%image_size==0:row_count += 1if i%image_size==0:col_count = image_size# 坐標(biāo)原點(diǎn)在左上角VoxelCentreCoordinate['x'].append((col_count-1)*pixel_spacing+pixel_spacing/2)VoxelCentreCoordinate['y'].append((row_count-1)*pixel_spacing+pixel_spacing/2) # 將坐標(biāo)原點(diǎn)移動(dòng)到中心(HFS) VoxelCentreCoordinate['x'] = np.array(VoxelCentreCoordinate['x']) - (image_size/2 * pixel_spacing) VoxelCentreCoordinate['y'] = np.array(VoxelCentreCoordinate['y']) - (image_size/2 * pixel_spacing) # 將坐標(biāo)原點(diǎn)移動(dòng)到中心(FFS) VoxelCentreCoordinate['x'] = -1*VoxelCentreCoordinate['x']'''讀入待投影的圖像''' image = readdat(r'D:\phantom.dat',[256,256]) image = image.T image = image.flatten() proj = np.zeros((channel_num,int(view_num/step)))SourceCoordinate = dict() temp_dict = {'x':[],'y':[]} SourceCoordinate.update(temp_dict) sino = np.zeros((channel_num,view_num)) '''生成系統(tǒng)矩陣''' for i in range(1,view_num,step):print(i)sys_m = dict()temp_dict = {'weight':dict(), 'index':dict()}sys_m.update(temp_dict)theta = (i-1)/view_num*360*math.pi/180 SourceCoordinate['x'].append(-SOD*math.cos(theta))SourceCoordinate['y'].append(-SOD*math.sin(theta))voxel_count = 0h_landpoint_temp = []h = VoxelCentreCoordinate['y']*math.cos(theta) - VoxelCentreCoordinate['x']*math.sin(theta) # (65536,)h_landpoint = h * SDD / np.sqrt((SourceCoordinate['x'][int((i-1)/step)]-VoxelCentreCoordinate['x'])**2+(SourceCoordinate['y'][int((i-1)/step)]-VoxelCentreCoordinate['y'])**2-h**2)idx1 = np.floor((h_landpoint-(Ray2DetCentre[0]-channel_spacing/2)) / channel_spacing).astype(np.int64) # (65536,)idx2 = (1 + idx1).astype(np.int64)weight1 = abs(h_landpoint-Ray2DetCentre[idx1])/channel_spacing # (65536,)weight2 = 1 - weight1for j in range(channel_num):if j in idx1:temp = idx1==jtemp = temp+0else:continueif j in idx2:temp2 = idx2==jtemp2 = temp2+0sino[j,i-1] = np.sum(temp*weight1*image + temp2*weight2*image)else:sino[j,i-1] = np.sum(temp*weight1*image)

結(jié)果展示

圖1 256*256大小的’Shepp-Logan’圖像

圖2 非并行扇形束正投結(jié)果（左），并行扇形束正投結(jié)果（右）

小結(jié)

本博客展示的代碼塊僅對(duì)二維圖像的像素進(jìn)行了并行處理，投影角度和探測(cè)器通道這兩個(gè)維度上應(yīng)該也可以進(jìn)行并行計(jì)算。歡迎大家討論指正。

總結(jié)

以上是生活随笔為你收集整理的扇形束CT重建快速生成系统矩阵（system matrix）的Python实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。