Implement Felix Indexer

9 files changed, 953 insertions, 9 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 715d26c0..52cab8a1 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -66,6 +66,12 @@ Invoke Mosflm.  To use this option, 'ipmosflm' must be in your shell's search pa
 .PD
 This is a implementation of the \fBdirax\fR algorithm, with some very small changes such as using a 1D Fourier transform for finding the lattice repeats.  This algorithm is implemented natively within CrystFEL meaning that no external software is required.
 
+.IP \fBfelix\fR
+.PD
+Invoke Felix, which will use your cell parameters to find multiple crystals in each pattern.
+.sp
+The Felix indexer has been developed by Soeren Schmidt <ssch@fysik.dtu.dk>. To use this option, 'Felix' must be in your shell's search path. This can be a link to the latest version of Felix. If you see the Felix version information when you run \fBFelix\fR on the command line, things are set up correctly.
+
 .IP \fBxds\fR
 .PD
 Invoke XDS, and use its REFIDX procedure to attempt to index the pattern.
@@ -119,7 +125,6 @@ It's risky to use \fBmosflm-nolatt\fR in conjunction with either \fB-comb\fR or
 
 If you don't know what to give for this option, try \fB--indexing=asdf,dirax-axes,mosflm-axes-latt,mosflm-axes-nolatt,xds\fR.
 
-
 .SH PEAK INTEGRATION
 If the pattern could be successfully indexed, peaks will be predicted in the pattern and their intensities measured.  You have a choice of integration methods, and you specify the method using \fB--integration\fR.  Choose from:
 
@@ -381,6 +386,118 @@ For a full explanation of how the internal layout of the data file can be  descr
 
 You can use \fBlist_events\fR to prepare a list of each event in one or more input files.  Note that you only need to do this if you need to perform some sorting or filtering on this list.  If you want to process every event in a file, simply specify the filename in the input file.
 
+.PD 0
+.IP \fB--felix-options\fR
+.PD
+Specify a comma-separated list of keyword arguments to change the default parameters passed to the felix indexer. For asthetics, this list can also be preceeded and followed by quotations - as in: \fB--felix-options="arg1=10,arg2=600"\fR. A list of parameters which can be modified through this option is detailed below.
+
+.PD 0
+.IP "\fBspacegroup=\fIn\fR"
+.PD
+Specify the spacegroup number of the crystal structure. For example: spacegroup=96 is P43212.
+
+.PD 0
+.IP "\fBn_voxels=\fIn\fR"
+.PD
+Specify the number of voxels used in the Rodrigues space search.
+(Default: 100)
+
+.PD 0
+.IP "\fBtest_fraction=\fIfract\fR"
+.PD
+Specify the fraction of space to test.
+(Default: 0.75)
+
+.PD 0
+.IP "\fBsigma_tth=\fIsig\fR"
+.PD
+Specify the sigma of the two-theta angle.
+(Default: 0.15)
+
+.PD 0
+.IP "\fBsigma_eta=\fIsig\fR"
+.PD
+Specify the sigma of the eta angle.
+(Default: 0.2)
+
+.PD 0
+.IP "\fBsigma_omega=\fIsig\fR"
+.PD
+Specify the sigma of the omega angle.
+(Default: 0.2)
+
+.PD 0
+.IP "\fBn_sigmas=\fIn\fR"
+.PD
+Specify the number of sigmas to use.
+(Default: 2)
+
+.PD 0
+.IP "\fBtthrange_min=\fImin\fR"
+.PD
+Specify the lower bound of the 2theta angle to consider when indexing (degrees).
+(Default: 0.0)
+
+.PD 0
+.IP "\fBtthrange_max=\fImax\fR"
+.PD
+Specify the upper bound of the 2theta angle to consider when indexing (degrees).
+(Default: 30.0)
+
+.PD 0
+.IP "\fBetarange_min=\fImin\fR"
+.PD
+Specify the lower bound of the eta angle to consider when indexing (degrees).
+(Default: 0.0)
+
+.PD 0
+.IP "\fBetarange_max=\fImax\fR"
+.PD
+Specify the upper bound of the eta angle to consider when indexing (degrees).
+(Default: 360.0)
+
+.PD 0
+.IP "\fBmin_measurements=\fIn\fR"
+.PD
+Specify a cutoff for how many Rodrigues lines cross in the first search.
+(Default: 15)
+
+.PD 0
+.IP "\fBmin_completeness=\fIfract\fR"
+.PD
+Specify a cutoff for the fraction of projected spots which are found in the pattern.
+(Default: 0.001)
+
+.PD 0
+.IP "\fBmin_uniqueness=\fIn\fR"
+.PD
+Specify a cutoff for the fraction of found spots which can belong to other crystallites.
+(Default: 0.5)
+
+.PD 0
+.IP "\fBforce4frustums\fR"
+.PD
+Tell Felix to use four frustums.
+
+.PD 0
+.IP "\fBorispace_octa\fR"
+.PD
+Tell Felix to use octahedral orispace as opposed to frustum.
+(Only possible with Felix.0.31)
+
+.PD 0
+.IP "\fBreadhkl=\fIfile\fR"
+.PD
+Specify a file from which to read in the ideal gvectors for the crystal.
+
+.PD 0
+.IP "\fBmaxtime=\fIn\fR"
+.PD
+Specify a maximum time which Felix is allowed to attempt to index the pattern (seconds)
+(Default: 30 s)
+
+
+
 .SH AUTHOR
 This page was written by Thomas White.
 
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index 77e5991b..f2760dac 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -10,7 +10,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.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/xds.c src/integration.c src/predict-refine.c \
-                         src/histogram.c src/events.c
+                         src/histogram.c src/events.c src/felix.c
 
 if HAVE_FFTW
 libcrystfel_la_SOURCES += src/asdf.c
