dribinski_compare.py

import numpy as np
import matplotlib.pyplot as plt
import sys, os
import abel
import scipy.ndimage
import bz2

transforms = [
  ("basex"         ,abel.basex.basex_transform),
  ("direct"        ,abel.direct.direct_transform),
  ("hansenlaw"     ,abel.hansenlaw.hansenlaw_transform),
  ("onion_bordas"  ,abel.onion_bordas.onion_bordas_transform),
  ("onion_peeling" ,abel.dasch.onion_peeling_transform),
  ("three_point"   ,abel.dasch.three_point_transform),
  ("two_point"     ,abel.dasch.two_point_transform),
  ("linbasex"      ,abel.linbasex.linbasex_transform)]


ntrans = len(transforms)  # number of transforms

infile = bz2.BZ2File('../data/O2-ANU1024.txt.bz2')
IM = np.loadtxt(infile)

IModd = abel.tools.center.center_image(IM, center="convolution", odd_size=True, square=True)
IModd = abel.tools.analytical.SampleImage(1024).image + 100*np.random.random((1024,1024))

Q = abel.tools.symmetry.get_image_quadrants(IModd, reorient=True)
Q0 = Q[1]

h, w = np.shape(IM)

fig, axs  = plt.subplots(4, 2, figsize=(3.487,7.3), sharex=True, sharey=True)
fig1,axs1 = plt.subplots(3, 1, figsize=(3.487,5))

for num, (ax, (label, transFunc)) in enumerate(zip(axs.ravel(), transforms)):
    print(label)
    
    # if label == 'Direct-P':
    #     targs = dict(backend='python', int_func=np.trapz, dr=0.1)
    #
    # elif label == 'direct_C':
    #     targs = dict(backend='c', dr=0.1)
        
    if label == 'linbasex':
        targs = dict(rcond=0.0001, smoothing=0.5)
    
    else:
        targs = dict()
    
    if label == 'basex':
        targs = dict(reg=200)
    
    # T = abel.Transform(IM, method=transFunc, center='convolution',
    #                    transform_options=targs,
    #                    center_options=dict(square=True),
    #                    angular_integration=True)

        
    
    trans = transFunc(Q0, direction="inverse", **targs)
    
    if label == 'linbasex':
        # trans *= np.average(prev/trans) * 1.
        pass
    r, inten = abel.tools.vmi.angular_integration(trans, origin=(0, 0), dr=0.1)

    inten /= 1e6
    
    im = ax.imshow(trans[::-1], cmap='gist_heat_r',
              origin='lower', aspect='auto', vmin=0,
              vmax=20)
    
    ax.set_title(label, fontsize=10, x=0.96, y=0.90, ha='right', va='top',
                 weight='bold', color='k')

    if label == 'linbasex':
        pargs = dict(alpha=0.3, zorder=0)
        
    else:
        pargs = dict()
        
        
    axs1[0].plot(r, inten, lw=1,             **pargs)
    axs1[1].plot(r, inten, lw=1,label=label, **pargs)
    axs1[2].plot(r, inten, lw=1,             **pargs)
    
    prev = trans
    

axc = fig.add_axes([0.12, 0.92, 0.85, 0.02])
cbar = plt.colorbar(im, orientation="horizontal", cax=axc, label='Intensity')
cbar.ax.xaxis.set_ticks_position('top')
cbar.ax.xaxis.set_label_position('top')

for label in cbar.ax.xaxis.get_ticklabels():
    label.set_fontsize(6)


for ax in axs.ravel():
    ax.set_xlim(0, 512)
    ax.set_ylim(0, 512)
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontsize(6)
    
fig.subplots_adjust(left=0.13, bottom=0.06, right=0.99, top=0.91, wspace=0.08,
                    hspace=0.08)

axs[-1,0].set_xlabel('$r$ (pixels)')
axs[-1,1].set_xlabel('$r$ (pixels)')

for ax in axs[:,0]:
    ax.set_ylabel('$z$ (pixels)')


for ax in axs1:
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontsize(8)
    
    ax.grid(color='k', alpha=0.1)

axs1[0].set_xlim(0, 512)
axs1[1].set_xlim(355, 385)
axs1[1].legend(fontsize=6)

axs1[2].set_xlim(100, 300)
axs1[2].set_ylim(-0.07, 0.20)

def place_letter(letter, ax, color='k', offset=(0,0)):
    ax.annotate(letter, xy=(0.02, 0.97), xytext=offset,
                xycoords='axes fraction', textcoords='offset points',
                color=color, ha='left', va='top', weight='bold')
    
for ax, letter in zip(axs.ravel(), 'abcdefgh'):
    place_letter(letter+')', ax)

for ax, letter in zip(axs1.ravel(), 'abcdefgh'):
    place_letter(letter+')', ax, color='k')

fig1.subplots_adjust(left=0.13, bottom=0.09, right=0.96, top=0.98, hspace=0.19)
axs1[-1].set_xlabel('$r$ (pixels)')

fig.savefig('dribinski.png', dpi=200)
fig1.savefig('dribinski.png', dpi=200)
plt.show()

