From e9e7f10da6797b5acae846564d3eb03aa7419a03 Mon Sep 17 00:00:00 2001 From: eddyem Date: Thu, 25 Jun 2015 11:43:41 +0300 Subject: [PATCH] Zernike works --- Zernike/Makefile | 16 +- Zernike/Z-BTA_test.c | 380 +++++++++--------------- Zernike/simple_list.c | 19 -- Zernike/simple_list.h | 1 + Zernike/spots.c | 595 ++++++++++++++++++++++++++++++++++++++ Zernike/spots.h | 104 +++++++ Zernike/zern_private.h | 3 +- Zernike/zernike_annular.c | 1 + 8 files changed, 843 insertions(+), 276 deletions(-) create mode 100644 Zernike/spots.c create mode 100644 Zernike/spots.h diff --git a/Zernike/Makefile b/Zernike/Makefile index ab301ce..725991a 100644 --- a/Zernike/Makefile +++ b/Zernike/Makefile @@ -1,23 +1,15 @@ LOADLIBES = -lm -lgsl -lgslcblas -SRCS = zernike.c zernikeR.c zernike_annular.c +SRCS = zernike.c zernikeR.c zernike_annular.c Z-BTA_test.c simple_list.c spots.c CC = gcc DEFINES = -D_GNU_SOURCE #-D_XOPEN_SOURCE=501 CXX = gcc CFLAGS = -Wall -Werror $(DEFINES) OBJS = $(SRCS:.c=.o) -all : zernike btatest -zernike : $(OBJS) test.o - $(CC) $(CFLAGS) test.o $(OBJS) $(LOADLIBES) -o zernike -btatest : $(OBJS) Z-BTA_test.o simple_list.o - $(CC) $(CFLAGS) Z-BTA_test.o simple_list.o $(OBJS) $(LOADLIBES) -o btatest + +all : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LOADLIBES) -o btatest clean: /bin/rm -f *.o *~ depend: $(CXX) -MM $(SRCS) - -### -# name1.o : header1.h header2.h ... -test.o zernike.o zernikeR.o zernike_annular.o Z-BTA_test.o : zernike.h -zernike.o zernikeR.o zernike_annular.o : zern_private.h -simple_list.o Z-BTA_test.o : simple_list.h diff --git a/Zernike/Z-BTA_test.c b/Zernike/Z-BTA_test.c index a3bc18d..87a981e 100644 --- a/Zernike/Z-BTA_test.c +++ b/Zernike/Z-BTA_test.c @@ -23,45 +23,14 @@ #include #include #include -#include -#include -#include // open/close/etc -#include -#include -#include // mmap #include "zernike.h" -#include "simple_list.h" +#include "spots.h" -#define _(...) __VA_ARGS__ extern char *__progname; -#define ERR(...) err(1, __VA_ARGS__) -#define ERRX(...) errx(1, __VA_ARGS__) - -// x0 = (*spot)->c.x - AxisX; -// y0 = AxisY - (*spot)->c.y; -double AxisX = 1523., AxisY = 1533.; // BUG from fitsview : xreal = x + AxisX, yreal = y - AxisY - -char *name0 = NULL, *name1 = NULL; // prefocal/postfocal filenames -double distance = -1.; // distance between images -double pixsize = 30.; // pixel size -double ImHeight = 500.; // half of image height in pixels - -// spots for spotlist -typedef struct{ - int id; // spot identificator - double x; // coordinates of center - double y; -} spot; - -// hartmannogram -typedef struct{ - char *filename; // spots-filename - int len; // amount of spots - List *spots; // spotlist -} hartmann; - +char *name0 = NULL, *name1 = NULL; // filenames with points of pre- and postfocal images coordinates +char *gradname = NULL; // file with pre-computed gradients /** * print usage with optional message & exit with error(1) @@ -76,14 +45,15 @@ void usage(char *fmt, ...){ printf("\n\n"); } va_end(ap); - // "Использование:\t%s опции\n" + // ":\t%s \n" printf(_("Usage:\t%s options\n"), __progname); - // "Опции:\n" + // ":\n" printf(_("Required options:\n")); printf("\t--prefocal, -0 \t\tfilename with spotslist in " RED "prefocal" OLDCOLOR " image\n"); printf("\t--postfocal,-1 \t\tfilename with spotslist in " RED "postfocal" OLDCOLOR " image\n"); printf("\t--distance, -D \tdistance between images in millimeters\n"); + printf("\t--gradients, -G \t\tfilename with precomputed gradients\n"); printf(_("Unnesessary options:\n")); printf("\t--pixsize, -p \tpixel size in microns (default: 30)\n"); exit(1); @@ -112,7 +82,7 @@ double *myatod(double *num, const char *str){ */ void parse_args(int argc, char **argv){ int i; - char short_options[] = "0:1:D:p:"; // all short equivalents + char short_options[] = "0:1:D:p:G:"; // struct option long_options[] = { /* { name, has_arg, flag, val }, where: * name - name of long parameter @@ -125,7 +95,8 @@ void parse_args(int argc, char **argv){ {"prefocal", 1, 0, '0'}, {"postfocal", 1, 0, '1'}, {"distance", 1, 0, 'D'}, - {"pixsize", 1, 0, '0'}, + {"pixsize", 1, 0, 'p'}, + {"gradients", 1, 0, 'G'}, { 0, 0, 0, 0 } }; if(argc == 1){ @@ -136,9 +107,6 @@ void parse_args(int argc, char **argv){ if((opt = getopt_long(argc, argv, short_options, long_options, NULL)) == -1) break; switch(opt){ - /*case 0: // only long option - // do something? - break;*/ case '0': name0 = strdup(optarg); break; @@ -152,6 +120,10 @@ void parse_args(int argc, char **argv){ case 'p': if(!myatod(&pixsize, optarg)) usage("Parameter should be an int or floating point number!"); + else pixsize *= 1e-3; + break; + case 'G': + gradname = strdup(optarg); break; default: usage(NULL); @@ -160,243 +132,165 @@ void parse_args(int argc, char **argv){ argc -= optind; argv += optind; if(argc > 0){ - // "Игнорирую аргумент[ы]:\n" + // " []:\n" printf(_("Ignore argument[s]:\n")); for (i = 0; i < argc; i++) printf("%s ", argv[i]); printf("\n"); } - if(!name0 || !name1) - usage("You should point to both spots-files"); + if((!name0 || !name1) && !gradname) + usage("You should point to both spots-files or file with gradients"); if(distance < 0.) usage("What is the distance between images?"); } -typedef struct{ - char *data; - size_t len; -} mmapbuf; -/** - * 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) - ERR(_("Can't open %s for reading"), filename); - if(fstat (fd, &statbuf) < 0) - ERR(_("Can't stat %s"), filename); - Mlen = statbuf.st_size; - if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) - ERR(_("Mmap error for input")); - 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)) - ERR(_("Can't munmap")); - FREE(b); -} - -/** - * Read spots-file and fill hartmann structure - * @param filename - name of spots-file - * @return dynamically allocated hartmanogram structure - */ -hartmann *read_spots(char *filename){ - assert(filename); - mmapbuf *M = NULL; - int L = 0; - List *spots = NULL, *curspot = NULL; - M = My_mmap(filename); - hartmann *H = MALLOC(hartmann, 1); - H->filename = strdup(filename); - char *pos = M->data, *epos = pos + M->len; - for(; pos && pos < epos; pos = strchr(pos+1, '\n')){ - spot *Spot = MALLOC(spot, 1); - double x, y; - if(3 != sscanf(pos, "%d %*s %*s %*s %*s %lf %lf", &Spot->id, &x, &y)) - continue; - Spot->x = x + AxisX - ImHeight; - Spot->y = y - AxisY + ImHeight; - L++; - if(spots) - curspot = list_add(&curspot, LIST_T(Spot)); - else - curspot = list_add(&spots, LIST_T(Spot)); - }; - H->len = L; - H->spots = spots; - My_munmap(M); - return H; -} - -void h_free(hartmann **H){ - listfree_function(free); - list_free(&((*H)->spots)); - listfree_function(NULL); - free((*H)->filename); - FREE(*H); -} - -// temporary structure for building of coordinates-gradients list -typedef struct{ - double x; - double y; - double Dx; - double Dy; -} CG; -/** - * - * @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; - assert(H); assert(H[0]); assert(H[1]); - List *S0 = H[0]->spots, *S1; - List *CG_list = NULL, *curCG = NULL; - //printf(RED "\nspots\n" OLDCOLOR "\n"); - double Scale = pixsize * 1e-6 / distance / 2., S = 0.; // tg(2a)=dx/D -> a \approx dx/(2D) - /* - * 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(; S0; S0 = S0->next){ - spot *Sp0 = (spot*)S0->data; - int Id0 = Sp0->id; - S1 = H[1]->spots; - for(; S1; S1 = S1->next){ - spot *Sp1 = (spot*)S1->data; - if(Sp1->id > Id0) break; // point with Id0 not found - if(Sp1->id == Id0){ - CG *cg = MALLOC(CG, 1); -/* printf("id=%d (%g, %g), dx=%g, dy=%g\n", Sp0->id, Sp0->x, Sp0->y, - (Sp1->x-Sp0->x)*Scale, (Sp1->y-Sp0->y)*Scale); -*/ cg->x = Sp0->x; cg->y = Sp0->y; - cg->Dx = -(Sp1->x-Sp0->x)*Scale; - cg->Dy = -(Sp1->y-Sp0->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; -} int main(int argc, char **argv){ - int _U_ i; - double scale; - hartmann _U_ *images[2]; + int i, j; + //double scale; + hartmann *images[2]; + mirror *mir; + size_t L; + polar *coords = NULL; + point *grads = NULL; parse_args(argc, argv); - images[0] = read_spots(name0); - images[1] = read_spots(name1); - polar *coords = NULL; point *grads = NULL; - size_t _U_ L = get_gradients(images, &coords, &grads, &scale); - h_free(&images[0]); - h_free(&images[1]); -/* printf(GREEN "\nSpots:\n" OLDCOLOR "\n\tr\ttheta\tDx\tDy\n"); + if(!gradname){ + images[0] = read_spots(name0,0); + images[1] = read_spots(name1,1); + mir = calc_mir_coordinates(images); + getQ(mir, images[0]); + calc_Hartmann_constant(mir, images[0]); + spot_diagram *spot_dia = calc_spot_diagram(mir, images[0], 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); + h_free(&images[0]); + h_free(&images[1]); + }else{ +// L = read_gradients(gradname, &coords, &grads, &scale); + } +/* + // 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, grads[i].y); + grads[i].x*1000., grads[i].y*1000.); } -*/ int Zsz, lastidx; - double *Zidxs = LS_gradZdecomposeR(15, L, coords, grads, &Zsz, &lastidx); + */ + // gradients decomposition (Less squares, Zhao polinomials) + int Zsz, lastidx; + double *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("\n" RED "GradZ decompose, coefficients (%d):" OLDCOLOR "\n", lastidx); - for(i = 0; i < lastidx; i++) printf("%5.3f, ", Zidxs[i]); + printf("\nGradZ decompose, coefficients (%d):\n", lastidx); + for(i = 0; i < lastidx; i++) printf("%g, ", Zidxs[i]); printf("\n\n"); - const int GridSize = 15; + int GridSize = 15; int G2 = GridSize * GridSize; polar *rect = MALLOC(polar, G2), *rptr = rect; - double *mirZ = MALLOC(double, G2); - #define ADD 0. - int j; double Stp = 2./((double)GridSize - 1.); + double Stp = 2./((double)GridSize - 1.); for(j = 0; j < GridSize; j++){ - double y = ((double) j + ADD) * Stp - 1.; + double y = ((double) j) * Stp - 1.; for(i = 0; i < GridSize; i++, rptr++){ - double x = ((double) i + ADD)* Stp - 1.; + double x = ((double) i)* Stp - 1.; double R2 = x*x + y*y; rptr->r = sqrt(R2); rptr->theta = atan2(y, x); - //printf("x=%g, y=%g, r=%g, t=%g\n", x,y,rptr->r, rptr->theta); - if(R2 > 1.) continue; - mirZ[j*GridSize+i] = R2 / 32.; // mirror surface, z = r^2/(4f), or (z/R) = (r/R)^2/(4[f/R]) - //mirZ[j*GridSize+i] = R2; } } - printf("\n\nCoeff: %g\n\n", Zidxs[4]*sqrt(3.)); -Zidxs[4] = 0.; Zidxs[1] = 0.; Zidxs[2] = 0.; + // build uniform grid for mirror profile recalculation double *comp = ZcomposeR(lastidx, Zidxs, G2, rect); - printf("\n"); - double SS = 0.; int SScntr = 0; + double zero_val = 0., N = 0.; for(j = GridSize-1; j > -1; j--){ - for(i = 0; i < GridSize; i++){ - int idx = j*GridSize+i; - double Diff = mirZ[idx] - comp[idx]; + int idx = j*GridSize; + for(i = 0; i < GridSize; i++,idx++){ + if(comp[idx] > 1e-9){ + zero_val += comp[idx]; + N++; + } + } + } + zero_val /= N; + printf("zero: %g\n", zero_val); + for(j = GridSize-1; j > -1; j--){ + int idx = j*GridSize; + for(i = 0; i < GridSize; i++,idx++){ + //int idx = j*GridSize+i; // printf("%7.3f", Diff); - printf("%7.3f", comp[idx]*1e3); - if(rect[idx].r < 1.){ SS += Diff; SScntr++;} - //printf("%7.3f", (comp[idx] + 0.03)/ mirZ[idx]); + //if(comp[idx] > 1e-15){ + if(rect[idx].r > 1. || rect[idx].r < 0.2){ + printf("%7.3f", 0.); + }else{ + //printf("%7.3f", comp[idx] * scale); + printf("%7.3f", comp[idx] * MIR_R); + } } printf("\n"); } - SS /= SScntr; - printf("\naver: %g\n", SS); - for(j = GridSize-1; j > -1; j--){ - for(i = 0; i < GridSize; i++){ - int idx = j*GridSize+i; - //printf("%7.3f", (comp[idx] + SS) / mirZ[idx]); - double Z = (fabs(comp[idx]) < 2*DBL_EPSILON) ? 0. : comp[idx] + SS - mirZ[idx]; - printf("%7.3f", Z * 1e3); + + // save matrix 100x100 with mirror deformations in Octave file format + FILE *f = fopen("file.out", "w"); + GridSize = 100; + G2 = GridSize * GridSize; + if(f){ + FREE(rect); + rect = MALLOC(polar, G2); rptr = rect; + Stp = 2./((double)GridSize - 1.); + for(j = 0; j < GridSize; j++){ + double y = ((double) j) * Stp - 1.; + for(i = 0; i < GridSize; i++, rptr++){ + double x = ((double) i)* Stp - 1.; + double R2 = x*x + y*y; + rptr->r = sqrt(R2); + rptr->theta = atan2(y, x); + } + } + fprintf(f, "# generated by Z-BTA_test\n" + "# name: dev_matrix\n# type: matrix\n# rows: %d\n# columns: %d\n", + GridSize, GridSize); + FREE(comp); + comp = ZcomposeR(lastidx, Zidxs, G2, rect); // 100x100 + 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 ", comp[idx] * MIR_R); + } + fprintf(f,"\n"); } - printf("\n"); } - FREE(comp); FREE(Zidxs); + fclose(f); + FREE(comp); + FREE(Zidxs); FREE(coords); FREE(grads); return 0; } - - diff --git a/Zernike/simple_list.c b/Zernike/simple_list.c index cd507d3..479906a 100644 --- a/Zernike/simple_list.c +++ b/Zernike/simple_list.c @@ -91,22 +91,3 @@ void list_free(List **root){ *root = NULL; } -#ifdef STANDALONE -int main(void) { - List *tp = NULL, *root_p = NULL; - int i, ins[] = {4,2,6,1,3,4,7}; - for(i = 0; i < 7; i++){ - if(!(tp = list_add(&tp, ins[i]))) - err(1, "Malloc error"); // can't insert - if(!root_p) root_p = tp; - } - tp = root_p; - i = 0; - do{ - printf("element %d = %d\n", i++, tp->data); - tp = tp->next; - }while(tp); - list_free(&root_p); - return 0; -} -#endif diff --git a/Zernike/simple_list.h b/Zernike/simple_list.h index 643ee4f..95cc2e6 100644 --- a/Zernike/simple_list.h +++ b/Zernike/simple_list.h @@ -18,6 +18,7 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, * MA 02110-1301, USA. */ + #pragma once #ifndef __SIMPLE_LIST_H__ #define __SIMPLE_LIST_H__ diff --git a/Zernike/spots.c b/Zernike/spots.c new file mode 100644 index 0000000..c76be9c --- /dev/null +++ b/Zernike/spots.c @@ -0,0 +1,595 @@ +/* + * 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" + +#define MM_TO_ARCSEC(x) (x*206265./FOCAL_R) + +const double _4F = 4. * FOCAL_R; +const double _2F = 2. * FOCAL_R; + + +#define ERR(...) err(1, __VA_ARGS__) +#define ERRX(...) errx(1, __VA_ARGS__) + +double pixsize = 30.e-3; // CCD pixel size in millimeters +double distance = -1.; // distanse between pre- and postfocal images in millimeters + +typedef struct{ + char *data; + size_t len; +} mmapbuf; +/** + * 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) + ERR(_("Can't open %s for reading"), filename); + if(fstat (fd, &statbuf) < 0) + ERR(_("Can't stat %s"), filename); + Mlen = statbuf.st_size; + if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) + ERR(_("Mmap error for input")); + 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)) + ERR(_("Can't munmap")); + FREE(b); +} + +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, int prefocal){ + 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 * pixsize; + H->center.y = yc * pixsize; + printf("get common center: (%.1f. %.1f)mm (pixsize=%g)\n", H->center.x, H->center.y, 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) * pixsize; + spots[i].y = (spots[i].y - yc) * 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, hartmann *H){ + uint8_t *mirgot = mir->got; + int i, j; + point *tans = mir->tans, *pre = H->spots; + 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 prefoc (i) - prefocal hartmannogram + */ +void getQ(mirror *mir, hartmann *prefoc){ + point *pre = prefoc->spots; + point *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, hartmann *prefoc, double z){ + int i; + spot_diagram *SD = MALLOC(spot_diagram, 1); + point *pre = prefoc->spots, *spots = SD->spots; + point *tans = mir->tans; + uint8_t *got = mir->got; + memcpy(SD->got, mir->got, 258); + SD->center.x = prefoc->center.x + mir->tanc.x * z; + SD->center.y = 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 = 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/Zernike/spots.h b/Zernike/spots.h new file mode 100644 index 0000000..b691924 --- /dev/null +++ b/Zernike/spots.h @@ -0,0 +1,104 @@ +/* + * 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" + +// 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) + +#define _(...) __VA_ARGS__ +#define MirRadius 3. // main mirror radius + +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 +} 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 + 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, int prefocal); +void h_free(hartmann **H); +mirror *calc_mir_coordinates(hartmann *H[]); +void getQ(mirror *mir, hartmann *prefoc); +double calc_Hartmann_constant(mirror *mir, hartmann *H); +spot_diagram *calc_spot_diagram( mirror *mir, hartmann *H, 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/Zernike/zern_private.h b/Zernike/zern_private.h index dd710be..e057326 100644 --- a/Zernike/zern_private.h +++ b/Zernike/zern_private.h @@ -18,6 +18,7 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, * MA 02110-1301, USA. */ + #pragma once #ifndef __ZERN_PRIVATE_H__ #define __ZERN_PRIVATE_H__ @@ -44,8 +45,6 @@ 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); diff --git a/Zernike/zernike_annular.c b/Zernike/zernike_annular.c index 10ff7ae..7ad2c67 100644 --- a/Zernike/zernike_annular.c +++ b/Zernike/zernike_annular.c @@ -19,6 +19,7 @@ * MA 02110-1301, USA. */ + /* * These polynomials realized according to formulae in: @ARTICLE{1981JOSA...71...75M,