[WIP] Improving solid angle correction

[meghdadyazdi]

Trying to improve solid angle corrections. Current method might be a bit off or too simplified. This PR is a first step to make it better.

Please drop any papers, notes, or explanations that justify the changes or highlight the issues with the current approach. Would be helpful for context and future reference.

Feel free to tag others who should take a look.
@JJBekx @SAnsell @zdemat @fgjorup

[meghdadyazdi]

The new calculation of solid angle follows the approach described here in Eq. (7)

[JJBekx]

The attached presentation shows, on slide 16, the derivation of the solid angle correction as implemented in this merge request. It also illustrates how the old version was essentially the same, normalized to the (would-be) value of the center pixel.
lec_arul.pdf

[zdemat]

@JJBekx, you mean PONI-pixel (Pixel Of Normal Incidence)?

It seems this original normalization is quite common. Its advantage is that it is at the scale of 1 in most cases. The disadvantage is that if comparison of data for different detector positions or different detectors is required, one has to know the calibration (PONI-detector distance and pixel dimensions) to rescale to absolute intensity values.

Is it not better we leave the choice on the user? I.e. we add an option normalization: "to PONI-pixel" and "absolute"?

[zdemat]

Hej @JJBekx , @SAnsell has the point that in some cases (extra low angles) the numerical precision in the denominator (D^2 + px^2 + py^2), where (px^2 + py^2) << D^2, is insufficient. Would it be possible to do and expansion of (1+r)^(3/2) for r<<1? Of course there are extreme cases where e.g. (D^2 + px^2) << py^2, but they may be less common.

[JJBekx]

@JJBekx, you mean PONI-pixel (Pixel Of Normal Incidence)?

It seems this original normalization is quite common. Its advantage is that it is at the scale of 1 in most cases. The disadvantage is that if comparison of data for different detector positions or different detectors is required, one has to know the calibration (PONI-detector distance and pixel dimensions) to rescale to absolute intensity values.

Is it not better we leave the choice on the user? I.e. we add an option normalization: "to PONI-pixel" and "absolute"?

I agree that leaving it as an option to the user could be implemented, and would ease the transition of older versions of azint to newer ones that have this solid-angle correction.
However, I am rather partial to only keeping the "full" (absolute) correction, because it holds more information. Thinking back on our discussion regarding (meta-)data added in the Nexus formats, a situation we'd want to avoid is providing users with azint-processed data that loses data the user might need later (such as PONI-pixel relative normalization).
In this particular case, it'd be less worrisome, because it is rather straightforward in post-processing to revert between the two normalizations. But in general I'd opt for providing as much information, or as complete, to the user, as concisely as possible.

[JJBekx]

Hej @JJBekx , @SAnsell has the point that in some cases (extra low angles) the numerical precision in the denominator (D^2 + px^2 + py^2), where (px^2 + py^2) << D^2, is insufficient. Would it be possible to do and expansion of (1+r)^(3/2) for r<<1? Of course there are extreme cases where e.g. (D^2 + px^2) << py^2, but they may be less common.

It'd be possible of course, but I don't see how a truncated Taylor series would provide a more exact answer.
In this case, it'd essentially be a workaround, and I agree it could be necessary and implemented, but in general should we not opt for something more general to treat underflow like decimal or Kahan summation? Or is that too general for a single necessary quick fix you'd think?

[fgjorup]

I am joining this discussion a bit late, so I might have missed something, but to me it looks like
solid_angle = poni.dist * poni.det.pixel1 * poni.det.pixel2 / (poni.dist**2 + p1*p1 + p2*p2)**(3./2.)
is only valid for a perfectly perpendicular detector? I.e, when all pixels have the same z coordinate == poni.dist.
I would suggest something like
$SA = \frac{d_{poni}^3}{|\vec{v}|^3} $, where $d_{poni}$ is the point-of-nearest-incidence distance and $\vec{v}$ is the vector describing the position of each (sub)pixel relative to the sample position. This would of course give the relative correction, but since the distance and pixel dimensions are already available in the file, it is easy for users to calculate the absolute scale if they want.

I have a not-so-pretty example here:

