diff options
Diffstat (limited to 'libcrystfel')
| -rw-r--r-- | libcrystfel/Makefile.am | 7 | ||||
| -rw-r--r-- | libcrystfel/src/grainspotter.c | 500 | ||||
| -rw-r--r-- | libcrystfel/src/grainspotter.h | 55 | ||||
| -rw-r--r-- | libcrystfel/src/index.c | 44 | ||||
| -rw-r--r-- | libcrystfel/src/index.h | 13 | ||||
| -rw-r--r-- | libcrystfel/src/reax.c | 1212 | ||||
| -rw-r--r-- | libcrystfel/src/reax.h | 75 |
7 files changed, 4 insertions, 1902 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 989ce434..be7a2556 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -9,11 +9,10 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/reflist-utils.c src/filters.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ - src/grainspotter.c src/xds.c src/integration.c \ + src/xds.c src/integration.c \ src/histogram.c src/events.c src/asdf.c if HAVE_FFTW -libcrystfel_la_SOURCES += src/reax.c endif libcrystfel_la_includedir=$(includedir)/crystfel/ @@ -26,9 +25,9 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \ src/geometry.h src/peaks.h src/stream.h \ src/render.h src/index.h src/image.h \ src/filters.h src/dirax.h src/mosflm.h \ - src/reax.h src/cell-utils.h \ + src/cell-utils.h \ src/integer_matrix.h src/crystal.h \ - src/grainspotter.h src/xds.h \ + src/xds.h \ src/integration.h src/histogram.h \ src/events.h src/asdf.h diff --git a/libcrystfel/src/grainspotter.c b/libcrystfel/src/grainspotter.c deleted file mode 100644 index 261756d6..00000000 --- a/libcrystfel/src/grainspotter.c +++ /dev/null @@ -1,500 +0,0 @@ -/* - * grainspotter.c - * - * Invoke GrainSpotter for multi-crystal autoindexing - * - * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2014 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - - -#include <stdlib.h> -#include <stdio.h> -#include <math.h> -#include <pty.h> -#include <string.h> -#include <unistd.h> -#include <sys/wait.h> -#include <fcntl.h> -#include <assert.h> -#include <sys/ioctl.h> -#include <errno.h> - -#ifdef HAVE_CLOCK_GETTIME -#include <time.h> -#else -#include <sys/time.h> -#endif - -#include "image.h" -#include "utils.h" -#include "peaks.h" -#include "cell.h" -#include "cell-utils.h" -#include "grainspotter.h" - - -#define GRAINSPOTTER_VERBOSE 0 - - -/* Global private data, prepared once */ -struct grainspotter_private -{ - IndexingMethod indm; - UnitCell *cell; -}; - - -/* Data needed during the run of Grainspotter */ -struct grainspotter_data { - - struct grainspotter_private *gp; - - /* Low-level stuff */ - int pty; - pid_t pid; - char *rbuffer; - int rbufpos; - int rbuflen; - -}; - - -static int read_matrix(struct grainspotter_private *gp, struct image *image, - char *filename) -{ - FILE *fh; - int d1; - float d2; - float ubi11, ubi12, ubi13; - float ubi21, ubi22, ubi23; - float ubi31, ubi32, ubi33; - char line[1024]; - int r; - - fh = fopen(filename, "r"); - if ( fh == NULL ) { - ERROR("Can't open '%s'\n", filename); - return 1; - } - - /* Read and discard first line */ - if ( fgets(line, 1024, fh) == NULL ) { - ERROR("Failed to read GFF file.\n"); - return 1; - } - - do { - - Crystal *cr; - UnitCell *cell; - - /* One line per grain */ - if ( fgets(line, 1024, fh) == NULL ) { - ERROR("Failed to read GFF file.\n"); - return 1; - } - - STATUS("'%s'\n", line); - - r = sscanf(line, "%i %f %f %f %f %f %f %f %f %f %f %f %f" - "%f %f %f %f %f %f %f %f %f %f %f", - &d1, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, - &d2, &d2, &d2, &d2, &d2, &d2, - &ubi11, &ubi12, &ubi13, - &ubi21, &ubi22, &ubi23, - &ubi31, &ubi32, &ubi33); - - if ( r != 24 ) { - ERROR("Only %i parameters in GFF file\n", r); - return 1; - } - - cell = cell_new(); - - cell_set_cartesian(cell, ubi12/1e10, ubi13/1e10, ubi11/1e10, - ubi22/1e10, ubi23/1e10, ubi21/1e10, - ubi32/1e10, ubi33/1e10, ubi31/1e10); - - cr = crystal_new(); - if ( cr == NULL ) { - ERROR("Failed to allocate crystal.\n"); - return 0; - } - - crystal_set_cell(cr, cell); - image_add_crystal(image, cr); - - } while ( !feof(fh) ); - - fclose(fh); - - if ( gp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, image->crystals, - image->n_crystals) ) - { - free_all_crystals(image); - return 0; - } - } - - return 0; -} - - -static void gs_parseline(char *line, struct image *image, - struct grainspotter_data *gs) -{ - #if GRAINSPOTTER_VERBOSE - STATUS("%s\n", line); - #endif -} - - -static int grainspotter_readable(struct image *image, - struct grainspotter_data *gs) -{ - int rval; - int no_string = 0; - - rval = read(gs->pty, gs->rbuffer+gs->rbufpos, gs->rbuflen-gs->rbufpos); - - if ( (rval == -1) || (rval == 0) ) return 1; - - gs->rbufpos += rval; - assert(gs->rbufpos <= gs->rbuflen); - - while ( (!no_string) && (gs->rbufpos > 0) ) { - - int i; - int block_ready = 0; - - /* See if there's a full line in the buffer yet */ - for ( i=0; i<gs->rbufpos-1; i++ ) { - /* Means the last value looked at is rbufpos-2 */ - - if ( (gs->rbuffer[i] == '\r') - && (gs->rbuffer[i+1] == '\n') ) { - block_ready = 1; - break; - } - - } - - if ( block_ready ) { - - unsigned int new_rbuflen; - unsigned int endbit_length; - char *block_buffer = NULL; - - block_buffer = malloc(i+1); - memcpy(block_buffer, gs->rbuffer, i); - block_buffer[i] = '\0'; - - if ( block_buffer[0] == '\r' ) { - memmove(block_buffer, block_buffer+1, i); - } - - gs_parseline(block_buffer, image, gs); - free(block_buffer); - endbit_length = i+2; - - /* Now the block's been parsed, it should be - * forgotten about */ - memmove(gs->rbuffer, - gs->rbuffer + endbit_length, - gs->rbuflen - endbit_length); - - /* Subtract the number of bytes removed */ - gs->rbufpos = gs->rbufpos - endbit_length; - new_rbuflen = gs->rbuflen - endbit_length; - if ( new_rbuflen == 0 ) new_rbuflen = 256; - gs->rbuffer = realloc(gs->rbuffer, new_rbuflen); - gs->rbuflen = new_rbuflen; - - } else { - - if ( gs->rbufpos==gs->rbuflen ) { - - /* More buffer space is needed */ - gs->rbuffer = realloc(gs->rbuffer, - gs->rbuflen + 256); - gs->rbuflen = gs->rbuflen + 256; - /* The new space gets used at the next - * read, shortly... */ - - } - no_string = 1; - - } - - } - - return 0; -} - - -static void write_gve(struct image *image, struct grainspotter_private *gp) -{ - FILE *fh; - int i; - char filename[1024]; - double a, b, c, al, be, ga; - - snprintf(filename, 1023, "xfel-%i.gve", image->id); - - fh = fopen(filename, "w"); - if ( !fh ) { - ERROR("Couldn't open temporary file '%s'\n", filename); - return; - } - - cell_get_parameters(gp->cell, &a, &b, &c, &al, &be, &ga); - fprintf(fh, "%.6f %.6f %.6f %.6f %.6f %.6f P\n", a*1e10, b*1e10, c*1e10, - rad2deg(al), rad2deg(be), rad2deg(ga)); - fprintf(fh, "# wavelength = %.6f\n", image->lambda*1e10); - fprintf(fh, "# wedge = 0.000000\n"); - fprintf(fh, "# ds h k l\n"); - fprintf(fh, "# xr yr zr xc yc ds eta omega\n"); - - for ( i=0; i<image_feature_count(image->features); i++ ) { - - struct imagefeature *f; - - f = image_get_feature(image->features, i); - if ( f == NULL ) continue; - - fprintf(fh, "%.6f %.6f %.6f 0 0 %.6f %.6f %.6f 0\n", - f->rz/1e10, f->rx/1e10, f->ry/1e10, - modulus(f->rx, f->ry, f->rz)/1e10, /* dstar */ - rad2deg(atan2(f->ry, f->rx)), 0.0); /* eta, omega */ - - } - fclose(fh); -} - - -static char *write_ini(struct image *image) -{ - FILE *fh; - char *filename; - double tt; - - filename = malloc(1024); - if ( filename == NULL ) return NULL; - - snprintf(filename, 1023, "xfel-%i.ini", image->id); - - fh = fopen(filename, "w"); - if ( !fh ) { - ERROR("Couldn't open temporary file '%s'\n", filename); - free(filename); - return NULL; - } - - get_q_for_panel(image->det->furthest_out_panel, - image->det->furthest_out_fs, - image->det->furthest_out_ss, - &tt, 1.0/image->lambda); - - fprintf(fh, "spacegroup 96\n"); - fprintf(fh, "!dsrange 0 1.3\n"); - fprintf(fh, "tthrange 0 %.2f\n", 20.0); - fprintf(fh, "etarange 0 360\n"); - fprintf(fh, "domega 1\n"); - fprintf(fh, "omegarange -0.5 0.5\n"); - fprintf(fh, "filespecs xfel-%i.gve xfel-%i.log\n", - image->id, image->id); - fprintf(fh, "cuts 3 0.1 0.5\n"); - fprintf(fh, "eulerstep 8\n"); - fprintf(fh, "uncertainties 0.1 0.2 .5\n"); - fprintf(fh, "nsigmas 2\n"); - fprintf(fh, "minfracg 0.95\n"); - fprintf(fh, "Nhkls_in_indexing 500\n"); - fprintf(fh, "random 10000\n"); - fprintf(fh, "!positionfit\n"); - fprintf(fh, "genhkl\n"); - fprintf(fh, "xfelmode\n"); - - fclose(fh); - - return filename; -} - - -int grainspotter_index(struct image *image, IndexingPrivate *ipriv) -{ - unsigned int opts; - int status; - int rval; - struct grainspotter_data *grainspotter; - struct grainspotter_private *gp = (struct grainspotter_private *)ipriv; - char *ini_filename; - char gff_filename[1024]; - - write_gve(image, gp); - ini_filename = write_ini(image); - - if ( ini_filename == NULL ) { - ERROR("Failed to write ini file for GrainSpotter.\n"); - return 0; - } - - grainspotter = malloc(sizeof(struct grainspotter_data)); - if ( grainspotter == NULL ) { - ERROR("Couldn't allocate memory for GrainSpotter data.\n"); - return 0; - } - - grainspotter->gp = gp; - - snprintf(gff_filename, 1023, "xfel-%i.gff", image->id); - remove(gff_filename); - - grainspotter->pid = forkpty(&grainspotter->pty, NULL, NULL, NULL); - if ( grainspotter->pid == -1 ) { - ERROR("Failed to fork for GrainSpotter: %s\n", strerror(errno)); - return 0; - } - if ( grainspotter->pid == 0 ) { - - /* Child process: invoke GrainSpotter */ - struct termios t; - - /* Turn echo off */ - tcgetattr(STDIN_FILENO, &t); - t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); - tcsetattr(STDIN_FILENO, TCSANOW, &t); - - STATUS("Running GrainSpotter.0.93 '%s'\n", ini_filename); - execlp("GrainSpotter.0.93", "", ini_filename, (char *)NULL); - ERROR("Failed to invoke GrainSpotter.\n"); - _exit(0); - - } - free(ini_filename); - - grainspotter->rbuffer = malloc(256); - grainspotter->rbuflen = 256; - grainspotter->rbufpos = 0; - - /* Set non-blocking */ - opts = fcntl(grainspotter->pty, F_GETFL); - fcntl(grainspotter->pty, F_SETFL, opts | O_NONBLOCK); - - do { - - fd_set fds; - struct timeval tv; - int sval; - - FD_ZERO(&fds); - FD_SET(grainspotter->pty, &fds); - - tv.tv_sec = 30; - tv.tv_usec = 0; - - sval = select(grainspotter->pty+1, &fds, NULL, NULL, &tv); - - if ( sval == -1 ) { - - const int err = errno; - - switch ( err ) { - - case EINTR: - STATUS("Restarting select()\n"); - rval = 0; - break; - - default: - ERROR("select() failed: %s\n", strerror(err)); - rval = 1; - - } - - } else if ( sval != 0 ) { - rval = grainspotter_readable(image, grainspotter); - } else { - ERROR("No response from GrainSpotter..\n"); - rval = 1; - } - - } while ( !rval ); - - close(grainspotter->pty); - free(grainspotter->rbuffer); - waitpid(grainspotter->pid, &status, 0); - - if ( status != 0 ) { - ERROR("GrainSpotter doesn't seem to be working properly.\n"); - free(grainspotter); - return 0; - } - - if ( read_matrix(gp, image, gff_filename) != 0 ) { - free(grainspotter); - return 0; - } - - /* Success! */ - free(grainspotter); - return 1; -} - - - -IndexingPrivate *grainspotter_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) -{ - struct grainspotter_private *gp; - - if ( cell == NULL ) { - ERROR("GrainSpotter needs a unit cell.\n"); - return NULL; - } - - gp = calloc(1, sizeof(*gp)); - if ( gp == NULL ) return NULL; - - /* Flags that GrainSpotter knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS - | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS; - - gp->cell = cell; - gp->indm = *indm; - - return (IndexingPrivate *)gp; -} - - -void grainspotter_cleanup(IndexingPrivate *pp) -{ - struct grainspotter_private *p; - - p = (struct grainspotter_private *)pp; - free(p); -} diff --git a/libcrystfel/src/grainspotter.h b/libcrystfel/src/grainspotter.h deleted file mode 100644 index 3ab823ad..00000000 --- a/libcrystfel/src/grainspotter.h +++ /dev/null @@ -1,55 +0,0 @@ -/* - * grainspotter.h - * - * Invoke GrainSpotter for multi-crystal autoindexing - * - * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2013 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifndef GRAINSPOTTER_H -#define GRAINSPOTTER_H - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "cell.h" - -#ifdef __cplusplus -extern "C" { -#endif - -extern IndexingPrivate *grainspotter_prepare(IndexingMethod *indm, - UnitCell *cell, - struct detector *det, - float *ltl); - -extern void grainspotter_cleanup(IndexingPrivate *pp); - -extern int grainspotter_index(struct image *image, IndexingPrivate *p); - -#ifdef __cplusplus -} -#endif - -#endif /* GRAINSPOTTER_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 6b36af3e..c3aef88e 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -49,11 +49,8 @@ #include "xds.h" #include "detector.h" #include "index.h" -#include "reax.h" -#include "grainspotter.h" #include "geometry.h" #include "cell-utils.h" -#include "grainspotter.h" static int debug_index(struct image *image) @@ -102,16 +99,7 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, break; case INDEXING_XDS : - iprivs[n] = xds_prepare(&indm[n], cell, det, ltl); - break; - - case INDEXING_REAX : - iprivs[n] = reax_prepare(&indm[n], cell, det, ltl); - break; - - case INDEXING_GRAINSPOTTER : - iprivs[n] = grainspotter_prepare(&indm[n], cell, - det, ltl); + iprivs[n] = xds_prepare(&indm[n], cell, det, ltl); break; case INDEXING_DEBUG : @@ -184,14 +172,6 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) xds_cleanup(privs[n]); break; - case INDEXING_REAX : - reax_cleanup(privs[n]); - break; - - case INDEXING_GRAINSPOTTER : - grainspotter_cleanup(privs[n]); - break; - case INDEXING_DEBUG : free(privs[n]); break; @@ -257,14 +237,6 @@ static int try_indexer(struct image *image, IndexingMethod indm, return run_xds(image, ipriv); break; - case INDEXING_REAX : - return reax_index(ipriv, image); - break; - - case INDEXING_GRAINSPOTTER : - return grainspotter_index(image, ipriv); - break; - case INDEXING_DEBUG : return debug_index(image); @@ -389,14 +361,6 @@ char *indexer_str(IndexingMethod indm) strcpy(str, "mosflm"); break; - case INDEXING_REAX : - strcpy(str, "reax"); - break; - - case INDEXING_GRAINSPOTTER : - strcpy(str, "grainspotter"); - break; - case INDEXING_XDS : strcpy(str, "xds"); break; @@ -469,15 +433,9 @@ IndexingMethod *build_indexer_list(const char *str) } else if ( strcmp(methods[i], "mosflm") == 0) { list[++nmeth] = INDEXING_DEFAULTS_MOSFLM; - } else if ( strcmp(methods[i], "grainspotter") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER; - } else if ( strcmp(methods[i], "xds") == 0) { list[++nmeth] = INDEXING_DEFAULTS_XDS; - } else if ( strcmp(methods[i], "reax") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_REAX; - } else if ( strcmp(methods[i], "none") == 0) { list[++nmeth] = INDEXING_NONE; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index ff1cd5df..12c13df8 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -48,15 +48,6 @@ | INDEXING_CHECK_CELL_COMBINATIONS \ | INDEXING_USE_LATTICE_TYPE) -#define INDEXING_DEFAULTS_REAX (INDEXING_REAX | INDEXING_USE_LATTICE_TYPE \ - | INDEXING_USE_CELL_PARAMETERS \ - | INDEXING_CHECK_PEAKS) - -#define INDEXING_DEFAULTS_GRAINSPOTTER (INDEXING_GRAINSPOTTER \ - | INDEXING_USE_LATTICE_TYPE \ - | INDEXING_USE_CELL_PARAMETERS \ - | INDEXING_CHECK_PEAKS) - /* Axis check is needed for XDS, because it likes to permute the axes */ #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ | INDEXING_USE_CELL_PARAMETERS \ @@ -68,8 +59,6 @@ * @INDEXING_NONE: No indexing to be performed * @INDEXING_DIRAX: Invoke DirAx * @INDEXING_MOSFLM: Invoke MOSFLM - * @INDEXING_REAX: DPS algorithm using known cell parameters - * @INDEXING_GRAINSPOTTER: Invoke GrainSpotter * @INDEXING_XDS: Invoke XDS * @INDEXING_SIMULATION: Dummy value * @INDEXING_DEBUG: Results injector for debugging @@ -95,8 +84,6 @@ typedef enum { /* The core indexing methods themselves */ INDEXING_DIRAX = 1, INDEXING_MOSFLM = 2, - INDEXING_REAX = 3, - INDEXING_GRAINSPOTTER = 4, INDEXING_XDS = 5, INDEXING_SIMULATION = 6, INDEXING_DEBUG = 7, diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c deleted file mode 100644 index e2dd249a..00000000 --- a/libcrystfel/src/reax.c +++ /dev/null @@ -1,1212 +0,0 @@ -/* - * reax.c - * - * A new auto-indexer - * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * Copyright © 2014 Takanori Nakane - * - * Authors: - * 2011-2014 Thomas White <taw@physics.org> - * 2014 Takanori Nakane <nakane.t@gmail.com> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - - -#include <stdlib.h> -#include <stdio.h> -#include <math.h> -#include <assert.h> -#include <fftw3.h> -#include <fenv.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_linalg.h> -#include <gsl/gsl_eigen.h> - -#include "image.h" -#include "utils.h" -#include "peaks.h" -#include "cell.h" -#include "cell-utils.h" -#include "index.h" -#include "reax.h" - - -/* Minimum number of standard deviations above the mean a peak must be in the - * 1D FT to qualify as a candidate vector */ -#define MIN_SIGMAS (3.0) - - -/* Maximum number of times the angular tolerance that vectors are permitted to - * be together before they get merged by squash_vectors() */ -#define INC_TOL_MULTIPLIER (3.0) - - -/* Maximum number of candidate vectors to find (we will take the best ones) */ -#define MAX_CANDIDATES (1024) - - -/* Choose the best solution from this many candidate cells */ -#define MAX_REAX_CELL_CANDIDATES (32) - -struct dvec -{ - double x; - double y; - double z; - double th; - double ph; -}; - - -struct reax_candidate -{ - struct dvec v; /* This is the vector for the candidate */ - double fom; -}; - - -struct reax_search_v -{ - unsigned int smin; - unsigned int smax; /* Search for vector in this range */ - - struct reax_candidate *cand; /* Candidate vectors go here */ - int n_cand; /* There are this many candidates */ - int max_warned; -}; - - -struct reax_search -{ - struct reax_search_v *search; /* Search for these vectors */ - int n_search; /* There are this many vectors to find */ - double pmax; /* The maximum feature resolution */ -}; - - -struct reax_private -{ - IndexingMethod indm; - struct dvec *directions; - int n_dir; - double angular_inc; - UnitCell *cell; - - double *fft_in; - fftw_complex *fft_out; - fftw_plan plan; - int nel; - - fftw_complex *r_fft_in; - fftw_complex *r_fft_out; - fftw_plan r_plan; - int ch; - int cw; -}; - - -static void fill_and_transform(struct dvec *dir, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - const struct rigid_group *rg, struct detector *det) -{ - int n, i; - - for ( i=0; i<nel; i++ ) { - fft_in[i] = 0.0; - } - - n = image_feature_count(flist); - for ( i=0; i<n; i++ ) { - - struct imagefeature *f; - double val; - int idx; - - f = image_get_feature(flist, i); - if ( f == NULL ) continue; - - if ( rg != NULL ) { - - struct panel *p; - - p = find_panel(det, f->fs, f->ss); - assert(p != NULL); - - if ( !panel_is_in_rigid_group(rg, p) ) continue; - - } - - val = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z; - - idx = nel/2 + nel*val/(2.0*pmax); - fft_in[idx]++; - - } - - fftw_execute_dft_r2c(plan, fft_in, fft_out); -} - - -static void add_candidate(struct reax_search_v *s, struct reax_candidate *c) -{ - int idx; - - if ( s->n_cand == MAX_CANDIDATES ) { - if ( !s->max_warned ) { - ERROR("WARNING: Too many candidates.\n"); - s->max_warned = 1; - } - return; - } - - idx = s->n_cand++; - - s->cand[idx].v = c->v; - s->cand[idx].fom = c->fom; -} - - -static double check_dir(struct dvec *dir, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - struct reax_search *s, - const struct rigid_group *rg, struct detector *det) -{ - int i; - double tot; - - fill_and_transform(dir, flist, nel, pmax, fft_in, fft_out, - plan, rg, det); - - tot = 0.0; - for ( i=0; i<s->n_search; i++ ) { - - double tot = 0.0; - double peak = 0.0; - double peak_mod = 0.0; - double mean; - double sd = 0.0; - int j; - int n = 0; - - for ( j=0; j<nel/2+1; j++ ) { - - double re, im, am; - - re = fft_out[j][0]; - im = fft_out[j][1]; - am = sqrt(re*re + im*im); - - tot += am; - n++; - - if ( ( j >= s->search[i].smin ) - && ( j <= s->search[i].smax ) ) { - if ( am > peak ) { - peak = am; - peak_mod = (double)j/(2.0*pmax); - } - } - - } - mean = tot/(double)n; - - for ( j=0; j<nel/2+1; j++ ) { - - double re, im, am; - - re = fft_out[j][0]; - im = fft_out[j][1]; - am = sqrt(re*re + im*im); - - sd += pow(am - mean, 2.0); - - } - sd = sqrt(sd/(double)n); - - /* If sufficiently strong, add to list of candidates */ - if ( peak > mean+MIN_SIGMAS*sd ) { - - struct reax_candidate c; - - c.v.x = dir->x * peak_mod; - c.v.y = dir->y * peak_mod; - c.v.z = dir->z * peak_mod; - c.fom = peak; - - add_candidate(&s->search[i], &c); - - } - - } - - return tot; -} - - -/* Refine a direct space vector. From Clegg (1984) - * with added iteration because more reflections might get included as the - * refinement proceeds. */ -static double iterate_refine_vector(double *x, double *y, double *z, - ImageFeatureList *flist) -{ - int fi, n, err; - gsl_matrix *C; - gsl_vector *A; - gsl_vector *t; - gsl_matrix *s_vec; - gsl_vector *s_val; - double dtmax; - - A = gsl_vector_calloc(3); - C = gsl_matrix_calloc(3, 3); - - n = image_feature_count(flist); - fesetround(1); - for ( fi=0; fi<n; fi++ ) { - - struct imagefeature *f; - double val; - double kn, kno; - double xv[3]; - int i, j; - - f = image_get_feature(flist, fi); - if ( f == NULL ) continue; - - kno = f->rx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */ - kn = nearbyint(kno); - if ( fabs(kn - kno) > 0.3 ) continue; - - xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz; - - for ( i=0; i<3; i++ ) { - - val = gsl_vector_get(A, i); - gsl_vector_set(A, i, val+xv[i]*kn); - - for ( j=0; j<3; j++ ) { - val = gsl_matrix_get(C, i, j); - gsl_matrix_set(C, i, j, val+xv[i]*xv[j]); - } - - } - - } - - s_val = gsl_vector_calloc(3); - s_vec = gsl_matrix_calloc(3, 3); - err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val); - if ( err ) { - ERROR("SVD failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(C); - gsl_vector_free(A); - return 0.0; - } - - t = gsl_vector_calloc(3); - err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t); - if ( err ) { - ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(C); - gsl_vector_free(A); - gsl_vector_free(t); - return 0.0; - } - - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - - dtmax = fabs(*x - gsl_vector_get(t, 0)); - dtmax += fabs(*y - gsl_vector_get(t, 1)); - dtmax += fabs(*z - gsl_vector_get(t, 2)); - - *x = gsl_vector_get(t, 0); - *y = gsl_vector_get(t, 1); - *z = gsl_vector_get(t, 2); - - gsl_matrix_free(C); - gsl_vector_free(A); - gsl_vector_free(t); - - return dtmax; -} - - -static void refine_vector(ImageFeatureList *flist, struct dvec *dir) -{ - int i; - double sm; - - i = 0; - do { - - sm = iterate_refine_vector(&dir->x, &dir->y, &dir->z, flist); - i++; - - } while ( (sm > 0.001e-9) && (i<10) ); -} - - -static void squash_vectors(struct reax_search *s, double tol) -{ - int i; - - for ( i=0; i<s->n_search; i++ ) { - - struct reax_search_v *sv; - struct reax_candidate *new; - int j, k; - int n_invalid = 0; - int n_copied; - - sv = &s->search[i]; - - for ( j=0; j<sv->n_cand; j++ ) { - for ( k=0; k<sv->n_cand; k++ ) { - - struct reax_candidate *v1, *v2; - - if ( j == k ) continue; - - v1 = &sv->cand[j]; - v2 = &sv->cand[k]; - - if ( angle_between(v1->v.x, v1->v.y, v1->v.z, - v2->v.x, v2->v.y, v2->v.z) < tol ) - { - if ( !isnan(v1->fom) && !isnan(v2->fom ) ) { - if ( v1->fom > v2->fom ) { - v2->fom = NAN; - } else { - v1->fom = NAN; - } - n_invalid++; - } - } - - } - } - - new = calloc(sv->n_cand - n_invalid, - sizeof(struct reax_candidate)); - if ( new == NULL ) { - ERROR("Failed to allocate memory for squashed" - " candidate list.\n"); - return; - } - - n_copied = 0; - for ( j=0; j<sv->n_cand; j++ ) { - if ( !isnan(sv->cand[j].fom) ) { - - new[n_copied] = sv->cand[j]; - n_copied++; - - } - } - assert(sv->n_cand - n_invalid == n_copied); - - free(sv->cand); - //STATUS("Search vector %i:", i); - //STATUS(" squashed %i candidates down to %i\n", - // sv->n_cand, n_copied); - sv->n_cand = n_copied; - sv->cand = new; - - } -} - - -static UNUSED void show_vectors(struct reax_search *s, const char *pre) -{ - int i; - - /* For each direction being searched for */ - for ( i=0; i<s->n_search; i++ ) { - - int j; - - for ( j=0; j<s->search[i].n_cand; j++ ) { - STATUS("%s: %i/%i: %+6.2f %+6.2f %+6.2f nm %.2f\n", pre, - i, j, s->search[i].cand[j].v.x*1e9, - s->search[i].cand[j].v.y*1e9, - s->search[i].cand[j].v.z*1e9, - s->search[i].cand[j].fom); - } - - } -} - - -static void find_candidates(struct reax_private *p, - ImageFeatureList *flist, double pmax, - double *fft_in, fftw_complex *fft_out, - struct reax_search *s, - const char *rg, struct detector *det) -{ - int i; - - for ( i=0; i<s->n_search; i++ ) { - s->search[i].cand = calloc(MAX_CANDIDATES, - sizeof(struct reax_candidate)); - s->search[i].n_cand = 0; - } - - for ( i=0; i<p->n_dir; i++ ) { - check_dir(&p->directions[i], flist, - p->nel, pmax, fft_in, fft_out, p->plan, - s, NULL, NULL); - } - - squash_vectors(s, INC_TOL_MULTIPLIER*p->angular_inc); - - //show_vectors(s, "BEFORE"); - - for ( i=0; i<s->n_search; i++ ) { - - struct reax_search_v *sv; - int j; - - sv = &s->search[i]; - - for ( j=0; j<sv->n_cand; j++ ) { - refine_vector(flist, &sv->cand[j].v); - } - - } - - //show_vectors(s, "FINAL"); -} - - -/* Set up search parameters to look for all three cell axes */ -static struct reax_search *search_all_axes(UnitCell *cell, double pmax) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double mod_a, mod_b, mod_c; - double amin, amax; - double bmin, bmax; - double cmin, cmax; - unsigned int smin, smax; - const double ltol = 10.0; /* Direct space axis length tolerance in % */ - struct reax_search *s; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - mod_a = modulus(ax, ay, az); - amin = mod_a * (1.0-ltol/100.0); - amax = mod_a * (1.0+ltol/100.0); - - mod_b = modulus(bx, by, bz); - bmin = mod_b * (1.0-ltol/100.0); - bmax = mod_b * (1.0+ltol/100.0); - - mod_c = modulus(cx, cy, cz); - cmin = mod_c * (1.0-ltol/100.0); - cmax = mod_c * (1.0+ltol/100.0); - - s = malloc(3*sizeof(*s)); - s->pmax = pmax; - s->n_search = 3; - s->search = malloc(3*sizeof(struct reax_search_v)); - smin = 2.0*pmax * amin; smax = 2.0*pmax * amax; - s->search[0].smin = smin; s->search[0].smax = smax; - s->search[0].max_warned = 0; - smin = 2.0*pmax * bmin; smax = 2.0*pmax * bmax; - s->search[1].smin = smin; s->search[1].smax = smax; - s->search[1].max_warned = 0; - smin = 2.0*pmax * cmin; smax = 2.0*pmax * cmax; - s->search[2].smin = smin; s->search[2].smax = smax; - s->search[2].max_warned = 0; - - return s; -} - - -static double get_model_phase(double x, double y, double z, ImageFeatureList *f, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, const struct rigid_group *rg, - struct detector *det) -{ - struct dvec dir; - int s, i; - double max; - double re, im; - - dir.x = x; dir.y = y; dir.z = z; - - fill_and_transform(&dir, f, nel, pmax, fft_in, fft_out, plan, rg, det); - - s = -1; - max = 0.0; - for ( i=smin; i<=smax; i++ ) { - - double re, im, m; - - re = fft_out[i][0]; - im = fft_out[i][1]; - m = sqrt(re*re + im*im); - if ( m > max ) { - max = m; - s = i; - } - - } - - re = fft_out[s][0]; - im = fft_out[s][1]; - - return atan2(im, re); -} - - -static void refine_rigid_group(struct image *image, UnitCell *cell, - const struct rigid_group *rg, double pmax, - double *fft_in, fftw_complex *fft_out, - fftw_plan plan, int smin, int smax, - struct detector *det, struct reax_private *pr) -{ - double ax, ay, az, ma; - double bx, by, bz, mb; - double cx, cy, cz, mc; - double pha, phb, phc; - int i, j; - fftw_complex *r_fft_in; - fftw_complex *r_fft_out; - double m2m; - signed int aix, aiy; - signed int bix, biy; - signed int cix, ciy; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - ma = modulus(ax, ay, az); - mb = modulus(bx, by, bz); - mc = modulus(cx, cy, cz); - - pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, - pr->nel, pmax, fft_in, fft_out, plan, - smin, smax, rg, det); - phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, - pr->nel, pmax, fft_in, fft_out, plan, - smin, smax, rg, det); - phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, - pr->nel, pmax, fft_in, fft_out, plan, - smin, smax, rg, det); - - r_fft_in = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); - r_fft_out = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); - for ( i=0; i<pr->cw; i++ ) { - for ( j=0; j<pr->ch; j++ ) { - r_fft_in[i+pr->cw*j][0] = 0.0; - r_fft_in[i+pr->cw*j][1] = 0.0; - } - } - - ma = modulus(ax, ay, 0.0); - mb = modulus(bx, by, 0.0); - mc = modulus(cx, cy, 0.0); - m2m = ma; - if ( mb > m2m ) m2m = mb; - if ( mc > m2m ) m2m = mc; - - aix = (pr->cw/2)*ax/m2m; aiy = (pr->ch/2)*ay/m2m; - bix = (pr->cw/2)*bx/m2m; biy = (pr->ch/2)*by/m2m; - cix = (pr->cw/2)*cx/m2m; ciy = (pr->ch/2)*cy/m2m; - - if ( aix < 0 ) aix += pr->cw/2; - if ( bix < 0 ) bix += pr->cw/2; - if ( cix < 0 ) cix += pr->cw/2; - - if ( aiy < 0 ) aiy += pr->ch/2; - if ( biy < 0 ) biy += pr->ch/2; - if ( ciy < 0 ) ciy += pr->ch/2; - - r_fft_in[aix + pr->cw*aiy][0] = cos(pha); - r_fft_in[aix + pr->cw*aiy][1] = sin(pha); - r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][0] = cos(pha); - r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][1] = -sin(pha); - - r_fft_in[bix + pr->cw*biy][0] = cos(phb); - r_fft_in[bix + pr->cw*biy][1] = sin(phb); - r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][0] = cos(phb); - r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][1] = -sin(phb); - - r_fft_in[cix + pr->cw*ciy][0] = cos(phc); - r_fft_in[cix + pr->cw*ciy][1] = sin(phc); - r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][0] = cos(phc); - r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][1] = -sin(phc); - - const int tidx = 1; - r_fft_in[tidx][0] = 1.0; - r_fft_in[tidx][1] = 0.0; - -// STATUS("%i %i\n", aix, aiy); -// STATUS("%i %i\n", bix, biy); -// STATUS("%i %i\n", cix, ciy); - - fftw_execute_dft(pr->r_plan, r_fft_in, r_fft_out); - -// max = 0.0; -// FILE *fh = fopen("centering.dat", "w"); -// for ( i=0; i<pr->cw; i++ ) { -// for ( j=0; j<pr->ch; j++ ) { -// -// double re, im, am, ph; -// -// re = r_fft_out[i + pr->cw*j][0]; -// im = r_fft_out[i + pr->cw*j][1]; -// am = sqrt(re*re + im*im); -// ph = atan2(im, re); -// -// if ( am > max ) { -// max = am; -// max_i = i; -// max_j = j; -// } -// -// fprintf(fh, "%f ", am); -// -// } -// fprintf(fh, "\n"); -// } -// STATUS("Max at %i, %i\n", max_i, max_j); -// fclose(fh); -// exit(1); - -// STATUS("Offsets for '%s': %.2f, %.2f pixels\n", rg, dx, dy); -} - - -static UNUSED void refine_all_rigid_groups(struct image *image, UnitCell *cell, - double pmax, - double *fft_in, - fftw_complex *fft_out, - fftw_plan plan, int smin, int smax, - struct detector *det, - struct reax_private *p) -{ - int i; - - for ( i=0; i<image->det->n_rigid_groups; i++ ) { - refine_rigid_group(image, cell, image->det->rigid_groups[i], - pmax, fft_in, fft_out, plan, smin, smax, - det, p); - } -} - - -static double max_feature_resolution(ImageFeatureList *flist) -{ - double pmax; - int i, n; - - pmax = 0.0; - n = image_feature_count(flist); - for ( i=0; i<n; i++ ) { - - struct imagefeature *f; - double val; - - f = image_get_feature(flist, i); - if ( f == NULL ) continue; - - val = modulus(f->rx, f->ry, f->rz); - if ( val > pmax ) pmax = val; - - } - - return pmax; -} - - -static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c) -{ - struct rvec aCb; - double aCb_dot_c; - - /* "a" cross "b" */ - aCb.u = a.v*b.w - a.w*b.v; - aCb.v = - (a.u*b.w - a.w*b.u); - aCb.w = a.u*b.v - a.v*b.u; - - /* "a cross b" dot "c" */ - aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; - - if ( aCb_dot_c > 0.0 ) return 1; - return 0; -} - - -struct cell_candidate -{ - UnitCell *cell; - double fom; -}; - - -struct cell_candidate_list -{ - struct cell_candidate *cand; - int n_cand; -}; - - -static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose) -{ - int i; - int n_dup; - const int n_trials = 40; - - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - - gsl_rng *rng; - - /* This is a rubbish RNG, but it serves for this purpose: nothing more - * than "I couldn't be bothered to think of n_trials sets of random - * indices". */ - rng = gsl_rng_alloc(gsl_rng_taus2); - - cell_get_reciprocal(c1, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - cell_get_cartesian(c2, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - n_dup = 0; - for ( i=0; i<n_trials; i++ ) { - - signed int h, k, l; - double h2, k2, l2; - double rx, ry, rz; - double dev; - signed int h2i, k2i, l2i; - - h = gsl_rng_uniform_int(rng, 10); - k = gsl_rng_uniform_int(rng, 10); - l = gsl_rng_uniform_int(rng, 10); - - /* Position of this (randomly selected) - * reciprocal lattice point */ - rx = h*asx + k*bsx + l*csx; - ry = h*asy + k*bsy + l*csy; - rz = h*asz + k*bsz + l*csz; - - /* Indices of this point in the basis of the other cell */ - h2 = rx*ax + ry*ay + rz*az; - k2 = rx*bx + ry*by + rz*bz; - l2 = rx*cx + ry*cy + rz*cz; - - h2i = lrint(h2); - k2i = lrint(k2); - l2i = lrint(l2); - - dev = pow(h2i-h2, 2.0) + pow(k2i-k2, 2.0) + pow(l2i-l2, 2.0); - - if ( verbose ) { - STATUS("%3i %3i %3i -> %5.2f %5.2f %5.2f -> " - "%3i %3i %3i -> %5.2f\n", h, k, l, - h2, k2, l2, h2i, k2i, l2i, dev); - } - - if ( dev < 0.1 ) { - n_dup++; - } - - } - - if ( verbose ) { - STATUS("%i duplicates.\n", n_dup); - } - - gsl_rng_free(rng); - - if ( n_dup > 10 ) return 1; - return 0; -} - - -/* Return true if "cnew" accounts for more than 25% of the peaks predicted by - * any of the "ncells" cells in "cells". */ -static int twinned(UnitCell *cnew, struct cell_candidate_list *cl) -{ - int i; - - for ( i=0; i<cl->n_cand; i++ ) { - if ( check_twinning(cnew, cl->cand[i].cell, 0) ) return 1; - } - - return 0; -} - - -static int check_vector_combination(struct dvec *vi, struct dvec *vj, - struct dvec *vk, UnitCell *cell) -{ - double ang; - double a, b, c, al, be, ga; - const double angtol = deg2rad(5.0); - - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - - ang = angle_between(vi->x, vi->y, vi->z, vj->x, vj->y, vj->z); - if ( fabs(ang-ga) > angtol ) return 0; - - ang = angle_between(vi->x, vi->y, vi->z, vk->x, vk->y, vk->z); - if ( fabs(ang-be) > angtol ) return 0; - - ang = angle_between(vj->x, vj->y, vj->z, vk->x, vk->y, vk->z); - if ( fabs(ang-al) > angtol ) return 0; - - return 1; -} - - -static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, - double fom) -{ - struct cell_candidate cshift; - int i, cpos; - - cpos = cl->n_cand; - for ( i=0; i<cl->n_cand; i++ ) { - if ( fom > cl->cand[i].fom ) { - cpos = i; - break; - } - } - - cshift.cell = cnew; - cshift.fom = fom; - - for ( i=cpos; i<cl->n_cand; i++ ) { - - struct cell_candidate cshift2; - cshift2 = cl->cand[i]; - cl->cand[i] = cshift; - cshift = cshift2; - - } - - if ( cl->n_cand >= MAX_REAX_CELL_CANDIDATES ) { - /* "cshift" just fell off the end of the list */ - } else { - cl->cand[cl->n_cand++] = cshift; - } -} - - -static int check_cell(struct reax_private *dp, struct image *image, - UnitCell *cell) -{ - UnitCell *out; - Crystal *cr; - - 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 ( dp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - crystal_free(cr); /* Frees the cell as well */ - return 0; - } - } - - image_add_crystal(image, cr); - - return 1; -} - - -static int assemble_cells_from_candidates(struct image *image, - struct reax_search *s, - UnitCell *cell, - struct reax_private *p) -{ - int i, j, k; - signed int ti, tj, tk; - int rval; - struct cell_candidate_list cl; - - cl.cand = calloc(MAX_REAX_CELL_CANDIDATES, - sizeof(struct cell_candidate)); - if ( cl.cand == NULL ) { - ERROR("Failed to allocate cell candidate list.\n"); - return 0; - } - cl.n_cand = 0; - - /* Find candidates for axes 0 and 1 which have the right angle */ - for ( i=0; i<s->search[0].n_cand; i++ ) { - for ( j=0; j<s->search[1].n_cand; j++ ) { - for ( k=0; k<s->search[2].n_cand; k++ ) { - for ( ti=-1; ti<=1; ti+=2 ) { - for ( tj=-1; tj<=1; tj+=2 ) { - for ( tk=-1; tk<=1; tk+=2 ) { - - struct dvec vi, vj, vk; - struct rvec ai, bi, ci; - UnitCell *cnew; - double fom; - - vi = s->search[0].cand[i].v; - vj = s->search[1].cand[j].v; - vk = s->search[2].cand[k].v; - - vi.x *= ti; vi.y *= ti; vi.z *= ti; - vj.x *= tj; vj.y *= tj; vj.z *= tj; - vk.x *= tk; vk.y *= tk; vk.z *= tk; - - if ( !check_vector_combination(&vi, &vj, &vk, cell) ) continue; - - ai.u = vi.x; ai.v = vi.y; ai.w = vi.z; - bi.u = vj.x; bi.v = vj.y; bi.w = vj.z; - ci.u = vk.x; ci.v = vk.y; ci.w = vk.z; - - if ( !right_handed_vec(ai, bi, ci) ) continue; - - /* We have three vectors with the right angles */ - cnew = cell_new_from_direct_axes(ai, bi, ci); - - if ( twinned(cnew, &cl) ) { - cell_free(cnew); - continue; - } - - fom = 1.0; - add_cell_candidate(&cl, cnew, fom); - - } - } - } - } - } - } - - for ( i=0; i<cl.n_cand; i++ ) { - double a, b, c, al, be, ga; - double aA, bA, cA, alA, beA, gaA; - int w = 0; -// STATUS("%i: %f\n", i, cl.cand[i].fom); - cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga); - cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA, - &alA, &beA, &gaA); - if ( fabs(a - aA) > aA/10.0 ) w = 1; - if ( fabs(b - bA) > bA/10.0 ) w = 1; - if ( fabs(c - cA) > cA/10.0 ) w = 1; - if ( fabs(al - alA) > deg2rad(5.0) ) w = 1; - if ( fabs(be - beA) > deg2rad(5.0) ) w = 1; - if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1; - if ( w ) { - STATUS("This cell is a long way from that sought:\n"); - cell_print(cl.cand[i].cell); - } - } - - rval = 0; - for ( i=0; i<cl.n_cand; i++ ) { - if ( check_cell(p, image, cl.cand[i].cell) ) { - rval = 1; - break; - } - } - - free(cl.cand); - return rval; -} - - -int reax_index(IndexingPrivate *pp, struct image *image) -{ - struct reax_private *p; - double *fft_in; - fftw_complex *fft_out; - double pmax; - struct reax_search *s; - int i; - int rval = 0; - - p = (struct reax_private *)pp; - - fft_in = fftw_malloc(p->nel*sizeof(double)); - fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - pmax = max_feature_resolution(image->features); - - /* Sanity check */ - if ( pmax < 1e4 ) { - fftw_free(fft_in); - fftw_free(fft_out); - return 0; - } - - s = search_all_axes(p->cell, pmax); - find_candidates(p, image->features, pmax, fft_in, fft_out, s, - NULL, image->det); - -// refine_all_rigid_groups(image, image->candidate_cells[0], pmax, -// fft_in, fft_out, p->plan, smin, smax, -// image->det, p); - - rval = assemble_cells_from_candidates(image, s, p->cell, p); - - for ( i=0; i<s->n_search; i++ ) { - free(s->search[i].cand); - } - free(s->search); - free(s); - fftw_free(fft_in); - fftw_free(fft_out); - return rval; -} - - -IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) -{ - struct reax_private *p; - int samp; - double th; - - if ( cell == NULL ) { - ERROR("ReAx needs a unit cell.\n"); - return NULL; - } - - p = calloc(1, sizeof(*p)); - if ( p == NULL ) return NULL; - - p->cell = cell; - - /* Flags that ReAx knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS; - - /* Flags that ReAx requires */ - *indm |= INDEXING_USE_LATTICE_TYPE; - *indm |= INDEXING_USE_CELL_PARAMETERS; - - p->angular_inc = deg2rad(1.0); - - /* Reserve memory, over-estimating the number of directions */ - samp = 2.0*M_PI / p->angular_inc; - p->directions = malloc(samp*samp*sizeof(struct dvec)); - if ( p == NULL) { - free(p); - return NULL; - } - STATUS("Allocated space for %i directions\n", samp*samp); - - /* Generate vectors for 1D Fourier transforms */ - fesetround(1); /* Round to nearest */ - p->n_dir = 0; - for ( th=0.0; th<M_PI_2; th+=p->angular_inc ) { - - double ph, phstep, n_phstep; - - n_phstep = 2.0*M_PI*sin(th)/p->angular_inc; - n_phstep = nearbyint(n_phstep); - phstep = 2.0*M_PI/n_phstep; - - for ( ph=0.0; ph<2.0*M_PI; ph+=phstep ) { - - struct dvec *dir; - - assert(p->n_dir<samp*samp); - - dir = &p->directions[p->n_dir++]; - - dir->x = cos(ph) * sin(th); - dir->y = sin(ph) * sin(th); - dir->z = cos(th); - dir->th = th; - dir->ph = ph; - - } - - } - STATUS("Generated %i directions (angular increment %.3f deg)\n", - p->n_dir, rad2deg(p->angular_inc)); - - p->nel = 1024; - - /* These arrays are not actually used */ - p->fft_in = fftw_malloc(p->nel*sizeof(double)); - p->fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - p->plan = fftw_plan_dft_r2c_1d(p->nel, p->fft_in, p->fft_out, - FFTW_MEASURE); - - p->cw = 128; p->ch = 128; - - /* Also not used */ - p->r_fft_in = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); - p->r_fft_out = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); - - p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, - 1, FFTW_MEASURE); - - p->indm = *indm; - - return (IndexingPrivate *)p; -} - - -void reax_cleanup(IndexingPrivate *pp) -{ - struct reax_private *p; - - p = (struct reax_private *)pp; - - free(p->directions); - - fftw_destroy_plan(p->plan); - fftw_free(p->fft_in); - fftw_free(p->fft_out); - - fftw_destroy_plan(p->r_plan); - fftw_free(p->r_fft_in); - fftw_free(p->r_fft_out); - - free(p); -} diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h deleted file mode 100644 index 0bc96d9e..00000000 --- a/libcrystfel/src/reax.h +++ /dev/null @@ -1,75 +0,0 @@ -/* - * reax.h - * - * A new auto-indexer - * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2011-2012 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifndef REAX_H -#define REAX_H - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "index.h" -#include "cell.h" -#include "detector.h" - -#ifdef __cplusplus -extern "C" { -#endif - -#ifdef HAVE_FFTW - -extern IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl); - -extern void reax_cleanup(IndexingPrivate *pp); - -extern int reax_index(IndexingPrivate *pp, struct image *image); - -#else /* HAVE_FFTW */ - -static IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) -{ - return NULL; -} - -static void reax_cleanup(IndexingPrivate *pp) -{ -} - -static int reax_index(IndexingPrivate *pp, struct image *image) -{ -} - -#endif /* HAVE_FFTW */ - -#ifdef __cplusplus -} -#endif - -#endif /* REAX_H */ |