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.

_images/backprop_from_mie_2d_cylinder_offcenter.jpg

backprop_from_mie_2d_cylinder_offcenter.py

  1import matplotlib.pylab as plt
  2import numpy as np
  3
  4import odtbrain as odt
  5
  6from example_helper import load_data
  7
  8
  9# simulation data
 10sino, angles, cfg = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 11                              f_sino_imag="sino_imag.txt",
 12                              f_sino_real="sino_real.txt",
 13                              f_angles="mie_angles.txt",
 14                              f_info="mie_info.txt")
 15A, size = sino.shape
 16
 17# background sinogram computed with Mie theory
 18# miefield.GetSinogramCylinderRotation(radius, nmed, nmed, lD, lC, size, A,res)
 19u0 = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 20               f_sino_imag="u0_imag.txt",
 21               f_sino_real="u0_real.txt")
 22# create 2d array
 23u0 = np.tile(u0, size).reshape(A, size).transpose()
 24
 25# background field necessary to compute initial born field
 26# u0_single = mie.GetFieldCylinder(radius, nmed, nmed, lD, size, res)
 27u0_single = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 28                      f_sino_imag="u0_single_imag.txt",
 29                      f_sino_real="u0_single_real.txt")
 30
 31print("Example: Backpropagation from 2D Mie simulations")
 32print("Refractive index of medium:", cfg["nmed"])
 33print("Measurement position from object center:", cfg["lD"])
 34print("Wavelength sampling:", cfg["res"])
 35print("Performing backpropagation.")
 36
 37# Set measurement parameters
 38# Compute scattered field from cylinder
 39radius = cfg["radius"]  # wavelengths
 40nmed = cfg["nmed"]
 41ncyl = cfg["ncyl"]
 42
 43lD = cfg["lD"]  # measurement distance in wavelengths
 44lC = cfg["lC"]  # displacement from center of image
 45size = cfg["size"]
 46res = cfg["res"]  # px/wavelengths
 47A = cfg["A"]  # number of projections
 48
 49x = np.arange(size) - size / 2
 50X, Y = np.meshgrid(x, x)
 51rad_px = radius * res
 52phantom = np.array(((Y - lC * res)**2 + X**2) < rad_px**2,
 53                   dtype=np.float) * (ncyl - nmed) + nmed
 54
 55# Born
 56u_sinB = (sino / u0 * u0_single - u0_single)  # fake born
 57fB = odt.backpropagate_2d(u_sinB, angles, res, nmed, lD * res)
 58nB = odt.odt_to_ri(fB, res, nmed)
 59
 60# Rytov
 61u_sinR = odt.sinogram_as_rytov(sino / u0)
 62fR = odt.backpropagate_2d(u_sinR, angles, res, nmed, lD * res)
 63nR = odt.odt_to_ri(fR, res, nmed)
 64
 65# Rytov 50
 66u_sinR50 = odt.sinogram_as_rytov((sino / u0)[::5, :])
 67fR50 = odt.backpropagate_2d(u_sinR50, angles[::5], res, nmed, lD * res)
 68nR50 = odt.odt_to_ri(fR50, res, nmed)
 69
 70# Plot sinogram phase and amplitude
 71ph = odt.sinogram_as_radon(sino / u0)
 72
 73am = np.abs(sino / u0)
 74
 75# prepare plot
 76vmin = np.min(np.array([phantom, nB.real, nR50.real, nR.real]))
 77vmax = np.max(np.array([phantom, nB.real, nR50.real, nR.real]))
 78
 79fig, axes = plt.subplots(2, 3, figsize=(8, 5))
 80axes = np.array(axes).flatten()
 81
 82phantommap = axes[0].imshow(phantom, vmin=vmin, vmax=vmax)
 83axes[0].set_title("phantom \n(non-centered cylinder)")
 84
 85amplmap = axes[1].imshow(am, cmap="gray")
 86axes[1].set_title("amplitude sinogram \n(background-corrected)")
 87
 88phasemap = axes[2].imshow(ph, cmap="coolwarm")
 89axes[2].set_title("phase sinogram [rad] \n(background-corrected)")
 90
 91axes[3].imshow(nB.real, vmin=vmin, vmax=vmax)
 92axes[3].set_title("reconstruction (Born) \n(250 projections)")
 93
 94axes[4].imshow(nR50.real, vmin=vmin, vmax=vmax)
 95axes[4].set_title("reconstruction (Rytov) \n(50 projections)")
 96
 97axes[5].imshow(nR.real, vmin=vmin, vmax=vmax)
 98axes[5].set_title("reconstruction (Rytov) \n(250 projections)")
 99