gistfile1.txt

import numpy as np
import matplotlib.pyplot as plt
import sys, os
import abel
import scipy.ndimage
import bz2

transforms = [
  ("basex"         ,abel.basex.basex_transform),
  ("direct"        ,abel.direct.direct_transform),
  ("hansenlaw"     ,abel.hansenlaw.hansenlaw_transform),
  ("onion_bordas"  ,abel.onion_bordas.onion_bordas_transform),
  ("onion_peeling" ,abel.dasch.onion_peeling_transform),
  ("three_point"   ,abel.dasch.three_point_transform),
  ("two_point"     ,abel.dasch.two_point_transform),
  ("linbasex"      ,abel.linbasex.linbasex_transform)]


ntrans = len(transforms)  # number of transforms

infile = bz2.BZ2File('../data/O2-ANU1024.txt.bz2')
IM = np.loadtxt(infile)

IModd = abel.tools.center.center_image(IM, center="convolution", odd_size=True, square=True)
IModd = abel.tools.analytical.SampleImage(1024).image + 100*np.random.random((1024,1024))

Q = abel.tools.symmetry.get_image_quadrants(IModd, reorient=True)
Q0 = Q[1]

h, w = np.shape(IM)

fig, axs  = plt.subplots(4, 2, figsize=(3.487,7.3), sharex=True, sharey=True)
fig1,axs1 = plt.subplots(3, 1, figsize=(3.487,5))

for num, (ax, (label, transFunc)) in enumerate(zip(axs.ravel(), transforms)):
    print(label)
    
    # if label == 'Direct-P':
    #     targs = dict(backend='python', int_func=np.trapz, dr=0.1)
    #
    # elif label == 'direct_C':
    #     targs = dict(backend='c', dr=0.1)
        
    if label == 'linbasex':
        targs = dict(rcond=0.0001, smoothing=0.5)
    
    else:
        targs = dict()
    
    if label == 'basex':
        targs = dict(reg=200)
    
    # T = abel.Transform(IM, method=transFunc, center='convolution',
    #                    transform_options=targs,
    #                    center_options=dict(square=True),
    #                    angular_integration=True)

        
    
    trans = transFunc(Q0, direction="inverse", **targs)
    
    if label == 'linbasex':
        # trans *= np.average(prev/trans) * 1.
        pass
    r, inten = abel.tools.vmi.angular_integration(trans, origin=(0, 0), dr=0.1)

    inten /= 1e6
    
    im = ax.imshow(trans[::-1], cmap='gist_heat_r',
              origin='lower', aspect='auto', vmin=0,
              vmax=20)
    
    ax.set_title(label, fontsize=10, x=0.96, y=0.90, ha='right', va='top',
                 weight='bold', color='k')

    if label == 'linbasex':
        pargs = dict(alpha=0.3, zorder=0)
        
    else:
        pargs = dict()
        
        
    axs1[0].plot(r, inten, lw=1,             **pargs)
    axs1[1].plot(r, inten, lw=1,label=label, **pargs)
    axs1[2].plot(r, inten, lw=1,             **pargs)
    
    prev = trans
    

axc = fig.add_axes([0.12, 0.92, 0.85, 0.02])
cbar = plt.colorbar(im, orientation="horizontal", cax=axc, label='Intensity')
cbar.ax.xaxis.set_ticks_position('top')
cbar.ax.xaxis.set_label_position('top')

for label in cbar.ax.xaxis.get_ticklabels():
    label.set_fontsize(6)


for ax in axs.ravel():
    ax.set_xlim(0, 512)
    ax.set_ylim(0, 512)
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontsize(6)
    
fig.subplots_adjust(left=0.13, bottom=0.06, right=0.99, top=0.91, wspace=0.08,
                    hspace=0.08)

axs[-1,0].set_xlabel('$r$ (pixels)')
axs[-1,1].set_xlabel('$r$ (pixels)')

for ax in axs[:,0]:
    ax.set_ylabel('$z$ (pixels)')


for ax in axs1:
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontsize(8)
    
    ax.grid(color='k', alpha=0.1)

axs1[0].set_xlim(0, 512)
axs1[1].set_xlim(355, 385)
axs1[1].legend(fontsize=6)

axs1[2].set_xlim(100, 300)
axs1[2].set_ylim(-0.07, 0.20)

def place_letter(letter, ax, color='k', offset=(0,0)):
    ax.annotate(letter, xy=(0.02, 0.97), xytext=offset,
                xycoords='axes fraction', textcoords='offset points',
                color=color, ha='left', va='top', weight='bold')
    
for ax, letter in zip(axs.ravel(), 'abcdefgh'):
    place_letter(letter+')', ax)

for ax, letter in zip(axs1.ravel(), 'abcdefgh'):
    place_letter(letter+')', ax, color='k')

fig1.subplots_adjust(left=0.13, bottom=0.09, right=0.96, top=0.98, hspace=0.19)
axs1[-1].set_xlabel('$r$ (pixels)')

fig.savefig('dribinski.png', dpi=200)
fig1.savefig('dribinski.png', dpi=200)
plt.show()