fix_direct.py

import numpy as np
import matplotlib.pyplot as plt
import sys
import abel
import scipy.integrate

transforms = [
  # ("BASEX"                ,  abel.basex.basex_transform),
  # ("linbasex"             , abel.linbasex.linbasex_transform),
  # must have higher symmetry for linbasex
  ("Direct-C"               ,  abel.direct.direct_transform),      
  ("Direct-P"               ,  abel.direct.direct_transform)]
  
  # ("Hansen-Law"           ,  abel.hansenlaw.hansenlaw_transform),
  # ("Onion-peeling (Bordas)", abel.onion_bordas.onion_bordas_transform),
  # ("Onion-peeling (Dasch)" , abel.dasch.onion_peeling_transform),
  # ("Three-point"           , abel.dasch.three_point_transform),
  # ("Two-point"             , abel.dasch.two_point_transform)]

ntrans = len(transforms)  # number of transforms

n     = 101

case = 'gaussian'
# case = 'circ'

if case == 'gaussian':
    r_max = n
    sigma = n*0.25
    
    ref  = abel.tools.analytical.GaussianAnalytical(n, r_max, sigma, symmetric=False)
    func = ref.func
    proj = ref.abel 
    

    r    = ref.r
    dr   = ref.dr
    
if case == 'circ':
    r_max = 1.0
    def a(n, x):
        return np.sqrt(n*n - x*x)
    
    r = np.linspace(0, r_max, n)
    func = np.ones_like(r)
    proj = 2*np.sqrt(1-r**2)
    dr   = r[1]-r[0]


## add noise:
#proj = proj + 0.1*np.abs(np.random.random(np.shape(ref.abel)))
    
fig, ax = plt.subplots(1, 1, figsize=(6,6), sharex=True, sharey=True)

ax.plot(r, func, label='Analytical', lw=1)

for row, (label, transFunc) in enumerate(transforms):

    alpha=1
    lw=1
    
    print dr
    if label == 'Direct-P':
        inverse = transFunc(np.copy(proj),dr=dr, direction='inverse', backend='python', correction=True)
    elif label == 'Direct-C':
        inverse = transFunc(np.copy(proj),dr=dr, direction='inverse', backend='c', correction=True)
        lw=4
        alpha=0.4
    else:
        inverse = transFunc(np.copy(proj),dr=dr, direction='inverse')
    
    ax.plot(r, inverse, label=label, ls='dashed', alpha=alpha, lw=lw)
    
    ax.plot(r, (inverse-func)*10, lw=lw, alpha=alpha, label=label+' error (x10)')
    ax.axhline(0, color='r', alpha=0.3, lw=1)
    
    ax.legend(loc='upper right', frameon=False, fontsize=9)
    
    ax.grid(ls='solid', alpha=0.05, color='k')
    ax.xaxis.set_tick_params(direction='in')
    ax.yaxis.set_tick_params(direction='in')
    
    
ax.set_xlabel("r (pixel)")
ax.set_ylabel('z')


ax.grid(ls='solid', alpha=0.05, color='k')
ax.xaxis.set_tick_params(direction='in')
ax.yaxis.set_tick_params(direction='in')
if case == 'gaussian':
    ax.set_ylim(-0.7,1.2)
    ax.set_xlim(0,n*0.75)
else:
    ax.set_xlim(0,1)


fig.subplots_adjust(left=0.08, bottom=0.07, right=0.98, top=0.99, hspace=0.03)
fig.savefig('gaussian.png', dpi=200)
plt.show()