100# color bars
101cbkwargs = {"fraction": 0.045}
102plt.colorbar(phantommap, ax=axes[0], **cbkwargs)
103plt.colorbar(amplmap, ax=axes[1], **cbkwargs)
104plt.colorbar(phasemap, ax=axes[2], **cbkwargs)
105plt.colorbar(phantommap, ax=axes[3], **cbkwargs)
106plt.colorbar(phantommap, ax=axes[4], **cbkwargs)
107plt.colorbar(phantommap, ax=axes[5], **cbkwargs)
108
109plt.tight_layout()
110plt.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.

_images/backprop_from_mie_2d_weights_angles.jpg

backprop_from_mie_2d_weights_angles.py

  1import matplotlib.pylab as plt
  2import numpy as np
  3import unwrap
  4
  5import odtbrain as odt
  6
  7from example_helper import load_data
  8
  9
 10sino, angles, cfg = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 11                              f_angles="mie_angles.txt",
 12                              f_sino_real="sino_real.txt",
 13                              f_sino_imag="sino_imag.txt",
 14                              f_info="mie_info.txt")
 15A, size = sino.shape
 16
 17# background sinogram computed with Mie theory
 18# miefield.GetSinogramCylinderRotation(radius, nmed, nmed, lD, lC, size, A,res)
 19u0 = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 20               f_sino_imag="u0_imag.txt",
 21               f_sino_real="u0_real.txt")
 22# create 2d array
 23u0 = np.tile(u0, size).reshape(A, size).transpose()
 24
 25# background field necessary to compute initial born field
 26# u0_single = mie.GetFieldCylinder(radius, nmed, nmed, lD, size, res)
 27u0_single = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 28                      f_sino_imag="u0_single_imag.txt",
 29                      f_sino_real="u0_single_real.txt")
 30
 31
 32print("Example: Backpropagation from 2D FDTD simulations")
 33print("Refractive index of medium:", cfg["nmed"])
 34print("Measurement position from object center:", cfg["lD"])
 35print("Wavelength sampling:", cfg["res"])
 36print("Performing backpropagation.")
 37
 38# Set measurement parameters
 39# Compute scattered field from cylinder
 40radius = cfg["radius"]  # wavelengths
 41nmed = cfg["nmed"]
 42ncyl = cfg["ncyl"]
 43
 44lD = cfg["lD"]  # measurement distance in wavelengths
 45lC = cfg["lC"]  # displacement from center of image
 46size = cfg["size"]
 47res = cfg["res"]  # px/wavelengths
 48A = cfg["A"]  # number of projections
 49
 50x = np.arange(size) - size / 2.0
 51X, Y = np.meshgrid(x, x)
 52rad_px = radius * res
 53phantom = np.array(((Y - lC * res)**2 + X**2) < rad_px **
 54                   2, dtype=np.float) * (ncyl - nmed) + nmed
 55
 56u_sinR = odt.sinogram_as_rytov(sino / u0)
 57
 58# Rytov 200 projections
 59# remove 50 projections from total of 250 projections
 60remove200 = np.argsort(angles % .0002)[:50]
 61angles200 = np.delete(angles, remove200, axis=0)
 62u_sinR200 = np.delete(u_sinR, remove200, axis=0)
 63ph200 = unwrap.unwrap(np.angle(sino / u0))
 64ph200[remove200] = 0
 65
 66fR200 = odt.backpropagate_2d(u_sinR200, angles200, res, nmed, lD*res)
 67nR200 = odt.odt_to_ri(fR200, res, nmed)
 68fR200nw = odt.backpropagate_2d(u_sinR200, angles200, res, nmed, lD*res,
 69                               weight_angles=False)
 70nR200nw = odt.odt_to_ri(fR200nw, res, nmed)
 71
 72# Rytov 50 projections
 73remove50 = np.argsort(angles % .0002)[:200]
 74angles50 = np.delete(angles, remove50, axis=0)
 75u_sinR50 = np.delete(u_sinR, remove50, axis=0)
 76ph50 = unwrap.unwrap(np.angle(sino / u0))
 77ph50[remove50] = 0
 78
 79fR50 = odt.backpropagate_2d(u_sinR50, angles50, res, nmed, lD*res)
 80nR50 = odt.odt_to_ri(fR50, res, nmed)
 81fR50nw = odt.backpropagate_2d(u_sinR50, angles50, res, nmed, lD*res,
 82                              weight_angles=False)
 83nR50nw = odt.odt_to_ri(fR50nw, res, nmed)
 84
 85# prepare plot
 86kw_ri = {"vmin": 1.330,
 87         "vmax": 1.340}
 88
 89kw_ph = {"vmin": np.min(np.array([ph200, ph50])),
 90         "vmax": np.max(np.array([ph200, ph50])),
 91         "cmap": "coolwarm"}
 92
 93fig, axes = plt.subplots(2, 3, figsize=(8, 4))
 94axes = np.array(axes).flatten()
 95
 96phmap = axes[0].imshow(ph200, **kw_ph)
 97axes[0].set_title("Phase sinogram (200 proj.)")
 98
 99rimap = axes[1].imshow(nR200nw.real, **kw_ri)
