path: root/libcrystfel/src/index.c

/*
 * index.c
 *
 * Perform indexing (somehow)
 *
 * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
 *                       a research centre of the Helmholtz Association.
 * Copyright © 2012 Lorenzo Galli
 *
 * Authors:
 *   2010-2013 Thomas White <taw@physics.org>
 *   2010-2011 Richard Kirian <rkirian@asu.edu>
 *   2012      Lorenzo Galli
 *
 * 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 <string.h>
#include <assert.h>

#include "image.h"
#include "utils.h"
#include "peaks.h"
#include "dirax.h"
#include "mosflm.h"
#include "detector.h"
#include "index.h"
#include "index-priv.h"
#include "reax.h"
#include "geometry.h"
#include "cell-utils.h"


IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
                                   const char *filename, struct detector *det,
                                   struct beam_params *beam, float *ltl)
{
	int n;
	int nm = 0;
	IndexingPrivate **iprivs;

	while ( indm[nm] != INDEXING_NONE ) nm++;
	iprivs = malloc((nm+1) * sizeof(IndexingPrivate *));

	for ( n=0; n<nm; n++ ) {

		switch ( indm[n] & INDEXING_METHOD_MASK ) {

			case INDEXING_DIRAX :
			iprivs[n] = dirax_prepare(indm[n], cell, filename, det,
			                          beam, ltl);
			break;

			case INDEXING_MOSFLM :
			iprivs[n] = mosflm_prepare(indm[n], cell, filename, det,
			                          beam, ltl);
			break;

			case INDEXING_REAX :
			iprivs[n] = reax_prepare(cell, filename, det, beam);
			break;

			default :
			ERROR("Don't know how to prepare indexing method %i\n",
			      indm[n]);
			break;

		}

	}
	iprivs[n] = NULL;

	return iprivs;
}


void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
{
	int n = 0;

	if ( indms == NULL ) return;  /* Nothing to do */
	if ( privs == NULL ) return;  /* Nothing to do */

	while ( indms[n] != INDEXING_NONE ) {

		switch ( indms[n] & INDEXING_METHOD_MASK ) {

			case INDEXING_NONE :
			case INDEXING_DIRAX :
			case INDEXING_MOSFLM :
			/* No cleanup */
			/* FIXME: Not true */
			break;

			case INDEXING_REAX :
			reax_cleanup(privs[n]);
			break;

			default :
			ERROR("Don't know how to clean up indexing method %i\n",
			      indms[n]);
			break;

		}

		n++;

	}
}


void map_all_peaks(struct image *image)
{
	int i, n;

	n = image_feature_count(image->features);

	/* Map positions to 3D */
	for ( i=0; i<n; i++ ) {

		struct imagefeature *f;
		struct rvec r;

		f = image_get_feature(image->features, i);
		if ( f == NULL ) continue;

		r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
		f->rx = r.u;  f->ry = r.v;  f->rz = r.w;

	}
}


/* Return non-zero for "success" */
static int try_indexer(struct image *image, IndexingMethod indm,
                        IndexingPrivate *ipriv)
{
	switch ( indm & INDEXING_METHOD_MASK ) {

		case INDEXING_NONE :
		break;

		case INDEXING_DIRAX :
		return run_dirax(image, ipriv);
		break;

		case INDEXING_MOSFLM :
		return run_mosflm(image, ipriv);
		break;

		case INDEXING_REAX :
		return reax_index(image, ipriv);
		break;

		default :
		ERROR("Unrecognised indexing method: %i\n", indm);
		break;

	}

	return 0;
}


void index_pattern(struct image *image,
                   IndexingMethod *indms, IndexingPrivate **iprivs)
{
	int n = 0;

	if ( indms == NULL ) return;

	map_all_peaks(image);
	image->crystals = NULL;
	image->n_crystals = 0;

	while ( indms[n] != INDEXING_NONE ) {

		if ( try_indexer(image, indms[n], iprivs[n]) ) break;
		n++;

	}
}


/* Set the default indexer flags.  May need tweaking depending on the method */
static IndexingMethod defaults(IndexingMethod a)
{
	return a | INDEXING_CHECK_CELL_COMBINATIONS | INDEXING_CHECK_PEAKS;
}


/* Set the indexer flags for "raw mode" ("--cell-reduction=none") */
static IndexingMethod set_raw(IndexingMethod a)
{
	/* Disable all unit cell checks */
	a &= ~(INDEXING_CHECK_CELL_COMBINATIONS | INDEXING_CHECK_CELL_AXES);
	return a;
}


/* Set the indexer flags for "bad mode" ("--insane) */
static IndexingMethod set_bad(IndexingMethod a)
{
	/* Disable the peak check */
	return a & ~INDEXING_CHECK_PEAKS;
}


/* Set the indexer flags for "axes mode" ("--cell-reduction=compare") */
static IndexingMethod set_axes(IndexingMethod a)
{
	return (a & ~INDEXING_CHECK_CELL_COMBINATIONS)
	          | INDEXING_CHECK_CELL_AXES;
}


IndexingMethod *build_indexer_list(const char *str)
{
	int n, i;
	char **methods;
	IndexingMethod *list;
	int nmeth = 0;

	n = assplode(str, ",-", &methods, ASSPLODE_NONE);
	list = malloc((n+1)*sizeof(IndexingMethod));

	nmeth = -1;  /* So that the first method is #0 */
	for ( i=0; i<n; i++ ) {

		if ( strcmp(methods[i], "dirax") == 0) {
			list[++nmeth] = defaults(INDEXING_DIRAX);

		} else if ( strcmp(methods[i], "mosflm") == 0) {
			list[++nmeth] = defaults(INDEXING_MOSFLM);

		} else if ( strcmp(methods[i], "reax") == 0) {
			/* ReAx doesn't need any cell check */
			list[++nmeth] = set_raw(defaults(INDEXING_REAX));

		} else if ( strcmp(methods[i], "raw") == 0) {
			list[nmeth] = set_raw(list[nmeth]);

		} else if ( strcmp(methods[i], "bad") == 0) {
			list[nmeth] = set_bad(list[nmeth]);

		} else if ( strcmp(methods[i], "axes") == 0) {
			list[nmeth] = set_axes(list[nmeth]);

		} else {
			ERROR("Bad list of indexing methods: '%s'\n", str);
			return NULL;
		}

		free(methods[i]);

	}
	free(methods);
	list[++nmeth] = INDEXING_NONE;

	return list;
}