aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rwxr-xr-xscripts/test-facetron5
-rw-r--r--src/Makefile.am2
-rw-r--r--src/Makefile.in5
-rw-r--r--src/facetron.c247
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;
}