Madagascrar之线性Radon变换
本文简要介绍如何在python中使用Madagascrar中的Radon变换模块。
命令行中输入sfradon可以查看Radon变换的相关参数，如下所示1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38NAME
sfradon
DESCRIPTION
High-resolution Radon transform.
SYNOPSIS
sfradon < in.rsf > out.rsf offset=offset.rsf adj=y inv=adj spk=inv verb=n np= dp= p0= nx= eps=1. ns=1 tol=1.e-6 perc=50.0 fact=0.5 ox= dx= parab=n x0=
1. niter=100
PARAMETERS
bool adj=y [y/n] if y, perform adjoint operation
float dp= p sampling (if adj=y)
float dx=
float eps=1.
float fact=0.5 percentage for sharpening
bool inv=adj [y/n] if y, perform inverse operation
int niter=100
int np= number of p values (if adj=y)
int ns=1 number of sharpening cycles
int nx= number of offsets (if adj=n)
string offset= auxiliary input file name
float ox=
float p0= p origin (if adj=y)
bool parab=n [y/n] if y, parabolic Radon transform
float perc=50.0 percentage for sharpening
bool spk=inv [y/n] if y, use spiking (hi-res) inversion
float tol=1.e-6 inversion tolerance
bool verb=n [y/n] verbosity flag
float x0=1. reference offset
USED IN
cwp/geo2006TimeShiftImagingCondition/flat
geo384s/hw5/pradon
geo384s/hw5/radon
milano/taupvel/cmp
rsf/school2012/demo
rsf/su/rsflab16
rsf/su/rsftaup
xjtu/test/myradon2
SOURCE
system/seismic/Mradon.
下面我们通过一个简单的实例来了解一下，如何使用线性Radon的正反变换。
首先，我们使用numpy生成一个包含两条水平直线、大小为(256, 256)的二维数组，两条直线分别位于网格深度64和192处。数据生成代码和图片如下所示：1
2
3
4
5nx = 256
nz = 256
image = np.zeros((nx, nz), dtype = np.float32)
image[64,:]=1
image[192,:]=1
接下来，使用m8r模块中的radon分别创建拉冬正变化和反变换的算子radon_op和iradon_op，其中我们必须要通过put方法来输入数据的网格间距d1、d2和起始位置o2信息，以及拉冬变换所需要的相关参数信息，比如变换后的采样数np（决定了变换后数据的网格大小）、p起始位置p0以及采样间隔dp，是否使用高精度变换spk等。
这里，我们将原始数据横向和纵向采样间隔均设置为1，并令$NP=55$，$dp=2*pmax/(NP-1)$。需要注意的是，在反变换中，我们需要将d2设置为正变换中的dp，即将数据第二个维度采样间隔设置为正变换中p的采样间隔，并将第二个维度的起始位置o2设置为正变换中p0(即此例中的-pmax)；然后设置radon变换中的nx为原数据第一个维度大小即可得到与原数据大小相同的二维数据。1
2
3
4
5
6
7
8
9
10
11d1 = 1
d2 = 1
o2 = 0
pmax=1
NP=55
dp=2*pmax/(NP-1)
radon_op = m8r.put(d1 = d1, d2 = d2, o2 = o2).radon(np=NP, p0=-pmax, dp=dp, spk=False)
iradon_op = m8r.put(d1 = d1, d2 = dp, o2 = -pmax).radon(adj=False,
nx=image.shape[0],
dx=d2, ox=0,
spk=False)
分别将正、反变换算子应用于原始数据和变换后的数据，即可相应地得到拉冬域数据和原始数据。但需要注意的是，由于python和C中行列优先级的不同，我们需要将numpy数组转置后再输入到算子中才可得到正确的正变换结果，相应的，反变换的结果转置后才为真正的原始数据排列方式。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19# 正变换
radon_m8r = radon_op.apply(image.T)
# 反变换
idata_m8r = iradon_op.apply(radon_m8r)
# 绘图
dx, dy = 0.5 * 180.0 / max(image.shape), 0.5 / radon_m8r.shape[0]
fig, ax = plt.subplots(1, 1)
ax.imshow(radon_m8r.T,
extent = (-dx, 180.0 + dx, -dy, radon_m8r.shape[0] + dy),
cmap = plt.cm.Greys, aspect ='auto')
ax.set_xlabel("Projection angle (deg)")
ax.set_ylabel("Projection position (pixels)")
plt.show()
plt.imshow(image, cmap=plt.cm.Greys)
plt.show()
plt.imshow(idata_m8r.T, cmap=plt.cm.Greys)
plt.show()
正变换结果：
反变换结果：