def _ellipsoid_phantom_3d(space, ellipsoids):
    """Create an ellipsoid phantom in 3d space.

    Parameters
    ----------
    space : `DiscreteLp`
        Space in which the phantom should be generated. If ``space.shape`` is
        1 in an axis, a corresponding slice of the phantom is created
        (instead of squashing the whole phantom into the slice).
    ellipsoids : list of lists
        Each row should contain the entries ::

            'value',
            'axis_1', 'axis_2', 'axis_3',
            'center_x', 'center_y', 'center_z',
            'rotation_phi', 'rotation_theta', 'rotation_psi'

        The provided ellipsoids need to be specified relative to the
        reference cube ``[-1, -1, -1] x [1, 1, 1]``. Angles are to be given
        in radians.

    Returns
    -------
    phantom : ``space`` element
        3D ellipsoid phantom in ``space``.

    See Also
    --------
    shepp_logan : The typical use-case for this function.
    """
    # Blank volume
    p = np.zeros(space.shape, dtype=space.dtype)

    minp = space.grid.min_pt
    maxp = space.grid.max_pt

    # Create the pixel grid
    grid_in = space.grid.meshgrid

    # Move points to [-1, 1]
    grid = []
    for i in range(3):
        mean_i = (minp[i] + maxp[i]) / 2.0
        # Where space.shape = 1, we have minp = maxp, so we set diff_i = 1
        # to avoid division by zero. Effectively, this allows constructing
        # a slice of a 3D phantom.
        diff_i = (maxp[i] - minp[i]) / 2.0 or 1.0
        grid.append((grid_in[i] - mean_i) / diff_i)

    for ellip in ellipsoids:
        assert len(ellip) == 10

        intensity = ellip[0]
        a_squared = ellip[1] ** 2
        b_squared = ellip[2] ** 2
        c_squared = ellip[3] ** 2
        x0 = ellip[4]
        y0 = ellip[5]
        z0 = ellip[6]
        phi = ellip[7]
        theta = ellip[8]
        psi = ellip[9]

        scales = [1 / a_squared, 1 / b_squared, 1 / c_squared]
        center = (np.array([x0, y0, z0]) + 1.0) / 2.0

        # Create the offset x,y and z values for the grid
        if any([phi, theta, psi]):
            # Rotate the points to the expected coordinate system.
            cphi = np.cos(phi)
            sphi = np.sin(phi)
            ctheta = np.cos(theta)
            stheta = np.sin(theta)
            cpsi = np.cos(psi)
            spsi = np.sin(psi)

            mat = np.array([[cpsi * cphi - ctheta * sphi * spsi,
                             cpsi * sphi + ctheta * cphi * spsi,
                             spsi * stheta],
                            [-spsi * cphi - ctheta * sphi * cpsi,
                             -spsi * sphi + ctheta * cphi * cpsi,
                             cpsi * stheta],
                            [stheta * sphi,
                             -stheta * cphi,
                             ctheta]])

            # Calculate the points that could possibly be inside the volume
            # Since the points are rotated, we cannot do anything directional
            # without more logic
            max_radius = np.sqrt(
                np.abs(mat).dot([a_squared, b_squared, c_squared]))
            idx, shapes = _getshapes_3d(center, max_radius, space.shape)

            subgrid = [g[idi] for g, idi in zip(grid, shapes)]
            offset_points = [vec * (xi - x0i)[..., None]
                             for xi, vec, x0i in zip(subgrid,
                                                     mat.T,
                                                     [x0, y0, z0])]
            rotated = offset_points[0] + offset_points[1] + offset_points[2]
            np.square(rotated, out=rotated)
            radius = np.dot(rotated, scales)
        else:
            # Calculate the points that could possibly be inside the volume
            max_radius = np.sqrt([a_squared, b_squared, c_squared])
            idx, shapes = _getshapes_3d(center, max_radius, space.shape)

            subgrid = [g[idi] for g, idi in zip(grid, shapes)]
            squared_dist = [ai * (xi - x0i) ** 2
                            for xi, ai, x0i in zip(subgrid,
                                                   scales,
                                                   [x0, y0, z0])]

            # Parentheses to get best order for broadcasting
            radius = squared_dist[0] + (squared_dist[1] + squared_dist[2])

        # Find the points within the ellipse
        inside = radius <= 1

        # Add the ellipse intensity to those points
        p[idx][inside] += intensity

    return space.element(p)
