И что я сразу о GSL не вспомнил?

[eddy_em]

Бился-бился я с ортогональными на кольце полиномами Цернике, но что-то сразу не подумал, что можно ведь другим путем пойти: ведь можно (пусть и неоднозначно) разложить функцию даже по базису не ортогональных полиномов! А для этого всего-то достаточно использовать метод наименьших квадратов.

И наткнулся я на эту идею в этой статье:

@ARTICLE{2010ApOpt..49.6489M,
   author = {{Mahajan}, V.~N. and {Aftab}, M.},
    title = "{Systematic comparison of the use of annular and Zernike circle polynomials for annular wavefronts}",
  journal = {\ao},
     year = 2010,
    month = nov,
   volume = 49,
    pages = {6489},
      doi = {10.1364/AO.49.006489},
   adsurl = {http://adsabs.harvard.edu/abs/2010ApOpt..49.6489M},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

Вот простенький тестовый код:

#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include <math.h>

void fillmatrix(gsl_matrix *m, double *x){
	int i, j;
	size_t W = m->size2, H = m->size1, T = m->tda;
	double *xptr = x;
	for(i = 0; i < H; i++, xptr++){
		double *ptr = &(m->data[i*T]);
		for(j = 0; j < W; j++, ptr++)
			*ptr = sin((double)(j+1) * (*xptr));
	}
}

int main(int argc, char **argv){
	double b_data[] = {0.37, 1.07, 0.95, 0.43, -0.27};
	double xx[] = {0.1, 0.34, 0.76, 0.93, 1.12};
	gsl_matrix *M = gsl_matrix_calloc(5, 4);
	fillmatrix(M, xx);
	gsl_vector *x = gsl_vector_calloc(4);
	gsl_vector_view b = gsl_vector_view_array (b_data, 5);
	gsl_vector *tau = gsl_vector_calloc(4); // min(size(M))
	gsl_linalg_QR_decomp(M, tau);
	gsl_vector *residual = gsl_vector_calloc(5);
	gsl_linalg_QR_lssolve(M, tau, &b.vector, x, residual);
	printf ("x = \n");
	gsl_vector_fprintf(stdout, x, "%g");
	printf ("\nresidual = \n");
	gsl_vector_fprintf(stdout, residual, "%g");
	gsl_matrix_free(M);
	gsl_vector_free(x);
	gsl_vector_free(tau);
	gsl_vector_free(residual);
	return 0;
}

Здесь я просто раскладываю функцию 1.25·sin(3·x) (округленную до двух знаков после запятой) на неравномерной сетке аргумента x по синусному базису (с первой по четвертую гармоники).
Результат получается вполне приличным:

gcc -lgsl -lgslcblas -lm G.c -o vect && ./vect 
x = 
0.130783
-0.143849
1.32437
-0.0124236

residual = 
-0.00101798
0.000479872
-0.00021396
0.00014092
-2.66255e-05

(если результаты округлять до четвертого знака, то получается более точно: 1.27 вместо 1.32). Понятно, что при увеличении количества точек раз в 100 погрешность раз в 10 уменьшится.
Эксперименты с разложением волнового фронта по полиномам Цернике на непрямоугольной сетке (точнее — для реальных координат отверстий в диафрагме Гартманна) оставлю на следующую неделю, т.к. после обеда надо помыть машину и ехать в командировку в Ставрополь, слушать защиту наших магистрантов (завтра с утра). Так что в понедельник-вторник выясню, подойдет ли этот метод для разложения волнового фронта на кольце (да еще и на неравномерной сетке), а также — для разложения градиента этого самого фронта. Если все будет ОК, можно будет приступить к анализу реальных данных (во всяком случае, по моделированию статью можно будет сразу написать, а применительно к реальным данным — чуть позже; все равно пока нет у нас ни одной нормальной пары гартманнограмм, все время с погодой не везло).