#!/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 phase retrieval.
This code is part of tomopy: https://github.com/tomopy/tomopy
It was adapted for being stand-alone, without job distribution.
See the license above for more information.
"""
import numpy as np
fft2 = np.fft.fft2
ifft2 = np.fft.ifft2
__author__ = "Doga Gursoy"
__credits__ = "Mark Rivers, Xianghui Xiao"
__copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC."
__docformat__ = 'restructuredtext en'
__all__ = ['retrieve_phase']
BOLTZMANN_CONSTANT = 1.3806488e-16 # [erg/k]
SPEED_OF_LIGHT = 299792458e+2 # [cm/s]
PI = 3.14159265359
PLANCK_CONSTANT = 6.58211928e-19 # [keV*s]
def _wavelength(energy):
return 2 * PI * PLANCK_CONSTANT * SPEED_OF_LIGHT / energy
[docs]
def retrieve_phase(
tomo, pixel_size=1e-4, dist=50, energy=20,
alpha=1e-3, pad=True, ncore=None, nchunk=None):
"""
Perform single-step phase retrieval from phase-contrast measurements.
Parameters
----------
tomo : ndarray
3D tomographic data.
pixel_size : float, optional
Detector pixel size in cm.
dist : float, optional
Propagation distance of the wavefront in cm.
energy : float, optional
Energy of incident wave in keV.
alpha : float, optional
Regularization parameter.
pad : bool, optional
If True, extend the size of the projections by padding with zeros.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Approximated 3D tomographic phase data.
"""
# New dimensions and pad value after padding.
py, pz, val = _calc_pad(tomo, pixel_size, dist, energy, pad)
# Compute the reciprocal grid.
dx, dy, dz = tomo.shape
w2 = _reciprocal_grid(pixel_size, dy + 2 * py, dz + 2 * pz)
# Filter in Fourier space.
phase_filter = np.fft.fftshift(
_paganin_filter_factor(energy, dist, alpha, w2))
prj = np.full((dy + 2 * py, dz + 2 * pz), val, dtype='float32')
arr = _retrieve_phase(tomo, phase_filter, py, pz, prj, pad)
return arr
def _retrieve_phase(tomo, phase_filter, px, py, prj, pad):
dx, dy, dz = tomo.shape
num_jobs = tomo.shape[0]
normalized_phase_filter = phase_filter / phase_filter.max()
for m in range(num_jobs):
# Padding "constant" with border value
# prj is initially filled with "val"
prj[px:dy + px, py:dz + py] = tomo[m]
prj[:px] = prj[px]
prj[-px:] = prj[-px-1]
prj[:, :py] = prj[:, py][:, np.newaxis]
prj[:, -py:] = prj[:, -py-1][:, np.newaxis]
fproj = fft2(prj)
fproj *= normalized_phase_filter
proj = np.real(ifft2(fproj))
if pad:
proj = proj[px:dy + px, py:dz + py]
tomo[m] = proj
return tomo
def _calc_pad(tomo, pixel_size, dist, energy, pad):
"""
Calculate new dimensions and pad value after padding.
Parameters
----------
tomo : ndarray
3D tomographic data.
pixel_size : float
Detector pixel size in cm.
dist : float
Propagation distance of the wavefront in cm.
energy : float
Energy of incident wave in keV.
pad : bool
If True, extend the size of the projections by padding with zeros.
Returns
-------
int
Pad amount in projection axis.
int
Pad amount in sinogram axis.
float
Pad value.
"""
dx, dy, dz = tomo.shape
wavelength = _wavelength(energy)
py, pz, val = 0, 0, 0
if pad:
val = _calc_pad_val(tomo)
py = _calc_pad_width(dy, pixel_size, wavelength, dist)
pz = _calc_pad_width(dz, pixel_size, wavelength, dist)
return py, pz, val
def _paganin_filter_factor(energy, dist, alpha, w2):
return 1 / (_wavelength(energy) * dist * w2 / (4 * PI) + alpha)
def _calc_pad_width(dim, pixel_size, wavelength, dist):
pad_pix = np.ceil(PI * wavelength * dist / pixel_size ** 2)
return int((pow(2, np.ceil(np.log2(dim + pad_pix))) - dim) * 0.5)
def _calc_pad_val(tomo):
# mean of [(column 0 of radio) + (column -1 of radio)]/2.
return np.mean((tomo[..., 0] + tomo[..., -1]) * 0.5)
def _reciprocal_grid(pixel_size, nx, ny):
"""
Calculate reciprocal grid.
Parameters
----------
pixel_size : float
Detector pixel size in cm.
nx, ny : int
Size of the reciprocal grid along x and y axes.
Returns
-------
ndarray
Grid coordinates.
"""
# Sampling in reciprocal space.
indx = _reciprocal_coord(pixel_size, nx)
indy = _reciprocal_coord(pixel_size, ny)
np.square(indx, out=indx)
np.square(indy, out=indy)
return np.add.outer(indx, indy)
def _reciprocal_coord(pixel_size, num_grid):
"""
Calculate reciprocal grid coordinates for a given pixel size
and discretization.
Parameters
----------
pixel_size : float
Detector pixel size in cm.
num_grid : int
Size of the reciprocal grid.
Returns
-------
ndarray
Grid coordinates.
"""
n = num_grid - 1
rc = np.arange(-n, num_grid, 2, dtype = np.float32)
rc *= 0.5 / (n * pixel_size)
return rc