Add wavefront decomposition

benchmark.c

@@ -91,7 +91,7 @@ wavefront *calc_surfaceR(int sz, double *Zidxs, int znum){
  * @param znum   - amount of Zidx
  */
 wavefront *calc_surfaceRS(int sz, polar *points, double *Zidxs, int znum){
-	if(!Zidxs || !points || znum < 1 || sz < 1) return NULL;
+	if(!Zidxs || !points || znum < 1 || sz < 1) ERRX(_("Bad arguments of calc_surfaceRS:%s"));
 	wavefront *F = MALLOC(wavefront, 1);
 	F->size = sz;
 	F->coordinates = MALLOC(polar, sz);
@@ -119,7 +119,7 @@ polar *calc_BTA_Hpoints(int *sz){
 	pt->r = r[2] / rmax; pt->theta = 24.5 * ray_step + theta0;
 	++pt;
 	pt->r = r[7] / rmax; pt->theta = 29.5 * ray_step + theta0;
-	if(sz) *sz = 256;
+	if(sz) *sz = 258;
 	return pts;
 }
 
@@ -131,6 +131,39 @@ static char *bench_names[] = {
 	N_("LS gradient decomposition (LS_gradZdecomposeR)")
 };
 
+/**
+ * calculate difference of two wavefronts
+ * @param WF1, WF2 (i) - wavefronts
+ * @param sz (i) - size of WF1 and WF2
+ * @param std  (o) - RMS of difference
+ * @param maxd (o) - maximal difference value
+ * @param maxd (o) - maximal relative difference value
+ */
+void WF_diff(double *WF1, double *WF2, int sz, double *std, double *maxd, double *maxdr){
+	int i;
+	double _maxd = 0., _maxdr = 0., sum = 0., sum2 = 0.;
+	for(i = 0; i < sz; ++i){
+		double z = WF1[i];
+		double diff = z - WF2[i];
+		sum += diff;
+		sum2 += diff*diff;
+		diff = fabs(diff);
+		if(_maxd < diff){
+			_maxd = diff;
+			z = fabs(z);
+			if(z > WF_EPSILON){
+				_maxdr = diff / z;
+			//	diff /= z;
+			//	if(_maxdr < diff) _maxdr = diff;
+			}
+		}
+	}
+	sum /= sz; sum2 /= sz;
+	if(std) *std = sqrt(sum2 - sum*sum);
+	if(maxd) *maxd = _maxd;
+	if(maxdr) *maxdr = _maxdr;
+}
+
 /**
  * calculate errors and fill WF_stat fields for single measurement
  * @param stat    - output stat structure
@@ -143,31 +176,21 @@ static char *bench_names[] = {
 static void fill_stat(WF_benchmark *bench, WF_bench_type type, double *Zidxs, double *Zidxs0, wavefront *WF0, int verbose){
 	double *zdata0 = WF0->zdata;
 	polar *P = WF0->coordinates;
-	size_t i, sz = WF0->size;
+	int sz = WF0->size;
 	WF_stat *stat = &bench->stat[type];
 	int znum = bench->Znum;
 	// test errors by wavefront
 	double *zdata = ZcomposeR(znum, Zidxs, sz, P);
-	double maxd = 0., sum = 0., sum2 = 0.;
-	for(i = 0; i < sz; ++i){
-		double z = zdata0[i];
-		double diff = z - zdata[i];
-		sum += diff;
-		sum2 += diff*diff;
-		if(fabs(z) > WF_EPSILON){
-			diff = fabs(diff / z);
-			if(maxd < diff) maxd = diff;
-		}
-	}
-	stat->max_dWF = maxd;
-	sum /= sz; sum2 /= sz;
-	stat->WF_std = sqrt(sum2 - sum*sum);
+	//for(i = 0; i < znum; ++i) printf("ZC_%d\t%10.3f\t%10.3f\n", i, Zidxs0[i], Zidxs[i]);
+	//for(i = 0; i < sz; ++i) printf("Point_%d\t%10.3f\t%10.3f\n", i, zdata0[i], zdata[i]);
+	double maxdr;
+	WF_diff(zdata0, zdata, sz, &(stat->WF_std), &(stat->max_dWF), &maxdr);
 	stat->ZC_sum = MALLOC(double, znum);
 	stat->ZC_sum2 = MALLOC(double, znum);
 	stat->max_dZC = MALLOC(double, znum);
 	if(verbose){
 		printf("%s\n", bench_names[type]);
-		printf("Wavefront error: STD=%g, |max err|=%.3f%%; Zernike:\n#\t   val0   \t    val   \tdiff%%\n", stat->WF_std, maxd);
+		printf("Wavefront error: STD=%g, |max err| = %.3f, |max err rel|=%.3f%%; Zernike:\n#\t   val0   \t    val   \tdiff%%\n", stat->WF_std, stat->max_dWF, maxdr);
 	}
 	for(int p = 0; p < znum; ++p){
 		double diff = Zidxs0[p] - Zidxs[p];
@@ -202,7 +225,6 @@ static WF_benchmark *benchmark_step(double *Zidxs0, int znum, int verbose){
 	MESG(_("1. Scattered points for BTA hartmannogram"));
 	P = calc_BTA_Hpoints(&sz);
 	WF = calc_surfaceRS(sz, P, Zidxs0, znum);
-	FREE(P);
 	zdata = WF->zdata;
 	// ZdecomposeR
 	Zidxs = ZdecomposeR(N, sz, WF->coordinates, zdata, &Zsz, NULL);
@@ -216,12 +238,23 @@ static WF_benchmark *benchmark_step(double *Zidxs0, int znum, int verbose){
 	Zidxs = QR_decompose(N, sz, WF->coordinates, zdata, &Zsz, NULL);
 	fill_stat(bench, Scatter_QR, Zidxs, Zidxs0, WF, verbose);
 	FREE(Zidxs);
-	// directGradZdecomposeR
+	free_wavefront(&WF);
 	Zidxs0[0] = 0.; // decomposition by gradients cannot know anything about zero coefficient
+	WF = calc_surfaceRS(sz, P, Zidxs0, znum); // recalculate new wavefront
+	// directGradZdecomposeR
 	point *grads = directGradZcomposeR(znum, Zidxs0, sz, WF->coordinates);
 	Zidxs = directGradZdecomposeR(N, sz, WF->coordinates, grads, &Zsz, NULL);
 	fill_stat(bench, Scatter_grad, Zidxs, Zidxs0, WF, verbose);
 	FREE(Zidxs);
+	// LS_gradZdecomposeR
+	Zidxs = LS_gradZdecomposeR(N, sz, WF->coordinates, grads, &Zsz, NULL);
+	fill_stat(bench, Scatter_LSgrad, Zidxs, Zidxs0, WF, verbose);
+	FREE(Zidxs);
+	// gradZdecomposeR
+	Zidxs = gradZdecomposeR(N, sz, WF->coordinates, grads, &Zsz, NULL);
+	fill_stat(bench, Scatter_Zhao, Zidxs, Zidxs0, WF, verbose);
+	FREE(grads);
+	FREE(P);
 	free_wavefront(&WF);
 	/*WF = calc_surface(G.imagesize, Zidxs0, znum);
 	P = WF->coordinates;

benchmark.h

@@ -26,8 +26,9 @@
 #include "zernike.h"
 #include "usefull_macros.h"
 
+// in lambdas
 #ifndef WF_EPSILON
-#define WF_EPSILON (1e-12)
+#define WF_EPSILON (1e-3)
 #endif
 
 typedef struct{
@@ -51,6 +52,7 @@ typedef enum{
 	,Scatter_QR              // QR_decompose
 	,Scatter_grad            // directGradZdecomposeR
 	,Scatter_LSgrad          // LS_gradZdecomposeR
+	,Scatter_Zhao            // gradZdecomposeR
 	,WF_bench_size
 } WF_bench_type;
 
@@ -68,4 +70,6 @@ wavefront *calc_surfaceRS(int sz, polar *points, double *Zidxs, int znum);
 polar *calc_BTA_Hpoints(int *sz);
 void free_bench(WF_benchmark **bench);
 void do_benchmark(double *Zidxs0, int znum0);
+void WF_diff(double *WF1, double *WF2, int sz, double *std, double *maxd, double *maxdr);
+
 #endif // __BENCHMARK_H__

cmdlnopts.c

@@ -59,12 +59,13 @@ myoption cmdlnopts[] = {
 	{"verbose",	NO_ARGS,	NULL,	'v',	arg_int,	APTR(&verbose),		N_("be more and more verbose")},
 	{"rewrite",	NO_ARGS,	NULL,	'r',	arg_none,	APTR(&rewrite_ifexists),N_("rewrite existant files")},
 	{"gradients",NEED_ARG,	NULL,	'g',	arg_string,	APTR(&G.gradname),	N_("file with pre-computed gradients")},
+	{"wavefront",NEED_ARG,	NULL,	'w',	arg_string,	APTR(&G.wf_fname),	N_("file with mirror aberrations data (X, Y, dZ in meters)")},
 	{"pixsize",	NEED_ARG,	NULL,	'p',	arg_double,	APTR(&G.pixsize),	N_("CCD pixel size (if differs from FITS header)")},
 	{"zern-gen",NO_ARGS,	NULL,	'z',	arg_none,	APTR(&G.zerngen),	N_("generate FITS file with wavefront by given Zernike coefficients, in benchmark mode - use this coefficients instead of random")},
 	{"output",	NEED_ARG,	NULL,	'o',	arg_string,	APTR(&G.outfile),	N_("output file name (\"wavefront\" by default")},
 	{"benchmark",NO_ARGS,	NULL,	'b',	arg_none,	APTR(&G.benchmark),	N_("benchmark: test accuracy of different methods")},
 	{"imsize",	NEED_ARG,	NULL,	's',	arg_int,	APTR(&G.imagesize),	N_("size of output image (rectangular)")},
-	{"bench-maxP",NEED_ARG,	NULL,	'm',	arg_int,	APTR(&G.bench_maxP),N_("maximum amount of Zernike coefficients for benchmark")},
+	{"bench-maxP",NEED_ARG,	NULL,	'm',	arg_int,	APTR(&G.bench_maxP),N_("maximum amount of Zernike coefficients in decomposition")},
 	{"bench-iter",NEED_ARG,	NULL,	'i',	arg_int,	APTR(&G.bench_iter),N_("amount of iterations for benchmark (if -z is absent)")},
 	{"zprec",	NEED_ARG,	NULL,	'Z',	arg_double,	APTR(&Z_prec),		N_("minimal value of non-zero Zernike coefficients")},
 	end_option

cmdlnopts.h

@@ -29,6 +29,7 @@
 typedef struct{
 	char **rest_pars;    // names of input files or zernike coefficients
 	int rest_pars_num;   // amount of files' names / coefficients
+	char *wf_fname;      // filename with wavefront data to restore
 	char *gradname;      // file with pre-computed gradients
 	double pixsize;      // CCD pixel size
 	int zerngen;         // generate FITS file with Zernike surface

main.c

@@ -47,7 +47,33 @@ int main(int argc, char **argv){
 	point *grads = NULL;
 	initial_setup();
 	parse_args(argc, argv);
-	if(G.zerngen){ // user give his zernike coefficients
+	if(G.wf_fname){ // just restore wavefront
+		double Rmax;
+		wavefront *wf = read_wavefront(G.wf_fname, &Rmax);
+		printf("Rmax: %g meters\n", Rmax);
+		int i, N, Zsz, lastidx;
+		convert_Zidx(G.bench_maxP, &N, &Zsz);
+		double *Zidxs = LS_decompose(N, wf->size, wf->coordinates, wf->zdata, &Zsz, &lastidx);
+		printf("LS_decompose. %d non-zero coefficients:\nIDX\tcoefficient\n", lastidx);
+		for(i = 0; i < lastidx; ++i) printf("%3d\t%g\n", i, Zidxs[i]);
+		double *zdata = ZcomposeR(lastidx, Zidxs, wf->size, wf->coordinates);
+		double std, maxd, maxdr;
+		WF_diff(wf->zdata, zdata, wf->size, &std, &maxd, &maxdr);
+		printf("Difference: std=%g, maxD=%g (%g%%)\n", std, maxd, maxdr*100.);
+		FREE(Zidxs);
+		Zidxs = ZdecomposeR(N, wf->size, wf->coordinates, wf->zdata, &Zsz, &lastidx);
+		printf("ZdecomposeR. %d non-zero coefficients:\nIDX\tcoefficient\n", lastidx);
+		for(i = 0; i < lastidx; ++i) printf("%3d\t%g\n", i, Zidxs[i]);
+		FREE(zdata);
+		zdata = ZcomposeR(lastidx, Zidxs, wf->size, wf->coordinates);
+		WF_diff(wf->zdata, zdata, wf->size, &std, &maxd, &maxdr);
+		printf("Difference: std=%g, maxD=%g (%g%%)\n", std, maxd, maxdr*100.);
+		FREE(Zidxs);
+		FREE(zdata);
+		free_wavefront(&wf);
+		return 0;
+	}
+	if(G.zerngen){ // user give Zernike coefficients
 		if(G.rest_pars_num < 1) ERRX(_("You should give at least one Zernike coefficient"));
 		Zidxs = MALLOC(double, G.rest_pars_num);
 		for(i = 0; i < G.rest_pars_num; ++i){

spots.c

@@ -441,6 +441,56 @@ void calc_gradients(mirror *mir, spot_diagram *foc_spots){
 	}
 }
 
+/**
+ * Readout of wavefront file calculated somewhere outside
+ *
+ * @param fname  (i) - input file name
+ * @param scale  (o) - scale on mirror (Rmax)
+ * @return - wavefront read
+ *
+ * Wavefront measured in lambdas (WAVELEN) and coordinates are normalized by `scale` (Rmax)
+ * 1 lambda on wavefront is lambda/2 on mirror, so we need to multiply data by 2
+ */
+wavefront *read_wavefront(char *fname, double *scale){
+	if(!fname) return NULL;
+	wavefront *F = MALLOC(wavefront, 1);
+	mmapbuf *M = NULL;
+	double Rmax = 0.;
+	int sz, max_sz = 512;
+	polar *coordinates = MALLOC(polar, max_sz);
+	double *zdata = MALLOC(double, max_sz);
+	M = My_mmap(fname);
+	char *pos = M->data, *epos = pos + M->len;
+	for(sz = -1; pos && pos < epos; pos = strchr(pos+1, '\n')){
+		double x, y, z, R;
+		if(3 != sscanf(pos, "%lf %lf %lf", &x, &y, &z))
+			continue;
+		if(++sz >= max_sz){
+			max_sz += 512;
+			double *dptr = realloc(F->zdata, max_sz * sizeof(double));
+			if(dptr) zdata = dptr;
+			else ERR("realloc");
+			polar *ptr = realloc(F->coordinates, max_sz * sizeof(polar));
+			if(ptr) coordinates = ptr;
+			else ERR("realloc");
+		}
+		R = sqrt(x*x + y*y);
+		if(R > Rmax) Rmax = R;
+		coordinates[sz].r = R;
+		coordinates[sz].theta = atan2(y, x);
+		//DBG("%3d. x:%g, y:%g, z:%g, R:%g, theta: %g\n", sz, x,y,z, R, coordinates[sz].theta);
+		zdata[sz] = 2.*z/WAVELEN;
+	}
+	DBG("Found %d points", sz+1);
+	F->size = sz + 1;
+	F->zdata = zdata;
+	F->coordinates = coordinates;
+	for(; sz > -1; --sz) F->coordinates[sz].r /= Rmax;
+	if(scale) *scale = Rmax;
+	My_munmap(M);
+	return F;
+}
+
 #if 0
 /**
  *

spots.h

@@ -26,6 +26,7 @@
 #include <stdint.h>
 #include "zernike.h"
 #include "simple_list.h"
+#include "benchmark.h"
 
 extern double pixsize;
 extern double distance;
@@ -91,6 +92,8 @@ double calc_Hartmann_constant(mirror *mir);
 spot_diagram *calc_spot_diagram(mirror *mir, double z);
 void calc_gradients(mirror *mir, spot_diagram *foc_spots);
 
+wavefront *read_wavefront(char *fname, double *scale);
+
 /*
 size_t get_gradients(hartmann *H[], polar **coords, point **grads, double *scale);
 size_t read_gradients(char *gradname, polar **coords, point **grads, double *scale);

zernike.h

@@ -32,6 +32,8 @@
 #define MIR_R       (3025.)
 // Distance from mirror to hartmann mask
 #define HARTMANN_Z  (20017.)
+// reference wavelength (meters): 650nm
+#define WAVELEN     (6.5e-7)
 
 /*************** Data structures & typedefs ***************/
 // point coordinates

zernikeR.c

@@ -388,7 +388,6 @@ double *QR_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int
 	return Zidxs_corr;
 }
 
-
 /**
  * Components of Zj gradient without constant factor
  * @param n,m      - orders of polynomial
@@ -491,7 +490,7 @@ point *zerngradR(int p, int Sz, polar *P, double *norm){
 }
 
 /**
- * Decomposition of image with normals to wavefront by Zhao's polynomials
+ * Decomposition of image with normals to wavefront by Zernike polynomials
  * by least squares method through LU-decomposition
  * @param Nmax (i)   - maximum power of Zernike polinomial for decomposition
  * @param Sz, P(i)   - size (Sz) of points array (P)