ODL使用笔记

最新推荐文章于 2023-11-20 19:42:44 发布

小屁孩__

最新推荐文章于 2023-11-20 19:42:44 发布

阅读量1.3k

点赞数 2

文章标签：笔记

本文链接：https://blog.csdn.net/TZ_li_gen/article/details/131513317

版权

ODL简介

Operator Discretization Library (ODL) is a Python library for fast prototyping focusing on (but not restricted to) inverse problems. ODL is being developed at KTH Royal Institute of Technology, Stockholm, and Centrum Wiskunde & Informatica (CWI), Amsterdam.

运算符离散化库（ODL）是一个用于快速原型设计的 Python 库，专注于（但不限于）逆问题。 ODL正在KTH Royal Institute of Technology, Stockholm, and Centrum Wiskunde & Informatica (CWI), Amsterdam.开发。

ODL的主要目的是使数学家和应用科学家能够对现实世界的问题使用不同的数值方法，而不必自下而上地实现所有必要的部分。这是通过封装所有特定于应用程序的部件的运算符结构以及通常期望运算符、数据和附加参数的求解器的高级公式来实现的。这种方法的主要优点是

通过简单地切换运算符和数据，可以用相同的方法（例如TV）解决不同的问题。
同样的问题可以通过简单地调用不同的求解器来解决不同的方法。
求解器和特定于应用程序的代码只需在一个位置编写一次，并且可以单独测试。
添加新的应用程序或解决方案方法变得容易得多。

ODL 实现了许多抽象的数学概念，如集合、向量空间和运算符。

ODL安装

conda install -c odlgroup odl matplotlib pytest scikit-image

额外的依赖项

以下包是可选的，它们扩展了 ODL 的功能。其中一些需要 pip 才能安装。请参阅安装点进一步说明。

使用 matplotlib 的图像和绘图显示功能：
```
conda install matplotlib
```
使用 FFTW 的更快 FFT 后端（目前不在主流 conda 中）：
- 安装 FFTW C 库版本 3（所有可能的精度）。使用 Linux 包管理器执行此任务，或分别参阅 Windows 或 MacOS 说明。
- 通过运行以下命令安装 python 后端 pyFFTW：
```
pip install pyfftw
```
使用 PyWavelet 的 Wavelet 变换（目前不在主流 conda 中）：
```
pip install pywavelets
```
使用scikit-image进行光线变换的简单后端：
```
conda install scikit-image
```
快速 ASTRA 射线转换后端：
```
conda install -c astra-toolbox astra-toolbox
```
如果这不起作用，或者您想要更新的版本，请查看下面的 ASTRA 用于 X 射线断层扫描的部分。
绑定到 ProxImaL 凸优化包，CVXPY 的扩展：
```
pip install proximal
```

安装ODL常见的错误

1.module ‘odl’ has no attribute ‘RectGrid’

pip install odl==0.7.0

2.numpy需要小于1.24.0以下的版本

ODL使用

ODL是一个可以快速求解逆问题的python库。而计算机断层成像的解析或迭代重建方法本质上也是对目标函数进行求逆的过程，odl库中的tomo模块可以很好的解决2D/3D CT图像重建问题且可以利用GPU加速。有关odl库的详细文档描述见：

https://odl.readthedocs.io/

调用odl.tomo进行正反投操作会调用到astra工具箱，有关astra的详细文档描述见：

http://www.astra-toolbox.com/docs/

有关astra工具箱调用代码示例的ppt可见：

https://www.astra-toolbox.com/files/misc/ICTMS2019/20190722_ICTMS_ASTRA_workshop_2d.pdf

1. FBP（2D）

生成平行束投影的几何：
若重建空间recon_space为2D数据则parallel_beam_geometry调用Parallel2dGeometry，若recon_space为3D数据则parallel_beam_geometry调用Parallel3dAxisGeometry。

geometry = odl.tomo.parallel_beam_geometry(reco_space)

使用anaconda3安装odl，parallel_beam_geometry函数在odl库中的位置：

C:\ProgramData\Anaconda3\Lib\site-packages\odl\tomo\geometry\parallel.py

生成扇形束投影的几何：
若重建空间recon_space为2D数据则cone_beam_geometry调用FanFlatGeometry，若recon_space为3D数据则cone_beam_geometry调用ConeFlatGeometry。
若不设置num_angle的值则会使用默认的足够对数据进行完全采样的数值。

geometry = odl.tomo.cone_beam_geometry(reco_space,src_radius=40, det_radius=80,num_angles=None)

cone_beam_geometry函数在odl库中的位置：

C:\ProgramData\Anaconda3\Lib\site-packages\odl\tomo\geometry\conebeam.py

import odl
import matplotlib.pyplot as plt
import numpy as np
reco_space = odl.uniform_discr(min_pt=[-20, -20], max_pt=[20, 20], shape=[512, 512])
# parallel beam
geometry = odl.tomo.parallel_beam_geometry(reco_space)
# fan beam
#geometry = odl.tomo.cone_beam_geometry(reco_space,src_radius=40, det_radius=80,num_angles=None)
projector = odl.tomo.RayTransform(reco_space,geometry, impl='astra_cuda')
vol = odl.phantom.shepp_logan(projector.domain, modified=True)
proj = projector(vol)
backproj = projector.adjoint(proj) # 直接反投影
FBP = odl.tomo.fbp_op(projector) 
img = FBP(proj) # 滤波反投影
plt.imshow(np.flipud(proj.data.T),'gray')
plt.figure()
plt.imshow(np.flipud(img.data.T),'gray')
plt.figure()
plt.imshow(np.flipud(backproj.data.T),'gray')
plt.show()

