aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/xds.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-02-18 15:23:05 +0100
committerThomas White <taw@physics.org>2013-02-18 15:23:05 +0100
commit277197d8f482229e29b05db1fa4adc866c72e1ee (patch)
tree0506cc3dd378500d9429f51ae11be35cb37fa52b /libcrystfel/src/xds.c
parent31ec28ea0258bbb2272afa4346b55b183a43410c (diff)
Update XDS for new indexing subsystem
Diffstat (limited to 'libcrystfel/src/xds.c')
-rw-r--r--libcrystfel/src/xds.c285
1 files changed, 189 insertions, 96 deletions
diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c
index 838b71b7..786e2ae0 100644
--- a/libcrystfel/src/xds.c
+++ b/libcrystfel/src/xds.c
@@ -5,6 +5,7 @@
*
* Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
+ * Copyright © 2013 Cornelius Gati
*
* Authors:
* 2010-2013 Thomas White <taw@physics.org>
@@ -60,6 +61,15 @@
#define xds_VERBOSE 0
+/* Global private data, prepared once */
+struct xds_private
+{
+ IndexingMethod indm;
+ float *ltl;
+ UnitCell *cell;
+};
+
+
struct xds_data {
/* Low-level stuff */
@@ -179,97 +189,139 @@ static int xds_readable(struct image *image, struct xds_data *xds)
}
+static int check_cell(struct xds_private *xp, struct image *image,
+ UnitCell *cell)
+{
+ UnitCell *out;
+ Crystal *cr;
+
+ if ( xp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) {
+
+ out = match_cell(cell, xp->cell, 0, xp->ltl, 1);
+ if ( out == NULL ) return 0;
+
+ } else if ( xp->indm & INDEXING_CHECK_CELL_AXES ) {
+
+ out = match_cell(cell, xp->cell, 0, xp->ltl, 0);
+ if ( out == NULL ) return 0;
+
+ } else {
+ out = cell_new_from_cell(cell);
+ }
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+
+ crystal_set_cell(cr, out);
+
+ if ( xp->indm & INDEXING_CHECK_PEAKS ) {
+ if ( !peak_sanity_check(image, &cr, 1) ) {
+ cell_free(out);
+ crystal_free(cr);
+ return 0;
+ }
+ }
+ image_add_crystal(image, cr);
+ return 1;
+}
-static int read_newmat(const char *filename, struct image *image)
+static int read_cell(const char *filename, struct image *image,
+ struct xds_private *xp)
{
FILE * fh;
- float axstar, aystar, azstar, bxstar, bystar, bzstar, cxstar, cystar, czstar;
- char asx[11], asy[11], asz[11];
+ float axstar, aystar, azstar;
+ float bxstar, bystar, bzstar;
+ float cxstar, cystar, czstar;
+ char asx[11], asy[11], asz[11];
char bsx[11], bsy[11], bsz[11];
char csx[11], csy[11], csz[11];
char *rval, line[1024];
int r;
-
+ UnitCell *cell;
+
fh = fopen("IDXREF.LP", "r");
if ( fh == NULL ) {
ERROR("Couldn't open 'IDXREF.LP'\n");
- return 1;
+ return 0;
}
do {
rval = fgets(line, 1023, fh);
- } while (strcmp(line, " # COORDINATES OF REC. BASIS VECTOR LENGTH 1/LENGTH\n") != 0 );
+ } while ( strcmp(line, " # COORDINATES OF REC. BASIS VECTOR"
+ " LENGTH 1/LENGTH\n") != 0 );
- /*free line after chunk*/
+ /* Free line after chunk */
+ rval = fgets(line, 1023, fh);
rval = fgets(line, 1023, fh);
+
+ memcpy(asx, line+7, 10); asx[10] = '\0';
+ memcpy(asy, line+17, 10); asy[10] = '\0';
+ memcpy(asz, line+27, 10); asz[10] = '\0';
+
rval = fgets(line, 1023, fh);
-
- memcpy(asx, line+7, 10); asx[10] = '\0';
- memcpy(asy, line+17, 10); asy[10] = '\0';
- memcpy(asz, line+27, 10); asz[10] = '\0';
-
- rval = fgets(line, 1023, fh);
-
- memcpy(bsx, line+7, 10); bsx[10] = '\0';
- memcpy(bsy, line+17, 10); bsy[10] = '\0';
- memcpy(bsz, line+27, 10); bsz[10] = '\0';
-
- rval = fgets(line, 1023, fh);
-
- memcpy(csx, line+7, 10); csx[10] = '\0';
- memcpy(csy, line+17, 10); csy[10] = '\0';
- memcpy(csz, line+27, 10); csz[10] = '\0';
-
- r = sscanf(asx, "%f", &cxstar);
- r += sscanf(asy, "%f", &cystar);
- r += sscanf(asz, "%f", &czstar);
- r += sscanf(bsx, "%f", &bxstar);
- r += sscanf(bsy, "%f", &bystar);
- r += sscanf(bsz, "%f", &bzstar);
- r += sscanf(csx, "%f", &axstar);
- r += sscanf(csy, "%f", &aystar);
- r += sscanf(csz, "%f", &azstar);
-
- if ( r != 9 ) {
- STATUS("Fewer than 9 parameters found in NEWMAT file.\n");
- return 1;}
-
- fclose(fh);
-
- //printf("asx='%s'\n", asx);
- //printf("asy='%s'\n", asy);
- //printf("asz='%s'\n", asz);
- //printf("cxstar='%f'\n", cxstar);
- //printf("cystar='%f'\n", cystar);
- //printf("czstar='%f'\n", czstar);
- //printf("bsx='%s'\n", bsx);
- //printf("bsy='%s'\n", bsy);
- //printf("bsz='%s'\n", bsz);
- //printf("axstar='%f'\n", axstar);
- //printf("aystar='%f'\n", aystar);
- //printf("azstar='%f'\n", azstar);
- //printf("csx='%s'\n", csx);
- //printf("csy='%s'\n", csy);
- //printf("csz='%s'\n", csz);
- //printf("bxstar='%f'\n", bxstar);
- //printf("bystar='%f'\n", bystar);
- //printf("bzstar='%f'\n", bzstar);
-
- image->candidate_cells[0] = cell_new();
-
- cell_set_reciprocal(image->candidate_cells[0],
+
+ memcpy(bsx, line+7, 10); bsx[10] = '\0';
+ memcpy(bsy, line+17, 10); bsy[10] = '\0';
+ memcpy(bsz, line+27, 10); bsz[10] = '\0';
+
+ rval = fgets(line, 1023, fh);
+
+ memcpy(csx, line+7, 10); csx[10] = '\0';
+ memcpy(csy, line+17, 10); csy[10] = '\0';
+ memcpy(csz, line+27, 10); csz[10] = '\0';
+
+ r = sscanf(asx, "%f", &cxstar);
+ r += sscanf(asy, "%f", &cystar);
+ r += sscanf(asz, "%f", &czstar);
+ r += sscanf(bsx, "%f", &bxstar);
+ r += sscanf(bsy, "%f", &bystar);
+ r += sscanf(bsz, "%f", &bzstar);
+ r += sscanf(csx, "%f", &axstar);
+ r += sscanf(csy, "%f", &aystar);
+ r += sscanf(csz, "%f", &azstar);
+
+ if ( r != 9 ) {
+ STATUS("Fewer than 9 parameters found in NEWMAT file.\n");
+ return 0;
+ }
+
+ fclose(fh);
+
+ //printf("asx='%s'\n", asx);
+ //printf("asy='%s'\n", asy);
+ //printf("asz='%s'\n", asz);
+ //printf("cxstar='%f'\n", cxstar);
+ //printf("cystar='%f'\n", cystar);
+ //printf("czstar='%f'\n", czstar);
+ //printf("bsx='%s'\n", bsx);
+ //printf("bsy='%s'\n", bsy);
+ //printf("bsz='%s'\n", bsz);
+ //printf("axstar='%f'\n", axstar);
+ //printf("aystar='%f'\n", aystar);
+ //printf("azstar='%f'\n", azstar);
+ //printf("csx='%s'\n", csx);
+ //printf("csy='%s'\n", csy);
+ //printf("csz='%s'\n", csz);
+ //printf("bxstar='%f'\n", bxstar);
+ //printf("bystar='%f'\n", bystar);
+ //printf("bzstar='%f'\n", bzstar);
+
+ cell = cell_new();
+ cell_set_reciprocal(cell,
axstar*10e9, aystar*10e9, azstar*10e9,
bxstar*10e9, bystar*10e9, bzstar*10e9,
-cxstar*10e9, -cystar*10e9, -czstar*10e9);
- image->ncells = 1;
+ r = check_cell(xp, image, cell);
+ cell_free(cell);
- cell_print(image->candidate_cells[0]);
-
- return 0;
+ return r;
}
@@ -279,7 +331,7 @@ static void write_spot(struct image *image, const char *filename)
int i;
double fclen = 99.0e-3; /* fake camera length in m */
int n;
-
+
fh = fopen("SPOT.XDS", "w");
if ( !fh ) {
ERROR("Couldn't open temporary file '%s'\n", "SPOT.XDS");
@@ -287,8 +339,8 @@ static void write_spot(struct image *image, const char *filename)
}
n = image_feature_count(image->features);
- for ( i=0; i<n; i++ )
-
+ for ( i=0; i<n; i++ )
+
{
struct imagefeature *f;
struct panel *p;
@@ -304,9 +356,9 @@ static void write_spot(struct image *image, const char *filename)
ys = (f->fs-p->min_fs)*p->fsy + (f->ss-p->min_ss)*p->ssy;
rx = ((xs + p->cnx) / p->res);
ry = ((ys + p->cny) / p->res);
-
- //printf("xs=%f ys=%f ----> rx=%f ry=%f\n", xs, ys, rx, ry);
-
+
+ //printf("xs=%f ys=%f ----> rx=%f ry=%f\n", xs, ys, rx, ry);
+
x = rx*fclen/p->clen;
y = ry*fclen/p->clen; /* Peak positions in m */
@@ -316,7 +368,7 @@ static void write_spot(struct image *image, const char *filename)
y = (ry / 70e-6) + 1500;
if (f->intensity <= 0) continue;
-
+
fprintf(fh, "%10.2f %10.2f %10.2f %10.0f.\n",
x, y, 0.5, f->intensity);
@@ -329,7 +381,7 @@ static char *write_inp(struct image *image)
{
FILE *fh;
char *filename;
-
+
filename = malloc(1024);
if ( filename == NULL ) return NULL;
@@ -351,7 +403,7 @@ static char *write_inp(struct image *image)
fprintf(fh, "NAME_TEMPLATE_OF_DATA_FRAMES=/home/dikay/insu/data_exp1/ins_ssad_1_???.img \n");
fprintf(fh, "DATA_RANGE=1 1\n");
fprintf(fh, "SPOT_RANGE=1 1\n");
- fprintf(fh, "SPACE_GROUP_NUMBER= 94\n"); //CatB 94
+ fprintf(fh, "SPACE_GROUP_NUMBER= 94\n"); //CatB 94
fprintf(fh, "UNIT_CELL_CONSTANTS= 126.2 126.2 54.2 90 90 90\n"); //CatB 125.4 125.4 54.56 90 90 90
//fprintf(fh, "SPACE_GROUP_NUMBER=194\n"); //PS1 194
//fprintf(fh, "UNIT_CELL_CONSTANTS=281 281 165.2 90 90 120\n"); //PS1 281 281 165.2 90 90 120
@@ -372,7 +424,7 @@ static char *write_inp(struct image *image)
fprintf(fh, "INDEX_ERROR= 0.4\n");
//fprintf(fh, "INDEX_QUALITY= 0.5\n");
//fprintf(fh, "REFINE(IDXREF)= ALL\n");
- //fprintf(fh, "MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT= 1\n");
+ //fprintf(fh, "MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT= 1\n");
//fprintf(fh, "MAXIMUM_ERROR_OF_SPOT_POSITION= 20.0\n");
fclose(fh);
@@ -380,35 +432,36 @@ static char *write_inp(struct image *image)
return filename;
}
-void run_xds(struct image *image, UnitCell *cell)
+
+int run_xds(struct image *image, IndexingPrivate *priv)
{
-
- unsigned int opts;
- int status;
- int rval;
- int n;
- struct xds_data *xds;
- char *inp_filename;
-
+ unsigned int opts;
+ int status;
+ int rval;
+ int n;
+ struct xds_data *xds;
+ char *inp_filename;
+ struct xds_private *xp = (struct xds_private *)priv;
+
xds = malloc(sizeof(struct xds_data));
if ( xds == NULL ) {
ERROR("Couldn't allocate memory for xds data.\n");
- return;
+ return 0;
}
- xds->target_cell = cell;
+ xds->target_cell = xp->cell;
write_inp(image);
inp_filename = write_inp(image);
if ( inp_filename == NULL ) {
ERROR("Failed to write XDS.INP file for XDS.\n");
- return;
+ return 0;
}
-
+
n= image_feature_count(image->features);
- if (n < 25) return;
-
+ if (n < 25) return 0;
+
snprintf(xds->spotfile, 127, "SPOT.XDS");
write_spot(image, xds->spotfile);
@@ -418,7 +471,7 @@ void run_xds(struct image *image, UnitCell *cell)
if ( xds->pid == -1 ) {
ERROR("Failed to fork for XDS\n");
free(xds);
- return;
+ return 0;
}
if ( xds->pid == 0 ) {
@@ -429,7 +482,7 @@ void run_xds(struct image *image, UnitCell *cell)
tcgetattr(STDIN_FILENO, &t);
t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
tcsetattr(STDIN_FILENO, TCSANOW, &t);
-
+
execlp("xds", "", (char *)NULL);
ERROR("Failed to invoke XDS.\n");
_exit(0);
@@ -474,6 +527,7 @@ void run_xds(struct image *image, UnitCell *cell)
default:
ERROR("select() failed: %s\n", strerror(err));
rval = 1;
+ break;
}
@@ -489,15 +543,54 @@ void run_xds(struct image *image, UnitCell *cell)
close(xds->pty);
free(xds->rbuffer);
waitpid(xds->pid, &status, 0);
-
+
//FIXME//
//if ( xds->finished_ok == 1 ) {
// ERROR("XDS doesn't seem to be working properly.\n");
//} else {
/* Read the XDS NEWMAT file and get cell if found */
- read_newmat(xds->newmatfile, image);
-
+ rval = read_cell(xds->newmatfile, image, xp);
+
//}
free(xds);
+
+ return rval;
+}
+
+
+IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
+ const char *filename,
+ struct detector *det,
+ struct beam_params *beam, float *ltl)
+{
+ struct xds_private *xp;
+
+ if ( cell == NULL ) {
+ ERROR("XDS needs a unit cell.\n");
+ return NULL;
+ }
+
+ xp = calloc(1, sizeof(*xp));
+ if ( xp == NULL ) return NULL;
+
+ /* Flags that XDS knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS
+ | INDEXING_CHECK_CELL_AXES | INDEXING_USE_LATTICE_TYPE
+ | INDEXING_CHECK_PEAKS;
+
+ xp->ltl = ltl;
+ xp->cell = cell;
+ xp->indm = *indm;
+
+ return (IndexingPrivate *)xp;
+}
+
+
+void xds_cleanup(IndexingPrivate *pp)
+{
+ struct xds_private *xp;
+
+ xp = (struct xds_private *)pp;
+ free(xp);
}