@@ -30,7 +30,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \
                                  src/integer_matrix.h src/crystal.h \
                                  src/xds.h src/predict-refine.h \
                                  src/integration.h src/histogram.h \
-                                 src/events.h src/asdf.h
+                                 src/events.h src/asdf.h src/felix.h
 
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
 AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@
diff --git a/libcrystfel/src/felix.c b/libcrystfel/src/felix.c
new file mode 100644
index 00000000..2f416693
--- /dev/null
+++ b/libcrystfel/src/felix.c
@@ -0,0 +1,724 @@
+/*
+ * felix.c
+ *
+ * Invoke Felix for multi-crystal autoindexing.
+ *
+ * Copyright © 2015 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2015 Thomas White <taw@physics.org>
+ *   2015 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *
+ * 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 <sys/types.h>
+#include <sys/wait.h>
+#include <unistd.h>
+#include <assert.h>
+#include <fcntl.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 "felix.h"
+
+
+#define FELIX_VERBOSE 0
+#define MAX_MEAN_IA 0.25
+#define MIN_NGV 15
+#define MIN_NGV_UNIQUE 10
+
+
+/* Global private data, prepared once */
+struct felix_private
+{
+	IndexingMethod indm;
+	UnitCell *cell;
+
+	/* Options specific to Felix */
+	int spacegroup;
+	float tthrange_min;
+	float tthrange_max;
+	float etarange_min;
+	float etarange_max;
+	float domega;
+	float omegarange_min;
+	float omegarange_max;
+	int min_measurements;   /* related to "cuts" */
+	float min_completeness; /* related to "cuts" */
+	float min_uniqueness;   /* related to "cuts" */
+	int n_voxels;           /* related to "frustsumsize" */
+	float test_fraction;    /* related to "frustsumsize" */
+	float sigma_tth;        /* related to "uncertainties" */
+	float sigma_eta;        /* related to "uncertainties" */
+	float sigma_omega;      /* related to "uncertainties" */
+	int n_sigmas;
+	int force4frustums;
+
+	/*Felix v0.3 options*/
+	int orispace_frustum;
+	int orispace_octa;
+	char *readhkl_file;
+	float maxtime;
+
+};
+
+
+/* Data needed to call Felix */
+struct felix_data {
+
+	struct felix_private *gp;
+
+	/* Low-level stuff */
+	int                     pty;
+	pid_t                   pid;
+	char                    *rbuffer;
+	int                     rbufpos;
+	int                     rbuflen;
+
+};
+
+
+static int read_felix(struct felix_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;
+	float mean_ia;
+	/* int ngv_unique;  /* Used in Felix 0.1 */
+	int ngv;
+	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 *.felix file.\n");
+		return 1;
+	}
+
+	do {
+
+		Crystal *cr;
+		UnitCell *cell;
+
+		/* One line per grain */
+		if ( fgets( line, 1024, fh ) == NULL ) {
+			break;
+		}
+		/* File format of the .felix files
+		 * version 0.1 - 0.2
+		 *
+		 * r = sscanf(line, "%i %f %i %i %f %f %f %f %f %f %f %f %f"
+		 *	         "%f %f %f %f %f %f %f %f %f %f %f %f",
+		 *	         &d1, &mean_ia, &ngv, &ngv_unique, &d2, &d2, &d2,
+		 *	         &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2,
+		 *         &ubi11, &ubi12, &ubi13,
+		 *	         &ubi21, &ubi22, &ubi23,
+		 *	         &ubi31, &ubi32, &ubi33);
+		 *
+		 * if ( r != 25 ) {
+		 *		ERROR("Only %i parameters in .felix file\n", r);
+		 *		return 1;
+		 * }
+		 */
+
+		/* version 0.3 - present */
+		r = sscanf(line, "%i %f %i %f %f %f %f %f %f"
+					 "%f %f %f %f %f %f %f %f %f %f %f %f",
+					 &d1, &mean_ia, &ngv, &d2, &d2,
+					 &d2, &d2, &d2, &d2, &d2, &d2, &d2,
+					 &ubi11, &ubi12, &ubi13,
+					 &ubi21, &ubi22, &ubi23,
+					 &ubi31, &ubi32, &ubi33);
+
+		if ( r != 21 ) {
+			ERROR("Only %i parameters in .felix file, check version and format.\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);
+
+		/* Poor indexing criterion for Felix v0.1
+		 *
+		 * if (mean_ia > MAX_MEAN_IA || ngv < MIN_NGV ||
+		 * 		ngv_unique < MIN_NGV_UNIQUE ){
+		 * 		crystal_set_user_flag(cr, 1);
+		 */
+
+		/* Poor indexing criterion for Felix v0.3 */
+
+		if ( mean_ia > MAX_MEAN_IA || ngv < MIN_NGV ){
+			crystal_set_user_flag(cr, 1);
+		}
+
+		/* All crystals are saved to the image,
+		 * but only good ones will be later written to the stream file.
+		 */
+
+		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 felix_data *gs)
+{
+	#if FELIX_VERBOSE
+	STATUS("%s\n", line);
+	#endif
+}
+
+
+static int felix_readable(struct image *image, struct felix_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 felix_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, struct felix_private *gp)
+{
+	FILE *fh;
+	char *filename;
+	char gveFilename[1024];
+	char logFilename[1024];
+	double tt;
+
+	filename = malloc(1024);
+	if ( filename == NULL ) return NULL;
+
+	snprintf(filename, 1023, "xfel-%i.ini", image->id);
+	snprintf(gveFilename, 1023, "xfel-%i.gve", image->id);
+	snprintf(logFilename, 1023, "xfel-%i.log", 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 %i\n", gp->spacegroup);
+	fprintf(fh, "tthrange %f %f\n", gp->tthrange_min, gp->tthrange_max);
+	fprintf(fh, "etarange %f %f\n", gp->etarange_min, gp->etarange_max);
+	fprintf(fh, "domega %f\n", gp->domega);
+	fprintf(fh, "omegarange %f %f\n", gp->omegarange_min, gp->omegarange_max);
+	fprintf(fh, "filespecs %s %s\n", gveFilename, logFilename);
+	fprintf(fh, "cuts %i %f %f\n", gp->min_measurements, gp->min_completeness,
+			gp->min_uniqueness);
+	fprintf(fh, "frustumsize %i %f\n", gp->n_voxels,
+			gp->test_fraction);
+	fprintf(fh, "uncertainties %f %f %f\n", gp->sigma_tth,
+			gp->sigma_eta, gp->sigma_omega);
+	fprintf(fh, "nsigmas %i\n", gp->n_sigmas);
+
+	if ( gp->force4frustums == 1 ){
+		fprintf(fh, "force4frustums\n");
+	}
+
+	if ( gp->orispace_frustum == 1 ){
+		fprintf(fh, "orispace frustum\n");
+	} else if ( gp->orispace_octa ==1 ){
+		fprintf(fh, "orispace octa\n");
+	} else{
+		ERROR("No felix supported orispace specified.\n");
+		free(filename);
+		filename = NULL;
+	}
+
+	/* If an hkl file is not specified, generate the peak list. */
+	if ( gp->readhkl_file != NULL ){
+		fprintf(fh, "readhkl %s\n", gp->readhkl_file);
+	} else{
+		fprintf(fh, "genhkl\n");
+	}
+
+	fprintf(fh, "maxtime %f\n", gp->maxtime);
+
+	fclose(fh);
+
+	return filename;
+}
+
+
+int felix_index(struct image *image, IndexingPrivate *ipriv)
+{
+	unsigned int opts;
+	int status;
+	int rval;
+	struct felix_data *felix;
+	struct felix_private *gp = (struct felix_private *) ipriv;
+	char *ini_filename;
+	char gff_filename[1024];
+
+	write_gve (image, gp);
+	ini_filename = write_ini (image, gp);
+
+	if ( ini_filename == NULL ) {
+		ERROR("Failed to write ini file for Felix.\n");
+		return 0;
+	}
+
+	felix = malloc(sizeof(struct felix_data));
+	if ( felix == NULL ) {
+		ERROR("Couldn't allocate memory for Felix data.\n");
+		return 0;
+	}
+
+	felix->gp = gp;
+
+	snprintf(gff_filename, 1023, "xfel-%i.felix", image->id);
+	remove(gff_filename);
+
+	felix->pid = forkpty(&felix->pty, NULL, NULL, NULL);
+	if ( felix->pid == -1 ) {
+		ERROR("Failed to fork for Felix: %s\n", strerror(errno));
+		return 0;
+	}
+	if ( felix->pid == 0 ) {
+
+		/* Child process: invoke Felix */
+		struct termios t;
+
+		/* Turn echo off */
+		tcgetattr(STDIN_FILENO, &t);
+		t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+		tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+		STATUS("Running Felix '%s'\n", ini_filename);
+		execlp("Felix", "", ini_filename, (char *)NULL);
+		ERROR("Failed to invoke Felix.\n");
+		_exit(0);
+
+	}
+
+	free(ini_filename);
+
+	felix->rbuffer = malloc(256);
+	felix->rbuflen = 256;
+	felix->rbufpos = 0;
+
+	/* Set non-blocking */
+	opts = fcntl(felix->pty, F_GETFL);
+	fcntl(felix->pty, F_SETFL, opts | O_NONBLOCK);
+
+	do {
+
+		fd_set fds;
+		struct timeval tv;
+		int sval;
+
+		FD_ZERO(&fds);
+		FD_SET(felix->pty, &fds);
+
+		tv.tv_sec = gp->maxtime + 1.0;
+		tv.tv_usec = 0;
+
+		sval = select(felix->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 = felix_readable(image, felix);
+		} else {
+			ERROR("No response from Felix..\n");
+			rval = 1;
+		}
+
+	} while ( !rval );
+
+	close(felix->pty);
+	free(felix->rbuffer);
+	waitpid(felix->pid, &status, 0);
+
+	if ( status != 0 ) {
+		ERROR("Felix either timed out, or is not working properly.\n");
+		free(felix);
+		return 0;
+	}
+
+	if ( read_felix(gp, image, gff_filename) != 0 ) {
+		free(felix);
+		return 0;
+	}
+
+	/* Success! */
+	free(felix);
+	return 1;
+
+}
+
+IndexingPrivate *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+                                      struct detector *det, float *ltl,
+                                      const char *options)
+{
+	struct felix_private *gp;
+	char *temp_options=NULL;
+	char *option=NULL;
+
+	if ( cell == NULL ) {
+		ERROR("Felix needs a unit cell.\n");
+		return NULL;
+	}
+
+	if ( options != NULL ){
+		temp_options = strdup(options);
+	}
+
+	gp = calloc(1, sizeof(*gp));
+	if ( gp == NULL ) return NULL;
+
+	/* Flags that Felix knows about */
+	*indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS
+	       | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS;
+
+	gp->cell = cell;
+	gp->indm = *indm;
+
+	/* Default values of felix options */
+	gp->spacegroup = 0;
+	gp->tthrange_min = 0;
+	gp->tthrange_max = 30.0;
+	gp->etarange_min = 0;
+	gp->etarange_max = 360;
+	gp->domega = 2;
+	gp->omegarange_min = -1.0;
+	gp->omegarange_max = 1.0;
+	gp->min_measurements = 15;
+	gp->min_completeness = 0.001;
+	gp->min_uniqueness = 0.5;
+	gp->n_voxels = 100;
+	gp->test_fraction = 0.75;
+	gp->sigma_tth = 0.15;
+	gp->sigma_eta = 0.2;
+	gp->sigma_omega = 0.2;
+	gp->n_sigmas = 2;
+	gp->force4frustums = 0;
+	gp->orispace_frustum = 1;
+	gp->orispace_octa = 0;
+	gp->readhkl_file = NULL;
+	gp->maxtime = 30.0;
+
+	/*Parse the options string and fill in the necessary private variables. */
+	if (temp_options != NULL){
+	while((option=strsep(&temp_options, ","))!=NULL) {
+
+		if ( strncmp(option, "spacegroup=", 11) == 0 ){
+			gp->spacegroup = atoi(option+11);
+		}
+
+		if ( strncmp(option, "tthrange_min=", 13) == 0 ){
+			gp->tthrange_min = atof(option+13);
+		}
+
+		if ( strncmp(option, "tthrange_max=", 13 ) == 0 ){
+			gp->tthrange_max = atof(option+13);
+		}
+
+		if ( strncmp(option, "etarange_min=", 13 ) == 0 ){
+			gp->etarange_min = atof(option+13);
+		}
+
+		if ( strncmp(option, "etarange_max=", 13 ) == 0 ){
+			gp->etarange_max = atof(option+13);
+		}
+
+		if ( strncmp(option, "domega=", 7) == 0 ){
+			gp->domega = atof(option+7);
+		}
+
+		if ( strncmp(option, "omegarange_min=", 15) == 0 ){
+			gp->omegarange_min = atof(option+15);
+		}
+
+		if ( strncmp(option, "omegarange_max=", 15) == 0 ){
+			gp->omegarange_max = atof(option+15);
+		}
+
+		if ( strncmp(option, "min_measurements=", 17) == 0 ){
+			gp->min_measurements = atoi(option+17);
+		}
+
+		if ( strncmp(option, "min_completeness=", 17) == 0 ){
+			gp->min_completeness = atof(option+17);
+		}
+
+		if ( strncmp(option, "min_uniqueness=", 15) == 0 ){
+			gp->min_uniqueness = atof(option+15);
+		}
+
+		if ( strncmp(option, "n_voxels=", 9) == 0 ){
+			gp->n_voxels = atoi(option+9);
+		}
+
+		if ( strncmp(option, "test_fraction=", 14) == 0 ){
+			gp->test_fraction = atof(option+14);
+		}
+
+		if ( strncmp(option, "sigma_tth=", 10) == 0 ){
+			gp->sigma_tth = atof(option+10);
+		}
+
+		if ( strncmp(option, "sigma_eta=", 10) == 0 ){
+			gp->sigma_eta = atof(option+10);
+		}
+
+		if ( strncmp(option, "sigma_omega=", 12) == 0 ){
+			gp->sigma_omega = atof(option+12);
+		}
+
+		if ( strncmp(option, "n_sigmas=", 9) == 0 ){
+			gp->n_sigmas = atoi(option+9);
+		}
+
+		if ( strncmp(option, "force4frustums", 14) == 0 ){
+			gp->force4frustums = 1;
+		}
+
+		/* Only allow one orispace command.
+		 * orispace frustum is the default, so have to turn off to use octa.
+		 */
+		if ( strncmp(option, "orispace_octa", 13) == 0 ){
+			gp->orispace_octa = 1;
+			gp->orispace_frustum = 0;
+		}
+
+		if ( strncmp(option, "readhkl=", 8) == 0 ){
+			gp->readhkl_file = strdup(option+8);
+		}
+
+		if ( strncmp(option, "maxtime=", 8) == 0 ){
+			gp->maxtime = atof(option+8);
+		}
+
+	}
+	}
+
+	/* Make sure that they at least specified the spacegroup number.*/
+
+	if ( gp->spacegroup == 0 ){
+		ERROR("Felix requires that you specify the spacegroup number.\n");
+		ERROR("You should use the argument --felix-options=spacegroup=xx\n");
+		return NULL;
+	}
+
+	free(temp_options);
+	free(option);
+
+	return (IndexingPrivate *) gp;
+}
+
+
+void felix_cleanup(IndexingPrivate *pp)
+{
+	struct felix_private *p;
+
+	p = (struct felix_private *) pp;
+	free(p->readhkl_file);
+	free(p);
+}
diff --git a/libcrystfel/src/felix.h b/libcrystfel/src/felix.h
new file mode 100644
index 00000000..40771933
--- /dev/null
+++ b/libcrystfel/src/felix.h
@@ -0,0 +1,50 @@
+/*
+ * felix.h
+ *
+ * Invoke Felix 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>
+ *   2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *
+ * 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 FELIX_H
+#define FELIX_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "cell.h"
+
+extern IndexingPrivate *felix_prepare(IndexingMethod *indm,
+                                             UnitCell *cell,
+                                             struct detector *det,
+                                             float *ltl,
+                                             const char *options);
+
+extern void felix_cleanup(IndexingPrivate *pp);
+
+extern int felix_index(struct image *image, IndexingPrivate *p);
+
+
+#endif	/* FELIX_H */
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 944c3f10..b409e80e 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -12,6 +12,7 @@
  *   2010-2011 Richard Kirian <rkirian@asu.edu>
  *   2012      Lorenzo Galli
  *   2013      Cornelius Gati <cornelius.gati@cfel.de>
+ *   2015	   Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -51,6 +52,7 @@
 #include "index.h"
 #include "geometry.h"
 #include "cell-utils.h"
+#include "felix.h"
 
 
 static int debug_index(struct image *image)
@@ -67,7 +69,8 @@ static int debug_index(struct image *image)
 
 
 IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
-                                   struct detector *det, float *ltl)
+                                   struct detector *det, float *ltl,
+                                   const char *options)
 {
 	int n;
 	int nm = 0;
@@ -106,6 +109,10 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
 			iprivs[n] = (IndexingPrivate *)strdup("Hello!");
 			break;
 
+			case INDEXING_FELIX :
+			iprivs[n] = felix_prepare(&indm[n], cell, det, ltl, options);
+			break;
+
 			default :
 			ERROR("Don't know how to prepare indexing method %i\n",
 			      indm[n]);
@@ -170,10 +177,14 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
 			mosflm_cleanup(privs[n]);
 			break;
 
-                        case INDEXING_XDS :
+			case INDEXING_XDS :
 			xds_cleanup(privs[n]);
 			break;
 
+			case INDEXING_FELIX :
+			felix_cleanup(privs[n]);
+			break;
+
 			case INDEXING_DEBUG :
 			free(privs[n]);
 			break;
@@ -244,6 +255,11 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 
 		case INDEXING_DEBUG :
 		return debug_index(image);
+		break;
+
+		case INDEXING_FELIX :
+		return felix_index(image, ipriv);
+		break;
 
 		default :
 		ERROR("Unrecognised indexing method: %i\n", indm);
@@ -366,6 +382,10 @@ char *indexer_str(IndexingMethod indm)
 		strcpy(str, "mosflm");
 		break;
 
+		case INDEXING_FELIX :
+		strcpy(str, "felix");
+		break;
+
 		case INDEXING_XDS :
 		strcpy(str, "xds");
 		break;
@@ -438,9 +458,12 @@ IndexingMethod *build_indexer_list(const char *str)
 		} else if ( strcmp(methods[i], "mosflm") == 0) {
 			list[++nmeth] = INDEXING_DEFAULTS_MOSFLM;
 
-        } else if ( strcmp(methods[i], "xds") == 0) {
+		} else if ( strcmp(methods[i], "xds") == 0) {
 			list[++nmeth] = INDEXING_DEFAULTS_XDS;
 
+		} else if ( strcmp(methods[i], "felix") == 0) {
+			list[++nmeth] = INDEXING_DEFAULTS_FELIX;
+
 		} else if ( strcmp(methods[i], "none") == 0) {
 			list[++nmeth] = INDEXING_NONE;
 
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index d835d003..0ce74e37 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -12,6 +12,7 @@
  *   2010-2015 Thomas White <taw@physics.org>
  *   2010      Richard Kirian
  *   2012      Lorenzo Galli
+ *   2015      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -49,6 +50,10 @@
                                      | INDEXING_USE_LATTICE_TYPE               \
                                      | INDEXING_USE_CELL_PARAMETERS)
 
