diff options
| -rwxr-xr-x | scripts/test-facetron | 5 | ||||
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/Makefile.in | 5 | ||||
| -rw-r--r-- | src/facetron.c | 247 |
4 files changed, 223 insertions, 36 deletions
diff --git a/scripts/test-facetron b/scripts/test-facetron index e10132d7..823fa72d 100755 --- a/scripts/test-facetron +++ b/scripts/test-facetron @@ -1,8 +1,7 @@ #!/bin/sh -echo $1 | indexamajig --indexing=dirax -i - \ +echo $1 | ~/crystfel/src/indexamajig --indexing=dirax -i - \ -g ~/crystfel/examples/lcls-dec.geom \ -p ~/lcls/1JB0.pdb \ --near-bragg --filter-cm \ - 2>/dev/null | grep 'star' | cut -d"=" -f 2 | cut -d" " -f 1-4 \ - | ~/crystfel/src/facetron + | ~/crystfel/src/facetron -g ~/crystfel/examples/lcls-dec.geom diff --git a/src/Makefile.am b/src/Makefile.am index bbf3bd21..be423094 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -55,7 +55,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \ sfac.c calibrate_detector_LDADD = @LIBS@ -facetron_SOURCES = facetron.c cell.c +facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c facetron_LDADD = @LIBS@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c diff --git a/src/Makefile.in b/src/Makefile.in index 6a911fa7..4fae1697 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -64,11 +64,12 @@ am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \ statistics.$(OBJEXT) symmetry.$(OBJEXT) compare_hkl_OBJECTS = $(am_compare_hkl_OBJECTS) compare_hkl_DEPENDENCIES = -am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) am_cubeit_OBJECTS = cubeit.$(OBJEXT) cell.$(OBJEXT) \ hdf5-file.$(OBJEXT) utils.$(OBJEXT) cubeit_OBJECTS = $(am_cubeit_OBJECTS) cubeit_DEPENDENCIES = +am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \ + hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) facetron_OBJECTS = $(am_facetron_OBJECTS) facetron_DEPENDENCIES = am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \ @@ -281,7 +282,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \ sfac.c calibrate_detector_LDADD = @LIBS@ -facetron_SOURCES = facetron.c cell.c +facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c facetron_LDADD = @LIBS@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c cubeit_LDADD = @LIBS@ diff --git a/src/facetron.c b/src/facetron.c index bcaaac48..7a1eec53 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -20,9 +20,15 @@ #include <string.h> #include <unistd.h> #include <getopt.h> +#include <errno.h> #include "image.h" #include "cell.h" +#include "hdf5-file.h" +#include "detector.h" + + +#define MAX_HITS (1024) static void show_help(const char *s) @@ -33,70 +39,208 @@ static void show_help(const char *s) "\n" " -h, --help Display this help message.\n" "\n" -" -i, --input=<filename> Specify the name of the image to work on.\n" -" -m, --matrix=<filename> Specify the file which contains the initial\n" -" orientation matrix. Can be '-' for stdin,\n" -" which is the default. Units are nm^-1.\n" +" -i, --input=<filename> Specify the name of the input stream.\n" +" Can be '-' for stdin.\n" +" -g. --geometry=<file> Get detector geometry from file.\n" ); } -static UnitCell *read_orientation_matrix(const char *filename) +/* Predict peaks */ +static int find_intersections(struct image *image, UnitCell *cell, + double divergence, double bandwidth) { - FILE *mfh; - float u, v, w; - struct rvec as, bs, cs; - UnitCell *cell; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + struct reflhit *hits; + int np = 0; + int hmax, kmax, lmax; + double mres; + signed int h, k, l; + + hits = malloc(sizeof(struct reflhit)*MAX_HITS); + if ( hits == NULL ) return 0; + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &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); + + 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 xda, yda; + int p; + int found = 0; + + 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 */ + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + + 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; + + } + + if ( !found ) continue; + + hits[np].h = h; + hits[np].k = k; + hits[np].l = l; + np++; + printf("%i %i %i 0.0 (at %f,%f) %i\n", h, k, l, xda, yda, p); - if ( (filename == NULL) || (strcmp(filename, "-") == 0) ) { - mfh = stdin; - } else { - mfh = fopen(filename, "r"); } - if ( mfh == NULL ) { - ERROR("Failed to open matrix file '%s'\n", filename); - return NULL; } + } + + free(hits); /* FIXME! */ - if ( fscanf(mfh, "%f %f %f", &u, &v, &w) != 3 ) { + return np; +} + + +static UnitCell *read_orientation_matrix(FILE *fh) +{ + float u, v, w; + struct rvec as, bs, cs; + UnitCell *cell; + char line[1024]; + + if ( fgets(line, 1023, fh) == NULL ) return NULL; + if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) { ERROR("Couldn't read a-star\n"); return NULL; } as.u = u*1e9; as.v = v*1e9; as.w = w*1e9; - if ( fscanf(mfh, "%f %f %f", &u, &v, &w) != 3 ) { + if ( fgets(line, 1023, fh) == NULL ) return NULL; + if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) { ERROR("Couldn't read b-star\n"); return NULL; } bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9; - if ( fscanf(mfh, "%f %f %f", &u, &v, &w) != 3 ) { + if ( fgets(line, 1023, fh) == NULL ) return NULL; + if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) { ERROR("Couldn't read c-star\n"); return NULL; } cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9; cell = cell_new_from_axes(as, bs, cs); - fclose(mfh); return cell; } +static int find_chunk(FILE *fh, UnitCell **cell, char **filename) +{ + char line[1024]; + char *rval = NULL; + + do { + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + + chomp(line); + + if ( strncmp(line, "Reflections from indexing", 25) != 0 ) { + continue; + } + + *filename = strdup(line+29); + + /* Skip two lines (while checking for errors) */ + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + + *cell = read_orientation_matrix(fh); + if ( *cell == NULL ) { + STATUS("Got filename but no cell for %s\n", *filename); + continue; + } + + return 0; + + } while ( rval != NULL ); + + return 1; +} + + int main(int argc, char *argv[]) { int c; char *infile = NULL; - char *matrix = NULL; + FILE *fh; UnitCell *cell; + char *filename; + struct detector *det; + char *geometry = NULL; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, - {"matrix", 1, NULL, 'm'}, + {"geometry", 1, NULL, 'g'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:m:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) { switch (c) { case 'h' : @@ -107,8 +251,8 @@ int main(int argc, char *argv[]) infile = strdup(optarg); break; - case 'm' : - matrix = strdup(optarg); + case 'g' : + geometry = strdup(optarg); break; case 0 : @@ -120,14 +264,57 @@ int main(int argc, char *argv[]) } - cell = read_orientation_matrix(matrix); - free(matrix); - if ( cell == NULL ) { - ERROR("Couldn't read initial orientation matrix.\n"); + if ( infile == NULL ) infile = strdup("-"); + if ( strcmp(infile, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(infile, "r"); + } + free(infile); + if ( fh == NULL ) { + ERROR("Couldn't open input stream '%s'\n", infile); + return ENOENT; + } + + det = get_detector_geometry(geometry); + if ( det == NULL ) { + ERROR("Failed to read detector geometry from '%s'\n", geometry); return 1; } + free(geometry); + + /* Loop over all successfully indexed patterns */ + while ( find_chunk(fh, &cell, &filename) == 0 ) { + + struct image image; + struct hdfile *hdfile; + + STATUS("Integrating intensities from '%s'\n", filename); + + image.det = det; + + hdfile = hdfile_open(filename); + if ( hdfile == NULL ) { + return ENOENT; + } else if ( hdfile_set_image(hdfile, "/data/data0") ) { + ERROR("Couldn't select path\n"); + return ENOENT; + } + + hdf5_read(hdfile, &image, 0); + find_intersections(&image, cell, 5.0e-3, 3.0/100.0); + + hdfile_close(hdfile); + free(cell); + free(filename); + free(image.data); + free(image.flags); + + } - cell_print(cell); + free(det->panels); + free(det); + fclose(fh); return 0; } |