add sinogram generation, some refactoring, add basic functionalities

cellsino/__init__.py

@@ -1,3 +1,3 @@
 from ._version import version as __version__  # noqa: F401
-
+from . import phantoms  # noqa: F401
 from .sinogram import Sinogram  # noqa: F401

cellsino/elements/base_element.py

@@ -51,25 +51,29 @@ def transform(self, x=0, y=0, z=0, rot_main=0, rot_in_plane=0,
           Third, the points are rotated about the z-axis (``rot_in_plane``,
           within the imaging plane).
         """
+        # In ODTbrain convention, x -> -y and y -> x.
+        alpha = -rot_main
+        beta = rot_perp_plane
+        gamma = rot_in_plane
         Rx = np.array([
-                      [1,                      0,                       0],
-                      [0, np.cos(rot_perp_plane), -np.sin(rot_perp_plane)],
-                      [0, np.sin(rot_perp_plane),  np.cos(rot_perp_plane)],
+                      [1,             0,              0],
+                      [0, np.cos(alpha), -np.sin(alpha)],
+                      [0, np.sin(alpha),  np.cos(alpha)],
                       ])
 
         Ry = np.array([
-                      [np.cos(rot_main),  0, np.sin(rot_main)],
-                      [0,                 1,                0],
-                      [-np.sin(rot_main), 0, np.cos(rot_main)],
+                      [np.cos(beta),  0, np.sin(beta)],
+                      [0,             1,            0],
+                      [-np.sin(beta), 0, np.cos(beta)],
                       ])
 
         Rz = np.array([
-                      [np.cos(rot_in_plane), -np.sin(rot_in_plane), 0],
-                      [np.sin(rot_in_plane),  np.cos(rot_in_plane), 0],
-                      [0,                                        0, 1]
+                      [np.cos(gamma), -np.sin(gamma), 0],
+                      [np.sin(gamma),  np.cos(gamma), 0],
+                      [0,          0,                 1],
                       ])
 
-        R = np.dot(np.dot(Rx, Rz), Ry)
+        R = np.dot(np.dot(Ry, Rz), Rx)
         rotated = np.dot(R, self.points.T).T
         rotated_pad = np.pad(rotated, ((0, 0), (0, 1)), mode="constant",
                              constant_values=1)

cellsino/elements/el_sphere.py

@@ -28,3 +28,18 @@ def __init__(self, object_index, medium_index, fl_brightness,
                                      medium_index=medium_index,
                                      fl_brightness=fl_brightness,
                                      points=points)
+
+    def draw(self, pixel_size, grid_size):
+        ri = self.medium_index * np.ones(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
+
+        volume = (cx-xx)**2 + (cy-yy)**2 + (cz-zz)**2
+        inside = volume <= self.radius**2
+        ri[inside] = self.object_index
+        return ri

cellsino/phantoms/base_phantom.py

@@ -1,6 +1,6 @@
 class BasePhantom(object):
-    def __init__(self, elements=[]):
-        self.elements = elements
+    def __init__(self):
+        self.elements = []
 
     def __iter__(self):
         for el in self.elements:

cellsino/phantoms/ph_simple_cell.py

@@ -20,37 +20,40 @@ def __init__(self,
             - 2 small nucleoli (full fluorescence)
         """
         nleoi1 = Sphere(radius=1.5e-6,
-                        position=(-.25, 2, 2),
-                        fl_intensity=nucleoli_fl-nucleus_fl,
-                        n_object=medium_index+(nucleoli_index-nucleus_index),
+                        center=(-.25e-6, 2e-6, 2e-6),
+                        fl_brightness=nucleoli_fl-nucleus_fl,
+                        object_index=medium_index +
+                        (nucleoli_index-nucleus_index),
                         medium_index=medium_index,
                         )
 
         nleoi2 = Sphere(radius=1.5e-6,
-                        position=(+.25, -1, 2),
-                        fl_intensity=nucleoli_fl-nucleus_fl,
-                        n_object=medium_index+(nucleoli_index-nucleus_index),
+                        center=(+.25e-6, -1e-6, 2e-6),
+                        fl_brightness=nucleoli_fl-nucleus_fl,
+                        object_index=medium_index +
+                        (nucleoli_index-nucleus_index),
                         medium_index=medium_index,
                         )
 
         nuclus = Sphere(radius=4e-6,
-                        position=(0, 1, 1),
-                        fl_intensity=nucleus_shell_fl,
-                        n_object=medium_index+(nucleus_index-cytoplasm_index),
+                        center=(0e-6, 1e-6, 1e-6),
+                        fl_brightness=nucleus_shell_fl,
+                        object_index=medium_index +
+                        (nucleus_index-cytoplasm_index),
                         medium_index=medium_index,
                         )
 
         nuclus_shell = Sphere(radius=3.8e-6,
-                              position=(0, 1, 1),
-                              fl_intensity=nucleus_fl-nucleus_shell_fl,
-                              n_object=medium_index,
+                              center=(0e-6, 1e-6, 1e-6),
+                              fl_brightness=nucleus_fl-nucleus_shell_fl,
+                              object_index=medium_index,
                               medium_index=medium_index,
                               )
 
         cyto = Sphere(radius=5.5e-6,
-                      position=(0, 0, 0),
-                      fl_intensity=cytoplasm_fl,
-                      n_object=cytoplasm_index,
+                      center=(0, 0, 0),
+                      fl_brightness=cytoplasm_fl,
+                      object_index=cytoplasm_index,
                       medium_index=medium_index,
                       )
         self.elements = [nleoi1, nleoi2, nuclus, nuclus_shell, cyto]

cellsino/propagators/__init__.py

@@ -1,5 +1,9 @@
-from .rytov import Rytov
+from .pp_rytov import Rytov
+from .pp_projection import Projection
 
-available = [Rytov]
+dictionary = {
+    "projection": Projection,
+    "rytov": Rytov,
+}
 
-dictionary = {"rytov": Rytov}
+available = sorted(dictionary.keys())

cellsino/propagators/base_propagator.py

@@ -1,6 +1,7 @@
 import abc
 
 import numpy as np
+import qpimage
 
 from ..elements import Sphere
 
@@ -21,14 +22,22 @@ def __init__(self, phantom, wavelength, pixel_size, grid_size):
         self.wavelength = wavelength
         self.pixel_size = pixel_size
         self.grid_size = grid_size
-        self.center = np.array(self.grid_size) / 2 - .5
+        gx, gy = grid_size
+        #: center of the volume used
+        self.center = np.array([gx, gy, 0]) / 2 - .5
 
     def propagate(self):
         field = np.ones(self.grid_size, dtype=np.complex256)
         for element in self.phantom:
             if isinstance(element, Sphere):
-                field *= self.propagate_sphere(element)
-        return field
+                field *= self.propagate_sphere(element).field
+        qpifull = qpimage.QPImage(data=field,
+                                  which_data="field",
+                                  meta_data={"wavelength": self.wavelength,
+                                             "pixel size": self.pixel_size,
+                                             }
+                                  )
+        return qpifull
 
     @abc.abstractmethod
     def propagate_sphere(self, sphere):

cellsino/propagators/pp_projection.py

@@ -0,0 +1,19 @@
+import qpsphere
+
+from .base_propagator import BasePropagator
+
+
+class Projection(BasePropagator):
+    """Projection approximation"""
+
+    def propagate_sphere(self, sphere):
+        center = self.center + sphere.points[0]/self.pixel_size
+        qpi = qpsphere.models.projection(radius=sphere.radius,
+                                         sphere_index=sphere.object_index,
+                                         medium_index=sphere.medium_index,
+                                         wavelength=self.wavelength,
+                                         pixel_size=self.pixel_size,
+                                         grid_size=self.grid_size,
+                                         center=center[:2],
+                                         )
+        return qpi

cellsino/propagators/rytov.py → cellsino/propagators/pp_rytov.py

@@ -1,3 +1,4 @@
+import numpy as np
 import qpsphere
 
 from .base_propagator import BasePropagator
@@ -6,13 +7,19 @@
 class Rytov(BasePropagator):
     """Rytov approximation"""
 
-    def propagate_sphere(self, sphere):
+    def propagate_sphere(self, sphere, grid_sampling=150):
         center = self.center + sphere.points[0]/self.pixel_size
+        # speed up computation for smaller spheres on large grid
+        n = np.max(self.grid_size) * self.pixel_size / sphere.radius
+        radius_sampling = grid_sampling / n
         qpi = qpsphere.models.rytov(radius=sphere.radius,
                                     sphere_index=sphere.object_index,
                                     medium_index=sphere.medium_index,
                                     wavelength=self.wavelength,
                                     pixel_size=self.pixel_size,
                                     grid_size=self.grid_size,
-                                    center=center)
-        return qpi.field
+                                    center=center[:2],
+                                    focus=-center[2]*self.pixel_size,
+                                    radius_sampling=radius_sampling,
+                                    )
+        return qpi

