Three Point

Introduction

The “Three Point” Abel transform method exploits the observation that the value of the Abel inverted data at any radial position r is primarily determined from changes in the projection data in the neighborhood of r. This technique was developed by Dasch [1].

How it works

The projection data (raw data \(\mathbf{P}\)) is expanded as a quadratic function of \(r - r_{j*}\) in the neighborhood of each data point in \(\mathbf{P}\). In other words, \(\mathbf{P}'(r) = dP/dr\) is estimated using a 3-point approximation (to the derivative), similar to central differencing. Doing so enables an analytical integration of the inverse Abel integral around each point \(r_j\). The result of this integration is expressed as a linear operator \(\mathbf{D}\), operating on the projection data \(\mathbf{P}\) to give the underlying radial distribution \(\mathbf{F}\).

When to use it

Dasch recommends this method based on its speed of implementation, robustness in the presence of sharp edges, and low noise. He also notes that this technique works best for cases where the real difference between adjacent projections is much greater than the noise in the projections. This is important, because if the projections are oversampled (raw data \(\mathbf{P}\) taken with data points very close to each other), the spacing between adjacent projections is decreased, and the real difference between them becomes comparable with the noise in the projections. In such situations, the deconvolution is highly inaccurate, and the projection data \(\mathbf{P}\) must be smoothed before this technique is used. (Consider smoothing with scipy.ndimage.gaussian_filter.)

How to use it

To complete the inverse transform of a full image with the three_point method, simply use the abel.Transform class:

abel.Transform(myImage, method='three_point', direction='inverse').transform

Note that the forward Three point transform is not yet implemented in PyAbel.

If you would like to access the Three Point algorithm directly (to transform a right-side half-image), you can use abel.dasch.three_point_transform().

Example

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

"""example_dasch_methods.py.
"""

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

# Dribinski sample image size 501x501
n = 501
IM = abel.tools.analytical.SampleImage(n).func

# split into quadrants
origQ = abel.tools.symmetry.get_image_quadrants(IM)

# speed distribution of original image
orig_speed = abel.tools.vmi.angular_integration_3D(origQ[0], origin=(-1, 0))
scale_factor = orig_speed[1].max()

plt.plot(orig_speed[0], orig_speed[1]/scale_factor, linestyle='dashed',
         label="Dribinski sample")


# forward Abel projection
fIM = abel.Transform(IM, direction="forward", method="hansenlaw").transform

# split projected image into quadrants
Q = abel.tools.symmetry.get_image_quadrants(fIM)

dasch_transform = {
    "two_point": abel.dasch.two_point_transform,
    "three_point": abel.dasch.three_point_transform,
    "onion_peeling": abel.dasch.onion_peeling_transform
}

for method in dasch_transform.keys():
    Q0 = Q[0].copy()
# method inverse Abel transform
    AQ0 = dasch_transform[method](Q0)
# speed distribution
    speed = abel.tools.vmi.angular_integration_3D(AQ0, origin=(-1, 0))

    plt.plot(speed[0], speed[1]*orig_speed[1][14]/speed[1][14]/scale_factor,
             label=method)

plt.title("Dasch methods for Dribinski sample image $n={:d}$".format(n))
plt.xlim((0, 250))
plt.legend(loc='upper center', bbox_to_anchor=(0.35, 1), frameon=False)
plt.tight_layout()
# plt.savefig("plot_example_dasch_methods.png",dpi=100)
plt.show()
../_images/example_dasch_methods.svg

Notes

The algorithm contained two typos in Eq. (7) in the original citation [1]. A corrected form of these equations is presented in Karl Martin’s 2002 PhD thesis [2]. PyAbel uses the corrected version of the algorithm.

For more information on the PyAbel implementation of the three-point algorithm, please see issue #61 and Pull Request #64.

Citation