100axes[1].set_title("RI without angular weights")
101
102axes[2].imshow(nR200.real, **kw_ri)
103axes[2].set_title("RI with angular weights")
104
105axes[3].imshow(ph50, **kw_ph)
106axes[3].set_title("Phase sinogram (50 proj.)")
107
108axes[4].imshow(nR50nw.real, **kw_ri)
109axes[4].set_title("RI without angular weights")
110
111axes[5].imshow(nR50.real, **kw_ri)
112axes[5].set_title("RI with angular weights")
113
114# color bars
115cbkwargs = {"fraction": 0.045,
116            "format": "%.3f"}
117plt.colorbar(phmap, ax=axes[0], **cbkwargs)
118plt.colorbar(phmap, ax=axes[3], **cbkwargs)
119plt.colorbar(rimap, ax=axes[1], **cbkwargs)
120plt.colorbar(rimap, ax=axes[2], **cbkwargs)
121plt.colorbar(rimap, ax=axes[5], **cbkwargs)
122plt.colorbar(rimap, ax=axes[4], **cbkwargs)
123
124plt.tight_layout()
125plt.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.

_images/backprop_from_mie_2d_incomplete_coverage.jpg

backprop_from_mie_2d_incomplete_coverage.py

  1import matplotlib.pylab as plt
  2import numpy as np
  3
  4import odtbrain as odt
  5
  6from example_helper import load_data
  7
  8sino, angles, cfg = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
  9                              f_angles="mie_angles.txt",
 10                              f_sino_real="sino_real.txt",
 11                              f_sino_imag="sino_imag.txt",
 12                              f_info="mie_info.txt")
 13A, size = sino.shape
 14
 15# background sinogram computed with Mie theory
 16# miefield.GetSinogramCylinderRotation(radius, nmed, nmed, lD, lC, size, A,res)
 17u0 = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 18               f_sino_imag="u0_imag.txt",
 19               f_sino_real="u0_real.txt")
 20# create 2d array
 21u0 = np.tile(u0, size).reshape(A, size).transpose()
 22
 23# background field necessary to compute initial born field
 24# u0_single = mie.GetFieldCylinder(radius, nmed, nmed, lD, size, res)
 25u0_single = load_data("mie_2d_noncentered_cylinder_A250_R2.zip",
 26                      f_sino_imag="u0_single_imag.txt",
 27                      f_sino_real="u0_single_real.txt")
 28
 29print("Example: Backpropagation from 2D FDTD simulations")
 30print("Refractive index of medium:", cfg["nmed"])
 31print("Measurement position from object center:", cfg["lD"])
 32print("Wavelength sampling:", cfg["res"])
 33print("Performing backpropagation.")
 34
 35# Set measurement parameters
 36# Compute scattered field from cylinder
 37radius = cfg["radius"]  # wavelengths
 38nmed = cfg["nmed"]
 39ncyl = cfg["ncyl"]
 40
 41lD = cfg["lD"]  # measurement distance in wavelengths
 42lC = cfg["lC"]  # displacement from center of image
 43size = cfg["size"]
 44res = cfg["res"]  # px/wavelengths
 45A = cfg["A"]  # number of projections
 46
 47x = np.arange(size) - size / 2.0
 48X, Y = np.meshgrid(x, x)
 49rad_px = radius * res
 50phantom = np.array(((Y - lC * res)**2 + X**2) < rad_px **
 51                   2, dtype=np.float) * (ncyl - nmed) + nmed
 52
 53u_sinR = odt.sinogram_as_rytov(sino / u0)
 54
 55# Rytov 100 projections evenly distributed
 56removeeven = np.argsort(angles % .002)[:150]
 57angleseven = np.delete(angles, removeeven, axis=0)
 58u_sinReven = np.delete(u_sinR, removeeven, axis=0)
 59pheven = odt.sinogram_as_radon(sino / u0)
 60pheven[removeeven] = 0
 61
 62fReven = odt.backpropagate_2d(u_sinReven, angleseven, res, nmed, lD * res)
 63nReven = odt.odt_to_ri(fReven, res, nmed)
 64fRevennw = odt.backpropagate_2d(
 65    u_sinReven, angleseven, res, nmed, lD * res, weight_angles=False)
 66nRevennw = odt.odt_to_ri(fRevennw, res, nmed)
 67
 68# Rytov 100 projections more than 180
 69removemiss = 249 - \
 70    np.concatenate((np.arange(100), 100 + np.arange(150)[::3]))
 71anglesmiss = np.delete(angles, removemiss, axis=0)
 72u_sinRmiss = np.delete(u_sinR, removemiss, axis=0)
 73phmiss = odt.sinogram_as_radon(sino / u0)
 74phmiss[removemiss] = 0
 75
 76fRmiss = odt.backpropagate_2d(u_sinRmiss, anglesmiss, res, nmed, lD * res)
 77nRmiss = odt.odt_to_ri(fRmiss, res, nmed)
 78fRmissnw = odt.backpropagate_2d(
 79    u_sinRmiss, anglesmiss, res, nmed, lD * res, weight_angles=False)
 80nRmissnw = odt.odt_to_ri(fRmissnw, res, nmed)
 81
 82# Rytov 100 projections less than 180
 83removebad = 249 - np.arange(150)
 84anglesbad = np.delete(angles, removebad, axis=0)
 85u_sinRbad = np.delete(u_sinR, removebad, axis=0)
 86phbad = odt.sinogram_as_radon(sino / u0)
 87phbad[removebad] = 0
 88
 89fRbad = odt.backpropagate_2d(u_sinRbad, anglesbad, res, nmed, lD * res)
 90nRbad = odt.odt_to_ri(fRbad, res, nmed)
 91fRbadnw = odt.backpropagate_2d(
 92    u_sinRbad, anglesbad, res, nmed, lD * res, weight_angles=False)
 93nRbadnw = odt.odt_to_ri(fRbadnw, res, nmed)
 94
 95# prepare plot
 96kw_ri = {"vmin": np.min(np.array([phantom, nRmiss.real, nReven.real])),
 97         "vmax": np.max(np.array([phantom, nRmiss.real, nReven.real]))}
 98
 99kw_ph = {"vmin": np.min(np.array([pheven, phmiss])),
100         "vmax": np.max(np.array([pheven, phmiss])),
101         "cmap": "coolwarm"}
102
103fig, axes = plt.subplots(3, 3, figsize=(8, 6.5))
104
105axes[0, 0].set_title("100% coverage ({} proj.)".format(angleseven.shape[0]))
106phmap = axes[0, 0].imshow(pheven, **kw_ph)
107
108axes[0, 1].set_title("RI without angular weights")
109rimap = axes[0, 1].imshow(nRevennw.real, **kw_ri)
110
111axes[0, 2].set_title("RI with angular weights")
112rimap = axes[0, 2].imshow(nReven.real, **kw_ri)
113
114axes[1, 0].set_title("60% coverage ({} proj.)".format(anglesmiss.shape[0]))
115axes[1, 0].imshow(phmiss, **kw_ph)
116
117axes[1, 1].set_title("RI without angular weights")
118axes[1, 1].imshow(nRmissnw.real, **kw_ri)
119
120axes[1, 2].set_title("RI with angular weights")
121axes[1, 2].imshow(nRmiss.real, **kw_ri)
122
123axes[2, 0].set_title("40% coverage ({} proj.)".format(anglesbad.shape[0]))
124axes[2, 0].imshow(phbad, **kw_ph)
125
126axes[2, 1].set_title("RI without angular weights")
127axes[2, 1].imshow(nRbadnw.real, **kw_ri)
128
129axes[2, 2].set_title("RI with angular weights")
130axes[2, 2].imshow(nRbad.real, **kw_ri)
131
132# color bars
133cbkwargs = {"fraction": 0.045,
134            "format": "%.3f"}
135plt.colorbar(phmap, ax=axes[0, 0], **cbkwargs)
136plt.colorbar(phmap, ax=axes[1, 0], **cbkwargs)
137plt.colorbar(phmap, ax=axes[2, 0], **cbkwargs)
138plt.colorbar(rimap, ax=axes[0, 1], **cbkwargs)
139plt.colorbar(rimap, ax=axes[1, 1], **cbkwargs)
140plt.colorbar(rimap, ax=axes[2, 1], **cbkwargs)
141plt.colorbar(rimap, ax=axes[0, 2], **cbkwargs)
142plt.colorbar(rimap, ax=axes[1, 2], **cbkwargs)
143plt.colorbar(rimap, ax=axes[2, 2], **cbkwargs)
144
145plt.tight_layout()
146plt.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.

