diff options
Diffstat (limited to 'libcrystfel/src/detector.c')
| -rw-r--r-- | libcrystfel/src/detector.c | 2362 |
1 files changed, 0 insertions, 2362 deletions
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c deleted file mode 100644 index 5653f0f2..00000000 --- a/libcrystfel/src/detector.c +++ /dev/null @@ -1,2362 +0,0 @@ -/* - * detector.c - * - * Detector properties - * - * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * Copyright © 2012 Richard Kirian - * - * Authors: - * 2009-2019 Thomas White <taw@physics.org> - * 2014 Valerio Mariani - * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de> - * 2011 Andrew Aquila - * 2011 Richard Kirian <rkirian@asu.edu> - * - * This file is part of CrystFEL. - * - * CrystFEL 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 3 of the License, or - * (at your option) any later version. - * - * CrystFEL 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 CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#include <stdlib.h> -#include <math.h> -#include <stdio.h> -#include <string.h> -#include <assert.h> -#include <ctype.h> -#include <sys/stat.h> - -#include "image.h" -#include "utils.h" -#include "detector.h" -#include "hdf5-file.h" - - -/** - * \file detector.h - */ - - -struct rg_definition { - char *name; - char *pns; -}; - - -struct rgc_definition { - char *name; - char *rgs; -}; - - -static int atob(const char *a) -{ - if ( strcasecmp(a, "true") == 0 ) return 1; - if ( strcasecmp(a, "false") == 0 ) return 0; - return atoi(a); -} - - -static int assplode_algebraic(const char *a_orig, char ***pbits) -{ - int len, i; - int nexp; - char **bits; - char *a; - int idx, istr; - - len = strlen(a_orig); - - /* Add plus at start if no sign there already */ - if ( (a_orig[0] != '+') && (a_orig[0] != '-') ) { - len += 1; - a = malloc(len+1); - snprintf(a, len+1, "+%s", a_orig); - a[len] = '\0'; - - } else { - a = strdup(a_orig); - } - - /* Count the expressions */ - nexp = 0; - for ( i=0; i<len; i++ ) { - if ( (a[i] == '+') || (a[i] == '-') ) nexp++; - } - - bits = calloc(nexp, sizeof(char *)); - - /* Break the string up */ - idx = -1; - istr = 0; - assert((a[0] == '+') || (a[0] == '-')); - for ( i=0; i<len; i++ ) { - - char ch; - - ch = a[i]; - - if ( (ch == '+') || (ch == '-') ) { - if ( idx >= 0 ) bits[idx][istr] = '\0'; - idx++; - bits[idx] = malloc(len+1); - istr = 0; - } - - if ( !isdigit(ch) && (ch != '.') && (ch != '+') && (ch != '-') - && (ch != 'x') && (ch != 'y') && (ch != 'z') ) - { - ERROR("Invalid character '%c' found.\n", ch); - return 0; - } - - assert(idx >= 0); - bits[idx][istr++] = ch; - - } - if ( idx >= 0 ) bits[idx][istr] = '\0'; - - *pbits = bits; - free(a); - - return nexp; -} - - -/* Parses the scan directions (accounting for possible rotation) - * Assumes all white spaces have been already removed */ -static int dir_conv(const char *a, double *sx, double *sy, double *sz) -{ - int n; - char **bits; - int i; - - *sx = 0.0; *sy = 0.0; *sz = 0.0; - - n = assplode_algebraic(a, &bits); - - if ( n == 0 ) { - ERROR("Invalid direction '%s'\n", a); - return 1; - } - - for ( i=0; i<n; i++ ) { - - int len; - double val; - char axis; - int j; - - len = strlen(bits[i]); - assert(len != 0); - axis = bits[i][len-1]; - if ( (axis != 'x') && (axis != 'y') && (axis != 'z') ) { - ERROR("Invalid symbol '%c' - must be x, y or z.\n", - axis); - return 1; - } - - /* Chop off the symbol now it's dealt with */ - bits[i][len-1] = '\0'; - - /* Check for anything that isn't part of a number */ - for ( j=0; j<strlen(bits[i]); j++ ) { - if ( isdigit(bits[i][j]) ) continue; - if ( bits[i][j] == '+' ) continue; - if ( bits[i][j] == '-' ) continue; - if ( bits[i][j] == '.' ) continue; - ERROR("Invalid coefficient '%s'\n", bits[i]); - } - - if ( strlen(bits[i]) == 0 ) { - val = 1.0; - } else { - val = atof(bits[i]); - } - if ( strlen(bits[i]) == 1 ) { - if ( bits[i][0] == '+' ) val = 1.0; - if ( bits[i][0] == '-' ) val = -1.0; - } - switch ( axis ) { - - case 'x' : - *sx += val; - break; - - case 'y' : - *sy += val; - break; - - case 'z' : - *sz += val; - break; - } - - free(bits[i]); - - } - free(bits); - - return 0; -} - - -static int count_trailing_spaces(const char *string) { - - int i; - - for ( i=0; i<strlen(string); i++ ) { - if ( !isspace(string[i]) ) { - return i; - } - } - - return -1; -} - - -static void build_output_line(const char *line, char *new_line, - const char *string_to_write) -{ - int nsp, i, w, ww; - int trailsp; - int n_bits; - char **bits; - - n_bits = assplode(line, "=", &bits, ASSPLODE_NONE); - trailsp = count_trailing_spaces(bits[1]); - - strcat(new_line, bits[0]); - strcat(new_line, "="); - - for ( nsp=0; nsp < trailsp; nsp++ ) { - strcat(new_line, " "); - } - strcat(new_line, string_to_write); - - for ( w=0; w<strlen(line); w++ ) { - if ( strncmp(&line[w],";",1) == 0 ) { - for ( ww=w-1; ww>=0; ww-- ) { - if ( !isspace(line[ww])) { - strcat(new_line, &line[ww]); - break; - } - } - break; - } - } - for ( i=0; i<n_bits; i++) free(bits[i]); - free(bits); -} - - -struct rvec get_q_for_panel(struct panel *p, double fs, double ss, - double *ttp, double k) -{ - struct rvec q; - double ctt, twotheta; - double xs, ys, zs; - double az; - - /* Calculate 3D position of given position, in m */ - xs = (p->cnx + fs*p->fsx + ss*p->ssx) / p->res; - ys = (p->cny + fs*p->fsy + ss*p->ssy) / p->res; - zs = p->clen + (fs*p->fsz + ss*p->ssz) / p->res; - - ctt = zs/sqrt(xs*xs + ys*ys + zs*zs); - twotheta = acos(ctt); - az = atan2(ys, xs); - if ( ttp != NULL ) *ttp = twotheta; - - q.u = k * sin(twotheta)*cos(az); - q.v = k * sin(twotheta)*sin(az); - q.w = k * (ctt - 1.0); - - return q; -} - - -int in_bad_region(struct detector *det, struct panel *p, double fs, double ss) -{ - double rx, ry; - double xs, ys; - int i; - - /* No panel found -> definitely bad! */ - if ( p == NULL ) return 1; - - /* Convert xs and ys, which are in fast scan/slow scan coordinates, - * to x and y */ - xs = fs*p->fsx + ss*p->ssx; - ys = fs*p->fsy + ss*p->ssy; - - rx = xs + p->cnx; - ry = ys + p->cny; - - for ( i=0; i<det->n_bad; i++ ) { - - struct badregion *b = &det->bad[i]; - - if ( (b->panel != NULL) - && (strcmp(b->panel, p->name) != 0) ) continue; - - if ( b->is_fsss ) { - - int nfs, nss; - - /* fs/ss bad regions are specified according to the - * original coordinates */ - nfs = fs + p->orig_min_fs; - nss = ss + p->orig_min_ss; - - if ( nfs < b->min_fs ) continue; - if ( nfs > b->max_fs ) continue; - if ( nss < b->min_ss ) continue; - if ( nss > b->max_ss ) continue; - - } else { - - if ( rx < b->min_x ) continue; - if ( rx > b->max_x ) continue; - if ( ry < b->min_y ) continue; - if ( ry > b->max_y ) continue; - - } - - return 1; - } - - return 0; -} - - -int detector_has_clen_references(struct detector *det) -{ - int i; - - for ( i=0; i<det->n_panels; i++ ) { - if ( det->panels[i].clen_from != NULL ) return 1; - } - - return 0; -} - - -static void record_panel(struct panel *p, float *dp, int do_poisson, - gsl_rng *rng, double ph_per_e, double background, - double lambda, - int *n_neg1, int *n_inf1, int *n_nan1, - int *n_neg2, int *n_inf2, int *n_nan2, - double *max_tt) -{ - int fs, ss; - - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - double counts; - double cf; - double intensity, sa; - double pix_area, Lsq; - double xs, ys, rx, ry; - double dsq, proj_area; - float dval; - double twotheta; - - intensity = (double)dp[fs + p->w*ss]; - if ( isinf(intensity) ) (*n_inf1)++; - if ( intensity < 0.0 ) (*n_neg1)++; - if ( isnan(intensity) ) (*n_nan1)++; - - /* Area of one pixel */ - pix_area = pow(1.0/p->res, 2.0); - Lsq = pow(p->clen, 2.0); - - /* Calculate distance from crystal to pixel */ - xs = fs*p->fsx + ss*p->ssx; - ys = ss*p->fsy + ss*p->ssy; - rx = (xs + p->cnx) / p->res; - ry = (ys + p->cny) / p->res; - dsq = pow(rx, 2.0) + pow(ry, 2.0); - twotheta = atan2(sqrt(dsq), p->clen); - - /* Area of pixel as seen from crystal (approximate) */ - proj_area = pix_area * cos(twotheta); - - /* Projected area of pixel divided by distance squared */ - sa = proj_area / (dsq + Lsq); - - if ( do_poisson ) { - counts = poisson_noise(rng, intensity * ph_per_e * sa); - } else { - cf = intensity * ph_per_e * sa; - counts = cf; - } - - /* Number of photons in pixel */ - dval = counts + poisson_noise(rng, background); - - /* Convert to ADU */ - dval *= p->adu_per_photon; - - /* Saturation */ - if ( dval > p->max_adu ) dval = p->max_adu; - - dp[fs + p->w*ss] = dval; - - /* Sanity checks */ - if ( isinf(dp[fs + p->w*ss]) ) n_inf2++; - if ( isnan(dp[fs + p->w*ss]) ) n_nan2++; - if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++; - - if ( twotheta > *max_tt ) *max_tt = twotheta; - - } - } -} - - -void record_image(struct image *image, int do_poisson, double background, - gsl_rng *rng, double beam_radius, double nphotons) -{ - double total_energy, energy_density; - double ph_per_e; - double area; - double max_tt = 0.0; - int n_inf1 = 0; - int n_neg1 = 0; - int n_nan1 = 0; - int n_inf2 = 0; - int n_neg2 = 0; - int n_nan2 = 0; - int pn; - - /* How many photons are scattered per electron? */ - area = M_PI*pow(beam_radius, 2.0); - total_energy = nphotons * ph_lambda_to_en(image->lambda); - energy_density = total_energy / area; - ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0); - STATUS("Fluence = %8.2e photons, " - "Energy density = %5.3f kJ/cm^2, " - "Total energy = %5.3f microJ\n", - nphotons, energy_density/1e7, total_energy*1e6); - - fill_in_adu(image); - - for ( pn=0; pn<image->det->n_panels; pn++ ) { - - record_panel(&image->det->panels[pn], image->dp[pn], - do_poisson, rng, ph_per_e, background, - image->lambda, - &n_neg1, &n_inf1, &n_nan1, - &n_neg2, &n_inf2, &n_nan2, &max_tt); - } - - STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n", - rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); - - STATUS("Halve the d values to get the voxel size for a synthesis.\n"); - - if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) { - ERROR("WARNING: The raw calculation produced %i negative" - " values, %i infinities and %i NaNs.\n", - n_neg1, n_inf1, n_nan1); - ERROR("WARNING: After processing, there were %i negative" - " values, %i infinities and %i NaNs.\n", - n_neg2, n_inf2, n_nan2); - } -} - - -signed int find_orig_panel_number(struct detector *det, double fs, double ss) -{ - int p; - - for ( p=0; p<det->n_panels; p++ ) { - if ( (fs >= det->panels[p].orig_min_fs) - && (fs < det->panels[p].orig_max_fs+1) - && (ss >= det->panels[p].orig_min_ss) - && (ss < det->panels[p].orig_max_ss+1) ) return p; - } - - return -1; -} - - -/* Like find_panel(), but uses the original panel bounds, i.e. referring to - * what's in the HDF5 file */ -struct panel *find_orig_panel(struct detector *det, double fs, double ss) -{ - signed int pn = find_orig_panel_number(det, fs, ss); - if ( pn == -1 ) return NULL; - return &det->panels[pn]; -} - - -int panel_number(const struct detector *det, const struct panel *p) -{ - int pn; - - for ( pn=0; pn<det->n_panels; pn++ ) { - if ( &det->panels[pn] == p ) return pn; - } - - return det->n_panels; -} - - -void adjust_centering_for_rail(struct panel *p) -{ - double offs; - - /* Offset in +z direction from calibrated clen to actual */ - offs = p->clen - p->clen_for_centering; - p->cnx += p->rail_x * offs; - p->cny += p->rail_y * offs; - p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs; -} - - -void fill_in_adu(struct image *image) -{ - int i; - - if ( image->det == NULL ) return; - - for ( i=0; i<image->det->n_panels; i++ ) { - - struct panel *p = &image->det->panels[i]; - - /* Already have ADU per photon? */ - if ( !isnan(p->adu_per_photon) ) continue; - - if ( isnan(p->adu_per_eV) ) { - ERROR("Neither adu_per_eV nor adu_per_photon set for " - "panel %s\n", p->name); - continue; - } - - /* Convert ADU per eV to ADU per photon */ - p->adu_per_photon = ph_lambda_to_eV(image->lambda) - * p->adu_per_eV; - } -} - - -int panel_is_in_rigid_group(const struct rigid_group *rg, struct panel *p) -{ - int i; - - for ( i=0; i<rg->n_panels; i++ ) { - if ( rg->panels[i] == p ) { - return 1; - } - } - - return 0; -} - - -int rigid_group_is_in_collection(struct rg_collection *c, - struct rigid_group *rg) -{ - int i; - - for ( i=0; i<c->n_rigid_groups; i++ ) { - if ( c->rigid_groups[i] == rg ) { - return 1; - } - } - - return 0; -} - - -struct rg_collection *find_rigid_group_collection_by_name(struct detector *det, - const char *name) -{ - int i; - - for ( i=0; i<det->n_rg_collections; i++ ) { - if ( strcmp(det->rigid_group_collections[i]->name, - name) == 0 ) { - return det->rigid_group_collections[i]; - } - } - - return NULL; -} - - -static struct panel *new_panel(struct detector *det, const char *name) -{ - struct panel *new; - - det->n_panels++; - det->panels = realloc(det->panels, det->n_panels*sizeof(struct panel)); - - new = &det->panels[det->n_panels-1]; - memcpy(new, &det->defaults, sizeof(struct panel)); - - strcpy(new->name, name); - - /* Create a new copy of the camera length location if needed */ - if ( new->clen_from != NULL ) { - new->clen_from = strdup(new->clen_from); - } - - /* Create a new copy of the data location if needed */ - if ( new->data != NULL ) { - new->data = strdup(new->data); - } - - /* Create a new copy of the dim_structure if needed */ - if ( new->dim_structure != NULL ) { - - struct dim_structure *dim_copy; - int di; - - dim_copy = initialize_dim_structure(); - dim_copy->num_dims = new->dim_structure->num_dims; - dim_copy->dims = malloc(dim_copy->num_dims*sizeof(int)); - for ( di=0; di<dim_copy->num_dims; di++ ) { - dim_copy->dims[di] = new->dim_structure->dims[di]; - } - - new->dim_structure = dim_copy; - } - - /* Create a new copy of the bad pixel mask location */ - if ( new->mask != NULL ) { - new->mask = strdup(new->mask); - } - - return new; -} - - -static struct badregion *new_bad_region(struct detector *det, const char *name) -{ - struct badregion *new; - - det->n_bad++; - det->bad = realloc(det->bad, det->n_bad*sizeof(struct badregion)); - - new = &det->bad[det->n_bad-1]; - new->min_x = NAN; - new->max_x = NAN; - new->min_y = NAN; - new->max_y = NAN; - new->min_fs = 0; - new->max_fs = 0; - new->min_ss = 0; - new->max_ss = 0; - new->is_fsss = 99; /* Slightly nasty: means "unassigned" */ - new->panel = NULL; - strcpy(new->name, name); - - return new; -} - - -struct panel *find_panel_by_name(struct detector *det, const char *name) -{ - int i; - - for ( i=0; i<det->n_panels; i++ ) { - if ( strcmp(det->panels[i].name, name) == 0 ) { - return &det->panels[i]; - } - } - - return NULL; -} - - -static struct badregion *find_bad_region_by_name(struct detector *det, - const char *name) -{ - int i; - - for ( i=0; i<det->n_bad; i++ ) { - if ( strcmp(det->bad[i].name, name) == 0 ) { - return &det->bad[i]; - } - } - - return NULL; -} - - -static struct rigid_group *find_or_add_rg(struct detector *det, - const char *name) -{ - int i; - struct rigid_group **new; - struct rigid_group *rg; - - for ( i=0; i<det->n_rigid_groups; i++ ) { - - if ( strcmp(det->rigid_groups[i]->name, name) == 0 ) { - return det->rigid_groups[i]; - } - - } - - new = realloc(det->rigid_groups, - (1+det->n_rigid_groups)*sizeof(struct rigid_group *)); - if ( new == NULL ) return NULL; - - det->rigid_groups = new; - - rg = malloc(sizeof(struct rigid_group)); - if ( rg == NULL ) return NULL; - - det->rigid_groups[det->n_rigid_groups++] = rg; - - rg->name = strdup(name); - rg->panels = NULL; - rg->n_panels = 0; - - return rg; -} - - -static struct rg_collection *find_or_add_rg_coll(struct detector *det, - const char *name) -{ - int i; - struct rg_collection **new; - struct rg_collection *rgc; - - for ( i=0; i<det->n_rg_collections; i++ ) { - if ( strcmp(det->rigid_group_collections[i]->name, name) == 0 ) - { - return det->rigid_group_collections[i]; - } - } - - new = realloc(det->rigid_group_collections, - (1+det->n_rg_collections)*sizeof(struct rg_collection *)); - if ( new == NULL ) return NULL; - - det->rigid_group_collections = new; - - rgc = malloc(sizeof(struct rg_collection)); - if ( rgc == NULL ) return NULL; - - det->rigid_group_collections[det->n_rg_collections++] = rgc; - - rgc->name = strdup(name); - rgc->rigid_groups = NULL; - rgc->n_rigid_groups = 0; - - return rgc; -} - - -static void add_to_rigid_group(struct rigid_group *rg, struct panel *p) -{ - struct panel **pn; - - pn = realloc(rg->panels, (1+rg->n_panels)*sizeof(struct panel *)); - if ( pn == NULL ) { - ERROR("Couldn't add panel to rigid group.\n"); - return; - } - - rg->panels = pn; - rg->panels[rg->n_panels++] = p; -} - - -static void add_to_rigid_group_coll(struct rg_collection *rgc, - struct rigid_group *rg) -{ - struct rigid_group **r; - - r = realloc(rgc->rigid_groups, (1+rgc->n_rigid_groups)* - sizeof(struct rigid_group *)); - if ( r == NULL ) { - ERROR("Couldn't add rigid group to collection.\n"); - return; - } - - rgc->rigid_groups = r; - rgc->rigid_groups[rgc->n_rigid_groups++] = rg; -} - - -/* Free all rigid groups in detector */ -static void free_all_rigid_groups(struct detector *det) -{ - int i; - - if ( det->rigid_groups == NULL ) return; - for ( i=0; i<det->n_rigid_groups; i++ ) { - free(det->rigid_groups[i]->name); - free(det->rigid_groups[i]->panels); - free(det->rigid_groups[i]); - } - free(det->rigid_groups); -} - - -/* Free all rigid groups in detector */ -static void free_all_rigid_group_collections(struct detector *det) -{ - int i; - - if ( det->rigid_group_collections == NULL ) return; - for ( i=0; i<det->n_rg_collections; i++ ) { - free(det->rigid_group_collections[i]->name); - free(det->rigid_group_collections[i]->rigid_groups); - free(det->rigid_group_collections[i]); - } - free(det->rigid_group_collections); -} - - -static struct rigid_group *find_rigid_group_by_name(struct detector *det, - char *name) -{ - int i; - - for ( i=0; i<det->n_rigid_groups; i++ ) { - if ( strcmp(det->rigid_groups[i]->name, name) == 0 ) { - return det->rigid_groups[i]; - } - } - - return NULL; -} - - -static int parse_field_for_panel(struct panel *panel, const char *key, - const char *val, struct detector *det) -{ - int reject = 0; - - if ( strcmp(key, "min_fs") == 0 ) { - panel->orig_min_fs = atof(val); - } else if ( strcmp(key, "max_fs") == 0 ) { - panel->orig_max_fs = atof(val); - } else if ( strcmp(key, "min_ss") == 0 ) { - panel->orig_min_ss = atof(val); - } else if ( strcmp(key, "max_ss") == 0 ) { - panel->orig_max_ss = atof(val); - } else if ( strcmp(key, "corner_x") == 0 ) { - panel->cnx = atof(val); - } else if ( strcmp(key, "corner_y") == 0 ) { - panel->cny = atof(val); - } else if ( strcmp(key, "rail_direction") == 0 ) { - if ( dir_conv(val, &panel->rail_x, - &panel->rail_y, - &panel->rail_z) ) - { - ERROR("Invalid rail direction '%s'\n", val); - reject = 1; - } - } else if ( strcmp(key, "clen_for_centering") == 0 ) { - panel->clen_for_centering = atof(val); - } else if ( strcmp(key, "adu_per_eV") == 0 ) { - panel->adu_per_eV = atof(val); - } else if ( strcmp(key, "adu_per_photon") == 0 ) { - panel->adu_per_photon = atof(val); - } else if ( strcmp(key, "rigid_group") == 0 ) { - add_to_rigid_group(find_or_add_rg(det, val), panel); - } else if ( strcmp(key, "clen") == 0 ) { - - char *end; - double v = strtod(val, &end); - if ( end == val ) { - /* This means "fill in later" */ - panel->clen = -1.0; - panel->clen_from = strdup(val); - } else { - panel->clen = v; - panel->clen_from = NULL; - } - - } else if ( strcmp(key, "data") == 0 ) { - if ( strncmp(val,"/",1) != 0 ) { - ERROR("Invalid data location '%s'\n", val); - reject = -1; - } - panel->data = strdup(val); - - } else if ( strcmp(key, "mask") == 0 ) { - if ( strncmp(val,"/",1) != 0 ) { - ERROR("Invalid mask location '%s'\n", val); - reject = -1; - } - panel->mask = strdup(val); - - } else if ( strcmp(key, "mask_file") == 0 ) { - panel->mask_file = strdup(val); - - } else if ( strcmp(key, "saturation_map") == 0 ) { - panel->satmap = strdup(val); - } else if ( strcmp(key, "saturation_map_file") == 0 ) { - panel->satmap_file = strdup(val); - - } else if ( strcmp(key, "coffset") == 0) { - panel->coffset = atof(val); - } else if ( strcmp(key, "res") == 0 ) { - panel->res = atof(val); - } else if ( strcmp(key, "max_adu") == 0 ) { - panel->max_adu = atof(val); - } else if ( strcmp(key, "badrow_direction") == 0 ) { - panel->badrow = val[0]; /* First character only */ - if ( (panel->badrow != 'x') && (panel->badrow != 'y') - && (panel->badrow != 'f') && (panel->badrow != 's') - && (panel->badrow != '-') ) { - ERROR("badrow_direction must be x, y, f, s or '-'\n"); - ERROR("Assuming '-'\n."); - panel->badrow = '-'; - } - if ( panel->badrow == 'x' ) panel->badrow = 'f'; - if ( panel->badrow == 'y' ) panel->badrow = 's'; - } else if ( strcmp(key, "no_index") == 0 ) { - panel->no_index = atob(val); - } else if ( strcmp(key, "fs") == 0 ) { - if ( dir_conv(val, &panel->fsx, &panel->fsy, &panel->fsz) != 0 ) - { - ERROR("Invalid fast scan direction '%s'\n", val); - reject = 1; - } - } else if ( strcmp(key, "ss") == 0 ) { - if ( dir_conv(val, &panel->ssx, &panel->ssy, &panel->ssz) != 0 ) - { - ERROR("Invalid slow scan direction '%s'\n", val); - reject = 1; - } - } else if ( strncmp(key, "dim", 3) == 0) { - int dim_entry; - char *endptr; - if ( key[3] != '\0' ) { - if ( panel->dim_structure == NULL ) { - panel->dim_structure = initialize_dim_structure(); - } - dim_entry = strtoul(key+3, &endptr, 10); - if ( endptr[0] != '\0' ) { - ERROR("Invalid dimension number %s\n", key+3); - } else { - if ( set_dim_structure_entry(panel->dim_structure, - dim_entry, val) ) - { - ERROR("Failed to set dim structure entry\n"); - } - } - } else { - ERROR("'dim' must be followed by a number, e.g. 'dim0'\n"); - } - } else { - ERROR("Unrecognised field '%s'\n", key); - } - - return reject; -} - - -static int check_badr_fsss(struct badregion *badr, int is_fsss) -{ - /* First assignment? */ - if ( badr->is_fsss == 99 ) { - badr->is_fsss = is_fsss; - return 0; - } - - if ( is_fsss != badr->is_fsss ) { - ERROR("You can't mix x/y and fs/ss in a bad region.\n"); - return 1; - } - - return 0; -} - - -static int parse_field_bad(struct badregion *badr, const char *key, - const char *val) -{ - int reject = 0; - - if ( strcmp(key, "min_x") == 0 ) { - badr->min_x = atof(val); - reject = check_badr_fsss(badr, 0); - } else if ( strcmp(key, "max_x") == 0 ) { - badr->max_x = atof(val); - reject = check_badr_fsss(badr, 0); - } else if ( strcmp(key, "min_y") == 0 ) { - badr->min_y = atof(val); - reject = check_badr_fsss(badr, 0); - } else if ( strcmp(key, "max_y") == 0 ) { - badr->max_y = atof(val); - reject = check_badr_fsss(badr, 0); - } else if ( strcmp(key, "min_fs") == 0 ) { - badr->min_fs = atof(val); - reject = check_badr_fsss(badr, 1); - } else if ( strcmp(key, "max_fs") == 0 ) { - badr->max_fs = atof(val); - reject = check_badr_fsss(badr, 1); - } else if ( strcmp(key, "min_ss") == 0 ) { - badr->min_ss = atof(val); - reject = check_badr_fsss(badr, 1); - } else if ( strcmp(key, "max_ss") == 0 ) { - badr->max_ss = atof(val); - reject = check_badr_fsss(badr, 1); - } else if ( strcmp(key, "panel") == 0 ) { - badr->panel = strdup(val); - } else { - ERROR("Unrecognised field '%s'\n", key); - } - - return reject; -} - - -static void parse_toplevel(struct detector *det, struct beam_params *beam, - const char *key, const char *val, - struct rg_definition ***rg_defl, - struct rgc_definition ***rgc_defl, int *n_rg_defs, - int *n_rgc_defs, char **hdf5_peak_path) -{ - - if ( strcmp(key, "mask_bad") == 0 ) { - - char *end; - double v = strtod(val, &end); - - if ( end != val ) { - det->mask_bad = v; - } - - } else if ( strcmp(key, "mask_good") == 0 ) { - - char *end; - double v = strtod(val, &end); - - if ( end != val ) { - det->mask_good = v; - } - - } else if ( strcmp(key, "coffset") == 0 ) { - det->defaults.coffset = atof(val); - - } else if ( strcmp(key, "photon_energy") == 0 ) { - if ( beam != NULL ) { - double v; - char *end; - v = strtod(val, &end); - if ( (val[0] != '\0') && (end[0] == '\0') ) { - beam->photon_energy = v; - beam->photon_energy_from = NULL; - } else { - beam->photon_energy = 0.0; - beam->photon_energy_from = strdup(val); - } - } - - } else if ( strcmp(key, "photon_energy_bandwidth") == 0 ) { - if ( beam != NULL ) { - double v; - char *end; - v = strtod(val, &end); - if ( (val[0] != '\0') && (end[0] == '\0') ) { - beam->bandwidth = v; - } else { - ERROR("Invalid value for " - "photon_energy_bandwidth\n"); - } - } - - } else if ( strcmp(key, "photon_energy_scale") == 0 ) { - if ( beam != NULL ) { - beam->photon_energy_scale = atof(val); - } - - } else if ( strcmp(key, "peak_info_location") == 0 ) { - if ( hdf5_peak_path != NULL ) { - *hdf5_peak_path = strdup(val); - } - - } else if (strncmp(key, "rigid_group", 11) == 0 - && strncmp(key, "rigid_group_collection", 22) != 0 ) { - - struct rg_definition **new; - - new = realloc(*rg_defl, - ((*n_rg_defs)+1)*sizeof(struct rg_definition*)); - *rg_defl = new; - - (*rg_defl)[*n_rg_defs] = malloc(sizeof(struct rg_definition)); - (*rg_defl)[*n_rg_defs]->name = strdup(key+12); - (*rg_defl)[*n_rg_defs]->pns = strdup(val); - *n_rg_defs = *n_rg_defs+1; - - } else if ( strncmp(key, "rigid_group_collection", 22) == 0 ) { - - struct rgc_definition **new; - - new = realloc(*rgc_defl, ((*n_rgc_defs)+1)* - sizeof(struct rgc_definition*)); - *rgc_defl = new; - - (*rgc_defl)[*n_rgc_defs] = - malloc(sizeof(struct rgc_definition)); - (*rgc_defl)[*n_rgc_defs]->name = strdup(key+23); - (*rgc_defl)[*n_rgc_defs]->rgs = strdup(val); - *n_rgc_defs = *n_rgc_defs+1; - - } else if ( parse_field_for_panel(&det->defaults, key, val, det) ) { - ERROR("Unrecognised top level field '%s'\n", key); - } -} - - -/* Test if fs,ss in panel "p" is further {out,in} than {*p_max_d,*p_min_d}, and - * if so update det->furthest_{out,in}_{panel,fs,ss}. */ -static void check_point(struct panel *p, double fs, double ss, - double *p_min_d, double *p_max_d, struct detector *det) -{ - double xs, ys, rx, ry, d; - - xs = fs*p->fsx + ss*p->ssx; - ys = fs*p->fsy + ss*p->ssy; - - rx = (xs + p->cnx) / p->res; - ry = (ys + p->cny) / p->res; - - d = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); - - if ( d > *p_max_d ) { - - det->furthest_out_panel = p; - det->furthest_out_fs = fs; - det->furthest_out_ss = ss; - *p_max_d = d; - - } else if ( d < *p_min_d ) { - - det->furthest_in_panel = p; - det->furthest_in_fs = fs; - det->furthest_in_ss = ss; - *p_min_d = d; - - } -} - - -static void find_min_max_d(struct detector *det) -{ - double max_d, min_d; - int i; - - min_d = +INFINITY; - max_d = 0.0; - for ( i=0; i<det->n_panels; i++ ) { - - struct panel *p; - - p = &det->panels[i]; - - check_point(p, 0, 0, &min_d, &max_d, det); - check_point(p, p->w, 0, &min_d, &max_d, det); - check_point(p, 0, p->h, &min_d, &max_d, det); - check_point(p, p->w, p->h, &min_d, &max_d, det); - - } -} - - -struct detector *get_detector_geometry(const char *filename, - struct beam_params *beam) -{ - return get_detector_geometry_2(filename, beam, NULL); -} - - -struct detector *get_detector_geometry_from_string(const char *string_in, - struct beam_params *beam, - char **hdf5_peak_path) -{ - struct detector *det; - char **bits; - int done = 0; - int i; - int rgi, rgci; - int reject = 0; - int path_dim, mask_path_dim; - int dim_dim; - int dim_reject = 0; - int dim_dim_reject = 0; - struct rg_definition **rg_defl = NULL; - struct rgc_definition **rgc_defl = NULL; - int n_rg_definitions = 0; - int n_rgc_definitions = 0; - char *string; - char *string_orig; - size_t len; - - det = calloc(1, sizeof(struct detector)); - if ( det == NULL ) return NULL; - - if ( beam != NULL ) { - beam->photon_energy = 0.0; - beam->photon_energy_from = NULL; - beam->photon_energy_scale = 1.0; - beam->bandwidth = 0.00000001; - } - - det->n_panels = 0; - det->panels = NULL; - det->n_bad = 0; - det->bad = NULL; - det->mask_good = 0; - det->mask_bad = 0; - det->n_rigid_groups = 0; - det->rigid_groups = NULL; - det->path_dim = 0; - det->dim_dim = 0; - det->n_rg_collections = 0; - det->rigid_group_collections = NULL; - - /* The default defaults... */ - det->defaults.orig_min_fs = -1; - det->defaults.orig_min_ss = -1; - det->defaults.orig_max_fs = -1; - det->defaults.orig_max_ss = -1; - det->defaults.cnx = NAN; - det->defaults.cny = NAN; - det->defaults.clen = NAN; - det->defaults.coffset = 0.0; - det->defaults.res = -1.0; - det->defaults.badrow = '-'; - det->defaults.no_index = 0; - det->defaults.fsx = 1.0; - det->defaults.fsy = 0.0; - det->defaults.fsz = 0.0; - det->defaults.ssx = 0.0; - det->defaults.ssy = 1.0; - det->defaults.ssz = 0.0; - det->defaults.rail_x = NAN; /* The actual default rail direction */ - det->defaults.rail_y = NAN; /* is below */ - det->defaults.rail_z = NAN; - det->defaults.clen_for_centering = NAN; - det->defaults.adu_per_eV = NAN; - det->defaults.adu_per_photon = NAN; - det->defaults.max_adu = +INFINITY; - det->defaults.mask = NULL; - det->defaults.mask_file = NULL; - det->defaults.satmap = NULL; - det->defaults.satmap_file = NULL; - det->defaults.data = NULL; - det->defaults.dim_structure = NULL; - strncpy(det->defaults.name, "", 1023); - - string = strdup(string_in); - if ( string == NULL ) return NULL; - len = strlen(string); - for ( i=0; i<len; i++ ) { - if ( string_in[i] == '\r' ) string[i] = '\n'; - } - - /* Because 'string' will get modified */ - string_orig = string; - - do { - - int n1, n2; - char **path; - char *line; - struct badregion *badregion = NULL; - struct panel *panel = NULL; - char wholeval[1024]; - - const char *nl = strchr(string, '\n'); - if ( nl != NULL ) { - size_t len = nl - string; - line = strndup(string, nl-string); - line[len] = '\0'; - string += len+1; - } else { - line = strdup(string); - done = 1; - } - - if ( line[0] == ';' ) { - free(line); - continue; - } - - n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); - if ( n1 < 3 ) { - for ( i=0; i<n1; i++ ) free(bits[i]); - free(bits); - free(line); - continue; - } - - /* Stitch the pieces of the "value" back together */ - wholeval[0] = '\0'; /* Empty string */ - for ( i=2; i<n1; i++ ) { - if ( bits[i][0] == ';' ) break; /* Stop on comment */ - strncat(wholeval, bits[i], 1023); - } - - if ( bits[1][0] != '=' ) { - for ( i=0; i<n1; i++ ) free(bits[i]); - free(bits); - free(line); - continue; - } - - n2 = assplode(bits[0], "/\\.", &path, ASSPLODE_NONE); - if ( n2 < 2 ) { - - /* This was a top-level option, not handled above. */ - parse_toplevel(det, beam, bits[0], wholeval, &rg_defl, - &rgc_defl, &n_rg_definitions, - &n_rgc_definitions, hdf5_peak_path); - for ( i=0; i<n1; i++ ) free(bits[i]); - free(bits); - for ( i=0; i<n2; i++ ) free(path[i]); - free(path); - free(line); - continue; - } - - if ( strncmp(path[0], "bad", 3) == 0 ) { - badregion = find_bad_region_by_name(det, path[0]); - if ( badregion == NULL ) { - badregion = new_bad_region(det, path[0]); - } - } else { - panel = find_panel_by_name(det, path[0]); - if ( panel == NULL ) { - panel = new_panel(det, path[0]); - } - } - - if ( panel != NULL ) { - if ( parse_field_for_panel(panel, path[1], - wholeval, det) ) - { - reject = 1; - } - } else { - if ( parse_field_bad(badregion, path[1], wholeval) ) { - reject = 1; - } - } - - for ( i=0; i<n1; i++ ) free(bits[i]); - for ( i=0; i<n2; i++ ) free(path[i]); - free(bits); - free(path); - free(line); - - } while ( !done ); - - if ( det->n_panels == -1 ) { - ERROR("No panel descriptions in geometry file.\n"); - free(det); - return NULL; - } - - path_dim = -1; - dim_reject = 0; - - for ( i=0; i<det->n_panels; i++ ) { - - int panel_dim = 0; - char *next_instance; - - next_instance = det->panels[i].data; - - while ( next_instance ) { - next_instance = strstr(next_instance, "%"); - if ( next_instance != NULL ) { - next_instance += 1*sizeof(char); - panel_dim += 1; - } - } - - if ( path_dim == -1 ) { - path_dim = panel_dim; - } else { - if ( panel_dim != path_dim ) { - dim_reject = 1; - } - } - - } - - mask_path_dim = -1; - for ( i=0; i<det->n_panels; i++ ) { - - int panel_mask_dim = 0; - char *next_instance; - - if ( det->panels[i].mask != NULL ) { - - next_instance = det->panels[i].mask; - - while ( next_instance ) { - next_instance = strstr(next_instance, "%"); - if ( next_instance != NULL ) { - next_instance += 1*sizeof(char); - panel_mask_dim += 1; - } - } - - if ( mask_path_dim == -1 ) { - mask_path_dim = panel_mask_dim; - } else { - if ( panel_mask_dim != mask_path_dim ) { - dim_reject = 1; - } - } - - } - } - - if ( dim_reject == 1 ) { - ERROR("All panels' data and mask entries must have the same " - "number of placeholders\n"); - reject = 1; - } - - if ( mask_path_dim > path_dim ) { - ERROR("Number of placeholders in mask cannot be larger than " - "for data\n"); - reject = 1; - } - - det->path_dim = path_dim; - - dim_dim_reject = 0; - dim_dim = -1; - - for ( i=0; i<det->n_panels; i++ ) { - - int di; - int found_ss = 0; - int found_fs = 0; - int panel_dim_dim = 0; - - if ( det->panels[i].dim_structure == NULL ) { - det->panels[i].dim_structure = default_dim_structure(); - } - - for ( di=0; di<det->panels[i].dim_structure->num_dims; di++ ) { - - if ( det->panels[i].dim_structure->dims[di] == - HYSL_UNDEFINED ) { - dim_dim_reject = 1; - ERROR("Dimension %i for panel %s is undefined.\n", - di, det->panels[i].name); - } - if ( det->panels[i].dim_structure->dims[di] == - HYSL_PLACEHOLDER ) { - panel_dim_dim += 1; - } - if ( det->panels[i].dim_structure->dims[di] == - HYSL_SS ) { - found_ss += 1; - } - if ( det->panels[i].dim_structure->dims[di] == - HYSL_FS ) { - found_fs += 1; - } - - } - - if ( found_ss != 1 ) { - ERROR("Exactly one slow scan dim coordinate is needed " - "(found %i for panel %s)\n", found_ss, - det->panels[i].name); - dim_dim_reject = 1; - } - - if ( found_fs != 1 ) { - ERROR("Exactly one fast scan dim coordinate is needed " - "(found %i for panel %s)\n", found_fs, - det->panels[i].name); - dim_dim_reject = 1; - } - - if ( panel_dim_dim > 1 ) { - ERROR("Maximum one placeholder dim coordinate is " - "allowed (found %i for panel %s)\n", - panel_dim_dim, det->panels[i].name); - dim_dim_reject = 1; - } - - if ( dim_dim == -1 ) { - dim_dim = panel_dim_dim; - } else { - if ( panel_dim_dim != dim_dim ) { - dim_dim_reject = 1; - } - } - - } - - if ( dim_dim_reject == 1) { - reject = 1; - } - - det->dim_dim = dim_dim; - - for ( i=0; i<det->n_panels; i++ ) { - - struct panel *p = &det->panels[i]; - - if ( p->orig_min_fs < 0 ) { - ERROR("Please specify the minimum FS coordinate for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( p->orig_max_fs < 0 ) { - ERROR("Please specify the maximum FS coordinate for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( p->orig_min_ss < 0 ) { - ERROR("Please specify the minimum SS coordinate for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( p->orig_max_ss < 0 ) { - ERROR("Please specify the maximum SS coordinate for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( isnan(p->cnx) ) { - ERROR("Please specify the corner X coordinate for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( isnan(p->cny) ) { - ERROR("Please specify the corner Y coordinate for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( isnan(p->clen) && (p->clen_from == NULL) ) { - ERROR("Please specify the camera length for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( p->res < 0 ) { - ERROR("Please specify the resolution for" - " panel %s\n", det->panels[i].name); - reject = 1; - } - if ( isnan(p->adu_per_eV) && isnan(p->adu_per_photon) ) { - ERROR("Please specify either adu_per_eV or " - "adu_per_photon for panel %s\n", - det->panels[i].name); - reject = 1; - } - - if ( isnan(p->clen_for_centering) && !isnan(p->rail_x) ) - { - ERROR("You must specify clen_for_centering if you " - "specify the rail direction (panel %s)\n", - p->name); - reject = 1; - } - - if ( (p->mask_file != NULL) && (p->mask == NULL) ) { - ERROR("You have specified 'mask_file' but not 'mask'. " - "'mask_file' will therefore have no effect. " - "(panel %s)\n", p->name); - reject = 1; - } - - /* The default rail direction */ - if ( isnan(p->rail_x) ) { - p->rail_x = 0.0; - p->rail_y = 0.0; - p->rail_z = 1.0; - } - if ( isnan(p->clen_for_centering) ) p->clen_for_centering = 0.0; - - /* It's OK if the badrow direction is '0' */ - /* It's not a problem if "no_index" is still zero */ - /* The default transformation matrix is at least valid */ - - det->panels[i].w = det->panels[i].orig_max_fs - - det->panels[i].orig_min_fs+1; - det->panels[i].h = det->panels[i].orig_max_ss - - det->panels[i].orig_min_ss+1; - - } - - for ( i=0; i<det->n_bad; i++ ) { - if ( det->bad[i].is_fsss == 99 ) { - ERROR("Please specify the coordinate ranges for" - " bad region %s\n", det->bad[i].name); - reject = 1; - } - } - - free(det->defaults.clen_from); - free(det->defaults.data); - free(det->defaults.mask); - - for ( rgi=0; rgi<n_rg_definitions; rgi++) { - - int pi, n1; - struct rigid_group *rigidgroup = NULL; - - rigidgroup = find_or_add_rg(det, rg_defl[rgi]->name); - - n1 = assplode(rg_defl[rgi]->pns, ",", &bits, ASSPLODE_NONE); - - for ( pi=0; pi<n1; pi++ ) { - - struct panel *p; - - p = find_panel_by_name(det, bits[pi]); - if ( p == NULL ) { - ERROR("Cannot add panel to rigid group\n"); - ERROR("Panel not found: %s\n", bits[pi]); - return NULL; - } - add_to_rigid_group(rigidgroup, p); - free(bits[pi]); - } - free(bits); - free(rg_defl[rgi]->name); - free(rg_defl[rgi]->pns); - free(rg_defl[rgi]); - } - free(rg_defl); - - for ( rgci=0; rgci<n_rgc_definitions; rgci++ ) { - - int rgi, n2; - struct rg_collection *rgcollection = NULL; - - rgcollection = find_or_add_rg_coll(det, rgc_defl[rgci]->name); - - n2 = assplode(rgc_defl[rgci]->rgs, ",", &bits, ASSPLODE_NONE); - - for ( rgi=0; rgi<n2; rgi++ ) { - - struct rigid_group *r; - - r = find_rigid_group_by_name(det, bits[rgi]); - if ( r == NULL ) { - ERROR("Cannot add rigid group to collection\n"); - ERROR("Rigid group not found: %s\n", bits[rgi]); - return NULL; - } - add_to_rigid_group_coll(rgcollection, r); - free(bits[rgi]); - } - free(bits); - free(rgc_defl[rgci]->name); - free(rgc_defl[rgci]->rgs); - free(rgc_defl[rgci]); - - } - free(rgc_defl); - - if ( n_rg_definitions == 0 ) { - - int pi; - - for ( pi=0; pi<det->n_panels; pi++ ) { - - struct rigid_group *rigidgroup = NULL; - - rigidgroup = find_or_add_rg(det, det->panels[pi].name); - add_to_rigid_group(rigidgroup, &det->panels[pi]); - - } - } - - if ( n_rgc_definitions == 0 ) { - - int rgi; - struct rg_collection *rgcollection = NULL; - - rgcollection = find_or_add_rg_coll(det, "default"); - - for ( rgi=0; rgi<det->n_rigid_groups; rgi++ ) { - - add_to_rigid_group_coll(rgcollection, - det->rigid_groups[rgi]); - - } - } - - /* Calculate matrix inverses and other stuff */ - for ( i=0; i<det->n_panels; i++ ) { - - struct panel *p; - double d; - - p = &det->panels[i]; - - if ( p->fsx*p->ssy == p->ssx*p->fsy ) { - ERROR("Panel %i transformation singular.\n", i); - } - - d = (double)p->fsx*p->ssy - p->ssx*p->fsy; - p->xfs = p->ssy / d; - p->yfs = -p->ssx / d; - p->xss = -p->fsy / d; - p->yss = p->fsx / d; - - p->w = p->orig_max_fs - p->orig_min_fs + 1; - p->h = p->orig_max_ss - p->orig_min_ss + 1; - - } - - find_min_max_d(det); - free(string_orig); - - if ( reject ) return NULL; - - return det; -} - - -struct detector *get_detector_geometry_2(const char *filename, - struct beam_params *beam, - char **hdf5_peak_path) -{ - char *contents; - struct detector *det; - - contents = load_entire_file(filename); - if ( contents == NULL ) { - ERROR("Failed to load geometry file '%s'\n", filename); - return NULL; - } - - det = get_detector_geometry_from_string(contents, beam, hdf5_peak_path); - free(contents); - return det; -} - - -void free_detector_geometry(struct detector *det) -{ - int i; - - free_all_rigid_groups(det); - free_all_rigid_group_collections(det); - - for ( i=0; i<det->n_panels; i++ ) { - free(det->panels[i].clen_from); - free_dim_structure(det->panels[i].dim_structure); - } - - free(det->panels); - free(det->bad); - free(det); -} - - -static int rg_number(const struct detector *det, const struct rigid_group *rg) -{ - int i; - for ( i=0; i<det->n_rigid_groups; i++ ) { - if ( det->rigid_groups[i] == rg ) return i; - } - return det->n_rigid_groups; -} - - -struct detector *copy_geom(const struct detector *in) -{ - struct detector *out; - int i; - - if ( in == NULL ) return NULL; - - out = malloc(sizeof(struct detector)); - if ( out == NULL ) return NULL; - - /* Copy everything */ - memcpy(out, in, sizeof(struct detector)); - - out->panels = malloc(out->n_panels * sizeof(struct panel)); - memcpy(out->panels, in->panels, out->n_panels * sizeof(struct panel)); - - out->bad = malloc(out->n_bad * sizeof(struct badregion)); - memcpy(out->bad, in->bad, out->n_bad * sizeof(struct badregion)); - - /* Copy the panels */ - for ( i=0; i<out->n_panels; i++ ) { - - struct panel *p; - - /* Copy all fields */ - p = &out->panels[i]; - - /* Now fix up everything involving pointers... */ - - if ( p->clen_from != NULL ) { - /* Make a copy of the clen_from fields unique to this - * copy of the structure. */ - p->clen_from = strdup(p->clen_from); - } - - if ( p->data != NULL ) { - /* Make a copy of the data fields unique to this - * copy of the structure. */ - p->data = strdup(p->data); - } - - if ( p->dim_structure != NULL ) { - /* Make a copy of the dim_structure fields unique to this - * copy of the structure. */ - - struct dim_structure *dim_new; - int di; - - dim_new = initialize_dim_structure(); - dim_new->num_dims = p->dim_structure->num_dims; - dim_new->dims = malloc(dim_new->num_dims*sizeof(int)); - for ( di=0; di<dim_new->num_dims; di++ ) { - dim_new->dims[di] = p->dim_structure->dims[di]; - } - - p->dim_structure = dim_new; - - } - - if ( &in->panels[i] == in->furthest_out_panel ) { - out->furthest_out_panel = &out->panels[i]; - } - if ( &in->panels[i] == in->furthest_in_panel ) { - out->furthest_in_panel = &out->panels[i]; - } - - } - - /* Copy all the rigid groups */ - out->rigid_groups = malloc(out->n_rigid_groups*sizeof(struct rigid_group)); - if ( out->rigid_groups == NULL ) return NULL; - for ( i=0; i<out->n_rigid_groups; i++ ) { - - struct rigid_group *inrg; - struct rigid_group *rg; - int j; - - rg = malloc(sizeof(struct rigid_group)); - if ( rg == NULL ) return NULL; - - out->rigid_groups[i] = rg; - - inrg = in->rigid_groups[i]; - - rg->name = strdup(inrg->name); - if ( rg->name == NULL ) return NULL; - - rg->n_panels = inrg->n_panels; - rg->panels = malloc(inrg->n_panels*sizeof(struct panel *)); - if ( rg->panels == NULL ) return NULL; - - for ( j=0; j<rg->n_panels; j++ ) { - int k = panel_number(in, inrg->panels[j]); - rg->panels[j] = &out->panels[k]; - } - - } - - /* Copy all the rigid group collections */ - out->rigid_group_collections = malloc(out->n_rg_collections*sizeof(struct rg_collection)); - if ( out->rigid_group_collections == NULL ) return NULL; - for ( i=0; i<out->n_rg_collections; i++ ) { - - struct rg_collection *inrgc; - struct rg_collection *rgc; - int j; - - rgc = malloc(sizeof(struct rg_collection)); - if ( rgc == NULL ) return NULL; - - out->rigid_group_collections[i] = rgc; - - inrgc = in->rigid_group_collections[i]; - - rgc->name = strdup(inrgc->name); - if ( rgc->name == NULL ) return NULL; - - rgc->n_rigid_groups = inrgc->n_rigid_groups; - rgc->rigid_groups = malloc(rgc->n_rigid_groups*sizeof(struct rg_collection)); - if ( rgc->rigid_groups == NULL ) return NULL; - - for ( j=0; j<rgc->n_rigid_groups; j++ ) { - int k = rg_number(in, inrgc->rigid_groups[j]); - if ( k == in->n_rigid_groups ) return NULL; - rgc->rigid_groups[j] = out->rigid_groups[k]; - } - - } - - return out; -} - - -struct detector *simple_geometry(const struct image *image, int w, int h) -{ - struct detector *geom; - - geom = calloc(1, sizeof(struct detector)); - - geom->n_panels = 1; - geom->panels = calloc(1, sizeof(struct panel)); - - geom->panels[0].orig_min_fs = 0; - geom->panels[0].orig_max_fs = w-1; - geom->panels[0].orig_min_ss = 0; - geom->panels[0].orig_max_ss = h-1; - geom->panels[0].cnx = -w / 2.0; - geom->panels[0].cny = -h / 2.0; - geom->panels[0].max_adu = INFINITY; - geom->panels[0].orig_min_fs = -1; - geom->panels[0].orig_max_fs = -1; - geom->panels[0].orig_min_ss = -1; - geom->panels[0].orig_max_ss = -1; - - geom->panels[0].fsx = 1; - geom->panels[0].fsy = 0; - geom->panels[0].fsz = 0; - geom->panels[0].ssx = 0; - geom->panels[0].ssy = 1; - geom->panels[0].ssz = 0; - - geom->panels[0].xfs = 1; - geom->panels[0].xss = 0; - geom->panels[0].yfs = 0; - geom->panels[0].yss = 1; - - geom->panels[0].w = w; - geom->panels[0].h = h; - - geom->panels[0].mask = NULL; - geom->panels[0].data = NULL; - - find_min_max_d(geom); - - return geom; -} - - -int reverse_2d_mapping(double x, double y, struct detector *det, - struct panel **pp, double *pfs, double *pss) -{ - int i; - - for ( i=0; i<det->n_panels; i++ ) { - - struct panel *p = &det->panels[i]; - double cx, cy, fs, ss; - - /* Get position relative to corner */ - cx = x - p->cnx; - cy = y - p->cny; - - /* Reverse the transformation matrix */ - fs = cx*p->xfs + cy*p->yfs; - ss = cx*p->xss + cy*p->yss; - - /* In range? */ - if ( fs < 0 ) continue; - if ( ss < 0 ) continue; - if ( fs > p->w ) continue; - if ( ss > p->h ) continue; - - *pfs = fs; - *pss = ss; - *pp = p; - return 0; - - } - - return 1; -} - - -static void check_extents(struct panel p, double *min_x, double *min_y, - double *max_x, double *max_y, double fs, double ss) -{ - double xs, ys, rx, ry; - - xs = fs*p.fsx + ss*p.ssx; - ys = fs*p.fsy + ss*p.ssy; - - rx = xs + p.cnx; - ry = ys + p.cny; - - if ( rx > *max_x ) *max_x = rx; - if ( ry > *max_y ) *max_y = ry; - if ( rx < *min_x ) *min_x = rx; - if ( ry < *min_y ) *min_y = ry; -} - - -static void rewrite_panel_fields(const struct panel *p, char *line, - FILE *fh, char **bits, - int write_panel_coffset) -{ - char new_line[1024]; - char string_to_write[512]; - - strcpy(new_line,"\0"); - strcpy(string_to_write,"\0"); - - if(strstr(bits[1], "fs") != NULL && - strstr(bits[1], "min_fs") == NULL && - strstr(bits[1], "max_fs") == NULL && - strstr(bits[1], "offset") == NULL ) { - - sprintf(string_to_write, "%+fx %+fy", - p->fsx, p->fsy); - build_output_line(line, new_line, - string_to_write); - fputs(new_line, fh); - fputs("\n", fh); - return; - - } else if ( strstr(bits[1], "ss") != NULL && - strstr(bits[1], "min_ss") == NULL && - strstr(bits[1], "max_ss") == NULL) { - - sprintf(string_to_write, "%+fx %+fy", - p->ssx, p->ssy); - build_output_line(line, new_line, - string_to_write); - fputs(new_line, fh); - fputs("\n", fh); - return; - - } else if ( strstr(bits[1], "corner_x") != NULL) { - - sprintf(string_to_write, "%g", - p->cnx); - build_output_line(line, new_line, - string_to_write); - fputs(new_line, fh); - fputs("\n", fh); - return; - - } else if ( strstr(bits[1], "corner_y") != NULL) { - - sprintf(string_to_write, "%g", - p->cny); - build_output_line(line, new_line, - string_to_write); - fputs(new_line, fh); - fputs("\n", fh); - return; - - } else if ( strstr(bits[1], "coffset") != NULL) { - - if ( write_panel_coffset ) { - return; - } else { - fputs(line, fh); - fputs("\n", fh); - return; - } - - } else { - fputs(line, fh); - fputs("\n", fh); - } -} - - -double largest_q(struct image *image) -{ - struct rvec q; - double tt; - - if ( image->det == NULL ) { - ERROR("No detector geometry. assuming detector is infinite!\n"); - return INFINITY; - } - - q = get_q_for_panel(image->det->furthest_out_panel, - image->det->furthest_out_fs, - image->det->furthest_out_ss, - &tt, 1.0/image->lambda); - - return modulus(q.u, q.v, q.w); -} - - -double smallest_q(struct image *image) -{ - struct rvec q; - double tt; - - q = get_q_for_panel(image->det->furthest_in_panel, - image->det->furthest_in_fs, - image->det->furthest_in_ss, - &tt, 1.0/image->lambda); - - return modulus(q.u, q.v, q.w); -} - - -void get_pixel_extents(struct detector *det, - double *min_x, double *min_y, - double *max_x, double *max_y) -{ - int i; - - *min_x = 0.0; - *max_x = 0.0; - *min_y = 0.0; - *max_y = 0.0; - - /* To determine the maximum extents of the detector, put all four - * corners of each panel through the transformations and watch for the - * biggest */ - - for ( i=0; i<det->n_panels; i++ ) { - - check_extents(det->panels[i], min_x, min_y, max_x, max_y, - 0.0, 0.0); - - check_extents(det->panels[i], min_x, min_y, max_x, max_y, - 0.0, det->panels[i].h+1); - - check_extents(det->panels[i], min_x, min_y, max_x, max_y, - det->panels[i].w+1, 0.0); - - check_extents(det->panels[i], min_x, min_y, max_x, max_y, - det->panels[i].w+1, det->panels[i].h+1); - - - } -} - - -int write_detector_geometry_3(const char *geometry_data, - const char *output_filename, struct detector *det, - const char *additional_comment, - int write_panel_coffset) -{ - FILE *fh; - int done = 0; - - if ( geometry_data == NULL ) return 2; - if ( output_filename == NULL ) return 2; - if ( det->n_panels < 1 ) return 3; - - fh = fopen(output_filename, "w"); - if ( fh == NULL ) return 1; - - if ( additional_comment != NULL ) { - fputs("; ", fh); - fputs(additional_comment, fh); - fputs("\n", fh); - } - - if ( write_panel_coffset ) { - fputs("; Optimized panel offsets can be found at the " - "end of the file\n", fh); - } - - do { - - int n_bits; - char **bits; - int i; - struct panel *p; - char *line; - const char *nl; - - /* Get the next line */ - nl = strchr(geometry_data, '\n'); - if ( nl != NULL ) { - size_t len = nl - geometry_data; - line = strndup(geometry_data, nl-geometry_data); - line[len] = '\0'; - geometry_data += len+1; - } else { - /* Last line might now have newline at end */ - line = strdup(geometry_data); - done = 1; - } - - n_bits = assplode(line, "/=", &bits, ASSPLODE_NONE); - - if ( n_bits != 3 ) { - if ( write_panel_coffset && (bits != NULL) - && (strstr(bits[0], "coffset" ) != NULL) ) continue; - fputs(line, fh); - fputs("\n", fh); - } else { - - p = find_panel_by_name(det, bits[0]); - - if ( p != NULL ) { - rewrite_panel_fields(p, line, fh, bits, - write_panel_coffset); - } else { - fputs(line, fh); - fputs("\n", fh); - } - } - - for ( i=0; i<n_bits; i++ ) free(bits[i]); - - } while ( !done ); - - if ( write_panel_coffset ) { - - int pi; - - fputs("\n", fh); - - for ( pi=0; pi<det->n_panels; pi++ ) { - fprintf(fh, "%s/coffset = %f\n", - det->panels[pi].name, det->panels[pi].coffset); - } - } - - return 0; -} - - -int write_detector_geometry_2(const char *geometry_filename, - const char *output_filename, struct detector *det, - const char *additional_comment, - int write_panel_coffset) -{ - int r; - char *geometry_data = load_entire_file(geometry_filename); - r = write_detector_geometry_3(geometry_data, output_filename, det, - additional_comment, write_panel_coffset); - free(geometry_data); - return r; -} - - -int write_detector_geometry(const char *geometry_filename, - const char *output_filename, struct detector *det) -{ - return write_detector_geometry_2(geometry_filename, output_filename, - det, NULL, 0); -} - - -/** - * \param image An image structure - * \param min Minimum value of 1/d to be marked as bad - * \param max Maximum value of 1/d to be marked as bad - * - * Flags, in the bad pixel mask for \p image, every pixel whose resolution is - * between \p min and \p max. - * - */ - -void mark_resolution_range_as_bad(struct image *image, - double min, double max) -{ - int i; - - if ( isinf(min) && isinf(max) ) return; /* nothing to do */ - - for ( i=0; i<image->det->n_panels; i++ ) { - - int fs, ss; - struct panel *p = &image->det->panels[i]; - - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - struct rvec q; - double r; - q = get_q_for_panel(p, fs, ss, NULL, 1.0/image->lambda); - r = modulus(q.u, q.v, q.w); - if ( (r >= min) && (r <= max) ) { - image->bad[i][fs+p->w*ss] = 1; - } - } - } - - } -} - - -static int safe_strcmp(const char *a, const char *b) -{ - /* If both are NULL, they count as equal */ - if ( (a == NULL) && (b == NULL) ) return 0; - - /* Otherwise, if either is NULL then they're different */ - if ( a == NULL ) return 1; - if ( b == NULL ) return 1; - - /* Otherwise, normal string comparison */ - return strcmp(a, b); -} - - -/** - * \param det A detector structure - * \param element If manually selected by the user, the HDF5 element being used. - * Otherwise NULL. - * - * \returns Non-zero if the combination of \p det and \p element mean that all the - * data comes from a single block. - */ -int single_panel_data_source(struct detector *det, const char *element) -{ - int pi; - const char *first_datafrom = NULL; - const char *curr_datafrom = NULL; - - if ( det->panels[0].data == NULL ) { - first_datafrom = element; /* Might be NULL */ - } else { - first_datafrom = det->panels[0].data; - } - - for ( pi=1; pi<det->n_panels; pi++ ) { - - if ( det->panels[pi].data == NULL ) { - curr_datafrom = element; /* Might be NULL */ - } else { - curr_datafrom = det->panels[pi].data; - } - - if ( safe_strcmp(curr_datafrom, first_datafrom) != 0 ) { - return 0; - } - - } - - return 1; -} - - -int multi_event_geometry(struct detector *det) -{ - return (det->path_dim != 0) || (det->dim_dim != 0); -} |