И что я сразу о GSL не вспомнил?
Бился-бился я с ортогональными на кольце полиномами Цернике, но что-то сразу не подумал, что можно ведь другим путем пойти: ведь можно (пусть и неоднозначно) разложить функцию даже по базису не ортогональных полиномов! А для этого всего-то достаточно использовать метод наименьших квадратов.
И наткнулся я на эту идею в этой статье:
@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 #include #include 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 уменьшится.
Эксперименты с разложением волнового фронта по полиномам Цернике на непрямоугольной сетке (точнее - для реальных координат отверстий в диафрагме Гартманна) оставлю на следующую неделю, т.к. после обеда надо помыть машину и ехать в командировку в Ставрополь, слушать защиту наших магистрантов (завтра с утра). Так что в понедельник-вторник выясню, подойдет ли этот метод для разложения волнового фронта на кольце (да еще и на неравномерной сетке), а также - для разложения градиента этого самого фронта. Если все будет ОК, можно будет приступить к анализу реальных данных (во всяком случае, по моделированию статью можно будет сразу написать, а применительно к реальным данным - чуть позже; все равно пока нет у нас ни одной нормальной пары гартманнограмм, все время с погодой не везло).
И наткнулся я на эту идею в этой статье:
@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 #include #include 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 уменьшится.
Эксперименты с разложением волнового фронта по полиномам Цернике на непрямоугольной сетке (точнее - для реальных координат отверстий в диафрагме Гартманна) оставлю на следующую неделю, т.к. после обеда надо помыть машину и ехать в командировку в Ставрополь, слушать защиту наших магистрантов (завтра с утра). Так что в понедельник-вторник выясню, подойдет ли этот метод для разложения волнового фронта на кольце (да еще и на неравномерной сетке), а также - для разложения градиента этого самого фронта. Если все будет ОК, можно будет приступить к анализу реальных данных (во всяком случае, по моделированию статью можно будет сразу написать, а применительно к реальным данным - чуть позже; все равно пока нет у нас ни одной нормальной пары гартманнограмм, все время с погодой не везло).