+#define INDEXING_DEFAULTS_FELIX (INDEXING_FELIX                  \
+                                     | INDEXING_USE_LATTICE_TYPE               \
+                                     | INDEXING_USE_CELL_PARAMETERS)
+
 /* 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            \
@@ -60,6 +65,7 @@
  * @INDEXING_NONE: No indexing to be performed
  * @INDEXING_DIRAX: Invoke DirAx
  * @INDEXING_MOSFLM: Invoke MOSFLM
+ * @INDEXING_FELIX: Invoke Felix
  * @INDEXING_XDS: Invoke XDS
  * @INDEXING_SIMULATION: Dummy value
  * @INDEXING_DEBUG: Results injector for debugging
@@ -86,6 +92,7 @@ typedef enum {
 	/* The core indexing methods themselves */
 	INDEXING_DIRAX  = 1,
 	INDEXING_MOSFLM = 2,
+	INDEXING_FELIX = 4,
 	INDEXING_XDS = 5,
 	INDEXING_SIMULATION = 6,
 	INDEXING_DEBUG = 7,
@@ -125,7 +132,8 @@ extern char *indexer_str(IndexingMethod indm);
 #include "image.h"
 
 extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
-                                          struct detector *det, float *ltl);
+                                          struct detector *det, float *ltl,
+                                          const char *options);
 
 extern void index_pattern(struct image *image,
                           IndexingMethod *indms, IndexingPrivate **iprivs);
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 48d5004e..503078a7 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -256,7 +256,6 @@ IntegerMatrix *get_symop(const SymOpList *ops, const SymOpMask *m, int idx)
 	return ops->ops[idx];
 }
 
