compare pbasex with PyAbel(basex)

test-pbasex.py

import numpy as np
import matplotlib.pyplot as plt
from pbasex import pbasex, loadG
from quadrant import foldQuadrant, resizeFolded
import abel

# Load images
im = np.loadtxt('O2-ANU1024.txt.bz2')
raw = abel.tools.center.center_image(im, square=True)
n = raw.shape[-1]
n2 = n//2

araw = abel.Transform(raw, method='basex',
                      transform_options=dict(reg=2000)).transform

# Fold images into quadrants
x0, y0 = n2, n2
quadrant_filter = [1,1,1,1]
folded = resizeFolded(foldQuadrant(raw, x0, y0, quadrant_filter), n2+1)

# Load inversion data
gData = loadG('G_r512_k128_l4.h5', 1)

# Apply the pBASEX algorithm
out = pbasex(folded, gData, make_images=True, alpha=4.1e-5)

r = np.arange(n2+1)
fig, ax = plt.subplots(2, 2)
ax[0, 0].set_title('pbasex')
ax[0, 0].imshow(out['inv'], extent=[-n2, n2, -n2, n2], vmin=0)
ax[0, 1].set_title('basex (PyAbel)')
ax[0, 1].imshow(araw, vmin=0, vmax=araw[100:].max()*0.6,
                extent=[-n2, n2, -n2, n2])
ax[1, 0].plot(out['inv'][n2-5:n2+5, n2:].sum(axis=0)*r/10000, lw=1)
ax[1, 0].set_xlabel('pixels')
ax[1, 1].plot(araw[n2-5:n2+5, n2:].sum(axis=0), lw=1)
ax[1, 1].set_xlabel('pixels')

plt.savefig('test-pbasex.png', dpi=100, bbox_inches='tight')
plt.show()