diff --git a/benchmark.c b/benchmark.c index 89f38dc..c13f4f6 100644 --- a/benchmark.c +++ b/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; diff --git a/benchmark.h b/benchmark.h index b27b09f..39ed617 100644 --- a/benchmark.h +++ b/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__ diff --git a/cmdlnopts.c b/cmdlnopts.c index 68f6c02..ae0dbf5 100644 --- a/cmdlnopts.c +++ b/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 diff --git a/cmdlnopts.h b/cmdlnopts.h index 3488c43..5426ced 100644 --- a/cmdlnopts.h +++ b/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 diff --git a/main.c b/main.c index 2193a33..80f563c 100644 --- a/main.c +++ b/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){ diff --git a/spots.c b/spots.c index 3393501..4980a77 100644 --- a/spots.c +++ b/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 /** * diff --git a/spots.h b/spots.h index 522a19d..03c5c53 100644 --- a/spots.h +++ b/spots.h @@ -26,6 +26,7 @@ #include #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); diff --git a/zernike.h b/zernike.h index 3463061..f182f87 100644 --- a/zernike.h +++ b/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 diff --git a/zernikeR.c b/zernikeR.c index 12adfcd..6c423d5 100644 --- a/zernikeR.c +++ b/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)