-
 static signed int *v(signed int h, signed int k, signed int i, signed int l)
 {
 	signed int *vec = malloc(3*sizeof(signed int));
diff --git a/src/indexamajig.c b/src/indexamajig.c
index e5b6e19b..164fa071 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -140,6 +140,11 @@ static void show_help(const char *s)
 "     --no-refls-in-stream Do not record integrated reflections in the stream.\n"
 "     --int-diag=<cond>    Show debugging information about reflections.\n"
 "     --no-refine          Skip the prediction refinement step.\n"
+"\nLow-level options for the felix indexer:\n\n"
+"     --felix-options      Change the default arguments passed to the indexer.\n"
+"                          Given as a list of comma separated list of \n"
+"                            indexer specific, key word arguments.\n"
+"                          Example: \"arg1=10,arg2=500\" \n"
 );
 }
 
@@ -192,6 +197,7 @@ int main(int argc, char *argv[])
 	char *geom_filename = NULL;
 	struct beam_params beam;
 	int have_push_res = 0;
+	int len;
 
 	/* Defaults */
 	iargs.cell = NULL;
@@ -236,6 +242,7 @@ int main(int argc, char *argv[])
 	iargs.fix_bandwidth = -1.0;
 	iargs.fix_divergence = -1.0;
 	iargs.predict_refine = 1;
