Commit 29e67fab authored by Thomas White's avatar Thomas White
Browse files

Get detector geometry from file

parent 1915044f
......@@ -4,6 +4,6 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h \
src/statistics.h src/displaywindow.h src/render.h src/hdfsee.h \
data/displaywindow.ui src/dirax.h src/peaks.h src/index.h \
src/filters.h src/diffraction-gpu.h src/cl-utils.h \
data/defs.h src/parameters-lcls.tmp src/geometry-lcls.tmp \
data/defs.h src/parameters-lcls.tmp \
data/diffraction.cl data/sfac src/likelihood.h
SUBDIRS = src data
......@@ -201,7 +201,7 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h \
src/statistics.h src/displaywindow.h src/render.h src/hdfsee.h \
data/displaywindow.ui src/dirax.h src/peaks.h src/index.h \
src/filters.h src/diffraction-gpu.h src/cl-utils.h \
data/defs.h src/parameters-lcls.tmp src/geometry-lcls.tmp \
data/defs.h src/parameters-lcls.tmp \
data/diffraction.cl data/sfac src/likelihood.h
SUBDIRS = src data
......
......@@ -33,7 +33,7 @@ int map_position(struct image *image, double dx, double dy,
double x = 0.0;
double y = 0.0;
p = find_panel(&image->det, dx, dy);
p = find_panel(image->det, dx, dy);
if ( p == NULL ) return 1;
x = ((double)dx - p->cx);
......@@ -94,7 +94,7 @@ void record_image(struct image *image, int do_poisson)
ERROR("NaN at %i,%i\n", x, y);
}
p = find_panel(&image->det, x, y);
p = find_panel(image->det, x, y);
/* Area of one pixel */
pix_area = pow(1.0/p->res, 2.0);
......@@ -158,3 +158,98 @@ struct panel *find_panel(struct detector *det, int x, int y)
return NULL;
}
struct detector *get_detector_geometry(const char *filename)
{
FILE *fh;
struct detector *det;
char *rval;
char **bits;
int n_panels = -1;
int i;
fh = fopen(filename, "r");
if ( fh == NULL ) return NULL;
det = malloc(sizeof(struct detector));
if ( det == NULL ) {
fclose(fh);
return NULL;
}
do {
int n1, n2;
char **path;
char line[1024];
int np;
rval = fgets(line, 1023, fh);
if ( rval == NULL ) break;
chomp(line);
n1 = assplode(line, " \t", &bits, ASSPLODE_NONE);
if ( n1 < 3 ) continue;
if ( bits[1][0] != '=' ) continue;
if ( strcmp(bits[0], "n_panels") == 0 ) {
if ( n_panels != -1 ) {
ERROR("Duplicate n_panels statement.\n");
fclose(fh);
free(det);
return NULL;
}
n_panels = atoi(bits[2]);
det->panels = malloc(n_panels * sizeof(struct panel));
continue;
}
n2 = assplode(bits[0], "/\\.", &path, ASSPLODE_NONE);
np = atoi(path[0]);
if ( strcmp(path[1], "min_x") == 0 ) {
det->panels[np].min_x = atof(bits[2]);
} else if ( strcmp(path[1], "max_x") == 0 ) {
det->panels[np].max_x = atof(bits[2]);
} else if ( strcmp(path[1], "min_y") == 0 ) {
det->panels[np].min_y = atof(bits[2]);
} else if ( strcmp(path[1], "max_y") == 0 ) {
det->panels[np].max_y = atof(bits[2]);
} else if ( strcmp(path[1], "cx") == 0 ) {
det->panels[np].cx = atof(bits[2]);
} else if ( strcmp(path[1], "cy") == 0 ) {
det->panels[np].cy = atof(bits[2]);
} else if ( strcmp(path[1], "clen") == 0 ) {
det->panels[np].clen = atof(bits[2]);
} else if ( strcmp(path[1], "res") == 0 ) {
det->panels[np].res = atof(bits[2]);
} else {
ERROR("Unrecognised field '%s'\n", path[1]);
}
} while ( rval != NULL );
if ( n_panels == -1 ) {
ERROR("No panel descriptions in geometry file.\n");
fclose(fh);
free(det->panels);
free(det);
return NULL;
}
for ( i=0; i<n_panels; i++ ) {
STATUS("Panel %i, min_x = %i\n", i, det->panels[i].min_x);
STATUS("Panel %i, max_x = %i\n", i, det->panels[i].max_x);
STATUS("Panel %i, min_y = %i\n", i, det->panels[i].min_y);
STATUS("Panel %i, max_y = %i\n", i, det->panels[i].max_y);
STATUS("Panel %i, cx = %f\n", i, det->panels[i].cx);
STATUS("Panel %i, cy = %f\n", i, det->panels[i].cy);
STATUS("Panel %i, clen = %f\n", i, det->panels[i].clen);
STATUS("Panel %i, res = %f\n", i, det->panels[i].res);
}
return det;
}
......@@ -47,4 +47,6 @@ extern void record_image(struct image *image, int do_poisson);
extern struct panel *find_panel(struct detector *det, int x, int y);
extern struct detector *get_detector_geometry(const char *filename);
#endif /* DETECTOR_H */
......@@ -348,7 +348,7 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
const unsigned int x = xs / sampling;
const unsigned int y = ys / sampling; /* Integer part only */
p = find_panel(&image->det, x, y);
p = find_panel(image->det, x, y);
assert(p != NULL);
rx = ((float)xs - (sampling*p->cx)) / (sampling * p->res);
......@@ -377,7 +377,7 @@ double get_tt(struct image *image, unsigned int xs, unsigned int ys)
const unsigned int x = xs;
const unsigned int y = ys; /* Integer part only */
p = find_panel(&image->det, x, y);
p = find_panel(image->det, x, y);
rx = ((float)xs - p->cx) / p->res;
ry = ((float)ys - p->cy) / p->res;
......
/* Set up detector configuration */
image.det.n_panels = 2;
image.det.panels = malloc(image.det.n_panels*sizeof(struct panel));
/* Upper panel */
image.det.panels[0].min_x = 0;
image.det.panels[0].max_x = 1023;
image.det.panels[0].min_y = 512;
image.det.panels[0].max_y = 1023;
image.det.panels[0].cx = 491.9;
image.det.panels[0].cy = 440.7;
image.det.panels[0].clen = 67.8e-3;
image.det.panels[0].res = 13333.3; /* 75 micron pixel size */
/* Lower panel */
image.det.panels[1].min_x = 0;
image.det.panels[1].max_x = 1023;
image.det.panels[1].min_y = 0;
image.det.panels[1].max_y = 511;
image.det.panels[1].cx = 492.0;
image.det.panels[1].cy = 779.7;
image.det.panels[1].clen = 70.8e-3;
image.det.panels[1].res = 13333.3; /* 75 micron pixel size */
/* Set up detector configuration */
image.det.n_panels = 2;
image.det.panels = malloc(image.det.n_panels*sizeof(struct panel));
/* Upper panel */
image.det.panels[0].min_x = 0;
image.det.panels[0].max_x = 1023;
image.det.panels[0].min_y = 512;
image.det.panels[0].max_y = 1023;
image.det.panels[0].cx = 512.0;
image.det.panels[0].cy = 502.0;
image.det.panels[0].clen = 50.0e-3;
image.det.panels[0].res = 13333.3; /* 75 micron pixel size */
/* Lower panel */
image.det.panels[1].min_x = 0;
image.det.panels[1].max_x = 1023;
image.det.panels[1].min_y = 0;
image.det.panels[1].max_y = 511;
image.det.panels[1].cx = 512.0;
image.det.panels[1].cy = 522.0;
image.det.panels[1].clen = 50.0e-3;
image.det.panels[1].res = 13333.3; /* 75 micron pixel size */
......@@ -86,7 +86,7 @@ struct image {
UnitCell *indexed_cell;
UnitCell *candidate_cells[MAX_CELL_CANDIDATES];
int ncells;
struct detector det;
struct detector *det;
const char *filename;
struct reflhit *hits;
int n_hits;
......
......@@ -59,6 +59,7 @@ struct process_args
int config_nomatch;
int config_unpolar;
int config_sanity;
struct detector *det;
IndexingMethod indm;
const double *intensities;
const unsigned int *counts;
......@@ -86,6 +87,7 @@ static void show_help(const char *s)
" --indexing=<method> Use 'method' for indexing. Choose from:\n"
" none : no indexing\n"
" dirax : invoke DirAx\n"
" -g. --geometry=<file> Get detector geometry from file.\n"
"\n\nWith just the above options, this program does not do much of practical "
"use.\nYou should also enable some of the following:\n\n"
" --near-bragg Output a list of reflection intensities to stdout.\n"
......@@ -159,7 +161,7 @@ static struct image *get_simage(struct image *template, int alternate)
* - not necessarily the same as the original. */
image->width = 1024;
image->height = 1024;
image->det.n_panels = 2;
image->det->n_panels = 2;
if ( alternate ) {
......@@ -183,12 +185,12 @@ static struct image *get_simage(struct image *template, int alternate)
panels[1].clen = 56.7e-2; /* 56.7 cm */
panels[1].res = 13333.3; /* 75 microns/pixel */
image->det.panels = panels;
image->det->panels = panels;
} else {
/* Copy pointer to old geometry */
image->det.panels = template->det.panels;
image->det->panels = template->det->panels;
}
......@@ -263,6 +265,7 @@ static void *process_image(void *pargsv)
image.filename = filename;
image.hits = NULL;
image.n_hits = 0;
image.det = pargs->det;
/* View head-on (unit cell is tilted) */
image.orientation.w = 1.0;
......@@ -285,8 +288,6 @@ static void *process_image(void *pargsv)
return result;
}
#include "geometry-lcls.tmp"
hdf5_read(hdfile, &image);
if ( config_cmfilter ) {
......@@ -372,7 +373,6 @@ static void *process_image(void *pargsv)
done:
free(image.data);
free(image.flags);
free(image.det.panels);
image_feature_list_free(image.features);
free(image.hits);
hdfile_close(hdfile);
......@@ -409,6 +409,8 @@ int main(int argc, char *argv[])
int config_alternate = 0;
int config_unpolar = 0;
int config_sanity = 0;
struct detector *det;
char *geometry = NULL;
IndexingMethod indm;
char *indm_str = NULL;
UnitCell *cell;
......@@ -435,6 +437,7 @@ int main(int argc, char *argv[])
{"near-bragg", 0, &config_nearbragg, 1},
{"write-drx", 0, &config_writedrx, 1},
{"indexing", 1, NULL, 'z'},
{"geometry", 1, NULL, 'g'},
{"simulate", 0, &config_simulate, 1},
{"filter-cm", 0, &config_cmfilter, 1},
{"filter-noise", 0, &config_noisefilter, 1},
......@@ -450,7 +453,8 @@ int main(int argc, char *argv[])
};
/* Short options */
while ((c = getopt_long(argc, argv, "hi:wp:j:x:", longopts, NULL)) != -1) {
while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:",
longopts, NULL)) != -1) {
switch (c) {
case 'h' : {
......@@ -488,6 +492,11 @@ int main(int argc, char *argv[])
break;
}
case 'g' : {
geometry = strdup(optarg);
break;
}
case 0 : {
break;
}
......@@ -550,6 +559,18 @@ int main(int argc, char *argv[])
}
free(indm_str);
if ( geometry == NULL ) {
ERROR("You need to specify a geometry file with --geometry\n");
return 1;
}
det = get_detector_geometry(geometry);
if ( det == NULL ) {
ERROR("Failed to read detector geometry from '%s'\n", geometry);
return 1;
}
free(geometry);
cell = load_cell_from_pdb(pdb);
if ( cell == NULL ) {
if ( pdb == NULL ) {
......@@ -603,6 +624,7 @@ int main(int argc, char *argv[])
pargs->config_unpolar = config_unpolar;
pargs->config_sanity = config_sanity;
pargs->cell = cell;
pargs->det = det;
pargs->indm = indm;
pargs->intensities = intensities;
pargs->counts = counts;
......@@ -695,6 +717,8 @@ int main(int argc, char *argv[])
}
free(prefix);
free(det->panels);
free(det);
free(cell);
fclose(fh);
......
......@@ -47,6 +47,7 @@ static void show_help(const char *s)
" intensities file)\n"
" --simulation-details Show technical details of the simulation.\n"
" --gpu Use the GPU to speed up the calculation.\n"
" -g. --geometry=<file> Get detector geometry from file.\n"
"\n"
" --near-bragg Output h,k,l,I near Bragg conditions.\n"
" -n, --number=<N> Generate N images. Default 1.\n"
......@@ -187,6 +188,7 @@ int main(int argc, char *argv[])
char *filename = NULL;
char *grad_str = NULL;
char *outfile = NULL;
char *geometry = NULL;
GradientMethod grad;
int ndone = 0; /* Number of simulations done (images or not) */
int number = 1; /* Number used for filename of image */
......@@ -208,14 +210,16 @@ int main(int argc, char *argv[])
{"no-noise", 0, &config_nonoise, 1},
{"intensities", 1, NULL, 'i'},
{"powder", 1, NULL, 'w'},
{"gradients", 1, NULL, 'g'},
{"gradients", 1, NULL, 't'},
{"pdb", 1, NULL, 'p'},
{"output", 1, NULL, 'o'},
{"geometry", 1, NULL, 'g'},
{0, 0, NULL, 0}
};
/* Short options */
while ((c = getopt_long(argc, argv, "hrn:i:g:p:o:", longopts, NULL)) != -1) {
while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:",
longopts, NULL)) != -1) {
switch (c) {
case 'h' : {
......@@ -238,7 +242,7 @@ int main(int argc, char *argv[])
break;
}
case 'g' : {
case 't' : {
grad_str = strdup(optarg);
break;
}
......@@ -258,6 +262,11 @@ int main(int argc, char *argv[])
break;
}
case 'g' : {
geometry = strdup(optarg);
break;
}
case 0 : {
break;
}
......@@ -329,6 +338,11 @@ int main(int argc, char *argv[])
return 1;
}
if ( geometry == NULL ) {
ERROR("You need to specify a geometry file with --geometry\n");
return 1;
}
/* Define image parameters */
image.width = 1024;
image.height = 1024;
......@@ -343,7 +357,12 @@ int main(int argc, char *argv[])
image.f0 = 1.0;
image.f0_available = 1;
#include "geometry-lcls.tmp"
image.det = get_detector_geometry(geometry);
if ( image.det == NULL ) {
ERROR("Failed to read detector geometry from '%s'\n", geometry);
return 1;
}
free(geometry);
powder = calloc(image.width*image.height, sizeof(*powder));
......@@ -461,7 +480,8 @@ skip:
cleanup_gpu(gctx);
}
free(image.det.panels);
free(image.det->panels);
free(image.det);
free(powder);
free(cell);
free(intensities);
......
......@@ -166,7 +166,7 @@ static int integrate_peak(struct image *image, int xp, int yp,
if ( !((flags & 0x01) && (flags & 0x04)) ) return 1;
}
p = find_panel(&image->det, x+xp, y+yp);
p = find_panel(image->det, x+xp, y+yp);
if ( p == NULL ) return 1;
/* Area of one pixel */
......
......@@ -35,6 +35,7 @@ static void show_help(const char *s)
"Compare intensity lists.\n"
"\n"
" -h, --help Display this help message.\n"
" -g. --geometry=<file> Get detector geometry from file.\n"
"\n");
}
......@@ -46,16 +47,18 @@ int main(int argc, char *argv[])
int x, y;
struct hdfile *hdfile;
char *filename = NULL;
char *geometry = NULL;
/* Long options */
const struct option longopts[] = {
{"help", 0, NULL, 'h'},
{"input", 1, NULL, 'i'},
{"geometry", 1, NULL, 'g'},
{0, 0, NULL, 0}
};
/* Short options */
while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) {
while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) {
switch (c) {
case 'h' : {
......@@ -72,6 +75,11 @@ int main(int argc, char *argv[])
break;
}
case 'g' : {
geometry = strdup(optarg);
break;
}
default : {
return 1;
}
......@@ -84,7 +92,17 @@ int main(int argc, char *argv[])
return 1;
}
#include "geometry-lcls.tmp"
if ( geometry == NULL ) {
ERROR("You need to specify a geometry file with --geometry\n");
return 1;
}
image.det = get_detector_geometry(geometry);
if ( image.det == NULL ) {
ERROR("Failed to read detector geometry from '%s'\n", geometry);
return 1;
}
free(geometry);
hdfile = hdfile_open(filename);
hdfile_set_image(hdfile, "/data/data");
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment