Разложение волнового фронта по полиномам Цернике в сях

[eddy_em]

У меня появилось желание добить-таки наконец свой fitsview-hartmann, а для этого нужно добавить туда декомпозицию нормалей к волновому фронту по ортонормированному на кольце базису векторных полиномов. Для этого сначала надо реализовать разложение/восстановление волнового фронта по полиномам Цернике. Кстати, это еще пригодится мне в модели зеркала (можно будет помимо получения маски искажений формы зеркала из файла-изображения добавить возможность указания величины этих искажений коэффициентами Цернике).

Все целиком приводить не буду для экономии места. Вскорости (после "причесывания") код можно будет увидеть в моих сниппетах.
Сами полиномы Цернике взяты отсюда (т.к. каждый представляет их по-своему, я решил решающим словом выбрать википедию). Перевод нотации Нолля в n/m сделан в соответствии с более распространенным вариантом: 0 -> 0,0; 1 -> 1,-1; 2 -> 1,1; 3 -> 2,-2; 4 -> 2,0; 5 -> 2,2 (т.е. по линейному возрастанию m, а не абсолютному возрастанию, как в википедии). Код для перевода совершенно простой:

/**
 * Convert polynomial order in Noll notation into n/m
 * @param p (i) - order of Zernike polynomial in Noll notation
 * @param N (o) - order of polynomial
 * @param M (o) - angular parameter
 */
void convert_Zidx(int p, int *N, int *M){
	int n = (int) floor((-1.+sqrt(1.+8.*p)) / 2.);
	*M = (int)(2.0*(p - n*(n+1.)/2. - 0.5*(double)n));
	*N = n;
}

В вычислениях нет ничего хитрого: это обычное решение квадратного уравнения через дискриминант. Дело в том, что для заданного p значение n есть ни что иное, как p·(p+1)/2 (формула суммы арифметической прогрессии). А арифметическая прогрессия взялась здесь оттуда, что для каждого n количество возможных значений m (а, следовательно, и полиномов степени n) равно n+1. Ну, а вычисление m, думаю, и так понятно.
Вычисление полиномов Цернике пока никак не оптимизировано. Вот, собственно, эта функция:

/**
 * Calculate value of Zernike polynomial on rectangle matrix WxH pixels
 * Center of matrix will be zero point
 * Scale will be set by max(W/2,H/2)
 * @param n - order of polynomial (max: 100!)
 * @param m - angular parameter of polynomial
 * @param W - width of output array
 * @param H - height of output array
 * @param norm (o) - (if !NULL) normalize factor
 * @return NULL-terminated array of Zernike polynomials on given matrix
 */
