ENH: Add FORBILD phantom in 2d

odlgroup · adler-j · Nov 18, 2016 · Oct 31, 2016 · Nov 18, 2016 · Oct 31, 2016
commit c6952fafe218a9e8a6260540d0d14b9d5fb1e5f3
diff --git a/odl/phantom/transmission.py b/odl/phantom/transmission.py
@@ -22,7 +22,9 @@
 from future import standard_library
 standard_library.install_aliases()
 
+from odl.discr import DiscreteLp
 from odl.phantom.geometric import ellipse_phantom
+import numpy as np
 
 
 __all__ = ('shepp_logan_ellipses', 'shepp_logan')
@@ -128,6 +130,7 @@ def shepp_logan(space, modified=False):
 
     See Also
     --------
+    forbild : Similar phantom but with more complexity. Only supports 2d.
     shepp_logan_ellipses : Get the parameters that define this phantom
     odl.phantom.geometric.ellipse_phantom :
         Function for creating arbitrary ellipse phantoms
@@ -141,6 +144,183 @@ def shepp_logan(space, modified=False):
     return ellipse_phantom(space, ellipses)
 
 
+def _analytical_forbild_phantom(resolution, ear):
+    """Analytical description of forbild phantom.
+
+    Parameters
+    ----------
+    resolution : `bool`
+        If true, insers a small resolution insert to the left.
+    ear : `bool`
+        If true, insers a ear-like structure to the right.
+    """
+    sha = 0.2 * np.sqrt(3)
+    y016b = -14.294530834372887
+    a16b = 0.443194085308632
+    b16b = 3.892760834372886
+
+    E = [[-4.7, 4.3, 1.79989, 1.79989, 0, 0.010, 0],  # 1
+         [4.7, 4.3, 1.79989, 1.79989, 0, 0.010, 0],  # 2
+         [-1.08, -9, 0.4, 0.4, 0, 0.0025, 0],  # 3
+         [1.08, -9, 0.4, 0.4, 0, -0.0025, 0],  # 4
+         [0, 0, 9.6, 12, 0, 1.800, 0],  # 5
+         [0, 8.4, 1.8, 3.0, 0, -1.050, 0],  # 7
+         [1.9, 5.4, 0.41633, 1.17425, -31.07698, 0.750, 0],  # 8
+         [-1.9, 5.4, 0.41633, 1.17425, 31.07698, 0.750, 0],  # 9
+         [-4.3, 6.8, 1.8, 0.24, -30, 0.750, 0],  # 10
+         [4.3, 6.8, 1.8, 0.24, 30, 0.750, 0],  # 11
+         [0, -3.6, 1.8, 3.6, 0, -0.005, 0],  # 12
+         [6.39395, -6.39395, 1.2, 0.42, 58.1, 0.005, 0],  # 13
+         [0, 3.6, 2, 2, 0, 0.750, 4],  # 14
+         [0, 9.6, 1.8, 3.0, 0, 1.800, 4],  # 15
+         [0, 0, 9.0, 11.4, 0, 0.750, 3],  # 16a
+         [0, y016b, a16b, b16b, 0, 0.750, 1],  # 16b
+         [0, 0, 9.0, 11.4, 0, -0.750, ear],  # 6
+         [9.1, 0, 4.2, 1.8, 0, 0.750, 1]]  # R_ear
+    E = np.array(E)
+
+    # generate the air cavities in the right ear
+    cavity1 = np.arange(8.8, 5.6, -0.4)[:, None]
+    cavity2 = np.zeros([9, 1])
+    cavity3_7 = np.ones([53, 1]) * [0.15, 0.15, 0, -1.800, 0]
+
+    for j in range(1, 4):
+        kj = 8 - 2 * int(np.floor(j / 3))
+        dj = 0.2 * int(np.mod(j, 2))
+
+        cavity1 = np.vstack((cavity1,
+                             cavity1[0:kj] - dj,
+                             cavity1[0:kj] - dj))
+        cavity2 = np.vstack((cavity2,
+                             j * sha * np.ones([kj, 1]),
+                             -j * sha * np.ones([kj, 1])))
+
+    E_cavity = np.hstack((cavity1, cavity2, cavity3_7))
+
+    # generate the left ear (resolution pattern)
+    x0 = -7.0
+    y0 = -1.0
+    d0_xy = 0.04
+
+    d_xy = [0.0357, 0.0312, 0.0278, 0.0250]
+    ab = 0.5 * np.ones([5, 1]) * d_xy
+    ab = ab.T.ravel()[:, None] * np.ones([1, 4])
+    abr = ab.T.ravel()[:, None]
+
+    leftear4_7 = np.hstack([abr, abr, np.ones([80, 1]) * [0, 0.75, 0]])
+
+    x00 = np.zeros([0, 1])
+    y00 = np.zeros([0, 1])
+    for i in range(1, 5):
+        y00 = np.vstack((y00,
+                         (y0 + np.arange(0, 5) * 2 * d_xy[i - 1])[:, None]))
+        x00 = np.vstack((x00,
+                         (x0 + 2 * (i - 1) * d0_xy) * np.ones([5, 1])))
+
+    x00 = x00 * np.ones([1, 4])
+    x00 = x00.T.ravel()[:, None]
+    y00 = np.vstack([y00, y00 + 12 * d0_xy,
+                     y00 + 24 * d0_xy, y00 + 36 * d0_xy])
+
+    leftear = np.hstack([x00, y00, leftear4_7])
+    C = [[1.2, 1.2, 0.27884, 0.27884, 0.60687, 0.60687, 0.2,
+          0.2, -2.605, -2.605, -10.71177, y016b + 10.71177, 8.88740, -0.21260],
+         [0, 180, 90, 270, 90, 270, 0,
+          180, 15, 165, 90, 270, 0, 0]]
+    C = np.array(C)
+
+    if not resolution and not ear:
+        phantomE = E[:17, :]
+        phantomC = C[:, :12]
+    elif not resolution and ear:
+        phantomE = np.vstack([E, E_cavity])
+        phantomC = C
+    elif resolution and not ear:
+        phantomE = np.vstack([leftear, E[:17, :]])
+        phantomC = C[:, :12]
+    else:
+        phantomE = np.vstack([leftear, E, E_cavity])
+        phantomC = C
+
+    return phantomE, phantomC
+
+
+def forbild(space, resolution=False, ear=True):
+    """Standard `FORBILD phantom` in 2 dimensions.
+
+    The forbild phantom is intended for testing CT-algorithms and is intended
+    to be similar to a human head.
+
+    Parameters
+    ----------
+    space : `DiscreteLp`
+        The space in which the phantom should be corrected. Needs to be two-
+        dimensional.
+    resolution : `bool`, optional
+        If true, insers a small resolution insert to the left.
+    ear : `bool`, optional
+        If true, insers a ear-like structure to the right.
+
+    See Also
+    --------
+    shepp_logan : A more simple phantom with similar uses, also works in 3d.
+
+    References
+    ----------
+    .. _FORBILD phantom: www.imp.uni-erlangen.de/phantoms/head/head.html
+    .. _algorithm: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3426508/
+    """
+    def transposeravel(arr):
+        """Implement matlabs ``transpose(arr(:))``."""
+        return arr.T.ravel()
+
+    if not isinstance(space, DiscreteLp):
+        raise TypeError('`space` must be a `DiscreteLp`')
+    if not space.ndim == 2:
+        raise TypeError('`space` must be two dimensional')
+
+    # Create analytic description of phantom
+    phantomE, phantomC = _analytical_forbild_phantom(resolution, ear)
+
+    # Rescale points to the default grid.
+    # The forbild phantom is defined on [-12.8, 12.8] x [-12.8, 12.8]
+    xcoord, ycoord = space.points().T
+    xcoord = (xcoord - np.min(xcoord)) / (np.max(xcoord) - np.min(xcoord))
+    xcoord = 25.8 * xcoord - 12.8
+    ycoord = (ycoord - np.min(ycoord)) / (np.max(ycoord) - np.min(ycoord))
+    ycoord = 25.8 * ycoord - 12.8
+
+    # Compute the phantom values in each voxel
+    image = np.zeros(space.size)
+    nclipinfo = 0
+    for k in range(phantomE.shape[0]):
+        # Handle elliptic bounds
+        Vx0 = np.array([transposeravel(xcoord) - phantomE[k, 0],
+                        transposeravel(ycoord) - phantomE[k, 1]])
+        D = np.array([[1 / phantomE[k, 2], 0],
+                      [0, 1 / phantomE[k, 3]]])
+        phi = phantomE[k, 4] * np.pi / 180
+        Q = np.array([[np.cos(phi), np.sin(phi)],
+                      [-np.sin(phi), np.cos(phi)]])
+        f = phantomE[k, 5]
+        nclip = int(phantomE[k, 6])
+        equation1 = np.sum(((D.dot(Q)).dot(Vx0))**2, axis=0)
+        i = (equation1 <= 1.0)
+
+        # Handle clipping surfaces
+        if nclip > 0:
+            for j in range(1, nclip + 1):
+                d = phantomC[0, nclipinfo]
+                psi = phantomC[1, nclipinfo] * np.pi / 180.0
+                equation2 = np.array([np.cos(psi), np.sin(psi)]).dot(Vx0)
+                i *= (equation2 < d)
+                nclipinfo += 1
+
+        image[i] += f
+
+    return space.element(image)
+
+
 if __name__ == '__main__':
     # Show the phantoms
     import odl
@@ -149,6 +329,7 @@ def shepp_logan(space, modified=False):
     discr = odl.uniform_discr([-1, -1], [1, 1], [1000, 1000])
     shepp_logan(discr, modified=True).show('shepp_logan 2d modified=True')
     shepp_logan(discr, modified=False).show('shepp_logan 2d modified=False')
+    forbild(discr).show('forbild 2d', clim=[1.035, 1.065])
 
     # 3D
     discr = odl.uniform_discr([-1, -1, -1], [1, 1, 1], [300, 300, 300])