From 277197d8f482229e29b05db1fa4adc866c72e1ee Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2013 15:23:05 +0100 Subject: Update XDS for new indexing subsystem --- libcrystfel/src/index.c | 13 ++- libcrystfel/src/index.h | 3 + libcrystfel/src/xds.c | 285 ++++++++++++++++++++++++++++++++---------------- libcrystfel/src/xds.h | 13 ++- 4 files changed, 211 insertions(+), 103 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 20e7cf24..0ab467a6 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -88,7 +88,8 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, break; case INDEXING_XDS : - iprivs[n] = indexing_private(indm[n]); + iprivs[n] = xds_prepare(&indm[n], cell, filename, + det, beam, ltl); break; case INDEXING_REAX : @@ -161,7 +162,7 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) break; case INDEXING_XDS : - free(priv[n]); + xds_cleanup(privs[n]); break; case INDEXING_REAX : @@ -225,7 +226,7 @@ static int try_indexer(struct image *image, IndexingMethod indm, break; case INDEXING_XDS : - run_XDS(image, cell); + return run_xds(image, ipriv); break; case INDEXING_REAX : @@ -348,6 +349,10 @@ char *indexer_str(IndexingMethod indm) strcpy(str, "grainspotter"); break; + case INDEXING_XDS : + strcpy(str, "xds"); + break; + default : ERROR("Unrecognised indexing method %i\n", indm & INDEXING_METHOD_MASK); @@ -401,7 +406,7 @@ IndexingMethod *build_indexer_list(const char *str) list[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER; } else if ( strcmp(methods[i], "xds") == 0) { - list[i] = INDEXING_XDS; + list[++nmeth] = INDEXING_DEFAULTS_XDS; } else if ( strcmp(methods[i], "reax") == 0) { list[++nmeth] = INDEXING_DEFAULTS_REAX; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index d550468b..94b0d457 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -58,6 +58,9 @@ | INDEXING_USE_LATTICE_TYPE \ | INDEXING_CHECK_PEAKS) +#define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_CHECK_PEAKS) + /** * IndexingMethod: * @INDEXING_NONE: No indexing to be performed 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 @@ -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; ifs-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); } diff --git a/libcrystfel/src/xds.h b/libcrystfel/src/xds.h index edf83b0b..7dd7c6a8 100644 --- a/libcrystfel/src/xds.h +++ b/libcrystfel/src/xds.h @@ -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 @@ -27,17 +28,23 @@ * */ -#ifndef xds_H -#define xds_H +#ifndef XDS_H +#define XDS_H #ifdef HAVE_CONFIG_H #include #endif #include "cell.h" +#include "index.h" -extern void run_xds(struct image *image, UnitCell *cell); +extern int run_xds(struct image *image, IndexingPrivate *ipriv); +extern IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl); + +extern void xds_cleanup(IndexingPrivate *pp); #endif /* XDS_H */ -- cgit v1.2.3