7.1.3.2 Two dimensional polynomial fitting ¶

In Least squares fitting an introduction was given for fitting in Gnuastro and one-dimensional fitting was reviewed in One dimensional polynomial fitting (which should be read before this in a first read). For a 2D polynomial fit, two independent variables are needed. The two independent variables can be given either as an image or as table columns. We’ll start with an image here, and end with using table columns.

Let’s use the Arithmetic program to generate an image from the following 2D polynomial: \(Y=5+10X_2+X_1X_2\). In a FITS image, \(X_1\) is the first FITS axis (horizontal) and \(X_2\) is the second FITS axis (vertical). For the details of the operators, see Arithmetic operators.

$ astarithmetic 100 100 2 makenew indexonly set-i \
                i 100 % 1 + set-X1 \
                i 100 / 1 + set-X2 \
                5 10 X2 x + X1 X2 x + f64 \
                --output=pixels.fits
$ astscript-fits-view pixels.fits

The gradient of the 2D polynomial is clearly visible in the image. You can directly give such an image to the statistics program for fitting in a very similar way One dimensional polynomial fitting as shown below (with some outputs printed):

$ aststatistics pixels.fits --fit=polynomial --fitmaxpower=2
Statistics (GNU Astronomy Utilities) 0.24
-------
Fitting results (remove extra info with '--quiet' or '-q)
  Input file:    pixels.fits (hdu: 1) with 10000 non-blank rows.
  X column: Horizontal position of non-NAN pixels.
  Y column: Vertical position of non-NAN pixels.

Fit function [2d, Y=f(X1,X2)]: Y = c0 + c1.X1 + c2.X2 + c3.X1^2
+ c4.X1.X2 + c5.X2^2 + c6.X1^3 + c7.X1^2.X2 + c8.X1.X2^2 + c9.X2^3
+ ... + cn.X1^(j).X2^(d-j)
  N:  2
  c0:  +5.000000000073776e+00
  c1:  +2.460254222569347e-13
  c2:  +9.999999999995385e+00
  c3:  -2.773822838086915e-15
  c4:  +1.000000000000002e+00
  c5:  +4.474025316891783e-14

Covariance matrix:
  +3.318964094598024e-24 -6.837076956939867e-26 -6.837076956940090e-26
  +3.786011000755381e-28 +4.497322158473032e-28 +3.786011000755588e-28
  ...

Except for the functional form of the the 2D fit and the number of constants, all the outputs are similar to what was introduced in One dimensional polynomial fitting. In this context, the extremely small values (\(\sim10^{-14}\)) are just floating point errors and can be considered equal to zero. We therefore see that the three constants that we had originally set and nicely determined. But we are rarely lucky enough to have a value for the full plane like above. Let’s simulate such a situation by adding a noisy dataset (n, with each pixel having a random value from a Gaussian distribution with \(\sigma=100\)). In the command below, we’ll use that noisy dataset to randomly remove (set to NaN/blank) all pixels that have a value below 50 in it.

$ astarithmetic 100 100 2 makenew indexonly set-i \
                i 100 % 1 + set-X1 \
                i 100 / 1 + set-X2 \
                5 10 X2 x + X1 X2 x + f64 set-p \
                i 0 constant 100 mknoise-sigma set-n \
                p n 50 lt nan where --output=pixels.fits

$ astscript-fits-view pixels.fits

$ aststatistics pixels.fits --fit=polynomial --fitmaxpower=2

With the second command, you can see that many of the pixels have become NaN and from the third command, you can see that the fit is still good. You can try decreasing the number of useful (non-NaN) pixels (by increasing the 50 threshold above) and see how robust the fit is to fewer pixels.

As mentioned at the start of this section, you can also feed your input data as table columns. First let’s generate some table columns from the pixels.fits image we made above: we will feed the resulting table to statistics instead of table. The top command below will generate a 3 column table with the same number of rows as pixels in the image (one row per pixel). The first and second columns will be the first (horizontal) and second (vertical) positions of the pixel and the third column will be the pixel value.

$ astarithmetic pixels.fits set-i \
                i 1 size    set-width \
                i indexonly set-ind \
                ind width % 1 + u32 set-X1 \
                ind width / 1 + u32 set-X2 \
                i to-1d X2 to-1d X1 to-1d \
                --writeall --output=pixels-tab-raw.fits

If you inspect this table, you will notice that many rows are NaN (as expected!). To remove the NaN rows, you can use Table’s --noblank option:

$ asttable pixels-tab-raw.fits --noblank=_all \
           --output=for-2d-fit.fits

Now that we have the two independent variables (horizontal and vertical positions) along with the function/pixel value, we are ready to call Statistics for a 2D polynomial fit:

$ aststatistics for-2d-fit.fits --fit=polynomial -c1,2,3 \
                --fitmaxpower=2

You will notice that the fitted parameters are the same as the case where you gave an image. In fact, when given an image, Statistics does a similar thing to what we did here to extract the positions and values of all the non-NaN pixels into a table and then feed them to Gnuastro’s internal fitting library.