Commit 533c1043 authored by Thomas White's avatar Thomas White
Browse files

facetron: Work on geometry

parent f47803a3
......@@ -59,7 +59,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
calibrate_detector_LDADD = @LIBS@
facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
image.c geometry.c reflections.c stream.c thread-pool.c
image.c geometry.c reflections.c stream.c thread-pool.c \
beam-parameters.c
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
......
......@@ -312,7 +312,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
calibrate_detector_LDADD = @LIBS@
facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
image.c geometry.c reflections.c stream.c thread-pool.c
image.c geometry.c reflections.c stream.c thread-pool.c \
beam-parameters.c
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
......
......@@ -31,6 +31,7 @@
#include "geometry.h"
#include "peaks.h"
#include "thread-pool.h"
#include "beam-parameters.h"
static void show_help(const char *s)
......@@ -45,6 +46,9 @@ static void show_help(const char *s)
" (must be a file, not e.g. stdin)\n"
" -o, --output=<filename> Output filename. Default: facetron.hkl.\n"
" -g. --geometry=<file> Get detector geometry from file.\n"
" -b, --beam=<file> Get beam parameters from file (provides initial\n"
" values for parameters, and nominal wavelengths\n"
" if no per-shot value is found in an HDF5 file.\n"
" -x, --prefix=<p> Prefix filenames from input file with <p>.\n"
" --basename Remove the directory parts of the filenames.\n"
" --no-check-prefix Don't attempt to correct the --prefix.\n"
......@@ -91,6 +95,7 @@ static void integrate_image(int mytask, void *tasks)
int j, n;
struct hdfile *hdfile;
struct image *image = pargs->image;
double nominal_photon_energy = pargs->image->beam->photon_energy;
hdfile = hdfile_open(image->filename);
if ( hdfile == NULL ) {
......@@ -102,34 +107,31 @@ static void integrate_image(int mytask, void *tasks)
return;
}
if ( hdf5_read(hdfile, pargs->image, 0) ) {
/* FIXME: Nominal photon energy */
if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
ERROR("Couldn't read '%s'\n", image->filename);
hdfile_close(hdfile);
return;
}
/* Figure out which spots should appear in this pattern */
spots = find_intersections(image, image->indexed_cell,
image->div, image->bw, &n, 0);
spots = find_intersections(image, image->indexed_cell, &n, 1);
/* For each reflection, estimate the partiality */
for ( j=0; j<n; j++ ) {
signed int h, k, l;
float i_partial;
double p;
float xc, yc;
float i_full_est;
h = spots[j].h;
k = spots[j].k;
l = spots[j].l;
/* Calculated partiality of this spot in this pattern */
p = partiality(image, h, k, l);
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( p < 0.1 ) continue;
if ( spots[j].p < 0.1 ) continue;
/* Actual measurement of this reflection from this
* pattern? */
......@@ -139,8 +141,10 @@ static void integrate_image(int mytask, void *tasks)
continue;
}
i_full_est = i_partial * spots[j].p;
pthread_mutex_lock(pargs->list_lock);
integrate_intensity(pargs->i_full, h, k, l, i_partial);
integrate_intensity(pargs->i_full, h, k, l, i_full_est);
integrate_count(pargs->cts, h, k, l, 1);
if ( !find_item(pargs->obs, h, k, l) ) {
add_item(pargs->obs, h, k, l);
......@@ -249,6 +253,7 @@ int main(int argc, char *argv[])
int n_total_patterns;
struct image *images;
int n_iter = 10;
struct beam_params *beam = NULL;
/* Long options */
const struct option longopts[] = {
......@@ -256,6 +261,7 @@ int main(int argc, char *argv[])
{"input", 1, NULL, 'i'},
{"output", 1, NULL, 'o'},
{"geometry", 1, NULL, 'g'},
{"beam", 1, NULL, 'b'},
{"prefix", 1, NULL, 'x'},
{"basename", 0, &config_basename, 1},
{"no-check-prefix", 0, &config_checkprefix, 0},
......@@ -265,7 +271,7 @@ int main(int argc, char *argv[])
};
/* Short options */
while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:",
while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
longopts, NULL)) != -1)
{
......@@ -302,6 +308,15 @@ int main(int argc, char *argv[])
n_iter = atoi(optarg);
break;
case 'b' :
beam = get_beam_parameters(optarg);
if ( beam == NULL ) {
ERROR("Failed to load beam parameters"
" from '%s'\n", optarg);
return 1;
}
break;
case 0 :
break;
......@@ -348,6 +363,11 @@ int main(int argc, char *argv[])
}
free(geomfile);
if ( beam == NULL ) {
ERROR("You must provide a beam parameters file.\n");
return 1;
}
/* Prepare for iteration */
i_full = new_list_intensity();
obs = new_items();
......@@ -388,13 +408,14 @@ int main(int argc, char *argv[])
snprintf(fnamereal, 1023, "%s%s", prefix, filename);
images[i].filename = fnamereal;
images[i].div = 0.5e-3;
images[i].bw = 0.001;
images[i].div = beam->divergence;
images[i].bw = beam->bandwidth;
images[i].orientation.w = 1.0;
images[i].orientation.x = 0.0;
images[i].orientation.y = 0.0;
images[i].orientation.z = 0.0;
images[i].det = det;
images[i].beam = beam;
/* Muppet proofing */
images[i].data = NULL;
......
......@@ -15,18 +15,93 @@
#include <stdlib.h>
#include <gsl/gsl_poly.h>
#include <assert.h>
#include "utils.h"
#include "cell.h"
#include "image.h"
#include "peaks.h"
#include "beam-parameters.h"
#define MAX_CPEAKS (1024)
#define MAX_CPEAKS (256 * 256)
static signed int locate_peak(double x, double y, double z, double lambda,
struct detector *det, double *xdap, double *ydap)
{
int p;
signed int found = -1;
const double den = 1.0/lambda + z;
for ( p=0; p<det->n_panels; p++ ) {
double xd, yd, cl;
double xda, yda;
/* Camera length for this panel */
cl = det->panels[p].clen;
/* Coordinates of peak relative to central beam, in m */
xd = cl * x / den;
yd = cl * y / den;
/* Convert to pixels */
xd *= det->panels[p].res;
yd *= det->panels[p].res;
/* Add the coordinates of the central beam */
xda = xd + det->panels[p].cx;
yda = yd + det->panels[p].cy;
/* Now, is this on this panel? */
if ( xda < det->panels[p].min_x ) continue;
if ( xda > det->panels[p].max_x ) continue;
if ( yda < det->panels[p].min_y ) continue;
if ( yda > det->panels[p].max_y ) continue;
/* If peak appears on multiple panels, reject it */
if ( found != -1 ) return -1;
/* Woohoo! */
found = p;
*xdap = xda;
*ydap = yda;
}
return found;
}
static double excitation_error(double xl, double yl, double zl,
double ds_sq, double lambda)
{
double tt;
double a, b, c;
double r1, r2;
int n;
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/lambda);
a = 1.0;
b = - cos(tt) * 2.0/lambda;
c = pow(lambda, -2.0) - ds_sq;
/* FIXME: I don't trust GSL's quadratic solver */
n = gsl_poly_solve_quadratic(a, b, c, &r1, &r2);
assert(n == 2);
r1 -= 1.0/lambda;
r2 -= 1.0/lambda;
if ( r1 > r2 ) return r2;
return r1;
}
struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double divergence, double bandwidth,
int *n, int output)
{
double asx, asy, asz;
......@@ -37,6 +112,11 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
int hmax, kmax, lmax;
double mres;
signed int h, k, l;
double bandwidth = image->bw;
double divergence = image->div;
double lambda = image->lambda;
const double profile_cutoff = 0.05e9; /* 0.1 nm^-1 */
double llow, lhigh; /* Wavelength */
cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
if ( cpeaks == NULL ) {
......@@ -48,103 +128,91 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
/* FIXME: Account for left-handed indexing */
asz = -asz; bsz = -bsz; csz = -csz;
mres = 1.0 / 8.0e-10; /* 8 Angstroms */
hmax = mres / modulus(asx, asy, asz);
kmax = mres / modulus(bsx, bsy, bsz);
lmax = mres / modulus(csx, csy, csz);
/* "low" gives the largest Ewald sphere,
* "high" gives the smallest Ewald sphere. */
llow = lambda - lambda*bandwidth/2.0;
lhigh = lambda + lambda*bandwidth/2.0;
for ( h=-hmax; h<hmax; h++ ) {
for ( k=-kmax; k<kmax; k++ ) {
for ( l=-lmax; l<lmax; l++ ) {
double xl, yl, zl;
double ds_sq, dps_sq;
double delta, divfact;
double llow, lhigh;
double xd, yd, cl;
double ds, ds_sq;
double rlow, rhigh; /* "Excitation error" */
signed int plow, phigh;
double xdalow, xdahigh;
double ydalow, ydahigh;
double xda, yda;
int p;
int found = 0;
/* Ignore central beam */
if ( (h==0) && (k==0) && (l==0) ) continue;
llow = image->lambda - image->lambda*bandwidth/2.0;
lhigh = image->lambda + image->lambda*bandwidth/2.0;
/* Get the coordinates of the reciprocal lattice point */
zl = h*asz + k*bsz + l*csz;
if ( zl < 0.0 ) continue; /* Do this check very early */
if ( zl > 0.0 ) continue; /* Throw out if it's "in front" */
xl = h*asx + k*bsx + l*csx;
yl = h*asy + k*bsy + l*csy;
/* Calculate reciprocal lattice point modulus (and square) */
ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
delta = divergence/image->lambda;
dps_sq = ds_sq + pow(delta, 2.0); /* d'*^2 */
/* In range? */
divfact = 2.0 * delta * sqrt(xl*xl + yl*yl);
if ( ds_sq - 2.0*zl/llow > 0.0 ) continue;
if ( ds_sq - 2.0*zl/lhigh < 0.0 ) continue;
/* Work out which panel this peak would fall on */
for ( p=0; p<image->det->n_panels; p++ ) {
/* Camera length for this panel */
cl = image->det->panels[p].clen;
/* Coordinates of peak relative to central beam, in m */
xd = cl*xl / (ds_sq/(2.0*zl) - zl);
yd = cl*yl / (ds_sq/(2.0*zl) - zl);
/* Convert to pixels */
xd *= image->det->panels[p].res;
yd *= image->det->panels[p].res;
/* Add the coordinates of the central beam */
xda = xd + image->det->panels[p].cx;
yda = yd + image->det->panels[p].cy;
/* Now, is this on this panel? */
if ( xda < image->det->panels[p].min_x ) continue;
if ( xda > image->det->panels[p].max_x ) continue;
if ( yda < image->det->panels[p].min_y ) continue;
if ( yda > image->det->panels[p].max_y ) continue;
/* Woohoo! */
found = 1;
break;
ds = sqrt(ds_sq);
if ( ds > mres ) continue; /* Outside resolution range */
/* Calculate excitation errors */
rlow = excitation_error(xl, yl, zl, ds_sq, llow);
rhigh = excitation_error(xl, yl, zl, ds_sq, lhigh);
if ( (fabs(rlow) > profile_cutoff)
&& (fabs(rhigh) > profile_cutoff) ) {
/* Reflection is not close to Bragg at either of
* the two extremes of wavelength, so skip it. */
continue;
}
if ( !found ) continue;
/* Locate peak on detector, and check it doesn't span panels */
plow = locate_peak(xl, yl, zl, llow, image->det,
&xdalow, &ydalow);
if ( plow == -1 ) continue;
phigh = locate_peak(xl, yl, zl, lhigh, image->det,
&xdahigh, &ydahigh);
if ( phigh == -1 ) continue;
if ( plow != phigh ) continue;
xda = xdahigh;
yda = ydahigh;
cpeaks[np].h = h;
cpeaks[np].k = k;
cpeaks[np].l = l;
cpeaks[np].x = xda;
cpeaks[np].y = yda;
cpeaks[np].x = xdahigh;
cpeaks[np].y = ydahigh;
np++;
if ( output ) {
printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda);
printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %9e %9e\n",
h, k, l, 0.0, xda, yda, rlow, rhigh);
}
if ( np == MAX_CPEAKS ) goto out;
}
}
}
out:
*n = np;
return cpeaks;
}
/* Return the partiality of this reflection in this image */
double partiality(struct image *image, signed int h, signed int k, signed int l)
{
return 1.0;
}
double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
{
double itot = 0.0;
......
......@@ -18,11 +18,8 @@
#endif
extern struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double divergence, double bandwidth,
int *n, int output);
extern double partiality(struct image *image, signed int h,
signed int k, signed int l);
extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n);
......
......@@ -53,11 +53,21 @@ typedef struct _imagefeaturelist ImageFeatureList;
/* This structure represents a predicted peak in an image */
struct cpeak {
struct cpeak
{
/* Indices */
signed int h;
signed int k;
signed int l;
double min_distance;
/* Partiality */
double r1;
double r2;
double p;
/* Location in image */
int x;
int y;
};
......@@ -73,7 +83,7 @@ struct image {
UnitCell *candidate_cells[MAX_CELL_CANDIDATES];
int ncells;
struct detector *det;
struct beam_params *beam;
struct beam_params *beam; /* The nominal beam parameters */
char *filename;
struct cpeak *cpeaks;
int n_cpeaks;
......@@ -86,7 +96,7 @@ struct image {
double m; /* Mosaicity in radians */
/* Information about the radiation */
/* Per-shot radiation values */
double lambda; /* Wavelength in m */
double div; /* Divergence in radians */
double bw; /* Bandwidth as a fraction */
......
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