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 --- src/cl-utils.c | 200 ------- src/cl-utils.h | 27 - src/diffraction-gpu.c | 529 ----------------- src/diffraction-gpu.h | 57 -- src/diffraction.c | 463 --------------- src/diffraction.h | 34 -- src/filters.c | 130 ----- src/filters.h | 25 - src/geometry.c | 341 ----------- src/geometry.h | 26 - src/peaks.c | 594 ------------------- src/peaks.h | 38 -- src/reflist-utils.c | 502 ----------------- src/reflist-utils.h | 52 -- src/statistics.c | 668 ---------------------- src/statistics.h | 45 -- src/stream.c | 487 ---------------- src/stream.h | 49 -- src/symmetry.c | 1503 ------------------------------------------------- src/symmetry.h | 63 --- 20 files changed, 5833 deletions(-) delete mode 100644 src/cl-utils.c delete mode 100644 src/cl-utils.h delete mode 100644 src/diffraction-gpu.c delete mode 100644 src/diffraction-gpu.h delete mode 100644 src/diffraction.c delete mode 100644 src/diffraction.h delete mode 100644 src/filters.c delete mode 100644 src/filters.h delete mode 100644 src/geometry.c delete mode 100644 src/geometry.h delete mode 100644 src/peaks.c delete mode 100644 src/peaks.h delete mode 100644 src/reflist-utils.c delete mode 100644 src/reflist-utils.h delete mode 100644 src/statistics.c delete mode 100644 src/statistics.h delete mode 100644 src/stream.c delete mode 100644 src/stream.h delete mode 100644 src/symmetry.c delete mode 100644 src/symmetry.h (limited to 'src') diff --git a/src/cl-utils.c b/src/cl-utils.c deleted file mode 100644 index 65a09363..00000000 --- a/src/cl-utils.c +++ /dev/null @@ -1,200 +0,0 @@ -/* - * 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/src/diffraction-gpu.c b/src/diffraction-gpu.c deleted file mode 100644 index 605b1514..00000000 --- a/src/diffraction-gpu.c +++ /dev/null @@ -1,529 +0,0 @@ -/* - * 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/src/diffraction-gpu.h b/src/diffraction-gpu.h deleted file mode 100644 index a3bde4e1..00000000 --- a/src/diffraction-gpu.h +++ /dev/null @@ -1,57 +0,0 @@ -/* - * 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/src/diffraction.c b/src/diffraction.c deleted file mode 100644 index 9532a6ce..00000000 --- a/src/diffraction.c +++ /dev/null @@ -1,463 +0,0 @@ -/* - * 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/src/diffraction.h b/src/diffraction.h deleted file mode 100644 index f71d3cce..00000000 --- a/src/diffraction.h +++ /dev/null @@ -1,34 +0,0 @@ -/* - * 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/src/filters.c b/src/filters.c deleted file mode 100644 index c4e409df..00000000 --- a/src/filters.c +++ /dev/null @@ -1,130 +0,0 @@ -/* - * 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/src/filters.h b/src/filters.h deleted file mode 100644 index 6b35d7e8..00000000 --- a/src/filters.h +++ /dev/null @@ -1,25 +0,0 @@ -/* - * 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/src/geometry.c b/src/geometry.c deleted file mode 100644 index 485abba3..00000000 --- a/src/geometry.c +++ /dev/null @@ -1,341 +0,0 @@ -/* - * 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/src/geometry.h b/src/geometry.h deleted file mode 100644 index ddf04b80..00000000 --- a/src/geometry.h +++ /dev/null @@ -1,26 +0,0 @@ -/* - * 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/src/peaks.c b/src/peaks.c deleted file mode 100644 index f1a58d23..00000000 --- a/src/peaks.c +++ /dev/null @@ -1,594 +0,0 @@ -/* - * 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 "index.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/src/peaks.h b/src/peaks.h deleted file mode 100644 index 9d475ea9..00000000 --- a/src/peaks.h +++ /dev/null @@ -1,38 +0,0 @@ -/* - * 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/src/reflist-utils.c b/src/reflist-utils.c deleted file mode 100644 index b64e9979..00000000 --- a/src/reflist-utils.c +++ /dev/null @@ -1,502 +0,0 @@ -/* - * 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/src/reflist-utils.h b/src/reflist-utils.h deleted file mode 100644 index d14e5f8e..00000000 --- a/src/reflist-utils.h +++ /dev/null @@ -1,52 +0,0 @@ -/* - * 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/src/statistics.c b/src/statistics.c deleted file mode 100644 index 7990672d..00000000 --- a/src/statistics.c +++ /dev/null @@ -1,668 +0,0 @@ -/* - * 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/src/statistics.h b/src/statistics.h deleted file mode 100644 index 8fa78ea6..00000000 --- a/src/statistics.h +++ /dev/null @@ -1,45 +0,0 @@ -/* - * 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/src/stream.c b/src/stream.c deleted file mode 100644 index a7cdc2d9..00000000 --- a/src/stream.c +++ /dev/null @@ -1,487 +0,0 @@ -/* - * 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/src/symmetry.c b/src/symmetry.c deleted file mode 100644 index f0b24146..00000000 --- a/src/symmetry.c +++ /dev/null @@ -1,1503 +0,0 @@ -/* - * 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/src/symmetry.h b/src/symmetry.h deleted file mode 100644 index 071ebbde..00000000 --- a/src/symmetry.h +++ /dev/null @@ -1,63 +0,0 @@ -/* - * 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