diff options
| author | Thomas White <taw@physics.org> | 2011-11-15 12:05:55 +0100 |
|---|---|---|
| committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:40 +0100 |
| commit | 38089071300b8e04ed42236dd08d9055094fb3b8 (patch) | |
| tree | 91e1487ac820eb549e7652750867cd4fec039097 /libcrystfel | |
| parent | 404c612223dbfa0210902ebc5c9226927335aa65 (diff) | |
Introduce "libcrystfel"
Diffstat (limited to 'libcrystfel')
| -rw-r--r-- | libcrystfel/Makefile.am | 5 | ||||
| -rw-r--r-- | libcrystfel/Makefile.in | 701 | ||||
| -rw-r--r-- | libcrystfel/src/beam-parameters.c | 133 | ||||
| -rw-r--r-- | libcrystfel/src/beam-parameters.h | 40 | ||||
| -rw-r--r-- | libcrystfel/src/cell.c | 1164 | ||||
| -rw-r--r-- | libcrystfel/src/cell.h | 107 | ||||
| -rw-r--r-- | libcrystfel/src/detector.c | 1233 | ||||
| -rw-r--r-- | libcrystfel/src/detector.h | 131 | ||||
| -rw-r--r-- | libcrystfel/src/hdf5-file.c | 817 | ||||
| -rw-r--r-- | libcrystfel/src/hdf5-file.h | 55 | ||||
| -rw-r--r-- | libcrystfel/src/image.c | 144 | ||||
| -rw-r--r-- | libcrystfel/src/image.h | 184 | ||||
| -rw-r--r-- | libcrystfel/src/list_tmp.h | 106 | ||||
| -rw-r--r-- | libcrystfel/src/reflist.c | 1048 | ||||
| -rw-r--r-- | libcrystfel/src/reflist.h | 112 | ||||
| -rw-r--r-- | libcrystfel/src/thread-pool.c | 283 | ||||
| -rw-r--r-- | libcrystfel/src/thread-pool.h | 71 | ||||
| -rw-r--r-- | libcrystfel/src/utils.c | 488 | ||||
| -rw-r--r-- | libcrystfel/src/utils.h | 236 |
19 files changed, 7058 insertions, 0 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am new file mode 100644 index 00000000..daf2ef33 --- /dev/null +++ b/libcrystfel/Makefile.am @@ -0,0 +1,5 @@ +lib_LTLIBRARIES = libcrystfel.la +libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ + src/thread-pool.c src/image.c src/hdf5-file.c \ + src/beam-parameters.c +INCLUDES = "-I$(top_srcdir)/data" diff --git a/libcrystfel/Makefile.in b/libcrystfel/Makefile.in new file mode 100644 index 00000000..3465345d --- /dev/null +++ b/libcrystfel/Makefile.in @@ -0,0 +1,701 @@ +# Makefile.in generated by automake 1.11.1 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, +# Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + +@SET_MAKE@ + +VPATH = @srcdir@ +pkgdatadir = $(datadir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkglibexecdir = $(libexecdir)/@PACKAGE@ +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +build_triplet = @build@ +host_triplet = @host@ +subdir = libcrystfel +DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/m4/00gnulib.m4 \ + $(top_srcdir)/m4/extensions.m4 \ + $(top_srcdir)/m4/gnulib-common.m4 \ + $(top_srcdir)/m4/gnulib-comp.m4 $(top_srcdir)/m4/gtk-doc.m4 \ + $(top_srcdir)/m4/include_next.m4 \ + $(top_srcdir)/m4/lib-ignore.m4 $(top_srcdir)/m4/libtool.m4 \ + $(top_srcdir)/m4/ltoptions.m4 $(top_srcdir)/m4/ltsugar.m4 \ + $(top_srcdir)/m4/ltversion.m4 $(top_srcdir)/m4/lt~obsolete.m4 \ + $(top_srcdir)/m4/memchr.m4 $(top_srcdir)/m4/mmap-anon.m4 \ + $(top_srcdir)/m4/onceonly.m4 $(top_srcdir)/m4/stddef_h.m4 \ + $(top_srcdir)/m4/string_h.m4 $(top_srcdir)/m4/strndup.m4 \ + $(top_srcdir)/m4/strnlen.m4 $(top_srcdir)/m4/warn-on-use.m4 \ + $(top_srcdir)/m4/wchar_t.m4 $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(install_sh) -d +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +CONFIG_CLEAN_VPATH_FILES = +am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; +am__vpath_adj = case $$p in \ + $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \ + *) f=$$p;; \ + esac; +am__strip_dir = f=`echo $$p | sed -e 's|^.*/||'`; +am__install_max = 40 +am__nobase_strip_setup = \ + srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*|]/\\\\&/g'` +am__nobase_strip = \ + for p in $$list; do echo "$$p"; done | sed -e "s|$$srcdirstrip/||" +am__nobase_list = $(am__nobase_strip_setup); \ + for p in $$list; do echo "$$p $$p"; done | \ + sed "s| $$srcdirstrip/| |;"' / .*\//!s/ .*/ ./; s,\( .*\)/[^/]*$$,\1,' | \ + $(AWK) 'BEGIN { files["."] = "" } { files[$$2] = files[$$2] " " $$1; \ + if (++n[$$2] == $(am__install_max)) \ + { print $$2, files[$$2]; n[$$2] = 0; files[$$2] = "" } } \ + END { for (dir in files) print dir, files[dir] }' +am__base_list = \ + sed '$$!N;$$!N;$$!N;$$!N;$$!N;$$!N;$$!N;s/\n/ /g' | \ + sed '$$!N;$$!N;$$!N;$$!N;s/\n/ /g' +am__installdirs = "$(DESTDIR)$(libdir)" +LTLIBRARIES = $(lib_LTLIBRARIES) +libcrystfel_la_LIBADD = +am__dirstamp = $(am__leading_dot)dirstamp +am_libcrystfel_la_OBJECTS = src/reflist.lo src/utils.lo src/cell.lo \ + src/detector.lo src/thread-pool.lo src/image.lo \ + src/hdf5-file.lo src/beam-parameters.lo +libcrystfel_la_OBJECTS = $(am_libcrystfel_la_OBJECTS) +AM_V_lt = $(am__v_lt_$(V)) +am__v_lt_ = $(am__v_lt_$(AM_DEFAULT_VERBOSITY)) +am__v_lt_0 = --silent +DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir) +depcomp = $(SHELL) $(top_srcdir)/build-aux/depcomp +am__depfiles_maybe = depfiles +am__mv = mv -f +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) \ + $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) \ + $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \ + $(AM_CFLAGS) $(CFLAGS) +AM_V_CC = $(am__v_CC_$(V)) +am__v_CC_ = $(am__v_CC_$(AM_DEFAULT_VERBOSITY)) +am__v_CC_0 = @echo " CC " $@; +AM_V_at = $(am__v_at_$(V)) +am__v_at_ = $(am__v_at_$(AM_DEFAULT_VERBOSITY)) +am__v_at_0 = @ +CCLD = $(CC) +LINK = $(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) \ + $(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \ + $(AM_LDFLAGS) $(LDFLAGS) -o $@ +AM_V_CCLD = $(am__v_CCLD_$(V)) +am__v_CCLD_ = $(am__v_CCLD_$(AM_DEFAULT_VERBOSITY)) +am__v_CCLD_0 = @echo " CCLD " $@; +AM_V_GEN = $(am__v_GEN_$(V)) +am__v_GEN_ = $(am__v_GEN_$(AM_DEFAULT_VERBOSITY)) +am__v_GEN_0 = @echo " GEN " $@; +SOURCES = $(libcrystfel_la_SOURCES) +DIST_SOURCES = $(libcrystfel_la_SOURCES) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AM_DEFAULT_VERBOSITY = @AM_DEFAULT_VERBOSITY@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CCDEPMODE = @CCDEPMODE@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +Cairo_CFLAGS = @Cairo_CFLAGS@ +Cairo_LIBS = @Cairo_LIBS@ +DEFS = @DEFS@ +DEPDIR = @DEPDIR@ +DSYMUTIL = @DSYMUTIL@ +DUMPBIN = @DUMPBIN@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +FFTW_CFLAGS = @FFTW_CFLAGS@ +FFTW_LIBS = @FFTW_LIBS@ +FGREP = @FGREP@ +GDK_pixbuf_2_CFLAGS = @GDK_pixbuf_2_CFLAGS@ +GDK_pixbuf_2_LIBS = @GDK_pixbuf_2_LIBS@ +GDK_pixbuf_CFLAGS = @GDK_pixbuf_CFLAGS@ +GDK_pixbuf_LIBS = @GDK_pixbuf_LIBS@ +GNULIB_MBSCASECMP = @GNULIB_MBSCASECMP@ +GNULIB_MBSCASESTR = @GNULIB_MBSCASESTR@ +GNULIB_MBSCHR = @GNULIB_MBSCHR@ +GNULIB_MBSCSPN = @GNULIB_MBSCSPN@ +GNULIB_MBSLEN = @GNULIB_MBSLEN@ +GNULIB_MBSNCASECMP = @GNULIB_MBSNCASECMP@ +GNULIB_MBSNLEN = @GNULIB_MBSNLEN@ +GNULIB_MBSPBRK = @GNULIB_MBSPBRK@ +GNULIB_MBSPCASECMP = @GNULIB_MBSPCASECMP@ +GNULIB_MBSRCHR = @GNULIB_MBSRCHR@ +GNULIB_MBSSEP = @GNULIB_MBSSEP@ +GNULIB_MBSSPN = @GNULIB_MBSSPN@ +GNULIB_MBSSTR = @GNULIB_MBSSTR@ +GNULIB_MBSTOK_R = @GNULIB_MBSTOK_R@ +GNULIB_MEMCHR = @GNULIB_MEMCHR@ +GNULIB_MEMMEM = @GNULIB_MEMMEM@ +GNULIB_MEMPCPY = @GNULIB_MEMPCPY@ +GNULIB_MEMRCHR = @GNULIB_MEMRCHR@ +GNULIB_RAWMEMCHR = @GNULIB_RAWMEMCHR@ +GNULIB_STPCPY = @GNULIB_STPCPY@ +GNULIB_STPNCPY = @GNULIB_STPNCPY@ +GNULIB_STRCASESTR = @GNULIB_STRCASESTR@ +GNULIB_STRCHRNUL = @GNULIB_STRCHRNUL@ +GNULIB_STRDUP = @GNULIB_STRDUP@ +GNULIB_STRERROR = @GNULIB_STRERROR@ +GNULIB_STRERROR_R = @GNULIB_STRERROR_R@ +GNULIB_STRNCAT = @GNULIB_STRNCAT@ +GNULIB_STRNDUP = @GNULIB_STRNDUP@ +GNULIB_STRNLEN = @GNULIB_STRNLEN@ +GNULIB_STRPBRK = @GNULIB_STRPBRK@ +GNULIB_STRSEP = @GNULIB_STRSEP@ +GNULIB_STRSIGNAL = @GNULIB_STRSIGNAL@ +GNULIB_STRSTR = @GNULIB_STRSTR@ +GNULIB_STRTOK_R = @GNULIB_STRTOK_R@ +GNULIB_STRVERSCMP = @GNULIB_STRVERSCMP@ +GREP = @GREP@ +GSL_CFLAGS = @GSL_CFLAGS@ +GSL_LIBS = @GSL_LIBS@ +GTKDOC_CHECK = @GTKDOC_CHECK@ +GTKDOC_MKPDF = @GTKDOC_MKPDF@ +GTKDOC_REBASE = @GTKDOC_REBASE@ +GTK_CFLAGS = @GTK_CFLAGS@ +GTK_LIBS = @GTK_LIBS@ +HAVE_DECL_MEMMEM = @HAVE_DECL_MEMMEM@ +HAVE_DECL_MEMRCHR = @HAVE_DECL_MEMRCHR@ +HAVE_DECL_STRDUP = @HAVE_DECL_STRDUP@ +HAVE_DECL_STRERROR_R = @HAVE_DECL_STRERROR_R@ +HAVE_DECL_STRNDUP = @HAVE_DECL_STRNDUP@ +HAVE_DECL_STRNLEN = @HAVE_DECL_STRNLEN@ +HAVE_DECL_STRSIGNAL = @HAVE_DECL_STRSIGNAL@ +HAVE_DECL_STRTOK_R = @HAVE_DECL_STRTOK_R@ +HAVE_MBSLEN = @HAVE_MBSLEN@ +HAVE_MEMCHR = @HAVE_MEMCHR@ +HAVE_MEMPCPY = @HAVE_MEMPCPY@ +HAVE_RAWMEMCHR = @HAVE_RAWMEMCHR@ +HAVE_STPCPY = @HAVE_STPCPY@ +HAVE_STPNCPY = @HAVE_STPNCPY@ +HAVE_STRCASESTR = @HAVE_STRCASESTR@ +HAVE_STRCHRNUL = @HAVE_STRCHRNUL@ +HAVE_STRPBRK = @HAVE_STRPBRK@ +HAVE_STRSEP = @HAVE_STRSEP@ +HAVE_STRVERSCMP = @HAVE_STRVERSCMP@ +HAVE_WCHAR_T = @HAVE_WCHAR_T@ +HTML_DIR = @HTML_DIR@ +IGNORE_UNUSED_LIBRARIES_CFLAGS = @IGNORE_UNUSED_LIBRARIES_CFLAGS@ +INCLUDE_NEXT = @INCLUDE_NEXT@ +INCLUDE_NEXT_AS_FIRST_DIRECTIVE = @INCLUDE_NEXT_AS_FIRST_DIRECTIVE@ +INSTALL = @INSTALL@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LD = @LD@ +LDFLAGS = @LDFLAGS@ +LIBGNU_LIBDEPS = @LIBGNU_LIBDEPS@ +LIBGNU_LTLIBDEPS = @LIBGNU_LTLIBDEPS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LIPO = @LIPO@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAKEINFO = @MAKEINFO@ +MKDIR_P = @MKDIR_P@ +NEXT_AS_FIRST_DIRECTIVE_STDDEF_H = @NEXT_AS_FIRST_DIRECTIVE_STDDEF_H@ +NEXT_AS_FIRST_DIRECTIVE_STRING_H = @NEXT_AS_FIRST_DIRECTIVE_STRING_H@ +NEXT_STDDEF_H = @NEXT_STDDEF_H@ +NEXT_STRING_H = @NEXT_STRING_H@ +NM = @NM@ +NMEDIT = @NMEDIT@ +OBJDUMP = @OBJDUMP@ +OBJEXT = @OBJEXT@ +OTOOL = @OTOOL@ +OTOOL64 = @OTOOL64@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_URL = @PACKAGE_URL@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +PKG_CONFIG = @PKG_CONFIG@ +PKG_CONFIG_LIBDIR = @PKG_CONFIG_LIBDIR@ +PKG_CONFIG_PATH = @PKG_CONFIG_PATH@ +PRAGMA_COLUMNS = @PRAGMA_COLUMNS@ +PRAGMA_SYSTEM_HEADER = @PRAGMA_SYSTEM_HEADER@ +RANLIB = @RANLIB@ +REPLACE_MEMCHR = @REPLACE_MEMCHR@ +REPLACE_MEMMEM = @REPLACE_MEMMEM@ +REPLACE_NULL = @REPLACE_NULL@ +REPLACE_STPNCPY = @REPLACE_STPNCPY@ +REPLACE_STRCASESTR = @REPLACE_STRCASESTR@ +REPLACE_STRDUP = @REPLACE_STRDUP@ +REPLACE_STRERROR = @REPLACE_STRERROR@ +REPLACE_STRERROR_R = @REPLACE_STRERROR_R@ +REPLACE_STRNCAT = @REPLACE_STRNCAT@ +REPLACE_STRNDUP = @REPLACE_STRNDUP@ +REPLACE_STRNLEN = @REPLACE_STRNLEN@ +REPLACE_STRSIGNAL = @REPLACE_STRSIGNAL@ +REPLACE_STRSTR = @REPLACE_STRSTR@ +REPLACE_STRTOK_R = @REPLACE_STRTOK_R@ +SED = @SED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STDDEF_H = @STDDEF_H@ +STRIP = @STRIP@ +UNDEFINE_STRTOK_R = @UNDEFINE_STRTOK_R@ +VERSION = @VERSION@ +abs_builddir = @abs_builddir@ +abs_srcdir = @abs_srcdir@ +abs_top_builddir = @abs_top_builddir@ +abs_top_srcdir = @abs_top_srcdir@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_DUMPBIN = @ac_ct_DUMPBIN@ +am__include = @am__include@ +am__leading_dot = @am__leading_dot@ +am__quote = @am__quote@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +builddir = @builddir@ +datadir = @datadir@ +datarootdir = @datarootdir@ +docdir = @docdir@ +dvidir = @dvidir@ +exec_prefix = @exec_prefix@ +gl_LIBOBJS = @gl_LIBOBJS@ +gl_LTLIBOBJS = @gl_LTLIBOBJS@ +gltests_LIBOBJS = @gltests_LIBOBJS@ +gltests_LTLIBOBJS = @gltests_LTLIBOBJS@ +gltests_WITNESS = @gltests_WITNESS@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +htmldir = @htmldir@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libPNG_CFLAGS = @libPNG_CFLAGS@ +libPNG_LIBS = @libPNG_LIBS@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localedir = @localedir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +pdfdir = @pdfdir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +psdir = @psdir@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +srcdir = @srcdir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +top_build_prefix = @top_build_prefix@ +top_builddir = @top_builddir@ +top_srcdir = @top_srcdir@ +lib_LTLIBRARIES = libcrystfel.la +libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ + src/thread-pool.c src/image.c src/hdf5-file.c \ + src/beam-parameters.c + +INCLUDES = "-I$(top_srcdir)/data" +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + ( cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh ) \ + && { if test -f $@; then exit 0; else break; fi; }; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu libcrystfel/Makefile'; \ + $(am__cd) $(top_srcdir) && \ + $(AUTOMAKE) --gnu libcrystfel/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(am__aclocal_m4_deps): +install-libLTLIBRARIES: $(lib_LTLIBRARIES) + @$(NORMAL_INSTALL) + test -z "$(libdir)" || $(MKDIR_P) "$(DESTDIR)$(libdir)" + @list='$(lib_LTLIBRARIES)'; test -n "$(libdir)" || list=; \ + list2=; for p in $$list; do \ + if test -f $$p; then \ + list2="$$list2 $$p"; \ + else :; fi; \ + done; \ + test -z "$$list2" || { \ + echo " $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=install $(INSTALL) $(INSTALL_STRIP_FLAG) $$list2 '$(DESTDIR)$(libdir)'"; \ + $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=install $(INSTALL) $(INSTALL_STRIP_FLAG) $$list2 "$(DESTDIR)$(libdir)"; \ + } + +uninstall-libLTLIBRARIES: + @$(NORMAL_UNINSTALL) + @list='$(lib_LTLIBRARIES)'; test -n "$(libdir)" || list=; \ + for p in $$list; do \ + $(am__strip_dir) \ + echo " $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=uninstall rm -f '$(DESTDIR)$(libdir)/$$f'"; \ + $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=uninstall rm -f "$(DESTDIR)$(libdir)/$$f"; \ + done + +clean-libLTLIBRARIES: + -test -z "$(lib_LTLIBRARIES)" || rm -f $(lib_LTLIBRARIES) + @list='$(lib_LTLIBRARIES)'; for p in $$list; do \ + dir="`echo $$p | sed -e 's|/[^/]*$$||'`"; \ + test "$$dir" != "$$p" || dir=.; \ + echo "rm -f \"$${dir}/so_locations\""; \ + rm -f "$${dir}/so_locations"; \ + done +src/$(am__dirstamp): + @$(MKDIR_P) src + @: > src/$(am__dirstamp) +src/$(DEPDIR)/$(am__dirstamp): + @$(MKDIR_P) src/$(DEPDIR) + @: > src/$(DEPDIR)/$(am__dirstamp) +src/reflist.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/utils.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/cell.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/detector.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/thread-pool.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/image.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/hdf5-file.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/beam-parameters.lo: src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +libcrystfel.la: $(libcrystfel_la_OBJECTS) $(libcrystfel_la_DEPENDENCIES) + $(AM_V_CCLD)$(LINK) -rpath $(libdir) $(libcrystfel_la_OBJECTS) $(libcrystfel_la_LIBADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + -rm -f src/beam-parameters.$(OBJEXT) + -rm -f src/beam-parameters.lo + -rm -f src/cell.$(OBJEXT) + -rm -f src/cell.lo + -rm -f src/detector.$(OBJEXT) + -rm -f src/detector.lo + -rm -f src/hdf5-file.$(OBJEXT) + -rm -f src/hdf5-file.lo + -rm -f src/image.$(OBJEXT) + -rm -f src/image.lo + -rm -f src/reflist.$(OBJEXT) + -rm -f src/reflist.lo + -rm -f src/thread-pool.$(OBJEXT) + -rm -f src/thread-pool.lo + -rm -f src/utils.$(OBJEXT) + -rm -f src/utils.lo + +distclean-compile: + -rm -f *.tab.c + +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/beam-parameters.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/cell.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/detector.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/hdf5-file.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/image.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/thread-pool.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/utils.Plo@am__quote@ + +.c.o: +@am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\ +@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\ +@am__fastdepCC_TRUE@ $(am__mv) $$depbase.Tpo $$depbase.Po +@am__fastdepCC_FALSE@ $(AM_V_CC) @AM_BACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(COMPILE) -c -o $@ $< + +.c.obj: +@am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.obj$$||'`;\ +@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ `$(CYGPATH_W) '$<'` &&\ +@am__fastdepCC_TRUE@ $(am__mv) $$depbase.Tpo $$depbase.Po +@am__fastdepCC_FALSE@ $(AM_V_CC) @AM_BACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(COMPILE) -c -o $@ `$(CYGPATH_W) '$<'` + +.c.lo: +@am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.lo$$||'`;\ +@am__fastdepCC_TRUE@ $(LTCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\ +@am__fastdepCC_TRUE@ $(am__mv) $$depbase.Tpo $$depbase.Plo +@am__fastdepCC_FALSE@ $(AM_V_CC) @AM_BACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + -rm -rf src/.libs src/_libs + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) '{ files[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in files) print i; }; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + set x; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) '{ files[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in files) print i; }; }'`; \ + shift; \ + if test -z "$(ETAGS_ARGS)$$*$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + if test $$# -gt 0; then \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + "$$@" $$unique; \ + else \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$unique; \ + fi; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) '{ files[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in files) print i; }; }'`; \ + test -z "$(CTAGS_ARGS)$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && $(am__cd) $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) "$$here" + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \ + list='$(DISTFILES)'; \ + dist_files=`for file in $$list; do echo $$file; done | \ + sed -e "s|^$$srcdirstrip/||;t" \ + -e "s|^$$topsrcdirstrip/|$(top_builddir)/|;t"`; \ + case $$dist_files in \ + */*) $(MKDIR_P) `echo "$$dist_files" | \ + sed '/\//!d;s|^|$(distdir)/|;s,/[^/]*$$,,' | \ + sort -u` ;; \ + esac; \ + for file in $$dist_files; do \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + if test -d $$d/$$file; then \ + dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test -d "$(distdir)/$$file"; then \ + find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \ + fi; \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -fpR $(srcdir)/$$file "$(distdir)$$dir" || exit 1; \ + find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \ + fi; \ + cp -fpR $$d/$$file "$(distdir)$$dir" || exit 1; \ + else \ + test -f "$(distdir)/$$file" \ + || cp -p $$d/$$file "$(distdir)/$$file" \ + || exit 1; \ + fi; \ + done +check-am: all-am +check: check-am +all-am: Makefile $(LTLIBRARIES) +installdirs: + for dir in "$(DESTDIR)$(libdir)"; do \ + test -z "$$dir" || $(MKDIR_P) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + -test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES) + -rm -f src/$(DEPDIR)/$(am__dirstamp) + -rm -f src/$(am__dirstamp) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-generic clean-libLTLIBRARIES clean-libtool \ + mostlyclean-am + +distclean: distclean-am + -rm -rf src/$(DEPDIR) + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +html-am: + +info: info-am + +info-am: + +install-data-am: + +install-dvi: install-dvi-am + +install-dvi-am: + +install-exec-am: install-libLTLIBRARIES + +install-html: install-html-am + +install-html-am: + +install-info: install-info-am + +install-info-am: + +install-man: + +install-pdf: install-pdf-am + +install-pdf-am: + +install-ps: install-ps-am + +install-ps-am: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -rf src/$(DEPDIR) + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-libLTLIBRARIES + +.MAKE: install-am install-strip + +.PHONY: CTAGS GTAGS all all-am check check-am clean clean-generic \ + clean-libLTLIBRARIES clean-libtool ctags distclean \ + distclean-compile distclean-generic distclean-libtool \ + distclean-tags distdir dvi dvi-am html html-am info info-am \ + install install-am install-data install-data-am install-dvi \ + install-dvi-am install-exec install-exec-am install-html \ + install-html-am install-info install-info-am \ + install-libLTLIBRARIES install-man install-pdf install-pdf-am \ + install-ps install-ps-am install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ + tags uninstall uninstall-am uninstall-libLTLIBRARIES + + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/libcrystfel/src/beam-parameters.c b/libcrystfel/src/beam-parameters.c new file mode 100644 index 00000000..ac30ad04 --- /dev/null +++ b/libcrystfel/src/beam-parameters.c @@ -0,0 +1,133 @@ +/* + * beam-parameters.c + * + * Beam parameters + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdio.h> +#include <stdlib.h> + +#include "beam-parameters.h" +#include "utils.h" + + +struct beam_params *get_beam_parameters(const char *filename) +{ + FILE *fh; + struct beam_params *b; + char *rval; + int reject; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return NULL; + + b = calloc(1, sizeof(struct beam_params)); + if ( b == NULL ) return NULL; + + b->fluence = -1.0; + b->beam_radius = -1.0; + b->photon_energy = -1.0; + b->bandwidth = -1.0; + b->divergence = -1.0; + b->dqe = -1.0; + b->adu_per_photon = -1.0; + + do { + + int n1; + char line[1024]; + int i; + char **bits; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + chomp(line); + + n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); + if ( n1 < 3 ) { + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + continue; + } + + if ( bits[1][0] != '=' ) { + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + continue; + } + + if ( strcmp(bits[0], "beam/fluence") == 0 ) { + b->fluence = atof(bits[2]); + } else if ( strcmp(bits[0], "beam/radius") == 0 ) { + b->beam_radius = atof(bits[2]); + } else if ( strcmp(bits[0], "beam/photon_energy") == 0 ) { + b->photon_energy = atof(bits[2]); + } else if ( strcmp(bits[0], "beam/bandwidth") == 0 ) { + b->bandwidth = atof(bits[2]); + } else if ( strcmp(bits[0], "beam/divergence") == 0 ) { + b->divergence = atof(bits[2]); + } else if ( strcmp(bits[0], "detector/dqe") == 0 ) { + b->dqe = atof(bits[2]); + } else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) { + b->adu_per_photon = atof(bits[2]); + } else { + ERROR("Unrecognised field '%s'\n", bits[0]); + } + + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + + } while ( rval != NULL ); + fclose(fh); + + reject = 0; + + if ( b->fluence < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/fluence'.\n"); + reject = 1; + } + if ( b->beam_radius < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/radius'.\n"); + reject = 1; + } + if ( b->photon_energy < 0.0 ) { + ERROR("Invalid or unspecified value for" + " 'beam/photon_energy'.\n"); + reject = 1; + } + if ( b->bandwidth < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/bandwidth'.\n"); + reject = 1; + } + if ( b->divergence < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/divergence'.\n"); + reject = 1; + } + if ( b->dqe < 0.0 ) { + ERROR("Invalid or unspecified value for 'detector/dqe'.\n"); + reject = 1; + } + if ( b->adu_per_photon < 0.0 ) { + ERROR("Invalid or unspecified value for" + " 'detector/adu_per_photon'.\n"); + reject = 1; + } + + if ( reject ) { + ERROR("Please fix the above problems with the beam" + " parameters file and try again.\n"); + return NULL; + } + + return b; +} diff --git a/libcrystfel/src/beam-parameters.h b/libcrystfel/src/beam-parameters.h new file mode 100644 index 00000000..411ab449 --- /dev/null +++ b/libcrystfel/src/beam-parameters.h @@ -0,0 +1,40 @@ +/* + * beam-parameters.h + * + * Beam parameters + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef BEAM_PARAMETERS_H +#define BEAM_PARAMETERS_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +struct beam_params +{ + double fluence; /* photons per pulse */ + double beam_radius; /* metres */ + double photon_energy; /* eV per photon */ + double bandwidth; /* FWHM(wavelength) over wavelength. + * Note: current simulation code just uses + * a rectangular distribution with this as + * its (full) width. */ + double divergence; /* divergence (radians) */ + + double dqe; /* Detector DQE (fraction) */ + double adu_per_photon; /* Detector "gain" */ +}; + + +extern struct beam_params *get_beam_parameters(const char *filename); + + +#endif /* BEAM_PARAMETERS_H */ diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c new file mode 100644 index 00000000..c31697bc --- /dev/null +++ b/libcrystfel/src/cell.c @@ -0,0 +1,1164 @@ +/* + * cell.c + * + * Unit Cell Calculations + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <math.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> + +#include "cell.h" +#include "utils.h" +#include "image.h" + + +/** + * SECTION:unitcell + * @short_description: Unit cell + * @title: UnitCell + * @section_id: + * @see_also: + * @include: "cell.h" + * @Image: + * + * This structure represents a unit cell. + */ + + + +/* Weighting factor of lengths relative to angles */ +#define LWEIGHT (10.0e-9) + +typedef enum { + CELL_REP_CRYST, + CELL_REP_CART, + CELL_REP_RECIP +} CellRepresentation; + +struct _unitcell { + + CellRepresentation rep; + + /* Crystallographic representation */ + double a; /* m */ + double b; /* m */ + double c; /* m */ + double alpha; /* Radians */ + double beta; /* Radians */ + double gamma; /* Radians */ + + /* Cartesian representation */ + double ax; double bx; double cx; + double ay; double by; double cy; + double az; double bz; double cz; + + /* Cartesian representation of reciprocal axes */ + double axs; double bxs; double cxs; + double ays; double bys; double cys; + double azs; double bzs; double czs; + + char *pointgroup; + char *spacegroup; +}; + + +/************************** Setters and Constructors **************************/ + + +/** + * cell_new: + * + * Create a new %UnitCell. + * + * Returns: the new unit cell, or NULL on failure. + * + */ +UnitCell *cell_new() +{ + UnitCell *cell; + + cell = malloc(sizeof(UnitCell)); + if ( cell == NULL ) return NULL; + + cell->a = 1.0; + cell->b = 1.0; + cell->c = 1.0; + cell->alpha = M_PI_2; + cell->beta = M_PI_2; + cell->gamma = M_PI_2; + + cell->rep = CELL_REP_CRYST; + + cell->pointgroup = strdup("1"); + cell->spacegroup = strdup("P 1"); + + return cell; +} + + +/** + * cell_free: + * @cell: A %UnitCell to free. + * + * Frees a %UnitCell, and all internal resources concerning that cell. + * + */ +void cell_free(UnitCell *cell) +{ + if ( cell == NULL ) return; + free(cell->pointgroup); + free(cell->spacegroup); + free(cell); +} + + +void cell_set_parameters(UnitCell *cell, double a, double b, double c, + double alpha, double beta, double gamma) +{ + if ( cell == NULL ) return; + + cell->a = a; + cell->b = b; + cell->c = c; + cell->alpha = alpha; + cell->beta = beta; + cell->gamma = gamma; + + cell->rep = CELL_REP_CRYST; +} + + +void cell_set_cartesian(UnitCell *cell, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz) +{ + if ( cell == NULL ) return; + + cell->ax = ax; cell->ay = ay; cell->az = az; + cell->bx = bx; cell->by = by; cell->bz = bz; + cell->cx = cx; cell->cy = cy; cell->cz = cz; + + cell->rep = CELL_REP_CART; +} + + +void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az) +{ + if ( cell == NULL ) return; + cell->ax = ax; cell->ay = ay; cell->az = az; + cell->rep = CELL_REP_CART; +} + + +void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz) +{ + if ( cell == NULL ) return; + cell->bx = bx; cell->by = by; cell->bz = bz; + cell->rep = CELL_REP_CART; +} + + +void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz) +{ + if ( cell == NULL ) return; + cell->cx = cx; cell->cy = cy; cell->cz = cz; + cell->rep = CELL_REP_CART; +} + + +UnitCell *cell_new_from_parameters(double a, double b, double c, + double alpha, double beta, double gamma) +{ + UnitCell *cell; + + cell = cell_new(); + if ( cell == NULL ) return NULL; + + cell_set_parameters(cell, a, b, c, alpha, beta, gamma); + + return cell; +} + + +UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, + struct rvec cs) +{ + UnitCell *cell; + + cell = cell_new(); + if ( cell == NULL ) return NULL; + + cell->axs = as.u; cell->ays = as.v; cell->azs = as.w; + cell->bxs = bs.u; cell->bys = bs.v; cell->bzs = bs.w; + cell->cxs = cs.u; cell->cys = cs.v; cell->czs = cs.w; + + cell->rep = CELL_REP_RECIP; + + return cell; +} + + +UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c) +{ + UnitCell *cell; + + cell = cell_new(); + if ( cell == NULL ) return NULL; + + cell->ax = a.u; cell->ay = a.v; cell->az = a.w; + cell->bx = b.u; cell->by = b.v; cell->bz = b.w; + cell->cx = c.u; cell->cy = c.v; cell->cz = c.w; + + cell->rep = CELL_REP_CART; + + return cell; +} + + +UnitCell *cell_new_from_cell(UnitCell *orig) +{ + UnitCell *new; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + new = cell_new(); + + cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz); + cell_set_pointgroup(new, orig->pointgroup); + cell_set_spacegroup(new, orig->spacegroup); + + return new; +} + + +void cell_set_reciprocal(UnitCell *cell, + double asx, double asy, double asz, + double bsx, double bsy, double bsz, + double csx, double csy, double csz) +{ + if ( cell == NULL ) return; + + cell->axs = asx; cell->ays = asy; cell->azs = asz; + cell->bxs = bsx; cell->bys = bsy; cell->bzs = bsz; + cell->cxs = csx; cell->cys = csy; cell->czs = csz; + + cell->rep = CELL_REP_RECIP; +} + + +void cell_set_spacegroup(UnitCell *cell, const char *sym) +{ + free(cell->spacegroup); + cell->spacegroup = strdup(sym); +} + + +void cell_set_pointgroup(UnitCell *cell, const char *sym) +{ + free(cell->pointgroup); + cell->pointgroup = strdup(sym); +} + + +/************************* Getter helper functions ****************************/ + +static int cell_crystallographic_to_cartesian(UnitCell *cell, + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz) +{ + double tmp, V, cosalphastar, cstar; + + /* Firstly: Get a in terms of x, y and z + * +a (cryst) is defined to lie along +x (cart) */ + *ax = cell->a; + *ay = 0.0; + *az = 0.0; + + /* b in terms of x, y and z + * b (cryst) is defined to lie in the xy (cart) plane */ + *bx = cell->b*cos(cell->gamma); + *by = cell->b*sin(cell->gamma); + *bz = 0.0; + + tmp = cos(cell->alpha)*cos(cell->alpha) + + cos(cell->beta)*cos(cell->beta) + + cos(cell->gamma)*cos(cell->gamma) + - 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma); + V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp); + + cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha); + cosalphastar /= sin(cell->beta)*sin(cell->gamma); + + cstar = (cell->a * cell->b * sin(cell->gamma))/V; + + /* c in terms of x, y and z */ + *cx = cell->c*cos(cell->beta); + *cy = -cell->c*sin(cell->beta)*cosalphastar; + *cz = 1.0/cstar; + + return 0; +} + + +/* Why yes, I do enjoy long argument lists...! */ +static int cell_invert(double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz) +{ + int s; + gsl_matrix *m; + gsl_matrix *inv; + gsl_permutation *perm; + + m = gsl_matrix_alloc(3, 3); + if ( m == NULL ) { + ERROR("Couldn't allocate memory for matrix\n"); + return 1; + } + gsl_matrix_set(m, 0, 0, ax); + gsl_matrix_set(m, 0, 1, bx); + gsl_matrix_set(m, 0, 2, cx); + gsl_matrix_set(m, 1, 0, ay); + gsl_matrix_set(m, 1, 1, by); + gsl_matrix_set(m, 1, 2, cy); + gsl_matrix_set(m, 2, 0, az); + gsl_matrix_set(m, 2, 1, bz); + gsl_matrix_set(m, 2, 2, cz); + + /* Invert */ + perm = gsl_permutation_alloc(m->size1); + if ( perm == NULL ) { + ERROR("Couldn't allocate permutation\n"); + gsl_matrix_free(m); + return 1; + } + inv = gsl_matrix_alloc(m->size1, m->size2); + if ( inv == NULL ) { + ERROR("Couldn't allocate inverse\n"); + gsl_matrix_free(m); + gsl_permutation_free(perm); + return 1; + } + if ( gsl_linalg_LU_decomp(m, perm, &s) ) { + ERROR("Couldn't decompose matrix\n"); + gsl_matrix_free(m); + gsl_permutation_free(perm); + return 1; + } + if ( gsl_linalg_LU_invert(m, perm, inv) ) { + ERROR("Couldn't invert matrix\n"); + gsl_matrix_free(m); + gsl_permutation_free(perm); + return 1; + } + gsl_permutation_free(perm); + gsl_matrix_free(m); + + /* Transpose */ + gsl_matrix_transpose(inv); + + *asx = gsl_matrix_get(inv, 0, 0); + *bsx = gsl_matrix_get(inv, 0, 1); + *csx = gsl_matrix_get(inv, 0, 2); + *asy = gsl_matrix_get(inv, 1, 0); + *bsy = gsl_matrix_get(inv, 1, 1); + *csy = gsl_matrix_get(inv, 1, 2); + *asz = gsl_matrix_get(inv, 2, 0); + *bsz = gsl_matrix_get(inv, 2, 1); + *csz = gsl_matrix_get(inv, 2, 2); + + gsl_matrix_free(inv); + + return 0; +} + + +/********************************** Getters ***********************************/ + +int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, + double *alpha, double *beta, double *gamma) +{ + double ax, ay, az, bx, by, bz, cx, cy, cz; + + if ( cell == NULL ) return 1; + + switch ( cell->rep ) { + + case CELL_REP_CRYST: + /* Direct response */ + *a = cell->a; + *b = cell->b; + *c = cell->c; + *alpha = cell->alpha; + *beta = cell->beta; + *gamma = cell->gamma; + return 0; + + case CELL_REP_CART: + /* Convert cartesian -> crystallographic */ + *a = modulus(cell->ax, cell->ay, cell->az); + *b = modulus(cell->bx, cell->by, cell->bz); + *c = modulus(cell->cx, cell->cy, cell->cz); + + *alpha = angle_between(cell->bx, cell->by, cell->bz, + cell->cx, cell->cy, cell->cz); + *beta = angle_between(cell->ax, cell->ay, cell->az, + cell->cx, cell->cy, cell->cz); + *gamma = angle_between(cell->ax, cell->ay, cell->az, + cell->bx, cell->by, cell->bz); + return 0; + + case CELL_REP_RECIP: + /* Convert reciprocal -> crystallographic. + * Start by converting reciprocal -> cartesian */ + cell_invert(cell->axs, cell->ays, cell->azs, + cell->bxs, cell->bys, cell->bzs, + cell->cxs, cell->cys, cell->czs, + &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + /* Now convert cartesian -> crystallographic */ + *a = modulus(ax, ay, az); + *b = modulus(bx, by, bz); + *c = modulus(cx, cy, cz); + + *alpha = angle_between(bx, by, bz, cx, cy, cz); + *beta = angle_between(ax, ay, az, cx, cy, cz); + *gamma = angle_between(ax, ay, az, bx, by, bz); + return 0; + } + + return 1; +} + + +int cell_get_cartesian(UnitCell *cell, + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz) +{ + if ( cell == NULL ) return 1; + + switch ( cell->rep ) { + + case CELL_REP_CRYST: + /* Convert crystallographic -> cartesian. */ + return cell_crystallographic_to_cartesian(cell, + ax, ay, az, + bx, by, bz, + cx, cy, cz); + + case CELL_REP_CART: + /* Direct response */ + *ax = cell->ax; *ay = cell->ay; *az = cell->az; + *bx = cell->bx; *by = cell->by; *bz = cell->bz; + *cx = cell->cx; *cy = cell->cy; *cz = cell->cz; + return 0; + + case CELL_REP_RECIP: + /* Convert reciprocal -> cartesian */ + return cell_invert(cell->axs, cell->ays, cell->azs, + cell->bxs, cell->bys, cell->bzs, + cell->cxs, cell->cys, cell->czs, + ax, ay, az, bx, by, bz, cx, cy, cz); + + } + + return 1; +} + + +int cell_get_reciprocal(UnitCell *cell, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz) +{ + int r; + double ax, ay, az, bx, by, bz, cx, cy, cz; + if ( cell == NULL ) return 1; + + switch ( cell->rep ) { + + case CELL_REP_CRYST: + /* Convert crystallographic -> reciprocal */ + r = cell_crystallographic_to_cartesian(cell, + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + if ( r ) return r; + return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz, + asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + case CELL_REP_CART: + /* Convert cartesian -> reciprocal */ + cell_invert(cell->ax, cell->ay, cell->az, + cell->bx, cell->by, cell->bz, + cell->cx, cell->cy, cell->cz, + asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + return 0; + + case CELL_REP_RECIP: + /* Direct response */ + *asx = cell->axs; *asy = cell->ays; *asz = cell->azs; + *bsx = cell->bxs; *bsy = cell->bys; *bsz = cell->bzs; + *csx = cell->cxs; *csy = cell->cys; *csz = cell->czs; + return 0; + + } + + return 1; +} + + +const char *cell_get_pointgroup(UnitCell *cell) +{ + return cell->pointgroup; +} + + +const char *cell_get_spacegroup(UnitCell *cell) +{ + return cell->spacegroup; +} + + + + + +/********************************* Utilities **********************************/ + +static const char *cell_rep(UnitCell *cell) +{ + switch ( cell->rep ) { + case CELL_REP_CRYST: + return "crystallographic, direct space"; + case CELL_REP_CART: + return "cartesian, direct space"; + case CELL_REP_RECIP: + return "cartesian, reciprocal space"; + } + + return "unknown"; +} + + +/** + * cell_rotate: + * @in: A %UnitCell to rotate + * @quat: A %quaternion + * + * Rotate a %UnitCell using a %quaternion. + * + * Returns: a newly allocated rotated copy of @in. + * + */ +UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) +{ + struct rvec a, b, c; + struct rvec an, bn, cn; + UnitCell *out = cell_new_from_cell(in); + + cell_get_cartesian(in, &a.u, &a.v, &a.w, + &b.u, &b.v, &b.w, + &c.u, &c.v, &c.w); + + an = quat_rot(a, quat); + bn = quat_rot(b, quat); + cn = quat_rot(c, quat); + + cell_set_cartesian(out, an.u, an.v, an.w, + bn.u, bn.v, bn.w, + cn.u, cn.v, cn.w); + + return out; +} + + +void cell_print(UnitCell *cell) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double a, b, c, alpha, beta, gamma; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + STATUS(" a b c alpha beta gamma\n"); + STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n", + a*1e9, b*1e9, c*1e9, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az); + STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz); + STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz); + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + asx, asy, asz, modulus(asx, asy, asz)); + STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + bsx, bsy, bsz, modulus(bsx, bsy, bsz)); + STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + csx, csy, csz, modulus(csx, csy, csz)); + + STATUS("Point group: %s\n", cell_get_pointgroup(cell)); + STATUS("Cell representation is %s.\n", cell_rep(cell)); +} + + +#define MAX_CAND (1024) + +static int right_handed(struct rvec a, struct rvec b, struct rvec c) +{ + struct rvec aCb; + double aCb_dot_c; + + /* "a" cross "b" */ + aCb.u = a.v*b.w - a.w*b.v; + aCb.v = - (a.u*b.w - a.w*b.u); + aCb.w = a.u*b.v - a.v*b.u; + + /* "a cross b" dot "c" */ + aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; + + if ( aCb_dot_c > 0.0 ) return 1; + return 0; +} + + +struct cvec { + struct rvec vec; + float na; + float nb; + float nc; + float fom; +}; + + +static int same_vector(struct cvec a, struct cvec b) +{ + if ( a.na != b.na ) return 0; + if ( a.nb != b.nb ) return 0; + if ( a.nc != b.nc ) return 0; + return 1; +} + + +/* Attempt to make 'cell' fit into 'template' somehow */ +UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, + int reduce) +{ + signed int n1l, n2l, n3l; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + int i, j; + double lengths[3]; + double angles[3]; + struct cvec *cand[3]; + UnitCell *new_cell = NULL; + float best_fom = +999999999.9; /* Large number.. */ + int ncand[3] = {0,0,0}; + signed int ilow, ihigh; + float ltl = 5.0; /* percent */ + float angtol = deg2rad(1.5); + + if ( cell_get_reciprocal(template, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell for template.\n"); + return NULL; + } + + lengths[0] = modulus(asx, asy, asz); + lengths[1] = modulus(bsx, bsy, bsz); + lengths[2] = modulus(csx, csy, csz); + + angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); + angles[1] = angle_between(asx, asy, asz, csx, csy, csz); + angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); + + cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); + cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); + cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); + + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell.\n"); + return NULL; + } + + if ( reduce ) { + ilow = -2; ihigh = 4; + } else { + ilow = 0; ihigh = 1; + } + + /* Negative values mean 1/n, positive means n, zero means zero */ + for ( n1l=ilow; n1l<=ihigh; n1l++ ) { + for ( n2l=ilow; n2l<=ihigh; n2l++ ) { + for ( n3l=ilow; n3l<=ihigh; n3l++ ) { + + float n1, n2, n3; + signed int b1, b2, b3; + + n1 = (n1l>=0) ? (n1l) : (1.0/n1l); + n2 = (n2l>=0) ? (n2l) : (1.0/n2l); + n3 = (n3l>=0) ? (n3l) : (1.0/n3l); + + if ( !reduce ) { + if ( n1l + n2l + n3l > 1 ) continue; + } + + /* 'bit' values can be +1 or -1 */ + for ( b1=-1; b1<=1; b1+=2 ) { + for ( b2=-1; b2<=1; b2+=2 ) { + for ( b3=-1; b3<=1; b3+=2 ) { + + double tx, ty, tz; + double tlen; + int i; + + n1 *= b1; n2 *= b2; n3 *= b3; + + tx = n1*asx + n2*bsx + n3*csx; + ty = n1*asy + n2*bsy + n3*csy; + tz = n1*asz + n2*bsz + n3*csz; + tlen = modulus(tx, ty, tz); + + /* Test modulus for agreement with moduli of template */ + for ( i=0; i<3; i++ ) { + + if ( !within_tolerance(lengths[i], tlen, ltl) ) + continue; + + cand[i][ncand[i]].vec.u = tx; + cand[i][ncand[i]].vec.v = ty; + cand[i][ncand[i]].vec.w = tz; + cand[i][ncand[i]].na = n1; + cand[i][ncand[i]].nb = n2; + cand[i][ncand[i]].nc = n3; + cand[i][ncand[i]].fom = fabs(lengths[i] - tlen); + if ( ncand[i] == MAX_CAND ) { + ERROR("Too many candidates\n"); + } else { + ncand[i]++; + } + } + + } + } + } + + } + } + } + + if ( verbose ) { + STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); + } + + for ( i=0; i<ncand[0]; i++ ) { + for ( j=0; j<ncand[1]; j++ ) { + + double ang; + int k; + float fom1; + + if ( same_vector(cand[0][i], cand[1][j]) ) continue; + + /* Measure the angle between the ith candidate for axis 0 + * and the jth candidate for axis 1 */ + ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v, + cand[0][i].vec.w, cand[1][j].vec.u, + cand[1][j].vec.v, cand[1][j].vec.w); + + /* Angle between axes 0 and 1 should be angle 2 */ + if ( fabs(ang - angles[2]) > angtol ) continue; + + fom1 = fabs(ang - angles[2]); + + for ( k=0; k<ncand[2]; k++ ) { + + float fom2, fom3; + + if ( same_vector(cand[1][j], cand[2][k]) ) continue; + + /* Measure the angle between the current candidate for + * axis 0 and the kth candidate for axis 2 */ + ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v, + cand[0][i].vec.w, cand[2][k].vec.u, + cand[2][k].vec.v, cand[2][k].vec.w); + + /* ... it should be angle 1 ... */ + if ( fabs(ang - angles[1]) > angtol ) continue; + + fom2 = fom1 + fabs(ang - angles[1]); + + /* Finally, the angle between the current candidate for + * axis 1 and the kth candidate for axis 2 */ + ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, + cand[1][j].vec.w, cand[2][k].vec.u, + cand[2][k].vec.v, cand[2][k].vec.w); + + /* ... it should be angle 0 ... */ + if ( fabs(ang - angles[0]) > angtol ) continue; + + /* Unit cell must be right-handed */ + if ( !right_handed(cand[0][i].vec, cand[1][j].vec, + cand[2][k].vec) ) continue; + + fom3 = fom2 + fabs(ang - angles[0]); + fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom + + cand[2][k].fom); + + if ( fom3 < best_fom ) { + if ( new_cell != NULL ) free(new_cell); + new_cell = cell_new_from_reciprocal_axes( + cand[0][i].vec, cand[1][j].vec, + cand[2][k].vec); + best_fom = fom3; + } + + } + + } + } + + free(cand[0]); + free(cand[1]); + free(cand[2]); + + return new_cell; +} + + + +UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + int i; + double lengths[3]; + int used[3]; + struct rvec real_a, real_b, real_c; + struct rvec params[3]; + double alen, blen, clen; + float ltl = 5.0; /* percent */ + int have_real_a; + int have_real_b; + int have_real_c; + + /* Get the lengths to match */ + if ( cell_get_cartesian(template, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) + { + ERROR("Couldn't get cell for template.\n"); + return NULL; + } + alen = modulus(ax, ay, az); + blen = modulus(bx, by, bz); + clen = modulus(cx, cy, cz); + + /* Get the lengths from the cell and turn them into anonymous vectors */ + if ( cell_get_cartesian(cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) + { + ERROR("Couldn't get cell.\n"); + return NULL; + } + lengths[0] = modulus(ax, ay, az); + lengths[1] = modulus(bx, by, bz); + lengths[2] = modulus(cx, cy, cz); + used[0] = 0; used[1] = 0; used[2] = 0; + params[0].u = ax; params[0].v = ay; params[0].w = az; + params[1].u = bx; params[1].v = by; params[1].w = bz; + params[2].u = cx; params[2].v = cy; params[2].w = cz; + + real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; + real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; + real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; + + /* Check each vector against a and b */ + have_real_a = 0; + have_real_b = 0; + for ( i=0; i<3; i++ ) { + if ( within_tolerance(lengths[i], alen, ltl) + && !used[i] && !have_real_a ) + { + used[i] = 1; + memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); + have_real_a = 1; + } + if ( within_tolerance(lengths[i], blen, ltl) + && !used[i] && !have_real_b ) + { + used[i] = 1; + memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); + have_real_b = 1; + } + } + + /* Have we matched both a and b? */ + if ( !(have_real_a && have_real_b) ) return NULL; + + /* "c" is "the other one" */ + have_real_c = 0; + for ( i=0; i<3; i++ ) { + if ( !used[i] ) { + memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); + have_real_c = 1; + } + } + + if ( !have_real_c ) { + ERROR("Huh? Couldn't find the third vector.\n"); + ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); + return NULL; + } + + /* Flip c if not right-handed */ + if ( !right_handed(real_a, real_b, real_c) ) { + real_c.u = -real_c.u; + real_c.v = -real_c.v; + real_c.w = -real_c.w; + } + + return cell_new_from_direct_axes(real_a, real_b, real_c); +} + + +/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ +double resolution(UnitCell *cell, signed int h, signed int k, signed int l) +{ + double a, b, c, alpha, beta, gamma; + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha) + - cos(beta)*cos(beta) + - cos(gamma)*cos(gamma) + + 2*cos(alpha)*cos(beta)*cos(gamma) ); + + const double S11 = b*b*c*c*sin(alpha)*sin(alpha); + const double S22 = a*a*c*c*sin(beta)*sin(beta); + const double S33 = a*a*b*b*sin(gamma)*sin(gamma); + const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma)); + const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha)); + const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta)); + + const double brackets = S11*h*h + S22*k*k + S33*l*l + + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l; + const double oneoverdsq = brackets / Vsq; + const double oneoverd = sqrt(oneoverdsq); + + return oneoverd / 2; +} + + +static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym) +{ + char *new = NULL; + + if ( strcmp(sym, "P 1") == 0 ) new = "1"; + if ( strcmp(sym, "P 63") == 0 ) new = "6"; + if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222"; + if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222"; + if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422"; + + if ( new != NULL ) { + cell_set_pointgroup(cell, new); + } else { + ERROR("Can't determine point group for '%s'\n", sym); + } +} + + +UnitCell *load_cell_from_pdb(const char *filename) +{ + FILE *fh; + char *rval; + UnitCell *cell = NULL; + + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Couldn't open '%s'\n", filename); + return NULL; + } + + do { + + char line[1024]; + + rval = fgets(line, 1023, fh); + + if ( strncmp(line, "CRYST1", 6) == 0 ) { + + float a, b, c, al, be, ga; + char as[10], bs[10], cs[10]; + char als[8], bes[8], gas[8]; + char *sym; + int r; + + memcpy(as, line+6, 9); as[9] = '\0'; + memcpy(bs, line+15, 9); bs[9] = '\0'; + memcpy(cs, line+24, 9); cs[9] = '\0'; + memcpy(als, line+33, 7); als[7] = '\0'; + memcpy(bes, line+40, 7); bes[7] = '\0'; + memcpy(gas, line+47, 7); gas[7] = '\0'; + + r = sscanf(as, "%f", &a); + r += sscanf(bs, "%f", &b); + r += sscanf(cs, "%f", &c); + r += sscanf(als, "%f", &al); + r += sscanf(bes, "%f", &be); + r += sscanf(gas, "%f", &ga); + + if ( r != 6 ) { + STATUS("Couldn't understand CRYST1 line.\n"); + continue; + } + + cell = cell_new_from_parameters(a*1e-10, + b*1e-10, c*1e-10, + deg2rad(al), + deg2rad(be), + deg2rad(ga)); + + if ( strlen(line) > 65 ) { + sym = strndup(line+55, 10); + notrail(sym); + cell_set_pointgroup_from_pdb(cell, sym); + cell_set_spacegroup(cell, sym); + free(sym); + } else { + cell_set_pointgroup_from_pdb(cell, "P 1"); + cell_set_spacegroup(cell, "P 1"); + ERROR("CRYST1 line without space group.\n"); + } + + break; /* Done */ + } + + } while ( rval != NULL ); + + fclose(fh); + + return cell; +} + + +#ifdef GSL_FUDGE +/* Force the linker to bring in CBLAS to make GSL happy */ +void cell_fudge_gslcblas() +{ + STATUS("%p\n", cblas_sgemm); +} +#endif + + +UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot) +{ + UnitCell *out; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double xnew, ynew, znew; + + cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, + &bsz, &csx, &csy, &csz); + + /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */ + xnew = asx*cos(omega) + asy*sin(omega); + ynew = -asx*sin(omega) + asy*cos(omega); + znew = asz; + asx = xnew; asy = ynew; asz = znew; + xnew = bsx*cos(omega) + bsy*sin(omega); + ynew = -bsx*sin(omega) + bsy*cos(omega); + znew = bsz; + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx*cos(omega) + csy*sin(omega); + ynew = -csx*sin(omega) + csy*cos(omega); + znew = csz; + csx = xnew; csy = ynew; csz = znew; + + /* Rotate by "phi" about +x (not parallel to anything specific) */ + xnew = asx; + ynew = asy*cos(phi) + asz*sin(phi); + znew = -asy*sin(phi) + asz*cos(phi); + asx = xnew; asy = ynew; asz = znew; + xnew = bsx; + ynew = bsy*cos(phi) + bsz*sin(phi); + znew = -bsy*sin(phi) + bsz*cos(phi); + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx; + ynew = csy*cos(phi) + csz*sin(phi); + znew = -csy*sin(phi) + csz*cos(phi); + csx = xnew; csy = ynew; csz = znew; + + /* Rotate by "rot" about the new +z (in-plane rotation) */ + xnew = asx*cos(rot) + asy*sin(rot); + ynew = -asx*sin(rot) + asy*cos(rot); + znew = asz; + asx = xnew; asy = ynew; asz = znew; + xnew = bsx*cos(rot) + bsy*sin(rot); + ynew = -bsx*sin(rot) + bsy*cos(rot); + znew = bsz; + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx*cos(rot) + csy*sin(rot); + ynew = -csx*sin(rot) + csy*cos(rot); + znew = csz; + csx = xnew; csy = ynew; csz = znew; + + out = cell_new_from_cell(in); + cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + return out; +} + + +int cell_is_sensible(UnitCell *cell) +{ + double a, b, c, al, be, ga; + + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + if ( al + be + ga >= 2.0*M_PI ) return 0; + if ( al + be - ga >= 2.0*M_PI ) return 0; + if ( al - be + ga >= 2.0*M_PI ) return 0; + if ( - al + be + ga >= 2.0*M_PI ) return 0; + if ( al + be + ga <= 0.0 ) return 0; + if ( al + be - ga <= 0.0 ) return 0; + if ( al - be + ga <= 0.0 ) return 0; + if ( - al + be + ga <= 0.0 ) return 0; + if ( isnan(al) ) return 0; + if ( isnan(be) ) return 0; + if ( isnan(ga) ) return 0; + return 1; +} diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h new file mode 100644 index 00000000..b5d31fc6 --- /dev/null +++ b/libcrystfel/src/cell.h @@ -0,0 +1,107 @@ +/* + * cell.h + * + * Unit Cell Calculations + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef CELL_H +#define CELL_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "utils.h" + +/* A 3D vector in reciprocal space */ +struct rvec +{ + double u; + double v; + double w; +}; + + +/** + * UnitCell: + * + * This data structure is opaque. You must use the available accessor functions + * to read and write its contents. + **/ +typedef struct _unitcell UnitCell; + +extern UnitCell *cell_new(void); +extern UnitCell *cell_new_from_cell(UnitCell *orig); +extern void cell_free(UnitCell *cell); + +/* Lengths in m, angles in radians */ +extern UnitCell *cell_new_from_parameters(double a, double b, double c, + double alpha, double beta, double gamma); + +extern UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, + struct rvec cs); + +extern UnitCell *cell_new_from_direct_axes(struct rvec as, struct rvec bs, + struct rvec cs); + +extern void cell_set_cartesian(UnitCell *cell, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz); + +extern void cell_set_parameters(UnitCell *cell, double a, double b, double c, + double alpha, double beta, double gamma); + +extern void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az); +extern void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz); +extern void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz); +extern void cell_set_spacegroup(UnitCell *cell, const char *sym); +extern void cell_set_pointgroup(UnitCell *cell, const char *sym); + + +extern int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, + double *alpha, double *beta, double *gamma); + +extern int cell_get_cartesian(UnitCell *cell, + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz); + +extern int cell_get_reciprocal(UnitCell *cell, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz); + +extern void cell_set_reciprocal(UnitCell *cell, + double asx, double asy, double asz, + double bsx, double bsy, double bsz, + double csx, double csy, double csz); + +extern const char *cell_get_pointgroup(UnitCell *cell); + +extern const char *cell_get_spacegroup(UnitCell *cell); + +extern double resolution(UnitCell *cell, + signed int h, signed int k, signed int l); + +extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat); +extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, + double rot); + +extern void cell_print(UnitCell *cell); + +extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, + int reduce); + +extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template); + +extern UnitCell *load_cell_from_pdb(const char *filename); + +extern int cell_is_sensible(UnitCell *cell); + +#endif /* CELL_H */ diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c new file mode 100644 index 00000000..54bae1d7 --- /dev/null +++ b/libcrystfel/src/detector.c @@ -0,0 +1,1233 @@ +/* + * detector.c + * + * Detector properties + * + * (c) 2006-2011 Thomas White <taw@physics.org> + * (c) 2011 Rick Kirian <rkirian@asu.edu> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include <stdlib.h> +#include <math.h> +#include <stdio.h> +#include <string.h> +#include <assert.h> +#include <ctype.h> + +#include "image.h" +#include "utils.h" +#include "detector.h" +#include "beam-parameters.h" +#include "hdf5-file.h" + + +/** + * SECTION:detector + * @short_description: Detector geometry + * @title: Detector + * @section_id: + * @see_also: + * @include: "detector.h" + * @Image: + * + * This structure represents the detector geometry + */ + + +static int atob(const char *a) +{ + if ( strcasecmp(a, "true") == 0 ) return 1; + if ( strcasecmp(a, "false") == 0 ) return 0; + return atoi(a); +} + + +static int assplode_algebraic(const char *a_orig, char ***pbits) +{ + int len, i; + int nexp; + char **bits; + char *a; + int idx, istr; + + len = strlen(a_orig); + + /* Add plus at start if no sign there already */ + if ( (a_orig[0] != '+') && (a_orig[0] != '-') ) { + len += 1; + a = malloc(len+1); + snprintf(a, len+1, "+%s", a_orig); + a[len] = '\0'; + + } else { + a = strdup(a_orig); + } + + /* Count the expressions */ + nexp = 0; + for ( i=0; i<len; i++ ) { + if ( (a[i] == '+') || (a[i] == '-') ) nexp++; + } + + bits = calloc(nexp, sizeof(char *)); + + /* Break the string up */ + idx = -1; + istr = 0; + assert((a[0] == '+') || (a[0] == '-')); + for ( i=0; i<len; i++ ) { + + char ch; + + ch = a[i]; + + if ( (ch == '+') || (ch == '-') ) { + if ( idx >= 0 ) bits[idx][istr] = '\0'; + idx++; + bits[idx] = malloc(len+1); + istr = 0; + } + + if ( !isdigit(ch) && (ch != '.') && (ch != 'x') && (ch != 'y') + && (ch != '+') && (ch != '-') ) + { + ERROR("Invalid character '%C' found.\n", ch); + return 0; + } + + assert(idx >= 0); + bits[idx][istr++] = ch; + + } + if ( idx >= 0 ) bits[idx][istr] = '\0'; + + *pbits = bits; + free(a); + + return nexp; +} + + +/* Parses the scan directions (accounting for possible rotation) + * Assumes all white spaces have been already removed */ +static int dir_conv(const char *a, double *sx, double *sy) +{ + int n; + char **bits; + int i; + + *sx = 0.0; *sy = 0.0; + + n = assplode_algebraic(a, &bits); + + if ( n == 0 ) { + ERROR("Invalid direction '%s'\n", a); + return 1; + } + + for ( i=0; i<n; i++ ) { + + int len; + double val; + char axis; + int j; + + len = strlen(bits[i]); + assert(len != 0); + axis = bits[i][len-1]; + if ( (axis != 'x') && (axis != 'y') ) { + ERROR("Invalid symbol '%C' - must be x or y.\n", axis); + return 1; + } + + /* Chop off the symbol now it's dealt with */ + bits[i][len-1] = '\0'; + + /* Check for anything that isn't part of a number */ + for ( j=0; j<strlen(bits[i]); j++ ) { + if ( isdigit(bits[i][j]) ) continue; + if ( bits[i][j] == '+' ) continue; + if ( bits[i][j] == '-' ) continue; + if ( bits[i][j] == '.' ) continue; + ERROR("Invalid coefficient '%s'\n", bits[i]); + } + + if ( strlen(bits[i]) == 0 ) { + val = 1.0; + } else { + val = atof(bits[i]); + } + if ( strlen(bits[i]) == 1 ) { + if ( bits[i][0] == '+' ) val = 1.0; + if ( bits[i][0] == '-' ) val = -1.0; + } + if ( axis == 'x' ) { + *sx += val; + } else if ( axis == 'y' ) { + *sy += val; + } + + free(bits[i]); + + } + free(bits); + + //STATUS("'%s' -> %5.2fx + %5.2fy\n", a, *sx, *sy); + + return 0; +} + + +struct rvec get_q_for_panel(struct panel *p, double fs, double ss, + double *ttp, double k) +{ + struct rvec q; + double twotheta, r, az; + double rx, ry; + double xs, ys; + + /* Convert xs and ys, which are in fast scan/slow scan coordinates, + * to x and y */ + xs = fs*p->fsx + ss*p->ssx; + ys = fs*p->fsy + ss*p->ssy; + + rx = (xs + p->cnx) / p->res; + ry = (ys + p->cny) / p->res; + + /* Calculate q-vector for this sub-pixel */ + r = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + twotheta = atan2(r, p->clen); + az = atan2(ry, rx); + if ( ttp != NULL ) *ttp = twotheta; + + q.u = k * sin(twotheta)*cos(az); + q.v = k * sin(twotheta)*sin(az); + q.w = k * (cos(twotheta) - 1.0); + + return q; +} + + +struct rvec get_q(struct image *image, double fs, double ss, + double *ttp, double k) +{ + struct panel *p; + const unsigned int fsi = fs; + const unsigned int ssi = ss; /* Explicit rounding */ + + /* Determine which panel to use */ + p = find_panel(image->det, fsi, ssi); + assert(p != NULL); + + return get_q_for_panel(p, fs-(double)p->min_fs, ss-(double)p->min_ss, + ttp, k); +} + + +int in_bad_region(struct detector *det, double fs, double ss) +{ + double rx, ry; + struct panel *p; + double xs, ys; + int i; + + /* Determine which panel to use */ + p = find_panel(det, fs, ss); + + /* No panel found -> definitely bad! */ + if ( p == NULL ) return 1; + + /* Convert xs and ys, which are in fast scan/slow scan coordinates, + * to x and y */ + xs = (fs-(double)p->min_fs)*p->fsx + (ss-(double)p->min_ss)*p->ssx; + ys = (fs-(double)p->min_fs)*p->fsy + (ss-(double)p->min_ss)*p->ssy; + + rx = xs + p->cnx; + ry = ys + p->cny; + + for ( i=0; i<det->n_bad; i++ ) { + struct badregion *b = &det->bad[i]; + if ( rx < b->min_x ) continue; + if ( rx > b->max_x ) continue; + if ( ry < b->min_y ) continue; + if ( ry > b->max_y ) continue; + return 1; + } + + return 0; +} + + +double get_tt(struct image *image, double fs, double ss) +{ + double r, rx, ry; + struct panel *p; + double xs, ys; + + p = find_panel(image->det, fs, ss); + + /* Convert xs and ys, which are in fast scan/slow scan coordinates, + * to x and y */ + xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx; + ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy; + + rx = (xs + p->cnx) / p->res; + ry = (ys + p->cny) / p->res; + + r = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + return atan2(r, p->clen); +} + + +void record_image(struct image *image, int do_poisson) +{ + int x, y; + double total_energy, energy_density; + double ph_per_e; + double area; + double max_tt = 0.0; + int n_inf1 = 0; + int n_neg1 = 0; + int n_nan1 = 0; + int n_inf2 = 0; + int n_neg2 = 0; + int n_nan2 = 0; + + /* How many photons are scattered per electron? */ + area = M_PI*pow(image->beam->beam_radius, 2.0); + total_energy = image->beam->fluence * ph_lambda_to_en(image->lambda); + energy_density = total_energy / area; + ph_per_e = (image->beam->fluence /area) * pow(THOMSON_LENGTH, 2.0); + STATUS("Fluence = %8.2e photons, " + "Energy density = %5.3f kJ/cm^2, " + "Total energy = %5.3f microJ\n", + image->beam->fluence, energy_density/1e7, total_energy*1e6); + + for ( x=0; x<image->width; x++ ) { + for ( y=0; y<image->height; y++ ) { + + double counts; + double cf; + double intensity, sa; + double pix_area, Lsq; + double xs, ys, rx, ry; + double dsq, proj_area; + struct panel *p; + + intensity = (double)image->data[x + image->width*y]; + if ( isinf(intensity) ) n_inf1++; + if ( intensity < 0.0 ) n_neg1++; + if ( isnan(intensity) ) n_nan1++; + + p = find_panel(image->det, x, y); + + /* Area of one pixel */ + pix_area = pow(1.0/p->res, 2.0); + Lsq = pow(p->clen, 2.0); + + /* Area of pixel as seen from crystal (approximate) */ + proj_area = pix_area * cos(image->twotheta[x + image->width*y]); + + /* Calculate distance from crystal to pixel */ + xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx; + ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy; + rx = (xs + p->cnx) / p->res; + ry = (ys + p->cny) / p->res; + dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + /* Projected area of pixel divided by distance squared */ + sa = proj_area / (dsq + Lsq); + + if ( do_poisson ) { + counts = poisson_noise(intensity * ph_per_e + * sa * image->beam->dqe ); + } else { + cf = intensity * ph_per_e * sa * image->beam->dqe; + counts = cf; + } + + image->data[x + image->width*y] = counts + * image->beam->adu_per_photon; + + /* Sanity checks */ + if ( isinf(image->data[x+image->width*y]) ) n_inf2++; + if ( isnan(image->data[x+image->width*y]) ) n_nan2++; + if ( image->data[x+image->width*y] < 0.0 ) n_neg2++; + + if ( image->twotheta[x + image->width*y] > max_tt ) { + max_tt = image->twotheta[x + image->width*y]; + } + + } + progress_bar(x, image->width-1, "Post-processing"); + } + + STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n", + rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); + + double tt_side = image->twotheta[(image->width/2)+image->width*0]; + STATUS("At middle of bottom edge: %.2f deg, min d = %.2f nm\n", + rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9); + + tt_side = image->twotheta[0+image->width*(image->height/2)]; + STATUS("At middle of left edge: %.2f deg, min d = %.2f nm\n", + rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9); + + STATUS("Halve the d values to get the voxel size for a synthesis.\n"); + + if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) { + ERROR("WARNING: The raw calculation produced %i negative" + " values, %i infinities and %i NaNs.\n", + n_neg1, n_inf1, n_nan1); + ERROR("WARNING: After processing, there were %i negative" + " values, %i infinities and %i NaNs.\n", + n_neg2, n_inf2, n_nan2); + } +} + + +struct panel *find_panel(struct detector *det, double fs, double ss) +{ + int p; + + for ( p=0; p<det->n_panels; p++ ) { + if ( (fs >= det->panels[p].min_fs) + && (fs <= det->panels[p].max_fs) + && (ss >= det->panels[p].min_ss) + && (ss <= det->panels[p].max_ss) ) { + return &det->panels[p]; + } + } + + return NULL; +} + + +int find_panel_number(struct detector *det, int fs, int ss) +{ + int p; + + for ( p=0; p<det->n_panels; p++ ) { + if ( (fs >= det->panels[p].min_fs) + && (fs <= det->panels[p].max_fs) + && (ss >= det->panels[p].min_ss) + && (ss <= det->panels[p].max_ss) ) return p; + } + + return -1; +} + + +void fill_in_values(struct detector *det, struct hdfile *f) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + if ( p->clen_from != NULL ) { + p->clen = get_value(f, p->clen_from) * 1.0e-3; + free(p->clen_from); + p->clen_from = NULL; + } + + p->clen += p->coffset; + + } +} + + +static struct panel *new_panel(struct detector *det, const char *name) +{ + struct panel *new; + + det->n_panels++; + det->panels = realloc(det->panels, det->n_panels*sizeof(struct panel)); + + new = &det->panels[det->n_panels-1]; + memcpy(new, &det->defaults, sizeof(struct panel)); + + /* Create a new copy of the camera length location if needed */ + if ( new->clen_from != NULL ) { + new->clen_from = strdup(new->clen_from); + } + strcpy(new->name, name); + + return new; +} + + +static struct badregion *new_bad_region(struct detector *det, const char *name) +{ + struct badregion *new; + + det->n_bad++; + det->bad = realloc(det->bad, det->n_bad*sizeof(struct badregion)); + + new = &det->bad[det->n_bad-1]; + new->min_x = NAN; + new->max_x = NAN; + new->min_y = NAN; + new->max_y = NAN; + strcpy(new->name, name); + + return new; +} + + +struct panel *find_panel_by_name(struct detector *det, const char *name) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + if ( strcmp(det->panels[i].name, name) == 0 ) { + return &det->panels[i]; + } + } + + return NULL; +} + + +static struct badregion *find_bad_region_by_name(struct detector *det, + const char *name) +{ + int i; + + for ( i=0; i<det->n_bad; i++ ) { + if ( strcmp(det->bad[i].name, name) == 0 ) { + return &det->bad[i]; + } + } + + return NULL; +} + + +static char *find_or_add_rg(struct detector *det, const char *name) +{ + int i; + char **new; + char *tmp; + + for ( i=0; i<det->num_rigid_groups; i++ ) { + + if ( strcmp(det->rigid_groups[i], name) == 0 ) { + return det->rigid_groups[i]; + } + + } + + new = realloc(det->rigid_groups, + (1+det->num_rigid_groups)*sizeof(char *)); + if ( new == NULL ) return NULL; + + det->rigid_groups = new; + + tmp = strdup(name); + det->rigid_groups[det->num_rigid_groups] = tmp; + + det->num_rigid_groups++; + + return tmp; +} + + +static int parse_field_for_panel(struct panel *panel, const char *key, + const char *val, struct detector *det) +{ + int reject = 0; + + if ( strcmp(key, "min_fs") == 0 ) { + panel->min_fs = atof(val); + } else if ( strcmp(key, "max_fs") == 0 ) { + panel->max_fs = atof(val); + } else if ( strcmp(key, "min_ss") == 0 ) { + panel->min_ss = atof(val); + } else if ( strcmp(key, "max_ss") == 0 ) { + panel->max_ss = atof(val); + } else if ( strcmp(key, "corner_x") == 0 ) { + panel->cnx = atof(val); + } else if ( strcmp(key, "corner_y") == 0 ) { + panel->cny = atof(val); + } else if ( strcmp(key, "rigid_group") == 0 ) { + panel->rigid_group = find_or_add_rg(det, val); + } else if ( strcmp(key, "clen") == 0 ) { + + char *end; + double v = strtod(val, &end); + if ( end == val ) { + /* This means "fill in later" */ + panel->clen = -1.0; + panel->clen_from = strdup(val); + } else { + panel->clen = v; + panel->clen_from = NULL; + } + + } else if ( strcmp(key, "coffset") == 0) { + panel->coffset = atof(val); + } else if ( strcmp(key, "res") == 0 ) { + panel->res = atof(val); + } else if ( strcmp(key, "peak_sep") == 0 ) { + panel->peak_sep = atof(val); + } else if ( strcmp(key, "integr_radius") == 0 ) { + panel->integr_radius = atof(val); + } else if ( strcmp(key, "badrow_direction") == 0 ) { + panel->badrow = val[0]; /* First character only */ + if ( (panel->badrow != 'x') && (panel->badrow != 'y') + && (panel->badrow != 'f') && (panel->badrow != 's') + && (panel->badrow != '-') ) { + ERROR("badrow_direction must be x, y, f, s or '-'\n"); + ERROR("Assuming '-'\n."); + panel->badrow = '-'; + } + if ( panel->badrow == 'x' ) panel->badrow = 'f'; + if ( panel->badrow == 'y' ) panel->badrow = 's'; + } else if ( strcmp(key, "no_index") == 0 ) { + panel->no_index = atob(val); + } else if ( strcmp(key, "fs") == 0 ) { + if ( dir_conv(val, &panel->fsx, &panel->fsy) != 0 ) { + ERROR("Invalid fast scan direction '%s'\n", val); + reject = 1; + } + } else if ( strcmp(key, "ss") == 0 ) { + if ( dir_conv(val, &panel->ssx, &panel->ssy) != 0 ) { + ERROR("Invalid slow scan direction '%s'\n", val); + reject = 1; + } + } else { + ERROR("Unrecognised field '%s'\n", key); + } + + return reject; +} + + +static int parse_field_bad(struct badregion *panel, const char *key, + const char *val) +{ + int reject = 0; + + if ( strcmp(key, "min_x") == 0 ) { + panel->min_x = atof(val); + } else if ( strcmp(key, "max_x") == 0 ) { + panel->max_x = atof(val); + } else if ( strcmp(key, "min_y") == 0 ) { + panel->min_y = atof(val); + } else if ( strcmp(key, "max_y") == 0 ) { + panel->max_y = atof(val); + } else { + ERROR("Unrecognised field '%s'\n", key); + } + + return reject; +} + + +static void parse_toplevel(struct detector *det, const char *key, + const char *val) +{ + if ( strcmp(key, "mask") == 0 ) { + + det->mask = strdup(val); + + } else if ( strcmp(key, "mask_bad") == 0 ) { + + char *end; + double v = strtod(val, &end); + + if ( end != val ) { + det->mask_bad = v; + } + + } else if ( strcmp(key, "mask_good") == 0 ) { + + char *end; + double v = strtod(val, &end); + + if ( end != val ) { + det->mask_good = v; + } + + } else if ( strcmp(key, "peak_sep") == 0 ) { + det->defaults.peak_sep = atof(val); + } else if ( strcmp(key, "integr_radius") == 0 ) { + det->defaults.integr_radius = atof(val); + } else if ( parse_field_for_panel(&det->defaults, key, val, det) ) { + ERROR("Unrecognised top level field '%s'\n", key); + } +} + + +struct detector *get_detector_geometry(const char *filename) +{ + FILE *fh; + struct detector *det; + char *rval; + char **bits; + int i; + int reject = 0; + int x, y, max_fs, max_ss; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return NULL; + + det = calloc(1, sizeof(struct detector)); + if ( det == NULL ) { + fclose(fh); + return NULL; + } + + det->n_panels = 0; + det->panels = NULL; + det->n_bad = 0; + det->bad = NULL; + det->mask_good = 0; + det->mask_bad = 0; + det->mask = NULL; + det->num_rigid_groups = 0; + det->rigid_groups = NULL; + + /* The default defaults... */ + det->defaults.min_fs = -1; + det->defaults.min_ss = -1; + det->defaults.max_fs = -1; + det->defaults.max_ss = -1; + det->defaults.cnx = NAN; + det->defaults.cny = NAN; + det->defaults.clen = -1.0; + det->defaults.coffset = 0.0; + det->defaults.res = -1.0; + det->defaults.badrow = '-'; + det->defaults.no_index = 0; + det->defaults.peak_sep = 50.0; + det->defaults.integr_radius = 3.0; + det->defaults.fsx = 1.0; + det->defaults.fsy = 0.0; + det->defaults.ssx = 0.0; + det->defaults.ssy = 1.0; + det->defaults.rigid_group = NULL; + strncpy(det->defaults.name, "", 1023); + + do { + + int n1, n2; + char **path; + char line[1024]; + struct badregion *badregion = NULL; + struct panel *panel = NULL; + char *key; + char wholeval[1024]; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + chomp(line); + + if ( line[0] == ';' ) continue; + + n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); + if ( n1 < 3 ) { + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + continue; + } + + /* Stitch the pieces of the "value" back together */ + wholeval[0] = '\0'; /* Empty string */ + for ( i=2; i<n1; i++ ) { + if ( bits[i][0] == ';' ) break; /* Stop on comment */ + strncat(wholeval, bits[i], 1023); + } + + if ( bits[1][0] != '=' ) { + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + continue; + } + + n2 = assplode(bits[0], "/\\.", &path, ASSPLODE_NONE); + if ( n2 < 2 ) { + /* This was a top-level option, not handled above. */ + parse_toplevel(det, bits[0], bits[2]); + for ( i=0; i<n1; i++ ) free(bits[i]); + free(bits); + for ( i=0; i<n2; i++ ) free(path[i]); + free(path); + continue; + } + + if ( strncmp(path[0], "bad", 3) == 0 ) { + badregion = find_bad_region_by_name(det, path[0]); + if ( badregion == NULL ) { + badregion = new_bad_region(det, path[0]); + } + } else { + panel = find_panel_by_name(det, path[0]); + if ( panel == NULL ) { + panel = new_panel(det, path[0]); + } + } + + key = path[1]; + + if ( panel != NULL ) { + if ( parse_field_for_panel(panel, path[1], + wholeval, det) ) + { + reject = 1; + } + } else { + if ( parse_field_bad(badregion, path[1], wholeval) ) { + reject = 1; + } + } + + for ( i=0; i<n1; i++ ) free(bits[i]); + for ( i=0; i<n2; i++ ) free(path[i]); + free(bits); + free(path); + + } while ( rval != NULL ); + + if ( det->n_panels == -1 ) { + ERROR("No panel descriptions in geometry file.\n"); + fclose(fh); + free(det); + return NULL; + } + + max_fs = 0; + max_ss = 0; + for ( i=0; i<det->n_panels; i++ ) { + + if ( det->panels[i].min_fs < 0 ) { + ERROR("Please specify the minimum FS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].max_fs < 0 ) { + ERROR("Please specify the maximum FS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].min_ss < 0 ) { + ERROR("Please specify the minimum SS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].max_ss < 0 ) { + ERROR("Please specify the maximum SS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( isnan(det->panels[i].cnx) ) { + ERROR("Please specify the corner X coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( isnan(det->panels[i].cny) ) { + ERROR("Please specify the corner Y coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( (det->panels[i].clen < 0.0) + && (det->panels[i].clen_from == NULL) ) { + ERROR("Please specify the camera length for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].res < 0 ) { + ERROR("Please specify the resolution for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + /* It's OK if the badrow direction is '0' */ + /* It's not a problem if "no_index" is still zero */ + /* The default peak_sep is OK (maybe) */ + /* The default transformation matrix is at least valid */ + + if ( det->panels[i].max_fs > max_fs ) { + max_fs = det->panels[i].max_fs; + } + if ( det->panels[i].max_ss > max_ss ) { + max_ss = det->panels[i].max_ss; + } + + } + + for ( i=0; i<det->n_bad; i++ ) { + + if ( isnan(det->bad[i].min_x) ) { + ERROR("Please specify the minimum x coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + if ( isnan(det->bad[i].min_y) ) { + ERROR("Please specify the minimum y coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + if ( isnan(det->bad[i].max_x) ) { + ERROR("Please specify the maximum x coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + if ( isnan(det->bad[i].max_y) ) { + ERROR("Please specify the maximum y coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + } + + for ( x=0; x<=max_fs; x++ ) { + for ( y=0; y<=max_ss; y++ ) { + if ( find_panel(det, x, y) == NULL ) { + ERROR("Detector geometry invalid: contains gaps.\n"); + reject = 1; + goto out; + } + } + } + +out: + det->max_fs = max_fs; + det->max_ss = max_ss; + + /* Calculate matrix inverse */ + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p; + double d; + + p = &det->panels[i]; + + if ( p->fsx*p->ssy == p->ssx*p->fsy ) { + ERROR("Panel %i transformation singular.\n", i); + reject = 1; + } + + d = (double)p->fsx*p->ssy - p->ssx*p->fsy; + p->xfs = p->ssy / d; + p->yfs = -p->ssx / d; + p->xss = -p->fsy / d; + p->yss = p->fsx / d; + + } + + if ( reject ) return NULL; + + fclose(fh); + + return det; +} + + +void free_detector_geometry(struct detector *det) +{ + int i; + + for ( i=0; i<det->num_rigid_groups; i++ ) { + free(det->rigid_groups[i]); + } + free(det->rigid_groups); + + free(det->panels); + free(det->bad); + free(det->mask); + free(det); +} + + +struct detector *copy_geom(const struct detector *in) +{ + struct detector *out; + int i; + + out = malloc(sizeof(struct detector)); + memcpy(out, in, sizeof(struct detector)); + + if ( in->mask != NULL ) { + out->mask = strdup(in->mask); + } else { + out->mask = NULL; /* = in->mask */ + } + + out->panels = malloc(out->n_panels * sizeof(struct panel)); + memcpy(out->panels, in->panels, out->n_panels * sizeof(struct panel)); + + out->bad = malloc(out->n_bad * sizeof(struct badregion)); + memcpy(out->bad, in->bad, out->n_bad * sizeof(struct badregion)); + + if ( in->rigid_groups != NULL ) { + + out->rigid_groups = malloc(out->num_rigid_groups*sizeof(char *)); + memcpy(out->rigid_groups, in->rigid_groups, + out->num_rigid_groups*sizeof(char *)); + + for ( i=0; i<in->num_rigid_groups; i++ ) { + out->rigid_groups[i] = strdup(in->rigid_groups[i]); + } + + } + + for ( i=0; i<out->n_panels; i++ ) { + + struct panel *p; + + p = &out->panels[i]; + + if ( p->clen_from != NULL ) { + /* Make a copy of the clen_from fields unique to this + * copy of the structure. */ + p->clen_from = strdup(p->clen_from); + } + + } + + for ( i=0; i<in->num_rigid_groups; i++ ) { + + int j; + char *rg = in->rigid_groups[i]; + char *rgn = out->rigid_groups[i]; + + for ( j=0; j<in->n_panels; j++ ) { + + if ( in->panels[j].rigid_group == rg ) { + out->panels[j].rigid_group = rgn; + } + + } + + } + + return out; +} + + +struct detector *simple_geometry(const struct image *image) +{ + struct detector *geom; + + geom = calloc(1, sizeof(struct detector)); + + geom->n_panels = 1; + geom->panels = calloc(1, sizeof(struct panel)); + + geom->panels[0].min_fs = 0; + geom->panels[0].max_fs = image->width-1; + geom->panels[0].min_ss = 0; + geom->panels[0].max_ss = image->height-1; + geom->panels[0].cnx = -image->width / 2.0; + geom->panels[0].cny = -image->height / 2.0; + geom->panels[0].rigid_group = NULL; + + geom->panels[0].fsx = 1; + geom->panels[0].fsy = 0; + geom->panels[0].ssx = 0; + geom->panels[0].ssy = 1; + + geom->panels[0].xfs = 1; + geom->panels[0].xss = 0; + geom->panels[0].yfs = 0; + geom->panels[0].yss = 1; + + return geom; +} + + +int reverse_2d_mapping(double x, double y, double *pfs, double *pss, + struct detector *det) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p = &det->panels[i]; + double cx, cy, fs, ss; + + /* Get position relative to corner */ + cx = x - p->cnx; + cy = y - p->cny; + + /* Reverse the transformation matrix */ + fs = cx*p->xfs + cy*p->yfs; + ss = cx*p->xss + cy*p->yss; + + /* In range? */ + if ( fs < 0 ) continue; + if ( ss < 0 ) continue; + if ( fs > (p->max_fs-p->min_fs+1) ) continue; + if ( ss > (p->max_ss-p->min_ss+1) ) continue; + + *pfs = fs + p->min_fs; + *pss = ss + p->min_ss; + return 0; + + } + + return 1; +} + + +static void check_extents(struct panel p, double *min_x, double *min_y, + double *max_x, double *max_y, double fs, double ss) +{ + double xs, ys, rx, ry; + + xs = fs*p.fsx + ss*p.ssx; + ys = fs*p.fsy + ss*p.ssy; + + rx = xs + p.cnx; + ry = ys + p.cny; + + if ( rx > *max_x ) *max_x = rx; + if ( ry > *max_y ) *max_y = ry; + if ( rx < *min_x ) *min_x = rx; + if ( ry < *min_y ) *min_y = ry; +} + + +double largest_q(struct image *image) +{ + int fs, ss; + double ttm = 0.0; + double qmax = 0.0; + + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; ss++ ) { + + struct rvec q; + double tt; + + q = get_q(image, fs, ss, &tt, 1.0/image->lambda); + + if ( tt > ttm ) { + qmax = modulus(q.u, q.v, q.w); + ttm = tt; + } + + } + } + + return qmax; +} + + +double smallest_q(struct image *image) +{ + int fs, ss; + double ttm = +INFINITY; + double qmin = +INFINITY; + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; ss++ ) { + + struct rvec q; + double tt; + + q = get_q(image, fs, ss, &tt, 1.0/image->lambda); + + if ( tt < ttm ) { + qmin = modulus(q.u, q.v, q.w); + ttm = tt; + } + + } + } + + return qmin; +} + + +void get_pixel_extents(struct detector *det, + double *min_x, double *min_y, + double *max_x, double *max_y) +{ + int i; + + *min_x = 0.0; + *max_x = 0.0; + *min_y = 0.0; + *max_y = 0.0; + + /* To determine the maximum extents of the detector, put all four + * corners of each panel through the transformations and watch for the + * biggest */ + + for ( i=0; i<det->n_panels; i++ ) { + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + 0.0, + 0.0); + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + 0.0, + det->panels[i].max_ss-det->panels[i].min_ss+1); + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + det->panels[i].max_fs-det->panels[i].min_fs+1, + 0.0); + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + det->panels[i].max_fs-det->panels[i].min_fs+1, + det->panels[i].max_ss-det->panels[i].min_ss+1); + + + } +} + + +int write_detector_geometry(const char *filename, struct detector *det) +{ + struct panel *p; + int pi; + FILE *fh; + + if ( filename == NULL ) return 2; + if ( det->n_panels < 1 ) return 3; + + fh = fopen(filename, "w"); + if ( fh == NULL ) return 1; + + for ( pi=0; pi<det->n_panels; pi++) { + + p = &(det->panels[pi]); + + if ( p == NULL ) return 4; + + if ( pi > 0 ) fprintf(fh, "\n"); + + fprintf(fh, "%s/min_fs = %d\n", p->name, p->min_fs); + fprintf(fh, "%s/min_ss = %d\n", p->name, p->min_ss); + fprintf(fh, "%s/max_fs = %d\n", p->name, p->max_fs); + fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss); + fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow); + fprintf(fh, "%s/res = %g\n", p->name, p->res); + fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep); + fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from); + fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy); + fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy); + fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx); + fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny); + + if ( p->no_index ) { + fprintf(fh, "%s/no_index = 1\n", p->name); + } /* else don't clutter up the file */ + + if ( p->rigid_group != NULL ) { + fprintf(fh, "%s/rigid_group = %s\n", + p->name, p->rigid_group); + } + + } + fclose(fh); + + return 0; +} diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h new file mode 100644 index 00000000..dd5dccec --- /dev/null +++ b/libcrystfel/src/detector.h @@ -0,0 +1,131 @@ +/* + * detector.h + * + * Detector properties + * + * (c) 2006-2011 Thomas White <taw@physics.org> + * (c) 2011 Rick Kirian <rkirian@asu.edu> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef DETECTOR_H +#define DETECTOR_H + +struct image; +struct hdfile; + +#include "hdf5-file.h" +#include "image.h" + + +struct panel +{ + char name[1024]; /* Name for this panel */ + + int min_fs; /* Smallest FS value considered to be in the panel */ + int max_fs; /* Largest FS value considered to be in this panel */ + int min_ss; /* ... and so on */ + int max_ss; + double cnx; /* Location of corner (min_fs,min_ss) in pixels */ + double cny; + double coffset; + double clen; /* Camera length in metres */ + char *clen_from; + double res; /* Resolution in pixels per metre */ + char badrow; /* 'x' or 'y' */ + int no_index; /* Don't index peaks in this panel if non-zero */ + double peak_sep; /* Characteristic peak separation */ + double integr_radius; /* Peak integration radius */ + char *rigid_group; /* Rigid group, or -1 for none */ + + double fsx; + double fsy; + double ssx; + double ssy; + + double xfs; + double yfs; + double xss; + double yss; +}; + + +struct badregion +{ + char name[1024]; + double min_x; + double max_x; + double min_y; + double max_y; +}; + + +struct detector +{ + struct panel *panels; + int n_panels; + + int max_fs; + int max_ss; /* Size of overall array needed, minus 1 */ + + struct badregion *bad; + int n_bad; + + char *mask; + unsigned int mask_bad; + unsigned int mask_good; + + char **rigid_groups; + int num_rigid_groups; + + struct panel defaults; +}; + + +extern struct rvec get_q(struct image *image, double fs, double ss, + double *ttp, double k); + +extern struct rvec get_q_for_panel(struct panel *p, double fs, double ss, + double *ttp, double k); + +extern double get_tt(struct image *image, double xs, double ys); + +extern int in_bad_region(struct detector *det, double fs, double ss); + +extern void record_image(struct image *image, int do_poisson); + +extern struct panel *find_panel(struct detector *det, double fs, double ss); + +extern int find_panel_number(struct detector *det, int fs, int ss); + +extern struct detector *get_detector_geometry(const char *filename); + +extern void free_detector_geometry(struct detector *det); + +extern struct detector *simple_geometry(const struct image *image); + +extern void get_pixel_extents(struct detector *det, + double *min_x, double *min_y, + double *max_x, double *max_y); + +extern void fill_in_values(struct detector *det, struct hdfile *f); + +extern struct detector *copy_geom(const struct detector *in); + +extern int reverse_2d_mapping(double x, double y, double *pfs, double *pss, + struct detector *det); + +extern double largest_q(struct image *image); + +extern double smallest_q(struct image *image); + +extern struct panel *find_panel_by_name(struct detector *det, const char *name); + + +#endif /* DETECTOR_H */ diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c new file mode 100644 index 00000000..355e97f1 --- /dev/null +++ b/libcrystfel/src/hdf5-file.c @@ -0,0 +1,817 @@ +/* + * hdf5-file.c + * + * Read/write HDF5 data files + * + * (c) 2006-2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <stdio.h> +#include <stdint.h> +#include <hdf5.h> +#include <assert.h> + +#include "image.h" +#include "hdf5-file.h" +#include "utils.h" + + +struct hdfile { + + const char *path; /* Current data path */ + + size_t nx; /* Image width */ + size_t ny; /* Image height */ + + hid_t fh; /* HDF file handle */ + hid_t dh; /* Dataset handle */ + + int data_open; /* True if dh is initialised */ +}; + + +struct hdfile *hdfile_open(const char *filename) +{ + struct hdfile *f; + + f = malloc(sizeof(struct hdfile)); + if ( f == NULL ) return NULL; + + /* Please stop spamming my terminal */ + H5Eset_auto2(H5E_DEFAULT, NULL, NULL); + + f->fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + if ( f->fh < 0 ) { + ERROR("Couldn't open file: %s\n", filename); + free(f); + return NULL; + } + + f->data_open = 0; + + return f; +} + + +int hdfile_set_image(struct hdfile *f, const char *path) +{ + hsize_t size[2]; + hsize_t max_size[2]; + hid_t sh; + + f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT); + if ( f->dh < 0 ) { + ERROR("Couldn't open dataset\n"); + return -1; + } + f->data_open = 1; + + sh = H5Dget_space(f->dh); + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + ERROR("Dataset is not two-dimensional\n"); + return -1; + } + H5Sget_simple_extent_dims(sh, size, max_size); + H5Sclose(sh); + + f->nx = size[0]; + f->ny = size[1]; + + return 0; +} + + +int get_peaks(struct image *image, struct hdfile *f, const char *p) +{ + hid_t dh, sh; + hsize_t size[2]; + hsize_t max_size[2]; + int i; + float *buf; + herr_t r; + int tw; + + dh = H5Dopen2(f->fh, p, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Peak list (%s) not found.\n", p); + return 1; + } + + sh = H5Dget_space(dh); + if ( sh < 0 ) { + H5Dclose(dh); + ERROR("Couldn't get dataspace for peak list.\n"); + return 1; + } + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + ERROR("Peak list has the wrong dimensionality (%i).\n", + H5Sget_simple_extent_ndims(sh)); + H5Sclose(sh); + H5Dclose(dh); + return 1; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + tw = size[1]; + if ( (tw != 3) && (tw != 4) ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Peak list has the wrong dimensions.\n"); + return 1; + } + + buf = malloc(sizeof(float)*size[0]*size[1]); + if ( buf == NULL ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Couldn't reserve memory for the peak list.\n"); + return 1; + } + r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read peak list.\n"); + free(buf); + return 1; + } + + if ( image->features != NULL ) { + image_feature_list_free(image->features); + } + image->features = image_feature_list_new(); + + for ( i=0; i<size[0]; i++ ) { + + float fs, ss, val; + struct panel *p; + + fs = buf[tw*i+0]; + ss = buf[tw*i+1]; + val = buf[tw*i+2]; + + p = find_panel(image->det, fs, ss); + if ( p == NULL ) continue; + if ( p->no_index ) continue; + + image_add_feature(image->features, fs, ss, image, val, NULL); + + } + + free(buf); + H5Sclose(sh); + H5Dclose(dh); + + return 0; +} + + +static void cleanup(hid_t fh) +{ + int n_ids, i; + hid_t ids[256]; + + n_ids = H5Fget_obj_ids(fh, H5F_OBJ_ALL, 256, ids); + for ( i=0; i<n_ids; i++ ) { + + hid_t id; + H5I_type_t type; + + id = ids[i]; + type = H5Iget_type(id); + + if ( type == H5I_GROUP ) H5Gclose(id); + if ( type == H5I_DATASET ) H5Dclose(id); + if ( type == H5I_DATATYPE ) H5Tclose(id); + if ( type == H5I_DATASPACE ) H5Sclose(id); + if ( type == H5I_ATTR ) H5Aclose(id); + + } +} + + +void hdfile_close(struct hdfile *f) +{ + if ( f->data_open ) { + H5Dclose(f->dh); + } + cleanup(f->fh); + + H5Fclose(f->fh); + free(f); +} + + +int hdf5_write(const char *filename, const void *data, + int width, int height, int type) +{ + hid_t fh, gh, sh, dh; /* File, group, dataspace and data handles */ + hid_t ph; /* Property list */ + herr_t r; + hsize_t size[2]; + hsize_t max_size[2]; + + fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + if ( fh < 0 ) { + ERROR("Couldn't create file: %s\n", filename); + return 1; + } + + gh = H5Gcreate2(fh, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + if ( gh < 0 ) { + ERROR("Couldn't create group\n"); + H5Fclose(fh); + return 1; + } + + /* Note the "swap" here, according to section 3.2.5, + * "C versus Fortran Dataspaces", of the HDF5 user's guide. */ + size[0] = height; + size[1] = width; + max_size[0] = height; + max_size[1] = width; + sh = H5Screate_simple(2, size, max_size); + + /* Set compression */ + ph = H5Pcreate(H5P_DATASET_CREATE); + H5Pset_chunk(ph, 2, size); + H5Pset_deflate(ph, 3); + + dh = H5Dcreate2(gh, "data", type, sh, + H5P_DEFAULT, ph, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't create dataset\n"); + H5Fclose(fh); + return 1; + } + + /* Muppet check */ + H5Sget_simple_extent_dims(sh, size, max_size); + + r = H5Dwrite(dh, type, H5S_ALL, + H5S_ALL, H5P_DEFAULT, data); + if ( r < 0 ) { + ERROR("Couldn't write data\n"); + H5Dclose(dh); + H5Fclose(fh); + return 1; + } + + H5Pclose(ph); + H5Gclose(gh); + H5Dclose(dh); + H5Fclose(fh); + + return 0; +} + + +static double get_wavelength(struct hdfile *f) +{ + herr_t r; + hid_t dh; + double lambda; + int nm = 1; + + dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_nm", H5P_DEFAULT); + if ( dh < 0 ) { + dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_A", H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't get photon wavelength from HDF5 file.\n"); + return -1.0; + } + nm = 0; + } + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &lambda); + H5Dclose(dh); + + if ( r < 0 ) return -1.0; + if ( isnan(lambda) ) return -1.0; + + /* Convert nm -> m */ + if ( nm ) return lambda / 1.0e9; + return lambda / 1.0e10; +} + + +static double get_i0(struct hdfile *f) +{ + herr_t r; + hid_t dh; + double i0; + + dh = H5Dopen2(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT); + if ( dh < 0 ) return -1.0; + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &i0); + H5Dclose(dh); + if ( r < 0 ) return -1.0; + + return i0; +} + + +static void debodge_saturation(struct hdfile *f, struct image *image) +{ + hid_t dh, sh; + hsize_t size[2]; + hsize_t max_size[2]; + int i; + float *buf; + herr_t r; + + dh = H5Dopen2(f->fh, "/processing/hitfinder/peakinfo_saturated", + H5P_DEFAULT); + + if ( dh < 0 ) { + /* This isn't an error */ + return; + } + + sh = H5Dget_space(dh); + if ( sh < 0 ) { + H5Dclose(dh); + ERROR("Couldn't get dataspace for saturation table.\n"); + return; + } + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + H5Sclose(sh); + H5Dclose(dh); + return; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + if ( size[1] != 3 ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Saturation table has the wrong dimensions.\n"); + return; + } + + buf = malloc(sizeof(float)*size[0]*size[1]); + if ( buf == NULL ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Couldn't reserve memory for saturation table.\n"); + return; + } + r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read saturation table.\n"); + free(buf); + return; + } + + for ( i=0; i<size[0]; i++ ) { + + unsigned int x, y; + float val; + + x = buf[3*i+0]; + y = buf[3*i+1]; + val = buf[3*i+2]; + + image->data[x+image->width*y] = val / 5.0; + image->data[x+1+image->width*y] = val / 5.0; + image->data[x-1+image->width*y] = val / 5.0; + image->data[x+image->width*(y+1)] = val / 5.0; + image->data[x+image->width*(y-1)] = val / 5.0; + + } + + free(buf); + H5Sclose(sh); + H5Dclose(dh); +} + + +int hdf5_read(struct hdfile *f, struct image *image, int satcorr) +{ + herr_t r; + float *buf; + uint16_t *flags; + hid_t mask_dh; + + /* Note the "swap" here, according to section 3.2.5, + * "C versus Fortran Dataspaces", of the HDF5 user's guide. */ + image->width = f->ny; + image->height = f->nx; + + buf = malloc(sizeof(float)*f->nx*f->ny); + + r = H5Dread(f->dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, + H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read data\n"); + free(buf); + return 1; + } + image->data = buf; + + if ( (image->det != NULL) && (image->det->mask != NULL) ) { + + mask_dh = H5Dopen2(f->fh, image->det->mask, H5P_DEFAULT); + if ( mask_dh <= 0 ) { + ERROR("Couldn't open flags\n"); + image->flags = NULL; + } else { + flags = malloc(sizeof(uint16_t)*f->nx*f->ny); + r = H5Dread(mask_dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL, + H5P_DEFAULT, flags); + if ( r < 0 ) { + ERROR("Couldn't read flags\n"); + free(flags); + image->flags = NULL; + } else { + image->flags = flags; + } + H5Dclose(mask_dh); + } + + } + + /* Read wavelength from file */ + image->lambda = get_wavelength(f); + + image->i0 = get_i0(f); + if ( image->i0 < 0.0 ) { + ERROR("Couldn't read incident intensity - using 1.0.\n"); + image->i0 = 1.0; + image->i0_available = 0; + } else { + image->i0_available = 1; + } + + if ( satcorr ) debodge_saturation(f, image); + + return 0; +} + + +static int looks_like_image(hid_t h) +{ + hid_t sh; + hsize_t size[2]; + hsize_t max_size[2]; + + sh = H5Dget_space(h); + if ( sh < 0 ) return 0; + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + return 0; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + if ( ( size[0] > 64 ) && ( size[1] > 64 ) ) return 1; + + return 0; +} + + +double get_value(struct hdfile *f, const char *name) +{ + hid_t dh; + hid_t sh; + hsize_t size; + hsize_t max_size; + hid_t type; + hid_t class; + herr_t r; + double buf; + + dh = H5Dopen2(f->fh, name, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't open data\n"); + return 0.0; + } + + type = H5Dget_type(dh); + class = H5Tget_class(type); + + if ( class != H5T_FLOAT ) { + ERROR("Not a floating point value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + sh = H5Dget_space(dh); + if ( H5Sget_simple_extent_ndims(sh) != 1 ) { + ERROR("Not a scalar value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + H5Sget_simple_extent_dims(sh, &size, &max_size); + if ( size != 1 ) { + ERROR("Not a scalar value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &buf); + if ( r < 0 ) { + ERROR("Couldn't read value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + return buf; +} + + +struct copy_hdf5_field +{ + char **fields; + int n_fields; + int max_fields; +}; + + +struct copy_hdf5_field *new_copy_hdf5_field_list() +{ + struct copy_hdf5_field *n; + + n = calloc(1, sizeof(struct copy_hdf5_field)); + if ( n == NULL ) return NULL; + + n->max_fields = 32; + n->fields = malloc(n->max_fields*sizeof(char *)); + if ( n->fields == NULL ) { + free(n); + return NULL; + } + + return n; +} + + +void free_copy_hdf5_field_list(struct copy_hdf5_field *n) +{ + int i; + for ( i=0; i<n->n_fields; i++ ) { + free(n->fields[i]); + } + free(n->fields); + free(n); +} + + +void add_copy_hdf5_field(struct copy_hdf5_field *copyme, + const char *name) +{ + assert(copyme->n_fields < copyme->max_fields); /* FIXME */ + + copyme->fields[copyme->n_fields] = strdup(name); + if ( copyme->fields[copyme->n_fields] == NULL ) { + ERROR("Failed to add field for copying '%s'\n", name); + return; + } + + copyme->n_fields++; +} + + +void copy_hdf5_fields(struct hdfile *f, const struct copy_hdf5_field *copyme, + FILE *fh) +{ + int i; + + if ( copyme == NULL ) return; + + for ( i=0; i<copyme->n_fields; i++ ) { + + char *val; + char *field; + + field = copyme->fields[i]; + val = hdfile_get_string_value(f, field); + + if ( field[0] == '/' ) { + fprintf(fh, "hdf5%s = %s\n", field, val); + } else { + fprintf(fh, "hdf5/%s = %s\n", field, val); + } + + free(val); + + } +} + + +char *hdfile_get_string_value(struct hdfile *f, const char *name) +{ + hid_t dh; + hid_t sh; + hsize_t size; + hsize_t max_size; + hid_t type; + hid_t class; + + dh = H5Dopen2(f->fh, name, H5P_DEFAULT); + if ( dh < 0 ) return NULL; + + type = H5Dget_type(dh); + class = H5Tget_class(type); + + if ( class == H5T_STRING ) { + + herr_t r; + char *tmp; + hid_t sh; + + size = H5Tget_size(type); + tmp = malloc(size+1); + + sh = H5Screate(H5S_SCALAR); + + r = H5Dread(dh, type, sh, sh, H5P_DEFAULT, tmp); + if ( r < 0 ) goto fail; + + /* Two possibilities: + * String is already zero-terminated + * String is not terminated. + * Make sure things are done properly... */ + tmp[size] = '\0'; + chomp(tmp); + + return tmp; + + } + + sh = H5Dget_space(dh); + if ( H5Sget_simple_extent_ndims(sh) != 1 ) goto fail; + + H5Sget_simple_extent_dims(sh, &size, &max_size); + if ( size != 1 ) { + H5Dclose(dh); + goto fail; + } + + switch ( class ) { + case H5T_FLOAT : { + + herr_t r; + double buf; + char *tmp; + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &buf); + if ( r < 0 ) goto fail; + + tmp = malloc(256); + snprintf(tmp, 255, "%f", buf); + + return tmp; + + } + case H5T_INTEGER : { + + herr_t r; + int buf; + char *tmp; + + r = H5Dread(dh, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &buf); + if ( r < 0 ) goto fail; + + tmp = malloc(256); + snprintf(tmp, 255, "%d", buf); + + return tmp; + } + default : { + goto fail; + } + } + +fail: + H5Tclose(type); + H5Dclose(dh); + return NULL; +} + + +char **hdfile_read_group(struct hdfile *f, int *n, const char *parent, + int **p_is_group, int **p_is_image) +{ + hid_t gh; + hsize_t num; + char **res; + int i; + int *is_group; + int *is_image; + H5G_info_t ginfo; + + gh = H5Gopen2(f->fh, parent, H5P_DEFAULT); + if ( gh < 0 ) { + *n = 0; + return NULL; + } + + if ( H5Gget_info(gh, &ginfo) < 0 ) { + /* Whoopsie */ + *n = 0; + return NULL; + } + num = ginfo.nlinks; + *n = num; + if ( num == 0 ) return NULL; + + res = malloc(num*sizeof(char *)); + is_image = malloc(num*sizeof(int)); + is_group = malloc(num*sizeof(int)); + *p_is_image = is_image; + *p_is_group = is_group; + + for ( i=0; i<num; i++ ) { + + char buf[256]; + hid_t dh; + H5I_type_t type; + + H5Lget_name_by_idx(gh, ".", H5_INDEX_NAME, H5_ITER_NATIVE, + i, buf, 255, H5P_DEFAULT); + res[i] = malloc(256); + if ( strlen(parent) > 1 ) { + snprintf(res[i], 255, "%s/%s", parent, buf); + } else { + snprintf(res[i], 255, "%s%s", parent, buf); + } /* ick */ + + is_image[i] = 0; + is_group[i] = 0; + dh = H5Oopen(gh, buf, H5P_DEFAULT); + if ( dh < 0 ) continue; + type = H5Iget_type(dh); + + if ( type == H5I_GROUP ) { + is_group[i] = 1; + } else if ( type == H5I_DATASET ) { + is_image[i] = looks_like_image(dh); + } + H5Oclose(dh); + + } + + return res; +} + + +int hdfile_set_first_image(struct hdfile *f, const char *group) +{ + char **names; + int *is_group; + int *is_image; + int n, i, j; + + names = hdfile_read_group(f, &n, group, &is_group, &is_image); + if ( n == 0 ) return 1; + + for ( i=0; i<n; i++ ) { + + if ( is_image[i] ) { + hdfile_set_image(f, names[i]); + for ( j=0; j<n; j++ ) free(names[j]); + free(is_image); + free(is_group); + free(names); + return 0; + } else if ( is_group[i] ) { + if ( !hdfile_set_first_image(f, names[i]) ) { + for ( j=0; j<n; j++ ) free(names[j]); + free(is_image); + free(is_group); + free(names); + return 0; + } + } + + } + + for ( j=0; j<n; j++ ) free(names[j]); + free(is_image); + free(is_group); + free(names); + + return 1; +} diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h new file mode 100644 index 00000000..385e919c --- /dev/null +++ b/libcrystfel/src/hdf5-file.h @@ -0,0 +1,55 @@ +/* + * hdf5-file.h + * + * Read/write HDF5 data files + * + * (c) 2006-2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef HDF5_H +#define HDF5_H + +#include <stdint.h> +#include <hdf5.h> + +#include "image.h" + +struct hdfile; + +struct copy_hdf5_field; + + +extern int hdf5_write(const char *filename, const void *data, + int width, int height, int type); + +extern int hdf5_read(struct hdfile *f, struct image *image, int satcorr); + +extern struct hdfile *hdfile_open(const char *filename); +extern int hdfile_set_image(struct hdfile *f, const char *path); +extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile, + int binning, int16_t *maxp); +extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent, + int **p_is_group, int **p_is_image); +extern int hdfile_set_first_image(struct hdfile *f, const char *group); +extern void hdfile_close(struct hdfile *f); + +extern char *hdfile_get_string_value(struct hdfile *f, const char *name); +extern int get_peaks(struct image *image, struct hdfile *f, const char *p); +extern double get_value(struct hdfile *f, const char *name); + +extern struct copy_hdf5_field *new_copy_hdf5_field_list(void); +extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f); +extern void copy_hdf5_fields(struct hdfile *f, + const struct copy_hdf5_field *copyme, FILE *fh); +extern void add_copy_hdf5_field(struct copy_hdf5_field *copyme, + const char *name); + + +#endif /* HDF5_H */ diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c new file mode 100644 index 00000000..75a8af0a --- /dev/null +++ b/libcrystfel/src/image.c @@ -0,0 +1,144 @@ +/* + * image.c + * + * Handle images and image features + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include <stdlib.h> +#include <assert.h> +#include <math.h> +#include <stdio.h> + +#include "image.h" +#include "utils.h" + +/** + * SECTION:image + * @short_description: Data structure representing an image + * @title: Image + * @section_id: + * @see_also: + * @include: "image.h" + * @Image: + * + * The <structname>image</structname> structure represents an image, usually one + * frame from a large series of diffraction patterns, which might be from the + * same or different crystals. + */ + + +struct _imagefeaturelist +{ + struct imagefeature *features; + int n_features; +}; + + +void image_add_feature(ImageFeatureList *flist, double fs, double ss, + struct image *parent, double intensity, const char *name) +{ + if ( flist->features ) { + flist->features = realloc(flist->features, + (flist->n_features+1) + *sizeof(struct imagefeature)); + } else { + assert(flist->n_features == 0); + flist->features = malloc(sizeof(struct imagefeature)); + } + + flist->features[flist->n_features].fs = fs; + flist->features[flist->n_features].ss = ss; + flist->features[flist->n_features].intensity = intensity; + flist->features[flist->n_features].parent = parent; + flist->features[flist->n_features].name = name; + flist->features[flist->n_features].valid = 1; + + flist->n_features++; + +} + + +ImageFeatureList *image_feature_list_new() +{ + ImageFeatureList *flist; + + flist = malloc(sizeof(ImageFeatureList)); + + flist->n_features = 0; + flist->features = NULL; + + return flist; +} + + +void image_feature_list_free(ImageFeatureList *flist) +{ + if ( !flist ) return; + + if ( flist->features ) free(flist->features); + free(flist); +} + + +struct imagefeature *image_feature_closest(ImageFeatureList *flist, + double fs, double ss, + double *d, int *idx) +{ + int i; + double dmin = +HUGE_VAL; + int closest = 0; + + for ( i=0; i<flist->n_features; i++ ) { + + double ds; + + ds = distance(flist->features[i].fs, flist->features[i].ss, + fs, ss); + + if ( ds < dmin ) { + dmin = ds; + closest = i; + } + + } + + if ( dmin < +HUGE_VAL ) { + *d = dmin; + *idx = closest; + return &flist->features[closest]; + } + + *d = +INFINITY; + return NULL; +} + + +int image_feature_count(ImageFeatureList *flist) +{ + if ( flist == NULL ) return 0; + return flist->n_features; +} + + +struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx) +{ + /* Sanity check */ + if ( flist == NULL ) return NULL; + if ( idx > flist->n_features ) return NULL; + + if ( flist->features[idx].valid == 0 ) return NULL; + + return &flist->features[idx]; +} + + +void image_remove_feature(ImageFeatureList *flist, int idx) +{ + flist->features[idx].valid = 0; +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h new file mode 100644 index 00000000..4d4ac2d7 --- /dev/null +++ b/libcrystfel/src/image.h @@ -0,0 +1,184 @@ +/* + * image.h + * + * Handle images and image features + * + * (c) 2006-2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef IMAGE_H +#define IMAGE_H + +#include <stdint.h> +#include <complex.h> +#include <sys/types.h> + +#include "utils.h" +#include "cell.h" +#include "detector.h" +#include "reflist.h" + + +#define MAX_CELL_CANDIDATES (32) + + +/* Structure describing a feature in an image */ +struct imagefeature { + + struct image *parent; + double fs; + double ss; + double intensity; + + /* Reciprocal space coordinates (m^-1 of course) of this feature */ + double rx; + double ry; + double rz; + + /* Internal use only */ + int valid; + + const char *name; +}; + +/* An opaque type representing a list of image features */ +typedef struct _imagefeaturelist ImageFeatureList; + + +/** + * image: + * + * <programlisting> + * struct image + * { + * float *data; + * uint16_t *flags; + * double *twotheta; + * + * UnitCell *indexed_cell; + * UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; + * int ncells; + + * struct detector *det; + * struct beam_params *beam; + * char *filename; + * const struct copy_hdf5_field *copyme; + * + * int id; + * + * double m; + * + * double lambda; + * double div; + * double bw; + * double i0; + * int i0_available; + * double osf; + * double profile_radius; + * int pr_dud; + * + * int width; + * int height; + * + * RefList *reflections; + * + * ImageFeatureList *features; + * }; + * </programlisting> + * + * The field <structfield>data</structfield> contains the raw image data, if it + * is currently available. The data might be available throughout the + * processing of an experimental pattern, but it might not be available when + * simulating, scaling or merging patterns. Similarly, + * <structfield>flags</structfield> contains an array of the same dimensions + * as <structfield>data</structfield> to contain the bad pixel flags. + * <structfield>twotheta</structfield> likewise contains an array of 2*theta + * (scattering angle) values in radians, since these values are generated as a + * by-product of the scattering vector calculation and can be used later for + * calculating intensities from differential scattering cross sections. + * + * <structfield>candidate_cells</structfield> is an array of unit cells directly + * returned by the low-level indexing system. <structfield>ncells</structfield> + * is the number of candidate unit cells which were found. The maximum number + * of cells which may be returned is <function>MAX_CELL_CANDIDATES</function>. + * <structfield>indexed_cell</structfield> contains the "correct" unit cell + * after cell reduction or matching has been performed. The job of the cell + * reduction is to convert the list of candidate cells into a single indexed + * cell, or <function>NULL</function> on failure. + * + * <structfield>copyme</structfield> represents a list of HDF5 fields to copy + * to the output stream. + **/ +struct image; + +struct image { + + float *data; + uint16_t *flags; + double *twotheta; + + UnitCell *indexed_cell; + UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; + int ncells; + + struct detector *det; + struct beam_params *beam; /* The nominal beam parameters */ + char *filename; + const struct copy_hdf5_field *copyme; + + int id; /* ID number of the thread + * handling this image */ + + /* Information about the crystal */ + double m; /* Mosaicity in radians */ + + /* Per-shot radiation values */ + double lambda; /* Wavelength in m */ + double div; /* Divergence in radians */ + double bw; /* Bandwidth as a fraction */ + double i0; /* Incident intensity */ + int i0_available; /* 0 if f0 wasn't available + * from the input. */ + double osf; /* Overall scaling factor */ + double profile_radius; /* Radius of reflection */ + int pr_dud; /* Post refinement failed */ + + int width; + int height; + + /* Integrated (or about-to-be-integrated) reflections */ + RefList *reflections; + + /* Detected peaks */ + ImageFeatureList *features; + +}; + + +/* Feature lists */ +extern ImageFeatureList *image_feature_list_new(void); + +extern void image_feature_list_free(ImageFeatureList *flist); + +extern void image_add_feature(ImageFeatureList *flist, double x, double y, + struct image *parent, double intensity, + const char *name); + +extern void image_remove_feature(ImageFeatureList *flist, int idx); + +extern struct imagefeature *image_feature_closest(ImageFeatureList *flist, + double fs, double ss, + double *d, int *idx); + +extern int image_feature_count(ImageFeatureList *flist); +extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx); + +#endif /* IMAGE_H */ diff --git a/libcrystfel/src/list_tmp.h b/libcrystfel/src/list_tmp.h new file mode 100644 index 00000000..a524b2f9 --- /dev/null +++ b/libcrystfel/src/list_tmp.h @@ -0,0 +1,106 @@ +/* + * Template for creating indexed 3D lists of a given type, usually indexed + * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX. + * + * These are used, for example, for: + * - a list of 'double complex' values for storing structure factors, + * - a list of 'double' values for storing reflection intensities, + * - a list of 'unsigned int' values for counts of some sort. + * + * When LABEL and TYPE are #defined appropriately, including this header + * defines functions such as: + * - new_list_<LABEL>(), for creating a new list of the given type, + * - set_<LABEL>(), for setting a value in a list, + * - lookup_<LABEL>(), for retrieving values from a list, + * ..and so on. + * + * See src/utils.h for more information. + * + */ + +#include <stdio.h> + +#define ERROR_T(...) fprintf(stderr, __VA_ARGS__) + +static inline void LABEL(integrate)(TYPE *ref, signed int h, + signed int k, signed int l, + TYPE i) +{ + int idx; + + if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { + ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l); + ERROR_T("You need to re-configure INDMAX and re-run.\n"); + exit(1); + } + + if ( h < 0 ) h += IDIM; + if ( k < 0 ) k += IDIM; + if ( l < 0 ) l += IDIM; + + idx = h + (IDIM*k) + (IDIM*IDIM*l); + ref[idx] += i; +} + + +static inline void LABEL(set)(TYPE *ref, signed int h, + signed int k, signed int l, + TYPE i) +{ + int idx; + + if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { + ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l); + ERROR_T("You need to re-configure INDMAX and re-run.\n"); + exit(1); + } + + if ( h < 0 ) h += IDIM; + if ( k < 0 ) k += IDIM; + if ( l < 0 ) l += IDIM; + + idx = h + (IDIM*k) + (IDIM*IDIM*l); + ref[idx] = i; +} + + +static inline TYPE LABEL(lookup)(const TYPE *ref, signed int h, + signed int k, signed int l) +{ + int idx; + + if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { + ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l); + ERROR_T("You need to re-configure INDMAX and re-run.\n"); + exit(1); + } + + if ( h < 0 ) h += IDIM; + if ( k < 0 ) k += IDIM; + if ( l < 0 ) l += IDIM; + + idx = h + (IDIM*k) + (IDIM*IDIM*l); + return ref[idx]; +} + + +static inline TYPE *LABEL(new_list)(void) +{ + TYPE *r; + size_t r_size; + r_size = IDIM*IDIM*IDIM*sizeof(TYPE); + r = malloc(r_size); + memset(r, 0, r_size); + return r; +} + + +static inline void LABEL(zero_list)(TYPE *ref) +{ + memset(ref, 0, IDIM*IDIM*IDIM*sizeof(TYPE)); +} + + +#undef LABEL +#undef TYPE +#undef ERROR_T diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c new file mode 100644 index 00000000..18856c67 --- /dev/null +++ b/libcrystfel/src/reflist.c @@ -0,0 +1,1048 @@ +/* + * reflist.c + * + * Fast reflection/peak list + * + * (c) 2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include <stdlib.h> +#include <assert.h> +#include <stdio.h> +#include <pthread.h> + +#include "reflist.h" +#include "utils.h" + +/** + * SECTION:reflist + * @short_description: The fast reflection list + * @title: RefList + * @section_id: + * @see_also: + * @include: "reflist.h" + * @Image: + * + * The fast reflection list stores reflections in an RB-tree indexed + * by the Miller indices h, k and l. Any reflection can be found in a + * length of time which scales logarithmically with the number of reflections in + * the list. + * + * A RefList can contain any number of reflections, and can store more than + * one reflection with a given set of indices, for example when two distinct + * reflections are to be stored according to their asymmetric indices. + * + * There are getters and setters which can be used to get and set values for an + * individual reflection. The reflection list does not calculate any values, + * only stores what it was given earlier. As such, you will need to carefully + * examine which fields your prior processing steps have filled in. + */ + + +struct _refldata { + + /* Symmetric indices (i.e. the "real" indices) */ + signed int hs; + signed int ks; + signed int ls; + + /* Partiality and related geometrical stuff */ + double r1; /* First excitation error */ + double r2; /* Second excitation error */ + double p; /* Partiality */ + int clamp1; /* Clamp status for r1 */ + int clamp2; /* Clamp status for r2 */ + + /* Location in image */ + double fs; + double ss; + + /* The distance from the exact Bragg position to the coordinates + * given above. */ + double excitation_error; + + /* Non-zero if this reflection can be used for scaling */ + int scalable; + + /* Non-zero if this reflection should be used as a "guide star" for + * post refinement */ + int refinable; + + /* Intensity */ + double intensity; + double esd_i; + + /* Phase */ + double phase; + int have_phase; + + /* Redundancy */ + int redundancy; + + /* User-specified temporary values */ + double temp1; + double temp2; +}; + + +enum _nodecol { + RED, + BLACK +}; + + +struct _reflection { + + /* Listy stuff */ + unsigned int serial; /* Unique serial number, key */ + struct _reflection *child[2]; /* Child nodes */ + struct _reflection *next; /* Next and previous in doubly linked */ + struct _reflection *prev; /* list of duplicate reflections */ + enum _nodecol col; /* Colour (red or black) */ + + /* Payload */ + pthread_mutex_t lock; /* Protects the contents of "data" */ + struct _refldata data; +}; + + +struct _reflist { + + struct _reflection *head; + struct _reflection *tail; + +}; + + +#define SERIAL(h, k, l) ((((h)+256)<<18) + (((k)+256)<<9) + ((l)+256)) +#define GET_H(serial) ((((serial) & 0xfffc0000)>>18)-256) +#define GET_K(serial) ((((serial) & 0x0003fe00)>>9)-256) +#define GET_L(serial) (((serial) & 0x000001ff)-256) + +/**************************** Creation / deletion *****************************/ + +static Reflection *new_node(unsigned int serial) +{ + Reflection *new; + + new = calloc(1, sizeof(struct _reflection)); + new->serial = serial; + new->next = NULL; + new->prev = NULL; + new->child[0] = NULL; + new->child[1] = NULL; + new->col = RED; + pthread_mutex_init(&new->lock, NULL); + + return new; +} + + +/** + * reflist_new: + * + * Creates a new reflection list. + * + * Returns: the new reflection list, or NULL on error. + */ +RefList *reflist_new() +{ + RefList *new; + + new = malloc(sizeof(struct _reflist)); + if ( new == NULL ) return NULL; + + new->head = NULL; + + return new; +} + + +/** + * reflection_new: + * @h: The h index of the new reflection + * @k: The k index of the new reflection + * @l: The l index of the new reflection + * + * Creates a new individual reflection. You'll probably want to use + * add_refl_to_list() at some later point. + */ +Reflection *reflection_new(signed int h, signed int k, signed int l) +{ + return new_node(SERIAL(h, k, l)); +} + + +/** + * reflection_free: + * @refl: The reflection to free. + * + * Destroys an individual reflection. + */ +void reflection_free(Reflection *refl) +{ + pthread_mutex_destroy(&refl->lock); + free(refl); +} + + +static void recursive_free(Reflection *refl) +{ + if ( refl->child[0] != NULL ) recursive_free(refl->child[0]); + if ( refl->child[1] != NULL ) recursive_free(refl->child[1]); + + while ( refl != NULL ) { + Reflection *next = refl->next; + reflection_free(refl); + refl = next; + } +} + + +/** + * reflist_free: + * @list: The reflection list to free. + * + * Destroys a reflection list. + */ +void reflist_free(RefList *list) +{ + if ( list == NULL ) return; + if ( list->head != NULL ) { + recursive_free(list->head); + } /* else empty list */ + free(list); +} + + +/********************************** Search ************************************/ + +/** + * find_refl: + * @list: The reflection list to search in + * @h: The 'h' index to search for + * @k: The 'k' index to search for + * @l: The 'l' index to search for + * + * This function finds the first reflection in 'list' with the given indices. + * + * Since a %RefList can contain multiple reflections with the same indices, you + * may need to use next_found_refl() to get the other reflections. + * + * Returns: The found reflection, or NULL if no reflection with the given + * indices could be found. + **/ +Reflection *find_refl(const RefList *list, + signed int h, signed int k, signed int l) +{ + unsigned int search = SERIAL(h, k, l); + Reflection *refl; + + if ( list->head == NULL ) return NULL; + + /* Indices greater than or equal to 256 are filtered out when + * reflections are added, so don't even bother looking. + * (also, looking for such reflections causes trouble because the search + * serial number would be invalid) */ + if ( abs(h) >= 256 ) return NULL; + if ( abs(k) >= 256 ) return NULL; + if ( abs(l) >= 256 ) return NULL; + + refl = list->head; + + while ( refl != NULL ) { + + if ( refl->serial == search ) { + + assert(search == refl->serial); + assert(h == GET_H(refl->serial)); + assert(k == GET_K(refl->serial)); + assert(l == GET_L(refl->serial)); + return refl; + + } else { + + int dir = search > refl->serial; + if ( refl->child[dir] != NULL ) { + refl = refl->child[dir]; + } else { + /* Hit the bottom of the tree */ + return NULL; + } + + } + + } + + return NULL; +} + + +/** + * next_found_refl: + * @refl: A reflection returned by find_refl() or next_found_refl() + * + * This function returns the next reflection in @refl's list with the same + * indices. + * + * Returns: The found reflection, or NULL if there are no more reflections with + * the same indices. + **/ +Reflection *next_found_refl(Reflection *refl) +{ + if ( refl->next != NULL ) assert(refl->serial == refl->next->serial); + + return refl->next; /* Well, that was easy... */ +} + + +/********************************** Getters ***********************************/ + +/** + * get_excitation_error: + * @refl: A %Reflection + * + * Returns: The excitation error for the reflection. + **/ +double get_excitation_error(const Reflection *refl) +{ + return refl->data.excitation_error; +} + + +/** + * get_detector_pos: + * @refl: A %Reflection + * @fs: Location at which to store the fast scan offset of the reflection + * @ss: Location at which to store the slow scan offset of the reflection + * + **/ +void get_detector_pos(const Reflection *refl, double *fs, double *ss) +{ + *fs = refl->data.fs; + *ss = refl->data.ss; +} + + +/** + * get_indices: + * @refl: A %Reflection + * @h: Location at which to store the 'h' index of the reflection + * @k: Location at which to store the 'k' index of the reflection + * @l: Location at which to store the 'l' index of the reflection + * + **/ +void get_indices(const Reflection *refl, + signed int *h, signed int *k, signed int *l) +{ + *h = GET_H(refl->serial); + *k = GET_K(refl->serial); + *l = GET_L(refl->serial); +} + + +/** + * get_symmetric_indices: + * @refl: A %Reflection + * @hs: Location at which to store the 'h' index of the reflection + * @ks: Location at which to store the 'k' index of the reflection + * @ls: Location at which to store the 'l' index of the reflection + * + * This function gives the symmetric indices, that is, the "real" indices before + * squashing down to the asymmetric reciprocal unit. This may be useful if the + * list is indexed according to the asymmetric indices, but you still need + * access to the symmetric version. This happens during post-refinement. + * + **/ +void get_symmetric_indices(const Reflection *refl, + signed int *hs, signed int *ks, + signed int *ls) +{ + *hs = refl->data.hs; + *ks = refl->data.ks; + *ls = refl->data.ls; +} + + +/** + * get_partiality: + * @refl: A %Reflection + * + * Returns: The partiality of the reflection. + **/ +double get_partiality(const Reflection *refl) +{ + return refl->data.p; +} + + +/** + * get_intensity: + * @refl: A %Reflection + * + * Returns: The intensity of the reflection. + **/ +double get_intensity(const Reflection *refl) +{ + return refl->data.intensity; +} + + +/** + * get_partial: + * @refl: A %Reflection + * @r1: Location at which to store the first excitation error + * @r2: Location at which to store the second excitation error + * @p: Location at which to store the partiality + * @clamp_low: Location at which to store the first clamp status + * @clamp_high: Location at which to store the second clamp status + * + * This function is used during post refinement (in conjunction with + * set_partial()) to get access to the details of the partiality calculation. + * + **/ +void get_partial(const Reflection *refl, double *r1, double *r2, double *p, + int *clamp_low, int *clamp_high) +{ + *r1 = refl->data.r1; + *r2 = refl->data.r2; + *p = get_partiality(refl); + *clamp_low = refl->data.clamp1; + *clamp_high = refl->data.clamp2; +} + + +/** + * get_scalable: + * @refl: A %Reflection + * + * Returns: non-zero if this reflection can be scaled. + * + **/ +int get_scalable(const Reflection *refl) +{ + return refl->data.scalable; +} + + +/** + * get_refinable: + * @refl: A %Reflection + * + * Returns: non-zero if this reflection can be used for post refinement. + * + **/ +int get_refinable(const Reflection *refl) +{ + return refl->data.refinable; +} + + +/** + * get_redundancy: + * @refl: A %Reflection + * + * The redundancy of the reflection is the number of measurements that have been + * made of it. Note that a redundancy of zero may have a special meaning, such + * as that the reflection was impossible to integrate. Note further that each + * reflection in the list has its own redundancy, even if there are multiple + * copies of the reflection in the list. The total number of reflection + * measurements should always be the sum of the redundancies in the entire list. + * + * Returns: the number of measurements of this reflection. + * + **/ +int get_redundancy(const Reflection *refl) +{ + return refl->data.redundancy; +} + + +/** + * get_esd_intensity: + * @refl: A %Reflection + * + * Returns: the standard error in the intensity measurement (as returned by + * get_intensity()) for this reflection. + * + **/ +double get_esd_intensity(const Reflection *refl) +{ + return refl->data.esd_i; +} + + +/** + * get_phase: + * @refl: A %Reflection + * @have_phase: Place to store a non-zero value if the phase is set, or NULL. + * + * Returns: the phase for this reflection. + * + **/ +double get_phase(const Reflection *refl, int *have_phase) +{ + if ( have_phase != NULL ) *have_phase = refl->data.have_phase; + return refl->data.phase; +} + + +/** + * get_temp1: + * @refl: A %Reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + * Returns: the first temporary value for this reflection. + * + **/ +double get_temp1(const Reflection *refl) +{ + return refl->data.temp1; +} + + +/** + * get_temp2: + * @refl: A %Reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + * Returns: the second temporary value for this reflection. + * + **/ +double get_temp2(const Reflection *refl) +{ + return refl->data.temp2; +} + + +/********************************** Setters ***********************************/ + +/** + * copy_data: + * @to: %Reflection to copy data into + * @from: %Reflection to copy data from + * + * This function is used to copy the data (which is everything listed above in + * the list of getters and setters, apart from the indices themselves) from one + * reflection to another. This might be used when creating a new list from an + * old one, perhaps using the asymmetric indices instead of the raw indices for + * the new list. + * + **/ +void copy_data(Reflection *to, const Reflection *from) +{ + memcpy(&to->data, &from->data, sizeof(struct _refldata)); +} + + +/** + * set_detector_pos: + * @refl: A %Reflection + * @exerr: The excitation error for this reflection + * @fs: The fast scan offset of the reflection + * @ss: The slow scan offset of the reflection + * + **/ +void set_detector_pos(Reflection *refl, double exerr, double fs, double ss) +{ + refl->data.excitation_error = exerr; + refl->data.fs = fs; + refl->data.ss = ss; +} + + +/** + * set_partial: + * @refl: A %Reflection + * @r1: The first excitation error + * @r2: The second excitation error + * @p: The partiality + * @clamp_low: The first clamp status + * @clamp_high: The second clamp status + * + * This function is used during post refinement (in conjunction with + * get_partial()) to get access to the details of the partiality calculation. + * + **/ +void set_partial(Reflection *refl, double r1, double r2, double p, + double clamp_low, double clamp_high) +{ + refl->data.r1 = r1; + refl->data.r2 = r2; + refl->data.p = p; + refl->data.clamp1 = clamp_low; + refl->data.clamp2 = clamp_high; +} + + +/** + * set_int: + * @refl: A %Reflection + * @intensity: The intensity for the reflection. + * + * Set the intensity for the reflection. Note that retrieval is done with + * get_intensity(). + **/ +void set_int(Reflection *refl, double intensity) +{ + refl->data.intensity = intensity; +} + + +/** + * set_scalable: + * @refl: A %Reflection + * @scalable: Non-zero if this reflection should be scaled. + * + **/ +void set_scalable(Reflection *refl, int scalable) +{ + refl->data.scalable = scalable; +} + + +/** + * set_refinable: + * @refl: A %Reflection + * @refinable: Non-zero if this reflection can be used for post refinement. + * + **/ +void set_refinable(Reflection *refl, int refinable) +{ + refl->data.refinable = refinable; +} + + +/** + * set_redundancy: + * @refl: A %Reflection + * @red: New redundancy for the reflection + * + * The redundancy of the reflection is the number of measurements that have been + * made of it. Note that a redundancy of zero may have a special meaning, such + * as that the reflection was impossible to integrate. Note further that each + * reflection in the list has its own redundancy, even if there are multiple + * copies of the reflection in the list. The total number of reflection + * measurements should always be the sum of the redundancies in the entire list. + * + **/ +void set_redundancy(Reflection *refl, int red) +{ + refl->data.redundancy = red; +} + + +/** + * set_esd_intensity: + * @refl: A %Reflection + * @esd: New standard error for this reflection's intensity measurement + * + **/ +void set_esd_intensity(Reflection *refl, double esd) +{ + refl->data.esd_i = esd; +} + + +/** + * set_ph: + * @refl: A %Reflection + * @phase: New phase for the reflection + * + **/ +void set_ph(Reflection *refl, double phase) +{ + refl->data.phase = phase; + refl->data.have_phase = 1; +} + + +/** + * set_symmetric_indices: + * @refl: A %Reflection + * @hs: The 'h' index of the reflection + * @ks: The 'k' index of the reflection + * @ls: The 'l' index of the reflection + * + * This function gives the symmetric indices, that is, the "real" indices before + * squashing down to the asymmetric reciprocal unit. This may be useful if the + * list is indexed according to the asymmetric indices, but you still need + * access to the symmetric version. This happens during post-refinement. + * + **/ +void set_symmetric_indices(Reflection *refl, + signed int hs, signed int ks, signed int ls) +{ + refl->data.hs = hs; + refl->data.ks = ks; + refl->data.ls = ls; +} + + +/** + * set_temp1 + * @refl: A %Reflection + * @temp: New temporary value for the reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + **/ +void set_temp1(Reflection *refl, double temp) +{ + refl->data.temp1 = temp; +} + + +/** + * set_temp2 + * @refl: A %Reflection + * @temp: New temporary value for the reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + **/ +void set_temp2(Reflection *refl, double temp) +{ + refl->data.temp2 = temp; +} + + +/********************************* Insertion **********************************/ + +static Reflection *rotate_once(Reflection *refl, int dir) +{ + Reflection *s = refl->child[!dir]; + + refl->child[!dir] = s->child[dir]; + s->child[dir] = refl; + + refl->col = RED; + s->col = BLACK; + + return s; +} + + +static Reflection *rotate_twice(Reflection *refl, int dir) +{ + refl->child[!dir] = rotate_once(refl->child[!dir], !dir); + return rotate_once(refl, dir); +} + + +static int is_red(Reflection *refl) +{ + return (refl != NULL) && (refl->col == RED); +} + + +static Reflection *insert_node(Reflection *refl, Reflection *new) +{ + if ( refl == NULL ) { + + refl = new; + + } else { + + int dir; + Reflection *rcd; + + assert(new->serial != refl->serial); + dir = new->serial > refl->serial; + refl->child[dir] = insert_node(refl->child[dir], new); + + rcd = refl->child[dir]; + if ( is_red(rcd) ) { + + if ( is_red(refl->child[!dir]) ) { + + refl->col = RED; + refl->child[0]->col = BLACK; + refl->child[1]->col = BLACK; + + } else { + + if ( is_red(rcd->child[dir] ) ) { + refl = rotate_once(refl, !dir); + } else if ( is_red(rcd->child[!dir] ) ) { + refl = rotate_twice(refl, !dir); + } + + } + } + + } + + return refl; +} + + +/** + * add_refl + * @list: A %RefList + * @h: The 'h' index of the reflection + * @k: The 'k' index of the reflection + * @l: The 'l' index of the reflection + * + * Adds a new reflection to @list. Note that the implementation allows there to + * be multiple reflections with the same indices in the list, so this function + * should succeed even if the given indices already feature in the list. + * + * Returns: The newly created reflection, or NULL on failure. + * + **/ +Reflection *add_refl(RefList *list, signed int h, signed int k, signed int l) +{ + Reflection *new; + Reflection *f; + + assert(abs(h)<256); + assert(abs(k)<256); + assert(abs(l)<256); + + new = new_node(SERIAL(h, k, l)); + if ( new == NULL ) return NULL; + + f = find_refl(list, h, k, l); + if ( f == NULL ) { + + list->head = insert_node(list->head, new); + list->head->col = BLACK; + + } else { + + /* New reflection is identical to a previous one */ + while ( f->next != NULL ) { + f = f->next; + } + f->next = new; + new->prev = f; + } + + return new; +} + + +/** + * add_refl_to_list + * @refl: A %Reflection + * @list: A %RefList + * + * Adds a reflection to @list. The reflection that actually gets added will be + * a newly created one, and all the data will be copied across. The original + * reflection will be destroyed and the new reflection returned. + * + * Returns: The newly created reflection, or NULL on failure. + * + **/ +Reflection *add_refl_to_list(Reflection *refl, RefList *list) +{ + signed int h, k, l; + Reflection *r_added; + + get_indices(refl, &h, &k, &l); + r_added = add_refl(list, h, k, l); + if ( r_added == NULL ) return NULL; + + copy_data(r_added, refl); + reflection_free(refl); + + return r_added; +} + + +/********************************* Iteration **********************************/ + +struct _reflistiterator { + + int stack_size; + int stack_ptr; + Reflection **stack; + +}; + + +/** + * first_refl: + * @list: A %RefList to iterate over + * @piter: Address at which to store a %RefListIterator + * + * This function sets up the state required for iteration over the entire list, + * and then returns the first reflection in the list. An iterator object will + * be created and its address stored at the location given in piter. + * + * Returns: the first reflection in the list. + * + **/ +Reflection *first_refl(RefList *list, RefListIterator **piter) +{ + RefListIterator *iter; + + iter = malloc(sizeof(struct _reflistiterator)); + iter->stack_size = 32; + iter->stack = malloc(iter->stack_size*sizeof(Reflection *)); + iter->stack_ptr = 0; + *piter = iter; + + Reflection *refl = list->head; + + do { + + if ( refl != NULL ) { + iter->stack[iter->stack_ptr++] = refl; + if ( iter->stack_ptr == iter->stack_size ) { + iter->stack_size += 32; + iter->stack = realloc(iter->stack, + iter->stack_size*sizeof(Reflection *)); + } + refl = refl->child[0]; + continue; + } + + if ( iter->stack_ptr == 0 ) { + free(iter->stack); + free(iter); + return NULL; + } + + refl = iter->stack[--iter->stack_ptr]; + + return refl; + + } while ( 1 ); +} + + +/** + * next_refl: + * @refl: A reflection + * @iter: A %RefListIterator + * + * This function looks up the next reflection in the list that was given earlier + * to first_refl(). + * + * Returns: the next reflection in the list, or NULL if no more. + * + **/ +Reflection *next_refl(Reflection *refl, RefListIterator *iter) +{ + int returned = 1; + + do { + + if ( returned ) refl = refl->child[1]; + returned = 0; + + if ( refl != NULL ) { + + iter->stack[iter->stack_ptr++] = refl; + if ( iter->stack_ptr == iter->stack_size ) { + iter->stack_size += 32; + iter->stack = realloc(iter->stack, + iter->stack_size*sizeof(Reflection *)); + } + refl = refl->child[0]; + continue; + + } + if ( iter->stack_ptr == 0 ) { + free(iter->stack); + free(iter); + return NULL; + } + + return iter->stack[--iter->stack_ptr]; + + } while ( 1 ); +} + + +/*********************************** Voodoo ***********************************/ + +static int recursive_depth(Reflection *refl) +{ + int depth_left, depth_right; + + if ( refl == NULL ) return 0; + + depth_left = recursive_depth(refl->child[0]); + depth_right = recursive_depth(refl->child[1]); + + return 1 + biggest(depth_left, depth_right); +} + + +static int recursive_count(Reflection *refl) +{ + int count_left, count_right; + + if ( refl == NULL ) return 0; + + count_left = recursive_count(refl->child[0]); + count_right = recursive_count(refl->child[1]); + + return 1 + count_left + count_right; +} + + +/** + * num_reflections: + * @list: A %RefList + * + * Returns: the number of reflections in @list. + * + **/ +int num_reflections(RefList *list) +{ + return recursive_count(list->head); +} + + +/** + * tree_depth: + * @list: A %RefList + * + * If the depth of the tree is more than about 20, access to the list will be + * slow. This should never happen. + * + * Returns: the depth of the RB-tree used internally to represent @list. + * + **/ +int tree_depth(RefList *list) +{ + return recursive_depth(list->head); +} + + +/** + * lock_reflection: + * @refl: A %Reflection + * + * Acquires a lock on the reflection. + */ +void lock_reflection(Reflection *refl) +{ + pthread_mutex_lock(&refl->lock); +} + + +/** + * unlock_reflection: + * @refl: A %Reflection + * + * Releases a lock on the reflection. + */ +void unlock_reflection(Reflection *refl) +{ + pthread_mutex_unlock(&refl->lock); +} diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h new file mode 100644 index 00000000..955db69b --- /dev/null +++ b/libcrystfel/src/reflist.h @@ -0,0 +1,112 @@ +/* + * reflist.h + * + * Fast reflection/peak list + * + * (c) 2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef REFLIST_H +#define REFLIST_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +/** + * RefList: + * + * A %RefList represents a list of Bragg reflections. + * + * This data structure is opaque. You must use the available accessor functions + * to read and write its contents. + * + **/ +typedef struct _reflist RefList; + +/** + * Reflection: + * + * A %Reflection represents a single Bragg reflection. + * + * This data structure is opaque. You must use the available accessor functions + * to read and write its contents. + * + **/ +typedef struct _reflection Reflection; + +/** + * RefListIterator: + * + * A %RefListIterator is an opaque data type used when iterating over a + * %RefList. + * + **/ +typedef struct _reflistiterator RefListIterator; + +/* Creation/deletion */ +extern RefList *reflist_new(void); +extern void reflist_free(RefList *list); +extern Reflection *reflection_new(signed int h, signed int k, signed int l); +extern void reflection_free(Reflection *refl); + +/* Search */ +extern Reflection *find_refl(const RefList *list, signed int h, signed int k, signed int l); +extern Reflection *next_found_refl(Reflection *refl); + +/* Get */ +extern double get_excitation_error(const Reflection *refl); +extern void get_detector_pos(const Reflection *refl, double *fs, double *ss); +extern double get_partiality(const Reflection *refl); +extern void get_indices(const Reflection *refl, + signed int *h, signed int *k, signed int *l); +extern void get_symmetric_indices(const Reflection *refl, + signed int *hs, signed int *ks, + signed int *ls); +extern double get_intensity(const Reflection *refl); +extern void get_partial(const Reflection *refl, double *r1, double *r2, + double *p, int *clamp_low, int *clamp_high); +extern int get_scalable(const Reflection *refl); +extern int get_refinable(const Reflection *refl); +extern int get_redundancy(const Reflection *refl); +extern double get_temp1(const Reflection *refl); +extern double get_temp2(const Reflection *refl); +extern double get_esd_intensity(const Reflection *refl); +extern double get_phase(const Reflection *refl, int *have_phase); + +/* Set */ +extern void copy_data(Reflection *to, const Reflection *from); +extern void set_detector_pos(Reflection *refl, double exerr, + double fs, double ss); +extern void set_partial(Reflection *refl, double r1, double r2, double p, + double clamp_low, double clamp_high); +extern void set_int(Reflection *refl, double intensity); +extern void set_scalable(Reflection *refl, int scalable); +extern void set_refinable(Reflection *refl, int refinable); +extern void set_redundancy(Reflection *refl, int red); +extern void set_temp1(Reflection *refl, double temp); +extern void set_temp2(Reflection *refl, double temp); +extern void set_esd_intensity(Reflection *refl, double esd); +extern void set_ph(Reflection *refl, double phase); +extern void set_symmetric_indices(Reflection *refl, + signed int hs, signed int ks, signed int ls); + +/* Insertion */ +extern Reflection *add_refl(RefList *list, + signed int h, signed int k, signed int l); +extern Reflection *add_refl_to_list(Reflection *refl, RefList *list); + +/* Iteration */ +extern Reflection *first_refl(RefList *list, RefListIterator **piter); +extern Reflection *next_refl(Reflection *refl, RefListIterator *iter); + +/* Misc */ +extern int num_reflections(RefList *list); +extern int tree_depth(RefList *list); +extern void lock_reflection(Reflection *refl); +extern void unlock_reflection(Reflection *refl); + +#endif /* REFLIST_H */ diff --git a/libcrystfel/src/thread-pool.c b/libcrystfel/src/thread-pool.c new file mode 100644 index 00000000..0ae26e17 --- /dev/null +++ b/libcrystfel/src/thread-pool.c @@ -0,0 +1,283 @@ +/* + * thread-pool.c + * + * A thread pool implementation + * + * (c) 2006-2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> +#include <pthread.h> +#include <assert.h> + +#ifdef HAVE_CPU_AFFINITY +#include <sched.h> +#endif + + +#include "utils.h" + + +/** + * SECTION:thread-pool + * @short_description: The thread pool + * @title: The thread pool + * @section_id: + * @see_also: + * @include: "thread-pool.h" + * @Image: + * + * The thread pool helps when running many tasks in parallel. It takes care of + * starting and stopping threads, and presents a relatively simple interface to + * the individual programs. + */ + +/* ------------------------------ CPU affinity ------------------------------ */ + +#ifdef HAVE_CPU_AFFINITY + +static void set_affinity(int n, int cpu_num, int cpu_groupsize, int cpu_offset) +{ + cpu_set_t c; + int group; + int n_cpu_groups; + int i; + + if ( cpu_num == 0 ) return; + + CPU_ZERO(&c); + + /* Work out which group this thread belongs to */ + group = (n / cpu_groupsize) + cpu_offset; + + /* Work out which CPUs should be used for this group */ + n_cpu_groups = cpu_num / cpu_groupsize; + group = group % n_cpu_groups; + + /* Set flags */ + for ( i=0; i<cpu_groupsize; i++ ) { + + int cpu = cpu_groupsize*group + i; + + CPU_SET(cpu, &c); + + } + + if ( sched_setaffinity(0, sizeof(cpu_set_t), &c) ) { + + /* Cannot use ERROR() just yet */ + fprintf(stderr, "%i: Failed to set CPU affinity.\n", n); + + } +} + +#else /* HAVE_CPU_AFFINITY */ + +static void set_affinity(int n, int cpu_num, int cpu_groupsize, int cpu_offset) +{ + /* Do absolutely nothing */ +} + +#endif /* HAVE_CPU_AFFINITY */ + + +/* --------------------------- Status label stuff --------------------------- */ + +static int use_status_labels = 0; +static pthread_key_t status_label_key; +pthread_mutex_t stderr_lock = PTHREAD_MUTEX_INITIALIZER; + +struct worker_args +{ + struct task_queue_range *tqr; + struct task_queue *tq; + int id; + int cpu_num; + int cpu_groupsize; + int cpu_offset; +}; + + +signed int get_status_label() +{ + int *cookie; + + if ( !use_status_labels ) { + return -1; + } + + cookie = pthread_getspecific(status_label_key); + return *cookie; +} + + +struct task_queue +{ + pthread_mutex_t lock; + + int n_started; + int n_completed; + int max; + + void *(*get_task)(void *); + void (*finalise)(void *, void *); + void *queue_args; + void (*work)(void *, int); +}; + + +static void *task_worker(void *pargsv) +{ + struct worker_args *w = pargsv; + struct task_queue *q = w->tq; + int *cookie; + + set_affinity(w->id, w->cpu_num, w->cpu_groupsize, w->cpu_offset); + + cookie = malloc(sizeof(int)); + *cookie = w->id; + pthread_setspecific(status_label_key, cookie); + + free(w); + + do { + + void *task; + int cookie; + + /* Get a task */ + pthread_mutex_lock(&q->lock); + if ( (q->max) && (q->n_started >= q->max) ) { + pthread_mutex_unlock(&q->lock); + break; + } + task = q->get_task(q->queue_args); + + /* No more tasks? */ + if ( task == NULL ) { + pthread_mutex_unlock(&q->lock); + break; + } + + q->n_started++; + pthread_mutex_unlock(&q->lock); + + cookie = *(int *)pthread_getspecific(status_label_key); + q->work(task, cookie); + + /* Update totals etc */ + pthread_mutex_lock(&q->lock); + q->n_completed++; + if ( q->finalise ) { + q->finalise(q->queue_args, task); + } + pthread_mutex_unlock(&q->lock); + + } while ( 1 ); + + free(cookie); + + return NULL; +} + + +/** + * run_threads: + * @n_threads: The number of threads to run in parallel + * @work: The function to be called to do the work + * @get_task: The function which will determine the next unassigned task + * @final: The function which will be called to clean up after a task + * @queue_args: A pointer to any data required to determine the next task + * @max: Stop calling get_task after starting this number of jobs + * @cpu_num: The number of CPUs in the system + * @cpu_groupsize: The group size into which the CPUs are grouped + * @cpu_offset: The CPU group number at which to start pinning threads + * + * 'get_task' will be called every time a worker is idle. It returns either + * NULL, indicating that no further work is available, or a pointer which will + * be passed to 'work'. + * + * 'final' will be called once per image, and will be given both queue_args + * and the last task pointer. + * + * 'get_task' and 'final' will be called only under lock, and so do NOT need to + * be re-entrant or otherwise thread safe. 'work', of course, needs to be + * thread safe. + * + * Work will stop after 'max' tasks have been processed whether get_task + * returned NULL or not. If "max" is zero, all tasks will be processed. + * + * Returns: The number of tasks completed. + **/ +int run_threads(int n_threads, TPWorkFunc work, + TPGetTaskFunc get_task, TPFinalFunc final, + void *queue_args, int max, + int cpu_num, int cpu_groupsize, int cpu_offset) +{ + pthread_t *workers; + int i; + struct task_queue q; + + pthread_key_create(&status_label_key, NULL); + + workers = malloc(n_threads * sizeof(pthread_t)); + + pthread_mutex_init(&q.lock, NULL); + q.work = work; + q.get_task = get_task; + q.finalise = final; + q.queue_args = queue_args; + q.n_started = 0; + q.n_completed = 0; + q.max = max; + + /* Now it's safe to start using the status labels */ + if ( n_threads > 1 ) use_status_labels = 1; + + /* Start threads */ + for ( i=0; i<n_threads; i++ ) { + + struct worker_args *w; + + w = malloc(sizeof(struct worker_args)); + + w->tq = &q; + w->tqr = NULL; + w->id = i; + w->cpu_num = cpu_num; + w->cpu_groupsize = cpu_groupsize; + w->cpu_offset = cpu_offset; + + if ( pthread_create(&workers[i], NULL, task_worker, w) ) { + /* Not ERROR() here */ + fprintf(stderr, "Couldn't start thread %i\n", i); + n_threads = i; + break; + } + + } + + /* Join threads */ + for ( i=0; i<n_threads; i++ ) { + pthread_join(workers[i], NULL); + } + + use_status_labels = 0; + + free(workers); + + return q.n_completed; +} diff --git a/libcrystfel/src/thread-pool.h b/libcrystfel/src/thread-pool.h new file mode 100644 index 00000000..a99a7ade --- /dev/null +++ b/libcrystfel/src/thread-pool.h @@ -0,0 +1,71 @@ +/* + * thread-pool.h + * + * A thread pool implementation + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef THREAD_POOL_H +#define THREAD_POOL_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <pthread.h> + +extern pthread_mutex_t stderr_lock; +extern signed int get_status_label(void); + + +/** + * TPGetTaskFunc: + * @qargs: The queue_args pointer which was given to run_threads(). + * Returns: A pointer which will be passed to the worker function. + * + * This function is called, non-reentrantly, to get a new work item to give to + * your work function. The stuff you need to generate the new work item should + * have been stored in @qargs which was passed to run_threads(). + * + **/ +typedef void *(*TPGetTaskFunc)(void *qargs); + + +/** + * TPWorkFunc: + * @work: The queue_args pointer which was given to run_threads(). + * @cookie: A small integral number which is guaranteed to be unique among all + * currently running threads. + * + * This function is called, reentrantly, for each work item. + * + **/ +typedef void (*TPWorkFunc)(void *work, int cookie); + + +/** + * TPFinalFunc: + * @qargs: The queue_args pointer which was given to run_threads(). + * @work: The pointer which was returned by your get_task function. + * + * This function is called, non-reentrantly, after each work item has been + * completed. A typical use might be to update some counters inside @qargs + * according to fields withing @work which were filled by your 'work' function. + * + **/ +typedef void (*TPFinalFunc)(void *qargs, void *work); + + +extern int run_threads(int n_threads, TPWorkFunc work, + TPGetTaskFunc get_task, TPFinalFunc final, + void *queue_args, int max, + int cpu_num, int cpu_groupsize, int cpu_offset); + + +#endif /* THREAD_POOL_H */ diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c new file mode 100644 index 00000000..6b29627f --- /dev/null +++ b/libcrystfel/src/utils.c @@ -0,0 +1,488 @@ +/* + * utils.c + * + * Utility stuff + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#include <libgen.h> +#include <math.h> +#include <string.h> +#include <stdio.h> +#include <unistd.h> +#include <sys/types.h> +#include <sys/stat.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> + +#include "utils.h" +#include "image.h" + + +/** + * SECTION:utils + * @short_description: Miscellaneous utilities + * @title: Utilities + * @section_id: + * @see_also: + * @include: "utils.h" + * @Image: + * + * Wibble + */ + +/** + * show_matrix_eqn: + * @M: A matrix + * @v: A vector + * @r: The number of elements in @v and the side length of @M. + * + * Displays a matrix equation of the form @M.a = @v. + * + * @M must be square. + **/ +void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r) +{ + int i, j; + + for ( i=0; i<r; i++ ) { + STATUS("[ "); + for ( j=0; j<r; j++ ) { + STATUS("%+9.3e ", gsl_matrix_get(M, i, j)); + } + STATUS("][ a%2i ] = [ %+9.3e ]\n", i, gsl_vector_get(v, i)); + } +} + + +size_t notrail(char *s) +{ + size_t i; + size_t munched = 0; + + for ( i=strlen(s)-1; i>=0; i-- ) { + if ( (s[i] == ' ') || (s[i] == '\t') ) { + s[i] = '\0'; + munched++; + } else { + return munched; + } + } + + return munched; +} + + +void chomp(char *s) +{ + size_t i; + + if ( !s ) return; + + for ( i=0; i<strlen(s); i++ ) { + if ( (s[i] == '\n') || (s[i] == '\r') ) { + s[i] = '\0'; + return; + } + } +} + + +void progress_bar(int val, int total, const char *text) +{ + double frac; + int n, i; + char s[1024]; + const int width = 50; + + if ( total == 0 ) return; + + if ( !isatty(STDERR_FILENO) ) return; + if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; + + frac = (double)val/total; + n = (int)(frac*width); + + for ( i=0; i<n; i++ ) s[i] = '='; + for ( i=n; i<width; i++ ) s[i] = '.'; + s[width] = '\0'; + + pthread_mutex_lock(&stderr_lock); + fprintf(stderr, "\r%s: |%s|", text, s); + if ( val == total ) fprintf(stderr, "\n"); + pthread_mutex_unlock(&stderr_lock); + + fflush(stdout); +} + + +double random_flat(double max) +{ + return max * (double)random()/RAND_MAX; +} + + +double flat_noise(double expected, double width) +{ + double noise = random_flat(2.0*width); + return expected+noise-width; +} + + +double gaussian_noise(double expected, double stddev) +{ + double x1, x2, noise; + + /* A uniformly distributed random number between 0 and 1 */ + x1 = ((double)random()/RAND_MAX); + x2 = ((double)random()/RAND_MAX); + + noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2); + + return expected + noise*stddev; +} + + +static int fake_poisson_noise(double expected) +{ + double rf = gaussian_noise(expected, sqrt(expected)); + return (int)rf; +} + + +int poisson_noise(double expected) +{ + double L; + int k = 0; + double p = 1.0; + + /* For large values of the mean, we get big problems with arithmetic. + * In such cases, fall back on a Gaussian with the right variance. */ + if ( expected > 100.0 ) return fake_poisson_noise(expected); + + L = exp(-expected); + + do { + + double r; + + k++; + r = (double)random()/(double)RAND_MAX; + p *= r; + + } while ( p > L ); + + return k - 1; +} + + +/** + * SECTION:quaternion + * @short_description: Simple quaternion handling + * @title: Quaternion + * @section_id: + * @see_also: + * @include: "utils.h" + * @Image: + * + * There is a simple quaternion structure in CrystFEL. At the moment, it is + * only used when simulating patterns, as an argument to cell_rotate() to + * orient the unit cell. + */ + +/** + * quaternion_modulus: + * @q: A %quaternion + * + * If a quaternion represents a pure rotation, its modulus should be unity. + * + * Returns: the modulus of the given quaternion. + **/ +double quaternion_modulus(struct quaternion q) +{ + return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); +} + + +/** + * normalise_quaternion: + * @q: A %quaternion + * + * Rescales the quaternion such that its modulus is unity. + * + * Returns: the normalised version of @q + **/ +struct quaternion normalise_quaternion(struct quaternion q) +{ + double mod; + struct quaternion r; + + mod = quaternion_modulus(q); + + r.w = q.w / mod; + r.x = q.x / mod; + r.y = q.y / mod; + r.z = q.z / mod; + + return r; +} + + +/** + * random_quaternion: + * + * Returns: a randomly generated, normalised, quaternion. + **/ +struct quaternion random_quaternion() +{ + struct quaternion q; + + q.w = 2.0*(double)random()/RAND_MAX - 1.0; + q.x = 2.0*(double)random()/RAND_MAX - 1.0; + q.y = 2.0*(double)random()/RAND_MAX - 1.0; + q.z = 2.0*(double)random()/RAND_MAX - 1.0; + q = normalise_quaternion(q); + + return q; +} + + +/** + * quaternion_valid: + * @q: A %quaternion + * + * Checks if the given quaternion is normalised. + * + * This function performs a nasty floating point comparison of the form + * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be + * relied upon to spot anything other than the most obvious input error. + * + * Returns: 1 if the quaternion is normalised, 0 if not. + **/ +int quaternion_valid(struct quaternion q) +{ + double qmod; + + qmod = quaternion_modulus(q); + + /* Modulus = 1 to within some tolerance? + * Nasty allowance for floating-point accuracy follows... */ + if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; + + return 0; +} + + +/** + * quat_rot + * @q: A vector (in the form of an %rvec) + * @z: A %quaternion + * + * Rotates a vector according to a quaternion. + * + * Returns: A rotated version of @p. + **/ +struct rvec quat_rot(struct rvec q, struct quaternion z) +{ + struct rvec res; + double t01, t02, t03, t11, t12, t13, t22, t23, t33; + + t01 = z.w*z.x; + t02 = z.w*z.y; + t03 = z.w*z.z; + t11 = z.x*z.x; + t12 = z.x*z.y; + t13 = z.x*z.z; + t22 = z.y*z.y; + t23 = z.y*z.z; + t33 = z.z*z.z; + + res.u = (1.0 - 2.0 * (t22 + t33)) * q.u + + (2.0 * (t12 + t03)) * q.v + + (2.0 * (t13 - t02)) * q.w; + + res.v = (2.0 * (t12 - t03)) * q.u + + (1.0 - 2.0 * (t11 + t33)) * q.v + + (2.0 * (t01 + t23)) * q.w; + + res.w = (2.0 * (t02 + t13)) * q.u + + (2.0 * (t23 - t01)) * q.v + + (1.0 - 2.0 * (t11 + t22)) * q.w; + + return res; +} + + +/* Return non-zero if c is in delims */ +static int assplode_isdelim(const char c, const char *delims) +{ + size_t i; + for ( i=0; i<strlen(delims); i++ ) { + if ( c == delims[i] ) return 1; + } + return 0; +} + + +static int assplode_extract(char ***pbits, int n, size_t n_captured, + size_t start, const char *a) +{ + char **bits = *pbits; + bits = realloc(bits, sizeof(char *)*(n+1)); + bits[n] = malloc(n_captured+1); + memcpy(bits[n], a+start, n_captured); + bits[n][n_captured] = '\0'; + n++; + *pbits = bits; + return n; +} + + +/* Split the string 'a' using 'delims' as a zero-terminated list of + * deliminators. + * Store each segment in bits[0...n] where n is the number of segments and is + * the return value. pbits = &bits + * Each segment needs to be freed with free() when finished with. + * The array of bits also needs to be freed with free() when finished with, + * unless n=0 in which case bits==NULL + */ +int assplode(const char *a, const char *delims, char ***pbits, + AssplodeFlag flags) +{ + size_t i, start, n_captured; + int n, last_was_delim; + char **bits; + + n = 0; + i = 0; + n_captured = 0; + start = 0; + last_was_delim = 0; + bits = NULL; + while ( i < strlen(a) ) { + + if ( assplode_isdelim(a[i], delims) ) { + + if ( n_captured > 0 ) { + /* This is a deliminator after a sequence of + * non-deliminator chars */ + n = assplode_extract(&bits, n, n_captured, + start, a); + } + + n_captured = 0; + if ( (flags & ASSPLODE_DUPS) && last_was_delim ) { + n = assplode_extract(&bits, n, 0, start, a); + } + last_was_delim = 1; + + } else { + + if ( n_captured == 0 ) { + /* No characters currently found, so this is the + * start */ + start = i; + } + n_captured++; + last_was_delim = 0; + + } + + i++; + + } + /* Left over characters at the end? */ + if ( n_captured > 0 ) { + n = assplode_extract(&bits, n, n_captured, start, a); + } + + *pbits = bits; + return n; +} + + +char *check_prefix(char *prefix) +{ + int r; + struct stat statbuf; + char *new; + size_t len; + + /* Is "prefix" a directory? */ + r = stat(prefix, &statbuf); + if ( r != 0 ) { + /* "prefix" probably doesn't exist. This is fine - assume + * the user knows what they're doing, and that "prefix" + * suffixed with the actual filename will produce something + * sensible. */ + return prefix; + } + + if ( !S_ISDIR(statbuf.st_mode) ) { + /* Also fine, as above. */ + return prefix; + } + + /* Does the prefix end in a slash? */ + if ( prefix[strlen(prefix)-1] == '/' ) { + /* This looks sensible. */ + return prefix; + } + + STATUS("Your prefix ('%s') is a directory, but doesn't end" + " with a slash. I'm going to add it for you.\n", prefix); + STATUS("If this isn't what you want, run with --no-check-prefix.\n"); + len = strlen(prefix)+2; + new = malloc(len); + snprintf(new, len, "%s/", prefix); + free(prefix); + return new; +} + + +char *safe_basename(const char *in) +{ + int i; + char *cpy; + char *res; + + cpy = strdup(in); + + /* Get rid of any trailing slashes */ + for ( i=strlen(cpy)-1; i>0; i-- ) { + if ( cpy[i] == '/' ) { + cpy[i] = '\0'; + } else { + break; + } + } + + /* Find the base name */ + for ( i=strlen(cpy)-1; i>0; i-- ) { + if ( cpy[i] == '/' ) { + i++; + break; + } + } + + res = strdup(cpy+i); + /* If we didn't find a previous slash, i==0 so res==cpy */ + + free(cpy); + + return res; +} + + +#ifdef GSL_FUDGE +/* Force the linker to bring in CBLAS to make GSL happy */ +void utils_fudge_gslcblas() +{ + STATUS("%p\n", cblas_sgemm); +} +#endif diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h new file mode 100644 index 00000000..1179a57f --- /dev/null +++ b/libcrystfel/src/utils.h @@ -0,0 +1,236 @@ +/* + * utils.h + * + * Utility stuff + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef UTILS_H +#define UTILS_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <math.h> +#include <complex.h> +#include <string.h> +#include <stdlib.h> +#include <pthread.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> + + +#include "thread-pool.h" + + +/* -------------------------- Fundamental constants ------------------------ */ + +/* Electron charge in C */ +#define ELECTRON_CHARGE (1.6021773e-19) + +/* Planck's constant (Js) */ +#define PLANCK (6.62606896e-34) + +/* Speed of light in vacuo (m/s) */ +#define C_VACUO (299792458) + +/* Thomson scattering length (m) */ +#define THOMSON_LENGTH (2.81794e-15) + + +/* ------------------------------ Quaternions ------------------------------- */ + +/** + * quaternion: + * + * <programlisting> + * struct quaternion + * { + * double w + * double x + * double y + * double z + * }; + * </programlisting> + * + * A structure representing a quaternion. + * + **/ +struct quaternion; + +struct quaternion { + double w; + double x; + double y; + double z; +}; + +extern struct quaternion normalise_quaternion(struct quaternion q); +extern double quaternion_modulus(struct quaternion q); +extern struct quaternion random_quaternion(void); +extern int quaternion_valid(struct quaternion q); +extern struct rvec quat_rot(struct rvec q, struct quaternion z); + + +/* --------------------------- Useful functions ----------------------------- */ + +extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r); +extern size_t notrail(char *s); +extern void chomp(char *s); + +typedef enum { + ASSPLODE_NONE = 0, + ASSPLODE_DUPS = 1<<0 +} AssplodeFlag; +extern int assplode(const char *a, const char *delims, char ***pbits, + AssplodeFlag flags); + +extern void progress_bar(int val, int total, const char *text); +extern double random_flat(double max); +extern double flat_noise(double expected, double width); +extern double gaussian_noise(double expected, double stddev); +extern int poisson_noise(double expected); + +/* Keep these ones inline, to avoid function call overhead */ +static inline struct quaternion invalid_quaternion(void) +{ + struct quaternion quat; + quat.w = 0.0; + quat.x = 0.0; + quat.y = 0.0; + quat.z = 0.0; + return quat; +} + +static inline double distance(double x1, double y1, double x2, double y2) +{ + return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); +} + +static inline double modulus(double x, double y, double z) +{ + return sqrt(x*x + y*y + z*z); +} + +static inline double modulus_squared(double x, double y, double z) { + return x*x + y*y + z*z; +} + +static inline double distance3d(double x1, double y1, double z1, + double x2, double y2, double z2) +{ + return modulus(x1-x2, y1-y2, z1-z2); +} + +/* Answer in radians */ +static inline double angle_between(double x1, double y1, double z1, + double x2, double y2, double z2) +{ + double mod1 = modulus(x1, y1, z1); + double mod2 = modulus(x2, y2, z2); + double cosine = (x1*x2 + y1*y2 + z1*z2) / (mod1*mod2); + + /* Fix domain if outside due to rounding */ + if ( cosine > 1.0 ) cosine = 1.0; + if ( cosine < -1.0 ) cosine = -1.0; + + return acos(cosine); +} + +static inline int within_tolerance(double a, double b, double percent) +{ + double tol = fabs(a) * (percent/100.0); + if ( fabs(b-a) < tol ) return 1; + return 0; +} + + +/* ----------------------------- Useful macros ------------------------------ */ + +#define rad2deg(a) ((a)*180/M_PI) +#define deg2rad(a) ((a)*M_PI/180) + +#define is_odd(a) ((a)%2==1) + +#define biggest(a,b) ((a>b) ? (a) : (b)) +#define smallest(a,b) ((a<b) ? (a) : (b)) + + +/* Photon energy (J) to wavelength (m) */ +#define ph_en_to_lambda(a) ((PLANCK*C_VACUO)/(a)) + +/* Photon wavelength (m) to energy (J) */ +#define ph_lambda_to_en(a) ((PLANCK*C_VACUO)/(a)) + +/* eV to Joules */ +#define eV_to_J(a) ((a)*ELECTRON_CHARGE) + +/* Joules to eV */ +#define J_to_eV(a) ((a)/ELECTRON_CHARGE) + +#define UNUSED __attribute__((unused)) + +/* -------------------- Indexed lists for specified types ------------------- */ + +#include "defs.h" + +#define LIST_SIZE (IDIM*IDIM*IDIM) + +/* Create functions for storing reflection intensities indexed as h,k,l */ +#define LABEL(x) x##_intensity +#define TYPE double +#include "list_tmp.h" + +/* CAs above, but for phase values */ +#define LABEL(x) x##_phase +#define TYPE double +#include "list_tmp.h" + +/* As above, but for (unsigned) integer counts */ +#define LABEL(x) x##_count +#define TYPE unsigned int +#include "list_tmp.h" + +/* As above, but for simple flags */ +#define LABEL(x) x##_flag +#define TYPE unsigned char +#include "list_tmp.h" + + +/* ------------------------------ Message macros ---------------------------- */ + +extern pthread_mutex_t stderr_lock; + +#define ERROR(...) { \ + int error_print_val = get_status_label(); \ + pthread_mutex_lock(&stderr_lock); \ + if ( error_print_val >= 0 ) { \ + fprintf(stderr, "%3i: ", error_print_val); \ + } \ + fprintf(stderr, __VA_ARGS__); \ + pthread_mutex_unlock(&stderr_lock); \ + } + +#define STATUS(...) { \ + int status_print_val = get_status_label(); \ + pthread_mutex_lock(&stderr_lock); \ + if ( status_print_val >= 0 ) { \ + fprintf(stderr, "%3i: ", status_print_val); \ + } \ + fprintf(stderr, __VA_ARGS__); \ + pthread_mutex_unlock(&stderr_lock); \ + } + + +/* ------------------------------ File handling ----------------------------- */ + +extern char *check_prefix(char *prefix); +extern char *safe_basename(const char *in); + + +#endif /* UTILS_H */ |