cellsino/sinogram.py

@@ -16,7 +16,7 @@ def __init__(self, phantom, wavelength, pixel_size, grid_size):
 
     def compute(self, angles, path=None, propagator="rytov"):
         if isinstance(angles, int):
-            angles = np.linspace(0, 2*np.pi, endpoint=False)
+            angles = np.linspace(0, 2*np.pi, angles, endpoint=False)
 
         if path:
             write = True
@@ -34,28 +34,23 @@ def compute(self, angles, path=None, propagator="rytov"):
                                    self.grid_size[1]),
                                   dtype=float)
 
-        qpi_meta_data = {"wavelength": self.wavelength,
-                         "pixel size": self.pixel_size,
-                         }
-
         for ii, ang in enumerate(angles):
+            print(ii)
             ph = self.phantom.transform(rot_main=ang)
+
             pp = pp_dict[propagator](phantom=ph,
                                      wavelength=self.wavelength,
                                      pixel_size=self.pixel_size,
                                      grid_size=self.grid_size)
 
-            field = pp.propagate()
+            qpi = pp.propagate()
 
             if write:
                 with h5py.File(path, "a") as h5:
                     qps = qpimage.QPSeries(h5file=h5.require_group("qpseries"))
-                    qpi = qpimage.QPImage(data=field,
-                                          which_data="field",
-                                          meta_data=qpi_meta_data)
                     qps.add_qpimage(qpi)
             else:
-                sino_fields[ii] = field
+                sino_fields[ii] = qpi.field
 
         if write:
             return path

setup.py

@@ -26,7 +26,7 @@
     long_description=open('README.rst').read() if exists('README.rst') else '',
     install_requires=["h5py>=2.7.0",
                       "qpimage",
-                      "qpsphere",
+                      "qpsphere>=0.5.0",
                       "numpy>=1.12.0",
                       ],
     setup_requires=['pytest-runner'],