feat: add possibility to add (random) displacements

CHANGELOG

@@ -1,5 +1,7 @@
 0.2.0
  - feat: add capability to track progress of sinogram generation
+ - feat: add possibility to add (random) displacements for each
+   sinogram slice
  - fix: medium index not stored in exported sinogram file
 0.1.0
  - add basic example to docs 

cellsino/fluorescence.py

@@ -5,7 +5,7 @@
 
 class Fluorescence(object):
 
-    def __init__(self, phantom, grid_size, pixel_size):
+    def __init__(self, phantom, grid_size, pixel_size, displacement=(0, 0)):
         """Fluorescence projector
 
         Notes
@@ -20,6 +20,8 @@ def __init__(self, phantom, grid_size, pixel_size):
         gx, gy = grid_size
         #: center of the volume used
         self.center = np.array([gx, gy, 0]) / 2 - .5
+        self.center[0] += displacement[0]
+        self.center[1] += displacement[1]
 
     def project(self):
         fluor = np.zeros(self.grid_size, dtype=float)

cellsino/propagators/base_propagator.py

@@ -9,7 +9,8 @@
 class BasePropagator(object):
     __metaclass__ = abc.ABCMeta
 
-    def __init__(self, phantom, grid_size, pixel_size, wavelength):
+    def __init__(self, phantom, grid_size, pixel_size, wavelength,
+                 displacement=(0, 0)):
         """Base propagator
 
         Notes
@@ -25,20 +26,23 @@ def __init__(self, phantom, grid_size, pixel_size, wavelength):
         gx, gy = grid_size
         #: center of the volume used
         self.center = np.array([gx, gy, 0]) / 2 - .5
+        self.center[0] += displacement[0]
+        self.center[1] += displacement[1]
 
     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).field
-        qpifull = qpimage.QPImage(data=field,
-                                  which_data="field",
-                                  meta_data={
-                                      "wavelength": self.wavelength,
-                                      "pixel size": self.pixel_size,
-                                      "medium index": self.phantom.medium_index,
-                                      }
-                                  )
+        qpifull = qpimage.QPImage(
+                    data=field,
+                    which_data="field",
+                    meta_data={
+                        "wavelength": self.wavelength,
+                        "pixel size": self.pixel_size,
+                        "medium index": self.phantom.medium_index,
+                        }
+                    )
         return qpifull
 
     @abc.abstractmethod

cellsino/sinogram.py

@@ -19,14 +19,48 @@ def __init__(self, phantom, wavelength, pixel_size, grid_size):
         self.pixel_size = pixel_size
         self.grid_size = grid_size
 
-    def compute(self, angles, path=None, propagator="rytov",
-                count=None, max_count=None):
+    def compute(self, angles, displacements=None, propagator="rytov",
+                path=None, count=None, max_count=None):
+        """Compute sinogram data
+
+        Parameters
+        ----------
+        angles: 1d array of size N or int
+            If an array, defines the angular positions for each frame
+            of the rotation. If an int, defines the number of angles
+            of a steady rotation from 0 to 2π.
+        displacements: 2d ndarray of shape (N, 2) or float
+            A float value indicates the standard deviation of a
+            Gaussian noise distribution (using
+            :func:`numpy.random.normal`) in pixels.
+            A 2d array directly specifies the x-y-displacement for each
+            frame in pixels.
+        propagator: str
+            The propagator to use. Must be in
+            :data:`cellsino.propagators.available`.
+        path: str or pathlib.Path
+            If not None, the data will be written to this file and
+            a :class:`pathlib.Path` object will be returned.
+        count: multiprocessing.Value
+            May be used for tracking progress. At each computed
+            angle `count.value` is incremented by one.
+        max_count: multiprocessing.Value
+            May be used for tracking progress. This value is
+            incremented by `N`.
+        """
         if max_count is not None:
             max_count.value += angles.size
 
         if isinstance(angles, int):
             angles = np.linspace(0, 2*np.pi, angles, endpoint=False)
 
+        if displacements is None:
+            displacements = np.zeros((angles.shape[0], 2))
+        elif isinstance(displacements, float):
+            np.random.seed(47)  # for reproducibility
+            displacements = np.random.normal(scale=displacements,
+                                             size=(angles.shape[0], 2))
+
         if path:
             write = True
             path = pathlib.Path(path)
@@ -49,11 +83,13 @@ def compute(self, angles, path=None, propagator="rytov",
             pp = prop_dict[propagator](phantom=ph,
                                        grid_size=self.grid_size,
                                        pixel_size=self.pixel_size,
-                                       wavelength=self.wavelength)
+                                       wavelength=self.wavelength,
+                                       displacement=displacements[ii])
             qpi = pp.propagate()
             fluor = Fluorescence(phantom=ph,
                                  grid_size=self.grid_size,
-                                 pixel_size=self.pixel_size).project()
+                                 pixel_size=self.pixel_size,
+                                 displacement=displacements[ii]).project()
 
             if write:
                 with h5py.File(path, "a") as h5:

tests/test_displacement.py

@@ -0,0 +1,34 @@
+import numpy as np
+
+from cellsino import sinogram
+
+
+def test_displacement():
+    kw = {"phantom": "simple cell",
+          "wavelength": 550e-9,
+          "pixel_size": 0.7e-6,
+          "grid_size": (25, 25),
+          }
+    sino = sinogram.Sinogram(**kw)
+
+    angles = np.array([0, .5, .6])
+    field1, fluor1 = sino.compute(angles=angles,
+                                  propagator="projection")
+    field2, fluor2 = sino.compute(angles=angles,
+                                  displacements=np.ones((3, 2)),
+                                  propagator="projection")
+    field3 = np.roll(field1, (1, 1), axis=(1, 2))
+    fluor3 = np.roll(fluor1, (1, 1), axis=(1, 2))
+
+    assert np.allclose(field2, field3, rtol=0, atol=1e-13)
+    assert np.allclose(fluor2, fluor3, rtol=0, atol=1e-13)
+    assert not np.allclose(field2, field1, rtol=0, atol=1e-13)
+    assert not np.allclose(fluor2, fluor1, rtol=0, atol=1e-13)
+
+
+if __name__ == "__main__":
+    # Run all tests
+    loc = locals()
+    for key in list(loc.keys()):
+        if key.startswith("test_") and hasattr(loc[key], "__call__"):
+            loc[key]()