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