Geometries in ASTRA

Volume geometry and volume data

vol_geom = astra.create_vol_geom(num_rows, num_cols)

all lengths and sizes defined relative to voxel size
voxel size defaults to 1 unit (square voxels)

rec_id = astra.data2d.create(`-vol`, vol_geom)

rec_id = astra.data2d.create(’-vol’, vol_geom) 

rec_id = astra.data2d.create(’-vol’, vol_geom, 0) 
astra.data2d.store(rec_id, 0) 

rec_id = astra.data2d.create(’-vol’, vol_geom, V) 
astra.data2d.store(rec_id, V) 

V = astra.data2d.get(rec_id) 

astra.data2d.delete(rec_id)

Projection geometry and data

trajectory of source and detector plane： sod

number of detector elements：探元个数

detector size：探测器大小

parallel projection geometry

请添加图片描述

angles = np.linspace(0, np.pi, 180, False)
proj_geom = astra.create_proj_geom(‘parallel’, det_width, det_count, angles)

fan-beam projection geometry

请添加图片描述

angles = np.linspace(0, 2*np.pi, 360, False)
proj_geom = astra.create_proj_geom(‘fanflat’, det_width, det_count, angles, source_origin_dist, origin_detector_dist)

projection data

proj_id = astra.data2d.create(’-sino’, proj_geom) 
proj_id = astra.data2d.create(’-sino’, proj_geom, 0) 
proj_id = astra.data2d.create(’-sino’, proj_geom, V)

Reconstruction – SIRT (CPU)

# Configure geometry
angles = np.linspace(0, 2*np.pi, nAngles, False)
proj_geom = astra.create_proj_geom(‘fanflat', 1.0, detectorCount, angles,
originToSource, originToDetector)
vol_geom = astra.create_vol_geom(vy, vx)
# Data objects for input, output
sino_id = astra.data2d.create('-sino', proj_geom, S)
rec_id = astra.data2d.create('-vol', vol_geom)
proj_id = astra.create_projector(‘strip_fanflat', proj_geom, vol_geom)
# Configure algorithm
cfg = astra.astra_dict('SIRT')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sino_id
cfg['ProjectorDataId'] = proj_id
alg_id = astra.algorithm.create(cfg)
# Run
astra.algorithm.run(alg_id, 100)
rec = astra.data2d.get(rec_id)

Reconstruction – SIRT (GPU)

# Configure geometry
angles = np.linspace(0, 2*np.pi, nAngles, False)
proj_geom = astra.create_proj_geom(‘fanflat', 1.0, detectorCount, angles,
originToSource, originToDetector)
vol_geom = astra.create_vol_geom(vy, vx)
# Data objects for input, output
proj_id = astra.data2d.create('-sino', proj_geom, S)
rec_id = astra.data2d.create('-vol', vol_geom)
# Configure algorithm
cfg = astra.astra_dict('SIRT_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = proj_id
alg_id = astra.algorithm.create(cfg)
# Run
astra.algorithm.run(alg_id, 100)
rec = astra.data2d.get(rec_id)

SciPy linear operator

vol_geom = astra.create_vol_geom(256, 256)
proj_geom = astra.create_proj_geom('parallel', 1.0, 256, np.linspace2(0, np.pi, 180, False))
proj_id = astra.create_projector('strip', proj_geom, vol_geom)
# Create OpTomo operator
W = astra.OpTomo(proj_id)
# Forward projection
s = W * P
s = s.reshape(astra.geom_size(proj_geom))
# Reconstruction using scipy’s lsqr
output = scipy.sparse.linal.lsqr(W, s.ravel(), iter_lim=100)
rec = output[0].reshape(astra.geom_size(vol_geom))

Flexible geometries Example uses (2D and 3D):

Three 2D parameters per projection: s, d, u
请添加图片描述

Fan beam – vector form

# One single projection angle
vectors = np.zeros((1, 6))
angle = 0.1
source_dist = 2000
# source
vectors[0,0] = np.sin(angle) * source_dist
vectors[0,1] = -np.cos(angle) * source_dist 
# center of detector
vectors[0,2] = 0
vectors[0,3] = 0
# vector from detector pixel 0 to 1
vectors[0,4] = np.cos(angle) * 1.0
vectors[0,5] = np.sin(angle) * 1.0
proj_geom = astra.create_proj_geom('fanflat_vec', 256, vectors)
# Parallel beam - vector form
proj_geom = astra.create_proj_geom('parallel_vec', 256, vectors)

vectors consists of a 6 element row vector per projection.

Three 2D parameters per projection:

r: Ray direction

d: detector center

u: detector pixel basis vector

Using utility function:
# Standard fan beam
proj_geom = astra.create_proj_geom('fanflat', 1.0, 256, angles, 2000, 0)
# Convert to vector form
proj_geom_vec = astra.geom_2vec(proj_geom)
# Center-of-rotation correction (by -3.5 pixels)
proj_geom_cor = astra.geom_postalignment(proj_geom, -3.5)

Using utility function:
# Standard fan beam
proj_geom = astra.create_proj_geom('fanflat', 1.0, 256, angles, 2000, 0)
# Convert to vector form
proj_geom_vec = astra.geom_2vec(proj_geom)
# Center-of-rotation correction (by -3.5 pixels)
proj_geom_cor = astra.geom_postalignment(proj_geom, -3.5)