mirror of
https://github.com/eddyem/eddys_snippets.git
synced 2026-03-22 01:31:16 +03:00
Zernike works
This commit is contained in:
eddyem
@@ -23,45 +23,14 @@
|
||||
#include <getopt.h>
|
||||
#include <stdarg.h>
|
||||
#include <string.h>
|
||||
#include <assert.h>
|
||||
#include <err.h>
|
||||
#include <sys/stat.h> // open/close/etc
|
||||
#include <fcntl.h>
|
||||
#include <unistd.h>
|
||||
#include <sys/mman.h> // 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 <file>\t\tfilename with spotslist in " RED "prefocal" OLDCOLOR " image\n");
|
||||
printf("\t--postfocal,-1 <file>\t\tfilename with spotslist in " RED "postfocal" OLDCOLOR " image\n");
|
||||
printf("\t--distance, -D <distance in mm>\tdistance between images in millimeters\n");
|
||||
printf("\t--gradients, -G <file>\t\tfilename with precomputed gradients\n");
|
||||
printf(_("Unnesessary options:\n"));
|
||||
printf("\t--pixsize, -p <size in mkm>\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 <pixsize> 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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user