Simple Imaging Polynomial, distortion.

[Zero Pointer]

Hi Dustin!

What is the mathematical or/and physical meaning of the coefficients A and B (A_p_q and B_p_q) in the struct sip_t? 

What model (of distortion) is used to describe distortion?

Can You recommend to read any book on this theme?

I`m diving in a reading of Greisen & Calabretta`s articles, an articte "The SIP Convention for Representing Distortion in FITS Image Headers", but I can not to understand it. 

I understand that  A_p_q and B_p_q are something fundamental and astrometry.net can itself compute both polynoms.

Thank You!

[Danko Dnevic]

Hi,

 I'll try to give an answer.

SIP is a simplified bivariate (x,y or u,v in pairs) polynomial, reduced to terms where p+q <= order (A order or B order). A and B coefficients are generated, for each axis separately (A for x and B for y, but using BOTH input coordinates for correction calculation), by least squares fit to expected true coordinates (yielded from catalog). This is achieved in AN by using libraries, like GSL and CBLAS, but anyone can do the same in Matlab/Octave, or Python+NumPy. Kudos to Dustin, for a really great implementation in AN.

Here is 2008. paper that describes the concept:

https://fits.gsfc.nasa.gov/registry/sip/SIP_distortion_v1_0.pdf

BR,

DD

[Danko Dnevic]

I mistakenly left link for the paper you already mention...

[Danko Dnevic]

Here is some Octave/Matlab code, you can use to play with, to see the purpose of A and B:

CD is 2x2 matrix

 crpix iand arval are vectors containing coordinates of the image center. and corresponding Ra0 and Dec0

For given A_order A matrix is (A_order+1)x(A_order+1) matrix, but populated just with relevant SIP coeffs. i.e. for A_order==2:

A=[A_0_0, A_0_1, A_0_2; ...

    A_1_0, A_1_1, 0; ...

    A_2_0, 0, 0]

Same for B...

========================================================================

function [Ra, Dec, x_c, y_c] = pix2RaDec(CD, crpix, crval, A, B, u, v)

    x_i = u - crpix(1);

    y_i = v - crpix(2);

    Ra0 = crval(1);

    Dec0 = crval(2);

    % A and B are the distortion matrices for x and y respectively

    x_c = polyval2d_internal(x_i, y_i, A);

    y_c = polyval2d_internal(x_i, y_i, B);

    x_d = x_i + x_c;

    y_d = y_i + y_c;

    stdXY = CD * [x_d; y_d];

    X = deg2rad(stdXY(1));

    Y = deg2rad(stdXY(2));

    Ra  = mod(Ra0 + atand( X / (cosd(Dec0) - Y * sind(Dec0) )),360);

    Dec = asind ( (sind(Dec0)+Y*cosd(Dec0))/sqrt(1+X*X+Y*Y) );

endfunction

function val = polyval2d_internal(x, y, coeffs)

    % Evaluate a 2D polynomial with given coefficients at points (x, y)

    val = 0;

    [rows, cols] = size(coeffs);

    for i = 1:rows

        for j = 1:cols

            val = val + coeffs(i,j) * (x.^(i-1)) .* (y.^(j-1));

        endfor

    endfor

endfunction

Dana utorak, 9. srpnja 2024. u 08:13:23 UTC+2 korisnik zeropo...@gmail.com napisao je:

[Zero Pointer]

Hello!

Thank You for an answer and Matlab/Octave code! 

You are prevented my question about axes and matrices A and B)

What is the "AN"? Is it an astrometry.net?

Have a nice day!

[Danko Dnevic]

Hi,

you are welcome. 

Yes, by AN I mean astrometry.net

If you have any further questions, just ask, I'll try to answer...

BR,

DD

[Zero Pointer]

Hi everybody, hi DD!

There are coefficients A and B in the struct sip_t. They are makes some forward transform. 

Transform`s result  is convergence of two arrays of points: points from  the star catalogue and points from the frame: such every lenght between point from the frame and point from  the cataloge tends to zero (see attachment to message: there are a picture with flag "--no-tweak", and a picture without flag in the attachment).

With a flag "--no-tweak" at running of solve-field:

Without flag at running of solve-field:

What is this transform?

This transform is dependments from quad and quad`s position in the frame. Transform is uniqle for every frame.

This is not a lens distortion (I think so, because we can not use once computed coefficients for other frames).

Thank You!

[Danko Dnevic]

Hi,

I hope that I will correctly describe some AN internals. If something is wrong @Dustin can correct me... 

Think about quads as AN tool to recognize sources on the image, used even before SIP is "tweaked". Once AN has a valid guess, where center 'CRPIX (plate coordinates)' is looking at 'CRVAL' (Ra0 and Dec0) and how image is scaled (degrees per pixel) and rotated (2x2 CD matrix), field is matched with sources from the catalog. Only then, if tweak is used, SIP is generated, to compensate differences from "apparent" pixel positions and tangent plane positions, calculated from the catalog. SIP doesn't compensate just lens distortion, but for all differences within its polynomial reach (order). That is good and bad (see why below).

For narrow fields, taken above 45 degrees elevation, A and B can be reused as good distortion approximation, because it only "move" pixels based on calculated plane positions, ignoring star coordinates. Look at the Matlab code, you can see that only after CD and standard coordinate X,Y transformation, using Ra0 and Dec0 offset, pixel coordinates are related to Ra and Dec.

I mentioned that SIP doesn't care about distortion nature, is it intrinsic or extrinsic. it is a hammer and every distortion is a nail. If there is atmospheric refraction difference, it will try to compensate, rendering these coefficients unusable, but for that camera position.

BR,

DD

[Zero Pointer]

Hi everybody!

DD, thank You for Your detailed response.

Thank for a metafor "a hammer and a nail"

Have a nice day!