From c28092f7908dc88d1e8c4fe3efa3924fbe192051 Mon Sep 17 00:00:00 2001 From: Rick Kirian Date: Fri, 25 Mar 2011 16:24:19 +0100 Subject: calibrate_detector.c: work in progress --- src/calibrate_detector.c | 425 ++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 349 insertions(+), 76 deletions(-) (limited to 'src/calibrate_detector.c') diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c index e8b5ed37..35e0672c 100644 --- a/src/calibrate_detector.c +++ b/src/calibrate_detector.c @@ -20,111 +20,94 @@ #include #include #include +#include +#include +#include #include "utils.h" #include "image.h" #include "detector.h" #include "index.h" #include "hdf5-file.h" +#include "stream.h" +#include "peaks.h" + static void show_help(const char *s) { printf("Syntax: %s [options] -i \n\n", s); printf( -"Optimise detector geometry.\n" +"Stream-based optimisation of detector geometry.\n" "\n" " -h, --help Display this help message.\n" " -g. --geometry= Get detector geometry from file.\n" " -i, --input= Input filename.\n" -" --split Break up input file into panels according to\n" -" the detector geometry description.\n" +" -m, --method= The calibration method. Possiblities are\n" +" xy Determine panel shifts in plane of detector\n" +" (i.e. ignore camera length calibration\n" +" -o, --output= Output results here" +" -n, --npeaks= Don't refine unless this many peaks observed\n" +" (in the whole stream, not a single shot)\n" "\n"); } -static void split_image(struct image *image) -{ - int i; - const struct detector *det = image->det; - - for ( i=0; in_panels; i++ ) { - - int x, y; - struct panel *p; - float *data; - int width, height; - char filename[1024]; - - p = &det->panels[i]; - - width = 1 + p->max_fs - p->min_fs; - height = 1 + p->max_ss - p->min_ss; - data = malloc(width*height*sizeof(float)); - snprintf(filename, 1023, "%s-%i.h5", image->filename, i); - - STATUS("Panel %i, %i by %i pixels -> %s\n", i, width, height, - filename); - - for ( x=0; xmin_fs; - im_y = y+p->min_ss; - - data[x+width*y] = image->data[im_x+image->width*im_y]; - - } - } - - hdf5_write(filename, data, width, height, H5T_NATIVE_FLOAT); - - free(data); - - } -} - int main(int argc, char *argv[]) { - int c; + + + char c; struct image image; - int x, y; - struct hdfile *hdfile; + struct panel p; char *filename = NULL; char *geometry = NULL; - int split = 0; + char *method = NULL; + char *outputfile = NULL; + FILE *fh = NULL; + FILE *outfh = NULL; + int nFeatures; + int i; + int fail; + int nChunks; + int minpeaks=0; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, {"geometry", 1, NULL, 'g'}, - {"split", 0, &split, 1}, + {"method", 1, NULL, 'm'}, + {"output", 1, NULL, 'o'}, + {"npeaks", 1, NULL, 'n'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hi:g:m:o:n:", longopts, NULL)) != -1) { switch (c) { case 'h' : show_help(argv[0]); return 0; - case 0 : break; - case 'i' : filename = strdup(optarg); break; - case 'g' : geometry = strdup(optarg); break; - + case 'm' : + method = strdup(optarg); + break; + case 'o' : + outputfile = strdup(optarg); + break; + case 'n' : + minpeaks = atoi(optarg); + break; default : return 1; } @@ -135,45 +118,335 @@ int main(int argc, char *argv[]) ERROR("You must specify the input filename with -i\n"); return 1; } - image.filename = filename; + fh = fopen(filename,"r"); + if ( fh == NULL ) { + ERROR("Problem opening file\n"); + return 1; + } + printf("Read stream file: %s\n",filename); if ( geometry == NULL ) { ERROR("You need to specify a geometry file with --geometry\n"); return 1; - } + } image.det = get_detector_geometry(geometry); if ( image.det == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", geometry); + ERROR("Failed to read detector geometry from %s\n", geometry); return 1; } + printf("Read geometry file: %s\n",geometry); + image.width = image.det->max_fs; + image.height = image.det->max_ss; free(geometry); - hdfile = hdfile_open(filename); - hdfile_set_image(hdfile, "/data/data"); - hdf5_read(hdfile, &image, 1); - - if ( split ) { - split_image(&image); - exit(0); + if ( !(outputfile == NULL) ) { + printf("Writing result to file: %s\n",outputfile); + outfh = fopen(outputfile,"w"); + } else { + ERROR("You did not specify an output file.\n"); + return 1; + } + if ( minpeaks == 0 ) { + minpeaks = 1; + printf("You did not specify minimum number of peaks." + " Setting default value of %d\n",minpeaks); + } + + if ( method == NULL ) { + printf("You did not specify a refinement method-- setting default.\n"); + method = strdup("xy"); } - for ( x=0; xn_panels); + weightedSumSS = (double *)calloc(sizeof(double), image.det->n_panels); + summedWeights = (double *)calloc(sizeof(double), image.det->n_panels); + peaksFound = (int *)calloc(sizeof(int), image.det->n_panels); + meanShiftFS = (double *)calloc(sizeof(double), image.det->n_panels); + meanShiftSS = (double *)calloc(sizeof(double), image.det->n_panels); + + /* initialize arrays (is there a standard function to do this?) */ + int pi; + for (pi=0; pin_panels; pi++) { + weightedSumFS[pi] = 0; + weightedSumSS[pi] = 0; + summedWeights[pi] = 0; + peaksFound[pi] = 0; + meanShiftFS[pi] = 0; + meanShiftSS[pi] = 0; + } - r = get_q(&image, x, y, NULL, 1.0/image.lambda); - q = modulus(r.u, r.v, r.w); - intensity = image.data[x + image.width*y]; + fesetround(1); /* Round towards nearest */ + + nChunks = 0; + while (1) { + + // check for peaks in next file + fail = read_chunk(fh, &image); + if ( fail == 1 ) { + // FIXME:should check if this is EOF, or broken file handle + break; + } + nChunks += 1; + + + // move on to the next chunk if no peaks found + if ( image.features == NULL ) { + continue; + } + + /* update reflection list since geometry may have changed + from what is in the data stream (e.g. if iterating this + calibration procedure) */ + //image.reflections = + // find_projected_peaks(&image,image.indexed_cell,0, 0.1); + + //printf("chunk %d\n",nChunks); + //cell_print(image.indexed_cell); + + // now loop through peaks to determine mean panel shift + nFeatures = image_feature_count(image.features); + + for (i=0; ifs; + ss = thisFeature->ss; + + p = find_panel(image.det, fs, ss); + if ( p == NULL ) { + continue; + } + if ( p->no_index ) continue; + + /* now determine the predicted peak position */ + q = get_q(&image, fs, ss, &twotheta, 1.0/image.lambda); + + /* miller indices of nearest Bragg reflection */ + cell_get_cartesian(image.indexed_cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + + /* now get scattering vector for reflectin [hkl] + * this means solving the equation U*q = h */ + + double U[] = {ax, ay, az, + bx, by, bz, + cx, cy, cz}; + + double hvec[] = {h,k,l}; + + gsl_matrix_view m + = gsl_matrix_view_array (U, 3, 3); + + gsl_vector_view b + = gsl_vector_view_array (hvec, 3); + + gsl_vector *xx = gsl_vector_alloc (3); + + int s; + + gsl_permutation * perm = gsl_permutation_alloc (3); + gsl_linalg_LU_decomp (&m.matrix, perm, &s); + gsl_linalg_LU_solve (&m.matrix, perm, &b.vector, xx); + + + // outgoing wavevector + double x = xx->data[0]; + double y = xx->data[1]; + double z = xx->data[2]; + + double kk; + double xd, yd, cl; + double plx, ply; + + kk = 1/image.lambda; + const double den = kk + z; + + /* Camera length for this panel */ + cl = p->clen; + + /* Coordinates of peak relative to central beam, in m */ + xd = cl * x / den; + yd = cl * y / den; + + /* Convert to pixels */ + xd *= p->res; + yd *= p->res; + + /* Convert to relative to the panel corner */ + plx = xd - p->cnx; + ply = yd - p->cny; + + pfs = p->xfs*plx + p->yfs*ply; + pss = p->xss*plx + p->yss*ply; + + pfs += p->min_fs; + pss += p->min_ss; + + /* Now, is this on this panel? */ + if ( fs < p->min_fs ) continue; + if ( fs > p->max_fs ) continue; + if ( ss < p->min_ss ) continue; + if ( ss > p->max_ss ) continue; + + /* Finally, we have the shift in position of this peak */ + dfs = pfs - fs; + dss = pss - ss; + + /* Add this shift to the weighted sum over shifts */ + pi = find_panel_number(image.det,fs,ss); + thisWeight = 1; // FIXME: use real weighting some day + weightedSumFS[pi] += thisWeight*dfs; + weightedSumSS[pi] += thisWeight*dss; + summedWeights[pi] += thisWeight; + peaksFound[pi] += 1; + + } /* end loop through image features */ + + } /* end loop through stream chunks */ + + + + /* calculate weighted average shift in peak positions */ + for (pi=0; pi < image.det->n_panels; pi++) { + meanShiftFS[pi] = weightedSumFS[pi]/summedWeights[pi]; + meanShiftSS[pi] = weightedSumSS[pi]/summedWeights[pi]; + } + + /* now generate a new geometry file */ + for (pi=0; pi < image.det->n_panels; pi++) { + + p = image.det->panels[pi]; + + xsh = 0; + ysh = 0; + + if ( peaksFound[pi] >= minpeaks ) { + + //printf("meanShift: %f %f\n",meanShiftFS[pi],meanShiftSS[pi]); + + /* convert shifts from raw coords to lab frame */ + xsh = meanShiftFS[pi]*p.fsx + meanShiftSS[pi]*p.ssx; + ysh = meanShiftFS[pi]*p.fsy + meanShiftSS[pi]*p.ssy; + + /* add shifts to original panel corner locations */ + cnx = p.cnx + xsh; + cny = p.cny + ysh; + + //printf("new panel shift: %f %f\n",cnx,cny); + + + } else { + + /* not refined: use original coordinates */ + cnx = p.cnx; + cny = p.cny; + + } + + printf("panel %s, # peaks: %10d, mean shifts: %f %f\n",p.name, peaksFound[pi],xsh,ysh); + + //FIXME: there should be a function in geometry.c to write + // these values to text file, since it will be useful on + // other places as well (e.g. writing the geometry to + // the data stream) + fprintf(outfh,"%s/min_fs = %d\n",p.name,p.min_fs); + fprintf(outfh,"%s/min_ss = %d\n",p.name,p.min_ss); + fprintf(outfh,"%s/max_fs = %d\n",p.name,p.max_fs); + fprintf(outfh,"%s/max_ss = %d\n",p.name,p.max_ss); + fprintf(outfh,"%s/badrow_direction = %C\n",p.name,p.badrow); + fprintf(outfh,"%s/res = %g\n",p.name,p.res); + fprintf(outfh,"%s/peak_sep = %g\n",p.name,p.peak_sep); + fprintf(outfh,"%s/clen = %s\n",p.name,p.clen_from); + //FIXME: the following is sketchy, but it will work for now. we need + // to generalise the parser in detector.c + char coord; + char sign; + if (p.fsx != 0){ + if (p.fsx>0){sign='+';}else{sign='-';} + coord = 'x'; + } else { + if (p.fsy>0){sign='+';}else{sign='-';} + coord = 'y'; + } + fprintf(outfh,"%s/fs = %C%C\n",p.name, sign, coord); + if (p.ssx != 0){ + if (p.ssx>0){sign='+';}else{sign='-';} + coord = 'x'; + } else { + if (p.ssy>0){sign='+';}else{sign='-';} + coord = 'y'; + } + fprintf(outfh,"%s/ss = %C%C\n",p.name, sign, coord); + fprintf(outfh,"%s/corner_x = %g\n",p.name,cnx); + fprintf(outfh,"%s/corner_y = %g\n",p.name,cny); + if ( peaksFound[pi] < minpeaks ) { + fprintf(outfh,"%s/no_index = %d\n",p.name,1); + } else { + fprintf(outfh,"%s/no_index = %d\n",p.name,p.no_index); + } + fprintf(outfh,"\n\n"); + + } - printf("%5i\t%5i\t%7.3f\t%7i\n", x, y, q/1.0e9, intensity); + } else { + + printf("Refinement method %s not recognized\n",method); + return 1; } - } + + if ( !(outputfile == NULL) ) { + fclose(outfh); + } + + + + printf("Done!\n"); return 0; } -- cgit v1.2.3