def P_SA(xyz, pf=.99, pyfai_convention=True):
    """
    Combined polarization and solid-angle correction
    Parameters
    ----------
    xyz - ndarray (3,h,w)
          detector pixel coordinates as a (3,h,w) array of vectors
    pf  - float
          polarization factor (-1 vertical, +1 horizontal)

    pyfai_convention - bool
                       use pyFAI convention (-1 vertical, +1 horizontal)
                       else use Als-Nielsen convention (0 vertical, +1 horizontal)
                       pf = pf_pyfai*2-1
    Return
    ---------
    P_SA - ndarray (h,w)
           the combined pixel-wise polarization and solid-angle correction
    """
    if pyfai_convention:
        pf = pf*2-1
    
    r = np.sqrt(np.sum(xyz**2,axis=0))
    v = xyz/r
    p_h = v[0]**2*pf
    p_v = v[1]**2*(1-pf)
    p = 1.-(p_h+p_v)
    r_3 = r**(-3)
    SA = r_3/r_3.max()
    return p*SA

[zdemat]

Hej @fgjorup

do not be afraid. Our formula is correct :-) (px,py) are points on the detector plane. The z = poni.dist = D by definition. The PONI (point/pixel of normal incidence) is zero (px=0,py=0). So everything is perfectly orthogonal. That is the idea of PONI and the detector coordinate system related to PONI. All flat detectors are "perpendicular to the sample" in the single point and this PONI :-) Only the beam passing the sample point may hit the detector plane at oblique angle, i.e. not in the PONI. And the solid angle correction is about sample-detector geometry relation, not about beam.

[JJBekx]

Hej @fgjorup

do not be afraid. Our formula is correct :-) (px,py) are points on the detector plane. The z = poni.dist = D by definition. The PONI (point/pixel of normal incidence) is zero (px=0,py=0). So everything is perfectly orthogonal. That is the idea of PONI and the detector coordinate system related to PONI. All flat detectors are "perpendicular to the sample" in the single point and this PONI :-) Only the beam passing the sample point may hit the detector plane at oblique angle, i.e. not in the PONI. And the solid angle correction is about sample-detector geometry relation, not about beam.

Hi @zdemat,
(Edit: the part enclosed by ~~ in this comment is wrong, but kept for clarity)
~~
Though the formula for the solid-angle correction is correct, @fgjorup is right in pointing out that the poni.dist**2 in the denominator is wrong. The R³ in the denominator = sqrt(px_x² + px_y² + px_z²)³, but where the pixel coordinates have undergone the poni-translation and rotations so as to properly "see" a possibly skewed detector from the point-of-view of the sample-origin.

So it should be: solid_angle = poni.dist * poni.det.pixel1 * poni.det.pixel2 / (p3*p3 + p1*p1 + p2*p2)**(3./2.)
But then p3 should be passed as an argument in the function setup_corrections.

Perhaps it is worthy of note that this is not related to the choice of abs vs rel SA, and that this error is present in the main branch as well.
~~

Just for crystal-clarity:
$SA_{rel} = (D/R)^3$
$SA_{abs} = (dx * dy / D^2) * SA_{rel}$
with $D$ the PONI distance and $R$ the location vector of a pixel in the lab frame.

[fgjorup]

Hej @fgjorup

do not be afraid. Our formula is correct :-) (px,py) are points on the detector plane. The z = poni.dist = D by definition. The PONI (point/pixel of normal incidence) is zero (px=0,py=0). So everything is perfectly orthogonal. That is the idea of PONI and the detector coordinate system related to PONI. All flat detectors are "perpendicular to the sample" in the single point and this PONI :-) Only the beam passing the sample point may hit the detector plane at oblique angle, i.e. not in the PONI. And the solid angle correction is about sample-detector geometry relation, not about beam.

I see, the coordinates are in the "detector frame". As you see in my example I also calculate the polarization correction, so in that case I need to stay in the "lab frame", hence the confusion.

[JJBekx]

Hej @fgjorup
do not be afraid. Our formula is correct :-) (px,py) are points on the detector plane. The z = poni.dist = D by definition. The PONI (point/pixel of normal incidence) is zero (px=0,py=0). So everything is perfectly orthogonal. That is the idea of PONI and the detector coordinate system related to PONI. All flat detectors are "perpendicular to the sample" in the single point and this PONI :-) Only the beam passing the sample point may hit the detector plane at oblique angle, i.e. not in the PONI. And the solid angle correction is about sample-detector geometry relation, not about beam.

I see, the coordinates are in the "detector frame". As you see in my example I also calculate the polarization correction, so in that case I need to stay in the "lab frame", hence the confusion.

Oh, @zdemat @fgjorup, you're right. p1 and p2 are pixel positions after the poni translation, but without the rotations applied.
Considering we want the pixel location $R$ in the denominator, this seems odd to me? Why shouldn't we use the pos[:] array as defined in the transform function?