pyabel_fh.py

# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
import numpy as np
import abel
from scipy.special import jn
from scipy.ndimage import zoom
import matplotlib.pyplot as plt

#############################################################################
#
# Fourier-Hankel method:    A = H F
#
# 2017-08-29 
#  see https://github.com/PyAbel/PyAbel/issues/24#issuecomment-325547436
#  Steven Murray - implementation
#  Stephen Gibson, Dan Hickstein - adapted for PyAbel
#
#############################################################################

def construct_dht_matrix(N, nu  = 0, b = 1):
    m = np.arange(0, N).astype('float')
    n = np.arange(0, N).astype('float')
    return b * jn(nu, np.outer(b*m, n/N))


def dht(X, d = 1, nu = 0, axis=-1, b = 1):
    N = X.shape[axis]
    F = construct_dht_matrix(N,nu,b)*np.arange(0,N)
    return d**2 * np.tensordot(F, X, axes=([1],[axis])) 


def dft(X, axis=-1):
    # discrete Fourier transform
    n = X.shape[axis]
        
    # Build a slicer to remove last element from flipped array
    slc = [slice(None)] * len(X.shape)
    slc[axis] = slice(None,-1)
    
    X = np.append(np.flip(X,axis)[slc], X, axis=axis)  # make symmetric

    return np.abs(np.fft.rfft(X, axis=axis)[:n])/n


def hankel_fourier_transform(X, d=1, nu = 0, inverse=True, axis=-1):
    n = X.shape[axis]
    
    if inverse:
        fx = dft(X, axis=axis)  # Fourier transform
        hf = dht(fx, d=1/(n*d), nu=nu, b=np.pi, axis=axis) * n / 2  # Hankel transform
    else:
        hx = dht(X, d=1/(n*d), nu=nu, b=np.pi, axis=axis) * 2
        hf = dft(hx)
    
    return hf


def fourier_hankel_transform(IM, dr=1, inverse=True, basis_dir=None):
    """
    Parameters
    ----------
    IM : 1D or 2D numpy array
        Right-side half-image (or quadrant).
    dr : float
        Sampling size (=1 for pixel images), used for Jacobian scaling.
    direction : str
        'inverse' or 'forward' Abel transform
    Returns
    -------
    trans_IM : 1D or 2D numpy array
        Inverse or forward Abel transform half-image, the same shape as IM.
    """
    IM = np.atleast_2d(IM)
    transform_IM = hankel_fourier_transform(IM, inverse=inverse)
    
    if transform_IM.shape[0] == 1:
        transform_IM = transform_IM[0]   # flatten to a vector

    return transform_IM
 
n = 257
sigma = 20*n/100.
IM = abel.tools.analytical.GaussianAnalytical(n, n, sigma=sigma,
                                              symmetric=False)

fhIM = fourier_hankel_transform(IM.func)
tpIM = abel.dasch.three_point_transform(IM.func)

fig, ax = plt.subplots()
ax.plot(IM.r, IM.func/np.sqrt(np.pi)/sigma,
        label=r'Gaussian$/\sigma\sqrt{\pi}$')
ax.plot(IM.r, fhIM, label=r'Fourier-Hankel')
ax.plot(IM.r, tpIM , label=r'3pt')
ax.set_title('Fourier-Hankel 1-d test: $n={}$'.format(n))
ax.axis(xmax=n*0.6)
ax.set_xlabel('pixels')
ax.set_ylabel('intensity')
ax.legend()

plt.savefig("Fourier-Hankel-1D-Gauss.png", dpi=75)
plt.show()