_images/backprop_from_fdtd_2d.jpg

backprop_from_fdtd_2d.py

 1import matplotlib.pylab as plt
 2import numpy as np
 3import odtbrain as odt
 4
 5from example_helper import load_data
 6
 7
 8sino, angles, phantom, cfg = load_data("fdtd_2d_sino_A100_R13.zip",
 9                                       f_angles="fdtd_angles.txt",
10                                       f_sino_imag="fdtd_imag.txt",
11                                       f_sino_real="fdtd_real.txt",
12                                       f_info="fdtd_info.txt",
13                                       f_phantom="fdtd_phantom.txt",
14                                       )
15
16print("Example: Backpropagation from 2D FDTD simulations")
17print("Refractive index of medium:", cfg["nm"])
18print("Measurement position from object center:", cfg["lD"])
19print("Wavelength sampling:", cfg["res"])
20print("Performing backpropagation.")
21
22# Apply the Rytov approximation
23sino_rytov = odt.sinogram_as_rytov(sino)
24
25# perform backpropagation to obtain object function f
26f = odt.backpropagate_2d(uSin=sino_rytov,
27                         angles=angles,
28                         res=cfg["res"],
29                         nm=cfg["nm"],
30                         lD=cfg["lD"] * cfg["res"]
31                         )
32
33# compute refractive index n from object function
34n = odt.odt_to_ri(f, res=cfg["res"], nm=cfg["nm"])
35
36# compare phantom and reconstruction in plot
37fig, axes = plt.subplots(1, 3, figsize=(8, 2.8))
38
39axes[0].set_title("FDTD phantom")
40axes[0].imshow(phantom, vmin=phantom.min(), vmax=phantom.max())
41sino_phase = np.unwrap(np.angle(sino), axis=1)
42
43axes[1].set_title("phase sinogram")
44axes[1].imshow(sino_phase, vmin=sino_phase.min(), vmax=sino_phase.max(),
45               aspect=sino.shape[1] / sino.shape[0],
46               cmap="coolwarm")
47axes[1].set_xlabel("detector")
48axes[1].set_ylabel("angle [rad]")
49
50axes[2].set_title("reconstructed image")
51axes[2].imshow(n.real, vmin=phantom.min(), vmax=phantom.max())
52
53# set y ticks for sinogram
54labels = np.linspace(0, 2 * np.pi, len(axes[1].get_yticks()))
55labels = ["{:.2f}".format(i) for i in labels]
56axes[1].set_yticks(np.linspace(0, len(angles), len(labels)))
57axes[1].set_yticklabels(labels)
58
59plt.tight_layout()
60plt.show()