+	iargs.felix_options = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -294,6 +301,7 @@ int main(int argc, char *argv[])
 		{"fix-profile-radius", 1, NULL,               22},
 		{"fix-bandwidth",      1, NULL,               23},
 		{"fix-divergence",     1, NULL,               24},
+		{"felix-options",      1, NULL,               25},
 
 		{0, 0, NULL, 0}
 	};
@@ -463,6 +471,20 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 25 :
+			/* Remove leading and trailing quotes */
+			len = strlen(optarg);
+			if ( optarg[len-1] == '\'' || optarg[len-1] == '\"' ){
+				optarg[len-1] = 0;
+			}
+			if ( optarg[0] == '\'' || optarg[0] == '\"' ){
+				iargs.felix_options = strdup( optarg+1 );
+			}
+			else {
+				iargs.felix_options = strdup( optarg );
+			}
+			break;
+
 			case 0 :
 			break;
 
@@ -712,7 +734,7 @@ int main(int argc, char *argv[])
 	/* Prepare the indexer */
 	if ( indm != NULL ) {
 		ipriv = prepare_indexing(indm, iargs.cell, iargs.det,
-		                         iargs.tols);
+		                         iargs.tols, iargs.felix_options);
 		if ( ipriv == NULL ) {
 			ERROR("Failed to prepare indexing.\n");
 			return 1;
diff --git a/src/process_image.h b/src/process_image.h
index 0d5ca1b7..60da7cb8 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -89,6 +89,7 @@ struct index_args
 	float fix_bandwidth;
 	float fix_divergence;
 	int predict_refine;
+	char *felix_options;
 };