def _ellipsoid_phantom_3d ( space , ellipsoids ) : p = np . zeros ( space . shape , dtype = space . dtype ) minp = space . grid . min_pt maxp = space . grid . max_pt grid_in = space . grid . meshgrid grid = [ ] for i in range ( 3 ) : mean_i = ( minp [ i ] + maxp [ i ] ) / 2.0 diff_i = ( maxp [ i ] - minp [ i ] ) / 2.0 or 1.0 grid . append ( ( grid_in [ i ] - mean_i ) / diff_i ) for ellip in ellipsoids : assert len ( ellip ) == 10 intensity = ellip [ 0 ] a_squared = ellip [ 1 ] ** 2 b_squared = ellip [ 2 ] ** 2 c_squared = ellip [ 3 ] ** 2 x0 = ellip [ 4 ] y0 = ellip [ 5 ] z0 = ellip [ 6 ] phi = ellip [ 7 ] theta = ellip [ 8 ] psi = ellip [ 9 ] scales = [ 1 / a_squared , 1 / b_squared , 1 / c_squared ] center = ( np . array ( [ x0 , y0 , z0 ] ) + 1.0 ) / 2.0 if any ( [ phi , theta , psi ] ) : cphi = np . cos ( phi ) sphi = np . sin ( phi ) ctheta = np . cos ( theta ) stheta = np . sin ( theta ) cpsi = np . cos ( psi ) spsi = np . sin ( psi ) mat = np . array ( [ [ cpsi * cphi - ctheta * sphi * spsi , cpsi * sphi + ctheta * cphi * spsi , spsi * stheta ] , [ - spsi * cphi - ctheta * sphi * cpsi , - spsi * sphi + ctheta * cphi * cpsi , cpsi * stheta ] , [ stheta * sphi , - stheta * cphi , ctheta ] ] ) max_radius = np . sqrt ( np . abs ( mat ) . dot ( [ a_squared , b_squared , c_squared ] ) ) idx , shapes = _getshapes_3d ( center , max_radius , space . shape ) subgrid = [ g [ idi ] for g , idi in zip ( grid , shapes ) ] offset_points = [ vec * ( xi - x0i ) [ ... , None ] for xi , vec , x0i in zip ( subgrid , mat . T , [ x0 , y0 , z0 ] ) ] rotated = offset_points [ 0 ] + offset_points [ 1 ] + offset_points [ 2 ] np . square ( rotated , out = rotated ) radius = np . dot ( rotated , scales ) else : max_radius = np . sqrt ( [ a_squared , b_squared , c_squared ] ) idx , shapes = _getshapes_3d ( center , max_radius , space . shape ) subgrid = [ g [ idi ] for g , idi in zip ( grid , shapes ) ] squared_dist = [ ai * ( xi - x0i ) ** 2 for xi , ai , x0i in zip ( subgrid , scales , [ x0 , y0 , z0 ] ) ] radius = squared_dist [ 0 ] + ( squared_dist [ 1 ] + squared_dist [ 2 ] ) inside = radius <= 1 p [ idx ] [ inside ] += intensity return space . element ( p )
Create an ellipsoid phantom in 3d space.

    Parameters
    ----------
    space : `DiscreteLp`
        Space in which the phantom should be generated. If ``space.shape`` is
        1 in an axis, a corresponding slice of the phantom is created
        (instead of squashing the whole phantom into the slice).
    ellipsoids : list of lists
        Each row should contain the entries ::

            'value',
            'axis_1', 'axis_2', 'axis_3',
            'center_x', 'center_y', 'center_z',
            'rotation_phi', 'rotation_theta', 'rotation_psi'

        The provided ellipsoids need to be specified relative to the
        reference cube ``[-1, -1, -1] x [1, 1, 1]``. Angles are to be given
        in radians.

    Returns
    -------
    phantom : ``space`` element
        3D ellipsoid phantom in ``space``.

    See Also
    --------
    shepp_logan : The typical use-case for this function.