check_all.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)

sigma = 10

# 2D Gaussian -----------
# pixel coordinates
n = 200 
rows, cols = n, n
r2 = rows//2
c2 = cols//2

x = np.linspace(-c2+0.5, c2-0.5, cols)
y = np.linspace(-r2+0.5, r2-0.5, rows)


X, Y = np.meshgrid(x, y)

IM = gauss(X, 0, sigma) # cylinderical Gaussian max located at pixel=0
                        # @DanHickstein see #155
Q0, Q1, Q2, Q3 = abel.tools.symmetry.get_image_quadrants(IM)
qspeed = abel.tools.vmi.angular_integration(Q0, origin=(0, 0))

# Abel transform
AQ0 = abel.onion_peeling.onion_peeling_transform(Q0, shift_grid=False)
ospeed = abel.tools.vmi.angular_integration(AQ0, origin=(0, 0))

plt.subplot(221)
plt.plot(*qspeed, label="orig.")
plt.plot(ospeed[0], ospeed[1]*qspeed[1][2]/ospeed[1][2], label="onion")
plt.legend()
plt.title("n=100", fontsize=10)
plt.axis(xmax=49)

n = 201 
rows, cols = n, n
r2 = rows//2
c2 = cols//2

x = np.linspace(-c2, c2, cols)
y = np.linspace(-r2, r2, rows)

X, Y = np.meshgrid(x, y)

IM = gauss(X, 0, sigma) # cylinderical Gaussian max located at pixel=0
                        # @DanHickstein see #155
#Q0 = IM[:r2, c2:] # quadrant, top-right
Q0, Q1, Q2, Q3 = abel.tools.symmetry.get_image_quadrants(IM)
Q0_copy = Q0.copy()
qspeed = abel.tools.vmi.angular_integration(Q0, origin=(0, 0))


# Plots
plt.suptitle("inverse Abel transforms 2D-Gaussian cylinder $\sigma={:g}$".format(sigma))

Q0 = Q0_copy
AQ0 = abel.hansenlaw.hansenlaw_transform(Q0)
hspeed = abel.tools.vmi.angular_integration(AQ0, origin=(0, 0))
plt.subplot(222)
plt.plot(*qspeed, label="orig.")
plt.plot(hspeed[0], hspeed[1]*qspeed[1][2]/hspeed[1][2], label="hansenlaw")
plt.legend()
plt.title("n=101", fontsize=10)
plt.axis(xmax=49)

Q0 = Q0_copy
AQ0 = abel.basex.basex_transform(Q0)
bspeed = abel.tools.vmi.angular_integration(AQ0, origin=(0, 0))
plt.subplot(223)
plt.plot(*qspeed, label="orig.")
plt.plot(bspeed[0], bspeed[1]*qspeed[1][2]/bspeed[1][2], label="basex")
plt.legend()
plt.legend()
plt.title("n=101", fontsize=10)
plt.axis(xmax=49)

Q0 = Q0_copy
AQ0 = abel.three_point.three_point_transform(Q0)
tspeed = abel.tools.vmi.angular_integration(AQ0, origin=(0, 0))
plt.subplot(224)
plt.plot(*qspeed, label="orig.")
plt.plot(tspeed[0], tspeed[1]*qspeed[1][2]/tspeed[1][2], label="three_point")
plt.legend()
plt.title("n=101", fontsize=10)
plt.axis(xmax=49)


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