double *zernfunR(int n, int m, int W, int H, double *norm){
	double Z = 0., *Zarr = NULL;
	bool erparm = false;
	if(W < 2 || H < 2)
		errx(1, "Sizes of matrix must be > 2!");
	if(n > 100)
		errx(1, "Order of Zernike polynomial must be <= 100!");
	if(n < 0) erparm = true;
	if(n < iabs(m)) erparm = true; // |m| must be <= n
	int d = n - m;
	if(d % 2) erparm = true; // n-m must differ by a prod of 2
	if(erparm)
		errx(1, "Wrong parameters of Zernike polynomial (%d, %d)", n, m);
	if(!FK) build_factorial();
	double Xc = (W - 1.) / 2., Yc = (H - 1.) / 2., Rnorm = fmax(Xc, Yc); // coordinate of circle's middle
	int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, N = n+1, w = (W+1)/2, h = (H+1)/2, S=w*h;
	double **Rpow = MALLOC(double*, N); // powers of R
	Rpow[0] = MALLOC(double, S);
	for(j = 0; j < S; j++) Rpow[0][j] = 1.; // zero's power
	double *R = MALLOC(double, S);
	// first - fill array of R
	double xstart = (W%2) ? 0. : 0.5, ystart = (H%2) ? 0. : 0.5;
	for(j = 0; j < h; j++){
		double *pt = &R[j*w], x, ww = w;
		for(x = xstart; x < ww; x++, pt++){
			double yy = (j + ystart)/Rnorm, xx = x/Rnorm;
			*pt = sqrt(xx*xx+yy*yy);
		}
	}
	for(i = 1; i < N; i++){ // Rpow - is quater I of cartesian coordinates ('cause other are fully simmetrical)
		Rpow[i] = MALLOC(double, S);
		double *rp = Rpow[i], *rpo = Rpow[i-1];
		for(j = 0; j < h; j++){
			int idx = j*w;
			double *pt = &rp[idx], *pto = &rpo[idx], *r = &R[idx];
			for(k = 0; k < w; k++, pt++, pto++, r++)
				*pt = (*pto) * (*r); // R^{n+1} = R^n * R
		}
	}
	int SS = W * H;
	// now fill output matrix
	Zarr = MALLOC(double, SS); // output matrix W*H pixels
	double ZSum = 0.;
	for(j = 0; j < H; j++){
		double *Zptr = &Zarr[j*W];
		double Ryd = fabs(j - Yc);
		int Ry = w * (int)Ryd; // Y coordinate on R matrix
		for(i = 0; i < W; i++, Zptr++){
			Z = 0.;
			double Rxd = fabs(i - Xc);
			int Ridx = Ry + (int)Rxd; // coordinate on R matrix
			if(R[Ridx] > 1.) continue; // throw out points with R>1
			// calculate R_n^m
			for(k = 0; k <= iup; k++){ // Sum
				double z = (1. - 2. * (k % 2)) * FK[n - k]        //       (-1)^k * (n-k)!
					/(//----------------------------------- -----   -------------------------------
						FK[k]*FK[(n+m_abs)/2-k]*FK[(n-m_abs)/2-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)!
					);
				Z += z * Rpow[n-2*k][Ridx];  // *R^{n-2k}
			}
			// normalize
			double eps_m = (m) ? 1. : 2.;
			Z *= sqrt((double)(n+1) / M_PI / eps_m);
			double theta = atan2(j - Yc, i - Xc);
			// multiply to angular function:
			if(m){
				if(m > 0)
					Z *= cos(theta*(double)m_abs);
				else
					Z *= sin(theta*(double)m_abs);
			}
			*Zptr = Z;
			ZSum += Z*Z;
		}
	}
	if(norm) *norm = ZSum;
	// free unneeded memory
	FREE(R);
	for(i = 0; i < N; i++) FREE(Rpow[i]);
	FREE(Rpow);
	return Zarr;
}

Я "прибил гвоздями" ограничение на степень полинома, равное 100 (да и то, это много: уже где-то за 60 получаются сплошные nan'ы). Данная функция строит по заданным n и m полином Цернике на равномерной прямоугольной сетке размером WxH (если сетка не квадратная, единицей считается максимум из W/2 и H/2). В параметр norm (если он !NULL) помещается сумма квадратов полинома на единичном круге внутри этой матрицы. Это необходимо для дальнейшей нормализации. Функция возвращает динамически выделенную память с заданным полиномом.

Пока я писал этот код, столкнулся вот с чем: "волновой фронт" упорно не хотел восстанавливаться моими полиномами. Не сразу до меня дошло, что проблема кроется в том, что оригинальные полиномы вычисляются для непрерывных функций, да и в интеграле там есть r. Вот и пришлось для простоты добавить вычисление нормирующего коэффициента — суммы квадратов полинома по узлам заданной сетки. Зато этот вариант позволит вычислять полиномы даже по неравномерной сетке.

Для вычисления полинома в нотации Нолля служит функция-обертка:

/**
 * Zernike polynomials in Noll notation for square matrix
 * @param p - index of polynomial
 * @other params are like in zernfunR
 * @return zernfun
 */
double *zernfunNR(int p, int W, int H, double *norm){
	int n, m;
	convert_Zidx(p, &n, &m);
	return zernfunR(n,m,W,H,norm);
}

Декомпозиция "волнового фронта" выполняется просто: в пределах заданной степени полиномов Цернике производятся итерации по вычислению очередного коэффициента. При этом произведение величины "фронта" в данной точки на величину полинома делится на нормирующий коэффициент (чтобы учесть дискретность). Здесь же вычисляется максимальный номер ненулевого коэффициента. Для того, чтобы упростить вычисления (т.е. уменьшить коэффициенты высших степеней), я использовал промежуточное изображение, в которое копировалось изначальное, а затем на каждой итерации из него вычиталось произведение текущего полинома на текущий коэффициент. Вот и сама функция:

/**
 * Zernike decomposition of image 'image' WxH pixels
 * @param Nmax (i) - maximum power of Zernike polinomial for decomposition
 * @param W, H (i) - size of image
 * @param image(i) - image itself
 * @param Zsz  (o) - size of Z coefficients array
 * @param lastIdx (o) - (if !NULL) last non-zero coefficient
 * @return array of Zernike coefficients
 */
double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx){
	int i, SS = W*H, pmax, maxIdx = 0;
	double *Zidxs = NULL, *icopy = NULL;
	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
	printf("pmax = %d\n", pmax);
	Zidxs = MALLOC(double, pmax);
	icopy = MALLOC(double, W*H);
	memcpy(icopy, image, W*H*sizeof(double)); // make image copy to leave it unchanged
	*Zsz = pmax;
	for(i = 0; i < pmax; i++){ // now we fill array
		double norm, *Zcoeff = zernfunNR(i, W, H, &norm);
		int j;
		double *iptr = icopy, *zptr = Zcoeff, K = 0.;
		for(j = 0; j < SS; j++, iptr++, zptr++)
			K += (*zptr) * (*iptr) / norm; // multiply matrixes to get coefficient
		Zidxs[i] = K;
		if(fabs(K) < Z_prec){
			Zidxs[i] = 0.;
			continue; // there's no need to substract values that are less than our precision
		}
		maxIdx = i;
		iptr = icopy; zptr = Zcoeff;
		for(j = 0; j < SS; j++, iptr++, zptr++)
			*iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values
		FREE(Zcoeff);
	}
	if(lastIdx) *lastIdx = maxIdx;
	FREE(icopy);
	return Zidxs;
}

