initial working state with example

CHANGELOG

@@ -1,2 +1,9 @@
+0.1.0
+ - drawing refractive index and fluorescence phantoms
+ - sinogram generation
+ - basic functionalities (sphere element):
+   - simple cell phantom
+   - propagators: Rytov and projection
+   - fluorescence projector
 0.0.1
  - initial setup

cellsino/elements/base_element.py

@@ -17,7 +17,7 @@ def __init__(self, object_index, medium_index, fl_brightness,
         fl_brightness: float
             Fluorescence brightness
         points: 2d ndarray
-            Coordinates of the element the element
+            Coordinates of the element [m]
 
         Notes
         -----
@@ -36,6 +36,16 @@ def __init__(self, object_index, medium_index, fl_brightness,
         #: the object).
         self.points = np.array(points)
 
+    def draw(self, grid_size, pixel_size):
+        ri = np.ones(grid_size, dtype=float) * self.medium_index
+        fl = np.zeros(grid_size, dtype=float)
+        for pp in self.points:
+            # ODTbrain convention
+            cy, cz, cx = np.array(pp/pixel_size, dtype=int)
+            ri[cx, cy, cz] = self.object_index
+            fl[cx, cy, cz] = self.fl_brightness
+        return ri, fl
+
     def transform(self, x=0, y=0, z=0, rot_main=0, rot_in_plane=0,
                   rot_perp_plane=0):
         """Rotate and translate self.points

cellsino/elements/el_sphere.py

@@ -29,17 +29,23 @@ def __init__(self, object_index, medium_index, fl_brightness,
                                      fl_brightness=fl_brightness,
                                      points=points)
 
-    def draw(self, pixel_size, grid_size):
-        ri = self.medium_index * np.ones(grid_size, dtype=float)
+    @property
+    def center(self):
+        return self.points[0]
+
+    def draw(self, grid_size, pixel_size):
+        ri = np.ones(grid_size, dtype=float) * self.medium_index
+        fl = np.zeros(grid_size, dtype=float)
         center = np.array(grid_size) / 2 - .5
         x = (np.arange(grid_size[0]) - center[0]) * pixel_size
         y = (np.arange(grid_size[1]) - center[1]) * pixel_size
         z = (np.arange(grid_size[2]) - center[2]) * pixel_size
 
-        xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")
-        cx, cy, cz = self.center
+        xx, yy, zz = np.meshgrid(x, y, z, indexing="ij", sparse=True)
+        cy, cz, cx = self.center
 
         volume = (cx-xx)**2 + (cy-yy)**2 + (cz-zz)**2
         inside = volume <= self.radius**2
         ri[inside] = self.object_index
-        return ri
+        fl[inside] = self.fl_brightness
+        return ri, fl

cellsino/fluorescence.py

@@ -0,0 +1,47 @@
+import numpy as np
+
+from .elements import Sphere
+
+
+class Fluorescence(object):
+
+    def __init__(self, phantom, grid_size, pixel_size):
+        """Fluorescence projector
+
+        Notes
+        -----
+        - The origin of the coordinate system is the center of the grid.
+          For even grid sizes, the origin is between two pixels. For odd
+          grid sizes the origin coincides with the center pixel.
+        """
+        self.phantom = phantom
+        self.pixel_size = pixel_size
+        self.grid_size = grid_size
+        gx, gy = grid_size
+        #: center of the volume used
+        self.center = np.array([gx, gy, 0]) / 2 - .5
+
+    def project(self):
+        fluor = np.zeros(self.grid_size, dtype=float)
+        for element in self.phantom:
+            if isinstance(element, Sphere):
+                fluor += self.project_sphere(element)
+
+        return fluor
+
+    def project_sphere(self, sphere):
+        center = self.center + sphere.center/self.pixel_size
+        # grid
+        x = np.arange(self.grid_size[0]).reshape(-1, 1)
+        y = np.arange(self.grid_size[1]).reshape(1, -1)
+        cx, cy, _ = center
+        # sphere location
+        rpx = sphere.radius / self.pixel_size
+        r = rpx**2 - (x - cx)**2 - (y - cy)**2
+        # distance
+        z = np.zeros_like(r)
+        rvalid = r > 0
+        z[rvalid] = 2 * np.sqrt(r[rvalid])
+
+        fluor = z * sphere.fl_brightness
+        return fluor

cellsino/phantoms/base_phantom.py

@@ -1,6 +1,10 @@
+import numpy as np
+
+
 class BasePhantom(object):
-    def __init__(self):
+    def __init__(self, medium_index):
         self.elements = []
+        self.medium_index = medium_index
 
     def __iter__(self):
         for el in self.elements:
@@ -9,9 +13,19 @@ def __iter__(self):
     def append(self, element):
         self.elements.append(element)
 
+    def draw(self, grid_size, pixel_size):
+        ri = np.ones(grid_size, dtype=float) * self.medium_index
+        fl = np.zeros(grid_size, dtype=float)
+
+        for el in self:
+            riel, flel = el.draw(grid_size, pixel_size)
+            ri += riel - el.medium_index
+            fl += flel
+        return ri, fl
+
     def transform(self, x=0, y=0, z=0, rot_main=0, rot_in_plane=0,
                   rot_perp_plane=0):
-        ph = BasePhantom()
+        ph = BasePhantom(medium_index=self.medium_index)
         for el in self:
             ph.append(el.transform(x=x,
                                    y=y,

cellsino/phantoms/ph_simple_cell.py

@@ -57,3 +57,4 @@ def __init__(self,
                       medium_index=medium_index,
                       )
         self.elements = [nleoi1, nleoi2, nuclus, nuclus_shell, cyto]
+        self.medium_index = medium_index

cellsino/propagators/base_propagator.py

@@ -9,7 +9,7 @@
 class BasePropagator(object):
     __metaclass__ = abc.ABCMeta
 
-    def __init__(self, phantom, wavelength, pixel_size, grid_size):
+    def __init__(self, phantom, grid_size, pixel_size, wavelength):
         """Base propagator
 
         Notes

cellsino/propagators/pp_projection.py

@@ -7,7 +7,7 @@ class Projection(BasePropagator):
     """Projection approximation"""
 
     def propagate_sphere(self, sphere):
-        center = self.center + sphere.points[0]/self.pixel_size
+        center = self.center + sphere.center/self.pixel_size
         qpi = qpsphere.models.projection(radius=sphere.radius,
                                          sphere_index=sphere.object_index,
                                          medium_index=sphere.medium_index,

cellsino/sinogram.py

@@ -4,6 +4,7 @@
 import numpy as np
 import qpimage
 
+from .fluorescence import Fluorescence
 from .propagators import dictionary as pp_dict
 
 
@@ -35,22 +36,27 @@ def compute(self, angles, path=None, propagator="rytov"):
                                   dtype=float)
 
         for ii, ang in enumerate(angles):
-            print(ii)
             ph = self.phantom.transform(rot_main=ang)
 
             pp = pp_dict[propagator](phantom=ph,
-                                     wavelength=self.wavelength,
+                                     grid_size=self.grid_size,
                                      pixel_size=self.pixel_size,
-                                     grid_size=self.grid_size)
-
+                                     wavelength=self.wavelength)
             qpi = pp.propagate()
+            fluor = Fluorescence(phantom=ph,
+                                 grid_size=self.grid_size,
+                                 pixel_size=self.pixel_size).project()
 
             if write:
                 with h5py.File(path, "a") as h5:
                     qps = qpimage.QPSeries(h5file=h5.require_group("qpseries"))
                     qps.add_qpimage(qpi)
+
+                    h5fl = h5.require_group("fluorescence")
+                    h5fl.create_dataset("{}".format(ii), data=fluor)
             else:
                 sino_fields[ii] = qpi.field
+                sino_fluor[ii] = fluor
 
         if write:
             return path

examples/generate_example_images.py

@@ -0,0 +1,33 @@
+import os
+import os.path as op
+import sys
+
+import matplotlib.pylab as plt
+
+thisdir = op.dirname(op.abspath(__file__))
+sys.path.insert(0, op.dirname(thisdir))
+
+DPI = 80
+
+
+if __name__ == "__main__":
+    # Do not display example plots
+    plt.show = lambda: None
+    files = os.listdir(thisdir)
+    files = [f for f in files if f.endswith(".py")]
+    files = [f for f in files if not f == op.basename(__file__)]
+    files = sorted([op.join(thisdir, f) for f in files])
+
+    for f in files:
+        fname = f[:-3] + ".png"
+        if not op.exists(fname):
+            exec_str = open(f).read()
+            if exec_str.count("plt.show()"):
+                exec(exec_str)
+                plt.savefig(fname, dpi=DPI)
+                print("Image created: '{}'".format(fname))
+            else:
+                print("No image: '{}'".format(fname))
+        else:
+            print("Image skipped (already exists): '{}'".format(fname))
+        plt.close()

examples/simple_cell.py

@@ -0,0 +1,100 @@
+"""Reconstruction of a simple cell phantom
+
+This example uses :ref:`ODTbrain<odtbrain:index>` and
+:ref:`radontea<radontea:index>` for the reconstruction of the refractive
+index and the fluorescence of the simple cell phantom.
+The reconstruction is compared to the ground truth of the cell phantom.
+"""
+import cellsino
+import matplotlib.pylab as plt
+import numpy as np
+import odtbrain as odt
+import radontea as rt
+
+
+# number of sinogram angles
+num_ang = 160
+# sinogram acquisition angles
+angles = np.linspace(0, 2*np.pi, num_ang, endpoint=False)
+# detector grid size
+grid_size = (250, 250)
+# vacuum wavelength [m]
+wavelength = 550e-9
+# pixel size [m]
+pixel_size = 0.08e-6
+# refractive index of the surrounding medium
+medium_index = 1.335
+
+# initialize cell phantom
+phantom = cellsino.phantoms.SimpleCell()
+
+# initialize sinogram with geometric parameters
+sino = cellsino.Sinogram(phantom=phantom,
+                         wavelength=wavelength,
+                         pixel_size=pixel_size,
+                         grid_size=grid_size)
+
+# compute sinogram (field according to Rytov approximation and fluorescence)
+sino_field, sino_fluor = sino.compute(angles=angles, propagator="rytov")
+
+# reconstruction of refractive index
+sino_rytov = odt.sinogram_as_rytov(sino_field)
+potential = odt.backpropagate_3d(uSin=sino_rytov,
+                                 angles=angles,
+                                 res=wavelength/pixel_size,
+                                 nm=medium_index)
+ri = odt.odt_to_ri(f=potential,
+                   res=wavelength/pixel_size,
+                   nm=medium_index)
+
+# reconstruction of fluorescence
+fl = rt.backproject_3d(sinogram=sino_fluor,
+                       angles=angles)
+
+# reference for comparison
+rimod, flmod = phantom.draw(grid_size=ri.shape,
+                            pixel_size=pixel_size)
+
+# plotting
+idx = 150
+
+plt.figure(figsize=(7, 5.5))
+
+plotkwri = {"vmax": ri.real.max(),
+            "vmin": ri.real.min(),
+            "interpolation": "none",
+            "cmap": "viridis",
+            }
+
+plotkwfl = {"vmax": fl.max(),
+            "vmin": 0,
+            "interpolation": "none",
+            "cmap": "Greens_r",
+            }
+
+ax1 = plt.subplot(221, title="refractive index (ground_truth)")
+mapper = ax1.imshow(rimod[idx].real, **plotkwri)
+plt.colorbar(mappable=mapper, ax=ax1)
+
+ax2 = plt.subplot(222, title="refractive index (reconstruction)")
+mapper = ax2.imshow(ri[idx].real, **plotkwri)
+plt.colorbar(mappable=mapper, ax=ax2)
+
+ax3 = plt.subplot(223, title="fluorescence (ground truth)")
+mapper = ax3.imshow(flmod[idx], **plotkwfl)
+plt.colorbar(mappable=mapper, ax=ax3)
+
+ax4 = plt.subplot(224, title="fluorescence (reconstruction)")
+mapper = ax4.imshow(fl[idx], **plotkwfl)
+plt.colorbar(mappable=mapper, ax=ax4)
+
+for ax in [ax1, ax2, ax3, ax4]:
+    ax.grid(color="w", alpha=.5)
+    ax.set_xticks(np.arange(0, grid_size[0], 50))
+    ax.set_yticks(np.arange(0, grid_size[0], 50))
+    ax.set_xticklabels([])
+    ax.set_yticklabels([])
+
+plt.tight_layout()
+
+plt.show()