Ух ты ж ē-моē!!!
Таки добрался я наконец до своих баранов: вернулся к полиномам Цернике. И что оказалось: ведь отлично справляется метод наименьших квадратов! Т.е. в принципе, наверное, мучиться с ортогональными на кольце полиномами и не нужно!!!
Итак, вот собственно код, занимающийся поиском коэффициентов разложения «волнового фронта» по полиномам Цернике методом наименьших квадратов на любом наборе пробных точек:
double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation /* typedef struct { size_t size1; // rows (height) size_t size2; // columns (width) size_t tda; // data width (aligned) - data stride double * data; // pointer to block->data (matrix data itself) gsl_block * block; // block with matrix data (block->size is size) int owner; // ==1 if matrix owns 'block' (then block will be freed with matrix) } gsl_matrix; */ // Now allocate matrix: its Nth row is equation for Nth data point, // Mth column is Z_M gsl_matrix *M = gsl_matrix_calloc(Sz, pmax); // fill matrix with coefficients int x,y; size_t T = M->tda; for(x = 0; x < pmax; x++){ double norm, *Zcoeff = zernfunNR(x, Sz, P, &norm), *Zptr = Zcoeff; double *ptr = &(M->data[x]); for(y = 0; y < Sz; y++, ptr+=T, Zptr++) // fill xth polynomial *ptr = (*Zptr);// / norm; FREE(Zcoeff); } gsl_vector *ans = gsl_vector_calloc(pmax); gsl_vector_view b = gsl_vector_view_array(heights, Sz); gsl_vector *tau = gsl_vector_calloc(pmax); // min(size(M)) gsl_linalg_QR_decomp(M, tau); gsl_vector *residual = gsl_vector_calloc(Sz); gsl_linalg_QR_lssolve(M, tau, &b.vector, ans, residual); double *ptr = ans->data; int maxIdx = 0; double *Zidxs = MALLOC(double, pmax); for(x = 0; x < pmax; x++, ptr++){ if(fabs(*ptr) < Z_prec) continue; Zidxs[x] = *ptr; maxIdx = x; } gsl_matrix_free(M); gsl_vector_free(ans); gsl_vector_free(tau); gsl_vector_free(residual); if(lastIdx) *lastIdx = maxIdx; return Zidxs; }
Считает замечательно: для опорного набора коэффициентов (по которому генерировался «волновой фронт»)
double IdxsOri[] = {2., // смещение 1.1, -0.8, // наклон 5.5, -3.2, 0., // астигматизм, дефокус, аст. 6.8, 5.5, 0., 0.,// трилистник, кома 0., 0., 3.3, 1.4, 8.}; после восстановления получилось вот что:
2.000, 1.100, -0.800, 5.500, -3.200, 0.000, 6.800, 5.500, 0.000, 0.000, 0.000, 0.000, 3.300, 1.400, 8.000,
Я, честно говоря, даже не ожидал такой точности на всего-то 256 точках. И, кстати, считает быстро. Тайминги мне добавлять лень, но разложение по первым 6 степеням полиномов Цернике «волнового фронта» на 256 точек выполняется чуть ли не мгновенно.
Однако, рано я радовался: если максимальной степенью сделать 20, результаты получаются отвратительными. Вот где сказывается неортогональность, приводящая к неоднозначности решения.
В общем, реализовать полиномы, ортогональные на диске, мне все-таки придется! Либо же придется модифицировать задачу, постепенно подбирая коэффициенты: скажем, «скармливая» решателю по 3-4 очередных степени полинома (а в качестве значений ВФ брать невязки от предыдущего решения).
Итак, вот собственно код, занимающийся поиском коэффициентов разложения «волнового фронта» по полиномам Цернике методом наименьших квадратов на любом наборе пробных точек:
double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation /* typedef struct { size_t size1; // rows (height) size_t size2; // columns (width) size_t tda; // data width (aligned) - data stride double * data; // pointer to block->data (matrix data itself) gsl_block * block; // block with matrix data (block->size is size) int owner; // ==1 if matrix owns 'block' (then block will be freed with matrix) } gsl_matrix; */ // Now allocate matrix: its Nth row is equation for Nth data point, // Mth column is Z_M gsl_matrix *M = gsl_matrix_calloc(Sz, pmax); // fill matrix with coefficients int x,y; size_t T = M->tda; for(x = 0; x < pmax; x++){ double norm, *Zcoeff = zernfunNR(x, Sz, P, &norm), *Zptr = Zcoeff; double *ptr = &(M->data[x]); for(y = 0; y < Sz; y++, ptr+=T, Zptr++) // fill xth polynomial *ptr = (*Zptr);// / norm; FREE(Zcoeff); } gsl_vector *ans = gsl_vector_calloc(pmax); gsl_vector_view b = gsl_vector_view_array(heights, Sz); gsl_vector *tau = gsl_vector_calloc(pmax); // min(size(M)) gsl_linalg_QR_decomp(M, tau); gsl_vector *residual = gsl_vector_calloc(Sz); gsl_linalg_QR_lssolve(M, tau, &b.vector, ans, residual); double *ptr = ans->data; int maxIdx = 0; double *Zidxs = MALLOC(double, pmax); for(x = 0; x < pmax; x++, ptr++){ if(fabs(*ptr) < Z_prec) continue; Zidxs[x] = *ptr; maxIdx = x; } gsl_matrix_free(M); gsl_vector_free(ans); gsl_vector_free(tau); gsl_vector_free(residual); if(lastIdx) *lastIdx = maxIdx; return Zidxs; }
Считает замечательно: для опорного набора коэффициентов (по которому генерировался «волновой фронт»)
double IdxsOri[] = {2., // смещение 1.1, -0.8, // наклон 5.5, -3.2, 0., // астигматизм, дефокус, аст. 6.8, 5.5, 0., 0.,// трилистник, кома 0., 0., 3.3, 1.4, 8.}; после восстановления получилось вот что:
2.000, 1.100, -0.800, 5.500, -3.200, 0.000, 6.800, 5.500, 0.000, 0.000, 0.000, 0.000, 3.300, 1.400, 8.000,
Я, честно говоря, даже не ожидал такой точности на всего-то 256 точках. И, кстати, считает быстро. Тайминги мне добавлять лень, но разложение по первым 6 степеням полиномов Цернике «волнового фронта» на 256 точек выполняется чуть ли не мгновенно.
Однако, рано я радовался: если максимальной степенью сделать 20, результаты получаются отвратительными. Вот где сказывается неортогональность, приводящая к неоднозначности решения.
В общем, реализовать полиномы, ортогональные на диске, мне все-таки придется! Либо же придется модифицировать задачу, постепенно подбирая коэффициенты: скажем, «скармливая» решателю по 3-4 очередных степени полинома (а в качестве значений ВФ брать невязки от предыдущего решения).