onion peeling 1d- and 2d-gaussian testing

check_onion.py

# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

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

gauss = lambda r, r0, sigma: np.exp(-(r-r0)**2/sigma**2)

n = 100
rows, cols = n, n
r2 = rows//2 + rows % 2
c2 = cols//2 + cols % 2

sigma = 20*n/100

# 1D Gaussian -----------
r = np.linspace(0, c2-1, c2)

orig = gauss(r, 0, sigma)
orig_copy = orig.copy()

recon = abel.onion_peeling.onion_peeling_transform(orig, shift_grid=False)

ratio_1d = orig_copy[10]/recon[10]
print("1d ratio = {:g}".format(ratio_1d))

# 2D Gaussian -----------
# pixel coordinates
x = np.linspace(-c2, c2, cols+1)
y = np.linspace(-r2, r2, rows+1)

X, Y = np.meshgrid(x, y)
R, THETA = abel.tools.polar.cart2polar(X, Y)

IM = gauss(R, 0, sigma) # Gaussian max located at pixel=0
Q0 = IM[:r2, c2:] # quadrant, top-right
Q0_copy = Q0.copy()

# onion_peeling inverse Abel transform
AQ0 = abel.onion_peeling.onion_peeling_transform(Q0, shift_grid=False)

profQ0 = Q0_copy[-10:,:].sum(axis=0)
profAQ0 = AQ0[-10:,:].sum(axis=0)
ratio_2d = profQ0[10]/profAQ0[10]
print("2d ratio = {:g}".format(ratio_2d))

# Plots
plt.suptitle("onion peeling inverse Abel transform, $n={:d},"
             " \sigma={:g}$".format(n, sigma))
plt.subplot(121)
plt.plot(orig_copy, label="orig.")

plt.plot(recon*ratio_1d, label="onion")
plt.title("1D Gaussian", fontsize="small")
plt.xlabel("pixel radius")
plt.legend()

plt.subplot(122)
plt.plot(profQ0, label="orig.")


plt.plot(profAQ0*ratio_2d, label="onion")
plt.title("2D Gaussian", fontsize="small")
plt.xlabel("pixel radius")
plt.legend()

plt.savefig("onion_1d_2d.png", dpi=100)
plt.show()