Для того, чтобы "не связываться" со слишком малыми коэффициентами, я ввел ограничение на их амплитуду: по умолчанию это 10^-6, но при помощи функции set_prec(int val) можно установить требуемую точность.
Восстановление фронта по значениям коэффициентов Цернике — вообще простейшая задача:

/**
 * Zernike restoration of image WxH pixels by Zernike polynomials coefficients
 * @param Zsz  (i) - number of elements in coefficients array
 * @param Zidxs(i) - array with Zernike coefficients
 * @param W, H (i) - size of image
 * @return restored image
 */
double *Zcompose(int Zsz, double *Zidxs, int W, int H){
	int i, SS = W*H;
	double *image = MALLOC(double, SS);
	for(i = 0; i < Zsz; i++){ // now we fill array
		double K = Zidxs[i];
		if(fabs(K) < Z_prec) continue;
		double *Zcoeff = zernfunNR(i, W, H, NULL);
		int j;
		double *iptr = image, *zptr = Zcoeff;
		for(j = 0; j < SS; j++, iptr++, zptr++)
			*iptr += K * (*zptr); // add next Zernike polynomial
		FREE(Zcoeff);
	}
	return image;
}

Comments

Полиномы Цернике.

[asm786.livejournal.com]

Привет. Имел опыт работы с данными в области дескрипторов для описания формы изображений, там несколько проще, так как существует связь между полиномами и центральными моментами.
Инвариант Цернике в общем виде:


Связь с моментами:

Пример для p,q \in [0,4]:

На картинке продемонстрирована матрица из инвариантов Цернике, порядка 4. Значения u - центральные моменты изображения, вычисленные в единичной окружности. Изображение не нужно масштабировать, чтоб удовлетворить условию "внутри единичной окружности" - достаточно нормировать матрицу центральных моментов.

Если интересно, могу поделиться скриптом символьного вычисления таких вот матриц. Сорри за форматирование :)

Re: Полиномы Цернике.

[eddy-em.livejournal.com]

Привет!
А как эта связь между моментами и коэффициентами Цернике работает, если координатная сетка опорных точек неравномерная?

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

А скрипт - это интересно. Если не жалко, поделись, пожалуйста.

[asm786.livejournal.com]

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

Скрипт на языке Mathematica, https://dl.dropboxusercontent.com/u/108372303/ZernikeEvaluation.nb

В формулах выше убрана нормализация (те для получения правильных величин нужно все элементы матрицы разделить на Pi), ну и используется факт, что для нормированных центральных моментов величины \mu_01 и \mu_10 равны нулю.

[eddy-em.livejournal.com]

С интегралами на неравномерной сетке похуже: площадь-то вокруг точек постоянно меняется.
За скрипт спасибо: огромный такой и страшный ☺ Посмотрю (правда, в Mathematica я полный нуль).