From abee079d267836c02b3cc5522ced896618ec70d2 Mon Sep 17 00:00:00 2001 From: eddyem Date: Thu, 14 Jul 2016 09:47:49 +0300 Subject: [PATCH] copy --- Makefile | 24 ++ Readme.md | 9 + benchmark.c | 307 +++++++++++++++++++++ benchmark.h | 71 +++++ cmdlnopts.c | 97 +++++++ cmdlnopts.h | 49 ++++ main.c | 201 ++++++++++++++ main.h | 0 parseargs.c | 497 +++++++++++++++++++++++++++++++++ parseargs.h | 124 +++++++++ saveimg.c | 294 ++++++++++++++++++++ saveimg.h | 33 +++ simple_list.c | 93 +++++++ simple_list.h | 39 +++ spots.c | 556 +++++++++++++++++++++++++++++++++++++ spots.h | 98 +++++++ usefull_macros.c | 360 ++++++++++++++++++++++++ usefull_macros.h | 134 +++++++++ zern_private.h | 46 ++++ zernike.c | 494 +++++++++++++++++++++++++++++++++ zernike.h | 92 +++++++ zernikeR.c | 687 ++++++++++++++++++++++++++++++++++++++++++++++ zernikeR.h | 0 zernike_annular.c | 388 ++++++++++++++++++++++++++ zernike_annular.h | 0 25 files changed, 4693 insertions(+) create mode 100644 Makefile create mode 100644 Readme.md create mode 100644 benchmark.c create mode 100644 benchmark.h create mode 100644 cmdlnopts.c create mode 100644 cmdlnopts.h create mode 100644 main.c create mode 100644 main.h create mode 100644 parseargs.c create mode 100644 parseargs.h create mode 100644 saveimg.c create mode 100644 saveimg.h create mode 100644 simple_list.c create mode 100644 simple_list.h create mode 100644 spots.c create mode 100644 spots.h create mode 100644 usefull_macros.c create mode 100644 usefull_macros.h create mode 100644 zern_private.h create mode 100644 zernike.c create mode 100644 zernike.h create mode 100644 zernikeR.c create mode 100644 zernikeR.h create mode 100644 zernike_annular.c create mode 100644 zernike_annular.h diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..60ef2d5 --- /dev/null +++ b/Makefile @@ -0,0 +1,24 @@ +# run `make DEF=...` to add extra defines +OUT = zerndeco +LDFLAGS = -lm -lgsl -lgslcblas -lpng -lcfitsio -fopenmp +SRCS = $(wildcard *.c) +CC = gcc +DEFINES = $(DEF) -D_GNU_SOURCE -DEBUG +CFLAGS = -std=gnu99 -fopenmp -Wall -Werror -Wextra $(DEFINES) +OBJS = $(SRCS:.c=.o) + +all : $(OUT) + +$(OUT) : $(OBJS) + $(CC) $(OBJS) $(LDFLAGS) -o $(OUT) + +$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS) + touch $@ + +%.h: ; + +clean: + /bin/rm -f *.o *~ + +depend: + $(CC) -MM $(SRCS) diff --git a/Readme.md b/Readme.md new file mode 100644 index 0000000..093ed33 --- /dev/null +++ b/Readme.md @@ -0,0 +1,9 @@ +Wavefront decomposition by Zernike polynomials +============================================== + + +Coordinate system - Cartesian (X from left to right, Y from bottom to top) +On images coordinate system the same, so first byte of image is its left bottom corner! + +make: +DEF=DWF_EPSILON=eps - set minimal value of vawefront that would mean zero to `eps` (default: 1e-12) diff --git a/benchmark.c b/benchmark.c new file mode 100644 index 0000000..89f38dc --- /dev/null +++ b/benchmark.c @@ -0,0 +1,307 @@ +/* + * benchmark.c - test of different methods + * + * Copyright 2016 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include +#include // for srand48(time(NULL)) +#include "benchmark.h" +#include "usefull_macros.h" +#include "cmdlnopts.h" + +void free_wavefront(wavefront **F){ + FREE((*F)->zdata); + FREE((*F)->coordinates); + FREE(*F); +} +static void free_stat(WF_stat *st){ + if(!st) return; + FREE(st->ZC_sum); + FREE(st->ZC_sum2); + FREE(st->max_dZC); +} + +void free_bench(WF_benchmark **bench){ + if(!bench || !*bench) return; + WF_benchmark *b = *bench; + for(int i = 0; i < WF_bench_size; ++i) + free_stat(&(b->stat[i])); + FREE(*bench); +} + +/** + * calculate Z-coordinate of surface with rectangle nodes by given Zernike coefficients + * @param sz - image size + * @param Zidxs - Zernike coefficients + * @param znum - amount of Zidx + */ +double *calc_surface(int sz, double *Zidxs, int znum){ + if(!Zidxs || znum < 1) return NULL; + return Zcompose(znum, Zidxs, sz, sz); +} + +/** + * calculate surface with rectangle nodes by given Zernike coefficients using `R`-functions + * @param sz - image size + * @param Zidxs - Zernike coefficients + * @param znum - amount of Zidx + */ +wavefront *calc_surfaceR(int sz, double *Zidxs, int znum){ + int G2 = sz * sz; + if(!Zidxs || znum < 1) return NULL; + wavefront *F = MALLOC(wavefront, 1); + F->size = G2; + F->coordinates = MALLOC(polar, G2); + polar *rptr = F->coordinates; + double Stp = 2./((double)sz - 1.); + for(int j = 0; j < sz; ++j){ + double y = ((double) j) * Stp - 1.; + for(int i = 0; i < sz; ++i, ++rptr){ + double x = ((double) i)* Stp - 1.; + double R2 = x*x + y*y; + rptr->r = sqrt(R2); + rptr->theta = atan2(y, x); + } + } + F->zdata = ZcomposeR(znum, Zidxs, G2, F->coordinates); + return F; +} + +/** + * calculate surface with scattered nodes by given Zernike coefficients using `R`-functions + * @param sz - amount of points + * @param points - given nodes coordinates + * @param Zidxs - Zernike coefficients + * @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; + wavefront *F = MALLOC(wavefront, 1); + F->size = sz; + F->coordinates = MALLOC(polar, sz); + memcpy(F->coordinates, points, sizeof(polar)*sz); + F->zdata = ZcomposeR(znum, Zidxs, sz, F->coordinates); + return F; +} + +/** + * calculate coordinates of wavefront for BTA Hartmann mask + */ +polar *calc_BTA_Hpoints(int *sz){ + // radius on Hartmann mask (in mm) - for normalize coordinates @ [-1, 1] + double rmax = MIR_R * (1. - HARTMANN_Z/FOCAL_R); + double r[] = {175., 247., 295., 340., 379., 414., 448., 478.}; // radii of rings + polar *pts = MALLOC(polar, 258), *pt = pts; + const double ray_step = M_PI/16., theta0 = M_PI/32; // angle between rays & first ray theta coordinate + for(int ring = 0; ring < 8; ++ring){ // main rings + for(int ray = 0; ray < 32; ++ray, ++pt){ + pt->r = r[ring] / rmax; + pt->theta = ray * ray_step + theta0; + } + } + // and markers: ring 2 ray 24.5 and ring 7 ray 29.5 + 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; + return pts; +} + +static char *bench_names[] = { + N_("Direct decomposition (ZdecomposeR)"), + N_("Less square decomposition (LS_decompose)"), + N_("QR-decomposition (QR_decompose)"), + N_("Direct gradient decomposition (directGradZdecomposeR)"), + N_("LS gradient decomposition (LS_gradZdecomposeR)") +}; + +/** + * calculate errors and fill WF_stat fields for single measurement + * @param stat - output stat structure + * @param Zidxs - calculated ZC + * @param Zidxs0 - original ZC + * @param znum - size of 'Zidxs' + * @param WF0 - original wavefront + * @param verbose - ==1 for messages in stdout + */ +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; + 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); + 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); + } + for(int p = 0; p < znum; ++p){ + double diff = Zidxs0[p] - Zidxs[p]; + stat->ZC_sum[p] = diff; + stat->ZC_sum2[p] = diff*diff; + if(fabs(Zidxs0[p]) > Z_prec) stat->max_dZC[p] = fabs(diff / Zidxs0[p]); + if(verbose){ + printf("%d\t%10f\t%10f\t", p, Zidxs0[p], Zidxs[p]); + if(fabs(Zidxs0[p]) > Z_prec) printf("%4.1f\n", stat->max_dZC[p]*100.); + else printf(" - \n"); + } + } +} + +/** + * Process a benchmark step with different methods of wavefront restoration + * @param Zidxs - Zernike coefficients + * @param znum - length of Zidxs + * @param verbose - ==1 for detailed information on stdout + */ +static WF_benchmark *benchmark_step(double *Zidxs0, int znum, int verbose){ + #define MESG(...) do{if(verbose){printf(__VA_ARGS__); printf("\n");}}while(0) + int sz; //G2 = G.imagesize * G.imagesize; + wavefront *WF; + polar *P; + double *Zidxs, *zdata; + int N, M, Zsz; + convert_Zidx(znum, &N, &M); // calculate max power of Zc + WF_benchmark *bench = MALLOC(WF_benchmark, 1); + bench->Znum = znum; + MESG(_("Direct decomposition of given wavefront")); + 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); + fill_stat(bench, Scatter_direct, Zidxs, Zidxs0, WF, verbose); + FREE(Zidxs); + // LS_decompose + Zidxs = LS_decompose(N, sz, WF->coordinates, zdata, &Zsz, NULL); + fill_stat(bench, Scatter_LS, Zidxs, Zidxs0, WF, verbose); + FREE(Zidxs); + // QR_decompose + Zidxs = QR_decompose(N, sz, WF->coordinates, zdata, &Zsz, NULL); + fill_stat(bench, Scatter_QR, Zidxs, Zidxs0, WF, verbose); + FREE(Zidxs); + // directGradZdecomposeR + Zidxs0[0] = 0.; // decomposition by gradients cannot know anything about zero coefficient + 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); + free_wavefront(&WF); + /*WF = calc_surface(G.imagesize, Zidxs0, znum); + P = WF->coordinates; + zdata = WF->zdata; + grads = gradZcomposeR(znum, Zidxs0, G2, P); + Zidxs = ZdecomposeR(N, G2, WF->coordinates;, zdata, &Zsz, &lastIdx); + FREE(grads); + */ + #undef MESG + return bench; +} + +/** + * Update values stored in `total` by values of `current` + */ +static void update_stat_values(WF_benchmark *total, WF_benchmark *current){ + ++total->measurements; + int znum = total->Znum; + for(int i = 0; i < WF_bench_size; ++i){ + WF_stat *out = &total->stat[i], *in = ¤t->stat[i]; + if(!out->ZC_sum || !in->ZC_sum) continue; // given field is empty - didn't run benchmark + out->WF_std += in->WF_std; + if(in->max_dWF > out->max_dWF) + out->max_dWF = in->max_dWF; + for(int j = 0; j < znum; ++j){ + out->ZC_sum[j] += in->ZC_sum[j]; + out->ZC_sum2[j] += in->ZC_sum2[j]; + if(out->max_dZC[i] < in->max_dZC[i]) + out->max_dZC[i] = in->max_dZC[i]; + } + } +} + +/** + * Process a benchmark with different methods of wavefront restoration + * @param Zidxs - user defined Zernike coefficients or NULL + * @param znum - length of Zidxs + */ +void do_benchmark(double *Zidxs0, int znum0){ + int znum = G.bench_maxP; + double *Zidxs = MALLOC(double, znum); + srand48(time(NULL)); + inline void gen_zidxs(){ // generate random Z coefficients in interval (-10..10)/p + for(int i = 0; i < znum; ++i) + Zidxs[i] = (20.*drand48() - 10.) / ((double)i+1); + } + WF_benchmark *total = NULL, *current = NULL; + if(Zidxs0){ // do only one iteration - for given coefficients + current = benchmark_step(Zidxs0, znum0, 1); + free_bench(¤t); + return; + } + if(G.bench_iter < 2) ERRX(_("Need at least two iterations for benchmark")); + gen_zidxs(); + total = benchmark_step(Zidxs, znum, 0); + total->measurements = 1; + for(int i = 1; i < G.bench_iter; ++i){ + gen_zidxs(); + current = benchmark_step(Zidxs, znum, 0); + update_stat_values(total, current); + free_bench(¤t); + } + // print stat values + double Nm = total->measurements; + for(int i = 0; i < WF_bench_size; ++i){ + WF_stat *stat = &total->stat[i]; + if(!stat->ZC_sum) continue; + printf("%s:\n = %g\n", bench_names[i], stat->WF_std/Nm); + printf("max_dWF=%g%%\n", stat->max_dWF); + printf(" p(n,m)\t Z_std \tmax_dZC(%%) \tsum2\tsum\n"); + for(int i = 0; i < znum; ++i){ + double sum2 = stat->ZC_sum2[i] / Nm; + double sum = stat->ZC_sum[i] / Nm; + int N,M; + convert_Zidx(i, &N, &M); + printf("%3d (%2d, %2d)\t%10g\t%10g\t%10g\t%10g\n", i, N, M, sqrt(sum2 - sum*sum), + stat->max_dZC[i]*100., sum2, sum); + } + } + free_bench(&total); + FREE(Zidxs); +} + diff --git a/benchmark.h b/benchmark.h new file mode 100644 index 0000000..b27b09f --- /dev/null +++ b/benchmark.h @@ -0,0 +1,71 @@ +/* + * benchmark.h + * + * Copyright 2016 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __BENCHMARK_H__ +#define __BENCHMARK_H__ + +#include "zernike.h" +#include "usefull_macros.h" + +#ifndef WF_EPSILON +#define WF_EPSILON (1e-12) +#endif + +typedef struct{ + size_t size; // size of arrays + double *zdata; // Z-coordinate of wavefront + polar *coordinates; // coordinates on WF (-1..1) +} wavefront; + +// errors of wavefront restoration method +typedef struct{ + double WF_std; // standard deviation for restored wavefront weighted by wavefront values + double max_dWF; // max{|WF - WF_0|/WF_0} if WF_0 > WF_EPSILON + double *ZC_sum; // data for std calculation: sum of differences + double *ZC_sum2; // and sum of their squares + double *max_dZC; // max|ZC/ZC_0 - 1| +} WF_stat; + +typedef enum{ + Scatter_direct = 0 // ZdecomposeR + ,Scatter_LS // LS_decompose + ,Scatter_QR // QR_decompose + ,Scatter_grad // directGradZdecomposeR + ,Scatter_LSgrad // LS_gradZdecomposeR + ,WF_bench_size +} WF_bench_type; + +// wavefront restoration methods (according to functions' names) +typedef struct{ + WF_stat stat[WF_bench_size]; + int Znum; // amount of Zernike coefficients + int measurements; // amount of measurements: std=sqrt(ZC_sum2/measurements-(ZC_sum/measurements)^2) +} WF_benchmark; + +void free_wavefront(wavefront **F); +double *calc_surface(int sz, double *Zidxs, int znum); +wavefront *calc_surfaceR(int sz, double *Zidxs, int znum); +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); +#endif // __BENCHMARK_H__ diff --git a/cmdlnopts.c b/cmdlnopts.c new file mode 100644 index 0000000..68f6c02 --- /dev/null +++ b/cmdlnopts.c @@ -0,0 +1,97 @@ +/* + * cmdlnopts.c - the only function that parse cmdln args and returns glob parameters + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include +#include +#include +#include +#include "zernike.h" +#include "cmdlnopts.h" +#include "usefull_macros.h" + +#define RAD 57.2957795130823 +#define D2R(x) ((x) / RAD) +#define R2D(x) ((x) * RAD) + +/* + * here are global parameters initialisation + */ +int help; +glob_pars G; + +int rewrite_ifexists = 0, // rewrite existing files == 0 or 1 + verbose = 0; // each -v increments this value, e.g. -vvv sets it to 3 +// DEFAULTS +// default global parameters +glob_pars const Gdefault = { + .pixsize = 30.e-3 + ,.outfile = "wavefront" + ,.imagesize = 1024 + ,.bench_maxP = 21 // Nmax = 5 + ,.bench_iter = 100 +}; + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +myoption cmdlnopts[] = { + // set 1 to param despite of its repeating number: + {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), N_("show this help")}, + {"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")}, + {"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-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 +}; + +/** + * Parse command line options and return dynamically allocated structure + * to global parameters + * @param argc - copy of argc from main + * @param argv - copy of argv from main + * @return allocated structure with global parameters + */ +glob_pars *parse_args(int argc, char **argv){ + int i; + void *ptr; + ptr = memcpy(&G, &Gdefault, sizeof(G)); assert(ptr); + // format of help: "Usage: progname [args]\n" + change_helpstring(_("Usage: %s [args] \n\n\tWhere args are:\n")); + // parse arguments + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(argc > 0){ + G.rest_pars_num = argc; + G.rest_pars = calloc(argc, sizeof(char*)); + for (i = 0; i < argc; i++) + G.rest_pars[i] = strdup(argv[i]); + } + return &G; +} + diff --git a/cmdlnopts.h b/cmdlnopts.h new file mode 100644 index 0000000..3488c43 --- /dev/null +++ b/cmdlnopts.h @@ -0,0 +1,49 @@ +/* + * cmdlnopts.h - comand line options for parceargs + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __CMDLNOPTS_H__ +#define __CMDLNOPTS_H__ + +#include "parseargs.h" + + +typedef struct{ + char **rest_pars; // names of input files or zernike coefficients + int rest_pars_num; // amount of files' names / coefficients + char *gradname; // file with pre-computed gradients + double pixsize; // CCD pixel size + int zerngen; // generate FITS file with Zernike surface + char *outfile; // output file name + int imagesize; // output image size (1024 by default) + int benchmark; // different methods benchmark + int bench_maxP; // max ZC number for behchmark (default - 20) + int bench_iter; // amount of iterations for benchmark (default - 100) +} glob_pars; + + +// default & global parameters +extern glob_pars const Gdefault; +extern int rewrite_ifexists, verbose; + +glob_pars *parse_args(int argc, char **argv); +extern glob_pars G; +#endif // __CMDLNOPTS_H__ diff --git a/main.c b/main.c new file mode 100644 index 0000000..2193a33 --- /dev/null +++ b/main.c @@ -0,0 +1,201 @@ +/* + * Z-BTA_test.c - simple test for models of hartmannograms + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include +#include +#include +#include +#include "usefull_macros.h" +#include "cmdlnopts.h" +#include "zernike.h" +#include "spots.h" +#include "saveimg.h" +#include "benchmark.h" + + +void signals(int signo){ + exit(signo); +} + +int main(int argc, char **argv){ + int i, j; + //double scale; + hartmann **images; + mirror *mir = NULL; + wavefront *comp; + double *Zidxs = NULL; + size_t L; + polar *coords = NULL; + point *grads = NULL; + initial_setup(); + parse_args(argc, argv); + if(G.zerngen){ // user give his 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){ + if(!str2double(&Zidxs[i], G.rest_pars[i])){ + ERRX(_("Bad double number: %s!"), G.rest_pars[i]); + } + } + } + if(G.benchmark){ + DBG("Benchmark"); + do_benchmark(Zidxs, G.rest_pars_num); + FREE(Zidxs); + return 0; + } + if(G.zerngen){ // generate Zernike wavefront + DBG("Generate Zernike surface"); + double *z = calc_surface(G.imagesize, Zidxs, G.rest_pars_num); + if(z){ + writeimg(G.outfile, G.imagesize, mir, z); + FREE(z); + } + return 0; + } + if(!G.gradname){ + DBG("Calculate wavefront by hartmannograms"); + if(G.rest_pars_num < 2) ERR(_("You should give at least two file names: pre- and postfocal!")); + images = MALLOC(hartmann*, G.rest_pars_num + 1); + for(i = 0; i < G.rest_pars_num; ++i) + images[i] = read_spots(G.rest_pars[i]); + mir = calc_mir_coordinates(images); + getQ(mir); + calc_Hartmann_constant(mir); + spot_diagram *spot_dia = calc_spot_diagram(mir, mir->z07); + //printf("\nmirror's coordinates should be corrected to (%g, %g)\n", + //spot_dia->center.x, spot_dia->center.y); + printf("Projection of center on mirror shifted by (%g, %g), image shifted by (%g, %g)\n", + images[1]->center.x*HARTMANN_Z/distance, images[1]->center.y*HARTMANN_Z/distance, + images[1]->center.x*(FOCAL_R-HARTMANN_Z)/distance, images[1]->center.y*(FOCAL_R-HARTMANN_Z)/distance); + double tr = sqrt(images[1]->center.x*images[1]->center.x+images[1]->center.y*images[1]->center.y); + printf("Beam tilt is %g''\n", tr/distance*206265.); + calc_gradients(mir, spot_dia); + coords = MALLOC(polar, mir->spotsnum); + grads = MALLOC(point, mir->spotsnum); + printf("\nGradients of aberrations (*1e-6):\nray# x y dx dy R phi\n"); + for(j = 0, i = 0; i < 258; ++i){ + if(!mir->got[i]) continue; + printf("%4d %8.2f %8.2f %10.6f %10.6f %10.2f %10.6f\n", + i, mir->spots[i].x, mir->spots[i].y, + mir->grads[i].x * 1e6, mir->grads[i].y * 1e6, + mir->pol_spots[i].r * MIR_R, mir->pol_spots[i].theta * 180. / M_PI); + memcpy(&coords[j], &mir->pol_spots[i], sizeof(polar)); + memcpy(&grads[j], &mir->grads[i], sizeof(point)); + grads[j].x *= 1e3; + grads[j].y *= 1e3; + j++; + } + L = mir->spotsnum; + //scale = mir->Rmax; + FREE(spot_dia); + FREE(mir); + for(i = 0; i < G.rest_pars_num; ++i) + h_free(&images[i]); + FREE(images); + }else{ +// L = read_gradients(gradname, &coords, &grads, &scale); + return 0; + } +/* + // spots information + printf(GREEN "\nSpots:\n" OLDCOLOR "\n r\ttheta\tDx(mm/m)\tDy(mm/m)\n"); + for(i = 0; i < L; i++){ + printf("%8.1f%8.4f%8.4f%8.4f\n", coords[i].r, coords[i].theta, + grads[i].x*1000., grads[i].y*1000.); + } + */ + // gradients decomposition (Less squares, Zhao polinomials) + int Zsz, lastidx; + Zidxs = gradZdecomposeR(10, L, coords, grads, &Zsz, &lastidx); + //double *Zidxs = LS_gradZdecomposeR(10, L, coords, grads, &Zsz, &lastidx); +// Zidxs[1] = 0.; +// Zidxs[2] = 0.; + lastidx++; + printf("\nGradZ decompose, coefficients (%d):\n", lastidx); + for(i = 0; i < lastidx; i++) printf("%g, ", Zidxs[i]); + printf("\n\n"); + + int GridSize = 15; + comp = calc_surfaceR(GridSize, Zidxs, lastidx); + if(comp){ + double zero_val = 0., N = 0., *zd = comp->zdata; + polar *rect = comp->coordinates; + for(j = GridSize-1; j > -1; j--){ + int idx = j*GridSize; + for(i = 0; i < GridSize; i++,idx++){ + if(zd[idx] > 1e-9){ + zero_val += zd[idx]; + N++; + } + } + } + zero_val /= N; + printf("zero: %g\nMirror surface:\n", zero_val); + for(j = GridSize-1; j > -1; j--){ + int idx = j*GridSize; + for(i = 0; i < GridSize; i++,idx++){ + if(rect[idx].r > 1. || rect[idx].r < 0.2){ + printf(" "); + }else{ + printf("%7.3f", zd[idx] * MIR_R); + } + } + printf("\n"); + } + free_wavefront(&comp); + } + // save matrix 100x100 with mirror deformations in Octave file format + FILE *f = fopen("file.out", "w"); + if(f){ + GridSize = 100; + comp = calc_surfaceR(GridSize, Zidxs, lastidx); + if(comp){ + polar *rect = comp->coordinates; + double *zd = comp->zdata; + fprintf(f, "# generated by Z-BTA_test\n" + "# name: dev_matrix\n# type: matrix\n# rows: %d\n# columns: %d\n", + GridSize, GridSize); + for(j = 0; j < GridSize; j++){ + int idx = j*GridSize; + for(i = 0; i < GridSize; i++,idx++){ + if(rect[idx].r > 1. || rect[idx].r < 0.2) + fprintf(f, "0. "); + else + fprintf(f, "%g ", zd[idx] * MIR_R); + } + fprintf(f,"\n"); + } + free_wavefront(&comp); + } + fclose(f); + } + // save image + double *z = calc_surface(G.imagesize, Zidxs, lastidx); + if(z){ + writeimg(G.outfile, G.imagesize, mir, z); + FREE(z); + } + FREE(Zidxs); + FREE(coords); + FREE(grads); + return 0; +} diff --git a/main.h b/main.h new file mode 100644 index 0000000..e69de29 diff --git a/parseargs.c b/parseargs.c new file mode 100644 index 0000000..64bfab1 --- /dev/null +++ b/parseargs.c @@ -0,0 +1,497 @@ +/* + * parseargs.c - parsing command line arguments & print help + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include // printf +#include // getopt_long +#include // calloc, exit, strtoll +#include // assert +#include // strdup, strchr, strlen +#include // strcasecmp +#include // INT_MAX & so on +#include // gettext +#include // isalpha +#include "parseargs.h" +#include "usefull_macros.h" + +char *helpstring = "%s\n"; + +/** + * Change standard help header + * MAY consist ONE "%s" for progname + * @param str (i) - new format + */ +void change_helpstring(char *s){ + int pcount = 0, scount = 0; + char *str = s; + // check `helpstring` and set it to default in case of error + for(; pcount < 2; str += 2){ + if(!(str = strchr(str, '%'))) break; + if(str[1] != '%') pcount++; // increment '%' counter if it isn't "%%" + else{ + str += 2; // pass next '%' + continue; + } + if(str[1] == 's') scount++; // increment "%s" counter + }; + if(pcount > 1 || pcount != scount){ // amount of pcount and/or scount wrong + /// "Неправильный формат строки помощи" + ERRX(_("Wrong helpstring!")); + } + helpstring = s; +} + +/** + * Carefull atoll/atoi + * @param num (o) - returning value (or NULL if you wish only check number) - allocated by user + * @param str (i) - string with number must not be NULL + * @param t (i) - T_INT for integer or T_LLONG for long long (if argtype would be wided, may add more) + * @return TRUE if conversion sone without errors, FALSE otherwise + */ +static bool myatoll(void *num, char *str, argtype t){ + long long tmp, *llptr; + int *iptr; + char *endptr; + assert(str); + assert(num); + tmp = strtoll(str, &endptr, 0); + if(endptr == str || *str == '\0' || *endptr != '\0') + return FALSE; + switch(t){ + case arg_longlong: + llptr = (long long*) num; + *llptr = tmp; + break; + case arg_int: + default: + if(tmp < INT_MIN || tmp > INT_MAX){ + /// "Целое вне допустимого диапазона" + WARNX(_("Integer out of range")); + return FALSE; + } + iptr = (int*)num; + *iptr = (int)tmp; + } + return TRUE; +} + +// the same as myatoll but for double +// There's no NAN & INF checking here (what if they would be needed?) +static bool myatod(void *num, const char *str, argtype t){ + double tmp, *dptr; + float *fptr; + char *endptr; + assert(str); + tmp = strtod(str, &endptr); + if(endptr == str || *str == '\0' || *endptr != '\0') + return FALSE; + switch(t){ + case arg_double: + dptr = (double *) num; + *dptr = tmp; + break; + case arg_float: + default: + fptr = (float *) num; + *fptr = (float)tmp; + break; + } + return TRUE; +} + +/** + * Get index of current option in array options + * @param opt (i) - returning val of getopt_long + * @param options (i) - array of options + * @return index in array + */ +static int get_optind(int opt, myoption *options){ + int oind; + myoption *opts = options; + assert(opts); + for(oind = 0; opts->name && opts->val != opt; oind++, opts++); + if(!opts->name || opts->val != opt) // no such parameter + showhelp(-1, options); + return oind; +} + +/** + * reallocate new value in array of multiple repeating arguments + * @arg paptr - address of pointer to array (**void) + * @arg type - its type (for realloc) + * @return pointer to new (next) value + */ +void *get_aptr(void *paptr, argtype type){ + int i = 1; + void **aptr = *((void***)paptr); + if(aptr){ // there's something in array + void **p = aptr; + while(*p++) ++i; + } + size_t sz = 0; + switch(type){ + default: + case arg_none: + /// "Не могу использовать несколько параметров без аргументов!" + ERRX("Can't use multiple args with arg_none!"); + break; + case arg_int: + sz = sizeof(int); + break; + case arg_longlong: + sz = sizeof(long long); + break; + case arg_double: + sz = sizeof(double); + break; + case arg_float: + sz = sizeof(float); + break; + case arg_string: + sz = 0; + break; + /* case arg_function: + sz = sizeof(argfn *); + break;*/ + } + aptr = realloc(aptr, (i + 1) * sizeof(void*)); + *((void***)paptr) = aptr; + aptr[i] = NULL; + if(sz){ + aptr[i - 1] = malloc(sz); + }else + aptr[i - 1] = &aptr[i - 1]; + return aptr[i - 1]; +} + + +/** + * Parse command line arguments + * ! If arg is string, then value will be strdup'ed! + * + * @param argc (io) - address of argc of main(), return value of argc stay after `getopt` + * @param argv (io) - address of argv of main(), return pointer to argv stay after `getopt` + * BE CAREFUL! if you wanna use full argc & argv, save their original values before + * calling this function + * @param options (i) - array of `myoption` for arguments parcing + * + * @exit: in case of error this function show help & make `exit(-1)` + */ +void parseargs(int *argc, char ***argv, myoption *options){ + char *short_options, *soptr; + struct option *long_options, *loptr; + size_t optsize, i; + myoption *opts = options; + // check whether there is at least one options + assert(opts); + assert(opts[0].name); + // first we count how much values are in opts + for(optsize = 0; opts->name; optsize++, opts++); + // now we can allocate memory + short_options = calloc(optsize * 3 + 1, 1); // multiply by three for '::' in case of args in opts + long_options = calloc(optsize + 1, sizeof(struct option)); + opts = options; loptr = long_options; soptr = short_options; + // in debug mode check the parameters are not repeated +#ifdef EBUG + char **longlist = MALLOC(char*, optsize); + char *shortlist = MALLOC(char, optsize); +#endif + // fill short/long parameters and make a simple checking + for(i = 0; i < optsize; i++, loptr++, opts++){ + // check + assert(opts->name); // check name +#ifdef EBUG + longlist[i] = strdup(opts->name); +#endif + if(opts->has_arg){ + assert(opts->type != arg_none); // check error with arg type + assert(opts->argptr); // check pointer + } + if(opts->type != arg_none) // if there is a flag without arg, check its pointer + assert(opts->argptr); + // fill long_options + // don't do memcmp: what if there would be different alignment? + loptr->name = opts->name; + loptr->has_arg = (opts->has_arg < MULT_PAR) ? opts->has_arg : 1; + loptr->flag = opts->flag; + loptr->val = opts->val; + // fill short options if they are: + if(!opts->flag && opts->val){ +#ifdef EBUG + shortlist[i] = (char) opts->val; +#endif + *soptr++ = opts->val; + if(loptr->has_arg) // add ':' if option has required argument + *soptr++ = ':'; + if(loptr->has_arg == 2) // add '::' if option has optional argument + *soptr++ = ':'; + } + } + // sort all lists & check for repeating +#ifdef EBUG + int cmpstringp(const void *p1, const void *p2){ + return strcmp(* (char * const *) p1, * (char * const *) p2); + } + int cmpcharp(const void *p1, const void *p2){ + return (int)(*(char * const)p1 - *(char *const)p2); + } + qsort(longlist, optsize, sizeof(char *), cmpstringp); + qsort(shortlist,optsize, sizeof(char), cmpcharp); + char *prevl = longlist[0], prevshrt = shortlist[0]; + for(i = 1; i < optsize; ++i){ + if(longlist[i]){ + if(prevl){ + if(strcmp(prevl, longlist[i]) == 0) ERRX("double long arguments: --%s", prevl); + } + prevl = longlist[i]; + } + if(shortlist[i]){ + if(prevshrt){ + if(prevshrt == shortlist[i]) ERRX("double short arguments: -%c", prevshrt); + } + prevshrt = shortlist[i]; + } + } +#endif + // now we have both long_options & short_options and can parse `getopt_long` + while(1){ + int opt; + int oindex = 0, optind = 0; // oindex - number of option in argv, optind - number in options[] + if((opt = getopt_long(*argc, *argv, short_options, long_options, &oindex)) == -1) break; + if(opt == '?'){ + opt = optopt; + optind = get_optind(opt, options); + if(options[optind].has_arg == NEED_ARG || options[optind].has_arg == MULT_PAR) + showhelp(optind, options); // need argument + } + else{ + if(opt == 0 || oindex > 0) optind = oindex; + else optind = get_optind(opt, options); + } + opts = &options[optind]; + if(opt == 0 && opts->has_arg == NO_ARGS) continue; // only long option changing integer flag + // now check option + if(opts->has_arg == NEED_ARG || opts->has_arg == MULT_PAR) + if(!optarg) showhelp(optind, options); // need argument + void *aptr; + if(opts->has_arg == MULT_PAR){ + aptr = get_aptr(opts->argptr, opts->type); + }else + aptr = opts->argptr; + bool result = TRUE; + // even if there is no argument, but argptr != NULL, think that optarg = "1" + if(!optarg) optarg = "1"; + switch(opts->type){ + default: + case arg_none: + if(opts->argptr) *((int*)aptr) += 1; // increment value + break; + case arg_int: + result = myatoll(aptr, optarg, arg_int); + break; + case arg_longlong: + result = myatoll(aptr, optarg, arg_longlong); + break; + case arg_double: + result = myatod(aptr, optarg, arg_double); + break; + case arg_float: + result = myatod(aptr, optarg, arg_float); + break; + case arg_string: + result = (*((void**)aptr) = (void*)strdup(optarg)); + break; + case arg_function: + result = ((argfn)aptr)(optarg); + break; + } + if(!result){ + showhelp(optind, options); + } + } + *argc -= optind; + *argv += optind; +} + +/** + * compare function for qsort + * first - sort by short options; second - sort arguments without sort opts (by long options) + */ +static int argsort(const void *a1, const void *a2){ + const myoption *o1 = (myoption*)a1, *o2 = (myoption*)a2; + const char *l1 = o1->name, *l2 = o2->name; + int s1 = o1->val, s2 = o2->val; + int *f1 = o1->flag, *f2 = o2->flag; + // check if both options has short arg + if(f1 == NULL && f2 == NULL){ // both have short arg + return (s1 - s2); + }else if(f1 != NULL && f2 != NULL){ // both don't have short arg - sort by long + return strcmp(l1, l2); + }else{ // only one have short arg -- return it + if(f2) return -1; // a1 have short - it is 'lesser' + else return 1; + } +} + +/** + * Show help information based on myoption->help values + * @param oindex (i) - if non-negative, show only help by myoption[oindex].help + * @param options (i) - array of `myoption` + * + * @exit: run `exit(-1)` !!! + */ +void showhelp(int oindex, myoption *options){ + int max_opt_len = 0; // max len of options substring - for right indentation + const int bufsz = 255; + char buf[bufsz+1]; + myoption *opts = options; + assert(opts); + assert(opts[0].name); // check whether there is at least one options + if(oindex > -1){ // print only one message + opts = &options[oindex]; + printf(" "); + if(!opts->flag && isalpha(opts->val)) printf("-%c, ", opts->val); + printf("--%s", opts->name); + if(opts->has_arg == 1) printf("=arg"); + else if(opts->has_arg == 2) printf("[=arg]"); + printf(" %s\n", _(opts->help)); + exit(-1); + } + // header, by default is just "progname\n" + printf("\n"); + if(strstr(helpstring, "%s")) // print progname + printf(helpstring, __progname); + else // only text + printf("%s", helpstring); + printf("\n"); + // count max_opt_len + do{ + int L = strlen(opts->name); + if(max_opt_len < L) max_opt_len = L; + }while((++opts)->name); + max_opt_len += 14; // format: '-S , --long[=arg]' - get addition 13 symbols + opts = options; + // count amount of options + int N; for(N = 0; opts->name; ++N, ++opts); + if(N == 0) exit(-2); + // Now print all help (sorted) + opts = options; + qsort(opts, N, sizeof(myoption), argsort); + do{ + int p = sprintf(buf, " "); // a little indent + if(!opts->flag && opts->val) // .val is short argument + p += snprintf(buf+p, bufsz-p, "-%c, ", opts->val); + p += snprintf(buf+p, bufsz-p, "--%s", opts->name); + if(opts->has_arg == 1) // required argument + p += snprintf(buf+p, bufsz-p, "=arg"); + else if(opts->has_arg == 2) // optional argument + p += snprintf(buf+p, bufsz-p, "[=arg]"); + assert(p < max_opt_len); // there would be magic if p >= max_opt_len + printf("%-*s%s\n", max_opt_len+1, buf, _(opts->help)); // write options & at least 2 spaces after + ++opts; + }while(--N); + printf("\n\n"); + exit(-1); +} + +/** + * get suboptions from parameter string + * @param str - parameter string + * @param opt - pointer to suboptions structure + * @return TRUE if all OK + */ +bool get_suboption(char *str, mysuboption *opt){ + int findsubopt(char *par, mysuboption *so){ + int idx = 0; + if(!par) return -1; + while(so[idx].name){ + if(strcasecmp(par, so[idx].name) == 0) return idx; + ++idx; + } + return -1; // badarg + } + bool opt_setarg(mysuboption *so, int idx, char *val){ + mysuboption *soptr = &so[idx]; + bool result = FALSE; + void *aptr = soptr->argptr; + switch(soptr->type){ + default: + case arg_none: + if(soptr->argptr) *((int*)aptr) += 1; // increment value + result = TRUE; + break; + case arg_int: + result = myatoll(aptr, val, arg_int); + break; + case arg_longlong: + result = myatoll(aptr, val, arg_longlong); + break; + case arg_double: + result = myatod(aptr, val, arg_double); + break; + case arg_float: + result = myatod(aptr, val, arg_float); + break; + case arg_string: + result = (*((void**)aptr) = (void*)strdup(val)); + break; + case arg_function: + result = ((argfn)aptr)(val); + break; + } + return result; + } + char *tok; + bool ret = FALSE; + char *tmpbuf; + tok = strtok_r(str, ":,", &tmpbuf); + do{ + char *val = strchr(tok, '='); + int noarg = 0; + if(val == NULL){ // no args + val = "1"; + noarg = 1; + }else{ + *val++ = '\0'; + if(!*val || *val == ':' || *val == ','){ // no argument - delimeter after = + val = "1"; noarg = 1; + } + } + int idx = findsubopt(tok, opt); + if(idx < 0){ + /// "Неправильный параметр: %s" + WARNX(_("Wrong parameter: %s"), tok); + goto returning; + } + if(noarg && opt[idx].has_arg == NEED_ARG){ + /// "%s: необходим аргумент!" + WARNX(_("%s: argument needed!"), tok); + goto returning; + } + if(!opt_setarg(opt, idx, val)){ + /// "Неправильный аргумент \"%s\" параметра \"%s\"" + WARNX(_("Wrong argument \"%s\" of parameter \"%s\""), val, tok); + goto returning; + } + }while((tok = strtok_r(NULL, ":,", &tmpbuf))); + ret = TRUE; +returning: + return ret; +} diff --git a/parseargs.h b/parseargs.h new file mode 100644 index 0000000..b8351d9 --- /dev/null +++ b/parseargs.h @@ -0,0 +1,124 @@ +/* + * parseargs.h - headers for parsing command line arguments + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __PARSEARGS_H__ +#define __PARSEARGS_H__ + +#include // bool +#include + +#ifndef TRUE + #define TRUE true +#endif + +#ifndef FALSE + #define FALSE false +#endif + +// macro for argptr +#define APTR(x) ((void*)x) + +// if argptr is a function: +typedef bool(*argfn)(void *arg); + +/* + * type of getopt's argument + * WARNING! + * My function change value of flags by pointer, so if you want to use another type + * make a latter conversion, example: + * char charg; + * int iarg; + * myoption opts[] = { + * {"value", 1, NULL, 'v', arg_int, &iarg, "char val"}, ..., end_option}; + * ..(parse args).. + * charg = (char) iarg; + */ +typedef enum { + arg_none = 0, // no arg + arg_int, // integer + arg_longlong, // long long + arg_double, // double + arg_float, // float + arg_string, // char * + arg_function // parse_args will run function `bool (*fn)(char *optarg, int N)` +} argtype; + +/* + * Structure for getopt_long & help + * BE CAREFUL: .argptr is pointer to data or pointer to function, + * conversion depends on .type + * + * ATTENTION: string `help` prints through macro PRNT(), bu default it is gettext, + * but you can redefine it before `#include "parseargs.h"` + * + * if arg is string, then value wil be strdup'ed like that: + * char *str; + * myoption opts[] = {{"string", 1, NULL, 's', arg_string, &str, "string val"}, ..., end_option}; + * *(opts[1].str) = strdup(optarg); + * in other cases argptr should be address of some variable (or pointer to allocated memory) + * + * NON-NULL argptr should be written inside macro APTR(argptr) or directly: (void*)argptr + * + * !!!LAST VALUE OF ARRAY SHOULD BE `end_option` or ZEROS !!! + * + */ +typedef enum{ + NO_ARGS = 0, // first three are the same as in getopt_long + NEED_ARG = 1, + OPT_ARG = 2, + MULT_PAR +} hasarg; + +typedef struct{ + // these are from struct option: + const char *name; // long option's name + hasarg has_arg; // 0 - no args, 1 - nesessary arg, 2 - optionally arg, 4 - need arg & key can repeat (args are stored in null-terminated array) + int *flag; // NULL to return val, pointer to int - to set its value of val (function returns 0) + int val; // short opt name (if flag == NULL) or flag's value + // and these are mine: + argtype type; // type of argument + void *argptr; // pointer to variable to assign optarg value or function `bool (*fn)(char *optarg, int N)` + const char *help; // help string which would be shown in function `showhelp` or NULL +} myoption; + +/* + * Suboptions structure, almost the same like myoption + * used in parse_subopts() + */ +typedef struct{ + const char *name; + hasarg has_arg; + argtype type; + void *argptr; +} mysuboption; + +// last string of array (all zeros) +#define end_option {0,0,0,0,0,0,0} +#define end_suboption {0,0,0,0} + +extern const char *__progname; + +void showhelp(int oindex, myoption *options); +void parseargs(int *argc, char ***argv, myoption *options); +void change_helpstring(char *s); +bool get_suboption(char *str, mysuboption *opt); + +#endif // __PARSEARGS_H__ diff --git a/saveimg.c b/saveimg.c new file mode 100644 index 0000000..8a8a7c8 --- /dev/null +++ b/saveimg.c @@ -0,0 +1,294 @@ +/* + * saveimg.c - functions to save data in png and FITS formats + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include "usefull_macros.h" +#include "cmdlnopts.h" // for flag "-r", which will tell to rewrite existing file +#include "saveimg.h" +#include +#include +#include + +#define Stringify(x) #x +#define OMP_FOR(x) _Pragma(Stringify(omp parallel for x)) + +#ifndef THREAD_NUMBER + #define THREAD_NUMBER 4 +#endif + +#ifndef OMP_NUM_THREADS + #define OMP_NUM_THREADS THREAD_NUMBER +#endif + + +#define TRYFITS(f, ...) \ +do{int status = 0; f(__VA_ARGS__, &status); \ + if (status){ ret = 0; \ + fits_report_error(stderr, status); \ + goto returning;} \ +}while(0) +#define WRITEKEY(...) \ +do{ int status = 0; \ + fits_write_key(fp, __VA_ARGS__, &status);\ + if(status) \ + fits_report_error(stderr, status); \ +}while(0) + +/** + * Create filename as outfile + number + "." + suffix + * number -- any number from 1 to 9999 + * This function simply returns "outfile.suffix" when "-r" option is set + * + * @param outfile - file name + * @param suffix - file suffix + * @return created filename or NULL + */ +char *createfilename(char* outfile, char* suffix){ + FNAME(); + struct stat filestat; + char buff[256], sfx[32]; + if(suffix) snprintf(sfx, 31, ".%s", suffix); + else sfx[0] = 0; // no suffix + if(rewrite_ifexists){ // there was key "-f": simply return copy of filename + if(snprintf(buff, 255, "%s%s", outfile, sfx) < 1){ + DBG("error"); + return NULL; + } + DBG("(-r present) filename: %s", buff); + return strdup(buff); + } + int num; + if(!outfile) outfile = ""; + for(num = 1; num < 10000; num++){ + if(snprintf(buff, 255, "%s_%04d%s", outfile, num, sfx) < 1){ + DBG("error"); + return NULL; + } + if(stat(buff, &filestat)){ // || !S_ISREG(filestat.st_mode)) // OK, file not exists + DBG("filename: %s", buff); + return strdup(buff); + } + } + DBG("n: %s\n", buff); + WARN("Oops! All numbers are busy or other error!"); + return NULL; +} + +typedef struct{ + double *image; + double min; + double max; + double avr; + double std; +} ImStat; + +static ImStat glob_stat; + +/** + * compute basics image statictics + * @param img - image data + * @param size - image size W*H + * @return + */ +void get_stat(double *img, size_t size){ + FNAME(); + if(glob_stat.image == img) return; + size_t i; + double pv, sum=0., sum2=0., sz=(double)size; + double max = -1., min = 1e15; + for(i = 0; i < size; i++){ + pv = (double) *img++; + sum += pv; + sum2 += (pv * pv); + if(max < pv) max = pv; + if(min > pv) min = pv; + } + glob_stat.image = img; + glob_stat.avr = sum/sz; + glob_stat.std = sqrt(fabs(sum2/sz - glob_stat.avr*glob_stat.avr)); + glob_stat.max = max; + glob_stat.min = min; + DBG("Image stat: max=%g, min=%g, avr=%g, std=%g", max, min, glob_stat.avr, glob_stat.std); +} + +/** + * Save data to fits file + * @param filename - filename to save to + * @param sz - image size: sz x sz + * @param imbox - image bounding box + * @data image data + * @return 0 if failed + */ +int writefits(char *filename, size_t sz, mirror *mir, double *data){ + FNAME(); + long naxes[2] = {sz, sz}; + static char* newname = NULL; + char buf[80]; + int ret = 1; + time_t savetime = time(NULL); + fitsfile *fp; + if(!filename) return 0; + newname = realloc(newname, strlen(filename + 2)); + sprintf(newname, "!%s", filename); // say cfitsio that file could be rewritten + TRYFITS(fits_create_file, &fp, newname); + TRYFITS(fits_create_img, fp, FLOAT_IMG, 2, naxes); + // FILE / Input file original name + WRITEKEY(TSTRING, "FILE", filename, "Input file original name"); + WRITEKEY(TSTRING, "DETECTOR", "Wavefront model", "Detector model"); + // IMAGETYP / object, flat, dark, bias, scan, eta, neon, push + WRITEKEY(TSTRING, "IMAGETYP", "object", "Image type"); + // DATAMAX, DATAMIN / Max,min pixel value + WRITEKEY(TDOUBLE, "DATAMAX", &glob_stat.max, "Max data value"); + WRITEKEY(TDOUBLE, "DATAMIN", &glob_stat.min, "Min data value"); + // Some Statistics + WRITEKEY(TDOUBLE, "DATAAVR", &glob_stat.avr, "Average data value"); + WRITEKEY(TDOUBLE, "DATASTD", &glob_stat.std, "Standart deviation of data value"); + // DATE / Creation date (YYYY-MM-DDThh:mm:ss, UTC) + strftime(buf, 79, "%Y-%m-%dT%H:%M:%S", gmtime(&savetime)); + WRITEKEY(TSTRING, "DATE", buf, "Creation date (YYYY-MM-DDThh:mm:ss, UTC)"); + // DATE-OBS / DATE OF OBS. + WRITEKEY(TSTRING, "DATE-OBS", buf, "DATE OF OBS. (YYYY-MM-DDThh:mm:ss, local)"); + // OBJECT / Object name + WRITEKEY(TSTRING, "OBJECT", "Modeled object", "Object name"); + // BINNING / Binning + WRITEKEY(TSTRING, "XBIN", "1", "Horizontal binning"); + WRITEKEY(TSTRING, "YBIN", "1", "Vertical binning"); + // PROG-ID / Observation program identifier + WRITEKEY(TSTRING, "PROG-ID", "Wavefront modeling", "Observation program identifier"); + // AUTHOR / Author of the program + WRITEKEY(TSTRING, "AUTHOR", "Edward V. Emelianov", "Author of the program"); + if(mir){ + WRITEKEY(TDOUBLE, "ZBESTFOC", &mir->zbestfoc, "Z coordinate of best focus"); + WRITEKEY(TDOUBLE, "Z70", &mir->z07, "Z coordinate of best PSF on 70\%"); + } + TRYFITS(fits_write_img, fp, TDOUBLE, 1, sz * sz, data); + TRYFITS(fits_close_file, fp); +returning: + return ret; +} + +static uint8_t *rowptr = NULL; +uint8_t *processRow(double *irow, size_t width, double min, double wd){ + FREE(rowptr); + //float umax = ((float)(UINT16_MAX-1)); + rowptr = MALLOC(uint8_t, width * 3); + OMP_FOR(shared(irow)) + for(size_t i = 0; i < width; i++){ + double gray = (irow[i] - min) / wd; + if(gray == 0.) continue; + int G = (int)(gray * 4.); + double x = 4.*gray - (double)G; + uint8_t *ptr = &rowptr[i*3]; + uint8_t r = 0, g = 0, b = 0; + switch(G){ + case 0: + g = (uint8_t)(255. * x + 0.5); + b = 255; + break; + case 1: + g = 255; + b = (uint8_t)(255. * (1. - x) + 0.5); + break; + case 2: + r = (uint8_t)(255. * x + 0.5); + g = 255; + break; + case 3: + r = 255; + g = (uint8_t)(255. * (1. - x) + 0.5); + break; + default: + r = 255; + } + ptr[0] = r; ptr[1] = g; ptr[2] = b; + //ptr[0] = ptr[1] = ptr[2] = gray*255; + } + return rowptr; +} + +int writepng(char *filename, size_t sz, double *data){ + FNAME(); + int ret = 1; + FILE *fp = NULL; + png_structp pngptr = NULL; + png_infop infoptr = NULL; + double min = glob_stat.min, wd = glob_stat.max - min; + double *row; + + if ((fp = fopen(filename, "w")) == NULL){ + perror("Can't open png file"); + ret = 0; + goto done; + } + if ((pngptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, + NULL, NULL, NULL)) == NULL){ + perror("Can't create png structure"); + ret = 0; + goto done; + } + if ((infoptr = png_create_info_struct(pngptr)) == NULL){ + perror("Can't create png info structure"); + ret = 0; + goto done; + } + png_init_io(pngptr, fp); + png_set_compression_level(pngptr, 1); + png_set_IHDR(pngptr, infoptr, sz, sz, 8, PNG_COLOR_TYPE_RGB,//16, PNG_COLOR_TYPE_GRAY, + PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, + PNG_FILTER_TYPE_DEFAULT); + png_write_info(pngptr, infoptr); + png_set_swap(pngptr); + size_t height = sz; + for(row = &data[sz*(height-1)]; height > 0; row -= sz, height--) + png_write_row(pngptr, (png_bytep)processRow(row, sz, min, wd)); + png_write_end(pngptr, infoptr); +done: + if(fp) fclose(fp); + if(pngptr) png_destroy_write_struct(&pngptr, &infoptr); + return ret; +} + +/** + * Save data to image file[s] with format t + * @param name - filename prefix or NULL to save to "outXXXX.format" + * @param t - image[s] type[s] + * @param width, height - image size + * @param imbox - image bounding box (for FITS header) + * @param mirror - mirror parameters (for FITS header) + * @param data image data + * @return number of saved images + */ +int writeimg(char *name, size_t sz, mirror *mir, double *data){ + FNAME(); + char *filename = NULL; + int ret = 0; + get_stat(data, sz*sz); + filename = createfilename(name, "fits"); + if(filename){ + ret = writefits(filename, sz, mir, data); + FREE(filename); + } + filename = createfilename(name, "png"); + if(filename){ + ret += writepng(filename, sz, data); + FREE(filename); + } + return ret; +} diff --git a/saveimg.h b/saveimg.h new file mode 100644 index 0000000..139a9ad --- /dev/null +++ b/saveimg.h @@ -0,0 +1,33 @@ +/* + * saveimg.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __SAVEIMG_H__ +#define __SAVEIMG_H__ + +#include +#include "spots.h" + +int writeimg(char *name, size_t sz, mirror *mir, double *data); + +char *createfilename(char* outfile, char* suffix); + +#endif // __SAVEIMG_H__ diff --git a/simple_list.c b/simple_list.c new file mode 100644 index 0000000..479906a --- /dev/null +++ b/simple_list.c @@ -0,0 +1,93 @@ +/* + * simple_list.c - simple one-direction list + * + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include +#include + +#include "simple_list.h" + +#define MALLOC alloc_simple_list +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +static void *alloc_simple_list(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} + +void (*listdata_free)(void *node_data) = NULL; // external function to free(listdata) + +#define FREE(arg) do{free(arg); arg = NULL;}while(0) + +/** + * add element v to list with root root (also this can be last element) + * @param root (io) - pointer to root (or last element) of list or NULL + * if *root == NULL, just created node will be placed there + * @param v - data inserted + * @return pointer to created node + */ +List *list_add(List **root, void *v){ + List *node, *last; + if((node = (List*) MALLOC(1, sizeof(List))) == 0) return NULL; // allocation error + node->data = v; // insert data + if(root){ + if(*root){ // there was root node - search last + last = *root; + while(last->next) last = last->next; + last->next = node; // insert pointer to new node into last element in list + } + if(!*root) *root = node; + } + return node; +} + +/** + * set listdata_free() + * @param fn - function that freec memory of (node) + */ +void listfree_function(void (*fn)(void *node)){ + listdata_free = fn; +} + +/** + * remove all nodes in list + * @param root - pointer to root node + */ +void list_free(List **root){ + List *node = *root, *next; + if(!root || !*root) return; + do{ + next = node->next; + if(listdata_free) + listdata_free(node->data); + free(node); + node = next; + }while(node); + *root = NULL; +} + diff --git a/simple_list.h b/simple_list.h new file mode 100644 index 0000000..95cc2e6 --- /dev/null +++ b/simple_list.h @@ -0,0 +1,39 @@ +/* + * simple_list.h - header file for simple list support + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __SIMPLE_LIST_H__ +#define __SIMPLE_LIST_H__ + +#define LIST_T(x) ((void*) x) + +typedef struct list_{ + void *data; + struct list_ *next; +} List; + +// add element v to list with root root (also this can be last element) +List *list_add(List **root, void *v); +// set listdata_free() +void listfree_function(void (*fn)(void *node)); +void list_free(List **root); + +#endif // __SIMPLE_LIST_H__ diff --git a/spots.c b/spots.c new file mode 100644 index 0000000..3393501 --- /dev/null +++ b/spots.c @@ -0,0 +1,556 @@ +/* + * spots.c + * + * Copyright 2015 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include // stat +#include // mmap +#include +#include +#include +#include +#include +#include +#include +#include "spots.h" +#include "cmdlnopts.h" +#include "usefull_macros.h" + +#define MM_TO_ARCSEC(x) (x*206265./FOCAL_R) + +const double _4F = 4. * FOCAL_R; +const double _2F = 2. * FOCAL_R; + +double distance = -1.; // distanse between pre- and postfocal images in millimeters + +void spots_free(List **spots){ + listfree_function(free); + list_free(spots); + listfree_function(NULL); +} + +/** + * Read spots-file, find center of hartmannogram & convert coordinates + * @param filename (i) - name of spots-file + * @return dynamically allocated hartmanogram structure + * COORDINATES ARE IN MILLIMETERS!!! + */ +hartmann *read_spots(char *filename){ + assert(filename); + mmapbuf *M = NULL; + M = My_mmap(filename); + hartmann *H = MALLOC(hartmann, 1); + point *spots = H->spots; + uint8_t *got = H->got; + H->filename = strdup(filename); + char *pos = M->data, *epos = pos + M->len; + // readout list of spots + for(; pos && pos < epos; pos = strchr(pos+1, '\n')){ + double x, y; + int id, a, b; + //if(3 != sscanf(pos, "%d %*s %*s %*s %*s %*s %*s %*s %*s %lf %lf", &id, &x, &y)) + if(3 != sscanf(pos, "%d %*s %*s %*s %*s %lf %lf", &id, &x, &y)) + continue; + a = id/100; b = id%100; + if(b < 32) id = a*32 + b; // main spots + else if(a == 2) id = 256; // inner marker + else id = 257; // outern marker + spots[id].x = x; +#ifdef MIR_Y + spots[id].y = -y; +#else + spots[id].y = y; +#endif + got[id] = 1; + }; + // get center: simply get center of each pair of opposite spots & calculate average + // IDs: xyy -- x(yy+16), yy=[0..15] + double xc = 0., yc = 0.; + int i, j, N = 0; + for(i = 0; i < 8; i++){ // circles + int id0 = i*32, id1 = id0 + 16; + for(j = 0; j < 16; j++, id0++, id1++){ // first half + if(!got[id0] || !got[id1]) continue; // there's no such pair + double c1, c2; + c1 = (spots[id0].x + spots[id1].x) / 2.; + c2 = (spots[id0].y + spots[id1].y) / 2.; + xc += c1; + yc += c2; + //printf("center: %g, %g\n", c1, c2); + ++N; + } + } + xc /= (double) N; + yc /= (double) N; + //printf("Calculated center: %g, %g\n", xc, yc); + H->center.x = xc * G.pixsize; + H->center.y = yc * G.pixsize; + printf("get common center: (%.1f. %.1f)mm (pixsize=%g)\n", H->center.x, H->center.y, G.pixsize); + // convert coordinates to center & fill spots array of H + for(i = 0; i < 258; i++){ + if(!got[i]) continue; + spots[i].x = (spots[i].x - xc) * G.pixsize; + spots[i].y = (spots[i].y - yc) * G.pixsize; + //printf("spot #%d: (%.1f, %.1f)mm\n", i, spots[i].x, spots[i].y); + } +#ifdef MIR_Y + const double stp = M_PI/16., an0 = -M_PI_2 + stp/2., _2pi = 2.*M_PI; +#else + const double stp = M_PI/16., an0 = +M_PI_2 - stp/2., _2pi = 2.*M_PI; +#endif + double dmean = 0.; +// printf("RAYS' theta:\n"); + for(i = 0; i < 32; i++){ // rays + int id = i, N=0; +#ifdef MIR_Y + double sum = 0., refang = an0 + stp * (double)i; +#else + double sum = 0., refang = an0 - stp * (double)i; +#endif +// printf("ray %d: ", i); + for(j = 0; j < 8; j++, id+=32){ // circles + if(!got[id]) continue; + double theta = atan2(spots[id].y, spots[id].x); +// printf("%.5f,", theta); + sum += theta; + ++N; + } + double meanang = sum/(double)N, delta = refang-meanang; + if(delta > _2pi) delta -= _2pi; + if(delta < 0.) delta += _2pi; + if(delta < 0.) delta += _2pi; +// printf("mean: %g, delta: %g\n", meanang, delta); + dmean += delta; + } + dmean /= 32.; + printf("MEAN delta: %g\n\n", dmean); + //if(!prefocal) dmean += M_PI; + double s,c; + sincos(dmean, &s, &c); +// printf("sin: %g, cos: %g\n",s,c); + // rotate data to ZERO + for(i = 0; i < 258; i++){ + if(!got[i]) continue; + double x = spots[i].x, y = spots[i].y; + spots[i].x = x*c+y*s; + spots[i].y = -x*s+y*c; + } + My_munmap(M); + return H; +} + +/** + * Calculate coordinates of points on mirror + * Modify coordinates of hartmannogram centers: + * !!! the center beam on prefocal hartmannogram will have zero coordinates !!! + * @param H (io) - array of pre- and postfocal H + * @param D - distance from mirror top to prefocal image in millimeters + * @return mirror structure with coordinates & tan parameters + */ +mirror *calc_mir_coordinates(hartmann **H){ + assert(H); assert(H[0]); assert(H[1]); + point *pre = H[0]->spots, *post = H[1]->spots, *prec = &H[0]->center, *postc = &H[1]->center; + mirror *mir = MALLOC(mirror, 1); + uint8_t *got_pre = H[0]->got, *got_post = H[1]->got, *mirgot = mir->got; + double *Z = MALLOC(double, 259); // 258 points + center beam + int i, iter; + point *tans = mir->tans, *spots = mir->spots, *mirc = &mir->center; + H[1]->center.x -= H[0]->center.x; + H[1]->center.y -= H[0]->center.y; + H[0]->center.x = 0.; + H[0]->center.y = 0.; + double dx = H[1]->center.x, dy = H[1]->center.y; + printf("H1 center: (%g, %g)\n", dx, dy); + double FCCnumerator = 0., FCCdenominator = 0.; // components of Fcc (Focus for minimal circle of confusion): + /* + * SUM (x_pre * tans_x + y_pre * tans_y) + * Fcc = ----------------------------------------- + * SUM (tanx_x^2 + tans_y^2) + */ + for(i = 0; i < 258; i++){ // check point pairs on pre/post + if(got_pre[i] && got_post[i]){ + double tx, ty; + mirgot[i] = 1; + ++mir->spotsnum; + //tx = (dx + post[i].x - pre[i].x)/distance; + tx = (post[i].x - pre[i].x)/distance; + tans[i].x = tx; + //ty = (dy + post[i].y - pre[i].y)/distance; + ty = (post[i].y - pre[i].y)/distance; + tans[i].y = ty; + FCCnumerator += pre[i].x*tx + pre[i].y*ty; + FCCdenominator += tx*tx + ty*ty; + } + } + mir->zbestfoc = -FCCnumerator/FCCdenominator; + double D = FOCAL_R - mir->zbestfoc; + printf("BEST (CC) focus: %g; D = %g\n", mir->zbestfoc, D); + printf("\nCalculate spots centers projections to mirror surface\n"); + point tanc; + tanc.x = (dx + postc->x - prec->x)/distance; + tanc.y = (dy + postc->y - prec->y)/distance; + memcpy(&mir->tanc, &tanc, sizeof(point)); + double Zerr = 10.; + /* + * X = x_pre + (Z-D)*tans_x + * Y = y_pre + (Z-D)*tans_y + * Z = (X^2 + Y^2) / (4F) + */ + for(iter = 0; iter < 10 && Zerr > 1e-6; iter++){ // calculate points position on mirror + Zerr = 0.; + double x, y; + for(i = 0; i < 258; i++){ + if(!mirgot[i]) continue; + x = pre[i].x + tans[i].x * (Z[i] - D); + y = pre[i].y + tans[i].y * (Z[i] - D); + spots[i].x = x; + spots[i].y = y; + double newZ = (x*x + y*y)/_4F; + double d = newZ - Z[i]; + Zerr += d*d; + Z[i] = newZ; + } + x = prec->x + tanc.x * (Z[258] - D); + y = prec->y + tanc.y * (Z[258] - D); + mirc->x = x; mirc->y = y; + Z[258] = (x*x + y*y)/_4F; + printf("iteration %d, sum(dZ^2) = %g\n", iter, Zerr); + } + // now calculate polar coordinates relative to calculated center + double xc = mirc->x, yc = mirc->y; + polar *rtheta = mir->pol_spots; +// double Rmax = 0.; + for(i = 0; i < 258; i++){ + if(!mirgot[i]) continue; + double x = spots[i].x - xc, y = spots[i].y - yc; + double r = sqrt(x*x + y*y); + rtheta[i].r = r; + rtheta[i].theta = atan2(y, x); +// if(Rmax < r) Rmax = r; + } +// Rmax = 3000.; + for(i = 0; i < 258; i++){ + if(mirgot[i]) + // rtheta[i].r /= Rmax; + rtheta[i].r /= MIR_R; + } +// mir->Rmax = Rmax; + FREE(Z); + return mir; +} + +/** + * Calculate Hartmann constant + * @param mir (i) - filled mirror structure + * @param prefoc (i) - prefocal hartmannogram + * @return constant value + * + * Hartmann constant is + * SUM(r_i^2*|F_i-F|) 200000 + * T = -------------------- * -------- + * SUM(r_i) F_m^2 + * + * where: + * F_m - mirror focus ratio + * r_i - radius of ith zone + * F_i - its focus (like in calc_mir_coordinates), counted from prefocal image + * SUM(r_i * F_i) + * F = ------------------ - mean focus (counted from prefocal image) + * SUM(r_i) + */ +double calc_Hartmann_constant(mirror *mir){ + uint8_t *mirgot = mir->got; + int i, j; + point *tans = mir->tans, *pre = mir->prefocal; + polar *p_spots = mir->pol_spots; + double foc_i[8], r_i[8], Rsum = 0.; + double FCCnumerator = 0., FCCdenominator = 0.; // components of Fcc (Focus for minimal circle of confusion): + for(j = 0; j < 8; ++j){ // cycle by circles + double Rj = 0.; + int Nj = 0, idx = j*32; + for(i = 0; i < 32; ++i, ++idx){ // run through points on same circle + if(!mirgot[idx]) continue; + ++Nj; + Rj += p_spots[idx].r; + double tx = tans[idx].x, ty = tans[idx].y; + FCCnumerator += pre[idx].x*tx + pre[idx].y*ty; + FCCdenominator += tx*tx + ty*ty; + } + foc_i[j] = -FCCnumerator/FCCdenominator; + r_i[j] = Rj / Nj; + Rsum += r_i[j]; + printf("focus on R = %g is %g\n", MIR_R * r_i[j], foc_i[j]); + } + double F = 0.; + for(j = 0; j < 8; ++j){ + F += foc_i[j] * r_i[j]; + } + F /= Rsum; + double numerator = 0.; + for(j = 0; j < 8; ++j){ + numerator += r_i[j]*r_i[j]*fabs(foc_i[j]-F); + } + printf("Mean focus is %g, numerator: %g, Rsum = %g\n", F, numerator, Rsum); + // multiply by MIR_R because r_i are normed by R + double T = MIR_R * 2e5/FOCAL_R/FOCAL_R*numerator/Rsum; + printf("\nHartmann value T = %g\n", T); + return T; +} + +static int cmpdbl(const void * a, const void * b){ + if (*(double*)a > *(double*)b) + return 1; + else return -1; +} + +/** + * Calculate energy in circle of confusion + * @param mir (i) - filled mirror structure + * @param hartm (i) - zero-terminated list of hartmannograms + */ +void getQ(mirror *mir){ + point *pre = mir->prefocal, *tans = mir->tans; + uint8_t *got = mir->got; + int i, j, N = mir->spotsnum, N03 = 0.3*N, N05 = 0.5*N, N07 = 0.7*N, N09 = 0.9*N; + double zstart = mir->zbestfoc - 3., z, zend = mir->zbestfoc + 3.01; + double *R = MALLOC(double, N); + double z03, z05, z07, z09, r03=1e3, r05=1e3, r07=1e3, r09=1e3; + printf("\nEnergy in circle of confusion\n"); + printf("z mean(R)'' std(R)'' R0.3'' R0.5'' R0.7'' R0.9'' Rmax''\n"); + for(z = zstart; z < zend; z += 0.1){ + j = 0; + double Rsum = 0., R2sum = 0.; + for(i = 0; i < 258; i++){ + if(!got[i]) continue; + double x, y, R2, R1; + x = pre[i].x + tans[i].x * z; + y = pre[i].y + tans[i].y * z; + R2 = x*x + y*y; + R1 = sqrt(R2); + R[j] = R1; + R2sum += R2; + Rsum += R1; + ++j; + } + qsort(R, N, sizeof(double), cmpdbl); + double R3 = R[N03], R5 = R[N05], R7 = R[N07], R9 = R[N09]; +/* for(i = N-1; i; --i) if(R[i] < R7) break; + R7 = R[i]; + for(i = N-1; i; --i) if(R[i] < R9) break; + R9 = R[i];*/ + printf("%6.2f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f\n", + z, MM_TO_ARCSEC(Rsum/N), MM_TO_ARCSEC(sqrt((R2sum - Rsum*Rsum/N)/(N-1.))), + MM_TO_ARCSEC(R3), MM_TO_ARCSEC(R5), + MM_TO_ARCSEC(R7), MM_TO_ARCSEC(R9), MM_TO_ARCSEC(R[N-1])); + if(r03 > R3){ + r03 = R3; z03 = z; + } + if(r05 > R5){ + r05 = R5; z05 = z; + } + if(r07 > R7){ + r07 = R7; z07 = z; + } + if(r09 > R9){ + r09 = R9; z09 = z; + } + } + mir->z07 = z07; + printf("\ngot best values: z03=%g (r=%g''), z05=%g (r=%g''), z07=%g (r=%g''), z09=%g (r=%g'')\n", + z03, MM_TO_ARCSEC(r03), z05, MM_TO_ARCSEC(r05), + z07, MM_TO_ARCSEC(r07), z09, MM_TO_ARCSEC(r09)); + printf("\nEnergy for z = %g\n R,'' q(r)\n", z07); + for(j = 0, i = 0; i < 258; ++i){ + if(!got[i]) continue; + double x, y; + x = pre[i].x + tans[i].x * z07; + y = pre[i].y + tans[i].y * z07; + R[j] = sqrt(x*x + y*y); + ++j; + } + qsort(R, N, sizeof(double), cmpdbl); + for(i = 0; i < N; ++i){ + printf("%8.6f %8.6f\n", MM_TO_ARCSEC(R[i]), (1.+(double)i)/N); + } + FREE(R); +} + +void h_free(hartmann **H){ + if(!H || !*H) return; + FREE((*H)->filename); + FREE(*H); +} + +/** + * Calculate spot diagram for given z value + * @param mir (i) - filled mirror structure + * @param prefoc (i) - prefocal hartmannogram + * @return allocated structure of spot diagram + */ +spot_diagram *calc_spot_diagram(mirror *mir, double z){ + int i; + spot_diagram *SD = MALLOC(spot_diagram, 1); + point *pre = mir->prefocal, *spots = SD->spots; + point *tans = mir->tans; + uint8_t *got = mir->got; + memcpy(SD->got, mir->got, 258); + SD->center.x = mir->prefoc_center.x + mir->tanc.x * z; + SD->center.y = mir->prefoc_center.y + mir->tanc.y * z; + printf("spots center: (%g, %g)\n", SD->center.x, SD->center.y); + printf("\nSpot diagram for z = %g (all in mm)\n ray# x y\n", z); + for(i = 0; i < 258; ++i){ + if(!got[i]) continue; + spots[i].x = pre[i].x + tans[i].x * z; + spots[i].y = pre[i].y + tans[i].y * z; + printf("%4d %10.6f %10.6f\n", i, spots[i].x, spots[i].y); + } + return SD; +} + +/** + * Calculate gradients of mirror surface aberrations + * @param mir (io) - mirror (mir->grads will be filled by gradients) + * @param foc_spots (i) - spot diagram for best focus + */ +void calc_gradients(mirror *mir, spot_diagram *foc_spots){ + int i; + point *grads = mir->grads, *spots = foc_spots->spots; + printf("Common tilt: d/dx = %g mkm/m, d/dy = %g mkm/m\n", + -foc_spots->center.x / _2F * 1e6, -foc_spots->center.y / _2F * 1e6); + uint8_t *got = mir->got; + for(i = 0; i < 258; ++i){ + if(!got[i]) continue; + grads[i].x = -(spots[i].x) / _2F; + grads[i].y = (spots[i].y) / _2F; + } +} + +#if 0 +/** + * + * @param H (i) - array of thwo hartmannograms (0 - prefocal, 1 - postfocal) + * @param coords (o) - array of gradients' coordinates on prefocal image (allocated here) - the same as H[0]->spots coordinates + * @param grads (o) - gradients' array (allocated here) + * @param scale (o) - scale of polar coordinate R (== Rmax) + * @return size of built arrays + */ +size_t get_gradients(hartmann *H[], polar **coords, point **grads, double *scale){ + size_t Sz = 0, i, j, L0, L1; + assert(H); assert(H[0]); assert(H[1]); + spot *S0 = H[0]->spots, *S1; + List *CG_list = NULL, *curCG = NULL; + double Scale = G.pixsize * 1e-6 / distance / 2., S = 0.; // tg(2a)=dx/D -> a \approx dx/(2D) + L0 = H[0]->len; + L1 = H[1]->len; + /* + * Both lists have sorted structure + * but they can miss some points - that's why we should find exact points. + * To store dinamycally data I use List + */ + for(i = 0; i < L0; ++i, ++S0){ + int Id0 = S0->id; + for(S1 = H[1]->spots, j=0; j < L1; ++j, ++S1){ + if(S1->id > Id0) break; // point with Id0 not found + if(S1->id == Id0){ + CG *cg = MALLOC(CG, 1); + cg->id = Id0; + cg->x = S0->x; cg->y = S0->y; + cg->Dx = -(S1->x - S0->x)*Scale; + cg->Dy = -(S1->y - S0->y)*Scale; + Sz++; + if(CG_list) + curCG = list_add(&curCG, cg); + else + curCG = list_add(&CG_list, cg); + break; + } + } + } + polar *C = MALLOC(polar, Sz), *cptr = C; + point *G = MALLOC(point, Sz), *gptr = G; + curCG = CG_list; + for(i = 0; i < Sz; i++, cptr++, gptr++, curCG = curCG->next){ + double x, y, length, R; + CG *cur = (CG*)curCG->data; + x = cur->x; y = cur->y; + R = sqrt(x*x + y*y); + cptr->r = R; + if(S < R) S = R; // find max R + cptr->theta = atan2(y, x); + x = cur->Dx; y = cur->Dy; + length = sqrt(1. + x*x + y*y); // length of vector for norm + gptr->x = x / length; + gptr->y = y / length; + } + cptr = C; + for(i = 0; i < Sz; i++, cptr++) + cptr->r /= S; + *scale = S; + *coords = C; *grads = G; + listfree_function(free); + list_free(&CG_list); + listfree_function(NULL); + return Sz; +} + +/** + * Readout of gradients file calculated somewhere outside + * + * @param coords (o) - array with coordinates on mirror (in meters), ALLOCATED HERE! + * @param grads (o) - array with gradients (meters per meters), ALLOCATED HERE! + * @param scale (o) - scale on mirror (Rmax) + * @return - amount of points readed + */ +size_t read_gradients(char *gradname, polar **coords, point **grads, double *scale){ + assert(gradname); + mmapbuf *M = NULL; + double Rmax = 0.; + M = My_mmap(gradname); + size_t L = 0; + char *pos = M->data, *epos = pos + M->len; + for(; pos && pos < epos; pos = strchr(pos+1, '\n')){ + double x, y, gx, gy, R; + if(4 != sscanf(pos, "%lf %lf %lf %lf", &x, &y, &gx, &gy)) + continue; + R = sqrt(x*x + y*y); + if(R > Rmax) Rmax = R; + L++; + } + printf("Found %zd points\n", L); + polar *C = MALLOC(polar, L), *cptr = C; + point *G = MALLOC(point, L), *gptr = G; + pos = M->data, epos = pos + M->len; + for(; pos && pos < epos; pos = strchr(pos+1, '\n')){ + double x, y, gx, gy, R; + if(4 != sscanf(pos, "%lf %lf %lf %lf", &x, &y, &gx, &gy)) + continue; + R = sqrt(x*x + y*y); + cptr->r = R / Rmax; + cptr->theta = atan2(y, x); + gptr->x = gx*1e6; + gptr->y = gy*1e6; + printf("%5.2f %5.2f %5.2f %5.2f\n",x, y, gx*1e6,gy*1e6); + cptr++, gptr++; + } + My_munmap(M); + *scale = Rmax; + *coords = C; *grads = G; + return L; +} +#endif diff --git a/spots.h b/spots.h new file mode 100644 index 0000000..522a19d --- /dev/null +++ b/spots.h @@ -0,0 +1,98 @@ +/* + * spots.h + * + * Copyright 2015 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __SPOTS_H__ +#define __SPOTS_H__ + +#include +#include "zernike.h" +#include "simple_list.h" + +extern double pixsize; +extern double distance; + +/* +// spot center +typedef struct{ + int id; // spot identificator + double x; // coordinates of center + double y; +} spot; +*/ + +// hartmannogram +typedef struct{ + char *filename; // spots-filename + uint8_t got[258]; // == 1 if there's this spot on image + point center; // coordinate of center + point spots[258]; // spotlist, markers have numbers 256 & 257 + double focus; // focus value for given image +} hartmann; + +// mirror +typedef struct{ + uint8_t got[258]; // == 1 if there's this spot on both pre-focal & post-focal images + int spotsnum; // total amount of spots used + point center; // x&y coordinates of center beam + point tanc; // tangent of center beam + point tans[258]; // tangents of beams == (r_post-r_pre)/d + point spots[258]; // x&y coordinates of spots on mirror surface (z=[x^2+y^2]/4f) + polar pol_spots[258]; // polar coordinates of spots on mirror surface according to center (norm by mirror R) + point grads[258]; // gradients to mirror surface + point prefocal[258];// coordinates of spots on prefocal image[s] (mean for many images) + point prefoc_center;// coordinate of prefocal hartmannogram[s] center + point postfocal[258]; // -//- on postfocal image[s] + double zbestfoc; // Z-coordinate (from pre-focal image) of best focus for minimal circle of confusion + double z07; // Z of best focus for q=0.7 +// double Rmax; // max R of spot (@ which pol_spots.r == 1.) +} mirror; + +// spot diagram +typedef struct{ + uint8_t got[258]; // the same as above + point center; // center beam's coordinate + point spots[258]; // spots + double z; // z from prefocal image +}spot_diagram; + +// gradients structure: point coordinates & gradient components +typedef struct{ + int id; + double x; + double y; + double Dx; + double Dy; +} CG; + +hartmann *read_spots(char *filename); +void h_free(hartmann **H); +mirror *calc_mir_coordinates(hartmann **H); +void getQ(mirror *mir); +double calc_Hartmann_constant(mirror *mir); +spot_diagram *calc_spot_diagram(mirror *mir, double z); +void calc_gradients(mirror *mir, spot_diagram *foc_spots); + +/* +size_t get_gradients(hartmann *H[], polar **coords, point **grads, double *scale); +size_t read_gradients(char *gradname, polar **coords, point **grads, double *scale); +*/ +#endif // __SPOTS_H__ diff --git a/usefull_macros.c b/usefull_macros.c new file mode 100644 index 0000000..1db9300 --- /dev/null +++ b/usefull_macros.c @@ -0,0 +1,360 @@ +/* + * usefull_macros.h - a set of usefull functions: memory, color etc + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include "usefull_macros.h" + +/** + * function for different purposes that need to know time intervals + * @return double value: time in seconds + */ +double dtime(){ + double t; + struct timeval tv; + gettimeofday(&tv, NULL); + t = tv.tv_sec + ((double)tv.tv_usec)/1e6; + return t; +} + +/******************************************************************************\ + * Coloured terminal +\******************************************************************************/ +int globErr = 0; // errno for WARN/ERR + +// pointers to coloured output printf +int (*red)(const char *fmt, ...); +int (*green)(const char *fmt, ...); +int (*_WARN)(const char *fmt, ...); + +/* + * format red / green messages + * name: r_pr_, g_pr_ + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int r_pr_(const char *fmt, ...){ + va_list ar; int i; + printf(RED); + va_start(ar, fmt); + i = vprintf(fmt, ar); + va_end(ar); + printf(OLDCOLOR); + return i; +} +int g_pr_(const char *fmt, ...){ + va_list ar; int i; + printf(GREEN); + va_start(ar, fmt); + i = vprintf(fmt, ar); + va_end(ar); + printf(OLDCOLOR); + return i; +} +/* + * print red error/warning messages (if output is a tty) + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int r_WARN(const char *fmt, ...){ + va_list ar; int i = 1; + fprintf(stderr, RED); + va_start(ar, fmt); + if(globErr){ + errno = globErr; + vwarn(fmt, ar); + errno = 0; + globErr = 0; + }else + i = vfprintf(stderr, fmt, ar); + va_end(ar); + i++; + fprintf(stderr, OLDCOLOR "\n"); + return i; +} + +static const char stars[] = "****************************************"; +/* + * notty variants of coloured printf + * name: s_WARN, r_pr_notty + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int s_WARN(const char *fmt, ...){ + va_list ar; int i; + i = fprintf(stderr, "\n%s\n", stars); + va_start(ar, fmt); + if(globErr){ + errno = globErr; + vwarn(fmt, ar); + errno = 0; + globErr = 0; + }else + i = +vfprintf(stderr, fmt, ar); + va_end(ar); + i += fprintf(stderr, "\n%s\n", stars); + i += fprintf(stderr, "\n"); + return i; +} +int r_pr_notty(const char *fmt, ...){ + va_list ar; int i; + i = printf("\n%s\n", stars); + va_start(ar, fmt); + i += vprintf(fmt, ar); + va_end(ar); + i += printf("\n%s\n", stars); + return i; +} + +/** + * Run this function in the beginning of main() to setup locale & coloured output + */ +void initial_setup(){ + // setup coloured output + if(isatty(STDOUT_FILENO)){ // make color output in tty + red = r_pr_; green = g_pr_; + }else{ // no colors in case of pipe + red = r_pr_notty; green = printf; + } + if(isatty(STDERR_FILENO)) _WARN = r_WARN; + else _WARN = s_WARN; + // Setup locale + setlocale(LC_ALL, ""); + setlocale(LC_NUMERIC, "C"); +#if defined GETTEXT_PACKAGE && defined LOCALEDIR + bindtextdomain(GETTEXT_PACKAGE, LOCALEDIR); + textdomain(GETTEXT_PACKAGE); +#endif +} + +/******************************************************************************\ + * Memory +\******************************************************************************/ +/* + * safe memory allocation for macro ALLOC + * @param N - number of elements to allocate + * @param S - size of single element (typically sizeof) + * @return pointer to allocated memory area + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) ERR("malloc"); + //assert(p); + return p; +} + +/** + * Mmap file to a memory area + * + * @param filename (i) - name of file to mmap + * @return stuct with mmap'ed file or die + */ +mmapbuf *My_mmap(char *filename){ + int fd; + char *ptr; + size_t Mlen; + struct stat statbuf; + /// "Не задано имя файла!" + if(!filename) ERRX(_("No filename given!")); + if((fd = open(filename, O_RDONLY)) < 0) + /// "Не могу открыть %s для чтения" + ERR(_("Can't open %s for reading"), filename); + if(fstat (fd, &statbuf) < 0) + /// "Не могу выполнить stat %s" + ERR(_("Can't stat %s"), filename); + Mlen = statbuf.st_size; + if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) + /// "Ошибка mmap" + ERR(_("Mmap error for input")); + /// "Не могу закрыть mmap'нутый файл" + if(close(fd)) ERR(_("Can't close mmap'ed file")); + mmapbuf *ret = MALLOC(mmapbuf, 1); + ret->data = ptr; + ret->len = Mlen; + return ret; +} + +void My_munmap(mmapbuf *b){ + if(munmap(b->data, b->len)){ + /// "Не могу munmap" + ERR(_("Can't munmap")); + } + FREE(b); +} + + +/******************************************************************************\ + * Terminal in no-echo mode +\******************************************************************************/ +static struct termios oldt, newt; // terminal flags +static int console_changed = 0; +// run on exit: +void restore_console(){ + if(console_changed) + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to previous state + console_changed = 0; +} + +// initial setup: +void setup_con(){ + if(console_changed) return; + tcgetattr(STDIN_FILENO, &oldt); + newt = oldt; + newt.c_lflag &= ~(ICANON | ECHO); + if(tcsetattr(STDIN_FILENO, TCSANOW, &newt) < 0){ + /// "Не могу настроить консоль" + WARN(_("Can't setup console")); + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); + signals(0); //quit? + } + console_changed = 1; +} + +/** + * Read character from console without echo + * @return char readed + */ +int read_console(){ + int rb; + struct timeval tv; + int retval; + fd_set rfds; + FD_ZERO(&rfds); + FD_SET(STDIN_FILENO, &rfds); + tv.tv_sec = 0; tv.tv_usec = 10000; + retval = select(1, &rfds, NULL, NULL, &tv); + if(!retval) rb = 0; + else { + if(FD_ISSET(STDIN_FILENO, &rfds)) rb = getchar(); + else rb = 0; + } + return rb; +} + +/** + * getchar() without echo + * wait until at least one character pressed + * @return character readed + */ +int mygetchar(){ // getchar() without need of pressing ENTER + int ret; + do ret = read_console(); + while(ret == 0); + return ret; +} + + +/******************************************************************************\ + * TTY with select() +\******************************************************************************/ +static struct termio oldtty, tty; // TTY flags +static int comfd = -1; // TTY fd + +// run on exit: +void restore_tty(){ + if(comfd == -1) return; + ioctl(comfd, TCSANOW, &oldtty ); // return TTY to previous state + close(comfd); + comfd = -1; +} + +#ifndef BAUD_RATE +#define BAUD_RATE B9600 +#endif +// init: +void tty_init(char *comdev){ + DBG("\nOpen port...\n"); + if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){ + WARN("Can't use port %s\n",comdev); + ioctl(comfd, TCSANOW, &oldtty); // return TTY to previous state + close(comfd); + signals(0); // quit? + } + DBG(" OK\nGet current settings... "); + if(ioctl(comfd,TCGETA,&oldtty) < 0){ // Get settings + /// "Не могу получить настройки" + WARN(_("Can't get settings")); + signals(0); + } + tty = oldtty; + tty.c_lflag = 0; // ~(ICANON | ECHO | ECHOE | ISIG) + tty.c_oflag = 0; + tty.c_cflag = BAUD_RATE|CS8|CREAD|CLOCAL; // 9.6k, 8N1, RW, ignore line ctrl + tty.c_cc[VMIN] = 0; // non-canonical mode + tty.c_cc[VTIME] = 5; + if(ioctl(comfd,TCSETA,&tty) < 0){ + /// "Не могу установить настройки" + WARN(_("Can't set settings")); + signals(0); + } + DBG(" OK\n"); +} +/** + * Read data from TTY + * @param buff (o) - buffer for data read + * @param length - buffer len + * @return amount of readed bytes + */ +size_t read_tty(uint8_t *buff, size_t length){ + ssize_t L = 0; + fd_set rfds; + struct timeval tv; + int retval; + FD_ZERO(&rfds); + FD_SET(comfd, &rfds); + tv.tv_sec = 0; tv.tv_usec = 50000; // wait for 50ms + retval = select(comfd + 1, &rfds, NULL, NULL, &tv); + if (!retval) return 0; + if(FD_ISSET(comfd, &rfds)){ + if((L = read(comfd, buff, length)) < 1) return 0; + } + return (size_t)L; +} + +int write_tty(uint8_t *buff, size_t length){ + ssize_t L = write(comfd, buff, length); + if((size_t)L != length){ + /// "Ошибка записи!" + WARN("Write error!"); + return 1; + } + return 0; +} + + +/** + * Safely convert data from string to double + * + * @param num (o) - double number read from string + * @param str (i) - input string + * @return 1 if success, 0 if fails + */ +int str2double(double *num, const char *str){ + double res; + char *endptr; + if(!str) return 0; + res = strtod(str, &endptr); + if(endptr == str || *str == '\0' || *endptr != '\0'){ + /// "Неправильный формат числа double!" + WARNX("Wrong double number format!"); + return 0; + } + if(num) *num = res; // you may run it like myatod(NULL, str) to test wether str is double number + return 1; +} diff --git a/usefull_macros.h b/usefull_macros.h new file mode 100644 index 0000000..3473e84 --- /dev/null +++ b/usefull_macros.h @@ -0,0 +1,134 @@ +/* + * usefull_macros.h - a set of usefull macros: memory, color etc + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __USEFULL_MACROS_H__ +#define __USEFULL_MACROS_H__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#if defined GETTEXT_PACKAGE && defined LOCALEDIR +/* + * GETTEXT + */ +#include +#define _(String) gettext(String) +#define gettext_noop(String) String +#define N_(String) gettext_noop(String) +#else +#define _(String) (String) +#define N_(String) (String) +#endif +#include +#include +#include +#include +#include +#include + + +// unused arguments with -Wall -Werror +#define _U_ __attribute__((__unused__)) + +/* + * Coloured messages output + */ +#define RED "\033[1;31;40m" +#define GREEN "\033[1;32;40m" +#define OLDCOLOR "\033[0;0;0m" + +#ifndef FALSE +#define FALSE (0) +#endif + +#ifndef TRUE +#define TRUE (1) +#endif + +/* + * ERROR/WARNING messages + */ +extern int globErr; +extern void signals(int sig); +#define ERR(...) do{globErr=errno; _WARN(__VA_ARGS__); signals(9);}while(0) +#define ERRX(...) do{globErr=0; _WARN(__VA_ARGS__); signals(9);}while(0) +#define WARN(...) do{globErr=errno; _WARN(__VA_ARGS__);}while(0) +#define WARNX(...) do{globErr=0; _WARN(__VA_ARGS__);}while(0) + +/* + * print function name, debug messages + * debug mode, -DEBUG + */ +#ifdef EBUG + #define FNAME() fprintf(stderr, "\n%s (%s, line %d)\n", __func__, __FILE__, __LINE__) + #define DBG(...) do{fprintf(stderr, "%s (%s, line %d): ", __func__, __FILE__, __LINE__); \ + fprintf(stderr, __VA_ARGS__); \ + fprintf(stderr, "\n");} while(0) +#else + #define FNAME() do{}while(0) + #define DBG(...) do{}while(0) +#endif //EBUG + +/* + * Memory allocation + */ +#define ALLOC(type, var, size) type * var = ((type *)my_alloc(size, sizeof(type))) +#define MALLOC(type, size) ((type *)my_alloc(size, sizeof(type))) +#define FREE(ptr) do{if(ptr){free(ptr); ptr = NULL;}}while(0) + +double dtime(); + +// functions for color output in tty & no-color in pipes +extern int (*red)(const char *fmt, ...); +extern int (*_WARN)(const char *fmt, ...); +extern int (*green)(const char *fmt, ...); +void * my_alloc(size_t N, size_t S); +void initial_setup(); + +// mmap file +typedef struct{ + char *data; + size_t len; +} mmapbuf; +mmapbuf *My_mmap(char *filename); +void My_munmap(mmapbuf *b); + +void restore_console(); +void setup_con(); +int read_console(); +int mygetchar(); + +void restore_tty(); +void tty_init(char *comdev); +size_t read_tty(uint8_t *buff, size_t length); +int write_tty(uint8_t *buff, size_t length); + +int str2double(double *num, const char *str); + +#endif // __USEFULL_MACROS_H__ diff --git a/zern_private.h b/zern_private.h new file mode 100644 index 0000000..626c944 --- /dev/null +++ b/zern_private.h @@ -0,0 +1,46 @@ +/* + * zern_private.h - private variables for zernike.c & zernikeR.c + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __ZERN_PRIVATE_H__ +#define __ZERN_PRIVATE_H__ + +#include +#include +#include +#include + +#ifndef iabs +#define iabs(a) (((a)<(0)) ? (-a) : (a)) +#endif + +extern double *FK; + +// zernike.c +void build_factorial(); +void free_rpow(double ***Rpow, int n); +void build_rpow(int W, int H, int n, double **Rad, double ***Rad_pow); +double **build_rpowR(int n, int Sz, polar *P); + +// zernike_annular.c +polar *conv_r(polar *r0, int Sz); + +#endif // __ZERN_PRIVATE_H__ diff --git a/zernike.c b/zernike.c new file mode 100644 index 0000000..a6302b9 --- /dev/null +++ b/zernike.c @@ -0,0 +1,494 @@ +/* + * zernike.c - wavefront decomposition using Zernike polynomials + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +/* + * Literature: +@ARTICLE{2007OExpr..1518014Z, + author = {{Zhao}, C. and {Burge}, J.~H.}, + title = "{Orthonormal vector polynomials in a unit circle Part I: basis set derived from gradients of Zernike polynomials}", + journal = {Optics Express}, + year = 2007, + volume = 15, + pages = {18014}, + doi = {10.1364/OE.15.018014}, + adsurl = {http://adsabs.harvard.edu/abs/2007OExpr..1518014Z}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} +@ARTICLE{2008OExpr..16.6586Z, + author = {{Zhao}, C. and {Burge}, J.~H.}, + title = "{Orthonormal vector polynomials in a unit circle, Part II : completing the basis set}", + journal = {Optics Express}, + year = 2008, + month = apr, + volume = 16, + pages = {6586}, + doi = {10.1364/OE.16.006586}, + adsurl = {http://adsabs.harvard.edu/abs/2008OExpr..16.6586Z}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + * + * !!!ATTENTION!!! Axe Y points upside-down! + */ + + +#include "zernike.h" +#include "zern_private.h" +#include "usefull_macros.h" + +double *FK = NULL; +/** + * Build pre-computed array of factorials from 1 to 100 + */ +void build_factorial(){ + double F = 1.; + int i; + if(FK) return; + FK = MALLOC(double, 100); + FK[0] = 1.; + for(i = 1; i < 100; i++) + FK[i] = (F *= (double)i); +} + +double Z_prec = 1e-6; // precision of Zernike coefficients + + +/** + * 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; +} + +/** + * Free array of R powers with power n + * @param Rpow (i) - array to free + * @param n - power of Zernike polinomial for that array (array size = n+1) + */ +void free_rpow(double ***Rpow, int n){ + int i, N = n+1; + for(i = 0; i < N; i++) FREE((*Rpow)[i]); + FREE(*Rpow); +} + +/** + * Build two arrays: with R and its powers (from 0 to n inclusive) + * for cartesian quoter I of matrix with size WxH + * @param W, H - size of initial matrix + * @param n - power of Zernike polinomial (array size = n+1) + * @param Rad (o) - NULL or array with R in quater I + * @param Rad_pow (o) - NULL or array with powers of R + */ +void build_rpow(int W, int H, int n, double **Rad, double ***Rad_pow){ + double Rnorm = fmax((W - 1.) / 2., (H - 1.) / 2.); + int i,j, k, 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 + } + } + if(Rad) *Rad = R; + else FREE(R); + if(Rad_pow) *Rad_pow = Rpow; + else free_rpow(&Rpow, n); +} + +/** + * Calculate value of Zernike polynomial on rectangular 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 array of Zernike polynomials on given matrix + */ +double *zernfun(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.; // coordinate of circle's middle + int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, w = (W+1)/2; + double *R, **Rpow; + build_rpow(W, H, n, &R, &Rpow); + int SS = W * H; + double ZSum = 0.; + // now fill output matrix + Zarr = MALLOC(double, SS); // output matrix W*H pixels + 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(2.*(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); + free_rpow(&Rpow, n); + return Zarr; +} + +/** + * Zernike polynomials in Noll notation for square matrix + * @param p - index of polynomial + * @other params are like in zernfun + * @return zernfun + */ +double *zernfunN(int p, int W, int H, double *norm){ + int n, m; + convert_Zidx(p, &n, &m); + return zernfun(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 + 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 = zernfunN(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; +} + +/** + * 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 = zernfunN(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; +} + + +/** + * Components of Zj gradient without constant factor + * all parameters are as in zernfun + * @return array of gradient's components + */ +point *gradZ(int n, int m, int W, int H, double *norm){ + point *gZ = NULL; + bool erparm = false; + if(W < 2 || H < 2) + errx(1, "Sizes of matrix must be > 2!"); + if(n > 100) + errx(1, "Order of gradient of Zernike polynomial must be <= 100!"); + if(n < 1) 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 gradient of Zernike polynomial (%d, %d)", n, m); + if(!FK) build_factorial(); + double Xc = (W - 1.) / 2., Yc = (H - 1.) / 2.; // coordinate of circle's middle + int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, isum = (n+m_abs)/2, w = (W+1)/2; + double *R, **Rpow; + build_rpow(W, H, n, &R, &Rpow); + int SS = W * H; + // now fill output matrix + gZ = MALLOC(point, SS); + double ZSum = 0.; + for(j = 0; j < H; j++){ + point *Zptr = &gZ[j*W]; + double Ryd = fabs(j - Yc); + int Ry = w * (int)Ryd; // Y coordinate on R matrix + for(i = 0; i < W; i++, Zptr++){ + double Rxd = fabs(i - Xc); + int Ridx = Ry + (int)Rxd; // coordinate on R matrix + double Rcurr = R[Ridx]; + if(Rcurr > 1. || fabs(Rcurr) < DBL_EPSILON) continue; // throw out points with R>1 + double theta = atan2(j - Yc, i - Xc); + // components of grad Zj: + + // 1. Theta_j + double Tj = 1., costm, sintm; + sincos(theta*(double)m_abs, &sintm, &costm); + if(m){ + if(m > 0) + Tj = costm; + else + Tj = sintm; + } + + // 2. dTheta_j/Dtheta + double dTj = 0.; + if(m){ + if(m < 0) + dTj = m_abs * costm; + else + dTj = -m_abs * sintm; + } + // 3. R_j & dR_j/dr + double Rj = 0., dRj = 0.; + for(k = 0; k <= iup; k++){ + double rj = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//----------------------------------- ----- ------------------------------- + FK[k]*FK[isum-k]*FK[iup-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)! + ); +//DBG("rj = %g (n=%d, k=%d)\n",rj, n, k); + int kidx = n - 2*k; + Rj += rj * Rpow[kidx][Ridx]; // *R^{n-2k} + if(kidx > 0) + dRj += rj * kidx * Rpow[kidx - 1][Ridx]; + } + // normalisation for Zernike + double eps_m = (m) ? 1. : 2., sq = sqrt(2.*(double)(n+1) / M_PI / eps_m); + Rj *= sq, dRj *= sq; + // 4. sin/cos + double sint, cost; + sincos(theta, &sint, &cost); + + // projections of gradZj + double TdR = Tj * dRj, RdT = Rj * dTj / Rcurr; + Zptr->x = TdR * cost - RdT * sint; + Zptr->y = TdR * sint + RdT * cost; + // norm factor + ZSum += Zptr->x * Zptr->x + Zptr->y * Zptr->y; + } + } + if(norm) *norm = ZSum; + // free unneeded memory + FREE(R); + free_rpow(&Rpow, n); + return gZ; +} + +/** + * Build components of vector polynomial Sj + * all parameters are as in zernfunN + * @return Sj(n,m) on image points + */ +point *zerngrad(int p, int W, int H, double *norm){ + int n, m; + convert_Zidx(p, &n, &m); + if(n < 1) errx(1, "Order of gradient Z must be > 0!"); + int m_abs = iabs(m); + int i,j; + double K = 1., K1 = 1.;///sqrt(2.*n*(n+1.)); + point *Sj = MALLOC(point, W*H); + point *Zj = gradZ(n, m, W, H, NULL); + double Zsum = 0.; + if(n == m_abs || n < 3){ // trivial case: n = |m| (in case of n=2,m=0 n'=0 -> grad(0,0)=0 + for(j = 0; j < H; j++){ + int idx = j*W; + point *Sptr = &Sj[idx], *Zptr = &Zj[idx]; + for(i = 0; i < W; i++, Sptr++, Zptr++){ + Sptr->x = K1*Zptr->x; + Sptr->y = K1*Zptr->y; + Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y; + } + } + }else{ + K = sqrt(((double)n+1.)/(n-1.)); + //K1 /= sqrt(2.); + // n != |m| + // I run gradZ() twice! But another variant (to make array of Zj) can meet memory lack + point *Zj_= gradZ(n-2, m, W, H, NULL); + for(j = 0; j < H; j++){ + int idx = j*W; + point *Sptr = &Sj[idx], *Zptr = &Zj[idx], *Z_ptr = &Zj_[idx]; + for(i = 0; i < W; i++, Sptr++, Zptr++, Z_ptr++){ + Sptr->x = K1*(Zptr->x - K * Z_ptr->x); + Sptr->y = K1*(Zptr->y - K * Z_ptr->y); + Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y; + } + } + FREE(Zj_); + } + FREE(Zj); + if(norm) *norm = Zsum; + return Sj; +} + +/** + * Decomposition of image with normals to wavefront by Zhao's polynomials + * all like Zdecompose + * @return array of coefficients + */ +double *gradZdecompose(int Nmax, int W, int H, point *image, int *Zsz, int *lastIdx){ + int i, SS = W*H, pmax, maxIdx = 0; + double *Zidxs = NULL; + point *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(point, SS); + memcpy(icopy, image, SS*sizeof(point)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 1; i < pmax; i++){ // now we fill array + double norm; + point *dZcoeff = zerngrad(i, W, H, &norm); + int j; + point *iptr = icopy, *zptr = dZcoeff; + double K = 0.; + for(j = 0; j < SS; j++, iptr++, zptr++) + K += (zptr->x*iptr->x + zptr->y*iptr->y) / norm; // multiply matrixes to get coefficient + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + iptr = icopy; zptr = dZcoeff; + for(j = 0; j < SS; j++, iptr++, zptr++){ + iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values + iptr->y -= K * zptr->y; + } + FREE(dZcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Restoration of image with normals Zhao's polynomials coefficients + * all like Zcompose + * @return restored image + */ +point *gradZcompose(int Zsz, double *Zidxs, int W, int H){ + int i, SS = W*H; + point *image = MALLOC(point, SS); + for(i = 1; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + point *Zcoeff = zerngrad(i, W, H, NULL); + int j; + point *iptr = image, *zptr = Zcoeff; + for(j = 0; j < SS; j++, iptr++, zptr++){ + iptr->x += K * zptr->x; + iptr->y += K * zptr->y; + } + FREE(Zcoeff); + } + return image; +} + +double *convGradIdxs(double *gradIdxs, int Zsz){ + double *idxNew = MALLOC(double, Zsz); + int i; + for(i = 1; i < Zsz; i++){ + int n,m; + convert_Zidx(i, &n, &m); + int j = ((n+2)*(n+4) + m) / 2; + if(j >= Zsz) j = 0; + idxNew[i] = (gradIdxs[i] - sqrt((n+3.)/(n+1.))*gradIdxs[j]); + } + return idxNew; +} diff --git a/zernike.h b/zernike.h new file mode 100644 index 0000000..3463061 --- /dev/null +++ b/zernike.h @@ -0,0 +1,92 @@ +/* + * zernike.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __ZERNIKE_H__ +#define __ZERNIKE_H__ + +#include +#include + +// focal ratio of BTA in millimeters +#define FOCAL_R (24024.) +// BTA primary mirror radius (in mm) +#define MIR_R (3025.) +// Distance from mirror to hartmann mask +#define HARTMANN_Z (20017.) + +/*************** Data structures & typedefs ***************/ +// point coordinates +typedef struct{ + double x,y; +} point; + +typedef struct{ + double r,theta; +} polar; + +// 2D array +typedef struct{ + double **data; + size_t len; // size of 1D arrays + size_t num; // number of 1D arrays +}_2D; + +extern double Z_prec; // precision of Zernike coefficients + +#ifndef DBL_EPSILON +#define DBL_EPSILON 2.2204460492503131e-16 +#endif + +/*************** Base functions ***************/ +void convert_Zidx(int p, int *N, int *M); +extern double Z_prec; + +/*************** Zernike on rectangular equidistant coordinate matrix ***************/ +double *zernfun(int n, int m, int W, int H, double *norm); +double *zernfunN(int p, int W, int H, double *norm); + +double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx); +double *Zcompose(int Zsz, double *Zidxs, int W, int H); + +double *gradZdecompose(int Nmax, int W, int H, point *image, int *Zsz, int *lastIdx); +point *gradZcompose(int Zsz, double *Zidxs, int W, int H); +double *convGradIdxs(double *gradIdxs, int Zsz); + +/*************** Zernike on a points set ***************/ +double *zernfunR(int n, int m, int Sz, polar *P, double *norm); +double *zernfunNR(int p, int Sz, polar *P, double *norm); +double *ZdecomposeR(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); +double *ZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P); +double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); +double *QR_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); +double *gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx); +double *LS_gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx); +point *gradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P); +point *directGradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P); +double *directGradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx); + +/*************** Annular Zernike ***************/ +_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm); +double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P); +double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); + +#endif // __ZERNIKE_H__ diff --git a/zernikeR.c b/zernikeR.c new file mode 100644 index 0000000..12adfcd --- /dev/null +++ b/zernikeR.c @@ -0,0 +1,687 @@ +/* + * zernikeR.c - Zernike decomposition for scattered points & annular aperture + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include + +#include "usefull_macros.h" +#include "zernike.h" +#include "zern_private.h" + +/** + * Build array with R powers (from 0 to n inclusive) + * @param n - power of Zernike polinomial (array size = n+1) + * @param Sz - size of P array + * @param P (i) - polar coordinates of points + */ +double **build_rpowR(int n, int Sz, polar *P){ + int i, j, N = n + 1; + double **Rpow = MALLOC(double*, N); + Rpow[0] = MALLOC(double, Sz); + for(i = 0; i < Sz; i++) Rpow[0][i] = 1.; // zero's power + for(i = 1; i < N; i++){ // Rpow - is quater I of cartesian coordinates ('cause other are fully simmetrical) + Rpow[i] = MALLOC(double, Sz); + double *rp = Rpow[i], *rpo = Rpow[i-1]; + polar *p = P; + for(j = 0; j < Sz; j++, rp++, rpo++, p++){ + *rp = (*rpo) * p->r; + } + } + return Rpow; +} + +bool check_parameters(int n, int m, int Sz, polar *P){ + bool erparm = false; + if(Sz < 3 || !P) + ERRX(_("Size of matrix must be > 2!")); + if(n > 100) + ERRX(_("Order of Zernike polynomial must be <= 100!")); + if(n < 0) erparm = true; + if(n < iabs(m)) erparm = true; // |m| must be <= n + if((n - m) % 2) erparm = true; // n-m must differ by a prod of 2 + if(erparm) + ERRX(_("Wrong parameters of Zernike polynomial (%d, %d)"), n, m); + else + if(!FK) build_factorial(); + return erparm; +} + +/** + * Zernike function for scattering data + * @param n,m - orders of polynomial + * @param Sz - number of points + * @param P(i) - array with points coordinates (polar, r<=1) + * @param norm(o) - (optional) norm coefficient + * @return dynamically allocated array with Z(n,m) for given array P + */ +double *zernfunR(int n, int m, int Sz, polar *P, double *norm){ + if(check_parameters(n, m, Sz, P)) return NULL; + int j, k, m_abs = iabs(m), iup = (n-m_abs)/2; + double **Rpow = build_rpowR(n, Sz, P); + double ZSum = 0.; + // now fill output matrix + double *Zarr = MALLOC(double, Sz); // output matrix + double *Zptr = Zarr; + polar *p = P; + for(j = 0; j < Sz; j++, p++, Zptr++){ + double Z = 0.; + if(p->r > 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][j]; // *R^{n-2k} + } + // normalize + double eps_m = (m) ? 1. : 2.; + Z *= sqrt(2.*(n+1.) / M_PI / eps_m ); + double m_theta = (double)m_abs * p->theta; + // multiply to angular function: + if(m){ + if(m > 0) + Z *= cos(m_theta); + else + Z *= sin(m_theta); + } + *Zptr = Z; + ZSum += Z*Z; + } + if(norm) *norm = ZSum; + // free unneeded memory + free_rpow(&Rpow, n); + return Zarr; +} + +/** + * Zernike polynomials in Noll notation for square matrix + * @param p - index of polynomial + * @other params are like in zernfunR + * @return zernfunR + */ +double *zernfunNR(int p, int Sz, polar *P, double *norm){ + int n, m; + convert_Zidx(p, &n, &m); + return zernfunR(n,m,Sz,P,norm); +} + +/** + * Zernike decomposition of WF with coordinates P + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param Sz, P(i) - size (Sz) of points array (P) + * @param heights(i) - wavefront walues in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of Zernike coefficients + */ +double *ZdecomposeR(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + double *Zidxs = NULL, *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(double, Sz); + memcpy(icopy, heights, Sz*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, Sz, P, &norm); + int j; + double *iptr = icopy, *zptr = Zcoeff, K = 0.; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += (*zptr) * (*iptr) / norm; // multiply matrixes to get coefficient + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; +// Without these 3 lines wavefront decomposition on scattered and/or small amount of points would be wrong ===================> + iptr = icopy; zptr = Zcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++) + *iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values +//<=================== */ + FREE(Zcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Restoration of image in points P by Zernike polynomials coefficients + * @param Zsz (i) - number of actual elements in coefficients array + * @param Zidxs(i) - array with Zernike coefficients + * @param Sz, P(i) - number (Sz) of points (P) + * @return restored image + */ +double *ZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){ + int i; + double *image = MALLOC(double, Sz); + for(i = 0; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + double *Zcoeff = zernfunNR(i, Sz, P, NULL); + int j; + double *iptr = image, *zptr = Zcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++) + *iptr += K * (*zptr); // add next Zernike polynomial + FREE(Zcoeff); + } + return image; +} + +/** + * Prints out GSL matrix + * @param M (i) - matrix to print + */ +void print_matrix(gsl_matrix *M){ + int x,y, H = M->size1, W = M->size2; + size_t T = M->tda; + printf("\nMatrix %dx%d:\n", W, H); + for(y = 0; y < H; y++){ + double *ptr = &(M->data[y*T]); + printf("str %6d", y); + for(x = 0; x < W; x++, ptr++) + printf("%6.1f ", *ptr); + printf("\n"); + } + printf("\n\n"); +} + +/** +To try less-squares fit I decide after reading of +@ARTICLE{2010ApOpt..49.6489M, + author = {{Mahajan}, V.~N. and {Aftab}, M.}, + title = "{Systematic comparison of the use of annular and Zernike circle polynomials for annular wavefronts}", + journal = {\ao}, + year = 2010, + month = nov, + volume = 49, + pages = {6489}, + doi = {10.1364/AO.49.006489}, + adsurl = {http://adsabs.harvard.edu/abs/2010ApOpt..49.6489M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} +*/ +/* + * n'th column of array m is polynomial of n'th degree, + * m'th row is m'th data point + * + * We fill matrix with polynomials by known datapoints coordinates, + * after that by less-square fit we get coefficients of decomposition + */ +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 + if(Zsz) *Zsz = pmax; +/* + (from GSL) + 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); + 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; +} + +/** + * Pseudo-annular Zernike polynomials + * They are just a linear composition of Zernike, made by Gram-Schmidt ortogonalisation + * Real Zernike koefficients restored by reverce transformation matrix + * + * Suppose we have a wavefront data in set of scattered points ${(x,y)}$, we want to + * find Zernike coefficients $z$ such that product of Zernike polynomials, $Z$, and + * $z$ give us that wavefront data with very little error (depending on $Z$ degree). + * + * Of cource, $Z$ isn't orthonormal on our basis, so we need to create some intermediate + * polynomials, $U$, which will be linear dependent from $Z$ (and this dependency + * should be strict and reversable, otherwise we wouldn't be able to reconstruct $z$ + * from $u$): $U = Zk$. So, we have: $W = Uu + \epsilon$ and $W = Zz + \epsilon$. + * + * $U$ is orthonormal, so $U^T = U^{-1}$ and (unlike to $Z$) this reverce matrix + * exists and is unique. This mean that $u = W^T U = U^T W$. + * Coefficients matrix, $k$ must be inversable, so $Uk^{-1} = Z$, this mean that + * $z = uk^{-1}$. + * + * Our main goal is to find this matrix, $k$. + * + * 1. QR-decomposition + * R. Navarro and J. Arines. Complete Modal Representation with Discrete Zernike + * Polynomials - Critical Sampling in Non Redundant Grids. INTECH Open Access + * Publisher, 2011. + * + * In this case non-orthogonal matrix $Z$ decomposed to orthogonal matrix $Q$ and + * right-triangular matrix $R$. In our case $Q$ is $U$ itself and $R$ is $k^{-1}$. + * QR-decomposition gives us an easy way to compute coefficient's matrix, $k$. + * Polynomials in $Q$ are linear-dependent from $Z$, each $n^{th}$ polynomial in $Q$ + * depends from first $n$ polynomials in $Z$. So, columns of $R$ are coefficients + * for making $U$ from $Z$; rows in $R$ are coefficients for making $z$ from $u$. + * + * 2. Cholesky decomposition + * In this case for any non-orthogonal matrix with real values we have: + * $A^{T}A = LL^{T}$, where $L$ is left-triangular matrix. + * Orthogonal basis made by equation $U = A(L^{-1})^T$. And, as $A = UL^T$, we + * can reconstruct coefficients matrix: $k = (L^{-1})^T$. + */ + +double *QR_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 + if(Sz < pmax) ERRX(_("Number of points must be >= number of polynomials!")); + if(Zsz) *Zsz = pmax; + int k, x,y; + + //make new polar + polar *Pn = conv_r(P, Sz); + // 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); + // Q-matrix (its first pmax columns is our basis) + gsl_matrix *Q = gsl_matrix_calloc(Sz, Sz); + // R-matrix (coefficients) + gsl_matrix *R = gsl_matrix_calloc(Sz, pmax); + // fill matrix with coefficients + size_t T = M->tda; + for(x = 0; x < pmax; x++){ + double norm, *Zcoeff = zernfunNR(x, Sz, Pn, &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 *tau = gsl_vector_calloc(pmax); // min(size(M)) + gsl_linalg_QR_decomp(M, tau); + + gsl_linalg_QR_unpack(M, tau, Q, R); +//print_matrix(R); + gsl_matrix_free(M); + gsl_vector_free(tau); + + double *Zidxs = MALLOC(double, pmax); + T = Q->tda; + size_t Tr = R->tda; + + for(k = 0; k < pmax; k++){ // cycle by powers + double sumD = 0.; + double *Qptr = &(Q->data[k]); + for(y = 0; y < Sz; y++, Qptr+=T){ // cycle by points + sumD += heights[y] * (*Qptr); + } + Zidxs[k] = sumD; + } + gsl_matrix_free(Q); + + // now restore Zernike coefficients + double *Zidxs_corr = MALLOC(double, pmax); + int maxIdx = 0; + for(k = 0; k < pmax; k++){ + double sum = 0.; + double *Rptr = &(R->data[k*(Tr+1)]), *Zptr = &(Zidxs[k]); + for(x = k; x < pmax; x++, Zptr++, Rptr++){ + sum += (*Zptr) * (*Rptr); + } + if(fabs(sum) < Z_prec) continue; + Zidxs_corr[k] = sum; + maxIdx = k; + } + gsl_matrix_free(R); + FREE(Zidxs); + FREE(Pn); + if(lastIdx) *lastIdx = maxIdx; + return Zidxs_corr; +} + + +/** + * Components of Zj gradient without constant factor + * @param n,m - orders of polynomial + * @param Sz - number of points + * @param P (i) - coordinates of reference points + * @param norm (o) - norm factor or NULL + * @return array of gradient's components + */ +point *gradZR(int n, int m, int Sz, polar *P, double *norm){ + if(check_parameters(n, m, Sz, P)) return NULL; + point *gZ = NULL; + int j, k, m_abs = iabs(m), iup = (n-m_abs)/2, isum = (n+m_abs)/2; + double **Rpow = build_rpowR(n, Sz, P); + // now fill output matrix + gZ = MALLOC(point, Sz); + point *Zptr = gZ; + double ZSum = 0.; + polar *p = P; + for(j = 0; j < Sz; j++, p++, Zptr++){ + if(p->r > 1.) continue; // throw out points with R>1 + double theta = p->theta; + // components of grad Zj: + + // 1. Theta_j; 2. dTheta_j/Dtheta + double Tj = 1., dTj = 0.; + if(m){ + double costm, sintm; + sincos(theta*(double)m_abs, &sintm, &costm); + if(m > 0){ + Tj = costm; + dTj = -m_abs * sintm; + }else{ + Tj = sintm; + dTj = m_abs * costm; + } + } + + // 3. R_j & dR_j/dr + double Rj = 0., dRj = 0.; + for(k = 0; k <= iup; k++){ + double rj = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//----------------------------------- ----- ------------------------------- + FK[k]*FK[isum-k]*FK[iup-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)! + ); +//DBG("rj = %g (n=%d, k=%d)\n",rj, n, k); + int kidx = n - 2*k; + Rj += rj * Rpow[kidx][j]; // *R^{n-2k} + if(kidx > 0) + dRj += rj * kidx * Rpow[kidx - 1][j]; // *(n-2k)*R^{n-2k-1} + } + // normalisation for Zernike + double eps_m = (m) ? 1. : 2., sq = sqrt(2.*(double)(n+1) / M_PI / eps_m); + Rj *= sq, dRj *= sq; + // 4. sin/cos + double sint, cost; + sincos(theta, &sint, &cost); + + // projections of gradZj + double TdR = Tj * dRj, RdT = Rj * dTj / p->r; + Zptr->x = TdR * cost - RdT * sint; + Zptr->y = TdR * sint + RdT * cost; + // norm factor + ZSum += Zptr->x * Zptr->x + Zptr->y * Zptr->y; + } + if(norm) *norm = ZSum; + // free unneeded memory + free_rpow(&Rpow, n); + return gZ; +} + +/** + * Build components of vector polynomial Sj + * @param p - index of polynomial + * @param Sz - number of points + * @param P (i) - coordinates of reference points + * @param norm (o) - norm factor or NULL + * @return Sj(n,m) on image points + */ +point *zerngradR(int p, int Sz, polar *P, double *norm){ + int n, m, i; + convert_Zidx(p, &n, &m); + if(n < 1) ERRX(_("Order of gradient Z must be > 0!")); + int m_abs = iabs(m); + double Zsum, K = 1.; + point *Sj = gradZR(n, m, Sz, P, &Zsum); + if(n != m_abs && n > 2){ // avoid trivial case: n = |m| (in case of n=2,m=0 n'=0 -> grad(0,0)=0 + K = sqrt(((double)n+1.)/(n-1.)); + Zsum = 0.; + point *Zj= gradZR(n-2, m, Sz, P, NULL); + point *Sptr = Sj, *Zptr = Zj; + for(i = 0; i < Sz; i++, Sptr++, Zptr++){ + Sptr->x -= K * Zptr->x; + Sptr->y -= K * Zptr->y; + Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y; + } + FREE(Zj); + } + if(norm) *norm = Zsum; + return Sj; +} + +/** + * Decomposition of image with normals to wavefront by Zhao's 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) + * @param grads(i) - wavefront gradients values in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of coefficients + */ +double *LS_gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx){ + int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + if(Zsz) *Zsz = pmax; + // Now allocate matrix: its Nth row is equation for Nth data point, + // Mth column is Z_M + int Sz2 = Sz*2, x, y; + gsl_matrix *M = gsl_matrix_calloc(Sz2, pmax); + // fill matrix with coefficients + size_t T = M->tda, S = T * Sz; // T is matrix stride, S - index of second coefficient + for(x = 1; x < pmax; x++){ + double norm; + int n, m; + convert_Zidx(x, &n, &m); + point *dZcoeff = gradZR(n, m, Sz, P, &norm), *dZptr = dZcoeff; + //point *dZcoeff = zerngradR(x, Sz, P, &norm), *dZptr = dZcoeff; + double *ptr = &(M->data[x]); + // X-component is top part of resulting matrix, Y is bottom part + for(y = 0; y < Sz; y++, ptr+=T, dZptr++){ // fill xth polynomial + *ptr = dZptr->x; + ptr[S] = dZptr->y; + } + FREE(dZcoeff); + } + + gsl_vector *ans = gsl_vector_calloc(pmax); + gsl_vector *b = gsl_vector_calloc(Sz2); + double *ptr = b->data; + for(x = 0; x < Sz; x++, ptr++, grads++){ + // fill components of vector b like components of matrix M + *ptr = grads->x; + ptr[Sz] = grads->y; + } + + gsl_vector *tau = gsl_vector_calloc(pmax); + gsl_linalg_QR_decomp(M, tau); + + gsl_vector *residual = gsl_vector_calloc(Sz2); + gsl_linalg_QR_lssolve(M, tau, b, ans, residual); + ptr = &ans->data[1]; + int maxIdx = 0; + double *Zidxs = MALLOC(double, pmax); + for(x = 1; 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(b); + gsl_vector_free(tau); + gsl_vector_free(residual); + + if(lastIdx) *lastIdx = maxIdx; + return Zidxs; +} + +/** + * Decomposition of image with normals to wavefront by Zhao's polynomials + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param Sz, P(i) - size (Sz) of points array (P) + * @param grads(i) - wavefront gradients values in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of coefficients + */ +double *gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + double *Zidxs = MALLOC(double, pmax); + point *icopy = MALLOC(point, Sz); + memcpy(icopy, grads, Sz*sizeof(point)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 1; i < pmax; i++){ // now we fill array + double norm; + point *dZcoeff = zerngradR(i, Sz, P, &norm); + int j; + point *iptr = icopy, *zptr = dZcoeff; + double K = 0.; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += zptr->x*iptr->x + zptr->y*iptr->y; // multiply matrixes to get coefficient + K /= norm; + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + iptr = icopy; zptr = dZcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values + iptr->y -= K * zptr->y; + } + FREE(dZcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Restoration of wavefront normals by given Zhao's polynomials coefficients + * all like Zcompose, but `Zidxs` are Zhao's (not direct Zernike)! + * @return restored image + */ +point *gradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){ + int i; + point *image = MALLOC(point, Sz); + for(i = 1; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + point *Zcoeff = zerngradR(i, Sz, P, NULL); + int j; + point *iptr = image, *zptr = Zcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + iptr->x += K * zptr->x; + iptr->y += K * zptr->y; + } + FREE(Zcoeff); + } + return image; +} + +/** + * Direct gradients of Zernike polynomials + */ +point *directGradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){ + point *image = MALLOC(point, Sz); + for(int p = 1; p < Zsz; ++p){ + int n, m; + double K = Zidxs[p]; + convert_Zidx(p, &n, &m); + point *Zcoeff = gradZR(n, m, Sz, P, NULL); + point *iptr = image, *zptr = Zcoeff; + for(int j = 0; j < Sz; ++j, ++iptr, ++zptr){ + iptr->x += K * zptr->x; + iptr->y += K * zptr->y; + } + FREE(Zcoeff); + } + return image; +} + +/** + * Direct decomposition of image with normals to wavefront by Zernike gradients + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param Sz, P(i) - size (Sz) of points array (P) + * @param grads(i) - wavefront gradients values in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of coefficients + */ +double *directGradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + double *Zidxs = MALLOC(double, pmax); + point *icopy = MALLOC(point, Sz); + memcpy(icopy, grads, Sz*sizeof(point)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 1; i < pmax; i++){ // now we fill array + double norm; + int n, m; + convert_Zidx(i, &n, &m); + point *dZcoeff = gradZR(n, m, Sz, P, &norm); + int j; + point *iptr = icopy, *zptr = dZcoeff; + double K = 0.; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += zptr->x*iptr->x + zptr->y*iptr->y; + K /= norm; + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + + iptr = icopy; zptr = dZcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values + iptr->y -= K * zptr->y; + } + + FREE(dZcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} diff --git a/zernikeR.h b/zernikeR.h new file mode 100644 index 0000000..e69de29 diff --git a/zernike_annular.c b/zernike_annular.c new file mode 100644 index 0000000..a61f275 --- /dev/null +++ b/zernike_annular.c @@ -0,0 +1,388 @@ +/* + * zernike_annular.c - annular Zernike polynomials + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + + +/* + * These polynomials realized according to formulae in: +@ARTICLE{1981JOSA...71...75M, + author = {{Mahajan}, V.~N.}, + title = "{Zernike annular polynomials for imaging systems with annular pupils.}", + journal = {Journal of the Optical Society of America (1917-1983)}, + keywords = {Astronomical Optics}, + year = 1981, + volume = 71, + pages = {75-85}, + adsurl = {http://adsabs.harvard.edu/abs/1981JOSA...71...75M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@ARTICLE{2006ApOpt..45.3442H, + author = {{Hou}, X. and {Wu}, F. and {Yang}, L. and {Wu}, S. and {Chen}, Q. + }, + title = "{Full-aperture wavefront reconstruction from annular subaperture interferometric data by use of Zernike annular polynomials and a matrix method for testing large aspheric surfaces}", + journal = {\ao}, + keywords = {Mathematics, Optical data processing, Metrology, Optical testing}, + year = 2006, + month = may, + volume = 45, + pages = {3442-3455}, + doi = {10.1364/AO.45.003442}, + adsurl = {http://adsabs.harvard.edu/abs/2006ApOpt..45.3442H}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + */ + +#include "zernike.h" +#include "zern_private.h" +#include "usefull_macros.h" + +// epsilon (==rmin) +double eps = -1., eps2 = -1.; // this global variable changes with each conv_r call! + +/** + * convert coordinates of points on ring into normalized coordinates in [0, 1] + * for further calculations of annular Zernike (expression for R_k^0) + * @param r0 (i) - original coordinates array + * @param Sz - size of r0 + * @return dinamically allocated array with new coordinates, which zero element is + * (0,0) + */ +polar *conv_r(polar *r0, int Sz){ + int x; + polar *Pn = MALLOC(polar, Sz); + double rmax = -1., rmin = 2.; + polar *p = r0, *p1 = Pn; + for(x = 0; x < Sz; x++, p++){ + if(p->r < rmin) rmin = p->r; + if(p->r > rmax) rmax = p->r; + } + if(rmin > rmax || fabs(rmin - rmax) <= DBL_EPSILON || rmax > 1. || rmin < 0.) + ERRX(_("polar coordinates should be in [0,1]!")); + eps = rmin; + eps2 = eps*eps; + p = r0; + //rmax = 1. - eps2; // 1 - eps^2 -- denominator + rmax = rmax*rmax - eps2; + for(x = 0; x < Sz; x++, p++, p1++){ + p1->theta = p->theta; + p1->r = sqrt((p->r*p->r - eps2) / rmax); // r = (r^2 - eps^2) / (1 - eps^2) + } + return Pn; +} + +/** + * Allocate 2D array (H arrays of length W) + * @param len - array's width (number of elements in each 1D array + * @param num - array's height (number of 1D arrays) + * @return dinamically allocated 2D array filled with zeros + */ +_2D *_2ALLOC(int len, int num){ + _2D *r = MALLOC(_2D, 1); + r->data = MALLOC(double*, num); + int x; + for(x = 0; x < num; x++) + r->data[x] = MALLOC(double, len); + r->len = (size_t)len; + r->num = (size_t)num; + return r; +} + +void _2FREE(_2D **a){ + int x, W = (*a)->num; + for(x = 0; x < W; x++) free((*a)->data[x]); + free((*a)->data); + FREE(*a); +} + +/** + * Zernike radial function R_n^0 + * @param n - order of polynomial + * @param Sz - number of points + * @param P (i) - modified polar coordinates: sqrt((r^2-eps^2)/(1-eps^2)) + * @param R (o) - array with size Sz to which polynomial will be stored + */ +void Rn0(int n, int Sz, polar *P, double *R){ + if(Sz < 1 || !P) + ERRX(_("Size of matrix must be > 0!")); + if(n < 0) + ERRX(_("Order of Zernike polynomial must be > 0!")); + if(!FK) build_factorial(); + int j, k, iup = n/2; + double **Rpow = build_rpowR(n, Sz, P); + // now fill output matrix + double *Zptr = R; + polar *p = P; + for(j = 0; j < Sz; j++, p++, Zptr++){ + double Z = 0.; + if(p->r > 1.) continue; // throw out points with R>1 + // calculate R_n^m + for(k = 0; k <= iup; k++){ // Sum + double F = FK[iup-k]; + double z = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//------------------------------------- --------------------- + FK[k]*F*F // k!((n/2-k)!)^2 + ); + Z += z * Rpow[n-2*k][j]; // *R^{n-2k} + } + *Zptr = Z; +//DBG("order %d, point (%g, %g): %g\n", n, p->r, p->theta, Z); + } + // free unneeded memory + free_rpow(&Rpow, n); +} + +/** + * Convert parameters n, m into index in Noll notation, p + * @param n, m - Zernike orders + * @return order in Noll notation + */ +int p_from_nm(int n, int m){ + return (n*n + 2*n + m)/2; +} + +/** + * Create array [p][pt] of annular Zernike polynomials + * @param pmax - max number of polynomial in Noll notation + * @param Sz - size of points array + * @param P (i) - array with points' coordinates + * @param Norm(o)- array [0..pmax] of sum^2 for each polynomial (dynamically allocated) or NULL + * @return + */ +//#undef DBG +//#define DBG(...) +_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm){ + int m, n, j, nmax, mmax, jmax, mW, jm; + convert_Zidx(pmax, &nmax, &mmax); + mmax = nmax; // for a case when p doesn't include a whole n powers + jmax = nmax / 2; + mW = mmax + 1; // width of row + jm = jmax + mW; // height of Q0 and h + polar *Pn = conv_r(P, Sz); + if(eps < 0. || eps2 < 0. || eps > 1. || eps2 > 1.) + ERRX(_("Wrong epsilon value!!!")); + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + _2D *h = _2ALLOC(mW, jm); // constants -> array [(jmax+1)+((mmax+1)-1)][mmax+1] + _2D *Q0= _2ALLOC(mW, jm); // Q(0) -> array [(jmax+1)+((mmax+1)-1)][mmax+1] + _2D *Q = _2ALLOC(Sz, (jmax + 1)*mW); // functions -> array [(jmax+1)*(mmax+1)][Sz] +DBG("Allocated Q&h, size of Q: %zdx%zd, size of H: %zdx%zd\n", Q->num, Q->len, h->num, h->len); +DBG("pmax = %d, nmax = mmax = %d, jmax = %d, jm = %d\n", pmax, nmax, jmax, jm); +DBG("eps=%g, eps2=%g\n", eps, eps2); + // in Q index of [j][m][point] is [j*mmax+m][point] + double **hd = h->data, **Qd = Q->data, **Q0d = Q0->data; + // fill h_j^0 & Q_j^0 + for(j = 0; j < jm; j++){ // main matrix + triangle +DBG("j = %d, m = 0\n", j); + polar Zero = {0., 0.}; + // h_j^0 = (1-eps^2)/(2(2j+1)) + hd[j][0] = (1.-eps2)/(2.*j+1.)/2.; + // Q_j^0(rho, eps) = R_{2j}^0(r) +DBG("Q,h [%d][0], h=%g\n", j, hd[j][0]); + if(j <= jmax) Rn0(2*j, Sz, Pn, Qd[mW*j]); + Rn0(2*j, 1, &Zero, &Q0d[j][0]); // fill also Q(0) + } + FREE(Pn); + // fill h_j^m & Q_j&m + int JMAX = jm-1; // max number of j + for(m = 1; m <= mmax; m++, JMAX--){ // m in outern cycle because of Q_{j+1}! +DBG("\n\nm=%d\n",m); + int mminusone = m - 1; + int idx = mminusone;// index Q_j^{m-1} + for(j = 0; j < JMAX; j++, idx+=mW){ +DBG("\nj=%d\n",j); + // 2(2j+2m-1) * h_j^{m-1} + // H_j^m = ------------------------------ + // (j+m)(1-eps^2) * Q_j^{m-1}(0) + double H = 2*(2*(j+m)-1)/(j+m)/(1-eps2) * hd[j][mminusone]/Q0d[j][mminusone]; + // h_j^m = -H_j^m * Q_{j+1}^{m-1}(0) + hd[j][m] = -H * Q0d[j+1][mminusone]; +DBG("hd[%d][%d] = %g (H = %g)\n", j, m, hd[j][m], H); + // j Q_i^{m-1}(0) * Q_i^{m-1}(r) + // Q_j^m(r) = H_j^m * Sum ----------------------------- + // i=0 h_i^{m-1} + if(j <= jmax){ + int pt; + for(pt = 0; pt < Sz; pt++){ // cycle by points + int i, idx1 = mminusone; + double S = 0; + for(i = 0; i <=j; i++, idx1+=mW){ // cycle Sum + S += Q0d[i][mminusone] * Qd[idx1][pt] / hd[i][mminusone]; + } + Qd[idx+1][pt] = H * S; // Q_j^m(rho^2) +//DBG("Q[%d][%d](%d) = %g\n", j, m, pt, H * S); + } + } + // and Q(0): + double S = 0.; + int i; + for(i = 0; i <=j; i++){ + double qim = Q0d[i][mminusone]; + S += qim*qim / hd[i][mminusone]; + } + Q0d[j][m] = H * S; // Q_j^m(0) + DBG("Q[%d][%d](0) = %g\n", j, m, H * S); + } + } + _2FREE(&Q0); + + // allocate memory for Z_p + _2D *Zarr = _2ALLOC(Sz, pmax + 1); // pmax is max number _including_ + double **Zd = Zarr->data; + // now we should calculate R_n^m +DBG("R_{2j}^0\n"); + // R_{2j}^0(rho) = Q_j^0(rho^2) + size_t S = Sz*sizeof(double); + for(j = 0; j <= jmax; j++){ + int pidx = p_from_nm(2*j, 0); + if(pidx > pmax) break; + DBG("R[%d]\n", pidx); + memcpy(Zd[pidx], Qd[j*mW], S); + } +DBG("R_n^n\n"); + // R_n^n(rho) = rho^n / sqrt(Sum_{i=0}^n[eps^{2i}]) + double **Rpow = build_rpowR(nmax, Sz, P); // powers of rho + double epssum = 1., epspow = 1.; + for(n = 1; n <= nmax; n++){ + int pidx = p_from_nm(n, -n); // R_n^{-n} + if(pidx > pmax) break; + DBG("R[%d] (%d, %d)\n", pidx, n, -n); + epspow *= eps2; // eps^{2n} + epssum += epspow;// sum_0^n eps^{2n} + double Sq = sqrt(epssum); // sqrt(sum...) + double *rptr = Zd[pidx], *pptr = Rpow[n]; + int pt; + for(pt = 0; pt < Sz; pt++, rptr++, pptr++) + *rptr = *pptr / Sq; // rho^n / sqrt(...) + int pidx1 = p_from_nm(n, n); // R_n^{n} + if(pidx1 > pmax) break; + DBG("R[%d] (%d, %d)\n", pidx1, n, n); + memcpy(Zd[pidx1], Zd[pidx], S); + } +DBG("R_n^m\n"); + /* / 1 - eps^2 \ + R_{2j+|m|}^{|m|}(rho, eps) = sqrt| ------------------ | * rho^m * Q_j^m(rho^2) + \ 2(2j+|m|+1)*h_j^m / */ + for(j = 1; j <= jmax; j++){ + for(m = 1; m <= mmax; m++){ + int N = 2*j + m, pt, pidx = p_from_nm(N, -m); + if(pidx > pmax) continue; + DBG("R[%d]\n", pidx); + double *rptr = Zd[pidx], *pptr = Rpow[m], *qptr = Qd[j*mW+m]; + double hjm = hd[j][m]; + for(pt = 0; pt < Sz; pt++, rptr++, pptr++, qptr++) + *rptr = sqrt((1-eps2)/2/(N+1)/hjm) * (*pptr) * (*qptr); + int pidx1 = p_from_nm(N, m); + if(pidx1 > pmax) continue; + DBG("R[%d]\n", pidx1); + memcpy(Zd[pidx1], Zd[pidx], S); + } + } + free_rpow(&Rpow, nmax); + _2FREE(&h); + _2FREE(&Q); + // last step: fill Z_n^m and find sum^2 + int p; + double *sum2 = MALLOC(double, pmax+1); + for(p = 0; p <= pmax; p++){ + convert_Zidx(p, &n, &m); + DBG("p=%d (n=%d, m=%d)\n", p, n, m); + int pt; + double (*fn)(double val); + double empty(double _U_ val){return 1.;} + if(m == 0) fn = empty; + else{if(m > 0) + fn = cos; + else + fn = sin;} + double m_abs = fabs((double)m), *rptr = Zd[p], *sptr = &sum2[p]; + polar *Pptr = P; + for(pt = 0; pt < Sz; pt++, rptr++, Pptr++){ + *rptr *= fn(m_abs * Pptr->theta); + //DBG("\tpt %d, val = %g\n", pt, *rptr); + *sptr += (*rptr) * (*rptr); + } + DBG("SUM p^2 = %g\n", *sptr); + } + if(Norm) *Norm = sum2; + else FREE(sum2); + return Zarr; +} +//#undef DBG +//#define DBG(...) do{fprintf(stderr, __VA_ARGS__); }while(0) + +/** + * Compose wavefront by koefficients of annular Zernike polynomials + * @params like for ZcomposeR + * @return composed wavefront (dynamically allocated) + */ +double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P){ + int i; + double *image = MALLOC(double, Sz); + _2D *Zannular = ann_Z(Zsz-1, Sz, P, NULL); + double **Zcoeff = Zannular->data; + for(i = 0; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + int j; + double *iptr = image, *zptr = Zcoeff[i]; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + *iptr += K * (*zptr); // add next Zernike polynomial + } + } + _2FREE(&Zannular); + return image; +} + + +/** + * Decompose wavefront by annular Zernike + * @params like ZdecomposeR + * @return array of Zernike coefficients + */ +double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + double *Zidxs = NULL, *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(double, Sz); + memcpy(icopy, heights, Sz*sizeof(double)); // make image copy to leave it unchanged + *Zsz = pmax; + double *norm; + _2D *Zannular = ann_Z(pmax-1, Sz, P, &norm); + double **Zcoeff = Zannular->data; + for(i = 0; i < pmax; i++){ // now we fill array + int j; + double *iptr = icopy, *zptr = Zcoeff[i], K = 0., norm_i = norm[i]; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += (*zptr) * (*iptr) / norm_i; // multiply matrixes to get coefficient + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + iptr = icopy; zptr = Zcoeff[i]; + for(j = 0; j < Sz; j++, iptr++, zptr++) + *iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values + } + if(lastIdx) *lastIdx = maxIdx; + _2FREE(&Zannular); + FREE(norm); + FREE(icopy); + return Zidxs; +} diff --git a/zernike_annular.h b/zernike_annular.h new file mode 100644 index 0000000..e69de29