aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2020-06-03 17:16:59 +0200
committerThomas White <taw@physics.org>2020-07-29 18:53:33 +0200
commit6e5d500a2786131c7304b6c001142e524791bc12 (patch)
tree0225c8f3dea1545332ef9294693ffc6b3263097b
parent2037afa4693d4032606d98eb67522fabb58a42ea (diff)
Remove "struct detector" completely, part I
record_image has been moved to pattern_sim.c
-rw-r--r--libcrystfel/CMakeLists.txt2
-rw-r--r--libcrystfel/src/detector.c2362
-rw-r--r--libcrystfel/src/detector.h333
-rw-r--r--libcrystfel/src/image.h23
-rw-r--r--src/pattern_sim.c126
5 files changed, 127 insertions, 2719 deletions
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index 32bf6fe3..b3953dbe 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -33,7 +33,6 @@ set(LIBCRYSTFEL_SOURCES
src/reflist.c
src/utils.c
src/cell.c
- src/detector.c
src/thread-pool.c
src/image.c
src/hdf5-file.c
@@ -83,7 +82,6 @@ set(LIBCRYSTFEL_HEADERS
src/reflist-utils.h
src/thread-pool.h
src/utils.h
- src/detector.h
src/geometry.h
src/peakfinder8.h
src/peaks.h
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);
-}
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
deleted file mode 100644
index 2176861f..00000000
--- a/libcrystfel/src/detector.h
+++ /dev/null
@@ -1,333 +0,0 @@
-/*
- * detector.h
- *
- * 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>
- * 2011-2012 Richard Kirian <rkirian@asu.edu>
- * 2014 Valerio Mariani
- * 2011 Andrew Aquila
- *
- * 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/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DETECTOR_H
-#define DETECTOR_H
-
-struct rigid_group;
-struct rg_collection;
-struct detector;
-struct panel;
-struct badregion;
-struct beam_params;
-struct hdfile;
-struct event;
-
-#include "hdf5-file.h"
-#include "image.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/**
- * \file detector.h
- * Detector geometry structure and related functions.
- */
-
-
-struct rigid_group
-{
- char *name;
- struct panel **panels;
- int n_panels;
-};
-
-
-struct rg_collection
-{
- char *name;
- struct rigid_group **rigid_groups;
- int n_rigid_groups;
-};
-
-
-/**
- * Represents one panel of a detector
- */
-struct panel
-{
- /** Text name for panel (fixed length array) */
- char name[1024];
-
- /** \name Location of corner in units of the pixel size of this panel */
- /**@{*/
- double cnx;
- double cny;
- /**@}*/
-
- /** The offset to be applied from \ref clen */
- double coffset;
-
- /** The distance from the interaction point to the corner of the
- * first pixel */
- double clen;
-
- /** Location to get \ref clen from, e.g. from HDF5 file */
- char *clen_from;
-
- /** Location of mask data */
- char *mask;
-
- /** Filename for mask data */
- char *mask_file;
-
- /** Location of per-pixel saturation map */
- char *satmap;
-
- /** Filename for saturation map */
- char *satmap_file;
-
- /** Resolution in pixels per metre */
- double res;
-
- /** Readout direction (for filtering out clusters of peaks)
- * ('x' or 'y') */
- char badrow;
-
- /** Non-zero if panel should be considered entirely bad */
- int no_index;
-
- /** Number of detector intensity units per photon */
- double adu_per_photon;
-
- /** Treat pixel as unreliable if higher than this */
- double max_adu;
-
- /** Location of data in file */
- char *data;
-
- /** Number of detector intensity units per eV of photon energy */
- double adu_per_eV;
-
- /** Dimension structure */
- struct dim_structure *dim_structure;
-
- /** \name Transformation matrix from pixel coordinates to lab frame */
- /*@{*/
- double fsx;
- double fsy;
- double fsz;
- double ssx;
- double ssy;
- double ssz;
- /*@}*/
-
- /** \name Rail direction */
- /*@{*/
- double rail_x;
- double rail_y;
- double rail_z;
- /*@}*/
-
- /* Value of clen (without coffset) at which beam is centered */
- double clen_for_centering;
-
- /** \name Inverse of 2D part of transformation matrix */
- /*@{*/
- double xfs;
- double yfs;
- double xss;
- double yss;
- /*@}*/
-
- /** \name Position of the panel in the data block in the file.
- * The panels may get moved around when the file is loaded (see
- * hdf5_read2()), especially if the panels come from different HDF5
- * elements. */
- /*@{*/
- int orig_min_fs;
- int orig_max_fs;
- int orig_min_ss;
- int orig_max_ss;
- /*@}*/
-
- /** Width, calculated as max_fs-min_fs+1 */
- int w;
-
- /*** Height, calculated as max_ss-min_ss+1 */
- int h;
-};
-
-
-struct badregion
-{
- char name[1024];
- int is_fsss;
- char *panel;
-
- double min_x;
- double max_x;
- double min_y;
- double max_y;
-
- /* Specified INCLUSIVELY */
- int min_fs;
- int max_fs;
- int min_ss;
- int max_ss;
-
-};
-
-
-struct detector
-{
- struct panel *panels;
- int n_panels;
-
- struct badregion *bad;
- int n_bad;
-
- unsigned int mask_bad;
- unsigned int mask_good;
-
- struct rigid_group **rigid_groups;
- int n_rigid_groups;
-
- struct rg_collection **rigid_group_collections;
- int n_rg_collections;
-
- /* Location of the pixel furthest away from the beam position, which
- * will have the largest value of 2theta regardless of camera length
- * and wavelength */
- struct panel *furthest_out_panel;
- double furthest_out_fs;
- double furthest_out_ss;
-
- /* As above, but for the smallest 2theta */
- struct panel *furthest_in_panel;
- double furthest_in_fs;
- double furthest_in_ss;
-
- int path_dim;
- int dim_dim;
-
- struct panel defaults;
-};
-
-
-extern struct rvec get_q_for_panel(struct panel *p, double fs, double ss,
- double *ttp, double k);
-
-extern double get_tt(struct image *image, double xs, double ys, int *err);
-
-extern int in_bad_region(struct detector *det, struct panel *p,
- double fs, double ss);
-
-extern void record_image(struct image *image, int do_poisson, double background,
- gsl_rng *rng, double beam_radius, double nphotons);
-
-extern struct panel *find_orig_panel(struct detector *det,
- double fs, double ss);
-
-extern signed int find_orig_panel_number(struct detector *det,
- double fs, double ss);
-
-extern int panel_number(const struct detector *det, const struct panel *p);
-
-extern struct detector *get_detector_geometry(const char *filename,
- struct beam_params *beam);
-
-extern struct detector *get_detector_geometry_2(const char *filename,
- struct beam_params *beam,
- char **hdf5_peak_path);
-
-extern struct detector *get_detector_geometry_from_string(const char *string,
- struct beam_params *beam,
- char **hdf5_peak_path);
-
-extern void free_detector_geometry(struct detector *det);
-
-extern struct detector *simple_geometry(const struct image *image, int w, int h);
-
-extern void get_pixel_extents(struct detector *det,
- double *min_x, double *min_y,
- double *max_x, double *max_y);
-
-extern void fill_in_adu(struct image *image);
-extern void adjust_centering_for_rail(struct panel *p);
-
-extern int panel_is_in_rigid_group(const struct rigid_group *rg,
- struct panel *p);
-
-extern int rigid_group_is_in_collection(struct rg_collection *c,
- struct rigid_group *rg);
-
-extern struct detector *copy_geom(const struct detector *in);
-
-extern int reverse_2d_mapping(double x, double y, struct detector *det,
- struct panel **pp, double *pfs, double *pss);
-
-extern double largest_q(struct image *image);
-
-extern double smallest_q(struct image *image);
-
-extern struct panel *find_panel_by_name(struct detector *det, const char *name);
-
-extern int write_detector_geometry_2(const char *geometry_filename,
- const char *output_filename,
- struct detector *det,
- const char *additional_comment,
- int write_panel_coffset);
-
-extern int write_detector_geometry_3(const char *geometry_data,
- const char *output_filename,
- struct detector *det,
- const char *additional_comment,
- int write_panel_coffset);
-
-extern int write_detector_geometry(const char *geometry_filename,
- const char *output_filename,
- struct detector *det);
-
-extern void mark_resolution_range_as_bad(struct image *image,
- double min, double max);
-
-
-extern int single_panel_data_source(struct detector *det, const char *element);
-
-struct rg_collection *find_rigid_group_collection_by_name(struct detector *det,
- const char *name);
-
-extern int detector_has_clen_references(struct detector *det);
-
-extern int multi_event_geometry(struct detector *det);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif /* DETECTOR_H */
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 44ce24ac..0235e0c2 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -35,19 +35,15 @@
#ifndef IMAGE_H
#define IMAGE_H
-struct detector;
-
#include <stdint.h>
#include <complex.h>
#include <sys/types.h>
struct imagefeature;
-struct sample;
struct image;
#include "utils.h"
#include "cell.h"
-#include "detector.h"
#include "reflist.h"
#include "crystal.h"
#include "index.h"
@@ -101,17 +97,6 @@ enum imagefile_type
typedef struct _imagefeaturelist ImageFeatureList;
-struct beam_params
-{
- double photon_energy; /**< eV per photon */
- char *photon_energy_from; /**< HDF5 dataset name */
- double photon_energy_scale; /**< Scale factor for photon energy, if it
- * comes from an image header */
- double bandwidth; /**< FWHM bandwidth as a fraction of
- * wavelength */
-};
-
-
struct image
{
/** The image data, by panel */
@@ -138,14 +123,8 @@ struct image
/** Number of times the indexer was tried before succeeding */
int n_indexing_tries;
- /** The detector structure
- * @{ */
+ /** The detector structure */
struct detgeom *detgeom;
- struct detector *det; /* FIXME: Deprecated */
- /** @} */
-
- /** The nominal beam parameters (or where to get them) */
- struct beam_params *beam; /* FIXME: Deprecated */
/** \name The filename and event ID for the image
* @{ */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 1fce0f60..9a144a36 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -116,6 +116,132 @@ static void show_help(const char *s)
}
+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);
+ }
+}
+
+
static double *intensities_from_list(RefList *list, SymOpList *sym)
{
Reflection *refl;