Equation for calculating pixel position from two theta and azimuth angles

[CPrescher]

Hi,

i am currently working on a distortion correction for large area detectors with holeplates for relatively small distortions (below 1 px). However, this requires a lot of different optimizations to figure out, the exact rotation or tilt of holeplate in respect to detector etc.

For this to be fast, it would be really nice to have an equation or function which can calculate a pixel position (with subpixel resolution) based on two theta and azimuth angles. I was searching through the pyFAI code but could not find it. Unfortunately, my math skills seem to be not good enough to reverse the equation.

Is there a solution for this?

Best wishes, Clemens

[vallsv]

Hi Clemens,

Jerome have designed an InvertGeometry.

It is used by our GUI. I can't let you it is really accurate, cause it was not designed for.

I am not sure, but here is few lines of code which maybe fit what you was searching for.

# tth/chi -> x/y
            invertGeometry = InvertGeometry(
                geometry.array_from_unit(typ="center", unit=pyFAI.units.TTH_RAD, scale=False),
                geometry.chiArray())
            pixel = invertGeometry(tthRad, chiRad, True)

# x/y -> tth/chi
            ax, ay = numpy.array([pixel[1]]), numpy.array([pixel[0]])
            tth = geometry.tth(ay, ax)[0]
            chi = geometry.chi(ay, ax)[0]

[kif]

Hi Clemens,

As Valentin commented, the transformation is not completely bijective and multiple pixels can have the same 2th/chi ... one example is tiled detectors like the "initial" Pilatus detector (different from the one sold by Dectris today) described in https://doi.org/10.1016/j.nima.2005.05.032

The InverGeometry class (you can find it in the pyFAI.ext.invert_geometry, it is written in Cython for performances) is doing the inversion in a rather coarse way. It initially searches the best pixel then performs an interpolation based on 2 neighboring pixels in each direction (hoping no module boundaries here). For GUI application is it good enough. For geometry correction, with sub-pixel precision, I would rather perform this refinement based on the corner-array which provides the corner coordinates of pixels. This would removes the burden of module boundaries.

[CPrescher]

Hi,

thank you for your answer. I just tried it out and found it is much slower than my current implementation for 5000 points (1h43min to 6.81s). I think it tries to create an interpolate for each point right?

The way I currently get the xy position is the following:

tth_interp2d= interp2d(grid_x, grid_y, geom.twoThetaArray(), kind='cubic')

chi_interp2d = interp2d(grid_x, grid_y, geom.chiArray(), kind='cubic')

xy_c = np.array([get_xy2(tth, chi, geom, plate_interp2d, chi_plate_interp2d) for tth, chi in zip(xy_tth, xy_chi)])

where get_xy2 is the following:

def get_xy2(tth, azi, geometry, tth_interp2d, chi_interp2d):
    """
    :param tth: tth value of point
    :param azi: azi value of point
    :param geometry: pyFAI geometry
    :type geometry: Geometry
    :param tth_interp2d: interpolator for the two theta array example:
                tth_interp2d = interp2d(grid_x, grid_y, geometry.twoThetaArray(img_shape), kind='cubic')
    :param chi_interp2d: interpolator for the two theta array
    :return: (x, y) tuple
    """

    # get start values from simple geometric consideration (without tilt calculation)
    fit2d_params = geometry.getFit2D()
    x_start = float(geometry.dist * np.tan(tth) * np.cos(azi) / geometry.pixel1 + fit2d_params["centerX"])
    y_start = float(geometry.dist * np.tan(tth) * np.sin(azi) / geometry.pixel1 + fit2d_params["centerY"])

    def min_func(x):
        tth_c = tth_interp2d(x[0], x[1])
        azi_c = chi_interp2d(x[0], x[1])
        chi_square = (tth - tth_c) ** 2 + (azi - azi_c) ** 2
        return chi_square * 1e5

    res = minimize(min_func, np.array([x_start, y_start]), method='L-BFGS-B')

    return res.x

this might be the fastest, when not having an exact equation.

Best wishes, Clemens

[kif]

Well, the result is not the same as out code parses all the pixel to search the nearest.

Your code assumes some "smooth" behavior in a larger neighborhood.

This runs in less than 10s:

import numpy, pyFAI, pyFAI.geometry
N = 5000
detector = pyFAI.detector_factory("Eiger1M")
geo = pyFAI.geometry.Geometry(dist=1e-3, detector=detector)
r = geo.array_from_unit(typ="center", unit=pyFAI.units.TTH_RAD, scale=False)
c = geo.chiArray()
ig = InvertGeometry(r, c)
x=numpy.random.random(size=N)*pi/2
y=numpy.random.random(size=N)*pi/2
res = [ig(a,b) for a,b in zip(x,y)]

Of course, it depends on the detector size and on the speed of the computer.