Ух ты ж ē-моē!!!

[eddy_em]

Таки добрался я наконец до своих баранов: вернулся к полиномам Цернике. И что оказалось: ведь отлично справляется метод наименьших квадратов! Т.е. в принципе, наверное, мучиться с ортогональными на кольце полиномами и не нужно!!!

Итак, вот собственно код, занимающийся поиском коэффициентов разложения «волнового фронта» по полиномам Цернике методом наименьших квадратов на любом наборе пробных точек:

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 очередных степени полинома (а в качестве значений ВФ брать невязки от предыдущего решения).