# Example: rBasex beam block¶

```from __future__ import division, print_function

import numpy as np
import matplotlib.pyplot as plt
from copy import copy

import abel
from abel.tools.analytical import SampleImage
from abel.tools.vmi import rharmonics
from abel.rbasex import rbasex_transform

# This example demonstrates analysis of velocity-map images with "damaged"
# areas, in this case, with some parts obstructed by a beam block (see
# https://aip.scitation.org/doi/10.1063/1.4921990 for a practical example).
# First, a general Abel-transform method is used naively to demonstrate
# artifacts produced in the reconstructed image and the extracted radial
# distributions.
# Second, radial distributions are extracted from the same reconstructed image,
# but with its artifacts masked, showing good agreement with the actual
# distributions.
# Third, the rBasex method is used to transform the initial damaged image with
# the experimental artifacts masked, yielding a correct and cleaner
# reconstructed image and correct reconstructed distributions.

R = 150  # image radius
N = 2 * R + 1  # full image width and height
block_r = 20  # beam-block disk radius
block_w = 5  # beam-block holder width

vlim = 3.6  # intensity limits for images
ylim = (-1.3, 2.7)  # limits for plots

# create source distribution and its profiles for reference
source = SampleImage(N).func / 100
r_src, P0_src, P2_src = rharmonics(source)

# simulate experimental image:
# projection
im, _ = rbasex_transform(source, direction='forward')
# Poissonian noise
im[im < 0] = 0
im = np.random.RandomState(0).poisson(im)
# image coordinates
im_x = np.arange(float(N)) - R
im_y = R - np.arange(float(N))[:, None]
im_r = np.sqrt(im_x**2 + im_y**2)
im = im / (1 + np.exp(-(im_r - block_r)))
im[:R] *= 1 / (1 + np.exp(-(np.abs(im_x) - block_w)))

# reconstruct "as is" by a general Abel-transform method
rec_abel = abel.Transform(im, method='two_point').transform
# extract profiles "as is"
r_abel, P0_abel, P2_abel = rharmonics(rec_abel)
# extract profiles from masked reconstruction

# reconstruct masked image with rBasex
r_rbasex, P0_rbasex, P2_rbasex = distr_rbasex.rharmonics()

# plotting...
plt.figure(figsize=(7, 7))

cmap_hot = copy(plt.cm.hot)
cmap_seismic = copy(plt.cm.seismic)

def plots(row,
pr_title, r, P0, P2):
# input image
if im is not None:
plt.subplot(3, 4, 4 * row + 1)
plt.title(im_title)
plt.axis('off')

# transformed image
plt.subplot(3, 4, 4 * row + 2)
plt.title(tr_title)
plt.axis('off')

# profiles
plt.subplot(3, 2, 2 * row + 2)
plt.title(pr_title)
plt.plot(r_src, P0_src, 'C0--', lw=1)
plt.plot(r_src, P2_src, 'C3--', lw=1)
plt.plot(r, P0, 'C0', lw=1, label='\$P_0(r)\$')
plt.plot(r, P2, 'C3', lw=1, label='\$P_2(r)\$')
plt.xlim((0, R))
plt.ylim(ylim)
plt.legend()

plots(0,
'Raw image', im, None,
'Two-point', rec_abel, None,
'Profiles', r_abel, P0_abel, P2_abel)

plots(1,
None, None, None,