From 469efb904b59f137ac9e85e5ff23edd0c113de5c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 12:17:59 +0100 Subject: Move a load more stuff into libcrystfel --- libcrystfel/Makefile.am | 8 +- libcrystfel/Makefile.in | 54 +- libcrystfel/src/cl-utils.c | 200 +++++ libcrystfel/src/cl-utils.h | 27 + libcrystfel/src/diffraction-gpu.c | 529 +++++++++++++ libcrystfel/src/diffraction-gpu.h | 57 ++ libcrystfel/src/diffraction.c | 463 ++++++++++++ libcrystfel/src/diffraction.h | 34 + libcrystfel/src/filters.c | 130 ++++ libcrystfel/src/filters.h | 25 + libcrystfel/src/geometry.c | 341 +++++++++ libcrystfel/src/geometry.h | 26 + libcrystfel/src/peaks.c | 593 +++++++++++++++ libcrystfel/src/peaks.h | 38 + libcrystfel/src/reflist-utils.c | 502 +++++++++++++ libcrystfel/src/reflist-utils.h | 52 ++ libcrystfel/src/statistics.c | 668 +++++++++++++++++ libcrystfel/src/statistics.h | 45 ++ libcrystfel/src/stream.c | 487 ++++++++++++ libcrystfel/src/stream.h | 49 ++ libcrystfel/src/symmetry.c | 1503 +++++++++++++++++++++++++++++++++++++ libcrystfel/src/symmetry.h | 63 ++ 22 files changed, 5891 insertions(+), 3 deletions(-) create mode 100644 libcrystfel/src/cl-utils.c create mode 100644 libcrystfel/src/cl-utils.h create mode 100644 libcrystfel/src/diffraction-gpu.c create mode 100644 libcrystfel/src/diffraction-gpu.h create mode 100644 libcrystfel/src/diffraction.c create mode 100644 libcrystfel/src/diffraction.h create mode 100644 libcrystfel/src/filters.c create mode 100644 libcrystfel/src/filters.h create mode 100644 libcrystfel/src/geometry.c create mode 100644 libcrystfel/src/geometry.h create mode 100644 libcrystfel/src/peaks.c create mode 100644 libcrystfel/src/peaks.h create mode 100644 libcrystfel/src/reflist-utils.c create mode 100644 libcrystfel/src/reflist-utils.h create mode 100644 libcrystfel/src/statistics.c create mode 100644 libcrystfel/src/statistics.h create mode 100644 libcrystfel/src/stream.c create mode 100644 libcrystfel/src/stream.h create mode 100644 libcrystfel/src/symmetry.c create mode 100644 libcrystfel/src/symmetry.h (limited to 'libcrystfel') diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index daf2ef33..304eb7bb 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -1,5 +1,11 @@ lib_LTLIBRARIES = libcrystfel.la libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/thread-pool.c src/image.c src/hdf5-file.c \ - src/beam-parameters.c + src/beam-parameters.c src/geometry.c src/statistics.c \ + src/symmetry.c src/stream.c src/peaks.c \ + src/reflist-utils.c src/filters.c src/diffraction.c \ + src/diffraction-gpu.c src/cl-utils.c + INCLUDES = "-I$(top_srcdir)/data" +AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib +AM_CPPFLAGS += -I$(top_srcdir)/lib diff --git a/libcrystfel/Makefile.in b/libcrystfel/Makefile.in index 3465345d..9464b1e6 100644 --- a/libcrystfel/Makefile.in +++ b/libcrystfel/Makefile.in @@ -83,7 +83,10 @@ libcrystfel_la_LIBADD = am__dirstamp = $(am__leading_dot)dirstamp am_libcrystfel_la_OBJECTS = src/reflist.lo src/utils.lo src/cell.lo \ src/detector.lo src/thread-pool.lo src/image.lo \ - src/hdf5-file.lo src/beam-parameters.lo + src/hdf5-file.lo src/beam-parameters.lo src/geometry.lo \ + src/statistics.lo src/symmetry.lo src/stream.lo src/peaks.lo \ + src/reflist-utils.lo src/filters.lo src/diffraction.lo \ + src/diffraction-gpu.lo src/cl-utils.lo libcrystfel_la_OBJECTS = $(am_libcrystfel_la_OBJECTS) AM_V_lt = $(am__v_lt_$(V)) am__v_lt_ = $(am__v_lt_$(AM_DEFAULT_VERBOSITY)) @@ -341,9 +344,14 @@ top_srcdir = @top_srcdir@ lib_LTLIBRARIES = libcrystfel.la libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/thread-pool.c src/image.c src/hdf5-file.c \ - src/beam-parameters.c + src/beam-parameters.c src/geometry.c src/statistics.c \ + src/symmetry.c src/stream.c src/peaks.c \ + src/reflist-utils.c src/filters.c src/diffraction.c \ + src/diffraction-gpu.c src/cl-utils.c INCLUDES = "-I$(top_srcdir)/data" +AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib \ + -I$(top_srcdir)/lib all: all-am .SUFFIXES: @@ -424,6 +432,18 @@ src/image.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/hdf5-file.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/beam-parameters.lo: src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/geometry.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/statistics.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/symmetry.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/stream.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/peaks.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/reflist-utils.lo: src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +src/filters.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/diffraction.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/diffraction-gpu.lo: src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +src/cl-utils.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) libcrystfel.la: $(libcrystfel_la_OBJECTS) $(libcrystfel_la_DEPENDENCIES) $(AM_V_CCLD)$(LINK) -rpath $(libdir) $(libcrystfel_la_OBJECTS) $(libcrystfel_la_LIBADD) $(LIBS) @@ -433,14 +453,34 @@ mostlyclean-compile: -rm -f src/beam-parameters.lo -rm -f src/cell.$(OBJEXT) -rm -f src/cell.lo + -rm -f src/cl-utils.$(OBJEXT) + -rm -f src/cl-utils.lo -rm -f src/detector.$(OBJEXT) -rm -f src/detector.lo + -rm -f src/diffraction-gpu.$(OBJEXT) + -rm -f src/diffraction-gpu.lo + -rm -f src/diffraction.$(OBJEXT) + -rm -f src/diffraction.lo + -rm -f src/filters.$(OBJEXT) + -rm -f src/filters.lo + -rm -f src/geometry.$(OBJEXT) + -rm -f src/geometry.lo -rm -f src/hdf5-file.$(OBJEXT) -rm -f src/hdf5-file.lo -rm -f src/image.$(OBJEXT) -rm -f src/image.lo + -rm -f src/peaks.$(OBJEXT) + -rm -f src/peaks.lo + -rm -f src/reflist-utils.$(OBJEXT) + -rm -f src/reflist-utils.lo -rm -f src/reflist.$(OBJEXT) -rm -f src/reflist.lo + -rm -f src/statistics.$(OBJEXT) + -rm -f src/statistics.lo + -rm -f src/stream.$(OBJEXT) + -rm -f src/stream.lo + -rm -f src/symmetry.$(OBJEXT) + -rm -f src/symmetry.lo -rm -f src/thread-pool.$(OBJEXT) -rm -f src/thread-pool.lo -rm -f src/utils.$(OBJEXT) @@ -451,10 +491,20 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/beam-parameters.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/cell.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/cl-utils.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/detector.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction-gpu.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/filters.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/geometry.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/hdf5-file.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/image.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/peaks.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist-utils.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/statistics.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/stream.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/symmetry.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/thread-pool.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/utils.Plo@am__quote@ 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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#ifdef HAVE_CL_CL_H +#include +#else +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef CLUTILS_H +#define CLUTILS_H + +#ifdef HAVE_CONFIG_H +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#ifdef HAVE_CL_CL_H +#include +#else +#include +#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; isinc_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; isinc_lut_ptrs[n-1][i] = 1.0; + } + } else { + for ( i=1; isinc_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; idet->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; fsmin_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; iintensities = 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; iflags = 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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include +#include +#include +#include +#include +#include +#include + +#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= 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; fswidth; fs++ ) { + for ( ss=0; ssheight; ss++ ) { + + int fs_step, ss_step, kstep; + int divxval, divyval; + int idx = fs + image->width*ss; + + for ( fs_step=0; fs_stepdata[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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#ifdef GSL_FUDGE +#include +#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; xwidth; x++ ) { + for ( y=0; yheight; 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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef FILTERS_H +#define FILTERS_H + +#ifdef HAVE_CONFIG_H +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include + +#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; in_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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef GEOMETRY_H +#define GEOMETRY_H + + +#ifdef HAVE_CONFIG_H +#include +#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 + * 2011 Andrew Martin + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include + +#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; ifeatures, 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; jfeatures, 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; jfeatures, 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; idet->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; idet->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; ifeatures); 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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef PEAKS_H +#define PEAKS_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include +#include + + +#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 1 ) { + + STATUS("Found %i %i %i: %i times:\n", h, k, l, found); + + for ( j=0; j *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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#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 + * (c) 2011 Rick Kirian + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include + +#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; ifeatures = 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; ifeatures); 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; jdet->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; kdet, 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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef STREAM_H +#define STREAM_H + +#ifdef HAVE_CONFIG_H +#include +#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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#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; in_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; in_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; imask[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; ih, 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; iops[i]; + + /* Apply op 'i' to all the current ops in the list */ + nj = num_ops(e); + for ( j=0; jorder-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; imask[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; imask[i] = 1; + for ( j=0; jmask[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 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; iops[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; iops[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; imask[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; imask[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; imask[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; imask[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; imask[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; imask[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; imask[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; imask[i] == 0 ) continue; + + for ( j=1; jops[k], + &src_reordered->ops[i], + &target->ops[j]) ) + { + used->mask[k] = 0; + } + } + + } + + } + + twins = new_symoplist(); + for ( i=0; imask[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; iops[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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef SYMMETRY_H +#define SYMMETRY_H + +#ifdef HAVE_CONFIG_H +#include +#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 */ -- cgit v1.2.3