abelloop.py

# from abel.analytical import GaussianAnalytical
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate
import time
import warnings
# warnings.filterwarnings("ignore")

r = np.linspace(0,10,51)
height = 5
sigma  = 3
A0     = 1

orig = A0 * np.exp(-r**2/sigma**2)  # analytical gaussian
proj = sigma * np.sqrt(np.pi) * A0*np.exp(-r**2/sigma**2) # analytical projection

orig2d = np.tile(orig, (height, 1))

def zero_loop_abel(IM,r):
    # in the following,
    # first dimension is z
    # second dimention is r
    # third dimension is y
    height,size = np.shape(IM)

    ORIG = np.zeros((height,size,size))
    ORIG = IM[:,:,None]

    plt.show()
    Z,R,Y = np.meshgrid(np.ones(height),r,r,indexing='ij') # set up grids

    dr = r[1]-r[0]

    # set up the integral: https://en.wikipedia.org/wiki/Abel_transform
    P = 2 * ORIG * R * dr/ np.sqrt(R**2-Y**2)
    P[R<=Y] = 0   # Set the limits of r from y -> inf

    F = scipy.integrate.simps(P,axis=1) # take the integral
    return F


def one_loop_abel(IM,r):
    size = np.shape(r)[0]
    dr = r[1]-r[0] # this assumes even spacing of r, but this is not strictly necessary
    res = np.zeros((height,size))

    for index,row in enumerate(IM): # loop over rows
        ORIG = row[:,None]
        R,Y = np.meshgrid(r,r,indexing='ij') # set up grids

        P = 2. * ORIG * R * dr / np.sqrt(R**2-Y**2) # set up the integral

        P[R<=Y] = 0   # set limits of integral
        F = scipy.integrate.simps(P,axis=0) # take the integral    
        res[index] = F

    return res


def two_loop_abel(IM,r):
    # start of Abel transform: 
    y_list = r # Here we use the same r-spacing for y, but this does not need to be the case
    dr = r[1]-r[0] # this assumes even spacing of r, but this is not strictly necessary
    res = np.zeros(np.shape(IM))

    for i,row in enumerate(IM): # loops over rows
        for j,y in enumerate(y_list):  # loop over y-values of final image
            if y != np.max(r):
                ind = (r>y)
                P = 2 * row[ind] * r[ind] * dr / np.sqrt(r[ind]**2-y**2)
                F = scipy.integrate.simps(P) # take the integral
                res[i,j] = F
            else:
                res[i,j] = 0

    return res


t = time.time()
zero = zero_loop_abel(orig2d,r)
zero = np.mean(zero,axis=0)
print 'zero loops: %.2f sec'%(time.time()-t)

t = time.time()
one = one_loop_abel(orig2d,r)
one = np.mean(one,axis=0)
print 'one loop: %.2f sec'%(time.time()-t)

t = time.time()
two = two_loop_abel(orig2d,r)
two = np.mean(two,axis=0)
print 'two loop: %.2f sec'%(time.time()-t)

plt.plot(r,orig,label='Original',color='k',alpha=0.5)
plt.plot(r,proj,label='Projection (analytical)',color='b',lw=2)
plt.plot(r,zero,label='Projection (zero-loop)',color='r',ls='dashed',lw=2)
plt.plot(r,one, label='Projection (one-loop)', color='m',ls='solid', lw=1)
plt.plot(r,two, label='Projection (two-loop)', color='r',ls='solid', lw=1)

plt.legend(loc=0)
plt.show()