Gaussian compare.py

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

transforms = [
  ("BASEX"                ,  abel.basex.basex_transform),   
  # ("linbasex"             , abel.linbasex.linbasex_transform),
  # must have higher symmetry for linbasex
  ("Direct"               ,  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     = 201

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 #+ 0.1*np.abs(np.random.random(np.shape(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]
    
    
fig, axs = plt.subplots(ntrans, 1, figsize=(4,7), sharex=True, sharey=True)

for row, (label, transFunc) in enumerate(transforms):
    axs[row].plot(r, func, label='Analytical', lw=1)

    inverse = transFunc(np.copy(proj),dr=dr, direction='inverse')
    axs[row].plot(r, inverse, label=label, ls='dashed', lw=1)
    
    # axs[row].plot(r, (inverse-func)*10, color='r', lw=1, label='Error (x10)')
    # axs[row].axhline(0, color='r', alpha=0.3, lw=1)
    
    axs[row].legend(loc='upper right', frameon=False, fontsize=9)
    
    
    axs[row].grid(ls='solid', alpha=0.05, color='k')
    axs[row].xaxis.set_tick_params(direction='in')
    axs[row].yaxis.set_tick_params(direction='in')
    
    
axs[-1].set_xlabel("x (pixel)")
axs[3].set_ylabel('y')

for ax in axs:
    ax.set_ylim(0.9,1.1)
    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_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()