Source code for tomopy.misc.phantom

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# #########################################################################
# Copyright (c) 2015-2019, UChicago Argonne, LLC. All rights reserved.    #
#                                                                         #
# Copyright 2015-2019. UChicago Argonne, LLC. This software was produced  #
# under U.S. Government contract DE-AC02-06CH11357 for Argonne National   #
# Laboratory (ANL), which is operated by UChicago Argonne, LLC for the    #
# U.S. Department of Energy. The U.S. Government has rights to use,       #
# reproduce, and distribute this software.  NEITHER THE GOVERNMENT NOR    #
# UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR        #
# ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.  If software is     #
# modified to produce derivative works, such modified software should     #
# be clearly marked, so as not to confuse it with the version available   #
# from ANL.                                                               #
#                                                                         #
# Additionally, redistribution and use in source and binary forms, with   #
# or without modification, are permitted provided that the following      #
# conditions are met:                                                     #
#                                                                         #
#     * Redistributions of source code must retain the above copyright    #
#       notice, this list of conditions and the following disclaimer.     #
#                                                                         #
#     * Redistributions in binary form must reproduce the above copyright #
#       notice, this list of conditions and the following disclaimer in   #
#       the documentation and/or other materials provided with the        #
#       distribution.                                                     #
#                                                                         #
#     * Neither the name of UChicago Argonne, LLC, Argonne National       #
#       Laboratory, ANL, the U.S. Government, nor the names of its        #
#       contributors may be used to endorse or promote products derived   #
#       from this software without specific prior written permission.     #
#                                                                         #
# THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS     #
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT       #
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS       #
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago     #
# Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,        #
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,    #
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;        #
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER        #
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT      #
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN       #
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE         #
# POSSIBILITY OF SUCH DAMAGE.                                             #
# #########################################################################

"""
Module for generating synthetic phantoms.
"""

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

import numpy as np
import skimage
import skimage.transform
import tifffile
import os.path
import logging

logger = logging.getLogger(__name__)


__author__ = "Doga Gursoy"
__copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC."
__docformat__ = 'restructuredtext en'
__all__ = ['baboon',
           'barbara',
           'cameraman',
           'checkerboard',
           'lena',
           'peppers',
           'shepp2d',
           'shepp3d',
           'phantom']


DATA_PATH = os.path.abspath(
    os.path.join(os.path.dirname(__file__), '..', 'data'))

try:
    resize_kwargs = {'anti_aliasing': False}
    ignore = skimage.transform.resize(np.zeros(5), [2], mode='constant',
                                      **resize_kwargs)
except TypeError:
    logger.debug("Determined that the anti_aliasing keyword is not needed.")
    resize_kwargs = dict()


