2D examples¶
These examples require raw data which are automatically
downloaded from the source repository by the script
example_helper.py
.
Please make sure that this script is present in the example
script folder.
Mie off-center cylinder¶
The in silico data set was created with the softare miefield. The data are 1D projections of an off-center cylinder of constant refractive index. The Born approximation is error-prone due to a relatively large radius of the cylinder (30 wavelengths) and a refractive index difference of 0.006 between cylinder and surrounding medium. The reconstruction of the refractive index with the Rytov approximation is in good agreement with the input data. When only 50 projections are used for the reconstruction, artifacts appear. These vanish when more projections are used for the reconstruction.
backprop_from_mie_2d_cylinder_offcenter.py
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 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | import matplotlib.pylab as plt
import numpy as np
import odtbrain as odt
from example_helper import load_data
# simulation data
sino, angles, cfg = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="sino_imag.txt",
f_sino_real="sino_real.txt",
f_angles="mie_angles.txt",
f_info="mie_info.txt")
A, size = sino.shape
# background sinogram computed with Mie theory
# miefield.GetSinogramCylinderRotation(radius, nmed, nmed, lD, lC, size, A,res)
u0 = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="u0_imag.txt",
f_sino_real="u0_real.txt")
# create 2d array
u0 = np.tile(u0, size).reshape(A, size).transpose()
# background field necessary to compute initial born field
# u0_single = mie.GetFieldCylinder(radius, nmed, nmed, lD, size, res)
u0_single = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="u0_single_imag.txt",
f_sino_real="u0_single_real.txt")
print("Example: Backpropagation from 2D Mie simulations")
print("Refractive index of medium:", cfg["nmed"])
print("Measurement position from object center:", cfg["lD"])
print("Wavelength sampling:", cfg["res"])
print("Performing backpropagation.")
# Set measurement parameters
# Compute scattered field from cylinder
radius = cfg["radius"] # wavelengths
nmed = cfg["nmed"]
ncyl = cfg["ncyl"]
lD = cfg["lD"] # measurement distance in wavelengths
lC = cfg["lC"] # displacement from center of image
size = cfg["size"]
res = cfg["res"] # px/wavelengths
A = cfg["A"] # number of projections
x = np.arange(size) - size / 2
X, Y = np.meshgrid(x, x)
rad_px = radius * res
phantom = np.array(((Y - lC * res)**2 + X**2) < rad_px**2,
dtype=np.float) * (ncyl - nmed) + nmed
# Born
u_sinB = (sino / u0 * u0_single - u0_single) # fake born
fB = odt.backpropagate_2d(u_sinB, angles, res, nmed, lD * res)
nB = odt.odt_to_ri(fB, res, nmed)
# Rytov
u_sinR = odt.sinogram_as_rytov(sino / u0)
fR = odt.backpropagate_2d(u_sinR, angles, res, nmed, lD * res)
nR = odt.odt_to_ri(fR, res, nmed)
# Rytov 50
u_sinR50 = odt.sinogram_as_rytov((sino / u0)[::5, :])
fR50 = odt.backpropagate_2d(u_sinR50, angles[::5], res, nmed, lD * res)
nR50 = odt.odt_to_ri(fR50, res, nmed)
# Plot sinogram phase and amplitude
ph = odt.sinogram_as_radon(sino / u0)
am = np.abs(sino / u0)
# prepare plot
vmin = np.min(np.array([phantom, nB.real, nR50.real, nR.real]))
vmax = np.max(np.array([phantom, nB.real, nR50.real, nR.real]))
fig, axes = plt.subplots(2, 3, figsize=(8, 5))
axes = np.array(axes).flatten()
phantommap = axes[0].imshow(phantom, vmin=vmin, vmax=vmax)
axes[0].set_title("phantom \n(non-centered cylinder)")
amplmap = axes[1].imshow(am, cmap="gray")
axes[1].set_title("amplitude sinogram \n(background-corrected)")
phasemap = axes[2].imshow(ph, cmap="coolwarm")
axes[2].set_title("phase sinogram [rad] \n(background-corrected)")
axes[3].imshow(nB.real, vmin=vmin, vmax=vmax)
axes[3].set_title("reconstruction (Born) \n(250 projections)")
axes[4].imshow(nR50.real, vmin=vmin, vmax=vmax)
axes[4].set_title("reconstruction (Rytov) \n(50 projections)")
axes[5].imshow(nR.real, vmin=vmin, vmax=vmax)
axes[5].set_title("reconstruction (Rytov) \n(250 projections)")
# color bars
cbkwargs = {"fraction": 0.045}
plt.colorbar(phantommap, ax=axes[0], **cbkwargs)
plt.colorbar(amplmap, ax=axes[1], **cbkwargs)
plt.colorbar(phasemap, ax=axes[2], **cbkwargs)
plt.colorbar(phantommap, ax=axes[3], **cbkwargs)
plt.colorbar(phantommap, ax=axes[4], **cbkwargs)
plt.colorbar(phantommap, ax=axes[5], **cbkwargs)
plt.tight_layout()
plt.show()
|
Mie cylinder with unevenly spaced angles¶
Angular weighting can significantly improve reconstruction quality when the angular projections are sampled at non-equidistant intervals [TPM81]. The in silico data set was created with the softare miefield. The data are 1D projections of a non-centered cylinder of constant refractive index 1.339 embedded in water with refractive index 1.333. The first column shows the used sinograms (missing angles are displayed as zeros) that were created from the original sinogram with 250 projections. The second column shows the reconstruction without angular weights and the third column shows the reconstruction with angular weights. The keyword argument weight_angles was introduced in version 0.1.1.
backprop_from_mie_2d_weights_angles.py
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 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 | import matplotlib.pylab as plt
import numpy as np
import unwrap
import odtbrain as odt
from example_helper import load_data
sino, angles, cfg = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_angles="mie_angles.txt",
f_sino_real="sino_real.txt",
f_sino_imag="sino_imag.txt",
f_info="mie_info.txt")
A, size = sino.shape
# background sinogram computed with Mie theory
# miefield.GetSinogramCylinderRotation(radius, nmed, nmed, lD, lC, size, A,res)
u0 = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="u0_imag.txt",
f_sino_real="u0_real.txt")
# create 2d array
u0 = np.tile(u0, size).reshape(A, size).transpose()
# background field necessary to compute initial born field
# u0_single = mie.GetFieldCylinder(radius, nmed, nmed, lD, size, res)
u0_single = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="u0_single_imag.txt",
f_sino_real="u0_single_real.txt")
print("Example: Backpropagation from 2D FDTD simulations")
print("Refractive index of medium:", cfg["nmed"])
print("Measurement position from object center:", cfg["lD"])
print("Wavelength sampling:", cfg["res"])
print("Performing backpropagation.")
# Set measurement parameters
# Compute scattered field from cylinder
radius = cfg["radius"] # wavelengths
nmed = cfg["nmed"]
ncyl = cfg["ncyl"]
lD = cfg["lD"] # measurement distance in wavelengths
lC = cfg["lC"] # displacement from center of image
size = cfg["size"]
res = cfg["res"] # px/wavelengths
A = cfg["A"] # number of projections
x = np.arange(size) - size / 2.0
X, Y = np.meshgrid(x, x)
rad_px = radius * res
phantom = np.array(((Y - lC * res)**2 + X**2) < rad_px **
2, dtype=np.float) * (ncyl - nmed) + nmed
u_sinR = odt.sinogram_as_rytov(sino / u0)
# Rytov 200 projections
# remove 50 projections from total of 250 projections
remove200 = np.argsort(angles % .0002)[:50]
angles200 = np.delete(angles, remove200, axis=0)
u_sinR200 = np.delete(u_sinR, remove200, axis=0)
ph200 = unwrap.unwrap(np.angle(sino / u0))
ph200[remove200] = 0
fR200 = odt.backpropagate_2d(u_sinR200, angles200, res, nmed, lD*res)
nR200 = odt.odt_to_ri(fR200, res, nmed)
fR200nw = odt.backpropagate_2d(u_sinR200, angles200, res, nmed, lD*res,
weight_angles=False)
nR200nw = odt.odt_to_ri(fR200nw, res, nmed)
# Rytov 50 projections
remove50 = np.argsort(angles % .0002)[:200]
angles50 = np.delete(angles, remove50, axis=0)
u_sinR50 = np.delete(u_sinR, remove50, axis=0)
ph50 = unwrap.unwrap(np.angle(sino / u0))
ph50[remove50] = 0
fR50 = odt.backpropagate_2d(u_sinR50, angles50, res, nmed, lD*res)
nR50 = odt.odt_to_ri(fR50, res, nmed)
fR50nw = odt.backpropagate_2d(u_sinR50, angles50, res, nmed, lD*res,
weight_angles=False)
nR50nw = odt.odt_to_ri(fR50nw, res, nmed)
# prepare plot
kw_ri = {"vmin": 1.330,
"vmax": 1.340}
kw_ph = {"vmin": np.min(np.array([ph200, ph50])),
"vmax": np.max(np.array([ph200, ph50])),
"cmap": "coolwarm"}
fig, axes = plt.subplots(2, 3, figsize=(8, 4))
axes = np.array(axes).flatten()
phmap = axes[0].imshow(ph200, **kw_ph)
axes[0].set_title("Phase sinogram (200 proj.)")
rimap = axes[1].imshow(nR200nw.real, **kw_ri)
axes[1].set_title("RI without angular weights")
axes[2].imshow(nR200.real, **kw_ri)
axes[2].set_title("RI with angular weights")
axes[3].imshow(ph50, **kw_ph)
axes[3].set_title("Phase sinogram (50 proj.)")
axes[4].imshow(nR50nw.real, **kw_ri)
axes[4].set_title("RI without angular weights")
axes[5].imshow(nR50.real, **kw_ri)
axes[5].set_title("RI with angular weights")
# color bars
cbkwargs = {"fraction": 0.045,
"format": "%.3f"}
plt.colorbar(phmap, ax=axes[0], **cbkwargs)
plt.colorbar(phmap, ax=axes[3], **cbkwargs)
plt.colorbar(rimap, ax=axes[1], **cbkwargs)
plt.colorbar(rimap, ax=axes[2], **cbkwargs)
plt.colorbar(rimap, ax=axes[5], **cbkwargs)
plt.colorbar(rimap, ax=axes[4], **cbkwargs)
plt.tight_layout()
plt.show()
|
Mie cylinder with incomplete angular coverage¶
This example illustrates how the backpropagation algorithm of ODTbrain handles incomplete angular coverage. All examples use 100 projections at 100%, 60%, and 40% total angular coverage. The keyword argument weight_angles that invokes angular weighting is set to True by default. The in silico data set was created with the softare miefield. The data are 1D projections of a non-centered cylinder of constant refractive index 1.339 embedded in water with refractive index 1.333. The first column shows the used sinograms (missing angles are displayed as zeros) that were created from the original sinogram with 250 projections. The second column shows the reconstruction without angular weights and the third column shows the reconstruction with angular weights. The keyword argument weight_angles was introduced in version 0.1.1.
A 180 degree coverage results in a good reconstruction of the object. Angular weighting as implemented in the backpropagation algorithm of ODTbrain automatically addresses uneven and incomplete angular coverage.
backprop_from_mie_2d_incomplete_coverage.py
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 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 | import matplotlib.pylab as plt
import numpy as np
import odtbrain as odt
from example_helper import load_data
sino, angles, cfg = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_angles="mie_angles.txt",
f_sino_real="sino_real.txt",
f_sino_imag="sino_imag.txt",
f_info="mie_info.txt")
A, size = sino.shape
# background sinogram computed with Mie theory
# miefield.GetSinogramCylinderRotation(radius, nmed, nmed, lD, lC, size, A,res)
u0 = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="u0_imag.txt",
f_sino_real="u0_real.txt")
# create 2d array
u0 = np.tile(u0, size).reshape(A, size).transpose()
# background field necessary to compute initial born field
# u0_single = mie.GetFieldCylinder(radius, nmed, nmed, lD, size, res)
u0_single = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
f_sino_imag="u0_single_imag.txt",
f_sino_real="u0_single_real.txt")
print("Example: Backpropagation from 2D FDTD simulations")
print("Refractive index of medium:", cfg["nmed"])
print("Measurement position from object center:", cfg["lD"])
print("Wavelength sampling:", cfg["res"])
print("Performing backpropagation.")
# Set measurement parameters
# Compute scattered field from cylinder
radius = cfg["radius"] # wavelengths
nmed = cfg["nmed"]
ncyl = cfg["ncyl"]
lD = cfg["lD"] # measurement distance in wavelengths
lC = cfg["lC"] # displacement from center of image
size = cfg["size"]
res = cfg["res"] # px/wavelengths
A = cfg["A"] # number of projections
x = np.arange(size) - size / 2.0
X, Y = np.meshgrid(x, x)
rad_px = radius * res
phantom = np.array(((Y - lC * res)**2 + X**2) < rad_px **
2, dtype=np.float) * (ncyl - nmed) + nmed
u_sinR = odt.sinogram_as_rytov(sino / u0)
# Rytov 100 projections evenly distributed
removeeven = np.argsort(angles % .002)[:150]
angleseven = np.delete(angles, removeeven, axis=0)
u_sinReven = np.delete(u_sinR, removeeven, axis=0)
pheven = odt.sinogram_as_radon(sino / u0)
pheven[removeeven] = 0
fReven = odt.backpropagate_2d(u_sinReven, angleseven, res, nmed, lD * res)
nReven = odt.odt_to_ri(fReven, res, nmed)
fRevennw = odt.backpropagate_2d(
u_sinReven, angleseven, res, nmed, lD * res, weight_angles=False)
nRevennw = odt.odt_to_ri(fRevennw, res, nmed)
# Rytov 100 projections more than 180
removemiss = 249 - \
np.concatenate((np.arange(100), 100 + np.arange(150)[::3]))
anglesmiss = np.delete(angles, removemiss, axis=0)
u_sinRmiss = np.delete(u_sinR, removemiss, axis=0)
phmiss = odt.sinogram_as_radon(sino / u0)
phmiss[removemiss] = 0
fRmiss = odt.backpropagate_2d(u_sinRmiss, anglesmiss, res, nmed, lD * res)
nRmiss = odt.odt_to_ri(fRmiss, res, nmed)
fRmissnw = odt.backpropagate_2d(
u_sinRmiss, anglesmiss, res, nmed, lD * res, weight_angles=False)
nRmissnw = odt.odt_to_ri(fRmissnw, res, nmed)
# Rytov 100 projections less than 180
removebad = 249 - np.arange(150)
anglesbad = np.delete(angles, removebad, axis=0)
u_sinRbad = np.delete(u_sinR, removebad, axis=0)
phbad = odt.sinogram_as_radon(sino / u0)
phbad[removebad] = 0
fRbad = odt.backpropagate_2d(u_sinRbad, anglesbad, res, nmed, lD * res)
nRbad = odt.odt_to_ri(fRbad, res, nmed)
fRbadnw = odt.backpropagate_2d(
u_sinRbad, anglesbad, res, nmed, lD * res, weight_angles=False)
nRbadnw = odt.odt_to_ri(fRbadnw, res, nmed)
# prepare plot
kw_ri = {"vmin": np.min(np.array([phantom, nRmiss.real, nReven.real])),
"vmax": np.max(np.array([phantom, nRmiss.real, nReven.real]))}
kw_ph = {"vmin": np.min(np.array([pheven, phmiss])),
"vmax": np.max(np.array([pheven, phmiss])),
"cmap": "coolwarm"}
fig, axes = plt.subplots(3, 3, figsize=(8, 6.5))
axes[0, 0].set_title("100% coverage ({} proj.)".format(angleseven.shape[0]))
phmap = axes[0, 0].imshow(pheven, **kw_ph)
axes[0, 1].set_title("RI without angular weights")
rimap = axes[0, 1].imshow(nRevennw.real, **kw_ri)
axes[0, 2].set_title("RI with angular weights")
rimap = axes[0, 2].imshow(nReven.real, **kw_ri)
axes[1, 0].set_title("60% coverage ({} proj.)".format(anglesmiss.shape[0]))
axes[1, 0].imshow(phmiss, **kw_ph)
axes[1, 1].set_title("RI without angular weights")
axes[1, 1].imshow(nRmissnw.real, **kw_ri)
axes[1, 2].set_title("RI with angular weights")
axes[1, 2].imshow(nRmiss.real, **kw_ri)
axes[2, 0].set_title("40% coverage ({} proj.)".format(anglesbad.shape[0]))
axes[2, 0].imshow(phbad, **kw_ph)
axes[2, 1].set_title("RI without angular weights")
axes[2, 1].imshow(nRbadnw.real, **kw_ri)
axes[2, 2].set_title("RI with angular weights")
axes[2, 2].imshow(nRbad.real, **kw_ri)
# color bars
cbkwargs = {"fraction": 0.045,
"format": "%.3f"}
plt.colorbar(phmap, ax=axes[0, 0], **cbkwargs)
plt.colorbar(phmap, ax=axes[1, 0], **cbkwargs)
plt.colorbar(phmap, ax=axes[2, 0], **cbkwargs)
plt.colorbar(rimap, ax=axes[0, 1], **cbkwargs)
plt.colorbar(rimap, ax=axes[1, 1], **cbkwargs)
plt.colorbar(rimap, ax=axes[2, 1], **cbkwargs)
plt.colorbar(rimap, ax=axes[0, 2], **cbkwargs)
plt.colorbar(rimap, ax=axes[1, 2], **cbkwargs)
plt.colorbar(rimap, ax=axes[2, 2], **cbkwargs)
plt.tight_layout()
plt.show()
|
FDTD cell phantom¶
The in silico data set was created with the FDTD software meep. The data are 1D projections of a 2D refractive index phantom. The reconstruction of the refractive index with the Rytov approximation is in good agreement with the phantom that was used in the simulation.
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 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 | import matplotlib.pylab as plt
import numpy as np
import odtbrain as odt
from example_helper import load_data
sino, angles, phantom, cfg = load_data("fdtd_2d_sino_A100_R13.zip",
f_angles="fdtd_angles.txt",
f_sino_imag="fdtd_imag.txt",
f_sino_real="fdtd_real.txt",
f_info="fdtd_info.txt",
f_phantom="fdtd_phantom.txt",
)
print("Example: Backpropagation from 2D FDTD simulations")
print("Refractive index of medium:", cfg["nm"])
print("Measurement position from object center:", cfg["lD"])
print("Wavelength sampling:", cfg["res"])
print("Performing backpropagation.")
# Apply the Rytov approximation
sino_rytov = odt.sinogram_as_rytov(sino)
# perform backpropagation to obtain object function f
f = odt.backpropagate_2d(uSin=sino_rytov,
angles=angles,
res=cfg["res"],
nm=cfg["nm"],
lD=cfg["lD"] * cfg["res"]
)
# compute refractive index n from object function
n = odt.odt_to_ri(f, res=cfg["res"], nm=cfg["nm"])
# compare phantom and reconstruction in plot
fig, axes = plt.subplots(1, 3, figsize=(8, 2.8))
axes[0].set_title("FDTD phantom")
axes[0].imshow(phantom, vmin=phantom.min(), vmax=phantom.max())
sino_phase = np.unwrap(np.angle(sino), axis=1)
axes[1].set_title("phase sinogram")
axes[1].imshow(sino_phase, vmin=sino_phase.min(), vmax=sino_phase.max(),
aspect=sino.shape[1] / sino.shape[0],
cmap="coolwarm")
axes[1].set_xlabel("detector")
axes[1].set_ylabel("angle [rad]")
axes[2].set_title("reconstructed image")
axes[2].imshow(n.real, vmin=phantom.min(), vmax=phantom.max())
# set y ticks for sinogram
labels = np.linspace(0, 2 * np.pi, len(axes[1].get_yticks()))
labels = ["{:.2f}".format(i) for i in labels]
axes[1].set_yticks(np.linspace(0, len(angles), len(labels)))
axes[1].set_yticklabels(labels)
plt.tight_layout()
plt.show()
|