aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/cl-utils.c200
-rw-r--r--libcrystfel/src/cl-utils.h27
-rw-r--r--libcrystfel/src/diffraction-gpu.c529
-rw-r--r--libcrystfel/src/diffraction-gpu.h57
-rw-r--r--libcrystfel/src/diffraction.c463
-rw-r--r--libcrystfel/src/diffraction.h34
-rw-r--r--libcrystfel/src/filters.c130
-rw-r--r--libcrystfel/src/filters.h25
-rw-r--r--libcrystfel/src/geometry.c341
-rw-r--r--libcrystfel/src/geometry.h26
-rw-r--r--libcrystfel/src/peaks.c593
-rw-r--r--libcrystfel/src/peaks.h38
-rw-r--r--libcrystfel/src/reflist-utils.c502
-rw-r--r--libcrystfel/src/reflist-utils.h52
-rw-r--r--libcrystfel/src/statistics.c668
-rw-r--r--libcrystfel/src/statistics.h45
-rw-r--r--libcrystfel/src/stream.c487
-rw-r--r--libcrystfel/src/stream.h49
-rw-r--r--libcrystfel/src/symmetry.c1503
-rw-r--r--libcrystfel/src/symmetry.h63
20 files changed, 5832 insertions, 0 deletions
diff --git a/libcrystfel/src/cl-utils.c b/libcrystfel/src/cl-utils.c
new file mode 100644
index 00000000..65a09363
--- /dev/null
+++ b/libcrystfel/src/cl-utils.c
@@ -0,0 +1,200 @@
+/*
+ * cl-utils.c
+ *
+ * OpenCL utility functions
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#ifdef HAVE_CL_CL_H
+#include <CL/cl.h>
+#else
+#include <cl.h>
+#endif
+
+#include "utils.h"
+
+
+const char *clError(cl_int err)
+{
+ switch ( err ) {
+ case CL_SUCCESS : return "no error";
+ case CL_INVALID_PLATFORM : return "invalid platform";
+ case CL_INVALID_KERNEL : return "invalid kernel";
+ case CL_INVALID_ARG_INDEX : return "invalid argument index";
+ case CL_INVALID_ARG_VALUE : return "invalid argument value";
+ case CL_INVALID_MEM_OBJECT : return "invalid memory object";
+ case CL_INVALID_SAMPLER : return "invalid sampler";
+ case CL_INVALID_ARG_SIZE : return "invalid argument size";
+ case CL_INVALID_COMMAND_QUEUE : return "invalid command queue";
+ case CL_INVALID_CONTEXT : return "invalid context";
+ case CL_INVALID_VALUE : return "invalid value";
+ case CL_INVALID_EVENT_WAIT_LIST : return "invalid wait list";
+ case CL_MAP_FAILURE : return "map failure";
+ case CL_MEM_OBJECT_ALLOCATION_FAILURE : return "object allocation failure";
+ case CL_OUT_OF_HOST_MEMORY : return "out of host memory";
+ case CL_OUT_OF_RESOURCES : return "out of resources";
+ case CL_INVALID_KERNEL_NAME : return "invalid kernel name";
+ case CL_INVALID_KERNEL_ARGS : return "invalid kernel arguments";
+ case CL_INVALID_WORK_GROUP_SIZE : return "invalid work group size";
+ case CL_IMAGE_FORMAT_NOT_SUPPORTED : return "image format not supported";
+ case CL_INVALID_WORK_DIMENSION : return "invalid work dimension";
+ default :
+ return "unknown error";
+ }
+}
+
+
+static char *get_device_string(cl_device_id dev, cl_device_info info)
+{
+ int r;
+ size_t size;
+ char *val;
+
+ r = clGetDeviceInfo(dev, info, 0, NULL, &size);
+ if ( r != CL_SUCCESS ) {
+ ERROR("Couldn't get device vendor size: %s\n",
+ clError(r));
+ return NULL;
+ }
+ val = malloc(size);
+ r = clGetDeviceInfo(dev, info, size, val, NULL);
+ if ( r != CL_SUCCESS ) {
+ ERROR("Couldn't get dev vendor: %s\n", clError(r));
+ return NULL;
+ }
+
+ return val;
+}
+
+
+cl_device_id get_cl_dev(cl_context ctx, int n)
+{
+ cl_device_id *dev;
+ cl_int r;
+ size_t size;
+ int i, num_devs;
+
+ /* Get the required size of the array */
+ r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, 0, NULL, &size);
+ if ( r != CL_SUCCESS ) {
+ ERROR("Couldn't get array size for devices: %s\n", clError(r));
+ return 0;
+ }
+
+ dev = malloc(size);
+ r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, size, dev, NULL);
+ if ( r != CL_SUCCESS ) {
+ ERROR("Couldn't get device: %s\n", clError(r));
+ return 0;
+ }
+ num_devs = size/sizeof(cl_device_id);
+
+ if ( n >= num_devs ) {
+ ERROR("Device ID out of range\n");
+ return 0;
+ }
+
+ if ( n < 0 ) {
+
+ STATUS("Available devices:\n");
+ for ( i=0; i<num_devs; i++ ) {
+
+ char *vendor;
+ char *name;
+
+ vendor = get_device_string(dev[i], CL_DEVICE_VENDOR);
+ name = get_device_string(dev[i], CL_DEVICE_NAME);
+
+ STATUS("Device %i: %s %s\n", i, vendor, name);
+
+ }
+ n = 0;
+
+ STATUS("Using device 0. Use --gpu-dev to choose another.\n");
+
+ } else {
+
+ char *vendor;
+ char *name;
+
+ vendor = get_device_string(dev[n], CL_DEVICE_VENDOR);
+ name = get_device_string(dev[n], CL_DEVICE_NAME);
+
+ STATUS("Using device %i: %s %s\n", n, vendor, name);
+
+ }
+
+ return dev[n];
+}
+
+
+static void show_build_log(cl_program prog, cl_device_id dev)
+{
+ cl_int r;
+ char log[4096];
+ size_t s;
+
+ r = clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, 4096, log,
+ &s);
+
+ STATUS("%s\n", log);
+}
+
+
+cl_program load_program(const char *filename, cl_context ctx,
+ cl_device_id dev, cl_int *err, const char *extra_cflags)
+{
+ FILE *fh;
+ cl_program prog;
+ char *source;
+ size_t len;
+ cl_int r;
+ char cflags[1024] = "";
+
+ fh = fopen(filename, "r");
+ if ( fh == NULL ) {
+ ERROR("Couldn't open '%s'\n", filename);
+ *err = CL_INVALID_PROGRAM;
+ return 0;
+ }
+ source = malloc(16384);
+ len = fread(source, 1, 16383, fh);
+ fclose(fh);
+ source[len] = '\0';
+
+ prog = clCreateProgramWithSource(ctx, 1, (const char **)&source,
+ NULL, err);
+ if ( *err != CL_SUCCESS ) {
+ ERROR("Couldn't load source\n");
+ return 0;
+ }
+
+ snprintf(cflags, 1023, "-Werror ");
+ strncat(cflags, "-I"DATADIR"/crystfel/ ", 1023-strlen(cflags));
+ strncat(cflags, "-cl-no-signed-zeros ", 1023-strlen(cflags));
+ strncat(cflags, extra_cflags, 1023-strlen(cflags));
+
+ r = clBuildProgram(prog, 0, NULL, cflags, NULL, NULL);
+ if ( r != CL_SUCCESS ) {
+ ERROR("Couldn't build program '%s'\n", filename);
+ show_build_log(prog, dev);
+ *err = r;
+ return 0;
+ }
+
+ free(source);
+ *err = CL_SUCCESS;
+ return prog;
+}
diff --git a/libcrystfel/src/cl-utils.h b/libcrystfel/src/cl-utils.h
new file mode 100644
index 00000000..21a7ecd2
--- /dev/null
+++ b/libcrystfel/src/cl-utils.h
@@ -0,0 +1,27 @@
+/*
+ * cl-utils.h
+ *
+ * OpenCL utility functions
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef CLUTILS_H
+#define CLUTILS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+extern const char *clError(cl_int err);
+extern cl_device_id get_cl_dev(cl_context ctx, int n);
+extern cl_program load_program(const char *filename, cl_context ctx,
+ cl_device_id dev, cl_int *err,
+ const char *extra_cflags);
+
+
+#endif /* CLUTILS_H */
diff --git a/libcrystfel/src/diffraction-gpu.c b/libcrystfel/src/diffraction-gpu.c
new file mode 100644
index 00000000..605b1514
--- /dev/null
+++ b/libcrystfel/src/diffraction-gpu.c
@@ -0,0 +1,529 @@
+/*
+ * diffraction-gpu.c
+ *
+ * Calculate diffraction patterns by Fourier methods (GPU version)
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include <complex.h>
+
+#ifdef HAVE_CL_CL_H
+#include <CL/cl.h>
+#else
+#include <cl.h>
+#endif
+
+#include "image.h"
+#include "utils.h"
+#include "cell.h"
+#include "diffraction.h"
+#include "cl-utils.h"
+#include "beam-parameters.h"
+
+
+#define SAMPLING (4)
+#define BWSAMPLING (10)
+#define DIVSAMPLING (1)
+#define SINC_LUT_ELEMENTS (4096)
+
+
+struct gpu_context
+{
+ cl_context ctx;
+ cl_command_queue cq;
+ cl_program prog;
+ cl_kernel kern;
+ cl_mem intensities;
+ cl_mem flags;
+
+ /* Array of sinc LUTs */
+ cl_mem *sinc_luts;
+ cl_float **sinc_lut_ptrs;
+ int max_sinc_lut; /* Number of LUTs, i.e. one greater than the maximum
+ * index. This equals the highest allowable "n". */
+};
+
+
+static void check_sinc_lut(struct gpu_context *gctx, int n)
+{
+ cl_int err;
+ cl_image_format fmt;
+ int i;
+
+ if ( n > gctx->max_sinc_lut ) {
+
+ gctx->sinc_luts = realloc(gctx->sinc_luts,
+ n*sizeof(*gctx->sinc_luts));
+ gctx->sinc_lut_ptrs = realloc(gctx->sinc_lut_ptrs,
+ n*sizeof(*gctx->sinc_lut_ptrs));
+
+ for ( i=gctx->max_sinc_lut; i<n; i++ ) {
+ gctx->sinc_lut_ptrs[i] = NULL;
+ }
+
+ gctx->max_sinc_lut = n;
+ }
+
+ if ( gctx->sinc_lut_ptrs[n-1] != NULL ) return;
+
+ /* Create a new sinc LUT */
+ gctx->sinc_lut_ptrs[n-1] = malloc(SINC_LUT_ELEMENTS*sizeof(cl_float));
+ gctx->sinc_lut_ptrs[n-1][0] = n;
+ if ( n == 1 ) {
+ for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
+ gctx->sinc_lut_ptrs[n-1][i] = 1.0;
+ }
+ } else {
+ for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
+ double x, val;
+ x = (double)i/SINC_LUT_ELEMENTS;
+ val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
+ gctx->sinc_lut_ptrs[n-1][i] = val;
+ }
+ }
+
+ fmt.image_channel_order = CL_INTENSITY;
+ fmt.image_channel_data_type = CL_FLOAT;
+
+ gctx->sinc_luts[n-1] = clCreateImage2D(gctx->ctx,
+ CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
+ &fmt, SINC_LUT_ELEMENTS, 1, 0,
+ gctx->sinc_lut_ptrs[n-1], &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't create LUT for %i\n", n);
+ return;
+ }
+}
+
+
+static int set_arg_float(struct gpu_context *gctx, int idx, float val)
+{
+ cl_int err;
+ err = clSetKernelArg(gctx->kern, idx, sizeof(cl_float), &val);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't set kernel argument %i: %s\n",
+ idx, clError(err));
+ return 1;
+ }
+
+ return 0;
+}
+
+
+static int set_arg_int(struct gpu_context *gctx, int idx, int val)
+{
+ cl_int err;
+
+ err = clSetKernelArg(gctx->kern, idx, sizeof(cl_int), &val);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't set kernel argument %i: %s\n",
+ idx, clError(err));
+ return 1;
+ }
+
+ return 0;
+}
+
+
+static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val)
+{
+ cl_int err;
+
+ err = clSetKernelArg(gctx->kern, idx, sizeof(cl_mem), &val);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't set kernel argument %i: %s\n",
+ idx, clError(err));
+ return 1;
+ }
+
+ return 0;
+}
+
+
+void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
+ int na, int nb, int nc, UnitCell *ucell)
+{
+ cl_int err;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ float klow, khigh;
+ int i;
+ cl_float16 cell;
+ cl_int4 ncells;
+ const int sampling = SAMPLING;
+ cl_float bwstep;
+ int n_inf = 0;
+ int n_neg = 0;
+ cl_float divxlow, divxstep;
+ cl_float divylow, divystep;
+ int n_nan = 0;
+ int sprod;
+
+ if ( gctx == NULL ) {
+ ERROR("GPU setup failed.\n");
+ return;
+ }
+
+ cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az;
+ cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz;
+ cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz;
+
+ /* Calculate wavelength */
+ klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0));
+ khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0));
+ bwstep = (khigh-klow) / BWSAMPLING;
+
+ /* Calculate divergence stuff */
+ divxlow = -image->beam->divergence/2.0;
+ divylow = -image->beam->divergence/2.0;
+ divxstep = image->beam->divergence / DIVSAMPLING;
+ divystep = image->beam->divergence / DIVSAMPLING;
+
+ ncells.s[0] = na;
+ ncells.s[1] = nb;
+ ncells.s[2] = nc;
+ ncells.s[3] = 0; /* unused */
+
+ /* Ensure all required LUTs are available */
+ check_sinc_lut(gctx, na);
+ check_sinc_lut(gctx, nb);
+ check_sinc_lut(gctx, nc);
+
+ if ( set_arg_float(gctx, 2, klow) ) return;
+ if ( set_arg_mem(gctx, 9, gctx->intensities) ) return;
+ if ( set_arg_int(gctx, 12, sampling) ) return;
+ if ( set_arg_float(gctx, 14, bwstep) ) return;
+ if ( set_arg_mem(gctx, 15, gctx->sinc_luts[na-1]) ) return;
+ if ( set_arg_mem(gctx, 16, gctx->sinc_luts[nb-1]) ) return;
+ if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return;
+ if ( set_arg_mem(gctx, 18, gctx->flags) ) return;
+ if ( set_arg_float(gctx, 23, divxlow) ) return;
+ if ( set_arg_float(gctx, 24, divxstep) ) return;
+ if ( set_arg_int(gctx, 25, DIVSAMPLING) ) return;
+ if ( set_arg_float(gctx, 26, divylow) ) return;
+ if ( set_arg_float(gctx, 27, divystep) ) return;
+ if ( set_arg_int(gctx, 28, DIVSAMPLING) ) return;
+
+ /* Unit cell */
+ err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't set unit cell: %s\n", clError(err));
+ return;
+ }
+
+ /* Local memory for reduction */
+ sprod = BWSAMPLING*SAMPLING*SAMPLING*DIVSAMPLING*DIVSAMPLING;
+ err = clSetKernelArg(gctx->kern, 13, sprod*sizeof(cl_float), NULL);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't set local memory: %s\n", clError(err));
+ return;
+ }
+
+ /* Allocate memory for the result */
+ image->data = calloc(image->width * image->height, sizeof(float));
+ image->twotheta = calloc(image->width * image->height, sizeof(double));
+
+ /* Iterate over panels */
+ for ( i=0; i<image->det->n_panels; i++ ) {
+
+ size_t dims[3];
+ size_t ldims[3] = {SAMPLING, SAMPLING,
+ BWSAMPLING * DIVSAMPLING * DIVSAMPLING};
+ struct panel *p;
+ cl_mem tt;
+ size_t tt_size;
+ cl_mem diff;
+ size_t diff_size;
+ float *diff_ptr;
+ float *tt_ptr;
+ int pan_width, pan_height;
+ int fs, ss;
+
+ p = &image->det->panels[i];
+
+ pan_width = 1 + p->max_fs - p->min_fs;
+ pan_height = 1 + p->max_ss - p->min_ss;
+
+ /* Buffer for the results of this panel */
+ diff_size = pan_width * pan_height * sizeof(cl_float);
+ diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY,
+ diff_size, NULL, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't allocate diffraction memory\n");
+ return;
+ }
+ tt_size = pan_width * pan_height * sizeof(cl_float);
+ tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size,
+ NULL, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't allocate twotheta memory\n");
+ return;
+ }
+
+ if ( set_arg_mem(gctx, 0, diff) ) return;
+ if ( set_arg_mem(gctx, 1, tt) ) return;
+ if ( set_arg_int(gctx, 3, pan_width) ) return;
+ if ( set_arg_float(gctx, 4, p->cnx) ) return;
+ if ( set_arg_float(gctx, 5, p->cny) ) return;
+ if ( set_arg_float(gctx, 6, p->res) ) return;
+ if ( set_arg_float(gctx, 7, p->clen) ) return;
+ if ( set_arg_int(gctx, 10, p->min_fs) ) return;
+ if ( set_arg_int(gctx, 11, p->min_ss) ) return;
+ if ( set_arg_float(gctx, 19, p->fsx) ) return;
+ if ( set_arg_float(gctx, 20, p->fsy) ) return;
+ if ( set_arg_float(gctx, 21, p->ssx) ) return;
+ if ( set_arg_float(gctx, 22, p->ssy) ) return;
+
+ dims[0] = pan_width * SAMPLING;
+ dims[1] = pan_height * SAMPLING;
+ dims[2] = BWSAMPLING * DIVSAMPLING * DIVSAMPLING;
+
+ err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 3, NULL,
+ dims, ldims, 0, NULL, NULL);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't enqueue diffraction kernel: %s\n",
+ clError(err));
+ return;
+ }
+
+ clFinish(gctx->cq);
+
+ diff_ptr = clEnqueueMapBuffer(gctx->cq, diff, CL_TRUE,
+ CL_MAP_READ, 0, diff_size,
+ 0, NULL, NULL, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't map diffraction buffer: %s\n",
+ clError(err));
+ return;
+ }
+ tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ,
+ 0, tt_size, 0, NULL, NULL, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't map tt buffer\n");
+ return;
+ }
+
+ for ( fs=0; fs<pan_width; fs++ ) {
+ for ( ss=0; ss<pan_height; ss++ ) {
+
+ float val, tt;
+ int tfs, tss;
+
+ val = diff_ptr[fs + pan_width*ss];
+ if ( isinf(val) ) n_inf++;
+ if ( val < 0.0 ) n_neg++;
+ if ( isnan(val) ) n_nan++;
+ tt = tt_ptr[fs + pan_width*ss];
+
+ tfs = p->min_fs + fs;
+ tss = p->min_ss + ss;
+ image->data[tfs + image->width*tss] = val;
+ image->twotheta[tfs + image->width*tss] = tt;
+
+ }
+ }
+
+ clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr,
+ 0, NULL, NULL);
+ clEnqueueUnmapMemObject(gctx->cq, tt, tt_ptr,
+ 0, NULL, NULL);
+
+ clReleaseMemObject(diff);
+ clReleaseMemObject(tt);
+
+ }
+
+
+ if ( n_neg + n_inf + n_nan ) {
+ ERROR("WARNING: The GPU calculation produced %i negative"
+ " values, %i infinities and %i NaNs.\n",
+ n_neg, n_inf, n_nan);
+ }
+
+}
+
+
+/* Setup the OpenCL stuff, create buffers, load the structure factor table */
+struct gpu_context *setup_gpu(int no_sfac,
+ const double *intensities, unsigned char *flags,
+ const char *sym, int dev_num)
+{
+ struct gpu_context *gctx;
+ cl_uint nplat;
+ cl_platform_id platforms[8];
+ cl_context_properties prop[3];
+ cl_int err;
+ cl_device_id dev;
+ size_t intensities_size;
+ float *intensities_ptr;
+ size_t flags_size;
+ float *flags_ptr;
+ size_t maxwgsize;
+ int i;
+ char cflags[512] = "";
+
+ STATUS("Setting up GPU...\n");
+
+ err = clGetPlatformIDs(8, platforms, &nplat);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't get platform IDs: %i\n", err);
+ return NULL;
+ }
+ if ( nplat == 0 ) {
+ ERROR("Couldn't find at least one platform!\n");
+ return NULL;
+ }
+ prop[0] = CL_CONTEXT_PLATFORM;
+ prop[1] = (cl_context_properties)platforms[0];
+ prop[2] = 0;
+
+ gctx = malloc(sizeof(*gctx));
+ gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
+ NULL, NULL, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't create OpenCL context: %i\n", err);
+ free(gctx);
+ return NULL;
+ }
+
+ dev = get_cl_dev(gctx->ctx, dev_num);
+
+ gctx->cq = clCreateCommandQueue(gctx->ctx, dev, 0, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't create OpenCL command queue\n");
+ free(gctx);
+ return NULL;
+ }
+
+ /* Create a single-precision version of the scattering factors */
+ intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float);
+ intensities_ptr = malloc(intensities_size);
+ if ( intensities != NULL ) {
+ for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+ intensities_ptr[i] = intensities[i];
+ }
+ } else {
+ for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+ intensities_ptr[i] = 1e5;
+ }
+ strncat(cflags, "-DFLAT_INTENSITIES ", 511-strlen(cflags));
+ }
+ gctx->intensities = clCreateBuffer(gctx->ctx,
+ CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
+ intensities_size, intensities_ptr, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't allocate intensities memory\n");
+ free(gctx);
+ return NULL;
+ }
+ free(intensities_ptr);
+
+ if ( sym != NULL ) {
+ if ( strcmp(sym, "1") == 0 ) {
+ strncat(cflags, "-DPG1 ", 511-strlen(cflags));
+ } else if ( strcmp(sym, "-1") == 0 ) {
+ strncat(cflags, "-DPG1BAR ", 511-strlen(cflags));
+ } else if ( strcmp(sym, "6/mmm") == 0 ) {
+ strncat(cflags, "-DPG6MMM ", 511-strlen(cflags));
+ } else if ( strcmp(sym, "6") == 0 ) {
+ strncat(cflags, "-DPG6 ", 511-strlen(cflags));
+ } else if ( strcmp(sym, "6/m") == 0 ) {
+ strncat(cflags, "-DPG6M ", 511-strlen(cflags));
+ } else {
+ ERROR("Sorry! Point group '%s' is not currently"
+ " supported on the GPU."
+ " I'm using '1' instead.\n", sym);
+ strncat(cflags, "-DPG1 ", 511-strlen(cflags));
+ }
+ } else {
+ if ( intensities != NULL ) {
+ ERROR("You gave me an intensities file but no point"
+ " group. I'm assuming '1'.\n");
+ strncat(cflags, "-DPG1 ", 511-strlen(cflags));
+ }
+ }
+
+ /* Create a flag array */
+ flags_size = IDIM*IDIM*IDIM*sizeof(cl_float);
+ flags_ptr = malloc(flags_size);
+ if ( flags != NULL ) {
+ for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+ flags_ptr[i] = flags[i];
+ }
+ } else {
+ for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+ flags_ptr[i] = 1.0;
+ }
+ }
+ gctx->flags = clCreateBuffer(gctx->ctx,
+ CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
+ flags_size, flags_ptr, &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't allocate flag buffer\n");
+ free(gctx);
+ return NULL;
+ }
+ free(flags_ptr);
+
+ gctx->prog = load_program(DATADIR"/crystfel/diffraction.cl", gctx->ctx,
+ dev, &err, cflags);
+ if ( err != CL_SUCCESS ) {
+ free(gctx);
+ return NULL;
+ }
+
+ gctx->kern = clCreateKernel(gctx->prog, "diffraction", &err);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't create kernel\n");
+ free(gctx);
+ return NULL;
+ }
+
+ gctx->max_sinc_lut = 0;
+ gctx->sinc_lut_ptrs = NULL;
+ gctx->sinc_luts = NULL;
+
+ clGetDeviceInfo(dev, CL_DEVICE_MAX_WORK_GROUP_SIZE,
+ sizeof(size_t), &maxwgsize, NULL);
+ STATUS("Maximum work group size = %lli\n", (long long int)maxwgsize);
+
+ return gctx;
+}
+
+
+void cleanup_gpu(struct gpu_context *gctx)
+{
+ int i;
+
+ clReleaseProgram(gctx->prog);
+ clReleaseMemObject(gctx->intensities);
+
+ /* Release LUTs */
+ for ( i=1; i<=gctx->max_sinc_lut; i++ ) {
+ if ( gctx->sinc_lut_ptrs[i-1] != NULL ) {
+ clReleaseMemObject(gctx->sinc_luts[i-1]);
+ free(gctx->sinc_lut_ptrs[i-1]);
+ }
+ }
+
+ free(gctx->sinc_luts);
+ free(gctx->sinc_lut_ptrs);
+
+ clReleaseCommandQueue(gctx->cq);
+ clReleaseContext(gctx->ctx);
+ free(gctx);
+}
diff --git a/libcrystfel/src/diffraction-gpu.h b/libcrystfel/src/diffraction-gpu.h
new file mode 100644
index 00000000..a3bde4e1
--- /dev/null
+++ b/libcrystfel/src/diffraction-gpu.h
@@ -0,0 +1,57 @@
+/*
+ * diffraction-gpu.h
+ *
+ * Calculate diffraction patterns by Fourier methods (GPU version)
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef DIFFRACTION_GPU_H
+#define DIFFRACTION_GPU_H
+
+#include "image.h"
+#include "cell.h"
+
+struct gpu_context;
+
+#if HAVE_OPENCL
+
+extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
+ int na, int nb, int nc, UnitCell *ucell);
+extern struct gpu_context *setup_gpu(int no_sfac,
+ const double *intensities,
+ const unsigned char *flags,
+ const char *sym, int dev_num);
+extern void cleanup_gpu(struct gpu_context *gctx);
+
+#else
+
+static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
+ int na, int nb, int nc, UnitCell *ucell)
+{
+ /* Do nothing */
+ ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n");
+}
+
+static struct gpu_context *setup_gpu(int no_sfac,
+ const double *intensities,
+ const unsigned char *flags,
+ const char *sym, int dev_num)
+{
+ return NULL;
+}
+
+static void cleanup_gpu(struct gpu_context *gctx)
+{
+}
+
+#endif
+
+#endif /* DIFFRACTION_GPU_H */
diff --git a/libcrystfel/src/diffraction.c b/libcrystfel/src/diffraction.c
new file mode 100644
index 00000000..9532a6ce
--- /dev/null
+++ b/libcrystfel/src/diffraction.c
@@ -0,0 +1,463 @@
+/*
+ * diffraction.c
+ *
+ * Calculate diffraction patterns by Fourier methods
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include <complex.h>
+#include <assert.h>
+#include <fenv.h>
+
+#include "image.h"
+#include "utils.h"
+#include "cell.h"
+#include "diffraction.h"
+#include "beam-parameters.h"
+#include "symmetry.h"
+
+
+#define SAMPLING (4)
+#define BWSAMPLING (10)
+#define DIVSAMPLING (1)
+#define SINC_LUT_ELEMENTS (4096)
+
+
+static double *get_sinc_lut(int n)
+{
+ int i;
+ double *lut;
+
+ lut = malloc(SINC_LUT_ELEMENTS*sizeof(double));
+ lut[0] = n;
+ if ( n == 1 ) {
+ for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
+ lut[i] = 1.0;
+ }
+ } else {
+ for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
+ double x, val;
+ x = (double)i/SINC_LUT_ELEMENTS;
+ val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
+ lut[i] = val;
+ }
+ }
+
+ return lut;
+}
+
+
+static double interpolate_lut(double *lut, double val)
+{
+ double i, pos, f;
+ unsigned int low, high;
+
+ pos = SINC_LUT_ELEMENTS * modf(fabs(val), &i);
+ low = (int)pos; /* Discard fractional part */
+ high = low + 1;
+ f = modf(pos, &i); /* Fraction */
+ if ( high == SINC_LUT_ELEMENTS ) high = 0;
+
+ return (1.0-f)*lut[low] + f*lut[high];
+}
+
+
+static double lattice_factor(struct rvec q, double ax, double ay, double az,
+ double bx, double by, double bz,
+ double cx, double cy, double cz,
+ double *lut_a, double *lut_b,
+ double *lut_c)
+{
+ struct rvec Udotq;
+ double f1, f2, f3;
+
+ Udotq.u = ax*q.u + ay*q.v + az*q.w;
+ Udotq.v = bx*q.u + by*q.v + bz*q.w;
+ Udotq.w = cx*q.u + cy*q.v + cz*q.w;
+
+ f1 = interpolate_lut(lut_a, Udotq.u);
+ f2 = interpolate_lut(lut_b, Udotq.v);
+ f3 = interpolate_lut(lut_c, Udotq.w);
+
+ return f1 * f2 * f3;
+}
+
+
+static double sym_lookup_intensity(const double *intensities,
+ const unsigned char *flags,
+ const SymOpList *sym,
+ signed int h, signed int k, signed int l)
+{
+ int i;
+ double ret = 0.0;
+
+ for ( i=0; i<num_equivs(sym, NULL); i++ ) {
+
+ signed int he;
+ signed int ke;
+ signed int le;
+ double f, val;
+
+ get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
+
+ f = (double)lookup_flag(flags, he, ke, le);
+ val = lookup_intensity(intensities, he, ke, le);
+
+ ret += f*val;
+
+ }
+
+ return ret;
+}
+
+
+static double sym_lookup_phase(const double *phases,
+ const unsigned char *flags, const SymOpList *sym,
+ signed int h, signed int k, signed int l)
+{
+ int i;
+ double ret = 0.0;
+
+ for ( i=0; i<num_equivs(sym, NULL); i++ ) {
+
+ signed int he;
+ signed int ke;
+ signed int le;
+ double f, val;
+
+ get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
+
+ f = (double)lookup_flag(flags, he, ke, le);
+ val = lookup_phase(phases, he, ke, le);
+
+ ret += f*val;
+
+ }
+
+ return ret;
+}
+
+
+static double interpolate_linear(const double *ref, const unsigned char *flags,
+ const SymOpList *sym, float hd,
+ signed int k, signed int l)
+{
+ signed int h;
+ double val1, val2;
+ float f;
+
+ h = (signed int)hd;
+ if ( hd < 0.0 ) h -= 1;
+ f = hd - (float)h;
+ assert(f >= 0.0);
+
+ val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
+ val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
+
+ val1 = val1;
+ val2 = val2;
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_bilinear(const double *ref,
+ const unsigned char *flags,
+ const SymOpList *sym,
+ float hd, float kd, signed int l)
+{
+ signed int k;
+ double val1, val2;
+ float f;
+
+ k = (signed int)kd;
+ if ( kd < 0.0 ) k -= 1;
+ f = kd - (float)k;
+ assert(f >= 0.0);
+
+ val1 = interpolate_linear(ref, flags, sym, hd, k, l);
+ val2 = interpolate_linear(ref, flags, sym, hd, k+1, l);
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_intensity(const double *ref,
+ const unsigned char *flags,
+ const SymOpList *sym,
+ float hd, float kd, float ld)
+{
+ signed int l;
+ double val1, val2;
+ float f;
+
+ l = (signed int)ld;
+ if ( ld < 0.0 ) l -= 1;
+ f = ld - (float)l;
+ assert(f >= 0.0);
+
+ val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l);
+ val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1);
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
+static double complex interpolate_phased_linear(const double *ref,
+ const double *phases,
+ const unsigned char *flags,
+ const SymOpList *sym,
+ float hd,
+ signed int k, signed int l)
+{
+ signed int h;
+ double val1, val2;
+ float f;
+ double ph1, ph2;
+ double re1, re2, im1, im2;
+ double re, im;
+
+ h = (signed int)hd;
+ if ( hd < 0.0 ) h -= 1;
+ f = hd - (float)h;
+ assert(f >= 0.0);
+
+ val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
+ val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
+ ph1 = sym_lookup_phase(phases, flags, sym, h, k, l);
+ ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l);
+
+ val1 = val1;
+ val2 = val2;
+
+ /* Calculate real and imaginary parts */
+ re1 = val1 * cos(ph1);
+ im1 = val1 * sin(ph1);
+ re2 = val2 * cos(ph2);
+ im2 = val2 * sin(ph2);
+
+ re = (1.0-f)*re1 + f*re2;
+ im = (1.0-f)*im1 + f*im2;
+
+ return re + im*I;
+}
+
+
+static double complex interpolate_phased_bilinear(const double *ref,
+ const double *phases,
+ const unsigned char *flags,
+ const SymOpList *sym,
+ float hd, float kd,
+ signed int l)
+{
+ signed int k;
+ double complex val1, val2;
+ float f;
+
+ k = (signed int)kd;
+ if ( kd < 0.0 ) k -= 1;
+ f = kd - (float)k;
+ assert(f >= 0.0);
+
+ val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l);
+ val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l);
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_phased_intensity(const double *ref,
+ const double *phases,
+ const unsigned char *flags,
+ const SymOpList *sym,
+ float hd, float kd, float ld)
+{
+ signed int l;
+ double complex val1, val2;
+ float f;
+
+ l = (signed int)ld;
+ if ( ld < 0.0 ) l -= 1;
+ f = ld - (float)l;
+ assert(f >= 0.0);
+
+ val1 = interpolate_phased_bilinear(ref, phases, flags, sym,
+ hd, kd, l);
+ val2 = interpolate_phased_bilinear(ref, phases, flags, sym,
+ hd, kd, l+1);
+
+ return cabs((1.0-f)*val1 + f*val2);
+}
+
+
+/* Look up the structure factor for the nearest Bragg condition */
+static double molecule_factor(const double *intensities, const double *phases,
+ const unsigned char *flags, struct rvec q,
+ double ax, double ay, double az,
+ double bx, double by, double bz,
+ double cx, double cy, double cz,
+ GradientMethod m, const SymOpList *sym)
+{
+ float hd, kd, ld;
+ signed int h, k, l;
+ double r;
+
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * cz;
+
+ /* No flags -> flat intensity distribution */
+ if ( flags == NULL ) return 1.0e5;
+
+ switch ( m ) {
+ case GRADIENT_MOSAIC :
+ fesetround(1); /* Round to nearest */
+ h = (signed int)rint(hd);
+ k = (signed int)rint(kd);
+ l = (signed int)rint(ld);
+ if ( abs(h) > INDMAX ) r = 0.0;
+ else if ( abs(k) > INDMAX ) r = 0.0;
+ else if ( abs(l) > INDMAX ) r = 0.0;
+ else r = sym_lookup_intensity(intensities, flags, sym, h, k, l);
+ break;
+ case GRADIENT_INTERPOLATE :
+ r = interpolate_intensity(intensities, flags, sym, hd, kd, ld);
+ break;
+ case GRADIENT_PHASED :
+ r = interpolate_phased_intensity(intensities, phases, flags,
+ sym, hd, kd, ld);
+ break;
+ default:
+ ERROR("This gradient method not implemented yet.\n");
+ exit(1);
+ }
+
+ return r;
+}
+
+
+void get_diffraction(struct image *image, int na, int nb, int nc,
+ const double *intensities, const double *phases,
+ const unsigned char *flags, UnitCell *cell,
+ GradientMethod m, const SymOpList *sym)
+{
+ unsigned int fs, ss;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ float klow, khigh, bwstep;
+ double *lut_a;
+ double *lut_b;
+ double *lut_c;
+ double divxlow, divylow, divxstep, divystep;
+
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ /* Allocate (and zero) the "diffraction array" */
+ image->data = calloc(image->width * image->height, sizeof(float));
+
+ /* Needed later for Lorentz calculation */
+ image->twotheta = malloc(image->width * image->height * sizeof(double));
+
+ klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0));
+ khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0));
+ bwstep = (khigh-klow) / BWSAMPLING;
+
+ divxlow = -image->beam->divergence/2.0;
+ divylow = -image->beam->divergence/2.0;
+ divxstep = image->beam->divergence / DIVSAMPLING;
+ divystep = image->beam->divergence / DIVSAMPLING;
+
+ lut_a = get_sinc_lut(na);
+ lut_b = get_sinc_lut(nb);
+ lut_c = get_sinc_lut(nc);
+
+ for ( fs=0; fs<image->width; fs++ ) {
+ for ( ss=0; ss<image->height; ss++ ) {
+
+ int fs_step, ss_step, kstep;
+ int divxval, divyval;
+ int idx = fs + image->width*ss;
+
+ for ( fs_step=0; fs_step<SAMPLING; fs_step++ ) {
+ for ( ss_step=0; ss_step<SAMPLING; ss_step++ ) {
+ for ( kstep=0; kstep<BWSAMPLING; kstep++ ) {
+ for ( divxval=0; divxval<DIVSAMPLING; divxval++ ) {
+ for ( divyval=0; divyval<DIVSAMPLING; divyval++ ) {
+
+ double k;
+ double intensity;
+ double f_lattice, I_lattice;
+ double I_molecule;
+ struct rvec q, qn;
+ double twotheta;
+ const double dfs = (double)fs
+ + ((double)fs_step / SAMPLING);
+ const double dss = (double)ss
+ + ((double)ss_step / SAMPLING);
+
+ double xdiv = divxlow + divxstep*(double)divxval;
+ double ydiv = divylow + divystep*(double)divyval;
+
+ /* Calculate k this time round */
+ k = klow + (double)kstep * bwstep;
+
+ qn = get_q(image, dfs, dss, &twotheta, k);
+
+ /* x divergence */
+ q.u = qn.u*cos(xdiv) +qn.w*sin(xdiv);
+ q.v = qn.v;
+ q.w = -qn.u*sin(xdiv) +qn.w*cos(xdiv);
+
+ qn = q;
+
+ /* y divergence */
+ q.v = qn.v*cos(ydiv) +qn.w*sin(ydiv);
+ q.w = -qn.v*sin(ydiv) +qn.w*cos(ydiv);
+
+ f_lattice = lattice_factor(q, ax, ay, az,
+ bx, by, bz,
+ cx, cy, cz,
+ lut_a, lut_b, lut_c);
+
+ I_molecule = molecule_factor(intensities,
+ phases, flags, q,
+ ax,ay,az,bx,by,bz,cx,cy,cz,
+ m, sym);
+
+ I_lattice = pow(f_lattice, 2.0);
+ intensity = I_lattice * I_molecule;
+
+ image->data[idx] += intensity;
+
+ if ( fs_step + ss_step + kstep == 0 ) {
+ image->twotheta[idx] = twotheta;
+ }
+
+ }
+ }
+ }
+ }
+ }
+
+ image->data[idx] /= (SAMPLING*SAMPLING*BWSAMPLING
+ *DIVSAMPLING*DIVSAMPLING);
+
+
+ }
+ progress_bar(fs, image->width-1, "Calculating diffraction");
+ }
+
+ free(lut_a);
+ free(lut_b);
+ free(lut_c);
+}
diff --git a/libcrystfel/src/diffraction.h b/libcrystfel/src/diffraction.h
new file mode 100644
index 00000000..f71d3cce
--- /dev/null
+++ b/libcrystfel/src/diffraction.h
@@ -0,0 +1,34 @@
+/*
+ * diffraction.h
+ *
+ * Calculate diffraction patterns by Fourier methods
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef DIFFRACTION_H
+#define DIFFRACTION_H
+
+#include "image.h"
+#include "cell.h"
+#include "symmetry.h"
+
+typedef enum {
+ GRADIENT_MOSAIC,
+ GRADIENT_INTERPOLATE,
+ GRADIENT_PHASED
+} GradientMethod;
+
+extern void get_diffraction(struct image *image, int na, int nb, int nc,
+ const double *intensities, const double *phases,
+ const unsigned char *flags, UnitCell *cell,
+ GradientMethod m, const SymOpList *sym);
+
+#endif /* DIFFRACTION_H */
diff --git a/libcrystfel/src/filters.c b/libcrystfel/src/filters.c
new file mode 100644
index 00000000..c4e409df
--- /dev/null
+++ b/libcrystfel/src/filters.c
@@ -0,0 +1,130 @@
+/*
+ * filters.c
+ *
+ * Image filtering
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <gsl/gsl_statistics_int.h>
+#ifdef GSL_FUDGE
+#include <gsl/gsl_blas.h>
+#endif
+
+#include "image.h"
+
+
+static int compare_vals(const void *ap, const void *bp)
+{
+ const signed int a = *(signed int *)ap;
+ const signed int b = *(signed int *)bp;
+
+ if ( a > b ) return 1;
+ if ( a < b ) return -1;
+ return 0;
+}
+
+
+static void clean_panel(struct image *image, int sx, int sy)
+{
+ int x, y;
+ const int s = sizeof(signed int);
+
+ for ( x=0; x<512; x++ ) {
+
+ signed int vals[128];
+ double m;
+
+ for ( y=0; y<128; y++ ) {
+ vals[y] = image->data[(x+sx)+(y+sy)*image->width];
+ }
+
+ qsort(&vals[0], 128, s, compare_vals);
+
+ m = gsl_stats_int_median_from_sorted_data(vals, 1, 128);
+
+ for ( y=0; y<128; y++ ) {
+ image->data[(x+sx)+(y+sy)*image->width] -= m;
+ }
+
+ }
+}
+
+
+/* Pre-processing to make life easier */
+void filter_cm(struct image *image)
+{
+ int px, py;
+
+ if ( (image->width != 1024) || (image->height != 1024) ) return;
+
+ for ( px=0; px<2; px++ ) {
+ for ( py=0; py<8; py++ ) {
+
+ clean_panel(image, 512*px, 128*py);
+
+ }
+ }
+
+}
+
+
+void filter_noise(struct image *image, float *old)
+{
+ int x, y;
+
+ for ( x=0; x<image->width; x++ ) {
+ for ( y=0; y<image->height; y++ ) {
+
+ int dx, dy;
+ int val = image->data[x+image->width*y];
+
+ if ( old != NULL ) old[x+image->width*y] = val;
+
+ /* FIXME: This isn't really the right thing to do
+ * at the edges. */
+ if ( (x==0) || (x==image->width-1)
+ || (y==0) || (y==image->height-1) ) {
+ if ( val < 0 ) val = 0;
+ continue;
+ }
+
+ for ( dx=-1; dx<=+1; dx++ ) {
+ for ( dy=-1; dy<=+1; dy++ ) {
+
+ int val2;
+
+ val2 = image->data[(x+dx)+image->width*(y+dy)];
+
+ if ( val2 < 0 ) val = 0;
+
+ }
+ }
+
+ image->data[x+image->width*y] = val;
+
+ }
+ }
+}
+
+
+#ifdef GSL_FUDGE
+/* Force the linker to bring in CBLAS to make GSL happy */
+void filters_fudge_gslcblas()
+{
+ STATUS("%p\n", cblas_sgemm);
+}
+#endif
diff --git a/libcrystfel/src/filters.h b/libcrystfel/src/filters.h
new file mode 100644
index 00000000..6b35d7e8
--- /dev/null
+++ b/libcrystfel/src/filters.h
@@ -0,0 +1,25 @@
+/*
+ * peaks.h
+ *
+ * Image filtering
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef FILTERS_H
+#define FILTERS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+extern void filter_cm(struct image *image);
+extern void filter_noise(struct image *image, float *old);
+
+
+#endif /* FILTERS_H */
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
new file mode 100644
index 00000000..485abba3
--- /dev/null
+++ b/libcrystfel/src/geometry.c
@@ -0,0 +1,341 @@
+/*
+ * geometry.c
+ *
+ * Geometry of diffraction
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <assert.h>
+
+#include "utils.h"
+#include "cell.h"
+#include "image.h"
+#include "peaks.h"
+#include "beam-parameters.h"
+#include "reflist.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+
+
+static signed int locate_peak(double x, double y, double z, double k,
+ struct detector *det, double *xdap, double *ydap)
+{
+ int i;
+ signed int found = -1;
+ const double den = k + z;
+
+ *xdap = -1; *ydap = -1;
+
+ for ( i=0; i<det->n_panels; i++ ) {
+
+ double xd, yd;
+ double fs, ss, plx, ply;
+ struct panel *p;
+
+ p = &det->panels[i];
+
+ /* Coordinates of peak relative to central beam, in m */
+ xd = p->clen * x / den;
+ yd = p->clen * y / den;
+
+ /* Convert to pixels */
+ xd *= p->res;
+ yd *= p->res;
+
+ /* Convert to relative to the panel corner */
+ plx = xd - p->cnx;
+ ply = yd - p->cny;
+
+ fs = p->xfs*plx + p->yfs*ply;
+ ss = p->xss*plx + p->yss*ply;
+
+ fs += p->min_fs;
+ ss += p->min_ss;
+
+ /* Now, is this on this panel? */
+ if ( fs < p->min_fs ) continue;
+ if ( fs > p->max_fs ) continue;
+ if ( ss < p->min_ss ) continue;
+ if ( ss > p->max_ss ) continue;
+
+ /* If peak appears on multiple panels, reject it */
+ if ( found != -1 ) return -1;
+
+ /* Woohoo! */
+ found = i;
+ *xdap = fs;
+ *ydap = ss;
+
+ }
+
+ return found;
+}
+
+
+static double excitation_error(double xl, double yl, double zl,
+ double ds, double k, double divergence,
+ double tt)
+{
+ double al;
+ double r;
+ double delta;
+
+ al = M_PI_2 - asin(-zl/ds);
+
+ r = ( ds * sin(al) / sin(tt) ) - k;
+
+ delta = sqrt(2.0 * pow(ds, 2.0) * (1.0-cos(divergence)));
+ if ( divergence > 0.0 ) {
+ r += delta;
+ } else {
+ r -= delta;
+ }
+
+ return r;
+}
+
+
+static double partiality(double r1, double r2, double r)
+{
+ double q1, q2;
+ double p1, p2;
+
+ /* Calculate degrees of penetration */
+ q1 = (r1 + r)/(2.0*r);
+ q2 = (r2 + r)/(2.0*r);
+
+ /* Convert to partiality */
+ p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
+ p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
+
+ return p2 - p1;
+}
+
+
+static Reflection *check_reflection(struct image *image,
+ signed int h, signed int k, signed int l,
+ double asx, double asy, double asz,
+ double bsx, double bsy, double bsz,
+ double csx, double csy, double csz)
+{
+ const int output = 0;
+ double xl, yl, zl;
+ double ds, ds_sq;
+ double rlow, rhigh; /* "Excitation error" */
+ signed int p; /* Panel number */
+ double xda, yda; /* Position on detector */
+ int close, inside;
+ double part; /* Partiality */
+ int clamp_low = 0;
+ int clamp_high = 0;
+ double bandwidth = image->bw;
+ double divergence = image->div;
+ double lambda = image->lambda;
+ double klow, kcen, khigh; /* Wavenumber */
+ Reflection *refl;
+ double tt;
+
+ /* "low" gives the largest Ewald sphere,
+ * "high" gives the smallest Ewald sphere. */
+ klow = 1.0/(lambda - lambda*bandwidth/2.0);
+ kcen = 1.0/lambda;
+ khigh = 1.0/(lambda + lambda*bandwidth/2.0);
+
+ /* Get the coordinates of the reciprocal lattice point */
+ zl = h*asz + k*bsz + l*csz;
+ /* Throw out if it's "in front". A tiny bit "in front" is OK. */
+ if ( zl > image->profile_radius ) return NULL;
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+
+ tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+kcen);
+ if ( tt > deg2rad(90.0) ) return NULL;
+
+ ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
+ ds = sqrt(ds_sq);
+
+ /* Calculate excitation errors */
+ rlow = excitation_error(xl, yl, zl, ds, klow, -divergence/2.0, tt);
+ rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence/2.0, tt);
+
+ /* Is the reciprocal lattice point close to either extreme of
+ * the sphere, maybe just outside the "Ewald volume"? */
+ close = (fabs(rlow) < image->profile_radius)
+ || (fabs(rhigh) < image->profile_radius);
+
+ /* Is the reciprocal lattice point somewhere between the
+ * extremes of the sphere, i.e. inside the "Ewald volume"? */
+ inside = signbit(rlow) ^ signbit(rhigh);
+
+ /* Can't be both inside and close */
+ if ( inside ) close = 0;
+
+ /* Neither? Skip it. */
+ if ( !(close || inside) ) return NULL;
+
+ /* If the "lower" Ewald sphere is a long way away, use the
+ * position at which the Ewald sphere would just touch the
+ * reflection. */
+ if ( rlow < -image->profile_radius ) {
+ rlow = -image->profile_radius;
+ clamp_low = -1;
+ }
+ if ( rlow > +image->profile_radius ) {
+ rlow = +image->profile_radius;
+ clamp_low = +1;
+ }
+ /* Likewise the "higher" Ewald sphere */
+ if ( rhigh < -image->profile_radius ) {
+ rhigh = -image->profile_radius;
+ clamp_high = -1;
+ }
+ if ( rhigh > +image->profile_radius ) {
+ rhigh = +image->profile_radius;
+ clamp_high = +1;
+ }
+ assert(clamp_low <= clamp_high);
+ /* The six possible combinations of clamp_{low,high} (including
+ * zero) correspond to the six situations in Table 3 of Rossmann
+ * et al. (1979). */
+
+ /* Calculate partiality */
+ part = partiality(rlow, rhigh, image->profile_radius);
+
+ /* Locate peak on detector. */
+ p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
+ if ( p == -1 ) return NULL;
+
+ /* Add peak to list */
+ refl = reflection_new(h, k, l);
+ set_detector_pos(refl, 0.0, xda, yda);
+ set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+ set_symmetric_indices(refl, h, k, l);
+ set_redundancy(refl, 1);
+
+ if ( output ) {
+ printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
+ h, k, l, 0.0, xda, yda, part);
+ }
+
+ return refl;
+}
+
+
+RefList *find_intersections(struct image *image, UnitCell *cell)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ RefList *reflections;
+ int hmax, kmax, lmax;
+ double mres;
+ signed int h, k, l;
+
+ reflections = reflist_new();
+
+ /* Cell angle check from Foadi and Evans (2011) */
+ if ( !cell_is_sensible(cell) ) {
+ ERROR("Invalid unit cell parameters given to"
+ " find_intersections()\n");
+ cell_print(cell);
+ return NULL;
+ }
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ /* We add a horrific 20% fudge factor because bandwidth, divergence
+ * and so on mean reflections appear beyond the largest q */
+ mres = 1.2 * largest_q(image);
+
+ hmax = mres / modulus(asx, asy, asz);
+ kmax = mres / modulus(bsx, bsy, bsz);
+ lmax = mres / modulus(csx, csy, csz);
+
+ if ( (hmax >= 256) || (kmax >= 256) || (lmax >= 256) ) {
+ ERROR("Unit cell is stupidly large.\n");
+ cell_print(cell);
+ if ( hmax >= 256 ) hmax = 255;
+ if ( kmax >= 256 ) kmax = 255;
+ if ( lmax >= 256 ) lmax = 255;
+ }
+
+ for ( h=-hmax; h<=hmax; h++ ) {
+ for ( k=-kmax; k<=kmax; k++ ) {
+ for ( l=-lmax; l<=lmax; l++ ) {
+
+ Reflection *refl;
+
+ refl = check_reflection(image, h, k, l,
+ asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+
+ if ( refl != NULL ) {
+ refl = add_refl_to_list(refl, reflections);
+ }
+
+ }
+ }
+ }
+
+ return reflections;
+}
+
+
+/* Calculate partialities and apply them to the image's reflections */
+void update_partialities(struct image *image)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ RefList *predicted;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz, &csx, &csy, &csz);
+
+ /* Scratch list to give check_reflection() something to add to */
+ predicted = reflist_new();
+
+ for ( refl = first_refl(image->reflections, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ Reflection *vals;
+ double r1, r2, p, x, y;
+ signed int h, k, l;
+ int clamp1, clamp2;
+
+ get_symmetric_indices(refl, &h, &k, &l);
+
+ vals = check_reflection(image, h, k, l,
+ asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+
+ if ( vals == NULL ) {
+ set_redundancy(refl, 0);
+ continue;
+ }
+ set_redundancy(refl, 1);
+
+ /* Transfer partiality stuff */
+ get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2);
+ set_partial(refl, r1, r2, p, clamp1, clamp2);
+
+ /* Transfer detector location */
+ get_detector_pos(vals, &x, &y);
+ set_detector_pos(refl, 0.0, x, y);
+ }
+
+ reflist_free(predicted);
+}
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
new file mode 100644
index 00000000..ddf04b80
--- /dev/null
+++ b/libcrystfel/src/geometry.h
@@ -0,0 +1,26 @@
+/*
+ * geometry.h
+ *
+ * Geometry of diffraction
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef GEOMETRY_H
+#define GEOMETRY_H
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "reflist.h"
+
+extern RefList *find_intersections(struct image *image, UnitCell *cell);
+
+extern void update_partialities(struct image *image);
+
+#endif /* GEOMETRY_H */
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
new file mode 100644
index 00000000..ad524c61
--- /dev/null
+++ b/libcrystfel/src/peaks.c
@@ -0,0 +1,593 @@
+/*
+ * peaks.c
+ *
+ * Peak search and other image analysis
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ * 2011 Andrew Martin
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <gsl/gsl_statistics_int.h>
+#include <pthread.h>
+#include <fenv.h>
+
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "detector.h"
+#include "filters.h"
+#include "diffraction.h"
+#include "reflist-utils.h"
+#include "beam-parameters.h"
+
+
+/* How close a peak must be to an indexed position to be considered "close"
+ * for the purposes of double hit detection and sanity checking. */
+#define PEAK_CLOSE (30.0)
+
+/* How close a peak must be to an indexed position to be considered "close"
+ * for the purposes of integration. */
+#define PEAK_REALLY_CLOSE (10.0)
+
+/* Degree of polarisation of X-ray beam */
+#define POL (1.0)
+
+static int cull_peaks_in_panel(struct image *image, struct panel *p)
+{
+ int i, n;
+ int nelim = 0;
+
+ n = image_feature_count(image->features);
+
+ for ( i=0; i<n; i++ ) {
+
+ struct imagefeature *f;
+ int j, ncol;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ if ( f->fs < p->min_fs ) continue;
+ if ( f->fs > p->max_fs ) continue;
+ if ( f->ss < p->min_ss ) continue;
+ if ( f->ss > p->max_ss ) continue;
+
+ /* How many peaks are in the same column? */
+ ncol = 0;
+ for ( j=0; j<n; j++ ) {
+
+ struct imagefeature *g;
+
+ if ( i==j ) continue;
+
+ g = image_get_feature(image->features, j);
+ if ( g == NULL ) continue;
+
+ if ( p->badrow == 'f' ) {
+ if ( fabs(f->ss - g->ss) < 2.0 ) ncol++;
+ } else if ( p->badrow == 's' ) {
+ if ( fabs(f->fs - g->fs) < 2.0 ) ncol++;
+ } /* else do nothing */
+
+ }
+
+ /* More than three? */
+ if ( ncol <= 3 ) continue;
+
+ /* Yes? Delete them all... */
+ nelim = 0;
+ for ( j=0; j<n; j++ ) {
+ struct imagefeature *g;
+ g = image_get_feature(image->features, j);
+ if ( g == NULL ) continue;
+ if ( p->badrow == 'f' ) {
+ if ( fabs(f->ss - g->ss) < 2.0 ) {
+ image_remove_feature(image->features,
+ j);
+ nelim++;
+ }
+ } else if ( p->badrow == 's' ) {
+ if ( fabs(f->fs - g->ss) < 2.0 ) {
+ image_remove_feature(image->features,
+ j);
+ nelim++;
+ }
+ } else {
+ ERROR("Invalid badrow direction.\n");
+ abort();
+ }
+
+ }
+
+ }
+
+ return nelim;
+}
+
+
+/* Post-processing of the peak list to remove noise */
+static int cull_peaks(struct image *image)
+{
+ int nelim = 0;
+ struct panel *p;
+ int i;
+
+ for ( i=0; i<image->det->n_panels; i++ ) {
+ p = &image->det->panels[i];
+ if ( p->badrow != '-' ) {
+ nelim += cull_peaks_in_panel(image, p);
+ }
+ }
+
+ return nelim;
+}
+
+
+/* Returns non-zero if peak has been vetoed.
+ * i.e. don't use result if return value is not zero. */
+int integrate_peak(struct image *image, int cfs, int css,
+ double *pfs, double *pss, double *intensity,
+ double *pbg, double *pmax, double *sigma,
+ int do_polar, int centroid, int bgsub)
+{
+ signed int fs, ss;
+ double lim, out_lim, mid_lim;
+ double lim_sq, out_lim_sq, mid_lim_sq;
+ double total = 0.0;
+ double fsct = 0.0;
+ double ssct = 0.0;
+ double noise = 0.0;
+ int noise_counts = 0;
+ double max = 0.0;
+ struct panel *p = NULL;
+ int pixel_counts = 0;
+ double noise_mean = 0.0;
+ double noise_meansq = 0.0;
+ struct beam_params *beam;
+ double aduph;
+
+ beam = image->beam;
+ if ( beam != NULL ) {
+ aduph = image->beam->adu_per_photon;
+ } else {
+ aduph = 1.0;
+ }
+
+ p = find_panel(image->det, cfs, css);
+ if ( p == NULL ) return 1;
+ if ( p->no_index ) return 1;
+
+ lim = p->integr_radius;
+ mid_lim = 3.0 + lim;
+ out_lim = 6.0 + lim;
+ lim_sq = pow(lim, 2.0);
+ mid_lim_sq = pow(mid_lim, 2.0);
+ out_lim_sq = pow(out_lim, 2.0);
+
+ for ( fs=-out_lim; fs<+out_lim; fs++ ) {
+ for ( ss=-out_lim; ss<+out_lim; ss++ ) {
+
+ double val;
+ double tt = 0.0;
+ double phi, pa, pb, pol;
+ uint16_t flags;
+ struct panel *p2;
+ int idx;
+
+ /* Outer mask radius */
+ if ( fs*fs + ss*ss > out_lim_sq ) continue;
+
+ if ( ((fs+cfs)>=image->width) || ((fs+cfs)<0) ) continue;
+ if ( ((ss+css)>=image->height) || ((ss+css)<0) ) continue;
+
+ /* Strayed off one panel? */
+ p2 = find_panel(image->det, fs+cfs, ss+css);
+ if ( p2 != p ) return 1;
+
+ idx = fs+cfs+image->width*(ss+css);
+
+ /* Veto this peak if we tried to integrate in a bad region */
+ if ( image->flags != NULL ) {
+
+ flags = image->flags[idx];
+
+ /* It must have all the "good" bits to be valid */
+ if ( !((flags & image->det->mask_good)
+ == image->det->mask_good) ) return 1;
+
+ /* If it has any of the "bad" bits, reject */
+ if ( flags & image->det->mask_bad ) return 1;
+
+ }
+
+ val = image->data[idx];
+
+ if ( do_polar ) {
+
+ tt = get_tt(image, fs+cfs, ss+css);
+
+ phi = atan2(ss+css, fs+cfs);
+ pa = pow(sin(phi)*sin(tt), 2.0);
+ pb = pow(cos(tt), 2.0);
+ pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
+
+ val /= pol;
+
+ }
+
+ if ( val > max ) max = val;
+
+ /* If outside inner mask, estimate noise from this region */
+ if ( fs*fs + ss*ss > mid_lim_sq ) {
+
+ /* Noise
+ * noise and noise_meansq are both in photons (^2) */
+ noise += val / image->beam->adu_per_photon;
+ noise_counts++;
+ noise_meansq += pow(val, 2.0);
+
+ } else if ( fs*fs + ss*ss < lim_sq ) {
+
+ /* Peak */
+ pixel_counts++;
+ total += val;
+ fsct += val*(cfs+fs);
+ ssct += val*(css+ss);
+
+ }
+
+ }
+ }
+
+ noise_mean = noise / noise_counts; /* photons */
+
+ /* The centroid is excitingly undefined if there is no intensity */
+ centroid = 0;
+ if ( centroid && (total != 0) ) {
+ *pfs = ((double)fsct / total) + 0.5;
+ *pss = ((double)ssct / total) + 0.5;
+ } else {
+ *pfs = (double)cfs + 0.5;
+ *pss = (double)css + 0.5;
+ }
+ if ( bgsub ) {
+ *intensity = total - aduph * pixel_counts*noise_mean; /* ADU */
+ } else {
+ *intensity = total; /* ADU */
+ }
+
+ if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
+
+ if ( sigma != NULL ) {
+
+ /* First term is standard deviation of background per pixel
+ * sqrt(pixel_counts) - increase of error for integrated value
+ * sqrt(2) - increase of error for background subtraction */
+ *sigma = sqrt(noise_meansq/noise_counts-(noise_mean*noise_mean))
+ * sqrt(2.0*pixel_counts) * aduph;
+
+ }
+
+ if ( pbg != NULL ) {
+ *pbg = aduph * (noise / noise_counts);
+ }
+ if ( pmax != NULL ) {
+ *pmax = max;
+ }
+
+ return 0;
+}
+
+
+static void search_peaks_in_panel(struct image *image, float threshold,
+ float min_gradient, float min_snr,
+ struct panel *p)
+{
+ int fs, ss, stride;
+ float *data;
+ double d;
+ int idx;
+ double f_fs = 0.0;
+ double f_ss = 0.0;
+ double intensity = 0.0;
+ double sigma = 0.0;
+ double pbg = 0.0;
+ double pmax = 0.0;
+ int nrej_dis = 0;
+ int nrej_pro = 0;
+ int nrej_fra = 0;
+ int nrej_bad = 0;
+ int nrej_snr = 0;
+ int nacc = 0;
+ int ncull;
+ const int pws = p->peak_sep/2;
+
+ data = image->data;
+ stride = image->width;
+
+ for ( fs = p->min_fs+1; fs <= p->max_fs-1; fs++ ) {
+ for ( ss = p->min_ss+1; ss <= p->max_ss-1; ss++ ) {
+
+ double dx1, dx2, dy1, dy2;
+ double dxs, dys;
+ double grad;
+ int mask_fs, mask_ss;
+ int s_fs, s_ss;
+ double max;
+ unsigned int did_something;
+ int r;
+
+ /* Overall threshold */
+ if ( data[fs+stride*ss] < threshold ) continue;
+
+ /* Get gradients */
+ dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss];
+ dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss];
+ dy1 = data[fs+stride*ss] - data[(fs+1)+stride*(ss+1)];
+ dy2 = data[fs+stride*(ss-1)] - data[fs+stride*ss];
+
+ /* Average gradient measurements from both sides */
+ dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
+ dys = ((dy1*dy1) + (dy2*dy2)) / 2;
+
+ /* Calculate overall gradient */
+ grad = dxs + dys;
+
+ if ( grad < min_gradient ) continue;
+
+ mask_fs = fs;
+ mask_ss = ss;
+
+ do {
+
+ max = data[mask_fs+stride*mask_ss];
+ did_something = 0;
+
+ for ( s_ss=biggest(mask_ss-pws/2,
+ p->min_ss);
+ s_ss<=smallest(mask_ss+pws/2,
+ p->max_ss);
+ s_ss++ ) {
+ for ( s_fs=biggest(mask_fs-pws/2,
+ p->min_fs);
+ s_fs<=smallest(mask_fs+pws/2,
+ p->max_fs);
+ s_fs++ ) {
+
+ if ( data[s_fs+stride*s_ss] > max ) {
+ max = data[s_fs+stride*s_ss];
+ mask_fs = s_fs;
+ mask_ss = s_ss;
+ did_something = 1;
+ }
+
+ }
+ }
+
+ /* Abort if drifted too far from the foot point */
+ if ( distance(mask_fs, mask_ss, fs, ss) >
+ p->peak_sep/2.0 )
+ {
+ break;
+ }
+
+ } while ( did_something );
+
+ /* Too far from foot point? */
+ if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep/2.0 ) {
+ nrej_dis++;
+ continue;
+ }
+
+ /* Should be enforced by bounds used above. Muppet check. */
+ assert(mask_fs <= p->max_fs);
+ assert(mask_ss <= p->max_ss);
+ assert(mask_fs >= p->min_fs);
+ assert(mask_ss >= p->min_ss);
+
+ /* Centroid peak and get better coordinates.
+ * Don't bother doing polarisation/SA correction, because the
+ * intensity of this peak is only an estimate at this stage. */
+ r = integrate_peak(image, mask_fs, mask_ss,
+ &f_fs, &f_ss, &intensity,
+ &pbg, &pmax, &sigma, 0, 1, 1);
+
+ if ( r ) {
+ /* Bad region - don't detect peak */
+ nrej_bad++;
+ continue;
+ }
+
+ /* It is possible for the centroid to fall outside the image */
+ if ( (f_fs < p->min_fs) || (f_fs > p->max_fs)
+ || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) {
+ nrej_fra++;
+ continue;
+ }
+
+ if (intensity/sigma < min_snr) {
+ nrej_snr++;
+ continue;
+ }
+
+ /* Check for a nearby feature */
+ image_feature_closest(image->features, f_fs, f_ss, &d, &idx);
+ if ( d < p->peak_sep/2.0 ) {
+ nrej_pro++;
+ continue;
+ }
+
+ /* Add using "better" coordinates */
+ image_add_feature(image->features, f_fs, f_ss, image, intensity,
+ NULL);
+ nacc++;
+
+ }
+ }
+
+ if ( image->det != NULL ) {
+ ncull = cull_peaks(image);
+ nacc -= ncull;
+ } else {
+ STATUS("Not culling peaks because I don't have a "
+ "detector geometry file.\n");
+ ncull = 0;
+ }
+
+// STATUS("%i accepted, %i box, %i proximity, %i outside panel, "
+// "%i in bad regions, %i with SNR < %g, %i badrow culled.\n",
+// nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
+// nrej_snr, min_snr, ncull);
+}
+
+
+void search_peaks(struct image *image, float threshold, float min_gradient,
+ float min_snr)
+{
+ int i;
+
+ if ( image->features != NULL ) {
+ image_feature_list_free(image->features);
+ }
+ image->features = image_feature_list_new();
+
+ for ( i=0; i<image->det->n_panels; i++ ) {
+
+ struct panel *p = &image->det->panels[i];
+
+ if ( p->no_index ) continue;
+ search_peaks_in_panel(image, threshold, min_gradient, min_snr, p);
+
+ }
+}
+
+
+int peak_sanity_check(struct image *image)
+{
+
+ int i;
+ int n_feat = 0;
+ int n_sane = 0;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ double min_dist = 0.25;
+
+ /* Round towards nearest */
+ fesetround(1);
+
+ /* Cell basis vectors for this image */
+ cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
+
+ /* Loop over peaks, checking proximity to nearest reflection */
+ for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+ struct imagefeature *f;
+ struct rvec q;
+ double h,k,l,hd,kd,ld;
+
+ /* Assume all image "features" are genuine peaks */
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+ n_feat++;
+
+ /* Reciprocal space position of found peak */
+ q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
+
+ /* Decimal and fractional Miller indices of nearest
+ * reciprocal lattice point */
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * cz;
+ h = lrint(hd);
+ k = lrint(kd);
+ l = lrint(ld);
+
+ /* Check distance */
+ if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist)
+ && (fabs(l - ld) < min_dist) )
+ {
+ n_sane++;
+ continue;
+ }
+
+ }
+
+ /* return 0 means fail test, return 1 means pass test */
+ // printf("%d out of %d peaks are \"sane\"\n",n_sane,n_feat);
+ if ( (float)n_sane / (float)n_feat < 0.5 ) return 0;
+
+ return 1;
+}
+
+
+/* Integrate the list of predicted reflections in "image" */
+void integrate_reflections(struct image *image, int polar, int use_closer,
+ int bgsub)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+
+ for ( refl = first_refl(image->reflections, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ double fs, ss, intensity;
+ double d;
+ int idx;
+ double bg, max;
+ double sigma;
+ double pfs, pss;
+ int r;
+
+ get_detector_pos(refl, &pfs, &pss);
+
+ /* Is there a really close feature which was detected? */
+ if ( use_closer ) {
+
+ struct imagefeature *f;
+
+ if ( image->features != NULL ) {
+ f = image_feature_closest(image->features,
+ pfs, pss, &d, &idx);
+ } else {
+ f = NULL;
+ }
+ if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) {
+
+ pfs = f->fs;
+ pss = f->ss;
+
+ }
+ }
+
+ r = integrate_peak(image, pfs, pss, &fs, &ss,
+ &intensity, &bg, &max, &sigma, polar, 0,
+ bgsub);
+
+ /* Record intensity and set redundancy to 1 on success */
+ if ( r == 0 ) {
+ set_int(refl, intensity);
+ set_esd_intensity(refl, sigma);
+ set_redundancy(refl, 1);
+ } else {
+ set_redundancy(refl, 0);
+ }
+
+ }
+}
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
new file mode 100644
index 00000000..9d475ea9
--- /dev/null
+++ b/libcrystfel/src/peaks.h
@@ -0,0 +1,38 @@
+/*
+ * peaks.h
+ *
+ * Peak search and other image analysis
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef PEAKS_H
+#define PEAKS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <pthread.h>
+
+#include "reflist.h"
+
+extern void search_peaks(struct image *image, float threshold,
+ float min_gradient, float min_snr);
+
+extern void integrate_reflections(struct image *image,
+ int polar, int use_closer, int bgsub);
+
+extern int peak_sanity_check(struct image * image);
+
+/* Exported so it can be poked by integration_check */
+extern int integrate_peak(struct image *image, int cfs, int css,
+ double *pfs, double *pss, double *intensity,
+ double *pbg, double *pmax, double *sigma,
+ int do_polar, int centroid, int bgsub);
+
+#endif /* PEAKS_H */
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
new file mode 100644
index 00000000..b64e9979
--- /dev/null
+++ b/libcrystfel/src/reflist-utils.c
@@ -0,0 +1,502 @@
+/*
+ * reflist-utils.c
+ *
+ * Utilities to complement the core reflist.c
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdio.h>
+#include <assert.h>
+
+
+#include "reflist.h"
+#include "cell.h"
+#include "utils.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+
+
+/**
+ * SECTION:reflist-utils
+ * @short_description: Reflection list utilities
+ * @title: RefList utilities
+ * @section_id:
+ * @see_also:
+ * @include: "reflist-utils.h"
+ * @Image:
+ *
+ * There are some utility functions associated with the core %RefList.
+ **/
+
+
+double *intensities_from_list(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double *out = new_list_intensity();
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ double intensity = get_intensity(refl);
+
+ get_indices(refl, &h, &k, &l);
+
+ set_intensity(out, h, k, l, intensity);
+
+ }
+
+ return out;
+}
+
+
+double *phases_from_list(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double *out = new_list_phase();
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ double phase = get_phase(refl, NULL);
+
+ get_indices(refl, &h, &k, &l);
+
+ set_phase(out, h, k, l, phase);
+
+ }
+
+ return out;
+
+}
+
+
+unsigned char *flags_from_list(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ unsigned char *out = new_list_flag();
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+
+ set_flag(out, h, k, l, 1);
+
+ }
+
+ return out;
+
+}
+
+
+int check_list_symmetry(RefList *list, const SymOpList *sym)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ SymOpMask *mask;
+
+ mask = new_symopmask(sym);
+ if ( mask == NULL ) {
+ ERROR("Couldn't create mask for list symmetry check.\n");
+ return 1;
+ }
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ int j;
+ int found = 0;
+ signed int h, k, l;
+ int n;
+
+ get_indices(refl, &h, &k, &l);
+
+ special_position(sym, mask, h, k, l);
+ n = num_equivs(sym, mask);
+
+ for ( j=0; j<n; j++ ) {
+
+ signed int he, ke, le;
+ Reflection *f;
+
+ get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
+
+ f = find_refl(list, he, ke, le);
+ if ( f != NULL ) found++;
+
+ }
+
+ assert(found != 0); /* That'd just be silly */
+ if ( found > 1 ) {
+
+ STATUS("Found %i %i %i: %i times:\n", h, k, l, found);
+
+ for ( j=0; j<n; j++ ) {
+
+ signed int he, ke, le;
+ Reflection *f;
+
+ get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
+
+ f = find_refl(list, he, ke, le);
+ if ( f != NULL ) {
+ STATUS("%3i %3i %3i\n", he, ke, le);
+ }
+
+ }
+ free_symopmask(mask);
+
+ return 1; /* Symmetry is wrong! */
+ }
+
+ }
+
+ free_symopmask(mask);
+
+ return 0;
+}
+
+
+int find_equiv_in_list(RefList *list, signed int h, signed int k,
+ signed int l, const SymOpList *sym, signed int *hu,
+ signed int *ku, signed int *lu)
+{
+ int i;
+ int found = 0;
+
+ for ( i=0; i<num_equivs(sym, NULL); i++ ) {
+
+ signed int he, ke, le;
+ Reflection *f;
+ get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
+ f = find_refl(list, he, ke, le);
+
+ /* There must only be one equivalent. If there are more, it
+ * indicates that the user lied about the input symmetry.
+ * This situation should have been checked for earlier by
+ * calling check_symmetry() with 'items' and 'mero'. */
+
+ if ( (f != NULL) && !found ) {
+ *hu = he; *ku = ke; *lu = le;
+ return 1;
+ }
+
+ }
+
+ return 0;
+}
+
+
+/**
+ * write_reflections_to_file:
+ * @fh: File handle to write to
+ * @list: The reflection list to write
+ * @cell: Unit cell to use for generating 1/d values, or NULL.
+ *
+ * This function writes the contents of @list to @fh, using @cell to generate
+ * 1/d values to ease later processing. If @cell is NULL, 1/d values will not
+ * be included ('-' will be written in their place).
+ *
+ * Reflections which have a redundancy of zero will not be written.
+ *
+ * The resulting list can be read back with read_reflections_from_file().
+ **/
+void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+
+ fprintf(fh, " h k l I phase sigma(I) "
+ " 1/d(nm^-1) counts fs/px ss/px\n");
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ double intensity, esd_i, s, ph;
+ int red;
+ double fs, ss;
+ char res[16];
+ char phs[16];
+ int have_phase;
+
+ get_indices(refl, &h, &k, &l);
+ get_detector_pos(refl, &fs, &ss);
+ intensity = get_intensity(refl);
+ esd_i = get_esd_intensity(refl);
+ red = get_redundancy(refl);
+ ph = get_phase(refl, &have_phase);
+
+ /* Reflections with redundancy = 0 are not written */
+ if ( red == 0 ) continue;
+
+ if ( cell != NULL ) {
+ s = 2.0 * resolution(cell, h, k, l);
+ snprintf(res, 16, "%10.2f", s/1e9);
+ } else {
+ strcpy(res, " -");
+ }
+
+ if ( have_phase ) {
+ snprintf(phs, 16, "%8.2f", rad2deg(ph));
+ } else {
+ strncpy(phs, " -", 15);
+ }
+
+ fprintf(fh,
+ "%3i %3i %3i %10.2f %s %10.2f %s %7i %6.1f %6.1f\n",
+ h, k, l, intensity, phs, esd_i, res, red,
+ fs, ss);
+
+ }
+}
+
+
+/**
+ * write_reflist:
+ * @filename: Filename
+ * @list: The reflection list to write
+ * @cell: Unit cell to use for generating 1/d values, or NULL.
+ *
+ * This function writes the contents of @list to @file, using @cell to generate
+ * 1/d values to ease later processing. If @cell is NULL, 1/d values will not
+ * be included ('-' will be written in their place).
+ *
+ * Reflections which have a redundancy of zero will not be written.
+ *
+ * The resulting list can be read back with read_reflections_from_file() or
+ * read_reflections().
+ *
+ * This is a convenience function which simply opens @filename and then calls
+ * write_reflections_to_file.
+ *
+ * Returns: zero on success, non-zero on failure.
+ **/
+int write_reflist(const char *filename, RefList *list, UnitCell *cell)
+{
+ FILE *fh;
+
+ if ( filename == NULL ) {
+ fh = stdout;
+ } else {
+ fh = fopen(filename, "w");
+ }
+
+ if ( fh == NULL ) {
+ ERROR("Couldn't open output file '%s'.\n", filename);
+ return 1;
+ }
+
+ write_reflections_to_file(fh, list, cell);
+ fprintf(fh, REFLECTION_END_MARKER"\n");
+
+ fclose(fh);
+
+ return 0;
+}
+
+
+RefList *read_reflections_from_file(FILE *fh)
+{
+ char *rval = NULL;
+ int first = 1;
+ RefList *out;
+
+ out = reflist_new();
+
+ do {
+
+ char line[1024];
+ signed int h, k, l;
+ float intensity, sigma, fs, ss;
+ char phs[1024];
+ char ress[1024];
+ int cts;
+ int r;
+ Reflection *refl;
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ chomp(line);
+
+ if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out;
+
+ r = sscanf(line, "%i %i %i %f %s %f %s %i %f %f",
+ &h, &k, &l, &intensity, phs, &sigma, ress, &cts,
+ &fs, &ss);
+ if ( (r != 10) && (!first) ) {
+ reflist_free(out);
+ return NULL;
+ }
+
+ first = 0;
+ if ( r == 10 ) {
+
+ double ph;
+ char *v;
+
+ refl = add_refl(out, h, k, l);
+ set_int(refl, intensity);
+ set_detector_pos(refl, 0.0, fs, ss);
+ set_esd_intensity(refl, sigma);
+ set_redundancy(refl, cts);
+
+ ph = strtod(phs, &v);
+ if ( v != NULL ) set_ph(refl, deg2rad(ph));
+
+ /* The 1/d value is actually ignored. */
+
+ }
+
+ } while ( rval != NULL );
+
+ /* Got read error of some kind before finding PEAK_LIST_END_MARKER */
+ return NULL;
+}
+
+
+RefList *read_reflections(const char *filename)
+{
+ FILE *fh;
+ RefList *out;
+
+ if ( filename == NULL ) {
+ fh = stdout;
+ } else {
+ fh = fopen(filename, "r");
+ }
+
+ if ( fh == NULL ) {
+ ERROR("Couldn't open input file '%s'.\n", filename);
+ return NULL;
+ }
+
+ out = read_reflections_from_file(fh);
+
+ fclose(fh);
+
+ return out;
+}
+
+
+/**
+ * asymmetric_indices:
+ * @in: A %RefList
+ * @sym: A %SymOpList
+ *
+ * This function creates a newly allocated copy of @in, but indexed using the
+ * asymmetric indices according to @sym instead of the original indices. The
+ * original indices are stored and can be retrieved using
+ * get_symmetric_indices() if required.
+ *
+ * Returns: the new %RefList, or NULL on failure.
+ **/
+RefList *asymmetric_indices(RefList *in, const SymOpList *sym)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ RefList *new;
+
+ new = reflist_new();
+ if ( new == NULL ) return NULL;
+
+ for ( refl = first_refl(in, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ signed int ha, ka, la;
+ Reflection *cr;
+
+ get_indices(refl, &h, &k, &l);
+
+ get_asymm(sym, h, k, l, &ha, &ka, &la);
+
+ cr = add_refl(new, ha, ka, la);
+ assert(cr != NULL);
+
+ copy_data(cr, refl);
+ set_symmetric_indices(cr, h, k, l);
+
+ }
+
+ return new;
+}
+
+
+/**
+ * resolution_limits:
+ * @list: A %RefList
+ * @cell: A %UnitCell
+ * @rmin: Place to store the minimum 1/d value
+ * @rmax: Place to store the maximum 1/d value
+ *
+ * This function calculates the minimum and maximum values of 1/d, where
+ * 2dsin(theta) = wavelength. The answers are in m^-1.
+ **/
+void resolution_limits(RefList *list, UnitCell *cell,
+ double *rmin, double *rmax)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+
+ *rmin = INFINITY;
+ *rmax = 0.0;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double r;
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+ r = 2.0 * resolution(cell, h, k, l);
+
+ if ( r > *rmax ) *rmax = r;
+ if ( r < *rmin ) *rmin = r;
+ }
+}
+
+
+/**
+ * max_intensity:
+ * @list: A %RefList
+ *
+ * Returns: The maximum intensity in @list.
+ **/
+double max_intensity(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double max;
+
+ max = -INFINITY;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double val = get_intensity(refl);
+ if ( val > max ) max = val;
+ }
+
+ return max;
+}
diff --git a/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h
new file mode 100644
index 00000000..d14e5f8e
--- /dev/null
+++ b/libcrystfel/src/reflist-utils.h
@@ -0,0 +1,52 @@
+/*
+ * reflist-utils.h
+ *
+ * Utilities to complement the core reflist.c
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef REFLIST_UTILS_H
+#define REFLIST_UTILS_H
+
+
+#include "reflist.h"
+#include "cell.h"
+#include "symmetry.h"
+
+
+#define REFLECTION_END_MARKER "End of reflections"
+
+
+extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell);
+
+extern int write_reflist(const char *filename, RefList *list, UnitCell *cell);
+
+extern RefList *read_reflections_from_file(FILE *fh);
+
+extern RefList *read_reflections(const char *filename);
+
+extern double *intensities_from_list(RefList *list);
+extern double *phases_from_list(RefList *list);
+extern unsigned char *flags_from_list(RefList *list);
+
+extern int check_list_symmetry(RefList *list, const SymOpList *sym);
+extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
+ signed int l, const SymOpList *sym, signed int *hu,
+ signed int *ku, signed int *lu);
+
+extern RefList *asymmetric_indices(RefList *in, const SymOpList *sym);
+
+extern void resolution_limits(RefList *list, UnitCell *cell,
+ double *rmin, double *rmax);
+
+extern double max_intensity(RefList *list);
+
+#endif /* REFLIST_UTILS_H */
diff --git a/libcrystfel/src/statistics.c b/libcrystfel/src/statistics.c
new file mode 100644
index 00000000..7990672d
--- /dev/null
+++ b/libcrystfel/src/statistics.c
@@ -0,0 +1,668 @@
+/*
+ * statistics.c
+ *
+ * Structure-factor statistics
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_min.h>
+#include <gsl/gsl_statistics.h>
+
+#include "statistics.h"
+#include "utils.h"
+
+/**
+ * SECTION:statistics
+ * @short_description: Intensity statistics and R-factors
+ * @title: Statistics
+ * @section_id:
+ * @see_also:
+ * @include: "statistics.h"
+ * @Image:
+ *
+ * These functions are for calculating various figures of merit.
+ */
+
+
+struct r_params {
+ RefList *list1;
+ RefList *list2;
+ int fom; /* Which FoM to use (see the enum just below) */
+};
+
+enum {
+ R_1_ZERO,
+ R_1_IGNORE,
+ R_2,
+ R_1_I,
+ R_DIFF_ZERO,
+ R_DIFF_IGNORE,
+ R_DIFF_INTENSITY,
+};
+
+
+/* Return the least squares optimal scaling factor when comparing intensities.
+ * list1,list2 are the two intensity lists to compare.
+ */
+double stat_scale_intensity(RefList *list1, RefList *list2)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ top += i1 * i2;
+ bot += i2 * i2;
+
+ }
+
+ return top/bot;
+}
+
+
+/* Return the least squares optimal scaling factor when comparing the square
+ * roots of the intensities (i.e. one approximation to the structure factor
+ * moduli).
+ * list1,list2 are the two intensity lists to compare (they contain intensities,
+ * not square rooted intensities).
+ */
+static double stat_scale_sqrti(RefList *list1, RefList *list2)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
+
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
+
+ top += f1 * f2;
+ bot += f2 * f2;
+
+ }
+
+ return top/bot;
+}
+
+
+static double internal_r1_ignorenegs(RefList *list1, RefList *list2,
+ double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
+
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
+ f2 *= scale;
+
+ top += fabs(f1 - f2);
+ bot += f1;
+
+ }
+
+ return top/bot;
+}
+
+
+static double internal_r1_negstozero(RefList *list1, RefList *list2,
+ double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
+
+ f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
+ f2 *= scale;
+
+ top += fabs(f1 - f2);
+ bot += f1;
+
+ }
+
+ return top/bot;
+}
+
+
+static double internal_r2(RefList *list1, RefList *list2, double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ i2 *= scale;
+
+ top += pow(i1 - i2, 2.0);
+ bot += pow(i1, 2.0);
+
+ }
+
+ return sqrt(top/bot);
+}
+
+
+static double internal_r_i(RefList *list1, RefList *list2, double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+ i2 *= scale;
+
+ top += fabs(i1-i2);
+ bot += fabs(i1);
+
+ }
+
+ return top/bot;
+}
+
+
+static double internal_rdiff_intensity(RefList *list1, RefList *list2,
+ double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+ i2 *= scale;
+
+ top += fabs(i1 - i2);
+ bot += i1 + i2;
+
+ }
+
+ return 2.0*top/bot;
+}
+
+
+static double internal_rdiff_negstozero(RefList *list1, RefList *list2,
+ double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
+
+ f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
+ f2 *= scale;
+
+ top += fabs(f1 - f2);
+ bot += f1 + f2;
+
+ }
+
+ return 2.0*top/bot;
+}
+
+
+static double internal_rdiff_ignorenegs(RefList *list1, RefList *list2,
+ double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
+
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
+ f2 *= scale;
+
+ top += fabs(f1 - f2);
+ bot += f1 + f2;
+
+ }
+
+ return 2.0*top/bot;
+}
+
+
+static double calc_r(double scale, void *params)
+{
+ struct r_params *rp = params;
+
+ switch ( rp->fom ) {
+ case R_1_ZERO :
+ return internal_r1_negstozero(rp->list1, rp->list2, scale);
+ case R_1_IGNORE :
+ return internal_r1_ignorenegs(rp->list1, rp->list2, scale);
+ case R_2 :
+ return internal_r2(rp->list1, rp->list2, scale);
+
+ case R_1_I :
+ return internal_r_i(rp->list1, rp->list2, scale);
+
+ case R_DIFF_ZERO :
+ return internal_rdiff_negstozero(rp->list1, rp->list2,scale);
+ case R_DIFF_IGNORE :
+ return internal_rdiff_ignorenegs(rp->list1, rp->list2, scale);
+ case R_DIFF_INTENSITY :
+ return internal_rdiff_intensity(rp->list1, rp->list2, scale);
+ }
+
+ ERROR("No such FoM!\n");
+ abort();
+}
+
+
+static double r_minimised(RefList *list1, RefList *list2, double *scalep, int fom,
+ int u)
+{
+ gsl_function F;
+ gsl_min_fminimizer *s;
+ int status;
+ double scale = 1.0;
+ struct r_params rp;
+ int iter = 0;
+
+ rp.list1 = list1;
+ rp.list2 = list2;
+ rp.fom = fom;
+
+ if ( u ) {
+
+ scale = 1.0;
+
+ } else {
+
+ F.function = &calc_r;
+ F.params = &rp;
+
+ s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);
+
+ /* Initial guess */
+ switch ( fom ) {
+ case R_1_ZERO :
+ case R_1_IGNORE :
+ case R_DIFF_ZERO :
+ case R_DIFF_IGNORE :
+ scale = stat_scale_sqrti(list1, list2);
+ break;
+ case R_2 :
+ case R_1_I :
+ case R_DIFF_INTENSITY :
+ scale = stat_scale_intensity(list1, list2);
+ break;
+ }
+ //STATUS("Initial scale factor estimate: %5.2e\n", scale);
+
+ /* Probably within an order of magnitude either side */
+ gsl_min_fminimizer_set(s, &F, scale, scale/10.0, scale*10.0);
+
+ do {
+
+ double lo, up;
+
+ /* Iterate */
+ if ( gsl_min_fminimizer_iterate(s) ) {
+ ERROR("Failed to find scale factor.\n");
+ return NAN;
+ }
+
+ /* Get the current estimate */
+ scale = gsl_min_fminimizer_x_minimum(s);
+ lo = gsl_min_fminimizer_x_lower(s);
+ up = gsl_min_fminimizer_x_upper(s);
+
+ /* Check for convergence */
+ status = gsl_min_test_interval(lo, up, 0.001, 0.0);
+
+ iter++;
+
+ } while ( status == GSL_CONTINUE );
+
+ if ( status != GSL_SUCCESS ) {
+ ERROR("Scale factor minimisation failed.\n");
+ }
+
+ gsl_min_fminimizer_free(s);
+
+ }
+
+ //STATUS("Final scale factor: %5.2e\n", scale);
+ *scalep = scale;
+ return calc_r(scale, &rp);
+}
+
+
+double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_1_IGNORE, u);
+}
+
+
+double stat_r1_zero(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_1_ZERO, u);
+}
+
+
+double stat_r2(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_2, u);
+}
+
+
+double stat_r1_i(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_1_I, u);
+}
+
+
+double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_DIFF_ZERO, u);
+}
+
+
+double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_DIFF_IGNORE, u);
+}
+
+
+double stat_rdiff_intensity(RefList *list1, RefList *list2, double *scalep, int u)
+{
+ return r_minimised(list1, list2, scalep, R_DIFF_INTENSITY, u);
+}
+
+
+double stat_pearson_i(RefList *list1, RefList *list2)
+{
+ double *vec1, *vec2;
+ int ni = num_reflections(list1);
+ double val;
+ int nacc = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ vec1 = malloc(ni*sizeof(double));
+ vec2 = malloc(ni*sizeof(double));
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ vec1[nacc] = i1;
+ vec2[nacc] = i2;
+ nacc++;
+ }
+
+ val = gsl_stats_correlation(vec1, 1, vec2, 1, nacc);
+ free(vec1);
+ free(vec2);
+
+ return val;
+}
+
+
+double stat_pearson_f_ignore(RefList *list1, RefList *list2)
+{
+ double *vec1, *vec2;
+ int ni = num_reflections(list1);
+ double val;
+ int nacc = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ vec1 = malloc(ni*sizeof(double));
+ vec2 = malloc(ni*sizeof(double));
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ if ( i1 < 0.0 ) continue;
+ if ( i2 < 0.0 ) continue;
+
+ f1 = sqrt(i1);
+ f2 = sqrt(i2);
+
+ vec1[nacc] = f1;
+ vec2[nacc] = f2;
+ nacc++;
+
+ }
+
+ val = gsl_stats_correlation(vec1, 1, vec2, 1, nacc);
+ free(vec1);
+ free(vec2);
+
+ return val;
+}
+
+
+double stat_pearson_f_zero(RefList *list1, RefList *list2)
+{
+ double *vec1, *vec2;
+ int ni = num_reflections(list1);
+ double val;
+ int nacc = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
+
+ vec1 = malloc(ni*sizeof(double));
+ vec2 = malloc(ni*sizeof(double));
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ double i1, i2;
+ double f1, f2;
+ signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
+
+ f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
+ f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
+
+ vec1[nacc] = f1;
+ vec2[nacc] = f2;
+ nacc++;
+
+ }
+
+ val = gsl_stats_correlation(vec1, 1, vec2, 1, nacc);
+ free(vec1);
+ free(vec2);
+
+ return val;
+}
diff --git a/libcrystfel/src/statistics.h b/libcrystfel/src/statistics.h
new file mode 100644
index 00000000..8fa78ea6
--- /dev/null
+++ b/libcrystfel/src/statistics.h
@@ -0,0 +1,45 @@
+/*
+ * statistics.h
+ *
+ * Structure-factor statistics
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef STATISTICS_H
+#define STATISTICS_H
+
+
+#include "reflist.h"
+
+extern double stat_scale_intensity(RefList *list1, RefList *list2);
+
+extern double stat_r1_zero(RefList *list1, RefList *list2,
+ double *scalep, int u);
+extern double stat_r1_ignore(RefList *list1, RefList *list2,
+ double *scalep, int u);
+
+extern double stat_r2(RefList *list1, RefList *list2, double *scalep, int u);
+
+extern double stat_r1_i(RefList *list1, RefList *list2, double *scalep, int u);
+
+extern double stat_rdiff_zero(RefList *list1, RefList *list2,
+ double *scalep, int u);
+extern double stat_rdiff_ignore(RefList *list1, RefList *list2,
+ double *scalep, int u);
+extern double stat_rdiff_intensity(RefList *list1, RefList *list2,
+ double *scalep, int u);
+
+extern double stat_pearson_i(RefList *list1, RefList *list2);
+extern double stat_pearson_f_zero(RefList *list1, RefList *list2);
+extern double stat_pearson_f_ignore(RefList *list1, RefList *list2);
+
+
+#endif /* STATISTICS_H */
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
new file mode 100644
index 00000000..a7cdc2d9
--- /dev/null
+++ b/libcrystfel/src/stream.c
@@ -0,0 +1,487 @@
+/*
+ * stream.c
+ *
+ * Stream tools
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ * (c) 2011 Rick Kirian <rkirian@asu.edu>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "cell.h"
+#include "utils.h"
+#include "image.h"
+#include "stream.h"
+#include "reflist.h"
+#include "reflist-utils.h"
+
+
+#define CHUNK_START_MARKER "----- Begin chunk -----"
+#define CHUNK_END_MARKER "----- End chunk -----"
+#define PEAK_LIST_START_MARKER "Peaks from peak search"
+#define PEAK_LIST_END_MARKER "End of peak list"
+#define REFLECTION_START_MARKER "Reflections measured after indexing"
+/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also
+ * used to terminate a standalone list of reflections */
+
+static void exclusive(const char *a, const char *b)
+{
+ ERROR("The stream options '%s' and '%s' are mutually exclusive.\n",
+ a, b);
+}
+
+
+int parse_stream_flags(const char *a)
+{
+ int n, i;
+ int ret = STREAM_NONE;
+ char **flags;
+
+ n = assplode(a, ",", &flags, ASSPLODE_NONE);
+
+ for ( i=0; i<n; i++ ) {
+
+ if ( strcmp(flags[i], "integrated") == 0) {
+
+ ret |= STREAM_INTEGRATED;
+
+ } else if ( strcmp(flags[i], "peaks") == 0) {
+ if ( ret & STREAM_PEAKS_IF_INDEXED ) {
+ exclusive("peaks", "peaksifindexed");
+ return -1;
+ }
+ if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) {
+ exclusive("peaks", "peaksifnotindexed");
+ return -1;
+ }
+ ret |= STREAM_PEAKS;
+
+ } else if ( strcmp(flags[i], "peaksifindexed") == 0) {
+ if ( ret & STREAM_PEAKS ) {
+ exclusive("peaks", "peaksifindexed");
+ return -1;
+ }
+ if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) {
+ exclusive("peaksifnotindexed",
+ "peaksifindexed");
+ return -1;
+ }
+ ret |= STREAM_PEAKS_IF_INDEXED;
+
+ } else if ( strcmp(flags[i], "peaksifnotindexed") == 0) {
+ if ( ret & STREAM_PEAKS ) {
+ exclusive("peaks", "peaksifnotindexed");
+ return -1;
+ }
+ if ( ret & STREAM_PEAKS_IF_INDEXED ) {
+ exclusive("peaksifnotindexed",
+ "peaksifindexed");
+ return -1;
+ }
+ ret |= STREAM_PEAKS_IF_NOT_INDEXED;
+
+ } else {
+ ERROR("Unrecognised stream flag '%s'\n", flags[i]);
+ return -1;
+ }
+
+ free(flags[i]);
+
+ }
+ free(flags);
+
+ return ret;
+}
+
+
+int count_patterns(FILE *fh)
+{
+ char *rval;
+
+ int n_total_patterns = 0;
+ do {
+ char line[1024];
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ chomp(line);
+ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++;
+
+ } while ( rval != NULL );
+
+ if ( ferror(fh) ) {
+ ERROR("Read error while counting patterns.\n");
+ return 0;
+ }
+
+ return n_total_patterns;
+}
+
+
+static int read_peaks(FILE *fh, struct image *image)
+{
+ char *rval = NULL;
+ int first = 1;
+
+ image->features = image_feature_list_new();
+
+ do {
+
+ char line[1024];
+ float x, y, d, intensity;
+ int r;
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ chomp(line);
+
+ if ( strcmp(line, PEAK_LIST_END_MARKER) == 0 ) return 0;
+
+ r = sscanf(line, "%f %f %f %f", &x, &y, &d, &intensity);
+ if ( (r != 4) && (!first) ) {
+ ERROR("Failed to parse peak list line.\n");
+ ERROR("The failed line was: '%s'\n", line);
+ return 1;
+ }
+
+ first = 0;
+ if ( r == 4 ) {
+ image_add_feature(image->features, x, y,
+ image, intensity, NULL);
+ }
+
+ } while ( rval != NULL );
+
+ /* Got read error of some kind before finding PEAK_LIST_END_MARKER */
+ return 1;
+}
+
+
+static void write_peaks(struct image *image, FILE *ofh)
+{
+ int i;
+
+ fprintf(ofh, PEAK_LIST_START_MARKER"\n");
+ fprintf(ofh, " fs/px ss/px (1/d)/nm^-1 Intensity\n");
+
+ for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+ struct imagefeature *f;
+ struct rvec r;
+ double q;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
+ q = modulus(r.u, r.v, r.w);
+
+ fprintf(ofh, "%7.2f %7.2f %10.2f %10.2f\n",
+ f->fs, f->ss, q/1.0e9, f->intensity);
+
+ }
+
+ fprintf(ofh, PEAK_LIST_END_MARKER"\n");
+}
+
+
+void write_chunk(FILE *ofh, struct image *i, struct hdfile *hdfile, int f)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double a, b, c, al, be, ga;
+
+ fprintf(ofh, CHUNK_START_MARKER"\n");
+
+ fprintf(ofh, "Image filename: %s\n", i->filename);
+
+ if ( i->indexed_cell != NULL ) {
+
+ cell_get_parameters(i->indexed_cell, &a, &b, &c,
+ &al, &be, &ga);
+ fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm,"
+ " %7.5f %7.5f %7.5f deg\n",
+ a*1.0e9, b*1.0e9, c*1.0e9,
+ rad2deg(al), rad2deg(be), rad2deg(ga));
+
+ cell_get_reciprocal(i->indexed_cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n",
+ asx/1e9, asy/1e9, asz/1e9);
+ fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
+ bsx/1e9, bsy/1e9, bsz/1e9);
+ fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
+ csx/1e9, csy/1e9, csz/1e9);
+
+ } else {
+
+ fprintf(ofh, "No unit cell from indexing.\n");
+
+ }
+
+ if ( i->i0_available ) {
+ fprintf(ofh, "I0 = %7.5f (arbitrary units)\n", i->i0);
+ } else {
+ fprintf(ofh, "I0 = invalid\n");
+ }
+
+ fprintf(ofh, "photon_energy_eV = %f\n",
+ J_to_eV(ph_lambda_to_en(i->lambda)));
+
+ if ( i->det != NULL ) {
+
+ int j;
+
+ for ( j=0; j<i->det->n_panels; j++ ) {
+ fprintf(ofh, "camera_length_%s = %f\n",
+ i->det->panels[j].name, i->det->panels[j].clen);
+ }
+
+ }
+
+ copy_hdf5_fields(hdfile, i->copyme, ofh);
+
+ if ( (f & STREAM_PEAKS)
+ || ((f & STREAM_PEAKS_IF_INDEXED) && (i->indexed_cell != NULL))
+ || ((f & STREAM_PEAKS_IF_NOT_INDEXED) && (i->indexed_cell == NULL)) )
+ {
+ fprintf(ofh, "\n");
+ write_peaks(i, ofh);
+ }
+
+ if ( f & STREAM_INTEGRATED ) {
+
+ fprintf(ofh, "\n");
+
+ if ( i->reflections != NULL ) {
+
+ fprintf(ofh, REFLECTION_START_MARKER"\n");
+ write_reflections_to_file(ofh, i->reflections,
+ i->indexed_cell);
+ fprintf(ofh, REFLECTION_END_MARKER"\n");
+
+ } else {
+
+ fprintf(ofh, "No integrated reflections.\n");
+
+ }
+ }
+
+ fprintf(ofh, CHUNK_END_MARKER"\n\n");
+}
+
+
+static int find_start_of_chunk(FILE *fh)
+{
+ char *rval = NULL;
+ char line[1024];
+
+ do {
+
+ rval = fgets(line, 1023, fh);
+
+ /* Trouble? */
+ if ( rval == NULL ) return 1;
+
+ chomp(line);
+
+ } while ( strcmp(line, CHUNK_START_MARKER) != 0 );
+
+ return 0;
+}
+
+
+/* Read the next chunk from a stream and fill in 'image' */
+int read_chunk(FILE *fh, struct image *image)
+{
+ char line[1024];
+ char *rval = NULL;
+ struct rvec as, bs, cs;
+ int have_as = 0;
+ int have_bs = 0;
+ int have_cs = 0;
+ int have_filename = 0;
+ int have_cell = 0;
+ int have_ev = 0;
+
+ if ( find_start_of_chunk(fh) ) return 1;
+
+ image->i0_available = 0;
+ image->i0 = 1.0;
+ image->lambda = -1.0;
+ image->features = NULL;
+ image->reflections = NULL;
+ image->indexed_cell = NULL;
+
+ do {
+
+ float u, v, w;
+
+ rval = fgets(line, 1023, fh);
+
+ /* Trouble? */
+ if ( rval == NULL ) break;
+
+ chomp(line);
+
+ if ( strncmp(line, "Image filename: ", 16) == 0 ) {
+ image->filename = strdup(line+16);
+ have_filename = 1;
+ }
+
+ if ( strncmp(line, "camera_length_", 14) == 0 ) {
+ if ( image->det != NULL ) {
+
+ int k;
+ char name[1024];
+ struct panel *p;
+
+ for ( k=0; k<strlen(line)-14; k++ ) {
+ char ch = line[k+14];
+ name[k] = ch;
+ if ( (ch == ' ') || (ch == '=') ) {
+ name[k] = '\0';
+ break;
+ }
+ }
+
+ p = find_panel_by_name(image->det, name);
+ if ( p == NULL ) {
+ ERROR("No panel '%s'\n", name);
+ } else {
+ p->clen = atof(line+14+k+3);
+ }
+
+ }
+ }
+
+ if ( strncmp(line, "I0 = ", 5) == 0 ) {
+ image->i0 = atof(line+5);
+ image->i0_available = 1;
+ }
+
+ if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) == 3 ) {
+ as.u = u*1e9; as.v = v*1e9; as.w = w*1e9;
+ have_as = 1;
+ }
+
+ if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) == 3 ) {
+ bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9;
+ have_bs = 1;
+ }
+
+ if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) == 3 ) {
+ cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9;
+ have_cs = 1;
+ }
+
+ if ( have_as && have_bs && have_cs ) {
+ if ( image->indexed_cell != NULL ) {
+ ERROR("Duplicate cell found in stream!\n");
+ cell_free(image->indexed_cell);
+ }
+ image->indexed_cell = cell_new_from_reciprocal_axes(as,
+ bs,
+ cs);
+ have_cell = 1;
+ have_as = 0; have_bs = 0; have_cs = 0;
+ }
+
+ if ( strncmp(line, "photon_energy_eV = ", 19) == 0 ) {
+ image->lambda = ph_en_to_lambda(eV_to_J(atof(line+19)));
+ have_ev = 1;
+ }
+
+ if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) {
+ if ( read_peaks(fh, image) ) {
+ ERROR("Failed while reading peaks\n");
+ return 1;
+ }
+ }
+
+ if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) {
+ image->reflections = read_reflections_from_file(fh);
+ if ( image->reflections == NULL ) {
+ ERROR("Failed while reading reflections\n");
+ return 1;
+ }
+ }
+
+ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) {
+ if ( have_filename && have_ev ) return 0;
+ ERROR("Incomplete chunk found in input file.\n");
+ return 1;
+ }
+
+ } while ( 1 );
+
+ return 1; /* Either error or EOF, don't care because we will complain
+ * on the terminal if it was an error. */
+}
+
+
+void write_stream_header(FILE *ofh, int argc, char *argv[])
+{
+ int i;
+
+ fprintf(ofh, "CrystFEL stream format 2.0\n");
+ fprintf(ofh, "Command line:");
+ for ( i=0; i<argc; i++ ) {
+ fprintf(ofh, " %s", argv[i]);
+ }
+ fprintf(ofh, "\n");
+ fflush(ofh);
+}
+
+
+int skip_some_files(FILE *fh, int n)
+{
+ char *rval = NULL;
+ int n_patterns = 0;
+
+ do {
+
+ char line[1024];
+
+ if ( n_patterns == n ) return 0;
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
+
+ } while ( rval != NULL );
+
+ return 1;
+}
+
+int is_stream(const char *filename) {
+ FILE *fh;
+ char line[1024];
+ char *rval = NULL;
+ fh = fopen(filename, "r");
+ rval = fgets(line, 1023, fh);
+ fclose(fh);
+ if ( rval == NULL ) {
+ return -1;
+ }
+ if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) {
+ return 1;
+ }
+ else {
+ return 0;
+ }
+ return -1;
+}
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
new file mode 100644
index 00000000..ba218fb9
--- /dev/null
+++ b/libcrystfel/src/stream.h
@@ -0,0 +1,49 @@
+/*
+ * stream.h
+ *
+ * Stream tools
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef STREAM_H
+#define STREAM_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+struct image;
+struct hdfile;
+
+/* Possible options dictating what goes into the output stream */
+enum
+{
+ STREAM_NONE = 0,
+ STREAM_INTEGRATED = 1<<0,
+ STREAM_PEAKS = 1<<2,
+ STREAM_PEAKS_IF_INDEXED = 1<<3,
+ STREAM_PEAKS_IF_NOT_INDEXED = 1<<4,
+};
+
+
+extern int count_patterns(FILE *fh);
+
+extern void write_stream_header(FILE *ofh, int argc, char *argv[]);
+
+extern void write_chunk(FILE *ofh, struct image *image, struct hdfile *hdfile,
+ int flags);
+
+extern int parse_stream_flags(const char *a);
+
+extern int read_chunk(FILE *fh, struct image *image);
+
+extern int skip_some_files(FILE *fh, int n);
+
+extern int is_stream(const char *filename);
+
+#endif /* STREAM_H */
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
new file mode 100644
index 00000000..f0b24146
--- /dev/null
+++ b/libcrystfel/src/symmetry.c
@@ -0,0 +1,1503 @@
+/*
+ * symmetry.c
+ *
+ * Symmetry
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#include "symmetry.h"
+#include "utils.h"
+
+
+/**
+ * SECTION:symmetry
+ * @short_description: Point symmetry handling
+ * @title: Symmetry
+ * @section_id:
+ * @see_also:
+ * @include: "symmetry.h"
+ * @Image:
+ *
+ * Routines to handle point symmetry.
+ */
+
+
+struct sym_op
+{
+ signed int *h;
+ signed int *k;
+ signed int *l; /* Contributions to h, k and l from h, k, i and l */
+ int order;
+};
+
+
+struct _symoplist
+{
+ struct sym_op *ops;
+ int n_ops;
+ int max_ops;
+ char *name;
+ int *divisors;
+ int num_equivs;
+};
+
+
+struct _symopmask
+{
+ const SymOpList *list;
+ int *mask;
+};
+
+
+
+static void alloc_ops(SymOpList *ops)
+{
+ ops->ops = realloc(ops->ops, ops->max_ops*sizeof(struct sym_op));
+ ops->divisors = realloc(ops->divisors, ops->max_ops*sizeof(int));
+}
+
+
+/**
+ * new_symopmask:
+ * @list: A %SymOpList
+ *
+ * Returns: a new %SymOpMask, which you can use when filtering out special
+ * reflections.
+ **/
+SymOpMask *new_symopmask(const SymOpList *list)
+{
+ SymOpMask *m;
+ int i;
+
+ m = malloc(sizeof(struct _symopmask));
+ if ( m == NULL ) return NULL;
+
+ m->list = list;
+ m->mask = malloc(sizeof(int)*list->n_ops);
+ if ( m->mask == NULL ) {
+ free(m);
+ return NULL;
+ }
+
+ for ( i=0; i<list->n_ops; i++ ) {
+ m->mask[i] = 1;
+ }
+
+ return m;
+}
+
+
+/* Creates a new SymOpList */
+static SymOpList *new_symoplist()
+{
+ SymOpList *new;
+ new = malloc(sizeof(SymOpList));
+ if ( new == NULL ) return NULL;
+ new->max_ops = 16;
+ new->n_ops = 0;
+ new->ops = NULL;
+ new->divisors = NULL;
+ new->name = NULL;
+ new->num_equivs = 1;
+ alloc_ops(new);
+ return new;
+}
+
+
+/**
+ * free_symoplist:
+ * @ops: A %SymOpList to free
+ *
+ * Frees a %SymOpList and all associated resources.
+ **/
+void free_symoplist(SymOpList *ops)
+{
+ int i;
+
+ if ( ops == NULL ) return;
+ for ( i=0; i<ops->n_ops; i++ ) {
+ free(ops->ops[i].h);
+ free(ops->ops[i].k);
+ free(ops->ops[i].l);
+ }
+ if ( ops->ops != NULL ) free(ops->ops);
+ if ( ops->name != NULL ) free(ops->name);
+ free(ops);
+}
+
+/**
+ * free_symopmask:
+ * @m: A %SymOpMask to free
+ *
+ * Frees a %SymOpMask and all associated resources.
+ **/
+void free_symopmask(SymOpMask *m)
+{
+ if ( m == NULL ) return;
+ free(m->mask);
+ free(m);
+}
+
+
+/* This returns the number of operations in "ops". This might be different
+ * to num_equivs() if the point group is being constructed. */
+static int num_ops(const SymOpList *ops)
+{
+ return ops->n_ops;
+}
+
+
+/* Add a operation to a SymOpList */
+static void add_symop(SymOpList *ops,
+ signed int *h, signed int *k, signed int *l,
+ int order)
+{
+ int n;
+
+ if ( ops->n_ops == ops->max_ops ) {
+ /* Pretty sure this never happens, but still... */
+ ops->max_ops += 16;
+ alloc_ops(ops);
+ }
+
+ n = ops->n_ops;
+ ops->ops[n].h = h;
+ ops->ops[n].k = k;
+ ops->ops[n].l = l;
+ ops->ops[n].order = order;
+ ops->n_ops++;
+}
+
+
+/* Add a operation to a SymOpList */
+static void add_copied_op(SymOpList *ops, struct sym_op *copyme)
+{
+ int n;
+ signed int *h, *k, *l;
+
+ if ( ops->n_ops == ops->max_ops ) {
+ ops->max_ops += 16;
+ alloc_ops(ops);
+ }
+
+ n = ops->n_ops;
+
+ h = malloc(3*sizeof(signed int));
+ k = malloc(3*sizeof(signed int));
+ l = malloc(3*sizeof(signed int));
+
+ memcpy(h, copyme->h, 3*sizeof(signed int));
+ memcpy(k, copyme->k, 3*sizeof(signed int));
+ memcpy(l, copyme->l, 3*sizeof(signed int));
+
+ ops->ops[n].h = h;
+ ops->ops[n].k = k;
+ ops->ops[n].l = l;
+ ops->ops[n].order = copyme->order;
+
+ ops->n_ops++;
+}
+
+
+/**
+ * num_equivs:
+ * @ops: A %SymOpList
+ * @m: A %SymOpMask, which has been shown to special_position()
+ *
+ * Returns: the number of equivalent reflections for a general reflection
+ * in point group "ops", which were not flagged by your call to
+ * special_position().
+ **/
+int num_equivs(const SymOpList *ops, const SymOpMask *m)
+{
+ int n = num_ops(ops);
+ int i;
+ int c;
+
+ if ( m == NULL ) return n;
+
+ c = 0;
+ for ( i=0; i<n; i++ ) {
+ if ( m->mask[i] ) c++;
+ }
+
+ return c;
+}
+
+
+static signed int *v(signed int h, signed int k, signed int i, signed int l)
+{
+ signed int *vec = malloc(3*sizeof(signed int));
+ /* Convert back to 3-index form now */
+ vec[0] = h-i; vec[1] = k-i; vec[2] = l;
+ return vec;
+}
+
+
+static void combine_ops(signed int *h1, signed int *k1, signed int *l1,
+ signed int *h2, signed int *k2, signed int *l2,
+ signed int *hnew, signed int *knew, signed int *lnew)
+{
+ /* Yay matrices */
+ hnew[0] = h1[0]*h2[0] + h1[1]*k2[0] + h1[2]*l2[0];
+ hnew[1] = h1[0]*h2[1] + h1[1]*k2[1] + h1[2]*l2[1];
+ hnew[2] = h1[0]*h2[2] + h1[1]*k2[2] + h1[2]*l2[2];
+
+ knew[0] = k1[0]*h2[0] + k1[1]*k2[0] + k1[2]*l2[0];
+ knew[1] = k1[0]*h2[1] + k1[1]*k2[1] + k1[2]*l2[1];
+ knew[2] = k1[0]*h2[2] + k1[1]*k2[2] + k1[2]*l2[2];
+
+ lnew[0] = l1[0]*h2[0] + l1[1]*k2[0] + l1[2]*l2[0];
+ lnew[1] = l1[0]*h2[1] + l1[1]*k2[1] + l1[2]*l2[1];
+ lnew[2] = l1[0]*h2[2] + l1[1]*k2[2] + l1[2]*l2[2];
+}
+
+
+static void combine_and_add_symop(struct sym_op *opi, int oi,
+ struct sym_op *opj,
+ SymOpList *s)
+{
+ int i;
+ signed int *h, *k, *l;
+
+ h = malloc(3*sizeof(signed int));
+ k = malloc(3*sizeof(signed int));
+ l = malloc(3*sizeof(signed int));
+ assert(h != NULL);
+ assert(k != NULL);
+ assert(l != NULL);
+
+ memcpy(h, opj->h, 3*sizeof(signed int));
+ memcpy(k, opj->k, 3*sizeof(signed int));
+ memcpy(l, opj->l, 3*sizeof(signed int));
+
+ for ( i=0; i<oi; i++ ) {
+
+ signed int hfs[3], kfs[3], lfs[3];
+
+ combine_ops(h, k, l, opi->h, opi->k, opi->l, hfs, kfs, lfs);
+
+ memcpy(h, hfs, 3*sizeof(signed int));
+ memcpy(k, kfs, 3*sizeof(signed int));
+ memcpy(l, lfs, 3*sizeof(signed int));
+
+ }
+
+// STATUS("Creating %3i %3i %3i\n", h[0], h[1], h[2]);
+// STATUS(" %3i %3i %3i\n", k[0], k[1], k[2]);
+// STATUS(" %3i %3i %3i\n", l[0], l[1], l[2]);
+
+ add_symop(s, h, k, l, 1);
+}
+
+
+/* Fill in the other operations for a point group starting from its
+ * generators */
+static SymOpList *expand_ops(SymOpList *s)
+{
+ int n, i;
+ SymOpList *e;
+
+ e = new_symoplist();
+ if ( e == NULL ) return NULL;
+ e->name = strdup(symmetry_name(s));
+
+ add_symop(e, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1), 1); /* I */
+
+ n = num_ops(s);
+ for ( i=0; i<n; i++ ) {
+
+ int j, nj;
+ struct sym_op *opi = &s->ops[i];
+
+ /* Apply op 'i' to all the current ops in the list */
+ nj = num_ops(e);
+ for ( j=0; j<nj; j++ ) {
+
+ int oi;
+
+ for ( oi=0; oi<opi->order-1; oi++ ) {
+ combine_and_add_symop(opi, oi+1, &e->ops[j], e);
+ }
+
+ }
+
+ }
+
+ free_symoplist(s);
+
+ return e;
+}
+
+
+/********************************* Triclinic **********************************/
+
+static SymOpList *make_1bar()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ new->name = strdup("-1");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_1()
+{
+ SymOpList *new = new_symoplist();
+ new->name = strdup("1");
+ return expand_ops(new);
+}
+
+
+/********************************* Monoclinic *********************************/
+
+static SymOpList *make_2m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
+ new->name = strdup("2/m");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_2_uab()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
+ new->name = strdup("2_uab");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_2()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
+ new->name = strdup("2");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
+ new->name = strdup("m");
+ return expand_ops(new);
+}
+
+
+/******************************** Orthorhombic ********************************/
+
+static SymOpList *make_mmm()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
+ add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m -| k */
+ new->name = strdup("mmm");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_222()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
+ new->name = strdup("222");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_mm2()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
+ add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m -| k */
+ new->name = strdup("mm2");
+ return expand_ops(new);
+}
+
+
+/********************************* Tetragonal *********************************/
+
+static SymOpList *make_4m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
+ new->name = strdup("4/m");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ new->name = strdup("4");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4mm()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| l */
+ new->name = strdup("4mm");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_422()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
+ new->name = strdup("422");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4bar()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
+ new->name = strdup("-4");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4bar2m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
+ new->name = strdup("-42m");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4barm2()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2); /* 2 // h+k */
+ new->name = strdup("-4m2");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4mmm()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| k */
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
+ new->name = strdup("4/mmm");
+ return expand_ops(new);
+}
+
+
+/************************** Trigonal (Rhombohedral) ***************************/
+
+static SymOpList *make_3_R()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* 3 // h+k+l */
+ new->name = strdup("3_R");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3bar_R()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* -3 // h+k+l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ new->name = strdup("-3_R");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_32_R()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* 3 // h+k+l */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1), 2); /* 2 -| 3 */
+ new->name = strdup("32_R");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3m_R()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* 3 // h+k+l */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2); /* m */
+ new->name = strdup("3m_R");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3barm_R()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* -3 // h+k+l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2); /* m */
+ new->name = strdup("-3m_R");
+ return expand_ops(new);
+}
+
+
+/*************************** Trigonal (Hexagonal) *****************************/
+
+static SymOpList *make_3_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ new->name = strdup("3_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3bar_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ new->name = strdup("-3_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_321_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2); /* 2 // h */
+ new->name = strdup("321_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_312_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1), 2); /* 2 // h+k */
+ new->name = strdup("312_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3m1_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
+ new->name = strdup("3m1_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_31m_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2); /* m -| (k+i) */
+ new->name = strdup("31m_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3barm1_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2); /* 2 // h */
+ new->name = strdup("-3m1_H");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_3bar1m_H()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1), 2); /* 2 // h+k */
+ new->name = strdup("-31m_H");
+ return expand_ops(new);
+}
+
+
+/********************************** Hexgonal **********************************/
+
+static SymOpList *make_6()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
+ new->name = strdup("6");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_6bar()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
+ new->name = strdup("-6");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_6m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
+ add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
+ new->name = strdup("6/m");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_622()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2); /* 2 // h */
+ new->name = strdup("622");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_6mm()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
+ new->name = strdup("6mm");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_6barm2()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
+ new->name = strdup("-6m2");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_6bar2m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
+ add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2); /* m -| (k+i) */
+ new->name = strdup("-62m");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_6mmm()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
+ add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ new->name = strdup("6/mmm");
+ return expand_ops(new);
+}
+
+
+/************************************ Cubic ***********************************/
+
+static SymOpList *make_23()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2// l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2// k */
+ add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3// h+k+l */
+ new->name = strdup("23");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_m3bar()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2// l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2// k */
+ add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3// h+k+l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ new->name = strdup("m-3");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_432()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);/* 2 // k */
+ add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3 // h+k+l */
+ new->name = strdup("432");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_4bar3m()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);/* 2 // k */
+ add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3 // h+k+l */
+ new->name = strdup("-43m");
+ return expand_ops(new);
+}
+
+
+static SymOpList *make_m3barm()
+{
+ SymOpList *new = new_symoplist();
+ add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
+ add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);/* 2 // k */
+ add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3 // h+k+l */
+ add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
+ new->name = strdup("m-3m");
+ return expand_ops(new);
+}
+
+
+/**
+ * get_pointgroup:
+ * @sym: A string representation of a point group
+ *
+ * This function parses @sym and returns the corresponding %SymOpList.
+ * In the string representation of the point group, use a preceding minus sign
+ * for any character which would have a "bar". Trigonal groups must be suffixed
+ * with either "_H" or "_R" for a hexagonal or rhombohedral lattice
+ * respectively.
+ *
+ * Examples: -1 1 2/m 2 m mmm 222 mm2 4/m 4 -4 4/mmm 422 -42m -4m2 4mm
+ * 3_R -3_R 32_R 3m_R -3m_R 3_H -3_H 321_H 312_H 3m1_H 31m_H -3m1_H -31m_H
+ * 6/m 6 -6 6/mmm 622 -62m -6m2 6mm 23 m-3 432 -43m m-3m.
+ **/
+SymOpList *get_pointgroup(const char *sym)
+{
+ /* Triclinic */
+ if ( strcmp(sym, "-1") == 0 ) return make_1bar();
+ if ( strcmp(sym, "1") == 0 ) return make_1();
+
+ /* Monoclinic */
+ if ( strcmp(sym, "2/m") == 0 ) return make_2m();
+ if ( strcmp(sym, "2") == 0 ) return make_2();
+ if ( strcmp(sym, "2_uab") == 0 ) return make_2_uab();
+ if ( strcmp(sym, "m") == 0 ) return make_m();
+
+ /* Orthorhombic */
+ if ( strcmp(sym, "mmm") == 0 ) return make_mmm();
+ if ( strcmp(sym, "222") == 0 ) return make_222();
+ if ( strcmp(sym, "mm2") == 0 ) return make_mm2();
+
+ /* Tetragonal */
+ if ( strcmp(sym, "4/m") == 0 ) return make_4m();
+ if ( strcmp(sym, "4") == 0 ) return make_4();
+ if ( strcmp(sym, "-4") == 0 ) return make_4bar();
+ if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm();
+ if ( strcmp(sym, "422") == 0 ) return make_422();
+ if ( strcmp(sym, "-42m") == 0 ) return make_4bar2m();
+ if ( strcmp(sym, "-4m2") == 0 ) return make_4barm2();
+ if ( strcmp(sym, "4mm") == 0 ) return make_4mm();
+
+ /* Trigonal (rhombohedral) */
+ if ( strcmp(sym, "3_R") == 0 ) return make_3_R();
+ if ( strcmp(sym, "-3_R") == 0 ) return make_3bar_R();
+ if ( strcmp(sym, "32_R") == 0 ) return make_32_R();
+ if ( strcmp(sym, "3m_R") == 0 ) return make_3m_R();
+ if ( strcmp(sym, "-3m_R") == 0 ) return make_3barm_R();
+
+ /* Trigonal (hexagonal) */
+ if ( strcmp(sym, "3_H") == 0 ) return make_3_H();
+ if ( strcmp(sym, "-3_H") == 0 ) return make_3bar_H();
+ if ( strcmp(sym, "321_H") == 0 ) return make_321_H();
+ if ( strcmp(sym, "312_H") == 0 ) return make_312_H();
+ if ( strcmp(sym, "3m1_H") == 0 ) return make_3m1_H();
+ if ( strcmp(sym, "31m_H") == 0 ) return make_31m_H();
+ if ( strcmp(sym, "-3m1_H") == 0 ) return make_3barm1_H();
+ if ( strcmp(sym, "-31m_H") == 0 ) return make_3bar1m_H();
+
+ /* Hexagonal */
+ if ( strcmp(sym, "6/m") == 0 ) return make_6m();
+ if ( strcmp(sym, "6") == 0 ) return make_6();
+ if ( strcmp(sym, "-6") == 0 ) return make_6bar();
+ if ( strcmp(sym, "6/mmm") == 0 ) return make_6mmm();
+ if ( strcmp(sym, "622") == 0 ) return make_622();
+ if ( strcmp(sym, "-62m") == 0 ) return make_6bar2m();
+ if ( strcmp(sym, "-6m2") == 0 ) return make_6barm2();
+ if ( strcmp(sym, "6mm") == 0 ) return make_6mm();
+
+ /* Cubic */
+ if ( strcmp(sym, "23") == 0 ) return make_23();
+ if ( strcmp(sym, "m-3") == 0 ) return make_m3bar();
+ if ( strcmp(sym, "432") == 0 ) return make_432();
+ if ( strcmp(sym, "-43m") == 0 ) return make_4bar3m();
+ if ( strcmp(sym, "m-3m") == 0 ) return make_m3barm();
+
+ ERROR("Unknown point group '%s'\n", sym);
+ return NULL;
+}
+
+
+static void do_op(const struct sym_op *op,
+ signed int h, signed int k, signed int l,
+ signed int *he, signed int *ke, signed int *le)
+{
+ *he = h*op->h[0] + k*op->h[1] + l*op->h[2];
+ *ke = h*op->k[0] + k*op->k[1] + l*op->k[2];
+ *le = h*op->l[0] + k*op->l[1] + l*op->l[2];
+}
+
+
+/**
+ * get_equiv:
+ * @ops: A %SymOpList
+ * @m: A %SymOpMask, which has been shown to special_position()
+ * @idx: Index of the operation to use
+ * @h: index of reflection
+ * @k: index of reflection
+ * @l: index of reflection
+ * @he: location to store h index of equivalent reflection
+ * @ke: location to store k index of equivalent reflection
+ * @le: location to store l index of equivalent reflection
+ *
+ * This function applies the @idx-th symmetry operation from @ops to the
+ * reflection @h, @k, @l, and stores the result at @he, @ke and @le.
+ *
+ * If you don't mind that the same equivalent might appear twice, simply call
+ * this function the number of times returned by num_ops(), using the actual
+ * point group. If repeating the same equivalent twice (for example, if the
+ * given reflection is a special high-symmetry one), call special_position()
+ * first to get a "specialised" SymOpList and use that instead.
+ **/
+void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
+ signed int h, signed int k, signed int l,
+ signed int *he, signed int *ke, signed int *le)
+{
+ const int n = num_ops(ops);
+
+ if ( m != NULL ) {
+
+ int i, c;
+
+ c = 0;
+ for ( i=0; i<n; i++ ) {
+
+ if ( (c == idx) && m->mask[i] ) {
+ do_op(&ops->ops[i], h, k, l, he, ke, le);
+ return;
+ }
+
+ if ( m->mask[i] ) {
+ c++;
+ }
+
+ }
+
+ ERROR("Index %i out of range for point group '%s' with"
+ " reflection %i %i %i\n",
+ idx, symmetry_name(ops), h, k, l);
+
+ *he = 0; *ke = 0; *le = 0;
+
+ return;
+
+ }
+
+
+
+ if ( idx >= n ) {
+
+ ERROR("Index %i out of range for point group '%s'\n", idx,
+ symmetry_name(ops));
+
+ *he = 0; *ke = 0; *le = 0;
+ return;
+
+ }
+
+ do_op(&ops->ops[idx], h, k, l, he, ke, le);
+}
+
+
+/**
+ * special_position:
+ * @ops: A %SymOpList, usually corresponding to a point group
+ * @m: A %SymOpMask created with new_symopmask()
+ * @h: index of a reflection
+ * @k: index of a reflection
+ * @l: index of a reflection
+ *
+ * This function determines which operations in @ops map the reflection @h, @k,
+ * @l onto itself, and uses @m to flag the operations in @ops which cause this.
+ *
+ **/
+void special_position(const SymOpList *ops, SymOpMask *m,
+ signed int h, signed int k, signed int l)
+{
+ int i, n;
+ signed int *htest;
+ signed int *ktest;
+ signed int *ltest;
+
+ assert(m->list = ops);
+
+ n = num_equivs(ops, NULL);
+ htest = malloc(n*sizeof(signed int));
+ ktest = malloc(n*sizeof(signed int));
+ ltest = malloc(n*sizeof(signed int));
+
+ for ( i=0; i<n; i++ ) {
+
+ signed int he, ke, le;
+ int j;
+
+ get_equiv(ops, NULL, i, h, k, l, &he, &ke, &le);
+
+ m->mask[i] = 1;
+ for ( j=0; j<i; j++ ) {
+ if ( (he==htest[j]) && (ke==ktest[j])
+ && (le==ltest[j]) )
+ {
+ m->mask[i] = 0;
+ break; /* Only need to find one */
+ }
+ }
+
+ htest[i] = he;
+ ktest[i] = ke;
+ ltest[i] = le;
+
+ }
+
+ free(htest);
+ free(ktest);
+ free(ltest);
+}
+
+
+static int any_negative(signed int h, signed int k, signed int l)
+{
+ if ( h < 0 ) return 1;
+ if ( k < 0 ) return 1;
+ if ( l < 0 ) return 1;
+ return 0;
+}
+
+
+/**
+ * get_asymm:
+ * @ops: A %SymOpList, usually corresponding to a point group
+ * @h: index of a reflection
+ * @k: index of a reflection
+ * @l: index of a reflection
+ * @hp: location for asymmetric index of reflection
+ * @kp: location for asymmetric index of reflection
+ * @lp: location for asymmetric index of reflection
+ *
+ * This function determines the asymmetric version of the reflection @h, @k, @l
+ * in symmetry group @ops, and puts the result in @hp, @kp, @lp.
+ *
+ * This is a relatively expensive operation because of its generality.
+ * Therefore, if you know you'll need to make repeated use of the asymmetric
+ * indices, consider creating a new %RefList indexed according to the asymmetric
+ * indices themselves with asymmetric_indices(). If you do that, you'll still
+ * be able to get the original versions of the indices with
+ * get_symmetric_indices().
+ *
+ **/
+void get_asymm(const SymOpList *ops,
+ signed int h, signed int k, signed int l,
+ signed int *hp, signed int *kp, signed int *lp)
+{
+ int nequiv;
+ int p;
+ signed int best_h, best_k, best_l;
+ int have_negs;
+
+ nequiv = num_equivs(ops, NULL);
+
+ best_h = h; best_k = k; best_l = l;
+ have_negs = any_negative(best_h, best_k, best_l);
+ for ( p=0; p<nequiv; p++ ) {
+
+ int will_have_negs;
+
+ get_equiv(ops, NULL, p, h, k, l, hp, kp, lp);
+
+ will_have_negs = any_negative(*hp, *kp, *lp);
+
+ /* Don't lose "no negs" status */
+ if ( !have_negs && will_have_negs ) continue;
+
+ if ( have_negs && !will_have_negs ) {
+ best_h = *hp; best_k = *kp; best_l = *lp;
+ have_negs = 0;
+ continue;
+ }
+
+ if ( *hp > best_h ) {
+ best_h = *hp; best_k = *kp; best_l = *lp;
+ have_negs = any_negative(best_h, best_k, best_l);
+ continue;
+ }
+ if ( *hp < best_h ) continue;
+
+ if ( *kp > best_k ) {
+ best_h = *hp; best_k = *kp; best_l = *lp;
+ have_negs = any_negative(best_h, best_k, best_l);
+ continue;
+ }
+ if ( *kp < best_k ) continue;
+
+ if ( *lp > best_l ) {
+ best_h = *hp; best_k = *kp; best_l = *lp;
+ have_negs = any_negative(best_h, best_k, best_l);
+ continue;
+ }
+
+ }
+
+ *hp = best_h; *kp = best_k; *lp = best_l;
+}
+
+
+static int is_inversion(const struct sym_op *op)
+{
+ if ( (op->h[0]!=-1) || (op->h[1]!=0) || (op->h[2]!=0) ) return 0;
+ if ( (op->k[0]!=0) || (op->k[1]!=-1) || (op->k[2]!=0) ) return 0;
+ if ( (op->l[0]!=0) || (op->l[1]!=0) || (op->l[2]!=-1) ) return 0;
+ return 1;
+}
+
+
+static int is_identity(const struct sym_op *op)
+{
+ if ( (op->h[0]!=1) || (op->h[1]!=0) || (op->h[2]!=0) ) return 0;
+ if ( (op->k[0]!=0) || (op->k[1]!=1) || (op->k[2]!=0) ) return 0;
+ if ( (op->l[0]!=0) || (op->l[1]!=0) || (op->l[2]!=1) ) return 0;
+ return 1;
+}
+
+
+static signed int determinant(const struct sym_op *op)
+{
+ signed int det = 0;
+
+ det += op->h[0] * (op->k[1]*op->l[2] - op->k[2]*op->l[1]);
+ det -= op->h[1] * (op->k[0]*op->l[2] - op->k[2]*op->l[0]);
+ det += op->h[2] * (op->k[0]*op->l[1] - op->k[1]*op->l[0]);
+
+ return det;
+}
+
+
+/**
+ * is_centrosymmetric:
+ * @s: A %SymOpList
+ *
+ * Returns: non-zero if @s contains an inversion operation
+ */
+int is_centrosymmetric(const SymOpList *s)
+{
+ int i, n;
+
+ n = num_ops(s);
+ for ( i=0; i<n; i++ ) {
+ if ( is_inversion(&s->ops[i]) ) return 1;
+ }
+
+ return 0;
+}
+
+
+static int ops_equal(const struct sym_op *op,
+ signed int *h, signed int *k, signed int *l)
+{
+ if ( (op->h[0]!=h[0]) || (op->h[1]!=h[1]) || (op->h[2]!=h[2]) )
+ return 0;
+ if ( (op->k[0]!=k[0]) || (op->k[1]!=k[1]) || (op->k[2]!=k[2]) )
+ return 0;
+ if ( (op->l[0]!=l[0]) || (op->l[1]!=l[1]) || (op->l[2]!=l[2]) )
+ return 0;
+ return 1;
+}
+
+
+static int struct_ops_equal(const struct sym_op *op1, const struct sym_op *op2)
+{
+ return ops_equal(op1, op2->h, op2->k, op2->l);
+}
+
+
+/* Return true if a*b = ans */
+static int check_mult(const struct sym_op *ans,
+ const struct sym_op *a, const struct sym_op *b)
+{
+ signed int *ans_h, *ans_k, *ans_l;
+ int val;
+
+ ans_h = malloc(3*sizeof(signed int));
+ ans_k = malloc(3*sizeof(signed int));
+ ans_l = malloc(3*sizeof(signed int));
+
+ combine_ops(a->h, a->k, a->l, b->h, b->k, b->l, ans_h, ans_k, ans_l);
+ val = ops_equal(ans, ans_h, ans_k, ans_l);
+
+ free(ans_h);
+ free(ans_k);
+ free(ans_l);
+
+ return val;
+}
+
+
+/**
+ * is_subgroup:
+ * @source: A %SymOpList
+ * @target: Another %SymOpList, which might be a subgroup of @source.
+ *
+ * Returns: non-zero if every operation in @target is also in @source.
+ **/
+int is_subgroup(const SymOpList *source, const SymOpList *target)
+{
+ int n_src, n_tgt;
+ int i;
+
+ n_src = num_ops(source);
+ n_tgt = num_ops(target);
+
+ for ( i=0; i<n_tgt; i++ ) {
+
+ int j;
+ int found = 0;
+
+ for ( j=0; j<n_src; j++ ) {
+
+ if ( struct_ops_equal(&target->ops[i],
+ &source->ops[j] ) )
+ {
+ found = 1;
+ break;
+ }
+
+ }
+
+ if ( !found ) return 0;
+
+ }
+
+ return 1;
+}
+
+
+/**
+ * get_ambiguities:
+ * @source: The "source" symmetry, a %SymOpList
+ * @target: The "target" symmetry, a %SymOpList
+
+ * Calculates twinning laws. Returns a %SymOpList containing the twinning
+ * operators, which are the symmetry operations which can be added to @target
+ * to generate @source. Only rotations are allowable - no mirrors nor
+ * inversions.
+ * To count the number of possibilities, use num_ops() on the result.
+ *
+ * Returns: A %SymOpList containing the twinning operators, or NULL if the
+ * source symmetry cannot be generated from that target symmetry without using
+ * mirror or inversion operations.
+ */
+SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target)
+{
+ int n_src, n_tgt;
+ int i;
+ SymOpList *twins;
+ SymOpList *src_reordered;
+ SymOpMask *used;
+ char *name;
+ int index;
+
+ n_src = num_ops(source);
+ n_tgt = num_ops(target);
+
+ if ( !is_subgroup(source, target) ) {
+ ERROR("'%s' is not a subgroup of '%s'\n",
+ symmetry_name(target), symmetry_name(source));
+ return NULL;
+ }
+
+ if ( n_src % n_tgt != 0 ) {
+ ERROR("Subgroup index would be fractional.\n");
+ return NULL;
+ }
+ index = n_src / n_tgt;
+
+ src_reordered = new_symoplist();
+ used = new_symopmask(source);
+
+ /* Find identity */
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( is_identity(&source->ops[i]) ) {
+ add_copied_op(src_reordered, &source->ops[i]);
+ used->mask[i] = 0;
+ }
+ }
+
+ /* Find binary options (order=2) of first kind (determinant positive) */
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( (source->ops[i].order == 2)
+ && (determinant(&source->ops[i]) > 0) ) {
+ add_copied_op(src_reordered, &source->ops[i]);
+ used->mask[i] = 0;
+ }
+ }
+
+ /* Find other operations of first kind (determinant positive) */
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( determinant(&source->ops[i]) > 0 ) {
+ add_copied_op(src_reordered, &source->ops[i]);
+ used->mask[i] = 0;
+ }
+ }
+
+ /* Find inversion */
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( is_inversion(&source->ops[i]) ) {
+ add_copied_op(src_reordered, &source->ops[i]);
+ used->mask[i] = 0;
+ }
+ }
+
+ /* Find binary options of second kind (determinant negative) */
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( (source->ops[i].order == 2)
+ && (determinant(&source->ops[i]) < 0) ) {
+ add_copied_op(src_reordered, &source->ops[i]);
+ used->mask[i] = 0;
+ }
+ }
+
+ /* Find other operations of second kind (determinant negative) */
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( determinant(&source->ops[i]) < 0 ) {
+ add_copied_op(src_reordered, &source->ops[i]);
+ used->mask[i] = 0;
+ }
+ }
+
+ int n_left_over = 0;
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ n_left_over++;
+ }
+ if ( n_left_over != 0 ) {
+ ERROR("%i operations left over after rearranging for"
+ " left coset decomposition.\n", n_left_over);
+ }
+
+ if ( num_ops(src_reordered) != num_ops(source) ) {
+ ERROR("%i ops went to %i after rearranging.\n",
+ num_ops(src_reordered), num_ops(source));
+ }
+
+ free_symopmask(used);
+ used = new_symopmask(src_reordered);
+
+ /* This is the first method from Flack (1987) */
+ for ( i=0; i<n_src; i++ ) {
+
+ int j;
+ if ( used->mask[i] == 0 ) continue;
+
+ for ( j=1; j<n_tgt; j++ ) {
+
+ int k;
+ for ( k=i+1; k<n_src; k++ ) {
+ if ( check_mult(&src_reordered->ops[k],
+ &src_reordered->ops[i],
+ &target->ops[j]) )
+ {
+ used->mask[k] = 0;
+ }
+ }
+
+ }
+
+ }
+
+ twins = new_symoplist();
+ for ( i=0; i<n_src; i++ ) {
+ if ( used->mask[i] == 0 ) continue;
+ if ( determinant(&src_reordered->ops[i]) < 0 ) {
+ /* A mirror or inversion turned up in the list.
+ * That means that no pure rotational ambiguity can
+ * account for this subgroup relationship. */
+ free_symoplist(twins);
+ free_symopmask(used);
+ free_symoplist(src_reordered);
+ return NULL;
+ }
+ add_copied_op(twins, &src_reordered->ops[i]);
+ }
+
+ free_symopmask(used);
+ free_symoplist(src_reordered);
+
+ name = malloc(64);
+ snprintf(name, 63, "%s -> %s", symmetry_name(source),
+ symmetry_name(target));
+ twins->name = name;
+
+ return twins;
+}
+
+
+static void add_chars(char *t, const char *s, int max_len)
+{
+ char *tmp;
+
+ tmp = strdup(t);
+
+ snprintf(t, max_len, "%s%s", tmp, s);
+ free(tmp);
+}
+
+
+static char *get_matrix_name(signed int *v)
+{
+ char *text;
+ const int max_len = 9;
+ int i;
+ int printed = 0;
+
+ text = malloc(max_len+1);
+ text[0] = '\0';
+
+ for ( i=0; i<3; i++ ) {
+
+ if ( v[i] == 0 ) continue;
+
+ if ( (i==0) && (v[0]==v[1]) ) {
+ if ( v[i]>0 ) add_chars(text, "-", max_len);
+ add_chars(text, "i", max_len);
+ v[1] -= v[0];
+ continue;
+ }
+
+ if ( printed ) add_chars(text, "+", max_len);
+
+ if ( v[i]<0 ) add_chars(text, "-", max_len);
+
+ if ( abs(v[i])>1 ) {
+ char num[3];
+ snprintf(num, 2, "%i", abs(v[i]));
+ add_chars(text, num, max_len);
+ }
+
+ switch ( i )
+ {
+ case 0 : add_chars(text, "h", max_len); break;
+ case 1 : add_chars(text, "k", max_len); break;
+ case 2 : add_chars(text, "l", max_len); break;
+ default : add_chars(text, "X", max_len); break;
+ }
+
+ printed = 1;
+
+ }
+
+ return text;
+}
+
+
+static char *name_equiv(const struct sym_op *op)
+{
+ char *h, *k, *l;
+ char *name;
+
+ h = get_matrix_name(op->h);
+ k = get_matrix_name(op->k);
+ l = get_matrix_name(op->l);
+ name = malloc(32);
+
+ snprintf(name, 31, "%s%s%s", h, k, l);
+ free(h);
+ free(k);
+ free(l);
+
+ return name;
+}
+
+
+/**
+ * describe_symmetry:
+ * @s: A %SymOpList
+ *
+ * Writes the name and a list of operations to stderr.
+ */
+void describe_symmetry(const SymOpList *s)
+{
+ int i, n;
+
+ n = num_equivs(s, NULL);
+
+ STATUS("%15s :", symmetry_name(s));
+
+ for ( i=0; i<n; i++ ) {
+ char *name = name_equiv(&s->ops[i]);
+ STATUS(" %6s", name);
+ free(name);
+ if ( (i!=0) && (i%8==0) ) STATUS("\n%15s ", "");
+ }
+ STATUS("\n");
+}
+
+
+/**
+ * symmetry_name:
+ * @ops: A %SymOpList
+ *
+ * Returns: a text description of @ops.
+ */
+const char *symmetry_name(const SymOpList *ops)
+{
+ return ops->name;
+}
diff --git a/libcrystfel/src/symmetry.h b/libcrystfel/src/symmetry.h
new file mode 100644
index 00000000..071ebbde
--- /dev/null
+++ b/libcrystfel/src/symmetry.h
@@ -0,0 +1,63 @@
+/*
+ * symmetry.h
+ *
+ * Symmetry
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef SYMMETRY_H
+#define SYMMETRY_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/**
+ * SymOpList
+ *
+ * The SymOpList is an opaque data structure containing a list of point symmetry
+ * operations. It could represent an point group or a list of indexing
+ * ambiguities (twin laws), or similar.
+ **/
+typedef struct _symoplist SymOpList;
+
+/**
+ * SymOpMask
+ *
+ * The SymOpMask is an opaque data structure containing a list of flags
+ * associated with point symmetry operations in a specific %SymOpList. It is
+ * used to filter the operations in the %SymOpList to avoid duplicating
+ * equivalent reflections when the reflection is somehow special (e.g. 'hk0').
+ **/
+typedef struct _symopmask SymOpMask;
+
+extern void free_symoplist(SymOpList *ops);
+extern SymOpList *get_pointgroup(const char *sym);
+
+extern SymOpMask *new_symopmask(const SymOpList *list);
+extern void free_symopmask(SymOpMask *m);
+
+extern void special_position(const SymOpList *ops, SymOpMask *m,
+ signed int h, signed int k, signed int l);
+extern void get_asymm(const SymOpList *ops,
+ signed int h, signed int k, signed int l,
+ signed int *hp, signed int *kp, signed int *lp);
+extern int num_equivs(const SymOpList *ops, const SymOpMask *m);
+extern void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
+ signed int h, signed int k, signed int l,
+ signed int *he, signed int *ke, signed int *le);
+
+extern SymOpList *get_ambiguities(const SymOpList *source,
+ const SymOpList *target);
+extern int is_subgroup(const SymOpList *source, const SymOpList *target);
+
+extern int is_centrosymmetric(const SymOpList *s);
+extern const char *symmetry_name(const SymOpList *ops);
+extern void describe_symmetry(const SymOpList *s);
+
+#endif /* SYMMETRY_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 18:00:14 +0000