From 25c3d29ed7701cadbb3813878f25b633a7cd7c2d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 13:59:17 +0100 Subject: Move indexing and rendering to libcrystfel --- Makefile.am | 12 +- Makefile.in | 46 +-- configure | 2 +- libcrystfel/Makefile.am | 4 +- libcrystfel/Makefile.in | 27 +- libcrystfel/src/dirax.c | 529 +++++++++++++++++++++++++++++++ libcrystfel/src/dirax.h | 26 ++ libcrystfel/src/index-priv.h | 29 ++ libcrystfel/src/index.c | 274 ++++++++++++++++ libcrystfel/src/index.h | 63 ++++ libcrystfel/src/mosflm.c | 609 ++++++++++++++++++++++++++++++++++++ libcrystfel/src/mosflm.h | 27 ++ libcrystfel/src/reax.c | 728 +++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/reax.h | 27 ++ libcrystfel/src/render.c | 442 ++++++++++++++++++++++++++ libcrystfel/src/render.h | 55 ++++ src/dirax.c | 529 ------------------------------- src/dirax.h | 26 -- src/index-priv.h | 29 -- src/index.c | 274 ---------------- src/index.h | 63 ---- src/mosflm.c | 609 ------------------------------------ src/mosflm.h | 27 -- src/reax.c | 728 ------------------------------------------- src/reax.h | 27 -- src/render.c | 442 -------------------------- src/render.h | 55 ---- 27 files changed, 2854 insertions(+), 2855 deletions(-) create mode 100644 libcrystfel/src/dirax.c create mode 100644 libcrystfel/src/dirax.h create mode 100644 libcrystfel/src/index-priv.h create mode 100644 libcrystfel/src/index.c create mode 100644 libcrystfel/src/index.h create mode 100644 libcrystfel/src/mosflm.c create mode 100644 libcrystfel/src/mosflm.h create mode 100644 libcrystfel/src/reax.c create mode 100644 libcrystfel/src/reax.h create mode 100644 libcrystfel/src/render.c create mode 100644 libcrystfel/src/render.h delete mode 100644 src/dirax.c delete mode 100644 src/dirax.h delete mode 100644 src/index-priv.h delete mode 100644 src/index.c delete mode 100644 src/index.h delete mode 100644 src/mosflm.c delete mode 100644 src/mosflm.h delete mode 100644 src/reax.c delete mode 100644 src/reax.h delete mode 100644 src/render.c delete mode 100644 src/render.h diff --git a/Makefile.am b/Makefile.am index b119e288..53321656 100644 --- a/Makefile.am +++ b/Makefile.am @@ -42,11 +42,10 @@ endif src_process_hkl_SOURCES = src/process_hkl.c -src_indexamajig_SOURCES = src/indexamajig.c src/index.c src/dirax.c \ - src/mosflm.c src/reax.c +src_indexamajig_SOURCES = src/indexamajig.c if BUILD_HDFSEE -src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c +src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c endif src_get_hkl_SOURCES = src/get_hkl.c @@ -57,7 +56,7 @@ src_check_hkl_SOURCES = src/check_hkl.c src_powder_plot_SOURCES = src/powder_plot.c -src_render_hkl_SOURCES = src/render_hkl.c src/povray.c src/render.c +src_render_hkl_SOURCES = src/render_hkl.c src/povray.c src_sum_stack_SOURCES = src/sum_stack.c @@ -82,10 +81,9 @@ tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \ INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data -EXTRA_DIST += src/dw-hdfsee.h src/render.h src/hdfsee.h src/dirax.h \ - src/mosflm.h src/index.h src/povray.h src/index-priv.h \ +EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/povray.h \ src/render_hkl.h src/post-refinement.h src/hrs-scaling.h \ - src/scaling-report.h src/reax.h + src/scaling-report.h crystfeldir = $(datadir)/crystfel crystfel_DATA = data/diffraction.cl data/defs.h data/hdfsee.ui diff --git a/Makefile.in b/Makefile.in index e054cd4b..bffe2b0d 100644 --- a/Makefile.in +++ b/Makefile.in @@ -114,18 +114,14 @@ src_get_hkl_OBJECTS = $(am_src_get_hkl_OBJECTS) src_get_hkl_LDADD = $(LDADD) src_get_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ $(top_builddir)/libcrystfel/libcrystfel.la -am__src_hdfsee_SOURCES_DIST = src/hdfsee.c src/dw-hdfsee.c \ - src/render.c +am__src_hdfsee_SOURCES_DIST = src/hdfsee.c src/dw-hdfsee.c @BUILD_HDFSEE_TRUE@am_src_hdfsee_OBJECTS = src/hdfsee.$(OBJEXT) \ -@BUILD_HDFSEE_TRUE@ src/dw-hdfsee.$(OBJEXT) \ -@BUILD_HDFSEE_TRUE@ src/render.$(OBJEXT) +@BUILD_HDFSEE_TRUE@ src/dw-hdfsee.$(OBJEXT) src_hdfsee_OBJECTS = $(am_src_hdfsee_OBJECTS) src_hdfsee_LDADD = $(LDADD) src_hdfsee_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ $(top_builddir)/libcrystfel/libcrystfel.la -am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ - src/index.$(OBJEXT) src/dirax.$(OBJEXT) src/mosflm.$(OBJEXT) \ - src/reax.$(OBJEXT) +am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS) src_indexamajig_LDADD = $(LDADD) src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ @@ -161,7 +157,7 @@ src_process_hkl_LDADD = $(LDADD) src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ $(top_builddir)/libcrystfel/libcrystfel.la am_src_render_hkl_OBJECTS = src/render_hkl.$(OBJEXT) \ - src/povray.$(OBJEXT) src/render.$(OBJEXT) + src/povray.$(OBJEXT) src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS) src_render_hkl_LDADD = $(LDADD) src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ @@ -551,11 +547,10 @@ EXTRA_DIST = configure m4/gnulib-cache.m4 tests/first_merge_check \ tests/second_merge_check tests/third_merge_check \ tests/fourth_merge_check tests/first_merge_check.hkl \ tests/third_merge_check.hkl tests/fourth_merge_check.hkl \ - src/dw-hdfsee.h src/render.h src/hdfsee.h src/dirax.h \ - src/mosflm.h src/index.h src/povray.h src/index-priv.h \ - src/render_hkl.h src/post-refinement.h src/hrs-scaling.h \ - src/scaling-report.h src/reax.h $(crystfel_DATA) \ - doc/twin-calculator.pdf doc/examples/lcls-dec.geom \ + src/dw-hdfsee.h src/hdfsee.h src/povray.h src/render_hkl.h \ + src/post-refinement.h src/hrs-scaling.h src/scaling-report.h \ + $(crystfel_DATA) doc/twin-calculator.pdf \ + doc/examples/lcls-dec.geom \ doc/examples/lcls-june-r0013-r0128.geom \ doc/examples/lcls-xpp-estimate.geom doc/examples/simple.geom \ doc/examples/lcls-dec.beam doc/examples/lcls-june.beam \ @@ -579,15 +574,13 @@ src_partial_sim_SOURCES = src/partial_sim.c src_pattern_sim_SOURCES = src/pattern_sim.c @HAVE_OPENCL_TRUE@tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src_process_hkl_SOURCES = src/process_hkl.c -src_indexamajig_SOURCES = src/indexamajig.c src/index.c src/dirax.c \ - src/mosflm.c src/reax.c - -@BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c +src_indexamajig_SOURCES = src/indexamajig.c +@BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src_get_hkl_SOURCES = src/get_hkl.c src_compare_hkl_SOURCES = src/compare_hkl.c src_check_hkl_SOURCES = src/check_hkl.c src_powder_plot_SOURCES = src/powder_plot.c -src_render_hkl_SOURCES = src/render_hkl.c src/povray.c src/render.c +src_render_hkl_SOURCES = src/render_hkl.c src/povray.c src_sum_stack_SOURCES = src/sum_stack.c src_calibrate_detector_SOURCES = src/calibrate_detector.c src_partialator_SOURCES = src/partialator.c src/post-refinement.c \ @@ -763,18 +756,11 @@ src/hdfsee.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/dw-hdfsee.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) -src/render.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) src/hdfsee$(EXEEXT): $(src_hdfsee_OBJECTS) $(src_hdfsee_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/hdfsee$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_hdfsee_OBJECTS) $(src_hdfsee_LDADD) $(LIBS) src/indexamajig.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) -src/index.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/dirax.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/mosflm.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) -src/reax.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/indexamajig$(EXEEXT): $(src_indexamajig_OBJECTS) $(src_indexamajig_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/indexamajig$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_indexamajig_OBJECTS) $(src_indexamajig_LDADD) $(LIBS) @@ -858,14 +844,11 @@ mostlyclean-compile: -rm -f src/calibrate_detector.$(OBJEXT) -rm -f src/check_hkl.$(OBJEXT) -rm -f src/compare_hkl.$(OBJEXT) - -rm -f src/dirax.$(OBJEXT) -rm -f src/dw-hdfsee.$(OBJEXT) -rm -f src/get_hkl.$(OBJEXT) -rm -f src/hdfsee.$(OBJEXT) -rm -f src/hrs-scaling.$(OBJEXT) - -rm -f src/index.$(OBJEXT) -rm -f src/indexamajig.$(OBJEXT) - -rm -f src/mosflm.$(OBJEXT) -rm -f src/partial_sim.$(OBJEXT) -rm -f src/partialator.$(OBJEXT) -rm -f src/pattern_sim.$(OBJEXT) @@ -873,8 +856,6 @@ mostlyclean-compile: -rm -f src/povray.$(OBJEXT) -rm -f src/powder_plot.$(OBJEXT) -rm -f src/process_hkl.$(OBJEXT) - -rm -f src/reax.$(OBJEXT) - -rm -f src/render.$(OBJEXT) -rm -f src/render_hkl.$(OBJEXT) -rm -f src/scaling-report.$(OBJEXT) -rm -f src/sum_stack.$(OBJEXT) @@ -890,14 +871,11 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/calibrate_detector.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/check_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/compare_hkl.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/dirax.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/dw-hdfsee.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/get_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/hdfsee.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/hrs-scaling.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/index.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/indexamajig.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/mosflm.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/partial_sim.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/partialator.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/pattern_sim.Po@am__quote@ @@ -905,8 +883,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/povray.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/powder_plot.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/process_hkl.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reax.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/scaling-report.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sum_stack.Po@am__quote@ diff --git a/configure b/configure index 7689717f..860ad20e 100755 --- a/configure +++ b/configure @@ -14231,7 +14231,7 @@ fi if test "${with_html_dir+set}" = set; then : withval=$with_html_dir; else - with_html_dir='${docdir}/reference' + with_html_dir='${datadir}/gtk-doc/html' fi HTML_DIR="$with_html_dir" diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 00428a1f..576044a6 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -4,7 +4,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/beam-parameters.c src/geometry.c src/statistics.c \ src/symmetry.c src/stream.c src/peaks.c \ src/reflist-utils.c src/filters.c src/diffraction.c \ - src/diffraction-gpu.c src/cl-utils.c + src/diffraction-gpu.c src/cl-utils.c src/render.c \ + src/index.c src/dirax.c src/mosflm.c src/reax.c libcrystfel_la_includedir=$(includedir)/crystfel/ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/diffraction-gpu.h \ @@ -14,6 +15,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/diffraction-gpu.h \ src/cl-utils.h src/filters.h src/list_tmp.h \ src/statistics.h src/utils.h src/detector.h \ src/geometry.h src/peaks.h src/stream.h \ + src/render.h src/index.h \ ../data/defs.h INCLUDES = "-I$(top_srcdir)/data" diff --git a/libcrystfel/Makefile.in b/libcrystfel/Makefile.in index 7d27bd59..a0b527e4 100644 --- a/libcrystfel/Makefile.in +++ b/libcrystfel/Makefile.in @@ -89,7 +89,8 @@ am_libcrystfel_la_OBJECTS = src/reflist.lo src/utils.lo src/cell.lo \ src/hdf5-file.lo src/beam-parameters.lo src/geometry.lo \ src/statistics.lo src/symmetry.lo src/stream.lo src/peaks.lo \ src/reflist-utils.lo src/filters.lo src/diffraction.lo \ - src/diffraction-gpu.lo src/cl-utils.lo + src/diffraction-gpu.lo src/cl-utils.lo src/render.lo \ + src/index.lo src/dirax.lo src/mosflm.lo src/reax.lo libcrystfel_la_OBJECTS = $(am_libcrystfel_la_OBJECTS) AM_V_lt = $(am__v_lt_$(V)) am__v_lt_ = $(am__v_lt_$(AM_DEFAULT_VERBOSITY)) @@ -351,7 +352,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/beam-parameters.c src/geometry.c src/statistics.c \ src/symmetry.c src/stream.c src/peaks.c \ src/reflist-utils.c src/filters.c src/diffraction.c \ - src/diffraction-gpu.c src/cl-utils.c + src/diffraction-gpu.c src/cl-utils.c src/render.c \ + src/index.c src/dirax.c src/mosflm.c src/reax.c libcrystfel_la_includedir = $(includedir)/crystfel/ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/diffraction-gpu.h \ @@ -361,6 +363,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/diffraction-gpu.h \ src/cl-utils.h src/filters.h src/list_tmp.h \ src/statistics.h src/utils.h src/detector.h \ src/geometry.h src/peaks.h src/stream.h \ + src/render.h src/index.h \ ../data/defs.h INCLUDES = "-I$(top_srcdir)/data" @@ -458,6 +461,11 @@ src/diffraction.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/diffraction-gpu.lo: src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/cl-utils.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/render.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/index.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/dirax.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/mosflm.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/reax.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) @@ -475,6 +483,8 @@ mostlyclean-compile: -rm -f src/diffraction-gpu.lo -rm -f src/diffraction.$(OBJEXT) -rm -f src/diffraction.lo + -rm -f src/dirax.$(OBJEXT) + -rm -f src/dirax.lo -rm -f src/filters.$(OBJEXT) -rm -f src/filters.lo -rm -f src/geometry.$(OBJEXT) @@ -483,12 +493,20 @@ mostlyclean-compile: -rm -f src/hdf5-file.lo -rm -f src/image.$(OBJEXT) -rm -f src/image.lo + -rm -f src/index.$(OBJEXT) + -rm -f src/index.lo + -rm -f src/mosflm.$(OBJEXT) + -rm -f src/mosflm.lo -rm -f src/peaks.$(OBJEXT) -rm -f src/peaks.lo + -rm -f src/reax.$(OBJEXT) + -rm -f src/reax.lo -rm -f src/reflist-utils.$(OBJEXT) -rm -f src/reflist-utils.lo -rm -f src/reflist.$(OBJEXT) -rm -f src/reflist.lo + -rm -f src/render.$(OBJEXT) + -rm -f src/render.lo -rm -f src/statistics.$(OBJEXT) -rm -f src/statistics.lo -rm -f src/stream.$(OBJEXT) @@ -509,13 +527,18 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/detector.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction-gpu.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/dirax.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/filters.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/geometry.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/hdf5-file.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/image.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/index.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/mosflm.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/peaks.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reax.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist-utils.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/statistics.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/stream.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/symmetry.Plo@am__quote@ diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c new file mode 100644 index 00000000..bb1cb808 --- /dev/null +++ b/libcrystfel/src/dirax.c @@ -0,0 +1,529 @@ +/* + * dirax.c + * + * Invoke the DirAx auto-indexing program + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if HAVE_FORKPTY_LINUX +#include +#elif HAVE_FORKPTY_BSD +#include +#endif + + +#include "image.h" +#include "dirax.h" +#include "utils.h" +#include "peaks.h" + + +#define DIRAX_VERBOSE 0 + +#define MAX_DIRAX_CELL_CANDIDATES (5) + + +typedef enum { + DIRAX_INPUT_NONE, + DIRAX_INPUT_LINE, + DIRAX_INPUT_PROMPT, + DIRAX_INPUT_ACL +} DirAxInputType; + + +struct dirax_data { + + /* DirAx auto-indexing low-level stuff */ + int pty; + pid_t pid; + char *rbuffer; + int rbufpos; + int rbuflen; + + /* DirAx auto-indexing high-level stuff */ + int step; + int finished_ok; + int read_cell; + int best_acl; + int best_acl_nh; + int acls_tried[MAX_CELL_CANDIDATES]; + int n_acls_tried; + +}; + + +static void dirax_parseline(const char *line, struct image *image, + struct dirax_data *dirax) +{ + int rf, i, di, acl, acl_nh; + float d; + + #if DIRAX_VERBOSE + char *copy; + + copy = strdup(line); + for ( i=0; iread_cell = 1; + image->candidate_cells[image->ncells] = cell_new(); + return; + } + i++; + } + + /* Parse unit cell vectors as appropriate */ + if ( dirax->read_cell == 1 ) { + /* First row of unit cell values */ + float ax, ay, az; + int r; + r = sscanf(line, "%f %f %f %f %f %f", + &d, &d, &d, &ax, &ay, &az); + if ( r != 6 ) { + ERROR("Couldn't understand cell line:\n"); + ERROR("'%s'\n", line); + dirax->read_cell = 0; + cell_free(image->candidate_cells[image->ncells]); + return; + } + cell_set_cartesian_a(image->candidate_cells[image->ncells], + ax*1e-10, ay*1e-10, az*1e-10); + dirax->read_cell++; + return; + } else if ( dirax->read_cell == 2 ) { + /* Second row of unit cell values */ + float bx, by, bz; + int r; + r = sscanf(line, "%f %f %f %f %f %f", + &d, &d, &d, &bx, &by, &bz); + if ( r != 6 ) { + ERROR("Couldn't understand cell line:\n"); + ERROR("'%s'\n", line); + dirax->read_cell = 0; + cell_free(image->candidate_cells[image->ncells]); + return; + } + cell_set_cartesian_b(image->candidate_cells[image->ncells], + bx*1e-10, by*1e-10, bz*1e-10); + dirax->read_cell++; + return; + } else if ( dirax->read_cell == 3 ) { + /* Third row of unit cell values */ + float cx, cy, cz; + int r; + r = sscanf(line, "%f %f %f %f %f %f", + &d, &d, &d, &cx, &cy, &cz); + if ( r != 6 ) { + ERROR("Couldn't understand cell line:\n"); + ERROR("'%s'\n", line); + dirax->read_cell = 0; + cell_free(image->candidate_cells[image->ncells]); + return; + } + cell_set_cartesian_c(image->candidate_cells[image->ncells++], + cx*1e-10, cy*1e-10, cz*1e-10); + dirax->read_cell = 0; + return; + } + + dirax->read_cell = 0; + + if ( sscanf(line, "%i %i %f %f %f %f %f %f %i", &acl, &acl_nh, + &d, &d, &d, &d, &d, &d, &di) == 9 ) { + if ( acl_nh > dirax->best_acl_nh ) { + + int i, found = 0; + + for ( i=0; in_acls_tried; i++ ) { + if ( dirax->acls_tried[i] == acl ) found = 1; + } + + if ( !found ) { + dirax->best_acl_nh = acl_nh; + dirax->best_acl = acl; + } + + } + } +} + + +static void dirax_sendline(const char *line, struct dirax_data *dirax) +{ + #if DIRAX_VERBOSE + char *copy; + int i; + + copy = strdup(line); + for ( i=0; ipty, line, strlen(line)); +} + + +static void dirax_send_next(struct image *image, struct dirax_data *dirax) +{ + char tmp[32]; + + switch ( dirax->step ) { + + case 1 : + dirax_sendline("\\echo off\n", dirax); + break; + + case 2 : + snprintf(tmp, 31, "read xfel-%i.drx\n", image->id); + dirax_sendline(tmp, dirax); + break; + + case 3 : + dirax_sendline("dmax 1000\n", dirax); + break; + + case 4 : + dirax_sendline("indexfit 2\n", dirax); + break; + + case 5 : + dirax_sendline("levelfit 1000\n", dirax); + break; + + case 6 : + dirax_sendline("go\n", dirax); + dirax->finished_ok = 1; + break; + + case 7 : + dirax_sendline("acl\n", dirax); + break; + + case 8 : + if ( dirax->best_acl_nh == 0 ) { + /* At this point, DirAx is presenting its ACL prompt + * and waiting for a single number. Use an extra + * newline to choose automatic ACL selection before + * exiting. */ + dirax_sendline("\nexit\n", dirax); + break; + } + snprintf(tmp, 31, "%i\n", dirax->best_acl); + dirax->acls_tried[dirax->n_acls_tried++] = dirax->best_acl; + dirax_sendline(tmp, dirax); + break; + + case 9 : + dirax_sendline("cell\n", dirax); + break; + + case 10 : + if ( dirax->n_acls_tried == MAX_DIRAX_CELL_CANDIDATES ) { + dirax_sendline("exit\n", dirax); + } else { + /* Go back round for another cell */ + dirax->best_acl_nh = 0; + dirax->step = 7; + dirax_sendline("acl\n", dirax); + } + break; + + default: + dirax_sendline("exit\n", dirax); + return; + + } + + dirax->step++; +} + + +static int dirax_readable(struct image *image, struct dirax_data *dirax) +{ + int rval; + int no_string = 0; + + rval = read(dirax->pty, dirax->rbuffer+dirax->rbufpos, + dirax->rbuflen-dirax->rbufpos); + + if ( (rval == -1) || (rval == 0) ) return 1; + + dirax->rbufpos += rval; + assert(dirax->rbufpos <= dirax->rbuflen); + + while ( (!no_string) && (dirax->rbufpos > 0) ) { + + int i; + int block_ready = 0; + DirAxInputType type = DIRAX_INPUT_NONE; + + /* See if there's a full line in the buffer yet */ + for ( i=0; irbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + + /* Is there a prompt in the buffer? */ + if ( (i+7 <= dirax->rbufpos) + && (!strncmp(dirax->rbuffer+i, "Dirax> ", 7)) ) { + block_ready = 1; + type = DIRAX_INPUT_PROMPT; + break; + } + + /* How about an ACL prompt? */ + if ( (i+10 <= dirax->rbufpos) + && (!strncmp(dirax->rbuffer+i, "acl/auto [", 10)) ) { + block_ready = 1; + type = DIRAX_INPUT_ACL; + break; + } + + if ( (dirax->rbuffer[i] == '\r') + && (dirax->rbuffer[i+1] == '\n') ) { + block_ready = 1; + type = DIRAX_INPUT_LINE; + break; + } + + } + + if ( block_ready ) { + + unsigned int new_rbuflen; + unsigned int endbit_length; + char *block_buffer = NULL; + + switch ( type ) { + + case DIRAX_INPUT_LINE : + + block_buffer = malloc(i+1); + memcpy(block_buffer, dirax->rbuffer, i); + block_buffer[i] = '\0'; + + if ( block_buffer[0] == '\r' ) { + memmove(block_buffer, block_buffer+1, i); + } + + dirax_parseline(block_buffer, image, dirax); + free(block_buffer); + endbit_length = i+2; + break; + + case DIRAX_INPUT_PROMPT : + + dirax_send_next(image, dirax); + endbit_length = i+7; + break; + + case DIRAX_INPUT_ACL : + + dirax_send_next(image,dirax ); + endbit_length = i+10; + break; + + default : + + /* Obviously, this never happens :) */ + ERROR("Unrecognised DirAx input mode! " + "I don't know how to understand DirAx\n"); + return 1; + + } + + /* Now the block's been parsed, it should be + * forgotten about */ + memmove(dirax->rbuffer, + dirax->rbuffer + endbit_length, + dirax->rbuflen - endbit_length); + + /* Subtract the number of bytes removed */ + dirax->rbufpos = dirax->rbufpos + - endbit_length; + new_rbuflen = dirax->rbuflen - endbit_length; + if ( new_rbuflen == 0 ) new_rbuflen = 256; + dirax->rbuffer = realloc(dirax->rbuffer, + new_rbuflen); + dirax->rbuflen = new_rbuflen; + + } else { + + if ( dirax->rbufpos==dirax->rbuflen ) { + + /* More buffer space is needed */ + dirax->rbuffer = realloc( + dirax->rbuffer, + dirax->rbuflen + 256); + dirax->rbuflen = dirax->rbuflen + + 256; + /* The new space gets used at the next + * read, shortly... */ + + } + no_string = 1; + + } + + } + + return 0; +} + + +static void write_drx(struct image *image) +{ + FILE *fh; + int i; + char filename[1024]; + + snprintf(filename, 1023, "xfel-%i.drx", image->id); + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + return; + } + fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. */ + + for ( i=0; ifeatures); i++ ) { + + struct imagefeature *f; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + fprintf(fh, "%10f %10f %10f %8f\n", + f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0); + + } + fclose(fh); +} + + +void run_dirax(struct image *image) +{ + unsigned int opts; + int status; + int rval; + struct dirax_data *dirax; + + write_drx(image); + + dirax = malloc(sizeof(struct dirax_data)); + if ( dirax == NULL ) { + ERROR("Couldn't allocate memory for DirAx data.\n"); + return; + } + + dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL); + if ( dirax->pid == -1 ) { + ERROR("Failed to fork for DirAx\n"); + return; + } + if ( dirax->pid == 0 ) { + + /* Child process: invoke DirAx */ + struct termios t; + + /* Turn echo off */ + tcgetattr(STDIN_FILENO, &t); + t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); + tcsetattr(STDIN_FILENO, TCSANOW, &t); + + execlp("dirax", "", (char *)NULL); + ERROR("Failed to invoke DirAx.\n"); + _exit(0); + + } + + dirax->rbuffer = malloc(256); + dirax->rbuflen = 256; + dirax->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(dirax->pty, F_GETFL); + fcntl(dirax->pty, F_SETFL, opts | O_NONBLOCK); + + dirax->step = 1; /* This starts the "initialisation" procedure */ + dirax->finished_ok = 0; + dirax->read_cell = 0; + dirax->n_acls_tried = 0; + dirax->best_acl_nh = 0; + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(dirax->pty, &fds); + + tv.tv_sec = 10; + tv.tv_usec = 0; + + sval = select(dirax->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + int err = errno; + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + } else if ( sval != 0 ) { + rval = dirax_readable(image, dirax); + } else { + ERROR("No response from DirAx..\n"); + rval = 1; + } + + } while ( !rval ); + + close(dirax->pty); + free(dirax->rbuffer); + waitpid(dirax->pid, &status, 0); + + if ( dirax->finished_ok == 0 ) { + ERROR("DirAx doesn't seem to be working properly.\n"); + } + + free(dirax); +} diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h new file mode 100644 index 00000000..8c429710 --- /dev/null +++ b/libcrystfel/src/dirax.h @@ -0,0 +1,26 @@ +/* + * dirax.h + * + * Invoke the DirAx auto-indexing program + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef DIRAX_H +#define DIRAX_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "utils.h" + + +extern void run_dirax(struct image *image); + + +#endif /* DIRAX_H */ diff --git a/libcrystfel/src/index-priv.h b/libcrystfel/src/index-priv.h new file mode 100644 index 00000000..3d8c8a22 --- /dev/null +++ b/libcrystfel/src/index-priv.h @@ -0,0 +1,29 @@ +/* + * index-priv.h + * + * Indexing private data + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef INDEXPRIV_H +#define INDEXPRIV_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include "index.h" + +struct _indexingprivate +{ + IndexingMethod indm; +}; + + +#endif /* INDEXPRIV_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c new file mode 100644 index 00000000..d5e76c50 --- /dev/null +++ b/libcrystfel/src/index.c @@ -0,0 +1,274 @@ +/* + * index.c + * + * Perform indexing (somehow) + * + * (c) 2006-2011 Thomas White + * (c) 2010 Richard Kirian + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" +#include "peaks.h" +#include "dirax.h" +#include "mosflm.h" +#include "detector.h" +#include "index.h" +#include "index-priv.h" +#include "reax.h" +#include "geometry.h" + + +/* Base class constructor for unspecialised indexing private data */ +IndexingPrivate *indexing_private(IndexingMethod indm) +{ + struct _indexingprivate *priv; + priv = calloc(1, sizeof(struct _indexingprivate)); + priv->indm = indm; + return priv; +} + + +static const char *maybes(int n) +{ + if ( n == 1 ) return ""; + return "s"; +} + + +IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + double nominal_photon_energy) +{ + int n; + int nm = 0; + IndexingPrivate **iprivs; + + while ( indm[nm] != INDEXING_NONE ) nm++; + STATUS("Preparing %i indexing method%s.\n", nm, maybes(nm)); + iprivs = malloc((nm+1) * sizeof(IndexingPrivate *)); + + for ( n=0; nindm ) { + case INDEXING_NONE : + free(priv[n]); + break; + case INDEXING_DIRAX : + free(priv[n]); + break; + case INDEXING_MOSFLM : + free(priv[n]); + break; + case INDEXING_REAX : + reax_cleanup(priv[n]); + break; + } + + n++; + + } +} + + +void map_all_peaks(struct image *image) +{ + int i, n; + + n = image_feature_count(image->features); + + /* Map positions to 3D */ + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); + f->rx = r.u; f->ry = r.v; f->rz = r.w; + + } +} + + +void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, + int cellr, int verbose, IndexingPrivate **ipriv, + int config_insane) +{ + int i; + int n = 0; + + if ( indm == NULL ) return; + + map_all_peaks(image); + image->indexed_cell = NULL; + + while ( indm[n] != INDEXING_NONE ) { + + image->ncells = 0; + + /* Index as appropriate */ + switch ( indm[n] ) { + case INDEXING_NONE : + return; + case INDEXING_DIRAX : + run_dirax(image); + break; + case INDEXING_MOSFLM : + run_mosflm(image, cell); + break; + case INDEXING_REAX : + reax_index(ipriv[n], image, cell); + break; + } + if ( image->ncells == 0 ) { + n++; + continue; + } + + for ( i=0; incells; i++ ) { + + UnitCell *new_cell = NULL; + UnitCell *cand = image->candidate_cells[i]; + + if ( verbose ) { + STATUS("--------------------\n"); + STATUS("Candidate cell %i (before matching):\n", + i); + cell_print(image->candidate_cells[i]); + STATUS("--------------------\n"); + } + + /* Match or reduce the cell as appropriate */ + switch ( cellr ) { + case CELLR_NONE : + new_cell = cell_new_from_cell(cand); + break; + case CELLR_REDUCE : + new_cell = match_cell(cand, cell, verbose, 1); + break; + case CELLR_COMPARE : + new_cell = match_cell(cand, cell, verbose, 0); + break; + case CELLR_COMPARE_AB : + new_cell = match_cell_ab(cand, cell); + break; + } + + /* No cell? Move on to the next candidate */ + if ( new_cell == NULL ) continue; + + /* Sanity check */ + image->reflections = find_intersections(image, new_cell); + image->indexed_cell = new_cell; + + if ( !config_insane && !peak_sanity_check(image) ) { + cell_free(new_cell); + image->indexed_cell = NULL; + continue; + } + + goto done; /* Success */ + + } + + for ( i=0; incells; i++ ) { + cell_free(image->candidate_cells[i]); + image->candidate_cells[i] = NULL; + } + + /* Move on to the next indexing method */ + n++; + + } + +done: + for ( i=0; incells; i++ ) { + /* May free(NULL) if all algorithms were tried and no success */ + cell_free(image->candidate_cells[i]); + } +} + + +IndexingMethod *build_indexer_list(const char *str, int *need_cell) +{ + int n, i; + char **methods; + IndexingMethod *list; + int tmp; + + if ( need_cell == NULL ) need_cell = &tmp; + *need_cell = 0; + + n = assplode(str, ",", &methods, ASSPLODE_NONE); + list = malloc((n+1)*sizeof(IndexingMethod)); + + for ( i=0; i + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef INDEX_H +#define INDEX_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include "cell.h" +#include "image.h" +#include "detector.h" + + +/* Indexing methods */ +typedef enum { + INDEXING_NONE, + INDEXING_DIRAX, + INDEXING_MOSFLM, + INDEXING_REAX, +} IndexingMethod; + + +/* Cell reduction methods */ +enum { + CELLR_NONE, + CELLR_REDUCE, + CELLR_COMPARE, + CELLR_COMPARE_AB, +}; + + +typedef struct _indexingprivate IndexingPrivate; + +extern IndexingPrivate *indexing_private(IndexingMethod indm); + +extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, + const char *filename, + struct detector *det, + double nominal_photon_energy); + +extern void map_all_peaks(struct image *image); + +extern void index_pattern(struct image *image, UnitCell *cell, + IndexingMethod *indm, int cellr, int verbose, + IndexingPrivate **priv, int config_insane); + +extern void cleanup_indexing(IndexingPrivate **priv); + +extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); + +#endif /* INDEX_H */ diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c new file mode 100644 index 00000000..5df5af21 --- /dev/null +++ b/libcrystfel/src/mosflm.c @@ -0,0 +1,609 @@ +/* + * mosflm.c + * + * Invoke the DPS auto-indexing algorithm through MOSFLM + * + * (c) 2010 Richard Kirian + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +/* TODO: + * + * Properly read the newmat file (don't use fscanf-- spaces between numers + * are not guaranteed) + * + * "Success" is indicated by existence of NEWMAT file written by mosflm. + * Better to interact with mosflm directly in order to somehow verify success. + * + * Investigate how these keywords affect mosflms behavior: + * + * MOSAICITY + * DISPERSION + * DIVERGENCE + * POLARISATION + * POSTREF BEAM + * POSTREF USEBEAM OFF + * PREREFINE ON + * EXTRA ON + * POSTREF ON + * + * These did not seem to affect the results by my (Rick's) experience, probably + * because they are only used conjunction with image intensity data, but it's + * worth another look at the documentation. + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if HAVE_FORKPTY_LINUX +#include +#elif HAVE_FORKPTY_BSD +#include +#endif + + +#include "image.h" +#include "mosflm.h" +#include "utils.h" +#include "peaks.h" + + +#define MOSFLM_VERBOSE 0 + + +typedef enum { + MOSFLM_INPUT_NONE, + MOSFLM_INPUT_LINE, + MOSFLM_INPUT_PROMPT +} MOSFLMInputType; + + +struct mosflm_data { + + /* MOSFLM auto-indexing low-level stuff */ + int pty; + pid_t pid; + char *rbuffer; + int rbufpos; + int rbuflen; + + /* MOSFLM high-level stuff */ + char newmatfile[128]; + char imagefile[128]; + char sptfile[128]; + int step; + int finished_ok; + UnitCell *target_cell; + +}; + + +static void mosflm_parseline(const char *line, struct image *image, + struct mosflm_data *dirax) +{ + #if MOSFLM_VERBOSE + char *copy; + int i; + + copy = strdup(line); + for ( i=0; ilambda; + + image->candidate_cells[0] = cell_new(); + + cell_set_reciprocal(image->candidate_cells[0], + asz*c, asy*c, asx*c, + bsz*c, bsy*c, bsx*c, + csz*c, csy*c, csx*c); + + image->ncells = 1; + + return 0; +} + + +/* Need to sort mosflm peaks by intensity... */ +struct sptline { + double x; /* x coordinate of peak */ + double y; /* y coordinate of peak */ + double h; /* height of peak */ + double s; /* sigma of peak */ +}; + + +static int compare_vals(const void *ap, const void *bp) +{ + const struct sptline a = *(struct sptline *)ap; + const struct sptline b = *(struct sptline *)bp; + + if ( a.h < b.h ) return 1; + if ( a.h > b.h ) return -1; + return 0; +} + + +/* Write .spt file for mosflm */ +static void write_spt(struct image *image, const char *filename) +{ + FILE *fh; + int i; + double fclen = 67.8; /* fake camera length in mm */ + double fpix = 0.075; /* fake pixel size in mm */ + double pix; + double height = 100.0; + double sigma = 1.0; + int nPeaks = image_feature_count(image->features); + struct sptline *sptlines; + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + return; + } + + fprintf(fh, "%10d %10d %10.8f %10.6f %10.6f\n", 1, 1, fpix, 1.0, 0.0); + fprintf(fh, "%10d %10d\n", 1, 1); + fprintf(fh, "%10.5f %10.5f\n", 0.0, 0.0); + + sptlines = malloc(sizeof(struct sptline)*nPeaks); + + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + p = find_panel(image->det, f->fs, f->ss); + if ( p == NULL ) continue; + + pix = 1000.0/p->res; /* pixel size in mm */ + height = f->intensity; + + xs = (f->fs-p->min_fs)*p->fsx + (f->ss-p->min_ss)*p->ssx; + ys = (f->ss-p->min_fs)*p->fsy + (f->ss-p->min_ss)*p->ssy; + rx = xs + p->cnx; + ry = ys + p->cny; + + sptlines[i].x = ry*pix*fclen/p->clen/1000.0; + sptlines[i].y = -rx*pix*fclen/p->clen/1000.0; + sptlines[i].h = height; + sptlines[i].s = sigma/1000.0; + + } + + qsort(sptlines, nPeaks, sizeof(struct sptline), compare_vals); + + for ( i=0; ipty, line, strlen(line)); +} + + +static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) +{ + char tmp[256]; + char symm[32]; + const char *sg; + double wavelength; + double a, b, c, alpha, beta, gamma; + int i, j; + + switch ( mosflm->step ) { + + case 1 : + mosflm_sendline("DETECTOR ROTATION HORIZONTAL" + " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" + " RECTANGULAR\n", mosflm); + break; + + case 2 : + if ( mosflm->target_cell != NULL ) { + cell_get_parameters(mosflm->target_cell, &a, &b, &c, + &alpha, &beta, &gamma); + snprintf(tmp, 255, + "CELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n", + a*1e10, b*1e10, c*1e10, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + mosflm_sendline(tmp, mosflm); + } else { + mosflm_sendline("# Do nothing\n", mosflm); + } + break; + + case 3 : + if ( mosflm->target_cell != NULL ) { + sg = cell_get_spacegroup(mosflm->target_cell); + /* Remove white space from space group */ + j = 0; + for ( i=0; ilambda*1e10; + snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength); + mosflm_sendline(tmp, mosflm); + break; + + case 7 : + snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile); + mosflm_sendline(tmp, mosflm); + break; + + case 8 : + snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile); + mosflm_sendline(tmp, mosflm); + break; + + case 9 : + snprintf(tmp, 255, "AUTOINDEX DPS FILE %s" + " IMAGE 1 MAXCELL 1000 REFINE\n", + mosflm->sptfile); + + /* "This option is only available if you e-mail Andrew Leslie + * and ask for it." - MOSFLM + * snprintf(tmp, 255, "AUTOINDEX NODISPLAY IMAGE 1 FILE %s\n", + * mosflm->sptfile); */ + mosflm_sendline(tmp, mosflm); + break; + + case 10 : + mosflm_sendline("GO\n", mosflm); + mosflm->finished_ok = 1; + break; + + default: + mosflm_sendline("exit\n", mosflm); + return; + + } + + mosflm->step++; +} + + +static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) +{ + int rval; + int no_string = 0; + + rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos, + mosflm->rbuflen-mosflm->rbufpos); + + if ( (rval == -1) || (rval == 0) ) return 1; + + mosflm->rbufpos += rval; + assert(mosflm->rbufpos <= mosflm->rbuflen); + + while ( (!no_string) && (mosflm->rbufpos > 0) ) { + + int i; + int block_ready = 0; + MOSFLMInputType type = MOSFLM_INPUT_NONE; + + /* See if there's a full line in the buffer yet */ + for ( i=0; irbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + + /* Is there a prompt in the buffer? */ + if ( (i+10 <= mosflm->rbufpos) + && (!strncmp(mosflm->rbuffer+i, "MOSFLM => ", 10)) ) { + block_ready = 1; + type = MOSFLM_INPUT_PROMPT; + break; + } + + if ( (mosflm->rbuffer[i] == '\r') + && (mosflm->rbuffer[i+1] == '\n') ) { + block_ready = 1; + type = MOSFLM_INPUT_LINE; + break; + } + + } + + if ( block_ready ) { + + unsigned int new_rbuflen; + unsigned int endbit_length; + char *block_buffer = NULL; + + switch ( type ) { + + case MOSFLM_INPUT_LINE : + + block_buffer = malloc(i+1); + memcpy(block_buffer, mosflm->rbuffer, i); + block_buffer[i] = '\0'; + + if ( block_buffer[0] == '\r' ) { + memmove(block_buffer, block_buffer+1, i); + } + + mosflm_parseline(block_buffer, image, mosflm); + free(block_buffer); + endbit_length = i+2; + break; + + case MOSFLM_INPUT_PROMPT : + + mosflm_send_next(image, mosflm); + endbit_length = i+7; + break; + + default : + + /* Obviously, this never happens :) */ + ERROR("Unrecognised MOSFLM input mode! " + "I don't know how to understand MOSFLM\n"); + return 1; + + } + + /* Now the block's been parsed, it should be + * forgotten about */ + memmove(mosflm->rbuffer, + mosflm->rbuffer + endbit_length, + mosflm->rbuflen - endbit_length); + + /* Subtract the number of bytes removed */ + mosflm->rbufpos = mosflm->rbufpos + - endbit_length; + new_rbuflen = mosflm->rbuflen - endbit_length; + if ( new_rbuflen == 0 ) new_rbuflen = 256; + mosflm->rbuffer = realloc(mosflm->rbuffer, + new_rbuflen); + mosflm->rbuflen = new_rbuflen; + + } else { + + if ( mosflm->rbufpos==mosflm->rbuflen ) { + + /* More buffer space is needed */ + mosflm->rbuffer = realloc( + mosflm->rbuffer, + mosflm->rbuflen + 256); + mosflm->rbuflen = mosflm->rbuflen + 256; + /* The new space gets used at the next + * read, shortly... */ + + } + no_string = 1; + + } + + } + + return 0; +} + + +void run_mosflm(struct image *image, UnitCell *cell) +{ + struct mosflm_data *mosflm; + unsigned int opts; + int status; + int rval; + + mosflm = malloc(sizeof(struct mosflm_data)); + if ( mosflm == NULL ) { + ERROR("Couldn't allocate memory for MOSFLM data.\n"); + return; + } + + mosflm->target_cell = cell; + + snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id); + write_img(image, mosflm->imagefile); /* Dummy image */ + + snprintf(mosflm->sptfile, 127, "xfel-%i_001.spt", image->id); + write_spt(image, mosflm->sptfile); + + snprintf(mosflm->newmatfile, 127, "xfel-%i.newmat", image->id); + remove(mosflm->newmatfile); + + mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL); + if ( mosflm->pid == -1 ) { + ERROR("Failed to fork for MOSFLM\n"); + free(mosflm); + return; + } + if ( mosflm->pid == 0 ) { + + /* Child process: invoke MOSFLM */ + struct termios t; + + /* Turn echo off */ + tcgetattr(STDIN_FILENO, &t); + t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); + tcsetattr(STDIN_FILENO, TCSANOW, &t); + + execlp("ipmosflm", "", (char *)NULL); + ERROR("Failed to invoke MOSFLM.\n"); + _exit(0); + + } + + mosflm->rbuffer = malloc(256); + mosflm->rbuflen = 256; + mosflm->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(mosflm->pty, F_GETFL); + fcntl(mosflm->pty, F_SETFL, opts | O_NONBLOCK); + + mosflm->step = 1; /* This starts the "initialisation" procedure */ + mosflm->finished_ok = 0; + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(mosflm->pty, &fds); + + tv.tv_sec = 10; + tv.tv_usec = 0; + + sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + int err = errno; + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + } else if ( sval != 0 ) { + rval = mosflm_readable(image, mosflm); + } else { + ERROR("No response from MOSFLM..\n"); + rval = 1; + } + + } while ( !rval ); + + close(mosflm->pty); + free(mosflm->rbuffer); + waitpid(mosflm->pid, &status, 0); + + if ( mosflm->finished_ok == 0 ) { + ERROR("MOSFLM doesn't seem to be working properly.\n"); + } else { + /* Read the mosflm NEWMAT file and get cell if found */ + read_newmat(mosflm->newmatfile, image); + } + + free(mosflm); +} diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h new file mode 100644 index 00000000..82b80a44 --- /dev/null +++ b/libcrystfel/src/mosflm.h @@ -0,0 +1,27 @@ +/* + * mosflm.h + * + * Invoke the DPS auto-indexing algorithm through MOSFLM + * + * (c) 2010 Richard Kirian + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef MOSFLM_H +#define MOSFLM_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "utils.h" + + +extern void run_mosflm(struct image *image, UnitCell *cell); + + +#endif /* MOSFLM_H */ diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c new file mode 100644 index 00000000..ff1ea582 --- /dev/null +++ b/libcrystfel/src/reax.c @@ -0,0 +1,728 @@ +/* + * reax.c + * + * A new auto-indexer + * + * (c) 2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" +#include "peaks.h" +#include "cell.h" +#include "index.h" +#include "index-priv.h" + + +struct dvec +{ + double x; + double y; + double z; + double th; + double ph; +}; + + +struct reax_private +{ + IndexingPrivate base; + struct dvec *directions; + int n_dir; + double angular_inc; + + double *fft_in; + fftw_complex *fft_out; + fftw_plan plan; + int nel; + + fftw_complex *r_fft_in; + fftw_complex *r_fft_out; + fftw_plan r_plan; + int ch; + int cw; +}; + + +static double check_dir(struct dvec *dir, ImageFeatureList *flist, + int nel, double pmax, double *fft_in, + fftw_complex *fft_out, fftw_plan plan, + int smin, int smax, + const char *rg, struct detector *det) +{ + int n, i; + double tot; + + for ( i=0; ifs, f->ss); + assert(p != NULL); + + if ( p->rigid_group != rg ) continue; + + } + + val = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z; + + idx = nel/2 + nel*val/(2.0*pmax); + fft_in[idx]++; + + } + + fftw_execute_dft_r2c(plan, fft_in, fft_out); + + tot = 0.0; + for ( i=smin; i<=smax; i++ ) { + double re, im; + re = fft_out[i][0]; + im = fft_out[i][1]; + tot += sqrt(re*re + im*im); + } + + return tot; +} + + +/* Refine a direct space vector. From Clegg (1984) */ +static double iterate_refine_vector(double *x, double *y, double *z, + ImageFeatureList *flist) +{ + int fi, n, err; + gsl_matrix *C; + gsl_vector *A; + gsl_vector *t; + gsl_matrix *s_vec; + gsl_vector *s_val; + double dtmax; + + A = gsl_vector_calloc(3); + C = gsl_matrix_calloc(3, 3); + + n = image_feature_count(flist); + fesetround(1); + for ( fi=0; firx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */ + kn = nearbyint(kno); + if ( kn - kno > 0.3 ) continue; + + xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz; + + for ( i=0; i<3; i++ ) { + + val = gsl_vector_get(A, i); + gsl_vector_set(A, i, val+xv[i]*kn); + + for ( j=0; j<3; j++ ) { + val = gsl_matrix_get(C, i, j); + gsl_matrix_set(C, i, j, val+xv[i]*xv[j]); + } + + } + + } + + s_val = gsl_vector_calloc(3); + s_vec = gsl_matrix_calloc(3, 3); + err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val); + if ( err ) { + ERROR("SVD failed: %s\n", gsl_strerror(err)); + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + gsl_matrix_free(C); + gsl_vector_free(A); + return 0.0; + } + + t = gsl_vector_calloc(3); + err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t); + if ( err ) { + ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + gsl_matrix_free(C); + gsl_vector_free(A); + gsl_vector_free(t); + return 0.0; + } + + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + + dtmax = fabs(*x - gsl_vector_get(t, 0)); + dtmax += fabs(*y - gsl_vector_get(t, 1)); + dtmax += fabs(*z - gsl_vector_get(t, 2)); + + *x = gsl_vector_get(t, 0); + *y = gsl_vector_get(t, 1); + *z = gsl_vector_get(t, 2); + + gsl_matrix_free(C); + gsl_vector_free(A); + + return dtmax; +} + + +static void refine_cell(struct image *image, UnitCell *cell, + ImageFeatureList *flist) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + int i; + double sm; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + i = 0; + do { + + sm = iterate_refine_vector(&ax, &ay, &az, flist); + sm += iterate_refine_vector(&bx, &by, &bz, flist); + sm += iterate_refine_vector(&cx, &cy, &cz, flist); + i++; + + } while ( (sm > 0.001e-9) && (i<10) ); + + cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz); + + if ( i == 10 ) { + cell_free(image->indexed_cell); + image->indexed_cell = NULL; + } +} + + +static void fine_search(struct reax_private *p, ImageFeatureList *flist, + int nel, double pmax, double *fft_in, + fftw_complex *fft_out, fftw_plan plan, + int smin, int smax, double th_cen, double ph_cen, + double *x, double *y, double *z) +{ + double fom = 0.0; + double th, ph; + double inc; + struct dvec dir; + int i, s; + double max, modv; + + inc = p->angular_inc / 100.0; + + for ( th=th_cen-p->angular_inc; th<=th_cen+p->angular_inc; th+=inc ) { + for ( ph=ph_cen-p->angular_inc; ph<=ph_cen+p->angular_inc; ph+=inc ) { + + double new_fom; + + dir.x = cos(ph) * sin(th); + dir.y = sin(ph) * sin(th); + dir.z = cos(th); + + new_fom = check_dir(&dir, flist, nel, pmax, fft_in, fft_out, + plan, smin, smax, NULL, NULL); + + if ( new_fom > fom ) { + fom = new_fom; + *x = dir.x; *y = dir.y; *z = dir.z; + } + + } + } + + s = -1; + max = 0.0; + for ( i=smin; i<=smax; i++ ) { + + double re, im, m; + + re = fft_out[i][0]; + im = fft_out[i][1]; + m = sqrt(re*re + im*im); + if ( m > max ) { + max = m; + s = i; + } + + } + assert(s>0); + + modv = (double)s / (2.0*pmax); + *x *= modv; *y *= modv; *z *= modv; +} + + +static double get_model_phase(double x, double y, double z, ImageFeatureList *f, + int nel, double pmax, double *fft_in, + fftw_complex *fft_out, fftw_plan plan, + int smin, int smax, const char *rg, + struct detector *det) +{ + struct dvec dir; + int s, i; + double max; + double re, im; + + dir.x = x; dir.y = y; dir.z = z; + + check_dir(&dir, f, nel, pmax, fft_in,fft_out, plan, smin, smax, + rg, det); + + s = -1; + max = 0.0; + for ( i=smin; i<=smax; i++ ) { + + double re, im, m; + + re = fft_out[i][0]; + im = fft_out[i][1]; + m = sqrt(re*re + im*im); + if ( m > max ) { + max = m; + s = i; + } + + } + + re = fft_out[s][0]; + im = fft_out[s][1]; + + return atan2(im, re); +} + + +static void refine_rigid_group(struct image *image, UnitCell *cell, + const char *rg, int nel, double pmax, + double *fft_in, fftw_complex *fft_out, + fftw_plan plan, int smin, int smax, + struct detector *det, struct reax_private *pr) +{ + double ax, ay, az, ma; + double bx, by, bz, mb; + double cx, cy, cz, mc; + double pha, phb, phc; + struct panel *p; + int i, j; + fftw_complex *r_fft_in; + fftw_complex *r_fft_out; + double m2m; + signed int aix, aiy; + signed int bix, biy; + signed int cix, ciy; + double max; + int max_i, max_j; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + ma = modulus(ax, ay, az); + mb = modulus(bx, by, bz); + mc = modulus(cx, cy, cz); + + pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, nel, pmax, + fft_in, fft_out, plan, smin, smax, rg, det); + phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, nel, pmax, + fft_in, fft_out, plan, smin, smax, rg, det); + phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, nel, pmax, + fft_in, fft_out, plan, smin, smax, rg, det); + + for ( i=0; in_panels; i++ ) { + if ( det->panels[i].rigid_group == rg ) { + p = &det->panels[i]; + break; + } + } + + r_fft_in = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); + r_fft_out = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); + for ( i=0; icw; i++ ) { + for ( j=0; jch; j++ ) { + r_fft_in[i+pr->cw*j][0] = 0.0; + r_fft_in[i+pr->cw*j][1] = 0.0; + } + } + + ma = modulus(ax, ay, 0.0); + mb = modulus(bx, by, 0.0); + mc = modulus(cx, cy, 0.0); + m2m = ma; + if ( mb > m2m ) m2m = mb; + if ( mc > m2m ) m2m = mc; + + aix = (pr->cw/2)*ax/m2m; aiy = (pr->ch/2)*ay/m2m; + bix = (pr->cw/2)*bx/m2m; biy = (pr->ch/2)*by/m2m; + cix = (pr->cw/2)*cx/m2m; ciy = (pr->ch/2)*cy/m2m; + + if ( aix < 0 ) aix += pr->cw/2; + if ( bix < 0 ) bix += pr->cw/2; + if ( cix < 0 ) cix += pr->cw/2; + + if ( aiy < 0 ) aiy += pr->ch/2; + if ( biy < 0 ) biy += pr->ch/2; + if ( ciy < 0 ) ciy += pr->ch/2; + + r_fft_in[aix + pr->cw*aiy][0] = cos(pha); + r_fft_in[aix + pr->cw*aiy][1] = sin(pha); + r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][0] = cos(pha); + r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][1] = -sin(pha); + + r_fft_in[bix + pr->cw*biy][0] = cos(phb); + r_fft_in[bix + pr->cw*biy][1] = sin(phb); + r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][0] = cos(phb); + r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][1] = -sin(phb); + + r_fft_in[cix + pr->cw*ciy][0] = cos(phc); + r_fft_in[cix + pr->cw*ciy][1] = sin(phc); + r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][0] = cos(phc); + r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][1] = -sin(phc); + + const int tidx = 1; + r_fft_in[tidx][0] = 1.0; + r_fft_in[tidx][1] = 0.0; + +// STATUS("%i %i\n", aix, aiy); +// STATUS("%i %i\n", bix, biy); +// STATUS("%i %i\n", cix, ciy); + + fftw_execute_dft(pr->r_plan, r_fft_in, r_fft_out); + +// max = 0.0; +// FILE *fh = fopen("centering.dat", "w"); +// for ( i=0; icw; i++ ) { +// for ( j=0; jch; j++ ) { +// +// double re, im, am, ph; +// +// re = r_fft_out[i + pr->cw*j][0]; +// im = r_fft_out[i + pr->cw*j][1]; +// am = sqrt(re*re + im*im); +// ph = atan2(im, re); +// +// if ( am > max ) { +// max = am; +// max_i = i; +// max_j = j; +// } +// +// fprintf(fh, "%f ", am); +// +// } +// fprintf(fh, "\n"); +// } +// STATUS("Max at %i, %i\n", max_i, max_j); +// fclose(fh); +// exit(1); + +// STATUS("Offsets for '%s': %.2f, %.2f pixels\n", rg, dx, dy); +} + + +static void refine_all_rigid_groups(struct image *image, UnitCell *cell, + int nel, double pmax, + double *fft_in, fftw_complex *fft_out, + fftw_plan plan, int smin, int smax, + struct detector *det, + struct reax_private *p) +{ + int i; + + for ( i=0; idet->num_rigid_groups; i++ ) { + refine_rigid_group(image, cell, image->det->rigid_groups[i], + nel, pmax, fft_in, fft_out, plan, smin, smax, + det, p); + } +} + + +void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) +{ + int i; + struct reax_private *p; + double fom; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + double mod_a, mod_b, mod_c; + double al, be, ga; + double th, ph; + double *fft_in; + fftw_complex *fft_out; + int smin, smax; + double amin, amax; + double bmin, bmax; + double cmin, cmax; + double pmax; + int n; + const double ltol = 5.0; /* Direct space axis length + * tolerance in percent */ + const double angtol = deg2rad(1.5); /* Direct space angle tolerance + * in radians */ + + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; + + fft_in = fftw_malloc(p->nel*sizeof(double)); + fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + mod_a = modulus(ax, ay, az); + amin = mod_a * (1.0-ltol/100.0); + amax = mod_a * (1.0+ltol/100.0); + + mod_b = modulus(bx, by, bz); + bmin = mod_b * (1.0-ltol/100.0); + bmax = mod_b * (1.0+ltol/100.0); + + mod_c = modulus(cx, cy, cz); + cmin = mod_c * (1.0-ltol/100.0); + cmax = mod_c * (1.0+ltol/100.0); + + al = angle_between(bx, by, bz, cx, cy, cz); + be = angle_between(ax, ay, az, cx, cy, cz); + ga = angle_between(ax, ay, az, bx, by, bz); + + pmax = 0.0; + n = image_feature_count(image->features); + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + val = modulus(f->rx, f->ry, f->rz); + if ( val > pmax ) pmax = val; + + } + + /* Sanity check */ + if ( pmax < 1e4 ) return; + + /* Search for a */ + smin = 2.0*pmax * amin; + smax = 2.0*pmax * amax; + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom; + + new_fom = check_dir(&p->directions[i], image->features, + p->nel, pmax, fft_in, fft_out, p->plan, + smin, smax, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, + p->plan, smin, smax, th, ph, &ax, &ay, &az); + + /* Search for b */ + smin = 2.0*pmax * bmin; + smax = 2.0*pmax * bmax; + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom, ang; + + ang = angle_between(p->directions[i].x, p->directions[i].y, + p->directions[i].z, ax, ay, az); + if ( fabs(ang-ga) > angtol ) continue; + + new_fom = check_dir(&p->directions[i], image->features, + p->nel, pmax, fft_in, fft_out, p->plan, + smin, smax, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, + p->plan, smin, smax, th, ph, &bx, &by, &bz); + + /* Search for c */ + smin = 2.0*pmax * cmin; + smax = 2.0*pmax * cmax; + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom, ang; + + ang = angle_between(p->directions[i].x, p->directions[i].y, + p->directions[i].z, ax, ay, az); + if ( fabs(ang-be) > angtol ) continue; + + ang = angle_between(p->directions[i].x, p->directions[i].y, + p->directions[i].z, bx, by, bz); + if ( fabs(ang-al) > angtol ) continue; + + new_fom = check_dir(&p->directions[i], image->features, + p->nel, pmax, fft_in, fft_out, p->plan, + smin, smax, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, + p->plan, smin, smax, th, ph, &cx, &cy, &cz); + + image->candidate_cells[0] = cell_new(); + cell_set_cartesian(image->candidate_cells[0], + ax, ay, az, bx, by, bz, cx, cy, cz); + + refine_all_rigid_groups(image, image->candidate_cells[0], p->nel, pmax, + fft_in, fft_out, p->plan, smin, smax, + image->det, p); + map_all_peaks(image); + refine_cell(image, image->candidate_cells[0], image->features); + + fftw_free(fft_in); + fftw_free(fft_out); + + image->ncells = 1; +} + + +IndexingPrivate *reax_prepare() +{ + struct reax_private *p; + int samp; + double th; + + p = calloc(1, sizeof(*p)); + if ( p == NULL ) return NULL; + + p->base.indm = INDEXING_REAX; + + p->angular_inc = deg2rad(1.7); /* From Steller (1997) */ + + /* Reserve memory, over-estimating the number of directions */ + samp = 2.0*M_PI / p->angular_inc; + p->directions = malloc(samp*samp*sizeof(struct dvec)); + if ( p == NULL) { + free(p); + return NULL; + } + STATUS("Allocated space for %i directions\n", samp*samp); + + /* Generate vectors for 1D Fourier transforms */ + fesetround(1); /* Round to nearest */ + p->n_dir = 0; + for ( th=0.0; thangular_inc ) { + + double ph, phstep, n_phstep; + + n_phstep = 2.0*M_PI*sin(th)/p->angular_inc; + n_phstep = nearbyint(n_phstep); + phstep = 2.0*M_PI/n_phstep; + + for ( ph=0.0; ph<2.0*M_PI; ph+=phstep ) { + + struct dvec *dir; + + assert(p->n_dirdirections[p->n_dir++]; + + dir->x = cos(ph) * sin(th); + dir->y = sin(ph) * sin(th); + dir->z = cos(th); + dir->th = th; + dir->ph = ph; + + } + + } + STATUS("Generated %i directions (angular increment %.3f deg)\n", + p->n_dir, rad2deg(p->angular_inc)); + + p->nel = 1024; + + /* These arrays are not actually used */ + p->fft_in = fftw_malloc(p->nel*sizeof(double)); + p->fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); + + p->plan = fftw_plan_dft_r2c_1d(p->nel, p->fft_in, p->fft_out, + FFTW_MEASURE); + + p->cw = 128; p->ch = 128; + + /* Also not used */ + p->r_fft_in = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); + p->r_fft_out = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); + + p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, + 1, FFTW_MEASURE); + + return (IndexingPrivate *)p; +} + + +void reax_cleanup(IndexingPrivate *pp) +{ + struct reax_private *p; + + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; + + fftw_destroy_plan(p->plan); + fftw_free(p->fft_in); + fftw_free(p->fft_out); + + fftw_destroy_plan(p->r_plan); + fftw_free(p->r_fft_in); + fftw_free(p->r_fft_out); + + free(p); +} diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h new file mode 100644 index 00000000..543cd0d5 --- /dev/null +++ b/libcrystfel/src/reax.h @@ -0,0 +1,27 @@ +/* + * reax.h + * + * A new auto-indexer + * + * (c) 2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef REAX_H +#define REAX_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "cell.h" + +extern IndexingPrivate *reax_prepare(void); +extern void reax_cleanup(IndexingPrivate *pp); + +extern void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell); + +#endif /* REAX_H */ diff --git a/libcrystfel/src/render.c b/libcrystfel/src/render.c new file mode 100644 index 00000000..dbaca2fd --- /dev/null +++ b/libcrystfel/src/render.c @@ -0,0 +1,442 @@ +/* + * render.c + * + * Render a high dynamic range buffer in some sensible way + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#ifdef HAVE_GTK +#include +#endif + +#include +#include +#include + +#ifdef HAVE_TIFF +#include +#endif + + +#include "hdf5-file.h" +#include "render.h" +#include "peaks.h" +#include "filters.h" +#include "utils.h" + + +static void render_rgb(double val, double max, + double *rp, double *gp, double *bp) +{ + int s; + double p; + double r, g, b; + + s = val / (max/6); + p = fmod(val, max/6.0); + p /= (max/6.0); + + r = 0.0; g = 0.0; b = 0.0; + + if ( (val < 0.0) ) { + s = 0; + p = 0; + } + if ( (val > max) ) { + s = 6; + } + switch ( s ) { + case 0 : { /* Black to blue */ + r = 0; g = 0; b = p; + break; + } + case 1 : { /* Blue to pink */ + r = p; g = 0; b = 1.0; + break; + } + case 2 : { /* Pink to red */ + r = 1.0; g = 0; b = (1.0-p)*1.0; + break; + } + case 3 : { /* Red to Orange */ + r = 1.0; g = 0.5*p; b = 0; + break; + } + case 4 : { /* Orange to Yellow */ + r = 1.0; g = 0.5 + 0.5*p; b = 0; + break; + } + case 5 : { /* Yellow to White */ + r = 1.0; g = 1.0; b = 1.0*p; + break; + } + case 6 : { /* Pixel has hit the maximum value */ + r = 1.0; g = 1.0; b = 1.0; + break; + } + } + + *rp = r; + *gp = g; + *bp = b; +} + + +static void render_ratio(double val, double max, + double *rp, double *gp, double *bp) +{ + if ( val <= 1.0 ) { + render_rgb(val, 2.0, rp, gp, bp); + } else { + /* Your homework is to simplify this expression */ + val = ((val-1.0)/(max-1.0)) * (max/2.0) + max/2.0; + render_rgb(val, max, rp, gp, bp); + } +} + + +static void render_mono(double val, double max, + double *rp, double *gp, double *bp) +{ + double p; + p = val / max; + if ( val < 0.0 ) p = 0.0; + if ( val > max ) p = 1.0; + *rp = p; + *gp = p; + *bp = p; +} + + +static void render_invmono(double val, double max, + double *rp, double *gp, double *bp) +{ + double p; + p = val / max; + p = 1.0 - p; + if ( val < 0.0 ) p = 1.0; + if ( val > max ) p = 0.0; + *rp = p; + *gp = p; + *bp = p; +} + + +void render_scale(double val, double max, int scale, + double *rp, double *gp, double *bp) +{ + switch ( scale ) { + case SCALE_COLOUR : + render_rgb(val, max, rp, gp, bp); + break; + case SCALE_MONO : + render_mono(val, max, rp, gp, bp); + break; + case SCALE_INVMONO : + render_invmono(val, max, rp, gp, bp); + break; + case SCALE_RATIO : + render_ratio(val, max, rp, gp, bp); + break; + } +} + + +#ifdef HAVE_GTK + +static float *get_binned_panel(struct image *image, int binning, + int min_fs, int max_fs, int min_ss, int max_ss) +{ + float *data; + int x, y; + int w, h; + int fw; + float *in; + + fw = image->width; + in = image->data; + + /* Some pixels might get discarded */ + w = (max_fs - min_fs + 1) / binning; + h = (max_ss - min_ss + 1) / binning; + + data = malloc(w*h*sizeof(float)); + + for ( x=0; xwidth*image->height; i++ ) { + if ( image->data[i] > max ) max = image->data[i]; + } + hdr = get_binned_panel(image, binning, min_fs, max_fs, min_ss, max_ss); + if ( hdr == NULL ) return NULL; + + /* Rendered (colourful) version */ + data = malloc(3*w*h); + if ( data == NULL ) { + free(hdr); + return NULL; + } + + max /= boost; + if ( max <= 6 ) { max = 10; } + /* These x,y coordinates are measured relative to the bottom-left + * corner */ + for ( y=0; ydet->n_panels; + GdkPixbuf **pixbufs; + + pixbufs = calloc(np, sizeof(GdkPixbuf*)); + if ( pixbufs == NULL ) { + *n_pixbufs = 0; + return NULL; + } + + for ( i=0; idet->panels[i].min_fs, + image->det->panels[i].max_fs, + image->det->panels[i].min_ss, + image->det->panels[i].max_ss); + + } + + *n_pixbufs = np; + + return pixbufs; +} + + +GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale) +{ + guchar *data; + size_t x, y; + int max; + + data = malloc(3*w*h); + if ( data == NULL ) return NULL; + + max = h; + + for ( y=0; ywidth); + TIFFSetField(th, TIFFTAG_IMAGELENGTH, image->height); + TIFFSetField(th, TIFFTAG_SAMPLESPERPIXEL, 1); + TIFFSetField(th, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP); + TIFFSetField(th, TIFFTAG_BITSPERSAMPLE, 32); + TIFFSetField(th, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); + TIFFSetField(th, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); + TIFFSetField(th, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); + TIFFSetField(th, TIFFTAG_ROWSPERSTRIP, + TIFFDefaultStripSize(th, image->width*4)); + + line = _TIFFmalloc(TIFFScanlineSize(th)); + for ( y=0; yheight; y++ ) { + memcpy(line, &image->data[(image->height-1-y)*image->width], + image->width*4); + TIFFWriteScanline(th, line, y, 0); + } + _TIFFfree(line); + + TIFFClose(th); + +#else + STATUS("No TIFF support.\n"); +#endif + return 0; +} + + +int render_tiff_int16(struct image *image, const char *filename, double boost) +{ +#ifdef HAVE_TIFF + TIFF *th; + int16_t *line; + int x, y; + double max; + + th = TIFFOpen(filename, "w"); + if ( th == NULL ) return 1; + + TIFFSetField(th, TIFFTAG_IMAGEWIDTH, image->width); + TIFFSetField(th, TIFFTAG_IMAGELENGTH, image->height); + TIFFSetField(th, TIFFTAG_SAMPLESPERPIXEL, 1); + TIFFSetField(th, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_INT); /* (signed) */ + TIFFSetField(th, TIFFTAG_BITSPERSAMPLE, 16); + TIFFSetField(th, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); + TIFFSetField(th, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); + TIFFSetField(th, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); + TIFFSetField(th, TIFFTAG_ROWSPERSTRIP, + TIFFDefaultStripSize(th, image->width*4)); + + line = _TIFFmalloc(TIFFScanlineSize(th)); + max = 0.0; + for ( y=0; yheight; y++ ) { + for ( x=0;xwidth; x++ ) { + double val; + val = image->data[x+image->height*y]; + if ( val > max ) max = val; + } + } + max /= 32767.0; + + for ( y=0; yheight; y++ ) { + for ( x=0;xwidth; x++ ) { + + double val; + + val = image->data[x+(image->height-1-y)*image->width]; + val *= ((double)boost/max); + + /* Clamp to 16-bit range, + * and work round inability of most readers to deal + * with signed integers. */ + val += 1000.0; + if ( val > 32767.0 ) val = 32767.0; + if ( val < 0.0 ) val = 0.0; + + line[x] = val; + } + + TIFFWriteScanline(th, line, y, 0); + } + _TIFFfree(line); + + TIFFClose(th); +#else + STATUS("No TIFF support.\n"); +#endif + return 0; +} + +#endif /* HAVE_GTK */ diff --git a/libcrystfel/src/render.h b/libcrystfel/src/render.h new file mode 100644 index 00000000..fbfbbdb2 --- /dev/null +++ b/libcrystfel/src/render.h @@ -0,0 +1,55 @@ +/* + * render.h + * + * Render a high dynamic range buffer in some sensible way + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef RENDER_H +#define RENDER_H + + +#include + +#include "image.h" + +enum { + SCALE_COLOUR, + SCALE_MONO, + SCALE_INVMONO, + SCALE_RATIO +}; + +/* Colour scale lookup */ +extern void render_scale(double val, double max, int scale, + double *rp, double *gp, double *bp); + + +#ifdef HAVE_GTK + +#include + +extern GdkPixbuf **render_panels(struct image *image, + int binning, int scale, double boost, + int *n_pixbufs); + +extern GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale); + +extern int render_tiff_fp(struct image *image, const char *filename); +extern int render_tiff_int16(struct image *image, const char *filename, + double boost); + +#endif /* HAVE_GTK */ + + + +#endif /* RENDER_H */ diff --git a/src/dirax.c b/src/dirax.c deleted file mode 100644 index bb1cb808..00000000 --- a/src/dirax.c +++ /dev/null @@ -1,529 +0,0 @@ -/* - * dirax.c - * - * Invoke the DirAx auto-indexing program - * - * (c) 2006-2010 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#if HAVE_FORKPTY_LINUX -#include -#elif HAVE_FORKPTY_BSD -#include -#endif - - -#include "image.h" -#include "dirax.h" -#include "utils.h" -#include "peaks.h" - - -#define DIRAX_VERBOSE 0 - -#define MAX_DIRAX_CELL_CANDIDATES (5) - - -typedef enum { - DIRAX_INPUT_NONE, - DIRAX_INPUT_LINE, - DIRAX_INPUT_PROMPT, - DIRAX_INPUT_ACL -} DirAxInputType; - - -struct dirax_data { - - /* DirAx auto-indexing low-level stuff */ - int pty; - pid_t pid; - char *rbuffer; - int rbufpos; - int rbuflen; - - /* DirAx auto-indexing high-level stuff */ - int step; - int finished_ok; - int read_cell; - int best_acl; - int best_acl_nh; - int acls_tried[MAX_CELL_CANDIDATES]; - int n_acls_tried; - -}; - - -static void dirax_parseline(const char *line, struct image *image, - struct dirax_data *dirax) -{ - int rf, i, di, acl, acl_nh; - float d; - - #if DIRAX_VERBOSE - char *copy; - - copy = strdup(line); - for ( i=0; iread_cell = 1; - image->candidate_cells[image->ncells] = cell_new(); - return; - } - i++; - } - - /* Parse unit cell vectors as appropriate */ - if ( dirax->read_cell == 1 ) { - /* First row of unit cell values */ - float ax, ay, az; - int r; - r = sscanf(line, "%f %f %f %f %f %f", - &d, &d, &d, &ax, &ay, &az); - if ( r != 6 ) { - ERROR("Couldn't understand cell line:\n"); - ERROR("'%s'\n", line); - dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); - return; - } - cell_set_cartesian_a(image->candidate_cells[image->ncells], - ax*1e-10, ay*1e-10, az*1e-10); - dirax->read_cell++; - return; - } else if ( dirax->read_cell == 2 ) { - /* Second row of unit cell values */ - float bx, by, bz; - int r; - r = sscanf(line, "%f %f %f %f %f %f", - &d, &d, &d, &bx, &by, &bz); - if ( r != 6 ) { - ERROR("Couldn't understand cell line:\n"); - ERROR("'%s'\n", line); - dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); - return; - } - cell_set_cartesian_b(image->candidate_cells[image->ncells], - bx*1e-10, by*1e-10, bz*1e-10); - dirax->read_cell++; - return; - } else if ( dirax->read_cell == 3 ) { - /* Third row of unit cell values */ - float cx, cy, cz; - int r; - r = sscanf(line, "%f %f %f %f %f %f", - &d, &d, &d, &cx, &cy, &cz); - if ( r != 6 ) { - ERROR("Couldn't understand cell line:\n"); - ERROR("'%s'\n", line); - dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); - return; - } - cell_set_cartesian_c(image->candidate_cells[image->ncells++], - cx*1e-10, cy*1e-10, cz*1e-10); - dirax->read_cell = 0; - return; - } - - dirax->read_cell = 0; - - if ( sscanf(line, "%i %i %f %f %f %f %f %f %i", &acl, &acl_nh, - &d, &d, &d, &d, &d, &d, &di) == 9 ) { - if ( acl_nh > dirax->best_acl_nh ) { - - int i, found = 0; - - for ( i=0; in_acls_tried; i++ ) { - if ( dirax->acls_tried[i] == acl ) found = 1; - } - - if ( !found ) { - dirax->best_acl_nh = acl_nh; - dirax->best_acl = acl; - } - - } - } -} - - -static void dirax_sendline(const char *line, struct dirax_data *dirax) -{ - #if DIRAX_VERBOSE - char *copy; - int i; - - copy = strdup(line); - for ( i=0; ipty, line, strlen(line)); -} - - -static void dirax_send_next(struct image *image, struct dirax_data *dirax) -{ - char tmp[32]; - - switch ( dirax->step ) { - - case 1 : - dirax_sendline("\\echo off\n", dirax); - break; - - case 2 : - snprintf(tmp, 31, "read xfel-%i.drx\n", image->id); - dirax_sendline(tmp, dirax); - break; - - case 3 : - dirax_sendline("dmax 1000\n", dirax); - break; - - case 4 : - dirax_sendline("indexfit 2\n", dirax); - break; - - case 5 : - dirax_sendline("levelfit 1000\n", dirax); - break; - - case 6 : - dirax_sendline("go\n", dirax); - dirax->finished_ok = 1; - break; - - case 7 : - dirax_sendline("acl\n", dirax); - break; - - case 8 : - if ( dirax->best_acl_nh == 0 ) { - /* At this point, DirAx is presenting its ACL prompt - * and waiting for a single number. Use an extra - * newline to choose automatic ACL selection before - * exiting. */ - dirax_sendline("\nexit\n", dirax); - break; - } - snprintf(tmp, 31, "%i\n", dirax->best_acl); - dirax->acls_tried[dirax->n_acls_tried++] = dirax->best_acl; - dirax_sendline(tmp, dirax); - break; - - case 9 : - dirax_sendline("cell\n", dirax); - break; - - case 10 : - if ( dirax->n_acls_tried == MAX_DIRAX_CELL_CANDIDATES ) { - dirax_sendline("exit\n", dirax); - } else { - /* Go back round for another cell */ - dirax->best_acl_nh = 0; - dirax->step = 7; - dirax_sendline("acl\n", dirax); - } - break; - - default: - dirax_sendline("exit\n", dirax); - return; - - } - - dirax->step++; -} - - -static int dirax_readable(struct image *image, struct dirax_data *dirax) -{ - int rval; - int no_string = 0; - - rval = read(dirax->pty, dirax->rbuffer+dirax->rbufpos, - dirax->rbuflen-dirax->rbufpos); - - if ( (rval == -1) || (rval == 0) ) return 1; - - dirax->rbufpos += rval; - assert(dirax->rbufpos <= dirax->rbuflen); - - while ( (!no_string) && (dirax->rbufpos > 0) ) { - - int i; - int block_ready = 0; - DirAxInputType type = DIRAX_INPUT_NONE; - - /* See if there's a full line in the buffer yet */ - for ( i=0; irbufpos-1; i++ ) { - /* Means the last value looked at is rbufpos-2 */ - - /* Is there a prompt in the buffer? */ - if ( (i+7 <= dirax->rbufpos) - && (!strncmp(dirax->rbuffer+i, "Dirax> ", 7)) ) { - block_ready = 1; - type = DIRAX_INPUT_PROMPT; - break; - } - - /* How about an ACL prompt? */ - if ( (i+10 <= dirax->rbufpos) - && (!strncmp(dirax->rbuffer+i, "acl/auto [", 10)) ) { - block_ready = 1; - type = DIRAX_INPUT_ACL; - break; - } - - if ( (dirax->rbuffer[i] == '\r') - && (dirax->rbuffer[i+1] == '\n') ) { - block_ready = 1; - type = DIRAX_INPUT_LINE; - break; - } - - } - - if ( block_ready ) { - - unsigned int new_rbuflen; - unsigned int endbit_length; - char *block_buffer = NULL; - - switch ( type ) { - - case DIRAX_INPUT_LINE : - - block_buffer = malloc(i+1); - memcpy(block_buffer, dirax->rbuffer, i); - block_buffer[i] = '\0'; - - if ( block_buffer[0] == '\r' ) { - memmove(block_buffer, block_buffer+1, i); - } - - dirax_parseline(block_buffer, image, dirax); - free(block_buffer); - endbit_length = i+2; - break; - - case DIRAX_INPUT_PROMPT : - - dirax_send_next(image, dirax); - endbit_length = i+7; - break; - - case DIRAX_INPUT_ACL : - - dirax_send_next(image,dirax ); - endbit_length = i+10; - break; - - default : - - /* Obviously, this never happens :) */ - ERROR("Unrecognised DirAx input mode! " - "I don't know how to understand DirAx\n"); - return 1; - - } - - /* Now the block's been parsed, it should be - * forgotten about */ - memmove(dirax->rbuffer, - dirax->rbuffer + endbit_length, - dirax->rbuflen - endbit_length); - - /* Subtract the number of bytes removed */ - dirax->rbufpos = dirax->rbufpos - - endbit_length; - new_rbuflen = dirax->rbuflen - endbit_length; - if ( new_rbuflen == 0 ) new_rbuflen = 256; - dirax->rbuffer = realloc(dirax->rbuffer, - new_rbuflen); - dirax->rbuflen = new_rbuflen; - - } else { - - if ( dirax->rbufpos==dirax->rbuflen ) { - - /* More buffer space is needed */ - dirax->rbuffer = realloc( - dirax->rbuffer, - dirax->rbuflen + 256); - dirax->rbuflen = dirax->rbuflen - + 256; - /* The new space gets used at the next - * read, shortly... */ - - } - no_string = 1; - - } - - } - - return 0; -} - - -static void write_drx(struct image *image) -{ - FILE *fh; - int i; - char filename[1024]; - - snprintf(filename, 1023, "xfel-%i.drx", image->id); - - fh = fopen(filename, "w"); - if ( !fh ) { - ERROR("Couldn't open temporary file '%s'\n", filename); - return; - } - fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. */ - - for ( i=0; ifeatures); i++ ) { - - struct imagefeature *f; - - f = image_get_feature(image->features, i); - if ( f == NULL ) continue; - - fprintf(fh, "%10f %10f %10f %8f\n", - f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0); - - } - fclose(fh); -} - - -void run_dirax(struct image *image) -{ - unsigned int opts; - int status; - int rval; - struct dirax_data *dirax; - - write_drx(image); - - dirax = malloc(sizeof(struct dirax_data)); - if ( dirax == NULL ) { - ERROR("Couldn't allocate memory for DirAx data.\n"); - return; - } - - dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL); - if ( dirax->pid == -1 ) { - ERROR("Failed to fork for DirAx\n"); - return; - } - if ( dirax->pid == 0 ) { - - /* Child process: invoke DirAx */ - struct termios t; - - /* Turn echo off */ - tcgetattr(STDIN_FILENO, &t); - t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); - tcsetattr(STDIN_FILENO, TCSANOW, &t); - - execlp("dirax", "", (char *)NULL); - ERROR("Failed to invoke DirAx.\n"); - _exit(0); - - } - - dirax->rbuffer = malloc(256); - dirax->rbuflen = 256; - dirax->rbufpos = 0; - - /* Set non-blocking */ - opts = fcntl(dirax->pty, F_GETFL); - fcntl(dirax->pty, F_SETFL, opts | O_NONBLOCK); - - dirax->step = 1; /* This starts the "initialisation" procedure */ - dirax->finished_ok = 0; - dirax->read_cell = 0; - dirax->n_acls_tried = 0; - dirax->best_acl_nh = 0; - - do { - - fd_set fds; - struct timeval tv; - int sval; - - FD_ZERO(&fds); - FD_SET(dirax->pty, &fds); - - tv.tv_sec = 10; - tv.tv_usec = 0; - - sval = select(dirax->pty+1, &fds, NULL, NULL, &tv); - - if ( sval == -1 ) { - int err = errno; - ERROR("select() failed: %s\n", strerror(err)); - rval = 1; - } else if ( sval != 0 ) { - rval = dirax_readable(image, dirax); - } else { - ERROR("No response from DirAx..\n"); - rval = 1; - } - - } while ( !rval ); - - close(dirax->pty); - free(dirax->rbuffer); - waitpid(dirax->pid, &status, 0); - - if ( dirax->finished_ok == 0 ) { - ERROR("DirAx doesn't seem to be working properly.\n"); - } - - free(dirax); -} diff --git a/src/dirax.h b/src/dirax.h deleted file mode 100644 index 8c429710..00000000 --- a/src/dirax.h +++ /dev/null @@ -1,26 +0,0 @@ -/* - * dirax.h - * - * Invoke the DirAx auto-indexing program - * - * (c) 2006-2010 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifndef DIRAX_H -#define DIRAX_H - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include "utils.h" - - -extern void run_dirax(struct image *image); - - -#endif /* DIRAX_H */ diff --git a/src/index-priv.h b/src/index-priv.h deleted file mode 100644 index 3d8c8a22..00000000 --- a/src/index-priv.h +++ /dev/null @@ -1,29 +0,0 @@ -/* - * index-priv.h - * - * Indexing private data - * - * (c) 2006-2010 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifndef INDEXPRIV_H -#define INDEXPRIV_H - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include "index.h" - -struct _indexingprivate -{ - IndexingMethod indm; -}; - - -#endif /* INDEXPRIV_H */ diff --git a/src/index.c b/src/index.c deleted file mode 100644 index d5e76c50..00000000 --- a/src/index.c +++ /dev/null @@ -1,274 +0,0 @@ -/* - * index.c - * - * Perform indexing (somehow) - * - * (c) 2006-2011 Thomas White - * (c) 2010 Richard Kirian - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - -#include "image.h" -#include "utils.h" -#include "peaks.h" -#include "dirax.h" -#include "mosflm.h" -#include "detector.h" -#include "index.h" -#include "index-priv.h" -#include "reax.h" -#include "geometry.h" - - -/* Base class constructor for unspecialised indexing private data */ -IndexingPrivate *indexing_private(IndexingMethod indm) -{ - struct _indexingprivate *priv; - priv = calloc(1, sizeof(struct _indexingprivate)); - priv->indm = indm; - return priv; -} - - -static const char *maybes(int n) -{ - if ( n == 1 ) return ""; - return "s"; -} - - -IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, - const char *filename, struct detector *det, - double nominal_photon_energy) -{ - int n; - int nm = 0; - IndexingPrivate **iprivs; - - while ( indm[nm] != INDEXING_NONE ) nm++; - STATUS("Preparing %i indexing method%s.\n", nm, maybes(nm)); - iprivs = malloc((nm+1) * sizeof(IndexingPrivate *)); - - for ( n=0; nindm ) { - case INDEXING_NONE : - free(priv[n]); - break; - case INDEXING_DIRAX : - free(priv[n]); - break; - case INDEXING_MOSFLM : - free(priv[n]); - break; - case INDEXING_REAX : - reax_cleanup(priv[n]); - break; - } - - n++; - - } -} - - -void map_all_peaks(struct image *image) -{ - int i, n; - - n = image_feature_count(image->features); - - /* Map positions to 3D */ - for ( i=0; ifeatures, i); - if ( f == NULL ) continue; - - r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); - f->rx = r.u; f->ry = r.v; f->rz = r.w; - - } -} - - -void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, - int cellr, int verbose, IndexingPrivate **ipriv, - int config_insane) -{ - int i; - int n = 0; - - if ( indm == NULL ) return; - - map_all_peaks(image); - image->indexed_cell = NULL; - - while ( indm[n] != INDEXING_NONE ) { - - image->ncells = 0; - - /* Index as appropriate */ - switch ( indm[n] ) { - case INDEXING_NONE : - return; - case INDEXING_DIRAX : - run_dirax(image); - break; - case INDEXING_MOSFLM : - run_mosflm(image, cell); - break; - case INDEXING_REAX : - reax_index(ipriv[n], image, cell); - break; - } - if ( image->ncells == 0 ) { - n++; - continue; - } - - for ( i=0; incells; i++ ) { - - UnitCell *new_cell = NULL; - UnitCell *cand = image->candidate_cells[i]; - - if ( verbose ) { - STATUS("--------------------\n"); - STATUS("Candidate cell %i (before matching):\n", - i); - cell_print(image->candidate_cells[i]); - STATUS("--------------------\n"); - } - - /* Match or reduce the cell as appropriate */ - switch ( cellr ) { - case CELLR_NONE : - new_cell = cell_new_from_cell(cand); - break; - case CELLR_REDUCE : - new_cell = match_cell(cand, cell, verbose, 1); - break; - case CELLR_COMPARE : - new_cell = match_cell(cand, cell, verbose, 0); - break; - case CELLR_COMPARE_AB : - new_cell = match_cell_ab(cand, cell); - break; - } - - /* No cell? Move on to the next candidate */ - if ( new_cell == NULL ) continue; - - /* Sanity check */ - image->reflections = find_intersections(image, new_cell); - image->indexed_cell = new_cell; - - if ( !config_insane && !peak_sanity_check(image) ) { - cell_free(new_cell); - image->indexed_cell = NULL; - continue; - } - - goto done; /* Success */ - - } - - for ( i=0; incells; i++ ) { - cell_free(image->candidate_cells[i]); - image->candidate_cells[i] = NULL; - } - - /* Move on to the next indexing method */ - n++; - - } - -done: - for ( i=0; incells; i++ ) { - /* May free(NULL) if all algorithms were tried and no success */ - cell_free(image->candidate_cells[i]); - } -} - - -IndexingMethod *build_indexer_list(const char *str, int *need_cell) -{ - int n, i; - char **methods; - IndexingMethod *list; - int tmp; - - if ( need_cell == NULL ) need_cell = &tmp; - *need_cell = 0; - - n = assplode(str, ",", &methods, ASSPLODE_NONE); - list = malloc((n+1)*sizeof(IndexingMethod)); - - for ( i=0; i - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifndef INDEX_H -#define INDEX_H - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include "cell.h" -#include "image.h" -#include "detector.h" - - -/* Indexing methods */ -typedef enum { - INDEXING_NONE, - INDEXING_DIRAX, - INDEXING_MOSFLM, - INDEXING_REAX, -} IndexingMethod; - - -/* Cell reduction methods */ -enum { - CELLR_NONE, - CELLR_REDUCE, - CELLR_COMPARE, - CELLR_COMPARE_AB, -}; - - -typedef struct _indexingprivate IndexingPrivate; - -extern IndexingPrivate *indexing_private(IndexingMethod indm); - -extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, - const char *filename, - struct detector *det, - double nominal_photon_energy); - -extern void map_all_peaks(struct image *image); - -extern void index_pattern(struct image *image, UnitCell *cell, - IndexingMethod *indm, int cellr, int verbose, - IndexingPrivate **priv, int config_insane); - -extern void cleanup_indexing(IndexingPrivate **priv); - -extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); - -#endif /* INDEX_H */ diff --git a/src/mosflm.c b/src/mosflm.c deleted file mode 100644 index 5df5af21..00000000 --- a/src/mosflm.c +++ /dev/null @@ -1,609 +0,0 @@ -/* - * mosflm.c - * - * Invoke the DPS auto-indexing algorithm through MOSFLM - * - * (c) 2010 Richard Kirian - * (c) 2006-2010 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -/* TODO: - * - * Properly read the newmat file (don't use fscanf-- spaces between numers - * are not guaranteed) - * - * "Success" is indicated by existence of NEWMAT file written by mosflm. - * Better to interact with mosflm directly in order to somehow verify success. - * - * Investigate how these keywords affect mosflms behavior: - * - * MOSAICITY - * DISPERSION - * DIVERGENCE - * POLARISATION - * POSTREF BEAM - * POSTREF USEBEAM OFF - * PREREFINE ON - * EXTRA ON - * POSTREF ON - * - * These did not seem to affect the results by my (Rick's) experience, probably - * because they are only used conjunction with image intensity data, but it's - * worth another look at the documentation. - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#if HAVE_FORKPTY_LINUX -#include -#elif HAVE_FORKPTY_BSD -#include -#endif - - -#include "image.h" -#include "mosflm.h" -#include "utils.h" -#include "peaks.h" - - -#define MOSFLM_VERBOSE 0 - - -typedef enum { - MOSFLM_INPUT_NONE, - MOSFLM_INPUT_LINE, - MOSFLM_INPUT_PROMPT -} MOSFLMInputType; - - -struct mosflm_data { - - /* MOSFLM auto-indexing low-level stuff */ - int pty; - pid_t pid; - char *rbuffer; - int rbufpos; - int rbuflen; - - /* MOSFLM high-level stuff */ - char newmatfile[128]; - char imagefile[128]; - char sptfile[128]; - int step; - int finished_ok; - UnitCell *target_cell; - -}; - - -static void mosflm_parseline(const char *line, struct image *image, - struct mosflm_data *dirax) -{ - #if MOSFLM_VERBOSE - char *copy; - int i; - - copy = strdup(line); - for ( i=0; ilambda; - - image->candidate_cells[0] = cell_new(); - - cell_set_reciprocal(image->candidate_cells[0], - asz*c, asy*c, asx*c, - bsz*c, bsy*c, bsx*c, - csz*c, csy*c, csx*c); - - image->ncells = 1; - - return 0; -} - - -/* Need to sort mosflm peaks by intensity... */ -struct sptline { - double x; /* x coordinate of peak */ - double y; /* y coordinate of peak */ - double h; /* height of peak */ - double s; /* sigma of peak */ -}; - - -static int compare_vals(const void *ap, const void *bp) -{ - const struct sptline a = *(struct sptline *)ap; - const struct sptline b = *(struct sptline *)bp; - - if ( a.h < b.h ) return 1; - if ( a.h > b.h ) return -1; - return 0; -} - - -/* Write .spt file for mosflm */ -static void write_spt(struct image *image, const char *filename) -{ - FILE *fh; - int i; - double fclen = 67.8; /* fake camera length in mm */ - double fpix = 0.075; /* fake pixel size in mm */ - double pix; - double height = 100.0; - double sigma = 1.0; - int nPeaks = image_feature_count(image->features); - struct sptline *sptlines; - - fh = fopen(filename, "w"); - if ( !fh ) { - ERROR("Couldn't open temporary file '%s'\n", filename); - return; - } - - fprintf(fh, "%10d %10d %10.8f %10.6f %10.6f\n", 1, 1, fpix, 1.0, 0.0); - fprintf(fh, "%10d %10d\n", 1, 1); - fprintf(fh, "%10.5f %10.5f\n", 0.0, 0.0); - - sptlines = malloc(sizeof(struct sptline)*nPeaks); - - for ( i=0; ifeatures, i); - if ( f == NULL ) continue; - - p = find_panel(image->det, f->fs, f->ss); - if ( p == NULL ) continue; - - pix = 1000.0/p->res; /* pixel size in mm */ - height = f->intensity; - - xs = (f->fs-p->min_fs)*p->fsx + (f->ss-p->min_ss)*p->ssx; - ys = (f->ss-p->min_fs)*p->fsy + (f->ss-p->min_ss)*p->ssy; - rx = xs + p->cnx; - ry = ys + p->cny; - - sptlines[i].x = ry*pix*fclen/p->clen/1000.0; - sptlines[i].y = -rx*pix*fclen/p->clen/1000.0; - sptlines[i].h = height; - sptlines[i].s = sigma/1000.0; - - } - - qsort(sptlines, nPeaks, sizeof(struct sptline), compare_vals); - - for ( i=0; ipty, line, strlen(line)); -} - - -static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) -{ - char tmp[256]; - char symm[32]; - const char *sg; - double wavelength; - double a, b, c, alpha, beta, gamma; - int i, j; - - switch ( mosflm->step ) { - - case 1 : - mosflm_sendline("DETECTOR ROTATION HORIZONTAL" - " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" - " RECTANGULAR\n", mosflm); - break; - - case 2 : - if ( mosflm->target_cell != NULL ) { - cell_get_parameters(mosflm->target_cell, &a, &b, &c, - &alpha, &beta, &gamma); - snprintf(tmp, 255, - "CELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n", - a*1e10, b*1e10, c*1e10, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); - mosflm_sendline(tmp, mosflm); - } else { - mosflm_sendline("# Do nothing\n", mosflm); - } - break; - - case 3 : - if ( mosflm->target_cell != NULL ) { - sg = cell_get_spacegroup(mosflm->target_cell); - /* Remove white space from space group */ - j = 0; - for ( i=0; ilambda*1e10; - snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength); - mosflm_sendline(tmp, mosflm); - break; - - case 7 : - snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile); - mosflm_sendline(tmp, mosflm); - break; - - case 8 : - snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile); - mosflm_sendline(tmp, mosflm); - break; - - case 9 : - snprintf(tmp, 255, "AUTOINDEX DPS FILE %s" - " IMAGE 1 MAXCELL 1000 REFINE\n", - mosflm->sptfile); - - /* "This option is only available if you e-mail Andrew Leslie - * and ask for it." - MOSFLM - * snprintf(tmp, 255, "AUTOINDEX NODISPLAY IMAGE 1 FILE %s\n", - * mosflm->sptfile); */ - mosflm_sendline(tmp, mosflm); - break; - - case 10 : - mosflm_sendline("GO\n", mosflm); - mosflm->finished_ok = 1; - break; - - default: - mosflm_sendline("exit\n", mosflm); - return; - - } - - mosflm->step++; -} - - -static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) -{ - int rval; - int no_string = 0; - - rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos, - mosflm->rbuflen-mosflm->rbufpos); - - if ( (rval == -1) || (rval == 0) ) return 1; - - mosflm->rbufpos += rval; - assert(mosflm->rbufpos <= mosflm->rbuflen); - - while ( (!no_string) && (mosflm->rbufpos > 0) ) { - - int i; - int block_ready = 0; - MOSFLMInputType type = MOSFLM_INPUT_NONE; - - /* See if there's a full line in the buffer yet */ - for ( i=0; irbufpos-1; i++ ) { - /* Means the last value looked at is rbufpos-2 */ - - /* Is there a prompt in the buffer? */ - if ( (i+10 <= mosflm->rbufpos) - && (!strncmp(mosflm->rbuffer+i, "MOSFLM => ", 10)) ) { - block_ready = 1; - type = MOSFLM_INPUT_PROMPT; - break; - } - - if ( (mosflm->rbuffer[i] == '\r') - && (mosflm->rbuffer[i+1] == '\n') ) { - block_ready = 1; - type = MOSFLM_INPUT_LINE; - break; - } - - } - - if ( block_ready ) { - - unsigned int new_rbuflen; - unsigned int endbit_length; - char *block_buffer = NULL; - - switch ( type ) { - - case MOSFLM_INPUT_LINE : - - block_buffer = malloc(i+1); - memcpy(block_buffer, mosflm->rbuffer, i); - block_buffer[i] = '\0'; - - if ( block_buffer[0] == '\r' ) { - memmove(block_buffer, block_buffer+1, i); - } - - mosflm_parseline(block_buffer, image, mosflm); - free(block_buffer); - endbit_length = i+2; - break; - - case MOSFLM_INPUT_PROMPT : - - mosflm_send_next(image, mosflm); - endbit_length = i+7; - break; - - default : - - /* Obviously, this never happens :) */ - ERROR("Unrecognised MOSFLM input mode! " - "I don't know how to understand MOSFLM\n"); - return 1; - - } - - /* Now the block's been parsed, it should be - * forgotten about */ - memmove(mosflm->rbuffer, - mosflm->rbuffer + endbit_length, - mosflm->rbuflen - endbit_length); - - /* Subtract the number of bytes removed */ - mosflm->rbufpos = mosflm->rbufpos - - endbit_length; - new_rbuflen = mosflm->rbuflen - endbit_length; - if ( new_rbuflen == 0 ) new_rbuflen = 256; - mosflm->rbuffer = realloc(mosflm->rbuffer, - new_rbuflen); - mosflm->rbuflen = new_rbuflen; - - } else { - - if ( mosflm->rbufpos==mosflm->rbuflen ) { - - /* More buffer space is needed */ - mosflm->rbuffer = realloc( - mosflm->rbuffer, - mosflm->rbuflen + 256); - mosflm->rbuflen = mosflm->rbuflen + 256; - /* The new space gets used at the next - * read, shortly... */ - - } - no_string = 1; - - } - - } - - return 0; -} - - -void run_mosflm(struct image *image, UnitCell *cell) -{ - struct mosflm_data *mosflm; - unsigned int opts; - int status; - int rval; - - mosflm = malloc(sizeof(struct mosflm_data)); - if ( mosflm == NULL ) { - ERROR("Couldn't allocate memory for MOSFLM data.\n"); - return; - } - - mosflm->target_cell = cell; - - snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id); - write_img(image, mosflm->imagefile); /* Dummy image */ - - snprintf(mosflm->sptfile, 127, "xfel-%i_001.spt", image->id); - write_spt(image, mosflm->sptfile); - - snprintf(mosflm->newmatfile, 127, "xfel-%i.newmat", image->id); - remove(mosflm->newmatfile); - - mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL); - if ( mosflm->pid == -1 ) { - ERROR("Failed to fork for MOSFLM\n"); - free(mosflm); - return; - } - if ( mosflm->pid == 0 ) { - - /* Child process: invoke MOSFLM */ - struct termios t; - - /* Turn echo off */ - tcgetattr(STDIN_FILENO, &t); - t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); - tcsetattr(STDIN_FILENO, TCSANOW, &t); - - execlp("ipmosflm", "", (char *)NULL); - ERROR("Failed to invoke MOSFLM.\n"); - _exit(0); - - } - - mosflm->rbuffer = malloc(256); - mosflm->rbuflen = 256; - mosflm->rbufpos = 0; - - /* Set non-blocking */ - opts = fcntl(mosflm->pty, F_GETFL); - fcntl(mosflm->pty, F_SETFL, opts | O_NONBLOCK); - - mosflm->step = 1; /* This starts the "initialisation" procedure */ - mosflm->finished_ok = 0; - - do { - - fd_set fds; - struct timeval tv; - int sval; - - FD_ZERO(&fds); - FD_SET(mosflm->pty, &fds); - - tv.tv_sec = 10; - tv.tv_usec = 0; - - sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv); - - if ( sval == -1 ) { - int err = errno; - ERROR("select() failed: %s\n", strerror(err)); - rval = 1; - } else if ( sval != 0 ) { - rval = mosflm_readable(image, mosflm); - } else { - ERROR("No response from MOSFLM..\n"); - rval = 1; - } - - } while ( !rval ); - - close(mosflm->pty); - free(mosflm->rbuffer); - waitpid(mosflm->pid, &status, 0); - - if ( mosflm->finished_ok == 0 ) { - ERROR("MOSFLM doesn't seem to be working properly.\n"); - } else { - /* Read the mosflm NEWMAT file and get cell if found */ - read_newmat(mosflm->newmatfile, image); - } - - free(mosflm); -} diff --git a/src/mosflm.h b/src/mosflm.h deleted file mode 100644 index 82b80a44..00000000 --- a/src/mosflm.h +++ /dev/null @@ -1,27 +0,0 @@ -/* - * mosflm.h - * - * Invoke the DPS auto-indexing algorithm through MOSFLM - * - * (c) 2010 Richard Kirian - * (c) 2006-2010 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifndef MOSFLM_H -#define MOSFLM_H - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include "utils.h" - - -extern void run_mosflm(struct image *image, UnitCell *cell); - - -#endif /* MOSFLM_H */ diff --git a/src/reax.c b/src/reax.c deleted file mode 100644 index ff1ea582..00000000 --- a/src/reax.c +++ /dev/null @@ -1,728 +0,0 @@ -/* - * reax.c - * - * A new auto-indexer - * - * (c) 2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "image.h" -#include "utils.h" -#include "peaks.h" -#include "cell.h" -#include "index.h" -#include "index-priv.h" - - -struct dvec -{ - double x; - double y; - double z; - double th; - double ph; -}; - - -struct reax_private -{ - IndexingPrivate base; - struct dvec *directions; - int n_dir; - double angular_inc; - - double *fft_in; - fftw_complex *fft_out; - fftw_plan plan; - int nel; - - fftw_complex *r_fft_in; - fftw_complex *r_fft_out; - fftw_plan r_plan; - int ch; - int cw; -}; - - -static double check_dir(struct dvec *dir, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, - const char *rg, struct detector *det) -{ - int n, i; - double tot; - - for ( i=0; ifs, f->ss); - assert(p != NULL); - - if ( p->rigid_group != rg ) continue; - - } - - val = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z; - - idx = nel/2 + nel*val/(2.0*pmax); - fft_in[idx]++; - - } - - fftw_execute_dft_r2c(plan, fft_in, fft_out); - - tot = 0.0; - for ( i=smin; i<=smax; i++ ) { - double re, im; - re = fft_out[i][0]; - im = fft_out[i][1]; - tot += sqrt(re*re + im*im); - } - - return tot; -} - - -/* Refine a direct space vector. From Clegg (1984) */ -static double iterate_refine_vector(double *x, double *y, double *z, - ImageFeatureList *flist) -{ - int fi, n, err; - gsl_matrix *C; - gsl_vector *A; - gsl_vector *t; - gsl_matrix *s_vec; - gsl_vector *s_val; - double dtmax; - - A = gsl_vector_calloc(3); - C = gsl_matrix_calloc(3, 3); - - n = image_feature_count(flist); - fesetround(1); - for ( fi=0; firx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */ - kn = nearbyint(kno); - if ( kn - kno > 0.3 ) continue; - - xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz; - - for ( i=0; i<3; i++ ) { - - val = gsl_vector_get(A, i); - gsl_vector_set(A, i, val+xv[i]*kn); - - for ( j=0; j<3; j++ ) { - val = gsl_matrix_get(C, i, j); - gsl_matrix_set(C, i, j, val+xv[i]*xv[j]); - } - - } - - } - - s_val = gsl_vector_calloc(3); - s_vec = gsl_matrix_calloc(3, 3); - err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val); - if ( err ) { - ERROR("SVD failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(C); - gsl_vector_free(A); - return 0.0; - } - - t = gsl_vector_calloc(3); - err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t); - if ( err ) { - ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(C); - gsl_vector_free(A); - gsl_vector_free(t); - return 0.0; - } - - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - - dtmax = fabs(*x - gsl_vector_get(t, 0)); - dtmax += fabs(*y - gsl_vector_get(t, 1)); - dtmax += fabs(*z - gsl_vector_get(t, 2)); - - *x = gsl_vector_get(t, 0); - *y = gsl_vector_get(t, 1); - *z = gsl_vector_get(t, 2); - - gsl_matrix_free(C); - gsl_vector_free(A); - - return dtmax; -} - - -static void refine_cell(struct image *image, UnitCell *cell, - ImageFeatureList *flist) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int i; - double sm; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - i = 0; - do { - - sm = iterate_refine_vector(&ax, &ay, &az, flist); - sm += iterate_refine_vector(&bx, &by, &bz, flist); - sm += iterate_refine_vector(&cx, &cy, &cz, flist); - i++; - - } while ( (sm > 0.001e-9) && (i<10) ); - - cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz); - - if ( i == 10 ) { - cell_free(image->indexed_cell); - image->indexed_cell = NULL; - } -} - - -static void fine_search(struct reax_private *p, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, double th_cen, double ph_cen, - double *x, double *y, double *z) -{ - double fom = 0.0; - double th, ph; - double inc; - struct dvec dir; - int i, s; - double max, modv; - - inc = p->angular_inc / 100.0; - - for ( th=th_cen-p->angular_inc; th<=th_cen+p->angular_inc; th+=inc ) { - for ( ph=ph_cen-p->angular_inc; ph<=ph_cen+p->angular_inc; ph+=inc ) { - - double new_fom; - - dir.x = cos(ph) * sin(th); - dir.y = sin(ph) * sin(th); - dir.z = cos(th); - - new_fom = check_dir(&dir, flist, nel, pmax, fft_in, fft_out, - plan, smin, smax, NULL, NULL); - - if ( new_fom > fom ) { - fom = new_fom; - *x = dir.x; *y = dir.y; *z = dir.z; - } - - } - } - - s = -1; - max = 0.0; - for ( i=smin; i<=smax; i++ ) { - - double re, im, m; - - re = fft_out[i][0]; - im = fft_out[i][1]; - m = sqrt(re*re + im*im); - if ( m > max ) { - max = m; - s = i; - } - - } - assert(s>0); - - modv = (double)s / (2.0*pmax); - *x *= modv; *y *= modv; *z *= modv; -} - - -static double get_model_phase(double x, double y, double z, ImageFeatureList *f, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, const char *rg, - struct detector *det) -{ - struct dvec dir; - int s, i; - double max; - double re, im; - - dir.x = x; dir.y = y; dir.z = z; - - check_dir(&dir, f, nel, pmax, fft_in,fft_out, plan, smin, smax, - rg, det); - - s = -1; - max = 0.0; - for ( i=smin; i<=smax; i++ ) { - - double re, im, m; - - re = fft_out[i][0]; - im = fft_out[i][1]; - m = sqrt(re*re + im*im); - if ( m > max ) { - max = m; - s = i; - } - - } - - re = fft_out[s][0]; - im = fft_out[s][1]; - - return atan2(im, re); -} - - -static void refine_rigid_group(struct image *image, UnitCell *cell, - const char *rg, int nel, double pmax, - double *fft_in, fftw_complex *fft_out, - fftw_plan plan, int smin, int smax, - struct detector *det, struct reax_private *pr) -{ - double ax, ay, az, ma; - double bx, by, bz, mb; - double cx, cy, cz, mc; - double pha, phb, phc; - struct panel *p; - int i, j; - fftw_complex *r_fft_in; - fftw_complex *r_fft_out; - double m2m; - signed int aix, aiy; - signed int bix, biy; - signed int cix, ciy; - double max; - int max_i, max_j; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - ma = modulus(ax, ay, az); - mb = modulus(bx, by, bz); - mc = modulus(cx, cy, cz); - - pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, nel, pmax, - fft_in, fft_out, plan, smin, smax, rg, det); - phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, nel, pmax, - fft_in, fft_out, plan, smin, smax, rg, det); - phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, nel, pmax, - fft_in, fft_out, plan, smin, smax, rg, det); - - for ( i=0; in_panels; i++ ) { - if ( det->panels[i].rigid_group == rg ) { - p = &det->panels[i]; - break; - } - } - - r_fft_in = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); - r_fft_out = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); - for ( i=0; icw; i++ ) { - for ( j=0; jch; j++ ) { - r_fft_in[i+pr->cw*j][0] = 0.0; - r_fft_in[i+pr->cw*j][1] = 0.0; - } - } - - ma = modulus(ax, ay, 0.0); - mb = modulus(bx, by, 0.0); - mc = modulus(cx, cy, 0.0); - m2m = ma; - if ( mb > m2m ) m2m = mb; - if ( mc > m2m ) m2m = mc; - - aix = (pr->cw/2)*ax/m2m; aiy = (pr->ch/2)*ay/m2m; - bix = (pr->cw/2)*bx/m2m; biy = (pr->ch/2)*by/m2m; - cix = (pr->cw/2)*cx/m2m; ciy = (pr->ch/2)*cy/m2m; - - if ( aix < 0 ) aix += pr->cw/2; - if ( bix < 0 ) bix += pr->cw/2; - if ( cix < 0 ) cix += pr->cw/2; - - if ( aiy < 0 ) aiy += pr->ch/2; - if ( biy < 0 ) biy += pr->ch/2; - if ( ciy < 0 ) ciy += pr->ch/2; - - r_fft_in[aix + pr->cw*aiy][0] = cos(pha); - r_fft_in[aix + pr->cw*aiy][1] = sin(pha); - r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][0] = cos(pha); - r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][1] = -sin(pha); - - r_fft_in[bix + pr->cw*biy][0] = cos(phb); - r_fft_in[bix + pr->cw*biy][1] = sin(phb); - r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][0] = cos(phb); - r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][1] = -sin(phb); - - r_fft_in[cix + pr->cw*ciy][0] = cos(phc); - r_fft_in[cix + pr->cw*ciy][1] = sin(phc); - r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][0] = cos(phc); - r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][1] = -sin(phc); - - const int tidx = 1; - r_fft_in[tidx][0] = 1.0; - r_fft_in[tidx][1] = 0.0; - -// STATUS("%i %i\n", aix, aiy); -// STATUS("%i %i\n", bix, biy); -// STATUS("%i %i\n", cix, ciy); - - fftw_execute_dft(pr->r_plan, r_fft_in, r_fft_out); - -// max = 0.0; -// FILE *fh = fopen("centering.dat", "w"); -// for ( i=0; icw; i++ ) { -// for ( j=0; jch; j++ ) { -// -// double re, im, am, ph; -// -// re = r_fft_out[i + pr->cw*j][0]; -// im = r_fft_out[i + pr->cw*j][1]; -// am = sqrt(re*re + im*im); -// ph = atan2(im, re); -// -// if ( am > max ) { -// max = am; -// max_i = i; -// max_j = j; -// } -// -// fprintf(fh, "%f ", am); -// -// } -// fprintf(fh, "\n"); -// } -// STATUS("Max at %i, %i\n", max_i, max_j); -// fclose(fh); -// exit(1); - -// STATUS("Offsets for '%s': %.2f, %.2f pixels\n", rg, dx, dy); -} - - -static void refine_all_rigid_groups(struct image *image, UnitCell *cell, - int nel, double pmax, - double *fft_in, fftw_complex *fft_out, - fftw_plan plan, int smin, int smax, - struct detector *det, - struct reax_private *p) -{ - int i; - - for ( i=0; idet->num_rigid_groups; i++ ) { - refine_rigid_group(image, cell, image->det->rigid_groups[i], - nel, pmax, fft_in, fft_out, plan, smin, smax, - det, p); - } -} - - -void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) -{ - int i; - struct reax_private *p; - double fom; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double mod_a, mod_b, mod_c; - double al, be, ga; - double th, ph; - double *fft_in; - fftw_complex *fft_out; - int smin, smax; - double amin, amax; - double bmin, bmax; - double cmin, cmax; - double pmax; - int n; - const double ltol = 5.0; /* Direct space axis length - * tolerance in percent */ - const double angtol = deg2rad(1.5); /* Direct space angle tolerance - * in radians */ - - assert(pp->indm == INDEXING_REAX); - p = (struct reax_private *)pp; - - fft_in = fftw_malloc(p->nel*sizeof(double)); - fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - mod_a = modulus(ax, ay, az); - amin = mod_a * (1.0-ltol/100.0); - amax = mod_a * (1.0+ltol/100.0); - - mod_b = modulus(bx, by, bz); - bmin = mod_b * (1.0-ltol/100.0); - bmax = mod_b * (1.0+ltol/100.0); - - mod_c = modulus(cx, cy, cz); - cmin = mod_c * (1.0-ltol/100.0); - cmax = mod_c * (1.0+ltol/100.0); - - al = angle_between(bx, by, bz, cx, cy, cz); - be = angle_between(ax, ay, az, cx, cy, cz); - ga = angle_between(ax, ay, az, bx, by, bz); - - pmax = 0.0; - n = image_feature_count(image->features); - for ( i=0; ifeatures, i); - if ( f == NULL ) continue; - - val = modulus(f->rx, f->ry, f->rz); - if ( val > pmax ) pmax = val; - - } - - /* Sanity check */ - if ( pmax < 1e4 ) return; - - /* Search for a */ - smin = 2.0*pmax * amin; - smax = 2.0*pmax * amax; - fom = 0.0; th = 0.0; ph = 0.0; - for ( i=0; in_dir; i++ ) { - - double new_fom; - - new_fom = check_dir(&p->directions[i], image->features, - p->nel, pmax, fft_in, fft_out, p->plan, - smin, smax, NULL, NULL); - if ( new_fom > fom ) { - fom = new_fom; - th = p->directions[i].th; - ph = p->directions[i].ph; - } - - } - fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, - p->plan, smin, smax, th, ph, &ax, &ay, &az); - - /* Search for b */ - smin = 2.0*pmax * bmin; - smax = 2.0*pmax * bmax; - fom = 0.0; th = 0.0; ph = 0.0; - for ( i=0; in_dir; i++ ) { - - double new_fom, ang; - - ang = angle_between(p->directions[i].x, p->directions[i].y, - p->directions[i].z, ax, ay, az); - if ( fabs(ang-ga) > angtol ) continue; - - new_fom = check_dir(&p->directions[i], image->features, - p->nel, pmax, fft_in, fft_out, p->plan, - smin, smax, NULL, NULL); - if ( new_fom > fom ) { - fom = new_fom; - th = p->directions[i].th; - ph = p->directions[i].ph; - } - - } - fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, - p->plan, smin, smax, th, ph, &bx, &by, &bz); - - /* Search for c */ - smin = 2.0*pmax * cmin; - smax = 2.0*pmax * cmax; - fom = 0.0; th = 0.0; ph = 0.0; - for ( i=0; in_dir; i++ ) { - - double new_fom, ang; - - ang = angle_between(p->directions[i].x, p->directions[i].y, - p->directions[i].z, ax, ay, az); - if ( fabs(ang-be) > angtol ) continue; - - ang = angle_between(p->directions[i].x, p->directions[i].y, - p->directions[i].z, bx, by, bz); - if ( fabs(ang-al) > angtol ) continue; - - new_fom = check_dir(&p->directions[i], image->features, - p->nel, pmax, fft_in, fft_out, p->plan, - smin, smax, NULL, NULL); - if ( new_fom > fom ) { - fom = new_fom; - th = p->directions[i].th; - ph = p->directions[i].ph; - } - - } - fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, - p->plan, smin, smax, th, ph, &cx, &cy, &cz); - - image->candidate_cells[0] = cell_new(); - cell_set_cartesian(image->candidate_cells[0], - ax, ay, az, bx, by, bz, cx, cy, cz); - - refine_all_rigid_groups(image, image->candidate_cells[0], p->nel, pmax, - fft_in, fft_out, p->plan, smin, smax, - image->det, p); - map_all_peaks(image); - refine_cell(image, image->candidate_cells[0], image->features); - - fftw_free(fft_in); - fftw_free(fft_out); - - image->ncells = 1; -} - - -IndexingPrivate *reax_prepare() -{ - struct reax_private *p; - int samp; - double th; - - p = calloc(1, sizeof(*p)); - if ( p == NULL ) return NULL; - - p->base.indm = INDEXING_REAX; - - p->angular_inc = deg2rad(1.7); /* From Steller (1997) */ - - /* Reserve memory, over-estimating the number of directions */ - samp = 2.0*M_PI / p->angular_inc; - p->directions = malloc(samp*samp*sizeof(struct dvec)); - if ( p == NULL) { - free(p); - return NULL; - } - STATUS("Allocated space for %i directions\n", samp*samp); - - /* Generate vectors for 1D Fourier transforms */ - fesetround(1); /* Round to nearest */ - p->n_dir = 0; - for ( th=0.0; thangular_inc ) { - - double ph, phstep, n_phstep; - - n_phstep = 2.0*M_PI*sin(th)/p->angular_inc; - n_phstep = nearbyint(n_phstep); - phstep = 2.0*M_PI/n_phstep; - - for ( ph=0.0; ph<2.0*M_PI; ph+=phstep ) { - - struct dvec *dir; - - assert(p->n_dirdirections[p->n_dir++]; - - dir->x = cos(ph) * sin(th); - dir->y = sin(ph) * sin(th); - dir->z = cos(th); - dir->th = th; - dir->ph = ph; - - } - - } - STATUS("Generated %i directions (angular increment %.3f deg)\n", - p->n_dir, rad2deg(p->angular_inc)); - - p->nel = 1024; - - /* These arrays are not actually used */ - p->fft_in = fftw_malloc(p->nel*sizeof(double)); - p->fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - p->plan = fftw_plan_dft_r2c_1d(p->nel, p->fft_in, p->fft_out, - FFTW_MEASURE); - - p->cw = 128; p->ch = 128; - - /* Also not used */ - p->r_fft_in = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); - p->r_fft_out = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); - - p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, - 1, FFTW_MEASURE); - - return (IndexingPrivate *)p; -} - - -void reax_cleanup(IndexingPrivate *pp) -{ - struct reax_private *p; - - assert(pp->indm == INDEXING_REAX); - p = (struct reax_private *)pp; - - fftw_destroy_plan(p->plan); - fftw_free(p->fft_in); - fftw_free(p->fft_out); - - fftw_destroy_plan(p->r_plan); - fftw_free(p->r_fft_in); - fftw_free(p->r_fft_out); - - free(p); -} diff --git a/src/reax.h b/src/reax.h deleted file mode 100644 index 543cd0d5..00000000 --- a/src/reax.h +++ /dev/null @@ -1,27 +0,0 @@ -/* - * reax.h - * - * A new auto-indexer - * - * (c) 2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifndef REAX_H -#define REAX_H - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include "cell.h" - -extern IndexingPrivate *reax_prepare(void); -extern void reax_cleanup(IndexingPrivate *pp); - -extern void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell); - -#endif /* REAX_H */ diff --git a/src/render.c b/src/render.c deleted file mode 100644 index dbaca2fd..00000000 --- a/src/render.c +++ /dev/null @@ -1,442 +0,0 @@ -/* - * render.c - * - * Render a high dynamic range buffer in some sensible way - * - * (c) 2006-2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#ifdef HAVE_GTK -#include -#endif - -#include -#include -#include - -#ifdef HAVE_TIFF -#include -#endif - - -#include "hdf5-file.h" -#include "render.h" -#include "peaks.h" -#include "filters.h" -#include "utils.h" - - -static void render_rgb(double val, double max, - double *rp, double *gp, double *bp) -{ - int s; - double p; - double r, g, b; - - s = val / (max/6); - p = fmod(val, max/6.0); - p /= (max/6.0); - - r = 0.0; g = 0.0; b = 0.0; - - if ( (val < 0.0) ) { - s = 0; - p = 0; - } - if ( (val > max) ) { - s = 6; - } - switch ( s ) { - case 0 : { /* Black to blue */ - r = 0; g = 0; b = p; - break; - } - case 1 : { /* Blue to pink */ - r = p; g = 0; b = 1.0; - break; - } - case 2 : { /* Pink to red */ - r = 1.0; g = 0; b = (1.0-p)*1.0; - break; - } - case 3 : { /* Red to Orange */ - r = 1.0; g = 0.5*p; b = 0; - break; - } - case 4 : { /* Orange to Yellow */ - r = 1.0; g = 0.5 + 0.5*p; b = 0; - break; - } - case 5 : { /* Yellow to White */ - r = 1.0; g = 1.0; b = 1.0*p; - break; - } - case 6 : { /* Pixel has hit the maximum value */ - r = 1.0; g = 1.0; b = 1.0; - break; - } - } - - *rp = r; - *gp = g; - *bp = b; -} - - -static void render_ratio(double val, double max, - double *rp, double *gp, double *bp) -{ - if ( val <= 1.0 ) { - render_rgb(val, 2.0, rp, gp, bp); - } else { - /* Your homework is to simplify this expression */ - val = ((val-1.0)/(max-1.0)) * (max/2.0) + max/2.0; - render_rgb(val, max, rp, gp, bp); - } -} - - -static void render_mono(double val, double max, - double *rp, double *gp, double *bp) -{ - double p; - p = val / max; - if ( val < 0.0 ) p = 0.0; - if ( val > max ) p = 1.0; - *rp = p; - *gp = p; - *bp = p; -} - - -static void render_invmono(double val, double max, - double *rp, double *gp, double *bp) -{ - double p; - p = val / max; - p = 1.0 - p; - if ( val < 0.0 ) p = 1.0; - if ( val > max ) p = 0.0; - *rp = p; - *gp = p; - *bp = p; -} - - -void render_scale(double val, double max, int scale, - double *rp, double *gp, double *bp) -{ - switch ( scale ) { - case SCALE_COLOUR : - render_rgb(val, max, rp, gp, bp); - break; - case SCALE_MONO : - render_mono(val, max, rp, gp, bp); - break; - case SCALE_INVMONO : - render_invmono(val, max, rp, gp, bp); - break; - case SCALE_RATIO : - render_ratio(val, max, rp, gp, bp); - break; - } -} - - -#ifdef HAVE_GTK - -static float *get_binned_panel(struct image *image, int binning, - int min_fs, int max_fs, int min_ss, int max_ss) -{ - float *data; - int x, y; - int w, h; - int fw; - float *in; - - fw = image->width; - in = image->data; - - /* Some pixels might get discarded */ - w = (max_fs - min_fs + 1) / binning; - h = (max_ss - min_ss + 1) / binning; - - data = malloc(w*h*sizeof(float)); - - for ( x=0; xwidth*image->height; i++ ) { - if ( image->data[i] > max ) max = image->data[i]; - } - hdr = get_binned_panel(image, binning, min_fs, max_fs, min_ss, max_ss); - if ( hdr == NULL ) return NULL; - - /* Rendered (colourful) version */ - data = malloc(3*w*h); - if ( data == NULL ) { - free(hdr); - return NULL; - } - - max /= boost; - if ( max <= 6 ) { max = 10; } - /* These x,y coordinates are measured relative to the bottom-left - * corner */ - for ( y=0; ydet->n_panels; - GdkPixbuf **pixbufs; - - pixbufs = calloc(np, sizeof(GdkPixbuf*)); - if ( pixbufs == NULL ) { - *n_pixbufs = 0; - return NULL; - } - - for ( i=0; idet->panels[i].min_fs, - image->det->panels[i].max_fs, - image->det->panels[i].min_ss, - image->det->panels[i].max_ss); - - } - - *n_pixbufs = np; - - return pixbufs; -} - - -GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale) -{ - guchar *data; - size_t x, y; - int max; - - data = malloc(3*w*h); - if ( data == NULL ) return NULL; - - max = h; - - for ( y=0; ywidth); - TIFFSetField(th, TIFFTAG_IMAGELENGTH, image->height); - TIFFSetField(th, TIFFTAG_SAMPLESPERPIXEL, 1); - TIFFSetField(th, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP); - TIFFSetField(th, TIFFTAG_BITSPERSAMPLE, 32); - TIFFSetField(th, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); - TIFFSetField(th, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); - TIFFSetField(th, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); - TIFFSetField(th, TIFFTAG_ROWSPERSTRIP, - TIFFDefaultStripSize(th, image->width*4)); - - line = _TIFFmalloc(TIFFScanlineSize(th)); - for ( y=0; yheight; y++ ) { - memcpy(line, &image->data[(image->height-1-y)*image->width], - image->width*4); - TIFFWriteScanline(th, line, y, 0); - } - _TIFFfree(line); - - TIFFClose(th); - -#else - STATUS("No TIFF support.\n"); -#endif - return 0; -} - - -int render_tiff_int16(struct image *image, const char *filename, double boost) -{ -#ifdef HAVE_TIFF - TIFF *th; - int16_t *line; - int x, y; - double max; - - th = TIFFOpen(filename, "w"); - if ( th == NULL ) return 1; - - TIFFSetField(th, TIFFTAG_IMAGEWIDTH, image->width); - TIFFSetField(th, TIFFTAG_IMAGELENGTH, image->height); - TIFFSetField(th, TIFFTAG_SAMPLESPERPIXEL, 1); - TIFFSetField(th, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_INT); /* (signed) */ - TIFFSetField(th, TIFFTAG_BITSPERSAMPLE, 16); - TIFFSetField(th, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); - TIFFSetField(th, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); - TIFFSetField(th, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); - TIFFSetField(th, TIFFTAG_ROWSPERSTRIP, - TIFFDefaultStripSize(th, image->width*4)); - - line = _TIFFmalloc(TIFFScanlineSize(th)); - max = 0.0; - for ( y=0; yheight; y++ ) { - for ( x=0;xwidth; x++ ) { - double val; - val = image->data[x+image->height*y]; - if ( val > max ) max = val; - } - } - max /= 32767.0; - - for ( y=0; yheight; y++ ) { - for ( x=0;xwidth; x++ ) { - - double val; - - val = image->data[x+(image->height-1-y)*image->width]; - val *= ((double)boost/max); - - /* Clamp to 16-bit range, - * and work round inability of most readers to deal - * with signed integers. */ - val += 1000.0; - if ( val > 32767.0 ) val = 32767.0; - if ( val < 0.0 ) val = 0.0; - - line[x] = val; - } - - TIFFWriteScanline(th, line, y, 0); - } - _TIFFfree(line); - - TIFFClose(th); -#else - STATUS("No TIFF support.\n"); -#endif - return 0; -} - -#endif /* HAVE_GTK */ diff --git a/src/render.h b/src/render.h deleted file mode 100644 index fbfbbdb2..00000000 --- a/src/render.h +++ /dev/null @@ -1,55 +0,0 @@ -/* - * render.h - * - * Render a high dynamic range buffer in some sensible way - * - * (c) 2006-2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifdef HAVE_CONFIG_H -#include -#endif - -#ifndef RENDER_H -#define RENDER_H - - -#include - -#include "image.h" - -enum { - SCALE_COLOUR, - SCALE_MONO, - SCALE_INVMONO, - SCALE_RATIO -}; - -/* Colour scale lookup */ -extern void render_scale(double val, double max, int scale, - double *rp, double *gp, double *bp); - - -#ifdef HAVE_GTK - -#include - -extern GdkPixbuf **render_panels(struct image *image, - int binning, int scale, double boost, - int *n_pixbufs); - -extern GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale); - -extern int render_tiff_fp(struct image *image, const char *filename); -extern int render_tiff_int16(struct image *image, const char *filename, - double boost); - -#endif /* HAVE_GTK */ - - - -#endif /* RENDER_H */ -- cgit v1.2.3