[docs]def baboon(size=512, dtype='float32'): """ Load test baboon image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'baboon.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) im = im.astype(dtype) return im
[docs]def barbara(size=512, dtype='float32'): """ Load test Barbara image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'barbara.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) return im.astype(dtype)
[docs]def cameraman(size=512, dtype='float32'): """ Load test cameraman image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'cameraman.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) return im.astype(dtype)
[docs]def checkerboard(size=512, dtype='float32'): """ Load test checkerboard image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'checkerboard.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) return im.astype(dtype)
[docs]def lena(size=512, dtype='float32'): """ Load test Lena image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'lena.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) return im.astype(dtype)
[docs]def peppers(size=512, dtype='float32'): """ Load test peppers image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'peppers.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) return im.astype(dtype)
[docs]def shepp2d(size=512, dtype='float32'): """ Load test Shepp-Logan image array. Parameters ---------- size : int or tuple of int, optional Size of the output image. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 2) fname = os.path.join(DATA_PATH, 'shepp2d.tif') im = tifffile.imread(fname) im = skimage.transform.resize(im, size, order=3, preserve_range=True, mode='constant', **resize_kwargs) im = np.expand_dims(im, 0) return im.astype(dtype)
def _totuple(size, dim): """ Converts size to tuple. """ if not isinstance(size, tuple): if dim == 2: size = (size, size) elif dim == 3: size = (size, size, size) return size
[docs]def shepp3d(size=128, dtype='float32'): """ Load 3D Shepp-Logan image array. Parameters ---------- size : int or tuple, optional Size of the 3D data. dtype : str, optional The desired data-type for the array. Returns ------- ndarray Output 3D test image. """ size = _totuple(size, 3) shepp_params = _array_to_params(_get_shepp_array()) return phantom(size, shepp_params, dtype).clip(0, np.inf)
[docs]def phantom(size, params, dtype='float32'): """ Generate a cube of given size using a list of ellipsoid parameters. Parameters ---------- size: tuple of int Size of the output cube. params: list of dict List of dictionaries with the parameters defining the ellipsoids to include in the cube. dtype: str, optional Data type of the output ndarray. Returns ------- ndarray 3D object filled with the specified ellipsoids. """ # instantiate ndarray cube obj = np.zeros(size, dtype=dtype) # define coords coords = _define_coords(size) # recursively add ellipsoids to cube for param in params: _ellipsoid(param, out=obj, coords=coords) return obj
def _ellipsoid(params, shape=None, out=None, coords=None): """ Generate a cube containing an ellipsoid defined by its parameters. If out is given, fills the given cube instead of creating a new one. """ # handle inputs if shape is None and out is None: raise ValueError("You need to set shape or out") if out is None: out = np.zeros(shape) if shape is None: shape = out.shape if len(shape) == 1: shape = shape, shape, shape elif len(shape) == 2: shape = shape[0], shape[1], 1 elif len(shape) > 3: raise ValueError("input shape must be lower or equal to 3") if coords is None: coords = _define_coords(shape) # rotate coords coords = _transform(coords, params) # recast as ndarray coords = np.asarray(coords) np.square(coords, out=coords) ellip_mask = coords.sum(axis=0) <= 1. ellip_mask.resize(shape) # fill ellipsoid with value out[ ellip_mask ] += params['A'] return out def _rotation_matrix(p): """ Defines an Euler rotation matrix from angles phi, theta and psi. """ cphi = np.cos(np.radians(p['phi'])) sphi = np.sin(np.radians(p['phi'])) ctheta = np.cos(np.radians(p['theta'])) stheta = np.sin(np.radians(p['theta'])) cpsi = np.cos(np.radians(p['psi'])) spsi = np.sin(np.radians(p['psi'])) alpha = [[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]] return np.asarray(alpha) def _define_coords(shape): """ Generate a tuple of coords in 3D with a given shape. """ mgrid = np.lib.index_tricks.nd_grid() cshape = np.asarray(1j) * shape x, y, z = mgrid[-1:1:cshape[0], -1:1:cshape[1], -1:1:cshape[2]] return x, y, z def _transform(coords, p): """ Apply rotation, translation and rescaling to a 3-tuple of coords. """ alpha = _rotation_matrix(p) out_coords = np.tensordot(alpha, coords, axes=1) _shape = (3,) + (1,) * ( out_coords.ndim - 1 ) _dt = out_coords.dtype M0 = np.array([p['x0'], p['y0'], p['z0']], dtype=_dt).reshape(_shape) sc = np.array([p['a'], p['b'], p['c']], dtype=_dt).reshape(_shape) out_coords -= M0 out_coords /= sc return out_coords def _get_shepp_array(): """ Returns the parameters for generating modified Shepp-Logan phantom. """ shepp_array = [ [1., .6900, .920, .810, 0., 0., 0., 90., 90., 90.], [-.8, .6624, .874, .780, 0., -.0184, 0., 90., 90., 90.], [-.2, .1100, .310, .220, .22, 0., 0., -108., 90., 100.], [-.2, .1600, .410, .280, -.22, 0., 0., 108., 90., 100.], [.1, .2100, .250, .410, 0., .35, -.15, 90., 90., 90.], [.1, .0460, .046, .050, 0., .1, .25, 90., 90., 90.], [.1, .0460, .046, .050, 0., -.1, .25, 90., 90., 90.], [.1, .0460, .023, .050, -.08, -.605, 0., 90., 90., 90.], [.1, .0230, .023, .020, 0., -.606, 0., 90., 90., 90.], [.1, .0230, .046, .020, .06, -.605, 0., 90., 90., 90.]] return shepp_array def _array_to_params(array): """ Converts list to a dictionary. """ # mandatory parameters to define an ellipsoid params_tuple = [ 'A', 'a', 'b', 'c', 'x0', 'y0', 'z0', 'phi', 'theta', 'psi'] array = np.asarray(array) out = [] for i in range(array.shape[0]): tmp = dict() for k, j in zip(params_tuple, list(range(array.shape[1]))): tmp[k] = array[i, j] out.append(tmp) return out