diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/Makefile.am | 6 | ||||
-rw-r--r-- | src/Makefile.in | 417 | ||||
-rw-r--r-- | src/control.h | 94 | ||||
-rw-r--r-- | src/displaywindow.c | 543 | ||||
-rw-r--r-- | src/displaywindow.h | 28 | ||||
-rw-r--r-- | src/imagedisplay.c | 244 | ||||
-rw-r--r-- | src/imagedisplay.h | 44 | ||||
-rw-r--r-- | src/ipr.c | 253 | ||||
-rw-r--r-- | src/ipr.h | 23 | ||||
-rw-r--r-- | src/itrans.c | 527 | ||||
-rw-r--r-- | src/itrans.h | 24 | ||||
-rw-r--r-- | src/main.c | 183 | ||||
-rw-r--r-- | src/main.h | 19 | ||||
-rw-r--r-- | src/mrc.c | 141 | ||||
-rw-r--r-- | src/mrc.h | 121 | ||||
-rw-r--r-- | src/qdrp.c | 254 | ||||
-rw-r--r-- | src/qdrp.h | 22 | ||||
-rw-r--r-- | src/readpng.c | 146 | ||||
-rw-r--r-- | src/readpng.h | 22 | ||||
-rw-r--r-- | src/reflections.c | 186 | ||||
-rw-r--r-- | src/reflections.h | 77 | ||||
-rw-r--r-- | src/trackball.c | 324 | ||||
-rw-r--r-- | src/trackball.h | 78 | ||||
-rw-r--r-- | src/utils.c | 44 | ||||
-rw-r--r-- | src/utils.h | 24 |
25 files changed, 3844 insertions, 0 deletions
diff --git a/src/Makefile.am b/src/Makefile.am new file mode 100644 index 0000000..1e532f1 --- /dev/null +++ b/src/Makefile.am @@ -0,0 +1,6 @@ +bin_PROGRAMS = dtr +dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c +dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas +AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ +AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" + diff --git a/src/Makefile.in b/src/Makefile.in new file mode 100644 index 0000000..d07e5be --- /dev/null +++ b/src/Makefile.in @@ -0,0 +1,417 @@ +# Makefile.in generated by automake 1.9.6 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005 Free Software Foundation, Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + +@SET_MAKE@ + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ +pkgdatadir = $(datadir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +top_builddir = .. +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +INSTALL = @INSTALL@ +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +bin_PROGRAMS = dtr$(EXEEXT) +subdir = src +DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(install_sh) -d +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +am__installdirs = "$(DESTDIR)$(bindir)" +binPROGRAMS_INSTALL = $(INSTALL_PROGRAM) +PROGRAMS = $(bin_PROGRAMS) +am_dtr_OBJECTS = main.$(OBJEXT) displaywindow.$(OBJEXT) \ + trackball.$(OBJEXT) reflections.$(OBJEXT) readpng.$(OBJEXT) \ + mrc.$(OBJEXT) imagedisplay.$(OBJEXT) utils.$(OBJEXT) \ + itrans.$(OBJEXT) qdrp.$(OBJEXT) ipr.$(OBJEXT) +dtr_OBJECTS = $(am_dtr_OBJECTS) +dtr_DEPENDENCIES = +DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir) +depcomp = $(SHELL) $(top_srcdir)/depcomp +am__depfiles_maybe = depfiles +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@ +SOURCES = $(dtr_SOURCES) +DIST_SOURCES = $(dtr_SOURCES) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMDEP_FALSE = @AMDEP_FALSE@ +AMDEP_TRUE = @AMDEP_TRUE@ +AMTAR = @AMTAR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CCDEPMODE = @CCDEPMODE@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +DEPDIR = @DEPDIR@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +GREP = @GREP@ +GTKGLEXT_CFLAGS = @GTKGLEXT_CFLAGS@ +GTKGLEXT_LIBS = @GTKGLEXT_LIBS@ +GTK_CFLAGS = @GTK_CFLAGS@ +GTK_LIBS = @GTK_LIBS@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAKEINFO = @MAKEINFO@ +OBJEXT = @OBJEXT@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +PKG_CONFIG = @PKG_CONFIG@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +ac_ct_CC = @ac_ct_CC@ +am__fastdepCC_FALSE = @am__fastdepCC_FALSE@ +am__fastdepCC_TRUE = @am__fastdepCC_TRUE@ +am__include = @am__include@ +am__leading_dot = @am__leading_dot@ +am__quote = @am__quote@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build_alias = @build_alias@ +datadir = @datadir@ +datarootdir = @datarootdir@ +docdir = @docdir@ +dvidir = @dvidir@ +exec_prefix = @exec_prefix@ +host_alias = @host_alias@ +htmldir = @htmldir@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localedir = @localedir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +pdfdir = @pdfdir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +psdir = @psdir@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c +dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas +AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ +AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .o .obj +$(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \ + && exit 0; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu src/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu src/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +install-binPROGRAMS: $(bin_PROGRAMS) + @$(NORMAL_INSTALL) + test -z "$(bindir)" || $(mkdir_p) "$(DESTDIR)$(bindir)" + @list='$(bin_PROGRAMS)'; for p in $$list; do \ + p1=`echo $$p|sed 's/$(EXEEXT)$$//'`; \ + if test -f $$p \ + ; then \ + f=`echo "$$p1" | sed 's,^.*/,,;$(transform);s/$$/$(EXEEXT)/'`; \ + echo " $(INSTALL_PROGRAM_ENV) $(binPROGRAMS_INSTALL) '$$p' '$(DESTDIR)$(bindir)/$$f'"; \ + $(INSTALL_PROGRAM_ENV) $(binPROGRAMS_INSTALL) "$$p" "$(DESTDIR)$(bindir)/$$f" || exit 1; \ + else :; fi; \ + done + +uninstall-binPROGRAMS: + @$(NORMAL_UNINSTALL) + @list='$(bin_PROGRAMS)'; for p in $$list; do \ + f=`echo "$$p" | sed 's,^.*/,,;s/$(EXEEXT)$$//;$(transform);s/$$/$(EXEEXT)/'`; \ + echo " rm -f '$(DESTDIR)$(bindir)/$$f'"; \ + rm -f "$(DESTDIR)$(bindir)/$$f"; \ + done + +clean-binPROGRAMS: + -test -z "$(bin_PROGRAMS)" || rm -f $(bin_PROGRAMS) +dtr$(EXEEXT): $(dtr_OBJECTS) $(dtr_DEPENDENCIES) + @rm -f dtr$(EXEEXT) + $(LINK) $(dtr_LDFLAGS) $(dtr_OBJECTS) $(dtr_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/displaywindow.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/imagedisplay.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ipr.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/itrans.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/main.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mrc.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/qdrp.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/readpng.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/reflections.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/trackball.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/utils.Po@am__quote@ + +.c.o: +@am__fastdepCC_TRUE@ if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ $<; \ +@am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/$*.Tpo" "$(DEPDIR)/$*.Po"; else rm -f "$(DEPDIR)/$*.Tpo"; exit 1; fi +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(COMPILE) -c $< + +.c.obj: +@am__fastdepCC_TRUE@ if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ `$(CYGPATH_W) '$<'`; \ +@am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/$*.Tpo" "$(DEPDIR)/$*.Po"; else rm -f "$(DEPDIR)/$*.Tpo"; exit 1; fi +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(COMPILE) -c `$(CYGPATH_W) '$<'` +uninstall-info-am: + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$tags $$unique; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(CTAGS_ARGS)$$tags$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$tags $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && cd $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) $$here + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \ + list='$(DISTFILES)'; for file in $$list; do \ + case $$file in \ + $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \ + $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \ + esac; \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test "$$dir" != "$$file" && test "$$dir" != "."; then \ + dir="/$$dir"; \ + $(mkdir_p) "$(distdir)$$dir"; \ + else \ + dir=''; \ + fi; \ + if test -d $$d/$$file; then \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \ + fi; \ + cp -pR $$d/$$file $(distdir)$$dir || exit 1; \ + else \ + test -f $(distdir)/$$file \ + || cp -p $$d/$$file $(distdir)/$$file \ + || exit 1; \ + fi; \ + done +check-am: all-am +check: check-am +all-am: Makefile $(PROGRAMS) +installdirs: + for dir in "$(DESTDIR)$(bindir)"; do \ + test -z "$$dir" || $(mkdir_p) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-binPROGRAMS clean-generic mostlyclean-am + +distclean: distclean-am + -rm -rf ./$(DEPDIR) + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +info: info-am + +info-am: + +install-data-am: + +install-exec-am: install-binPROGRAMS + +install-info: install-info-am + +install-man: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -rf ./$(DEPDIR) + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-binPROGRAMS uninstall-info-am + +.PHONY: CTAGS GTAGS all all-am check check-am clean clean-binPROGRAMS \ + clean-generic ctags distclean distclean-compile \ + distclean-generic distclean-tags distdir dvi dvi-am html \ + html-am info info-am install install-am install-binPROGRAMS \ + install-data install-data-am install-exec install-exec-am \ + install-info install-info-am install-man install-strip \ + installcheck installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic pdf pdf-am ps ps-am tags uninstall \ + uninstall-am uninstall-binPROGRAMS uninstall-info-am + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/src/control.h b/src/control.h new file mode 100644 index 0000000..0493757 --- /dev/null +++ b/src/control.h @@ -0,0 +1,94 @@ +/* + * control.h + * + * Common control structure + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef CONTROL_H +#define CONTROL_H + +#include <gtk/gtk.h> + +typedef enum ift_enum { + INPUT_NONE, + INPUT_QDRP, + INPUT_MRC +} InputFileType; + +typedef enum { + FORMULATION_CLEN, + FORMULATION_PIXELSIZE +} FormulationMode; + +typedef enum { + PEAKSEARCH_THRESHOLD, + PEAKSEARCH_ADAPTIVE_THRESHOLD, + PEAKSEARCH_LSQ, + PEAKSEARCH_ZAEFFERER +} PeakSearchMode; + +typedef enum { + RECONSTRUCTION_MAPPING, + RECONSTRUCTION_PREDICTION +} ReconstructionMode; + +typedef struct cctx_struct { + + /* Modes */ + InputFileType inputfiletype; + FormulationMode fmode; + ReconstructionMode rmode; + PeakSearchMode psmode; + + /* Input filename */ + char *filename; + + /* Basic parameters */ + double camera_length; + double omega; + double resolution; + double lambda; + double pixel_size; + + /* (QDRP) Parser flags */ + unsigned int started; + unsigned int camera_length_set; + unsigned int omega_set; + unsigned int max_d_set; + unsigned int resolution_set; + unsigned int lambda_set; + + /* Miscellaneous modes and operations */ + unsigned int max_d; + unsigned int first_image; + + /* Size of input images - assumed the same throughout. */ + unsigned int width; + unsigned int height; + unsigned int x_centre; + unsigned int y_centre; + + /* Information about the input images */ + unsigned int n_images; + + + /* Output */ + struct rctx_struct *reflectionctx; + + /* GTK bits */ + GtkWidget *combo_algorithm; + GtkWidget *combo_peaksearch; + +} ControlContext; + +extern ControlContext *control_read(const char *filename); + +#endif /* CONTROL_H */ diff --git a/src/displaywindow.c b/src/displaywindow.c new file mode 100644 index 0000000..58c9e10 --- /dev/null +++ b/src/displaywindow.c @@ -0,0 +1,543 @@ +/* + * displaywindow.c + * + * The display window + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#define _GNU_SOURCE +#include <string.h> +#include <gtk/gtk.h> +#include <stdlib.h> +#include <math.h> +#include <gdk/gdkgl.h> +#include <gtk/gtkgl.h> +#include <GL/gl.h> +#include <GL/glu.h> + +#include "trackball.h" +#include "reflections.h" +#include "main.h" + +enum { + DW_ORTHO, + DW_PERSPECTIVE +}; + +GtkUIManager *displaywindow_ui; +GtkActionGroup *displaywindow_action_group; +GtkWidget *displaywindow_window; +GtkWidget *displaywindow_bigvbox; +GtkWidget *displaywindow_status_bar; +int displaywindow_view = DW_ORTHO; +GtkWidget *displaywindow_drawing_area; +GLfloat displaywindow_distance = 150; +GLfloat displaywindow_x_pos = 0; +GLfloat displaywindow_y_pos = 0; + +void error_report(const char *message) { + + GtkWidget *window; + + window = gtk_message_dialog_new(GTK_WINDOW(displaywindow_window), GTK_DIALOG_DESTROY_WITH_PARENT, + GTK_MESSAGE_WARNING, GTK_BUTTONS_CLOSE, message); + + g_signal_connect_swapped(window, "response", G_CALLBACK(gtk_widget_destroy), window); + gtk_widget_show(window); + +} + +static void displaywindow_about() { + + GtkWidget *window; + + const gchar *authors[] = { + "Thomas White <taw27@cam.ac.uk>", + "Gordon Ball <gfb21@cam.ac.uk>", + "Paul Midgley <pam33@cam.ac.uk>", + NULL + }; + + window = gtk_about_dialog_new(); + + gtk_about_dialog_set_name(GTK_ABOUT_DIALOG(window), PACKAGE_NAME); + gtk_about_dialog_set_version(GTK_ABOUT_DIALOG(window), PACKAGE_VERSION); + gtk_about_dialog_set_copyright(GTK_ABOUT_DIALOG(window), "(c) 2006-2007 Thomas White <taw27@cam.ac.uk>"); + gtk_about_dialog_set_comments(GTK_ABOUT_DIALOG(window), "Diffraction Tomography Reconstruction"); + gtk_about_dialog_set_license(GTK_ABOUT_DIALOG(window), "(c) 2006-2007 Thomas White <taw27@cam.ac.uk>\n" + "Virtual trackball (c) Copyright 1993, 1994, Silicon Graphics, Inc.\n" + "\n" + "Research funded by:\n" + "The Engineering and Physical Sciences Research Council\n" + "FEI Electron Optics B.V."); + gtk_about_dialog_set_website(GTK_ABOUT_DIALOG(window), "http://www-hrem.msm.cam.ac.uk/"); + gtk_about_dialog_set_authors(GTK_ABOUT_DIALOG(window), authors); + + g_signal_connect(window, "response", G_CALLBACK(gtk_widget_destroy), NULL); + + gtk_widget_show_all(window); + +} + +static void displaywindow_close() { + gtk_exit(0); +} + +static void displaywindow_gl_set_ortho(GLfloat w, GLfloat h) { + + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + if ( w > h ) { + GLfloat aspect = w/h; + glOrtho(-aspect*(displaywindow_distance/2), aspect*(displaywindow_distance/2), -(displaywindow_distance/2), (displaywindow_distance/2), 2.0, 1000.0); + } else { + GLfloat aspect = h/w; + glOrtho(-(displaywindow_distance/2), (displaywindow_distance/2), -aspect*(displaywindow_distance/2), aspect*(displaywindow_distance/2), 2.0, 1000.0); + } + glMatrixMode(GL_MODELVIEW); + +} + +static void displaywindow_gl_set_perspective(GLfloat w, GLfloat h) { + + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + if ( w > h ) { + GLfloat aspect = w/h; + glFrustum(-aspect, aspect, -1.0, 1.0, 2.0, 1000.0); + } else { + GLfloat aspect = h/w; + glFrustum(-1.0, 1.0, -aspect, aspect, 2.0, 1000.0); + } + glMatrixMode(GL_MODELVIEW); + +} + +static gint displaywindow_changeview(GtkWidget *widget, GtkRadioAction *action) { + + GLfloat w = displaywindow_drawing_area->allocation.width; + GLfloat h = displaywindow_drawing_area->allocation.height; + displaywindow_view = gtk_radio_action_get_current_value(action); + + if ( displaywindow_view == DW_ORTHO ) { + displaywindow_gl_set_ortho(w, h); + } else { + displaywindow_gl_set_perspective(w, h); + } + + return 0; + +} + +static void displaywindow_addui_callback(GtkUIManager *ui, GtkWidget *widget, GtkContainer *container) { + + gtk_box_pack_start(GTK_BOX(container), widget, FALSE, FALSE, 0); + + /* Enable overflow menu if this is a toolbar */ + if ( GTK_IS_TOOLBAR(widget) ) { + gtk_toolbar_set_show_arrow(GTK_TOOLBAR(widget), TRUE); + } + +} + +static void displaywindow_addmenubar() { + + GtkActionEntry entries[] = { + + { "FileAction", NULL, "_File", NULL, NULL, NULL }, + { "SaveAction", GTK_STOCK_SAVE, "_Save Image...", NULL, NULL, G_CALLBACK(NULL) }, + { "CloseAction", GTK_STOCK_CLOSE, "_Close", NULL, NULL, G_CALLBACK(displaywindow_close) }, + + { "ViewAction", NULL, "_View", NULL, NULL, NULL }, + + { "ToolsAction", NULL, "_Tools", NULL, NULL, NULL }, + + { "HelpAction", NULL, "_Help", NULL, NULL, NULL }, + { "AboutAction", GTK_STOCK_ABOUT, "_About DTR...", NULL, NULL, G_CALLBACK(displaywindow_about) }, + + }; + guint n_entries = G_N_ELEMENTS(entries); + GtkRadioActionEntry toggles[] = { + { "OrthoAction", NULL, "_Orthographic", NULL, NULL, DW_ORTHO }, + { "PerspectiveAction", NULL, "_Perspective", NULL, NULL, DW_PERSPECTIVE }, + }; + guint n_toggles = G_N_ELEMENTS(toggles); + + GError *error = NULL; + + displaywindow_action_group = gtk_action_group_new("dtrdisplaywindow"); + gtk_action_group_add_actions(displaywindow_action_group, entries, n_entries, displaywindow_window); + gtk_action_group_add_radio_actions(displaywindow_action_group, toggles, n_toggles, -1, G_CALLBACK(displaywindow_changeview), NULL); + + displaywindow_ui = gtk_ui_manager_new(); + gtk_ui_manager_insert_action_group(displaywindow_ui, displaywindow_action_group, 0); + g_signal_connect(displaywindow_ui, "add_widget", G_CALLBACK(displaywindow_addui_callback), displaywindow_bigvbox); + if ( gtk_ui_manager_add_ui_from_file(displaywindow_ui, DATADIR"/dtr/displaywindow.ui", &error) == 0 ) { + fprintf(stderr, "Error loading message window menu bar: %s\n", error->message); + return; + } + + gtk_window_add_accel_group(GTK_WINDOW(displaywindow_window), gtk_ui_manager_get_accel_group(displaywindow_ui)); + gtk_ui_manager_ensure_update(displaywindow_ui); + +} + +static float displaywindow_x_start = 0.0; +static float displaywindow_y_start = 0.0; + +static float view_quat[4] = { 0.0, 0.0, 0.0, 1.0 }; + +static gboolean displaywindow_gl_button_press(GtkWidget *widget, GdkEventButton *event, gpointer data) { + displaywindow_x_start = event->x; + displaywindow_y_start = event->y; + return FALSE; +} + +static gint displaywindow_gl_motion_notify(GtkWidget *widget, GdkEventMotion *event, gpointer data) { + + float w = widget->allocation.width; + float h = widget->allocation.height; + float x = event->x; + float y = event->y; + float d_quat[4]; + + if ( event->state & GDK_CONTROL_MASK ) { + + /* Control-click changes 'zoom' */ + displaywindow_distance += event->y - displaywindow_y_start; + if ( displaywindow_view == DW_ORTHO ) { + displaywindow_gl_set_ortho(w, h); + } + + } else if ( event->state & GDK_SHIFT_MASK ) { + + /* Shift-click translates */ + displaywindow_x_pos += (event->x - displaywindow_x_start)/5; + displaywindow_y_pos += (event->y - displaywindow_y_start)/5; + + } else { + + /* Click rotates */ + trackball(d_quat, (2.0*displaywindow_x_start - w)/w, (h-2.0*displaywindow_y_start)/h, (2.0*x-w)/w, (h-2.0*y)/h); + add_quats(d_quat, view_quat, view_quat); + + } + + displaywindow_x_start = x; + displaywindow_y_start = y; + + gdk_window_invalidate_rect(widget->window, &widget->allocation, FALSE); + return TRUE; +} + +static gint displaywindow_gl_expose(GtkWidget *widget, GdkEventExpose *event, ControlContext *ctx) { + + GdkGLContext *glcontext = gtk_widget_get_gl_context(widget); + GdkGLDrawable *gldrawable = gtk_widget_get_gl_drawable(widget); + GLfloat light0_position[] = { 0.0, 0.0, 100.0, 0.0 }; + GLfloat light0_diffuse[] = { 1.0, 1.0, 1.0, 1.0 }; + GLfloat light0_specular[] = { 0.5, 0.5, 0.5, 1.0 }; + GLfloat light1_position[] = { 0.0, 0.0, -100.0, 0.0 }; + GLfloat light1_diffuse[] = { 0.8, 0.8, 0.8, 1.0 }; + GLfloat light1_specular[] = { 0.5, 0.5, 0.5, 1.0 }; + GLfloat green[] = { 0.0, 1.0, 0.0, 1.0 }; + GLfloat blue[] = { 0.0, 0.0, 1.0, 1.0 }; + GLfloat red[] = { 1.0, 0.0, 0.0, 1.0 }; + GLfloat gold[] = { 0.7, 0.7, 0.0, 1.0 }; + GLfloat yellow[] = { 1.0, 1.0, 0.0, 1.0 }; + GLfloat yellow_glass[] = { 1.0, 1.0, 0.0, 000.1 }; + float m[4][4]; +// GLfloat fog_density = 0.03; + Reflection *reflection; + GLUquadricObj *quadric; + + if ( !gdk_gl_drawable_gl_begin(gldrawable, glcontext) ) { + return 0; + } + + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + + glEnable(GL_DEPTH_TEST); +// glEnable(GL_FOG); +// glFogfv(GL_FOG_DENSITY, &fog_density); + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + + glLoadIdentity(); + glTranslatef(displaywindow_x_pos, -displaywindow_y_pos, -displaywindow_distance); + build_rotmatrix(m, view_quat); + glMultMatrixf(&m[0][0]); + + /* Set up lighting */ + glEnable(GL_LIGHTING); + glEnable(GL_LIGHT0); + glLightfv(GL_LIGHT0, GL_POSITION, light0_position); + glLightfv(GL_LIGHT0, GL_DIFFUSE, light0_diffuse); + glLightfv(GL_LIGHT0, GL_SPECULAR, light0_specular); + glEnable(GL_LIGHT1); + glLightfv(GL_LIGHT1, GL_POSITION, light1_position); + glLightfv(GL_LIGHT1, GL_DIFFUSE, light1_diffuse); + glLightfv(GL_LIGHT1, GL_SPECULAR, light1_specular); + + /* Bounding cube: 100 nm^-1 side length */ + glBegin(GL_LINE_LOOP); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red); + glVertex3f(50, 50, 50); + glVertex3f(-50, 50, 50); + glVertex3f(-50, -50, 50); + glVertex3f(50, -50, 50); + glEnd(); + glBegin(GL_LINE_LOOP); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red); + glVertex3f(50, 50, -50); + glVertex3f(-50, 50, -50); + glVertex3f(-50, -50, -50); + glVertex3f(50, -50, -50); + glEnd(); + glBegin(GL_LINES); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red); + glVertex3f(50, 50, 50); + glVertex3f(50, 50, -50); + glVertex3f(-50, 50, 50); + glVertex3f(-50, 50, -50); + glVertex3f(-50, -50, 50); + glVertex3f(-50, -50, -50); + glVertex3f(50, -50, 50); + glVertex3f(50, -50, -50); + glEnd(); + + /* Tilt axis */ + glBegin(GL_LINES); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow); + glVertex3f(50, 0.0, 0.0); + glVertex3f(-50, 0.0, 0.0); + glEnd(); + glBegin(GL_TRIANGLE_FAN); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow); + glVertex3f(50, 0.0, 0.0); + glVertex3f(50-5, 1.0, 1.0); + glVertex3f(50-5, -1.0, 1.0); + glVertex3f(50-5, -1.0, -1.0); + glVertex3f(50-5, 1.0, -1.0); + glVertex3f(50-5, 1.0, 1.0); + glEnd(); + + /* Plot all the normal (measured) reflections */ + reflection = ctx->reflectionctx->reflections; + glBegin(GL_POINTS); + do { + if ( reflection->type == REFLECTION_NORMAL ) { + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green); + glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9); + } + reflection = reflection->next; + } while ( reflection != NULL ); + glEnd(); + + /* Plot the other reflections */ + reflection = ctx->reflectionctx->reflections; + quadric = gluNewQuadric(); + do { + if ( reflection->type == REFLECTION_CENTRAL ) { + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue); + glPushMatrix(); + glTranslatef(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9); + gluSphere(quadric, 0.2, 32, 32); + glPopMatrix(); + } + if ( reflection->type == REFLECTION_GENERATED ) { + + /* Generated reflections are displayed by index, not coordinates. + xyz field in Reflection will be used later for "measured" location. */ + double x, y, z; + signed int h, k, l; + + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, gold); + glPushMatrix(); + + /* Change hkl to xyz */ + h = reflection->h; k = reflection->k; l = reflection->l; + x = h*ctx->reflectionctx->basis->a.x + k*ctx->reflectionctx->basis->b.x + l*ctx->reflectionctx->basis->c.x; + y = h*ctx->reflectionctx->basis->a.y + k*ctx->reflectionctx->basis->b.y + l*ctx->reflectionctx->basis->c.y; + z = h*ctx->reflectionctx->basis->a.z + k*ctx->reflectionctx->basis->b.z + l*ctx->reflectionctx->basis->c.z; + glTranslatef(x/1e9, y/1e9, z/1e9); + + gluSphere(quadric, 0.2, 32, 32); + glPopMatrix(); + + } + reflection = reflection->next; + } while ( reflection != NULL ); + + /* If this is an iterative prediction-refinement reconstruction, draw the unit cell */ + if ( ctx->rmode == RECONSTRUCTION_PREDICTION ) { + glBegin(GL_LINE_STRIP); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue); + glVertex3f(0.0, 0.0, 0.0); + glVertex3f(ctx->reflectionctx->basis->a.x/1e9, ctx->reflectionctx->basis->a.y/1e9, ctx->reflectionctx->basis->a.z/1e9); + glVertex3f(ctx->reflectionctx->basis->a.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->a.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->a.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->b.x/1e9, ctx->reflectionctx->basis->b.y/1e9, ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(0.0, 0.0, 0.0); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9, ctx->reflectionctx->basis->c.y/1e9, ctx->reflectionctx->basis->c.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->a.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->a.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->a.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->a.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->a.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->a.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9, ctx->reflectionctx->basis->c.y/1e9, ctx->reflectionctx->basis->c.z/1e9); + + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->b.x/1e9, ctx->reflectionctx->basis->b.y/1e9, ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->a.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->a.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->a.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->a.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->a.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->a.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->a.x/1e9 + ctx->reflectionctx->basis->b.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->a.y/1e9 + ctx->reflectionctx->basis->b.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->a.z/1e9 + ctx->reflectionctx->basis->b.z/1e9); + glVertex3f(ctx->reflectionctx->basis->c.x/1e9 + ctx->reflectionctx->basis->a.x/1e9, + ctx->reflectionctx->basis->c.y/1e9 + ctx->reflectionctx->basis->a.y/1e9, + ctx->reflectionctx->basis->c.z/1e9 + ctx->reflectionctx->basis->a.z/1e9); + glVertex3f(ctx->reflectionctx->basis->a.x/1e9, ctx->reflectionctx->basis->a.y/1e9, ctx->reflectionctx->basis->a.z/1e9); + glEnd(); + } + + /* Zero plane */ + glBegin(GL_QUADS); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow_glass); + glVertex3f(50, 50, 0.0); + glVertex3f(50, -50, 0.0); + glVertex3f(-50, -50, 0.0); + glVertex3f(-50, 50, 0.0); + glEnd(); + + if ( gdk_gl_drawable_is_double_buffered(gldrawable) ) { + gdk_gl_drawable_swap_buffers(gldrawable); + } else { + glFlush(); + } + + gdk_gl_drawable_gl_end(gldrawable); + + return 0; + +} + +static gint displaywindow_gl_realise(GtkWidget *widget, gpointer data) { + + GdkGLContext *glcontext = gtk_widget_get_gl_context(widget); + GdkGLDrawable *gldrawable = gtk_widget_get_gl_drawable(widget); + GLfloat w = widget->allocation.width; + GLfloat h = widget->allocation.height; + + if ( !gdk_gl_drawable_gl_begin(gldrawable, glcontext) ) { + return 0; + } + + glClearColor(0.002, 0.0, 0.1, 1.0); + glClearDepth(1.0); + displaywindow_gl_set_ortho(w, h); + + gdk_gl_drawable_gl_end(gldrawable); + + return 0; + +} + +static gboolean displaywindow_gl_configure(GtkWidget *widget, GdkEventConfigure *event, gpointer data) { + + GdkGLContext *glcontext = gtk_widget_get_gl_context(widget); + GdkGLDrawable *gldrawable = gtk_widget_get_gl_drawable(widget); + GLfloat w = widget->allocation.width; + GLfloat h = widget->allocation.height; + + /* Set viewport */ + if ( !gdk_gl_drawable_gl_begin(gldrawable, glcontext) ) { + return FALSE; + } + glViewport(0, 0, w, h); + gdk_gl_drawable_gl_end(gldrawable); + + /* Nudge the projection matrix routines to preserve the aspect ratio */ + if ( displaywindow_view == DW_ORTHO ) { + displaywindow_gl_set_ortho(w, h); + } else { + displaywindow_gl_set_perspective(w, h); + } + + return FALSE; + +} + +static gint displaywindow_closedown(GtkWidget *widget, gpointer data) { + + return 0; + +} + +void displaywindow_open(ControlContext *ctx) { + + const char *filename; + char *title; + GdkGLConfig *glconfig; + + filename = basename(ctx->filename); + title = malloc(10+strlen(filename)); + strcpy(title, "dtr: "); + strcat(title, filename); + + displaywindow_window = gtk_window_new(GTK_WINDOW_TOPLEVEL); + gtk_window_set_title(GTK_WINDOW(displaywindow_window), title); + free(title); + displaywindow_bigvbox = gtk_vbox_new(FALSE, FALSE); + gtk_container_add(GTK_CONTAINER(displaywindow_window), displaywindow_bigvbox); + displaywindow_addmenubar(); + + displaywindow_status_bar = gtk_statusbar_new(); + gtk_statusbar_set_has_resize_grip(GTK_STATUSBAR(displaywindow_status_bar), FALSE); + gtk_box_pack_end(GTK_BOX(displaywindow_bigvbox), displaywindow_status_bar, FALSE, FALSE, 0); + + g_signal_connect(GTK_OBJECT(displaywindow_window), "destroy", G_CALLBACK(displaywindow_closedown), NULL); + g_signal_connect_after(GTK_OBJECT(displaywindow_window), "destroy", G_CALLBACK(gtk_main_quit), NULL); + + /* GL stuff */ + glconfig = gdk_gl_config_new_by_mode(GDK_GL_MODE_RGB | GDK_GL_MODE_DEPTH | GDK_GL_MODE_DOUBLE); + if ( glconfig == NULL ) { + fprintf(stderr, "Can't find double-buffered visual.\n"); + exit(1); + } + gtk_container_set_reallocate_redraws(GTK_CONTAINER(displaywindow_window), TRUE); + displaywindow_drawing_area = gtk_drawing_area_new(); + gtk_widget_set_size_request(displaywindow_drawing_area, 640, 640); + gtk_widget_set_gl_capability(displaywindow_drawing_area, glconfig, NULL, TRUE, GDK_GL_RGBA_TYPE); + gtk_box_pack_start(GTK_BOX(displaywindow_bigvbox), displaywindow_drawing_area, TRUE, TRUE, 0); + gtk_widget_add_events(displaywindow_drawing_area, GDK_BUTTON1_MOTION_MASK | GDK_BUTTON_PRESS_MASK | GDK_VISIBILITY_NOTIFY_MASK); + g_signal_connect(GTK_OBJECT(displaywindow_drawing_area), "configure_event", G_CALLBACK(displaywindow_gl_configure), NULL); + g_signal_connect(GTK_OBJECT(displaywindow_drawing_area), "realize", G_CALLBACK(displaywindow_gl_realise), NULL); + g_signal_connect(GTK_OBJECT(displaywindow_drawing_area), "expose_event", G_CALLBACK(displaywindow_gl_expose), ctx); + g_signal_connect(GTK_OBJECT(displaywindow_drawing_area), "button_press_event", G_CALLBACK(displaywindow_gl_button_press), NULL); + g_signal_connect(GTK_OBJECT(displaywindow_drawing_area), "motion_notify_event", G_CALLBACK(displaywindow_gl_motion_notify), NULL); + + gtk_widget_show_all(displaywindow_window); + +} diff --git a/src/displaywindow.h b/src/displaywindow.h new file mode 100644 index 0000000..6eae540 --- /dev/null +++ b/src/displaywindow.h @@ -0,0 +1,28 @@ +/* + * displaywindow.h + * + * The display window + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef DISPLAYWINDOW_H +#define DISPLAYWINDOW_H + +#include <gtk/gtk.h> + +#include "control.h" + +extern GtkWidget *displaywindow_gtkwindow(void); +extern void displaywindow_open(ControlContext *ctx); +extern void displaywindow_statusbar(const char *message); +extern void error_report(const char *message); + +#endif /* DISPLAYWINDOW_H */ diff --git a/src/imagedisplay.c b/src/imagedisplay.c new file mode 100644 index 0000000..aaa9bf9 --- /dev/null +++ b/src/imagedisplay.c @@ -0,0 +1,244 @@ +/* + * imagedisplay.c + * + * Show raw and processed images + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <gtk/gtk.h> +#include <gdk-pixbuf/gdk-pixbuf.h> +#include <stdint.h> +#include <stdlib.h> +#include <string.h> +#include <stdio.h> +#include <math.h> + +#include "imagedisplay.h" +#include "utils.h" + +/* Free pixbuf data when reference count drops to zero */ +static void imagedisplay_free_data(guchar *image_eightbit, ImageDisplay *imagedisplay) { + free(image_eightbit); +} + +static void imagedisplay_update(ImageDisplay *imagedisplay) { + + unsigned int w, h, v_w, v_h; + float aspect_image, aspect_window; + + w = imagedisplay->width; + h = imagedisplay->height; + v_w = imagedisplay->view_width; + v_h = imagedisplay->view_height; + + /* Preserve aspect ratio */ + aspect_image = (float)w/h; + aspect_window = (float)v_w/v_h; + if ( aspect_window > aspect_image ) { + v_w = aspect_image*v_h; + } else { + v_h = v_w/aspect_image; + } + + /* Create the scaled pixbuf from the 8-bit display data */ + imagedisplay->pixbuf_scaled = gdk_pixbuf_scale_simple(imagedisplay->pixbuf, v_w, v_h, GDK_INTERP_BILINEAR); + + /* Create the image with the scaled pixbuf */ + if ( !imagedisplay->image ) { + imagedisplay->image = gtk_image_new_from_pixbuf(imagedisplay->pixbuf_scaled); + gtk_container_add(GTK_CONTAINER(imagedisplay->window), imagedisplay->image); + gtk_widget_show(imagedisplay->image); + } else { + gtk_image_set_from_pixbuf(GTK_IMAGE(imagedisplay->image), GDK_PIXBUF(imagedisplay->pixbuf_scaled)); + } + +} + +static gboolean imagedisplay_configure_event(GtkWidget *widget, GdkEventConfigure *event, ImageDisplay *imagedisplay) { + + imagedisplay->view_width = event->width; + imagedisplay->view_height = event->height; + imagedisplay_update(imagedisplay); + + return FALSE; + +} + +static void imagedisplay_put_data(ImageDisplay *imagedisplay, int16_t *image16) { + + unsigned int x, y; + unsigned int w, h; + + h = imagedisplay->height; + w = imagedisplay->width; + + /* Turn 16-bit image data into 8-bit display data */ + imagedisplay->data = malloc(3*w*h); + for ( y=0; y<h; y++ ) { + for ( x=0; x<w; x++ ) { + + int16_t val16 = image16[x+w*y]; + if ( val16 > 255 ) { + imagedisplay->data[3*( x+w*(h-1-y) )] = 255; + imagedisplay->data[3*( x+w*(h-1-y) )+1] = 255; + imagedisplay->data[3*( x+w*(h-1-y) )+2] = 255; + } else if ( val16 < 0 ) { + imagedisplay->data[3*( x+w*(h-1-y) )] = 0; + imagedisplay->data[3*( x+w*(h-1-y) )+1] = 0; + imagedisplay->data[3*( x+w*(h-1-y) )+2] = 0; + } else { + imagedisplay->data[3*( x+w*(h-1-y) )] = val16; + imagedisplay->data[3*( x+w*(h-1-y) )+1] = val16; + imagedisplay->data[3*( x+w*(h-1-y) )+2] = val16; + } + + } + } + + /* Create the pixbuf from the 8-bit display data */ + imagedisplay->pixbuf = gdk_pixbuf_new_from_data(imagedisplay->data, GDK_COLORSPACE_RGB, FALSE, 8, w, h, h*3, + (GdkPixbufDestroyNotify)imagedisplay_free_data, imagedisplay); + + imagedisplay_update(imagedisplay); + +} + +static void imagedisplay_close(GtkWidget *widget, ImageDisplay *imagedisplay) { + free(imagedisplay); +} + +/* Display an image */ +ImageDisplay *imagedisplay_open(int16_t *image, unsigned int width, unsigned int height, const char *title) { + + ImageDisplay *imagedisplay; + GdkGeometry geom; + + imagedisplay = malloc(sizeof(ImageDisplay)); + imagedisplay->image = NULL; + imagedisplay->width = width; + imagedisplay->height = height; + imagedisplay->view_width = 512; + imagedisplay->view_height = 512; + imagedisplay->title = strdup(title); + + imagedisplay->window = gtk_window_new(GTK_WINDOW_TOPLEVEL); + gtk_window_set_title(GTK_WINDOW(imagedisplay->window), imagedisplay->title); + + imagedisplay_put_data(imagedisplay, image); + + g_signal_connect(GTK_OBJECT(imagedisplay->window), "destroy", G_CALLBACK(imagedisplay_close), imagedisplay); + g_signal_connect(GTK_OBJECT(imagedisplay->window), "configure_event", G_CALLBACK(imagedisplay_configure_event), imagedisplay); + + geom.min_width = 128; + geom.min_height = 128; + gtk_window_set_geometry_hints(GTK_WINDOW(imagedisplay->window), GTK_WIDGET(imagedisplay->image), &geom, GDK_HINT_MIN_SIZE); + + gtk_window_set_default_size(GTK_WINDOW(imagedisplay->window), 512, 512); + + gtk_widget_show_all(imagedisplay->window); + + return imagedisplay; + +} + +void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx, double omega) { + + guchar *image_eightbit; + int w, h; + double gradient; + + if ( !imagedisplay->pixbuf ) return; + + w = imagedisplay->width; + h = imagedisplay->height; /* Size of pixbuf */ + + g_object_get(G_OBJECT(imagedisplay->pixbuf), "pixels", &image_eightbit, NULL); + + gradient = tan(M_PI*omega/180); + + if ( gradient > 1 ) { + double xs; + signed int x, y; + gradient = 1/gradient; + /* Start at the centre and draw a line out in each direction until it hits an edge. + This makes the whole thing a lot easier. */ + xs = ctx->x_centre; y = ctx->y_centre; + do { + x = xs; + image_eightbit[3*(x+w*(h-1-y))+0] = 255; + image_eightbit[3*(x+w*(h-1-y))+1] = 255; + image_eightbit[3*(x+w*(h-1-y))+2] = 0; + y++; + xs += gradient; + } while ( (xs<w) && (y<h) && (xs>=0) && (y>=0) ); + xs = ctx->x_centre; y = ctx->y_centre; + do { + x = xs; + image_eightbit[3*(x+w*(h-1-y))+0] = 255; + image_eightbit[3*(x+w*(h-1-y))+1] = 255; + image_eightbit[3*(x+w*(h-1-y))+2] = 0; + y--; + xs -= gradient; + } while ( (xs<w) && (y<h) && (xs>=0) && (y>=0) ); + } else { + double ys; + signed int x, y; + x = ctx->x_centre; ys = ctx->y_centre; + do { + y = ys; + image_eightbit[3*(x+w*(h-1-y))+0] = 255; + image_eightbit[3*(x+w*(h-1-y))+1] = 255; + image_eightbit[3*(x+w*(h-1-y))+2] = 0; + x++; + ys += gradient; + } while ( (x<w) && (ys<h) && (x>=0) && (ys>=0) ); + + x = ctx->x_centre; ys = ctx->y_centre; + do { + y = ys; + image_eightbit[3*(x+w*(h-1-y))+0] = 255; + image_eightbit[3*(x+w*(h-1-y))+1] = 255; + image_eightbit[3*(x+w*(h-1-y))+2] = 0; + x--; + ys -= gradient; + } while ( (x<w) && (ys<h) && (x>=0) && (ys>=0) ); + } + +} + +void imagedisplay_mark_point(ImageDisplay *imagedisplay, unsigned int x, unsigned int y) { + + guchar *image_eightbit; + int xd, yd; + int w, h; + + if ( !imagedisplay->pixbuf ) return; + + w = imagedisplay->width; + h = imagedisplay->height; + + g_object_get(G_OBJECT(imagedisplay->pixbuf), "pixels", &image_eightbit, NULL); + + if ( (y >= 0) && (y < h) ) { + for ( xd=biggest(x-3, 0); xd<=smallest(x+3, w-1); xd++ ) { + imagedisplay->data[3*( xd+w*(h-1-y) )] = 255; + imagedisplay->data[3*( xd+w*(h-1-y) )+1] = 0; + imagedisplay->data[3*( xd+w*(h-1-y) )+2] = 0; + } + } + if ( (x >= 0) && (x < w) ) { + for ( yd=biggest(y-3, 0); yd<=smallest(y+3, h-1); yd++ ) { + imagedisplay->data[3*( x+w*(h-1-yd) )] = 255; + imagedisplay->data[3*( x+w*(h-1-yd) )+1] = 0; + imagedisplay->data[3*( x+w*(h-1-yd) )+2] = 0; + } + } + +} diff --git a/src/imagedisplay.h b/src/imagedisplay.h new file mode 100644 index 0000000..0aa140f --- /dev/null +++ b/src/imagedisplay.h @@ -0,0 +1,44 @@ +/* + * imagedisplay.h + * + * Show raw and processed images + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef IMAGEDISPLAY_H +#define IMAGEDISPLAY_H + +#include <stdint.h> +#include <gtk/gtk.h> + +typedef struct struct_imagedisplay { + + unsigned int width; + unsigned int height; /* Of underlying image */ + char *title; + + guchar *data; + + GtkWidget *window; + GdkPixbuf *pixbuf; + GdkPixbuf *pixbuf_scaled; + GtkWidget *image; + + unsigned int view_width; + unsigned int view_height; /* Of window */ + +} ImageDisplay; + +extern ImageDisplay *imagedisplay_open(int16_t *image, unsigned int width, unsigned int height, const char *title); +extern void imagedisplay_mark_point(ImageDisplay *imagedisplay, unsigned int x, unsigned int y); +#include "control.h" +extern void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx, double omega); + +#endif /* IMAGEDISPLAY_H */ diff --git a/src/ipr.c b/src/ipr.c new file mode 100644 index 0000000..365ff75 --- /dev/null +++ b/src/ipr.c @@ -0,0 +1,253 @@ +/* + * ipr.c + * + * Iterative Prediction-Refinement Reconstruction + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> + +#include "control.h" +#include "reflections.h" +#include "itrans.h" +#include "utils.h" + +static int ipr_choose_max(Basis *basis) { + + if ( (basis->a.modulus >= basis->b.modulus) && (basis->a.modulus >= basis->c.modulus) ) return 1; + if ( (basis->b.modulus >= basis->a.modulus) && (basis->b.modulus >= basis->c.modulus) ) return 2; + if ( (basis->c.modulus >= basis->a.modulus) && (basis->c.modulus >= basis->b.modulus) ) return 3; + + /* Never reaches here... */ + fprintf(stderr, "IP: Couldn't determine maximum basis vector modulus\n"); + abort(); + +} + +/* Try to turn the basis into a more sensible one */ +static void ipr_try_compact(Basis *basis) { + + int i; + + for ( i=1; i<5; i++ ) { + + if ( modulus(basis->a.x - basis->b.x, basis->a.y - basis->b.y, basis->a.z - basis->b.z) < basis->a.modulus ) { + basis->a.x = basis->a.x - basis->b.x; basis->a.y = basis->a.y - basis->b.y; basis->a.y = basis->a.y - basis->b.y; + basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z); + } + if ( modulus(basis->a.x - basis->c.x, basis->a.z - basis->c.z, basis->a.z - basis->c.z) < basis->a.modulus ) { + basis->a.x = basis->a.x - basis->c.x; basis->a.y = basis->a.y - basis->c.y; basis->a.y = basis->a.y - basis->c.y; + basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z); + } + if ( modulus(basis->a.x + basis->b.x, basis->a.y + basis->b.y, basis->a.z + basis->b.z) < basis->a.modulus ) { + basis->a.x = basis->a.x + basis->b.x; basis->a.y = basis->a.y + basis->b.y; basis->a.y = basis->a.y + basis->b.y; + basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z); + } + if ( modulus(basis->a.x + basis->c.x, basis->a.z + basis->c.z, basis->a.z + basis->c.z) < basis->a.modulus ) { + basis->a.x = basis->a.x + basis->c.x; basis->a.y = basis->a.y + basis->c.y; basis->a.y = basis->a.y + basis->c.y; + basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z); + } + + if ( modulus(basis->b.x - basis->a.x, basis->b.y - basis->a.y, basis->b.z - basis->a.z) < basis->b.modulus ) { + basis->b.x = basis->b.x - basis->a.x; basis->b.y = basis->b.y - basis->a.y; basis->b.y = basis->b.y - basis->a.y; + basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z); + } + if ( modulus(basis->b.x - basis->c.x, basis->b.y - basis->c.y, basis->b.z - basis->c.z) < basis->b.modulus ) { + basis->b.x = basis->b.x - basis->c.x; basis->b.y = basis->b.y - basis->c.y; basis->b.y = basis->b.y - basis->c.y; + basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z); + } + if ( modulus(basis->b.x + basis->a.x, basis->b.y + basis->a.y, basis->b.z + basis->a.z) < basis->b.modulus ) { + basis->b.x = basis->b.x + basis->a.x; basis->b.y = basis->b.y + basis->a.y; basis->b.y = basis->b.y + basis->a.y; + basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z); + } + if ( modulus(basis->b.x + basis->c.x, basis->b.y + basis->c.y, basis->b.z + basis->c.z) < basis->b.modulus ) { + basis->b.x = basis->b.x + basis->c.x; basis->b.y = basis->b.y + basis->c.y; basis->b.y = basis->b.y + basis->c.y; + basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z); + } + + if ( modulus(basis->c.x - basis->a.x, basis->c.y - basis->a.z, basis->c.z - basis->a.z) < basis->c.modulus ) { + basis->c.x = basis->c.x - basis->a.x; basis->c.y = basis->c.y - basis->a.y; basis->c.y = basis->c.y - basis->a.y; + basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z); + } + if ( modulus(basis->c.x - basis->b.x, basis->c.y - basis->b.y, basis->c.z - basis->b.z) < basis->c.modulus ) { + basis->c.x = basis->c.x - basis->b.x; basis->c.y = basis->c.y - basis->b.y; basis->c.y = basis->c.y - basis->b.y; + basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z); + } + if ( modulus(basis->c.x + basis->a.x, basis->c.y + basis->a.z, basis->c.z + basis->a.z) < basis->c.modulus ) { + basis->c.x = basis->c.x + basis->a.x; basis->c.y = basis->c.y + basis->a.y; basis->c.y = basis->c.y + basis->a.y; + basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z); + } + if ( modulus(basis->c.x + basis->b.x, basis->c.y + basis->b.y, basis->c.z + basis->b.z) < basis->c.modulus ) { + basis->c.x = basis->c.x + basis->b.x; basis->c.y = basis->c.y + basis->b.y; basis->c.y = basis->c.y + basis->b.y; + basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z); + } + + } + +} + +/* Choose an initial set of three vectors for the expansion. + * This would probably be just as easy for the user to do... */ +static Basis *ipr_choose_initial_basis(ControlContext *ctx) { + + Basis *basis; + Reflection *reflection; + int cur_max; + + basis = malloc(sizeof(Basis)); + + /* Determine a vaguely sensible starting basis */ + basis->a.x = 0; basis->a.y = 0; basis->a.z = 0; basis->a.modulus = 1e30; + basis->b.x = 0; basis->b.y = 0; basis->b.z = 0; basis->b.modulus = 1e30; + basis->c.x = 0; basis->c.y = 0; basis->c.z = 0; basis->c.modulus = 1e30; + reflection = ctx->reflectionctx->reflections->next; /* First reflection, skip initial placeholder */ + cur_max = 1; + do { + double mod; + unsigned int angle_works; + + mod = modulus(reflection->x, reflection->y, reflection->z); + + angle_works = 1; + if ( basis->a.modulus != 1e30 ) { + double angle = angle_between(reflection->x, reflection->y, reflection->z, basis->a.x, basis->a.y, basis->a.z); + if ( angle < 20 ) angle_works = 0; + if ( angle > 160 ) angle_works = 0; + } + if ( basis->b.modulus != 1e30 ) { + double angle = angle_between(reflection->x, reflection->y, reflection->z, basis->b.x, basis->b.y, basis->b.z); + if ( angle < 20 ) angle_works = 0; + if ( angle > 160 ) angle_works = 0; + } + if ( basis->c.modulus != 1e30 ) { + double angle = angle_between(reflection->x, reflection->y, reflection->z, basis->c.x, basis->c.y, basis->c.z); + if ( angle < 20 ) angle_works = 0; + if ( angle > 160 ) angle_works = 0; + } + if ( (mod > 1e9) && angle_works ) { + cur_max = ipr_choose_max(basis); + if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 1) ) { + basis->a.x = reflection->x; + basis->a.y = reflection->y; + basis->a.z = reflection->z; + basis->a.modulus = mod; + } + if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 2) ) { + basis->b.x = reflection->x; + basis->b.y = reflection->y; + basis->b.z = reflection->z; + basis->b.modulus = mod; + } + if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 3) ) { + basis->c.x = reflection->x; + basis->c.y = reflection->y; + basis->c.z = reflection->z; + basis->c.modulus = mod; + } + } + reflection = reflection->next; + + } while ( reflection || (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) ); + + printf("IP: Initial basis:\n"); + printf("IP: a = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->a.x, basis->a.y, basis->a.z, basis->a.modulus); + printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus); + printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus); + + reflection_clear(ctx->reflectionctx); + reflection_add(ctx->reflectionctx, basis->a.x, basis->a.y, basis->a.z, 1, REFLECTION_GENERATED); + reflection_add(ctx->reflectionctx, basis->b.x, basis->b.y, basis->b.z, 1, REFLECTION_GENERATED); + reflection_add(ctx->reflectionctx, basis->c.x, basis->c.y, basis->c.z, 1, REFLECTION_GENERATED); + + ipr_try_compact(basis); + + return basis; + +} + +/* Start with the shortest vector + * Add it to a list, find out which of the three vectors, when added to the current length, + * produces the minimum length resultant. Repeat. Hence determine the order in which + * reflections should be measured. */ +static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) { + + ReflectionContext *ordered; + double max_res; + int max_order_a, max_order_b, max_order_c; + signed int h, k, l; + + /* NB This assumes the diffraction pattern is "vaguely" centered... */ + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) * ctx->pixel_size; + } else { + max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) / (ctx->lambda * ctx->camera_length); + } + max_res = max_res / 4; + + max_order_a = max_res/basis->a.modulus; + max_order_b = max_res/basis->b.modulus; + max_order_c = max_res/basis->c.modulus; + + ordered = reflection_init(); + + for ( h=-max_order_a; h<=max_order_a; h++ ) { + for ( k=-max_order_b; k<=max_order_b; k++ ) { + for ( l=-max_order_c; l<=max_order_c; l++ ) { + /* IPR-generated reflections are stored as hkl because their + g-vectors are going to be messed about with later. */ + reflection_add_index(ordered, h, k, l, 1, REFLECTION_GENERATED); + } + } + } + + return ordered; + +} + +/* Sort the list of reflections into order of increasing g */ +static void ipr_sort(ReflectionContext *rctx) { + + Reflection *prev; + Reflection *this; + Reflection *next; + + this = rctx->reflections; + next = this->next; + prev = NULL; + + while ( next ) { + + next = next->next; /* ! */ + prev = prev->next; + this = this->next; + + } + +} + +void ipr_reduce(ControlContext *ctx) { + + Basis *basis; + + basis = ipr_choose_initial_basis(ctx); + ctx->reflectionctx->basis = basis; + + /* Get rid of the original list and replace it with the prediction list */ + reflection_clear(ctx->reflectionctx); + free(ctx->reflectionctx); + ctx->reflectionctx = ipr_generate(ctx, basis); + ctx->reflectionctx->basis = basis; + + /* Sort the reflections into order of increasing g */ + ipr_sort(ctx->reflectionctx); + +} diff --git a/src/ipr.h b/src/ipr.h new file mode 100644 index 0000000..149e523 --- /dev/null +++ b/src/ipr.h @@ -0,0 +1,23 @@ +/* + * ipr.h + * + * Iterative prediction-refinement reconstruction + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef IPR_H +#define IPR_H + +#include "control.h" + +extern void ipr_reduce(ControlContext *ctx); + +#endif /* IPR_H */ diff --git a/src/itrans.c b/src/itrans.c new file mode 100644 index 0000000..0480f3f --- /dev/null +++ b/src/itrans.c @@ -0,0 +1,527 @@ +/* + * itrans.c + * + * Parameterise features in an image for reconstruction + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#include <stdlib.h> +#include <stdio.h> +#include <stdint.h> +#include <gsl/gsl_multifit_nlin.h> + +#include "reflections.h" +#include "control.h" +#include "imagedisplay.h" +#include "utils.h" + +#define MAX_LSQ_ITERATIONS 100 +#define PEAK_WINDOW_SIZE 20 + +typedef struct struct_sampledpoints { + gsl_matrix *coordinates; + gsl_vector *values; + unsigned int n_samples; +} SampledPoints; + +static unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) { + + int x, y; + unsigned int n_reflections; + + n_reflections = 0; + for ( x=1; x<ctx->width-1; x++ ) { + for ( y=1; y<ctx->height-1; y++ ) { + + double dx1, dx2, dy1, dy2; + double dxs, dys; + double grad; + + if ( ctx->max_d ) { + if ( (x-ctx->x_centre)*(x-ctx->x_centre) + (y-ctx->y_centre)*(y-ctx->y_centre) > ctx->max_d*ctx->max_d ) { + continue; + } + } + + /* Get gradients */ + dx1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*y]; + dx2 = image[(x-1)+ctx->width*y] - image[x+ctx->width*y]; + dy1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*(y+1)]; + dy2 = image[x+ctx->width*(y-1)] - image[x+ctx->width*y]; + + /* Average gradient measurements from both sides */ + dxs = ((dx1*dx1) + (dx2*dx2)) / 2; + dys = ((dy1*dy1) + (dy2*dy2)) / 2; + + /* Calculate overall gradient */ + grad = dxs + dys; + + if ( grad > 400 ) { + + unsigned int mask_x, mask_y; + unsigned int sx, sy; + double max; + unsigned int did_something = 1; + + mask_x = x; + mask_y = y; + + while ( (did_something) && (distance(mask_x, mask_y, x, y)<50) ) { + max = image[mask_x+ctx->width*mask_y]; + did_something = 0; + for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0); sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, ctx->height); sy++ ) { + for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0); sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, ctx->width); sx++ ) { + if ( image[sx+ctx->width*sy] > max ) { + max = image[sx+ctx->width*sy]; + mask_x = sx; + mask_y = sy; + did_something = 1; + } + } + } + } + + if ( !did_something ) { + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + reflection_add_from_reciprocal(ctx, (signed)(mask_x-ctx->width/2)*ctx->pixel_size, + (signed)(mask_y-ctx->height/2)*ctx->pixel_size, + tilt_degrees, image[mask_x + ctx->width*mask_y]); + } else { + reflection_add_from_dp(ctx, (signed)(mask_x-ctx->width/2), (signed)(mask_y-ctx->height/2), + tilt_degrees, image[mask_x + ctx->width*mask_y]); + } + if ( ctx->first_image ) { + imagedisplay_mark_point(imagedisplay, mask_x, mask_y); + } + n_reflections++; + } + + } + } + } + + return n_reflections; + +} + +static int itrans_peaksearch_lsq_f(const gsl_vector *gaussian, void *params, gsl_vector *residuals) { + + double cal; + double obs; + SampledPoints *samples = params; + double a, b, c, d, e, f; + unsigned int sample; + + a = gsl_vector_get(gaussian, 0); b = gsl_vector_get(gaussian, 1); c = gsl_vector_get(gaussian, 2); + d = gsl_vector_get(gaussian, 3); e = gsl_vector_get(gaussian, 4); f = gsl_vector_get(gaussian, 5); + + //printf("Trial parameters: a=%f b=%f c=%f d=%f e=%f f=%f\n", a, b, c, d, e, f); + + for ( sample=0; sample<samples->n_samples; sample++ ) { + int dx = gsl_matrix_get(samples->coordinates, sample, 0); + int dy = gsl_matrix_get(samples->coordinates, sample, 1); + cal = exp(a + b*dx + c*dy + d*dx*dx + e*dy*dy + f*dx*dy); + obs = gsl_vector_get(samples->values, sample); + gsl_vector_set(residuals, sample, (cal-obs)); + //printf("At %3i,%3i: residual=%f - %f = %f\n", dx, dy, cal, obs, cal-obs); + } + + return GSL_SUCCESS; + +} + +static int itrans_peaksearch_lsq_df(const gsl_vector *gaussian, void *params, gsl_matrix *J) { + + double a, b, c, d, e, f; + unsigned int sample; + SampledPoints *samples = params; + + a = gsl_vector_get(gaussian, 0); b = gsl_vector_get(gaussian, 1); c = gsl_vector_get(gaussian, 2); + d = gsl_vector_get(gaussian, 3); e = gsl_vector_get(gaussian, 4); f = gsl_vector_get(gaussian, 5); + + //printf("------a------------b------------c------------d------------e------------f-------\n"); + for ( sample=0; sample<samples->n_samples; sample++ ) { + + double ex; + int dx, dy; + double derivative; + + dx = gsl_matrix_get(samples->coordinates, sample, 0); + dy = gsl_matrix_get(samples->coordinates, sample, 1); + ex = exp(a + b*dx + c*dy + d*dx*dx + e*dy*dy + f*dx*dy); + + derivative = ex; /* wrt a */ + gsl_matrix_set(J, sample, 0, derivative); + derivative = dx * ex; /* wrt b */ + gsl_matrix_set(J, sample, 1, derivative); + derivative = dy * ex; /* wrt c */ + gsl_matrix_set(J, sample, 2, derivative); + derivative = dx * dx * ex; /* wrt d */ + gsl_matrix_set(J, sample, 3, derivative); + derivative = dy * dy * ex; /* wrt e */ + gsl_matrix_set(J, sample, 4, derivative); + derivative = dx * dy * ex; /* wrt f */ + gsl_matrix_set(J, sample, 5, derivative); + + /*if ( (dx == 0) && (dy == 0) ) { + printf("> "); + } else { + printf("| "); + } + printf("%10.2e | %10.2e | %10.2e | %10.2e | %10.2e | %10.2e", gsl_matrix_get(J, sample, 0), gsl_matrix_get(J, sample, 1), + gsl_matrix_get(J, sample, 2), gsl_matrix_get(J, sample, 3), gsl_matrix_get(J, sample, 4), gsl_matrix_get(J, sample, 5)); + if ( (dx == 0) && (dy == 0) ) { + printf(" <\n"); + } else { + printf(" |\n"); + }*/ + + } + //printf("-------------------------------------------------------------------------------\n"); + + return GSL_SUCCESS; + +} + +static int itrans_peaksearch_lsq_fdf(const gsl_vector *gaussian, void *params, gsl_vector *f, gsl_matrix *J) { + + itrans_peaksearch_lsq_f(gaussian, params, f); + itrans_peaksearch_lsq_df(gaussian, params, J); + + return GSL_SUCCESS; + +} + +static void itrans_interpolate(int16_t *image, int width, int x, int y) { + + int a, b, c, d; + double av_horiz, av_vert, av; + printf("Interpolating...\n"); + a = image[(x-1) + width*y]; b = image[(x+1) + width*y]; + c = image[x + width*(y-1)]; d = image[x + width*(y+1)]; + av_horiz = (a+b)/2; av_vert = (c+d)/2; + av = (av_horiz+av_vert) / 2; + image[x + width*y] = av; + +} + +static unsigned int itrans_peaksearch_threshold(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) { + + int width, height; + int x, y; + unsigned int n_reflections = 0; + + width = ctx->width; height = ctx->height; + + /* Simple Thresholding */ + for ( y=0; y<height; y++ ) { + for ( x=0; x<width; x++ ) { + + if ( image[x + width*y] > 10 ) { + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + reflection_add_from_reciprocal(ctx, (signed)(x-ctx->x_centre)*ctx->pixel_size, (signed)(y-ctx->y_centre)*ctx->pixel_size, + tilt_degrees, image[x + width*y]); + } else { + reflection_add_from_dp(ctx, (signed)(x-ctx->x_centre), (signed)(y-ctx->y_centre), tilt_degrees, image[x + width*y]); + } + if ( ctx->first_image ) { + imagedisplay_mark_point(imagedisplay, x, y); + } + n_reflections++; + } + + } + } + + return n_reflections; + +} + +static unsigned int itrans_peaksearch_adaptive_threshold(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) { + + int16_t max_val = 0; + int width, height; + unsigned int n_reflections = 0; + + width = ctx->width; height = ctx->height; + + /* Adaptive Thresholding */ + do { + + int max_x = 0; + int max_y = 0;; + int x, y; + + /* Locate the highest point */ + max_val = 0; + for ( y=0; y<height; y++ ) { + for ( x=0; x<width; x++ ) { + + if ( image[x + width*y] > max_val ) { + max_val = image[x + width*y]; + max_x = x; max_y = y; + } + + } + } + + if ( max_val > 50 ) { + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + reflection_add_from_reciprocal(ctx, (signed)(max_x-ctx->x_centre)*ctx->pixel_size, (signed)(max_y-ctx->y_centre)*ctx->pixel_size, + tilt_degrees, max_val); + } else { + reflection_add_from_dp(ctx, (signed)(max_x-ctx->x_centre), (signed)(max_y-ctx->y_centre), tilt_degrees, max_val); + } + if ( ctx->first_image ) { + imagedisplay_mark_point(imagedisplay, max_x, max_y); + } + n_reflections++; + + /* Remove it and its surroundings */ + for ( y=((max_y-10>0)?max_y-10:0); y<((max_y+10)<height?max_y+10:height); y++ ) { + for ( x=((max_x-10>0)?max_x-10:0); x<((max_x+10)<width?max_x+10:width); x++ ) { + image[x + width*y] = 0; + } + } + } + + } while ( max_val > 50 ); + + return n_reflections; + +} + +static unsigned int itrans_peaksearch_lsq(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) { + + int16_t max_val = 0; + int width, height; + unsigned int n_reflections = 0; + + width = ctx->width; height = ctx->height; + + /* Least-Squares Craziness. NB Doesn't quite work... */ + do { + + int max_x = 0; + int max_y = 0;; + int x, y; + gsl_multifit_fdfsolver *s; + gsl_multifit_function_fdf f; + unsigned int iter; + gsl_vector *gaussian; + int iter_status, conv_status; + gsl_matrix *covar; + SampledPoints samples; + int sample; + double ga, gb, gc, gd, ge, gf; + gboolean sanity; + int dx, dy; + int bb_lh, bb_rh, bb_top, bb_bot; + int16_t ival; + double sx, sy; + + /* Locate the highest point */ + max_val = 0; + for ( y=10; y<height-10; y++ ) { + for ( x=10; x<width-10; x++ ) { + + if ( image[x + width*y] > max_val ) { + max_val = image[x + width*y]; + max_x = x; max_y = y; + } + + } + } + x = max_x; y = max_y; + + /* Try fitting a Gaussian to this region and see what happens... */ + f.p = 6; + + /* Set initial estimate of the fit parameters. I = exp(a + bx + cy + dx^2 + ey^2 + fxy) */ + gaussian = gsl_vector_alloc(f.p); + gsl_vector_set(gaussian, 0, log(max_val)); /* a */ + gsl_vector_set(gaussian, 1, -0.01); /* b */ + gsl_vector_set(gaussian, 2, -0.01); /* c */ + gsl_vector_set(gaussian, 3, -0.01); /* d */ + gsl_vector_set(gaussian, 4, -0.01); /* e */ + gsl_vector_set(gaussian, 5, -0.01); /* f */ + + /* Determine the bounding box which contains most of the peak */ + /* Left */ dx = 0; do { ival = image[(x+dx) + width*(y+0)]; dx--; } while ( (ival>max_val/300) && (x+dx>0) ); bb_lh = dx; + /* Right */ dx = 0; do { ival = image[(x+dx) + width*(y+0)]; dx++; } while ( (ival>max_val/300) && (x+dx<width) ); bb_rh = dx; + /* Top */ dy = 0; do { ival = image[(x+0) + width*(y+dy)]; dy++; } while ( (ival>max_val/300) && (y+dy<height) ); bb_top = dy; + /* Bottom */ dy = 0; do { ival = image[(x+0) + width*(y+dy)]; dy--; } while ( (ival>max_val/300) && (y+dy>0) ); bb_bot = dy; + printf("Peak %i,%i: bounding box %i<dx<%i, %i<dy<%i\n", x, y, bb_lh, bb_rh, bb_bot, bb_top); + sanity = TRUE; + if ( (bb_rh-bb_lh) > 300 ) sanity = FALSE; + if ( (bb_top-bb_bot) > 300 ) sanity = FALSE; + if ( sanity ) { + + /* Choose which points inside this bounding box to use for curve fitting */ + samples.n_samples = ((bb_top-bb_bot)+1)*((bb_rh-bb_lh)+1); + samples.coordinates = gsl_matrix_alloc(samples.n_samples, 2); + sample = 0; + for ( sx=bb_lh; sx<=bb_rh; sx++ ) { + for ( sy=bb_bot; sy<=bb_top; sy++ ) { + gsl_matrix_set(samples.coordinates, sample, 0, sx); gsl_matrix_set(samples.coordinates, sample, 1, sy); + sample++; + } + } + + samples.values = gsl_vector_alloc(samples.n_samples); + f.n = samples.n_samples; + for ( sample=0; sample<samples.n_samples; sample++ ) { + int dx = gsl_matrix_get(samples.coordinates, sample, 0); + int dy = gsl_matrix_get(samples.coordinates, sample, 1); + gsl_vector_set(samples.values, sample, image[(x+dx) + width*(y+dy)]); + //printf("At %3i,%3i: value=%i\n", dx, dy, image[(x+dx) + width*(y+dy)]); + } + f.params = &samples; + + /* Execute the LSQ fitting procedure */ + s = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, f.n, f.p); + f.f = &itrans_peaksearch_lsq_f; f.df = &itrans_peaksearch_lsq_df; f.fdf = &itrans_peaksearch_lsq_fdf; + gsl_multifit_fdfsolver_set(s, &f, gaussian); + iter = 0; conv_status = GSL_CONTINUE; + do { + + /* Iterate */ + iter_status = gsl_multifit_fdfsolver_iterate(s); + iter++; + + /* Check for error */ + if ( iter_status ) { + break; + } + + /* Test for convergence */ + conv_status = gsl_multifit_test_delta(s->dx, s->x, 1e-6, 1e-6); + + } while ( (conv_status == GSL_CONTINUE) && (iter < MAX_LSQ_ITERATIONS) ); + + /* See how good the fit is */ + covar = gsl_matrix_alloc(f.p, f.p); + gsl_multifit_covar(s->J, 0.0, covar); + ga = gsl_vector_get(s->x, 0); gb = gsl_vector_get(s->x, 1); gc = gsl_vector_get(s->x, 2); + gd = gsl_vector_get(s->x, 3); ge = gsl_vector_get(s->x, 4); gf = gsl_vector_get(s->x, 5); + if ( fabs(exp(ga + gb*100 + gc*100 + gd*100*100 + ge*100*100 + gf*100*100)) > 10 ) { + printf("Failed sanity check: %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n", x, y, ga, gb, gc, gd, ge, gf); + sanity = FALSE; + } else { + sanity = TRUE; + } + if ( (conv_status == GSL_SUCCESS) && (!iter_status) && (sanity) ) { + + /* Good fit! */ + int dx, dy; + int bb_top, bb_bot, bb_lh, bb_rh; + double frac; + + printf("Fit converged after %i iterations: Centre %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n", iter, x, y, ga, gb, gc, gd, ge, gf); + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + reflection_add_from_reciprocal(ctx, (x-ctx->x_centre)*ctx->pixel_size, (y-ctx->y_centre)*ctx->pixel_size, + tilt_degrees, image[x + width*y]); + } else { + reflection_add_from_dp(ctx, (x-ctx->x_centre), (y-ctx->y_centre), tilt_degrees, image[x + width*y]); + } + if ( ctx->first_image ) { + imagedisplay_mark_point(imagedisplay, x, y); + } + n_reflections++; + + /* Remove this peak from the image */ + /* Find right-hand edge of bounding box */ + dx = 0; do { + double pval, ival; + pval = exp(ga + gb*dx + gc*0 + gd*dx*dx + ge*0*0 + gf*dx*-0); + ival = image[(x+dx) + width*(y+0)]; + frac = pval / ival; + dx++; + } while ( (frac > 0.1) && (x+dx<width) ); + bb_rh = dx; + /* Find left-hand edge of bounding box */ + dx = 0; do { + double pval, ival; + pval = exp(ga + gb*dx + gc*0 + gd*dx*dx + ge*0*0 + gf*dx*-0); + ival = image[(x+dx) + width*(y+0)]; + frac = pval / ival; + dx--; + } while ( (frac > 0.1) && (x+dx>0) ); + bb_lh = dx; + /* Find top edge of bounding box */ + dy = 0; do { + double pval, ival; + pval = exp(ga + gb*0 + gc*dy + gd*0*0 + ge*dy*dy + gf*0*dy); + ival = image[(x+0) + width*(y+dy)]; + frac = pval / ival; + dy++; + } while ( (frac > 0.1) && (y+dy<height) ); + bb_top = dy; + /* Find bottom edge of bounding box */ + dy = 0; do { + double pval, ival; + pval = exp(ga + gb*0 + gc*dy + gd*0*0 + ge*dy*dy + gf*0*dy); + ival = image[(x+0) + width*(y+dy)]; + frac = pval / ival; + dy--; + } while ( (frac > 0.1) && (y+dy>0) ); + bb_bot = dy; + printf("Fitted peak bounding box %i<dx<%i, %i<dy<%i\n", bb_lh, bb_rh, bb_bot, bb_top); + for ( dx=bb_lh; dx<bb_rh; dx++ ) { + for ( dy=bb_bot; dy<bb_top; dy++ ) { + double pval; + pval = exp(ga + gb*dx + gc*dy + gd*dx*dx + ge*dy*dy + gf*dx*dy); + image[(x+dx) + width*(y+dy)] -= pval; + } + } + + } else { + itrans_interpolate(image, width, x, y); + } + + gsl_matrix_free(covar); + gsl_multifit_fdfsolver_free(s); + gsl_vector_free(gaussian); + gsl_matrix_free(samples.coordinates); + gsl_vector_free(samples.values); + + } else { + printf("Failed bounding box sanity check\n"); + itrans_interpolate(image, width, x, y); + } + + + } while ( max_val > 50 ); + + return n_reflections; +} + +void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees) { + + unsigned int n_reflections; + ImageDisplay *imagedisplay = NULL; + + if ( ctx->first_image ) { + imagedisplay = imagedisplay_open(image, ctx->width, ctx->height, "Image Display"); + imagedisplay_add_tilt_axis(imagedisplay, ctx, ctx->omega); + } + + switch( ctx->psmode ) { + case PEAKSEARCH_THRESHOLD : n_reflections = itrans_peaksearch_threshold(image, ctx, tilt_degrees, imagedisplay); break; + case PEAKSEARCH_ADAPTIVE_THRESHOLD : n_reflections = itrans_peaksearch_adaptive_threshold(image, ctx, tilt_degrees, imagedisplay); break; + case PEAKSEARCH_LSQ : n_reflections = itrans_peaksearch_lsq(image, ctx, tilt_degrees, imagedisplay); break; + case PEAKSEARCH_ZAEFFERER : n_reflections = itrans_peaksearch_zaefferer(image, ctx, tilt_degrees, imagedisplay); break; + default: n_reflections = 0; + } + + if ( ctx->rmode == RECONSTRUCTION_PREDICTION ) { + + } + + ctx->first_image = 0; + printf(" ..... %i peaks found\n", n_reflections); + +} diff --git a/src/itrans.h b/src/itrans.h new file mode 100644 index 0000000..37f128b --- /dev/null +++ b/src/itrans.h @@ -0,0 +1,24 @@ +/* + * itrans.h + * + * Parameterise features in an image for reconstruction + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef ITRANS_H +#define ITRANS_H + +#include <stdint.h> + +#include "control.h" + +extern void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees); + +#endif /* ITRANS_H */ diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..272c9a8 --- /dev/null +++ b/src/main.c @@ -0,0 +1,183 @@ +/* + * main.c + * + * The Top Level Source File + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#define _GNU_SOURCE +#include <string.h> +#include <stdarg.h> +#include <stdlib.h> +#include <gtk/gtk.h> +#include <math.h> +#include <gdk/gdkgl.h> +#include <gtk/gtkgl.h> +#include <sys/stat.h> + +#include "displaywindow.h" +#include "reflections.h" +#include "mrc.h" +#include "qdrp.h" +#include "ipr.h" + +static gint main_method_window_response(GtkWidget *method_window, gint response, ControlContext *ctx) { + + if ( response == GTK_RESPONSE_OK ) { + + int val = -1; + + switch ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_algorithm)) ) { + case 0 : ctx->rmode = RECONSTRUCTION_MAPPING; break; + case 1 : ctx->rmode = RECONSTRUCTION_PREDICTION; break; + default: abort(); + } + + switch ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_peaksearch)) ) { + case 0 : ctx->psmode = PEAKSEARCH_THRESHOLD; break; + case 1 : ctx->psmode = PEAKSEARCH_ADAPTIVE_THRESHOLD; break; + case 2 : ctx->psmode = PEAKSEARCH_LSQ; break; + case 3 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break; + default: abort(); + } + + gtk_widget_destroy(method_window); + while ( gtk_events_pending() ) gtk_main_iteration(); + + if ( ctx->inputfiletype == INPUT_QDRP ) { + val = qdrp_read(ctx); + } else if ( ctx->inputfiletype == INPUT_MRC ) { + val = mrc_read(ctx); + } + + if ( ctx->rmode == RECONSTRUCTION_PREDICTION ) { + ipr_reduce(ctx); + } + + if ( val == 0 ) { + displaywindow_open(ctx); + } else { + fprintf(stderr, "Reconstruction failed.\n"); + gtk_exit(0); + } + + } else { + gtk_exit(0); + } + + return 0; + +} + +void main_method_dialog_open(ControlContext *ctx) { + + GtkWidget *method_window; + GtkWidget *vbox; + GtkWidget *hbox; + GtkWidget *table; + GtkWidget *method_label; + GtkWidget *peaksearch_label; + + method_window = gtk_dialog_new_with_buttons("Reconstruction Parameters", NULL, + GTK_DIALOG_DESTROY_WITH_PARENT, GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE, GTK_STOCK_OK, GTK_RESPONSE_OK, NULL); + + vbox = gtk_vbox_new(FALSE, 0); + hbox = gtk_hbox_new(TRUE, 0); + gtk_box_pack_start(GTK_BOX(GTK_DIALOG(method_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); + + table = gtk_table_new(3, 2, FALSE); + + method_label = gtk_label_new("Reconstruction Algorithm: "); + gtk_table_attach_defaults(GTK_TABLE(table), method_label, 1, 2, 1, 2); + gtk_misc_set_alignment(GTK_MISC(method_label), 1, 0.5); + ctx->combo_algorithm = gtk_combo_box_new_text(); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_algorithm), "Feature Detection and Mapping"); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_algorithm), "Iterative Prediction and Refinement"); + gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_algorithm), 1); + gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_algorithm, 2, 3, 1, 2); + gtk_table_set_row_spacing(GTK_TABLE(table), 1, 5); + + peaksearch_label = gtk_label_new("Peak Search: "); + gtk_table_attach_defaults(GTK_TABLE(table), peaksearch_label, 1, 2, 2, 3); + gtk_misc_set_alignment(GTK_MISC(peaksearch_label), 1, 0.5); + ctx->combo_peaksearch = gtk_combo_box_new_text(); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Simple Thresholding"); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Adaptive Thresholding"); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Least-Squares Fit"); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Zaefferer Gradient Search"); + gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 3); + gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 2, 3); + + gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), TRUE, TRUE, 5); + + g_signal_connect(G_OBJECT(method_window), "response", G_CALLBACK(main_method_window_response), ctx); + + gtk_widget_show_all(method_window); + +} + + +int main(int argc, char *argv[]) { + + char *filename; + ControlContext *ctx; + InputFileType type; + struct stat stat_buffer; + + gtk_init(&argc, &argv); + + if ( gtk_gl_init_check(&argc, &argv) == FALSE ) { + fprintf(stderr, "Could not initialise gtkglext\n"); + return 1; + } + + if ( gdk_gl_query_extension() == FALSE ) { + fprintf(stderr, "OpenGL not supported\n"); + return 1; + } + + if ( argc != 2 ) { + fprintf(stderr, "Syntax: %s [<MRC file> | qdrp.rc]\n", argv[0]); + return 1; + } + + filename = basename(argv[1]); + ctx = malloc(sizeof(ControlContext)); + type = INPUT_NONE; + + if ( strcmp(filename, "qdrp.rc") == 0 ) { + printf("QDRP input file detected.\n"); + ctx->inputfiletype = INPUT_QDRP; + } else if ( strcmp(filename+(strlen(filename)-4), ".mrc") == 0 ) { + printf("MRC tomography file detected.\n"); + ctx->inputfiletype = INPUT_MRC; + } else { + fprintf(stderr, "Input file must either be an MRC tomography file or a QDRP-style control file.\n"); + return 1; + } + ctx->filename = strdup(argv[1]); + ctx->max_d = 0; + + if ( stat(argv[1], &stat_buffer) == -1 ) { + fprintf(stderr, "File '%s' not found\n", argv[1]); + return 1; + } + + main_method_dialog_open(ctx); + + gtk_main(); + + free(ctx->filename); + free(ctx); + + return 0; + +} diff --git a/src/main.h b/src/main.h new file mode 100644 index 0000000..1cf20d1 --- /dev/null +++ b/src/main.h @@ -0,0 +1,19 @@ +/* + * main.h + * + * The Top Level Source File + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef MAIN_H +#define MAIN_H + +#endif /* MAIN_H */ diff --git a/src/mrc.c b/src/mrc.c new file mode 100644 index 0000000..4815626 --- /dev/null +++ b/src/mrc.c @@ -0,0 +1,141 @@ +/* + * mrc.c + * + * Read the MRC tomography format + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <stdio.h> +#include <stdint.h> + +#include "mrc.h" +#include "control.h" +#include "imagedisplay.h" +#include "itrans.h" +#include "reflections.h" + +int mrc_read(ControlContext *ctx) { + + FILE *fh; + MRCHeader mrc; + MRCExtHeader ext[1024]; + unsigned int i; + unsigned int extsize; + double sx, sy, sz; + unsigned int nx, ny, nz; + double pixel_size; + int16_t *image_total; + unsigned int x, y; + int max_x, max_y; + int16_t max_val; + ImageDisplay *sum_id; + + fh = fopen(ctx->filename, "rb"); + + /* Read primary header */ + fread(&mrc, sizeof(MRCHeader), 1, fh); + printf("%i images in series\n", mrc.nz); + if ( mrc.mode != 1 ) { + fprintf(stderr, "MR: Unknown MRC image mode\n"); + fclose(fh); + return -1; + } + if ( (mrc.nxstart != 0) || (mrc.nystart != 0) ) { + fclose(fh); + fprintf(stderr, "MR: Image origin must be at zero: found at %i,%i\n", mrc.nxstart, mrc.nystart); + return -1; + } + + /* Read all extended headers, one by one */ + extsize = 4*mrc.numintegers + 4*mrc.numfloats; + if ( extsize > sizeof(MRCExtHeader) ) { + fclose(fh); + fprintf(stderr, "MR: MRC extended header is too big - tweak mrc.h\n"); + return -1; + } + for ( i=0; i<mrc.nz; i++ ) { + fread(&ext[i], extsize, 1, fh); + } + + pixel_size = ext[0].pixel_size; + nx = mrc.nx; + ny = mrc.ny; + nz = (mrc.nx > mrc.ny)?mrc.nx:mrc.ny; + sx = nx * pixel_size; + sy = ny * pixel_size; + sz = nz * pixel_size; + ctx->reflectionctx = reflection_init(); + ctx->fmode = FORMULATION_PIXELSIZE; + ctx->first_image = 1; + ctx->width = mrc.nx; + ctx->height = mrc.ny; + + printf("Judging centre...\n"); + image_total = malloc(mrc.ny * mrc.nx * sizeof(int16_t)); + for ( y=0; y<mrc.ny; y++ ) { + for ( x=0; x<mrc.nx; x++ ) { + image_total[x + mrc.nx*y] = 0; + } + } + for ( i=0; i<mrc.nz; i++ ) { + int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t)); + fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET); + fread(image, mrc.nx*mrc.ny*2, 1, fh); + for ( y=0; y<mrc.ny; y++ ) { + for ( x=0; x<mrc.nx; x++ ) { + image_total[x + mrc.nx*y] += image[x + mrc.nx*y]/mrc.nz; + } + } + } + sum_id = imagedisplay_open(image_total, mrc.nx, mrc.ny, "Sum of All Images"); + /* Locate the highest point */ + max_val = 0; max_x = 0; max_y = 0; + for ( y=0; y<ctx->height; y++ ) { + for ( x=0; x<ctx->width; x++ ) { + if ( image_total[x + ctx->width*y] > max_val ) { + max_val = image_total[x + ctx->width*y]; + max_x = x; max_y = y; + } + } + } + imagedisplay_mark_point(sum_id, max_x, max_y); + ctx->x_centre = max_x; + ctx->y_centre = max_y; + + for ( i=0; i<mrc.nz; i++ ) { + + int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t)); + + printf("Image #%3i: tilt=%f omega=%f L=%f\t", i, ext[i].a_tilt, ext[i].tilt_axis, ext[i].magnification); + ctx->camera_length = ext[i].magnification; + ctx->lambda = 2.51e-12; /* 200kV. Fudged until Max puts the HT voltage in the MRC headers */ + ctx->omega = ext[i].tilt_axis; + ctx->pixel_size = ext[i].pixel_size; + + fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET); + fread(image, mrc.nx*mrc.ny*2, 1, fh); + + for ( y=0; y<mrc.ny; y++ ) { + for ( x=0; x<mrc.nx; x++ ) { + image_total[x + mrc.nx*y] += image[x + mrc.nx*y]/mrc.nz; + } + } + + itrans_process_image(image, ctx, ext[i].a_tilt); + free(image); + + } + + fclose(fh); + + return 0; + +} diff --git a/src/mrc.h b/src/mrc.h new file mode 100644 index 0000000..8633f05 --- /dev/null +++ b/src/mrc.h @@ -0,0 +1,121 @@ +/* + * mrc.h + * + * Read the MRC tomography format + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef MRC_H +#define MRC_H + +#include <stdint.h> + +#include "control.h" + +typedef struct struct_mrcheader { + + int32_t nx; + int32_t ny; + int32_t nz; + int32_t mode; + int32_t nxstart; + int32_t nystart; + int32_t nzstart; + int32_t mx; + int32_t my; + int32_t mz; + float xlen; + float ylen; + float zlen; + float alpha; + float beta; + float gamma; + int32_t mapc; + int32_t mapr; + int32_t maps; + float amin; + float amax; + float amean; + int16_t ispg; /* Space group number */ + int16_t nsymbt; + int32_t next; + int16_t dvid; + char extra[30]; + int16_t numintegers; + int16_t numfloats; + int16_t sub; + int16_t zfac; + float min2; + float max2; + float min3; + float max3; + float min4; + float max4; + int16_t idtype; + int16_t lens; + int16_t nd1; + int16_t nd2; + int16_t vd1; + int16_t vd2; + float tiltangles[9]; + float zorg; + float xorg; + float yorg; + int32_t nlabl; + char data[10][80]; + +} MRCHeader; + +typedef struct struct_mrcextheader { + + float a_tilt; + float b_tilt; + float x_stage; + float y_stage; + float z_stage; + float x_shift; + float y_shift; + float defocus; + float exp_time; + float mean_int; + float tilt_axis; + float pixel_size; + float magnification; + float mic_type; + float gun_type; + float d_number; + float voltage; + float focus_spread; + float mtf; + float df_start; + float focus_step; + float dac_setting; + float cs; + float semi_conv; + float info_limit; + float num_images; + float num_in_series; + float coma1; + float coma2; + float astig21; + float astig22; + float astig31; + float astig32; + float cam_type; + float cam_pos; + float padding[64]; /* Need to guarantee that reading an extended header (of any size) + will never write beyond the boundary of this structure... */ + +} MRCExtHeader; + +extern int mrc_read(ControlContext *ctx); + +#endif /* MRC_H */ diff --git a/src/qdrp.c b/src/qdrp.c new file mode 100644 index 0000000..42be13f --- /dev/null +++ b/src/qdrp.c @@ -0,0 +1,254 @@ +/* + * qdrp.c + * + * Handle QDRP-style control files + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#define _GNU_SOURCE +#include <stdlib.h> +#include <string.h> +#include <stdio.h> + +#include "readpng.h" +#include "control.h" +#include "reflections.h" + +static void qdrp_chomp(char *str) { + + size_t i; + + for ( i=0; i<strlen(str); i++ ) { + if ( str[i] == '\r' ) { + str[i] = '\0'; + return; + } + if ( str[i] == '\n' ) { + str[i] = '\0'; + return; + } + } +} + +#define skip_whitespace \ + while ( (( line[pos] == ' ' ) || ( line[pos] == '\t' )) && ( pos < strlen(line) ) ) { \ + pos++; \ + } + +#define skip_chars \ + while ( ( line[pos] != ' ' ) && ( line[pos] != '\t' ) && ( pos < strlen(line) ) ) { \ + pos++; \ + } + +static int qdrp_parsefileline(ControlContext *ctx, const char *line) { + + size_t pos = 0; + size_t mark = 0; + char *tilt_s; + char *file; + double tilt; + + skip_whitespace; + mark = pos; + skip_chars; + tilt_s = strndup(line+mark, pos-mark); + + skip_whitespace; + mark = pos; + skip_chars; + file = strndup(line+mark, pos-mark); + + tilt = strtod(tilt_s, NULL); + free(tilt_s); + printf("QD: Reading file: Tilt = %f deg, File='%s'\n", tilt, file); + + if ( readpng_read(file, tilt, ctx) ) { + printf("Reconstruction failed.\n"); + return 1; + } + + return 0; + +} + +static int qdrp_parseline(ControlContext *ctx, const char *line) { + + if ( ctx->started ) { + return qdrp_parsefileline(ctx, line); + } + + if ( ( line[0] == '-' ) && ( line[1] == '-') ) { + + if ( !ctx->camera_length_set ) { + printf("QD: Parameter 'camera-length' not specified!\n"); + return 1; + } + + if ( !ctx->lambda_set ) { + printf("QD: Parameter 'lambda' not specified!\n"); + return 1; + } + + if ( !ctx->resolution_set ) { + printf("QD: Parameter 'resolution' not specified!\n"); + return 1; + } + + if ( !ctx->omega_set ) { + printf("QD: Parameter 'omega' not specified, and not using tilt-axis-search.\n"); + return 1; + } + + ctx->started = 1; + ctx->reflectionctx = reflection_init(); + + } + + if ( !ctx->started ) { + + size_t pos = 0; + size_t mark = 0; + + skip_whitespace; + + if ( strncasecmp(line+pos, "resolution", 10) == 0 ) { + + char *resolution_s; + + skip_chars; + skip_whitespace; + mark = pos; + skip_chars; + resolution_s = strndup(line+mark, pos-mark); + ctx->resolution = strtod(resolution_s, NULL); + free(resolution_s); + ctx->resolution_set = 1; + printf("QD: resolution = %f pixels/m\n", ctx->resolution); + + } + + if ( strncasecmp(line+pos, "lambda", 6) == 0 ) { + + char *lambda_s; + + skip_chars; + skip_whitespace; + mark = pos; + skip_chars; + lambda_s = strndup(line+mark, pos-mark); + ctx->lambda = strtod(lambda_s, NULL); + free(lambda_s); + ctx->lambda_set = 1; + printf("QD: lambda = %e m\n", ctx->lambda); + + } + + if ( strncasecmp(line+pos, "camera-length", 13) == 0 ) { + + char *camera_length_s; + + skip_chars; + skip_whitespace; + mark = pos; + skip_chars; + camera_length_s = strndup(line+mark, pos-mark); + ctx->camera_length = strtod(camera_length_s, NULL); + free(camera_length_s); + ctx->camera_length_set = 1; + printf("QD: camera-length = %f m\n", ctx->camera_length); + + } + + if ( strncasecmp(line+pos, "max_d", 5) == 0 ) { + + char *max_d_s; + + skip_chars; + skip_whitespace; + mark = pos; + skip_chars; + max_d_s = strndup(line+mark, pos-mark); + ctx->max_d = strtod(max_d_s, NULL); + free(max_d_s); + ctx->max_d_set = 1; + if ( ctx->max_d == 0 ) { + printf("QD: max_d=0 is Not Allowed!\n"); + return 1; + } + printf("QD: max_d = %i pixels\n", ctx->max_d); + + } + + if ( strncasecmp(line+pos, "omega", 5) == 0 ) { + + char *omega_s; + + skip_chars; + skip_whitespace; + mark = pos; + skip_chars; + omega_s = strndup(line+mark, pos-mark); + ctx->omega = strtod(omega_s, NULL); + free(omega_s); + ctx->omega_set = 1; + printf("CO: omega = %f deg\n", ctx->omega); + + } + + } + + return 0; + +} + +int qdrp_read(ControlContext *ctx) { + + char *line; + + line = malloc(256); + + ctx->camera_length_set = 0; + ctx->omega_set = 0; + ctx->max_d_set = 0; + ctx->resolution_set = 0; + ctx->lambda_set = 0; + ctx->max_d = 0; + ctx->started = 0; + ctx->first_image = 1; + + FILE *fh; + + fh = fopen(ctx->filename, "r"); + if ( !fh ) { + printf("QD: Couldn't open control file '%s'\n", ctx->filename); + return -1; + } + + while ( !feof(fh) ) { + fgets(line, 256, fh); + if ( !feof(fh) ) { + qdrp_chomp(line); + if ( line[0] != '#' ) { + if ( qdrp_parseline(ctx, line) ) { + fclose(fh); + free(line); + return -1; + } + } + } + } + + fclose(fh); + + free(line); + + return 0; /* Success */ + +} diff --git a/src/qdrp.h b/src/qdrp.h new file mode 100644 index 0000000..3bbfaa0 --- /dev/null +++ b/src/qdrp.h @@ -0,0 +1,22 @@ +/* + * qdrp.h + * + * Handle QDRP-style control files + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef QDRP_H +#define QDRP_H + +#include "control.h" + +extern int qdrp_read(ControlContext *ctx); + +#endif /* QDRP_H */ diff --git a/src/readpng.c b/src/readpng.c new file mode 100644 index 0000000..bb5a60b --- /dev/null +++ b/src/readpng.c @@ -0,0 +1,146 @@ +/* + * readpng.c + * + * Read images + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <string.h> +#include <stdio.h> +#include <png.h> +#include <math.h> + +#include "control.h" +#include "reflections.h" +#include "itrans.h" + +int readpng_read(const char *filename, double tilt, ControlContext *ctx) { + + FILE *fh; + png_bytep header; + png_structp png_ptr; + png_infop info_ptr; + png_infop end_info; + unsigned int width; + unsigned int height; + unsigned int bit_depth; + unsigned int channels; + png_bytep *row_pointers; + unsigned int x; + unsigned int y; + unsigned int xorig; + unsigned int yorig; + int16_t *image; + + /* Open file */ + fh = fopen(filename, "rb"); + if ( !fh ) { + printf("RI: Couldn't open file '%s'\n", filename); + return -1; + } + + /* Check it's actually a PNG file */ + header = malloc(8); + fread(header, 1, 8, fh); + if ( png_sig_cmp(header, 0, 8)) { + printf("RI: Not a PNG file.\n"); + free(header); + fclose(fh); + return -1; + } + free(header); + + png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if ( !png_ptr ) { + printf("RI: Couldn't create PNG read structure.\n"); + fclose(fh); + return -1; + } + + info_ptr = png_create_info_struct(png_ptr); + if ( !info_ptr ) { + png_destroy_read_struct(&png_ptr, (png_infopp)NULL, (png_infopp)NULL); + printf("RI: Couldn't create PNG info structure.\n"); + fclose(fh); + return -1; + } + + end_info = png_create_info_struct(png_ptr); + if ( !end_info ) { + png_destroy_read_struct(&png_ptr, &info_ptr, (png_infopp)NULL); + printf("RI: Couldn't create PNG end info structure.\n"); + fclose(fh); + return -1; + } + + if ( setjmp(png_jmpbuf(png_ptr)) ) { + png_destroy_read_struct(&png_ptr, &info_ptr, &end_info); + fclose(fh); + printf("RI: PNG read failed.\n"); + return -1; + } + + png_init_io(png_ptr, fh); + png_set_sig_bytes(png_ptr, 8); + + /* Read! */ + png_read_png(png_ptr, info_ptr, PNG_TRANSFORM_STRIP_ALPHA, NULL); + + width = png_get_image_width(png_ptr, info_ptr); + height = png_get_image_height(png_ptr, info_ptr); + bit_depth = png_get_bit_depth(png_ptr, info_ptr); + channels = png_get_channels(png_ptr, info_ptr); +// printf("RI: width=%i, height=%i, depth=%i, channels=%i\n", width, height, bit_depth, channels); + if ( (bit_depth != 16) && (bit_depth != 8) ) { + printf("RI: Whoops! Can't handle images with other than 8 or 16 bpp yet...\n"); + png_destroy_read_struct(&png_ptr, &info_ptr, &end_info); + fclose(fh); + return -1; + } + + /* Get image data */ + row_pointers = png_get_rows(png_ptr, info_ptr); + + xorig = 320; + yorig = 320; /* For now */ + image = malloc(height * width * sizeof(int16_t)); + ctx->width = width; ctx->height = height; + ctx->fmode = FORMULATION_CLEN; + + for ( y=0; y<height; y++ ) { + for ( x=0; x<width; x++ ) { + + int val = 0; + + if ( (x-xorig)*(x-xorig) + (y-yorig)*(y-yorig) < ctx->max_d*ctx->max_d ) { + + if ( bit_depth == 16 ) { + val = row_pointers[y][channels*x] << 8; + val += row_pointers[y][(channels*x)+1]; + } + if ( bit_depth == 8 ) { + val = row_pointers[y][channels*x]; + } + + image[x + width*y] = val; + + } + } + } + + png_destroy_read_struct(&png_ptr, &info_ptr, &end_info); + fclose(fh); + + itrans_process_image(image, ctx, tilt); + + return 0; + +} diff --git a/src/readpng.h b/src/readpng.h new file mode 100644 index 0000000..ae76385 --- /dev/null +++ b/src/readpng.h @@ -0,0 +1,22 @@ +/* + * readpng.h + * + * Read images + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef READPNG_H +#define READPNG_H + +#include "control.h" +extern int readpng_read(const char *filename, double tilt, ControlContext *ctx); + +#endif /* READPNG_H */ diff --git a/src/reflections.c b/src/reflections.c new file mode 100644 index 0000000..dbe2a7d --- /dev/null +++ b/src/reflections.c @@ -0,0 +1,186 @@ +/* + * reflections.c + * + * Data structures in 3D space + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <math.h> +#include <stdio.h> + +#include "reflections.h" + +static void reflection_addfirst(ReflectionContext *reflectionctx) { + /* Create first items on lists - saves faffing later. Corresponds to a central marker. + Reflections are only stored if they have non-zero value. */ + reflectionctx->reflections = malloc(sizeof(Reflection)); + reflectionctx->reflections->next = NULL; + reflectionctx->reflections->x = 0; + reflectionctx->reflections->y = 0; + reflectionctx->reflections->z = 0; + reflectionctx->reflections->type = REFLECTION_CENTRAL; + reflectionctx->last_reflection = reflectionctx->reflections; +} + +/* Size of the volume (arbitrary units - m-1 in DTR), number of reflections across. The centre is defined to be at nx/2, ny/2, nz/2. + This defines the three-dimensional grid of reflections */ +ReflectionContext *reflection_init() { + + ReflectionContext *reflectionctx = malloc(sizeof(ReflectionContext)); + reflection_addfirst(reflectionctx); + + return reflectionctx; + +} + +void reflection_clear(ReflectionContext *reflectionctx) { + + Reflection *reflection = reflectionctx->reflections; + do { + Reflection *next = reflection->next; + free(reflection); + reflection = next; + } while ( reflection ); + + reflection_addfirst(reflectionctx); + +} + +void reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) { + + Reflection *new_reflection; + + new_reflection = malloc(sizeof(Reflection)); + new_reflection->next = NULL; + new_reflection->x = x; + new_reflection->y = y; + new_reflection->z = z; + new_reflection->intensity = intensity; + new_reflection->type = type; + + reflectionctx->last_reflection->next = new_reflection; + reflectionctx->last_reflection = new_reflection; + +} + +void reflection_add_index(ReflectionContext *reflectionctx, signed int h, signed int k, signed int l, double intensity, ReflectionType type) { + + Reflection *new_reflection; + + new_reflection = malloc(sizeof(Reflection)); + new_reflection->next = NULL; + new_reflection->h = h; + new_reflection->k = k; + new_reflection->l = l; + new_reflection->intensity = intensity; + new_reflection->type = type; + + reflectionctx->last_reflection->next = new_reflection; + reflectionctx->last_reflection = new_reflection; + +} + +/* x and y in metres (in real space!) */ +void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) { + + double nx, ny, nz; + + /* Real space */ + double L; + double d; + + /* Shared */ + double theta; + double psi; + + /* Reciprocal space */ + double r; + double tilt; + double omega; + double x_temp, y_temp, z_temp; + + tilt = 2*M_PI*(tilt_degrees/360); /* Convert to Radians */ + omega = 2*M_PI*(ctx->omega/360); /* Likewise */ + + d = sqrt((x*x) + (y*y)); + L = ctx->camera_length; + theta = atan2(d, L); + psi = atan2(y, x); + + r = 1/ctx->lambda; + x_temp = r*sin(theta)*cos(psi); + y_temp = -r*sin(theta)*sin(psi); /* Minus sign to define axes as y going upwards */ + z_temp = r- r*cos(theta); + + /* Apply the rotations... + First: rotate image clockwise until tilt axis is aligned horizontally. */ + nx = x_temp*cos(omega) + y_temp*-sin(omega); + ny = x_temp*sin(omega) + y_temp*cos(omega); + nz = z_temp; + + /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way. + This is because the crystal is rotated in the experiment, not the Ewald sphere. */ + x_temp = nx; y_temp = ny; z_temp = nz; + nx = x_temp; + ny = cos(tilt)*y_temp - sin(tilt)*z_temp; + nz = -sin(tilt)*y_temp - cos(tilt)*z_temp; + + reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); + +} + +/* Alternative formulation for when the input data is already in reciprocal space units + x and y in metres^-1 (in reciprocal space */ +void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) { + + double nx, ny, nz; + + /* Input (reciprocal space) */ + double dr; + + /* Shared */ + double theta; + double psi; + double r; + + /* Reciprocal space */ + double tilt; + double omega; + double x_temp, y_temp, z_temp; + + tilt = M_PI*(tilt_degrees/180); /* Convert to Radians */ + omega = M_PI*(ctx->omega/180); /* Likewise */ + + dr = sqrt((x*x) + (y*y)); + r = 1/ctx->lambda; + theta = atan2(dr, r); + psi = atan2(y, x); + + x_temp = r*sin(theta)*cos(psi); + y_temp = -r*sin(theta)*sin(psi); /* Minus sign to define axes as y going upwards */ + z_temp = r - r*cos(theta); + + /* Apply the rotations... + First: rotate image clockwise until tilt axis is aligned horizontally. */ + nx = x_temp*cos(omega) + y_temp*-sin(omega); + ny = x_temp*sin(omega) + y_temp*cos(omega); + nz = z_temp; + + /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way. + This is because the crystal is rotated in the experiment, not the Ewald sphere. */ + x_temp = nx; y_temp = ny; z_temp = nz; + nx = x_temp; + ny = cos(tilt)*y_temp - sin(tilt)*z_temp; + nz = -sin(tilt)*y_temp - cos(tilt)*z_temp; + + reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); + +} diff --git a/src/reflections.h b/src/reflections.h new file mode 100644 index 0000000..7104a27 --- /dev/null +++ b/src/reflections.h @@ -0,0 +1,77 @@ +/* + * reflections.h + * + * Data structures in 3D space + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef REFLECTION_H +#define REFLECTION_H + +typedef enum { + REFLECTION_NORMAL, /* Measured - expressed as x, y, z position */ + REFLECTION_CENTRAL, /* The central beam */ + REFLECTION_GENERATED /* Generated and intensity-measured - expressed as h, k, l index */ +} ReflectionType; + +typedef struct reflection_struct { + + double x; + double y; + double z; + double intensity; + + signed int h; + signed int k; + signed int l; + + ReflectionType type; + + struct reflection_struct *next; + +} Reflection; + +typedef struct { + + double x; + double y; + double z; + + double modulus; /* Convenience */ + +} Vector; + +typedef struct { + Vector a; + Vector b; + Vector c; +} Basis; + +typedef struct rctx_struct { + + Reflection *reflections; + Reflection *last_reflection; + + Basis *basis; + +} ReflectionContext; + +extern ReflectionContext *reflection_init(void); +void reflection_clear(ReflectionContext *reflectionctx); + +extern void reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type); +extern void reflection_add_index(ReflectionContext *reflectionctx, signed int h, signed int k, signed int l, double intensity, ReflectionType type); + +#include "control.h" + +extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity); +extern void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity); + +#endif /* REFLECTION_H */ diff --git a/src/trackball.c b/src/trackball.c new file mode 100644 index 0000000..f23d3db --- /dev/null +++ b/src/trackball.c @@ -0,0 +1,324 @@ +/* + * (c) Copyright 1993, 1994, Silicon Graphics, Inc. + * ALL RIGHTS RESERVED + * Permission to use, copy, modify, and distribute this software for + * any purpose and without fee is hereby granted, provided that the above + * copyright notice appear in all copies and that both the copyright notice + * and this permission notice appear in supporting documentation, and that + * the name of Silicon Graphics, Inc. not be used in advertising + * or publicity pertaining to distribution of the software without specific, + * written prior permission. + * + * THE MATERIAL EMBODIED ON THIS SOFTWARE IS PROVIDED TO YOU "AS-IS" + * AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE, + * INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR + * FITNESS FOR A PARTICULAR PURPOSE. IN NO EVENT SHALL SILICON + * GRAPHICS, INC. BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT, + * SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY + * KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION, + * LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF + * THIRD PARTIES, WHETHER OR NOT SILICON GRAPHICS, INC. HAS BEEN + * ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON + * ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE + * POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE. + * + * US Government Users Restricted Rights + * Use, duplication, or disclosure by the Government is subject to + * restrictions set forth in FAR 52.227.19(c)(2) or subparagraph + * (c)(1)(ii) of the Rights in Technical Data and Computer Software + * clause at DFARS 252.227-7013 and/or in similar or successor + * clauses in the FAR or the DOD or NASA FAR Supplement. + * Unpublished-- rights reserved under the copyright laws of the + * United States. Contractor/manufacturer is Silicon Graphics, + * Inc., 2011 N. Shoreline Blvd., Mountain View, CA 94039-7311. + * + * OpenGL(TM) is a trademark of Silicon Graphics, Inc. + */ +/* + * Trackball code: + * + * Implementation of a virtual trackball. + * Implemented by Gavin Bell, lots of ideas from Thant Tessman and + * the August '88 issue of Siggraph's "Computer Graphics," pp. 121-129. + * + * Vector manip code: + * + * Original code from: + * David M. Ciemiewicz, Mark Grossman, Henry Moreton, and Paul Haeberli + * + * Much mucking with by: + * Gavin Bell + */ +#include <math.h> +#include "trackball.h" + +/* + * This size should really be based on the distance from the center of + * rotation to the point on the object underneath the mouse. That + * point would then track the mouse as closely as possible. This is a + * simple example, though, so that is left as an Exercise for the + * Programmer. + */ +#define TRACKBALLSIZE (0.8) + +/* + * Local function prototypes (not defined in trackball.h) + */ +static float tb_project_to_sphere(float, float, float); +static void normalize_quat(float [4]); + +void +vzero(float *v) +{ + v[0] = 0.0; + v[1] = 0.0; + v[2] = 0.0; +} + +void +vset(float *v, float x, float y, float z) +{ + v[0] = x; + v[1] = y; + v[2] = z; +} + +void +vsub(const float *src1, const float *src2, float *dst) +{ + dst[0] = src1[0] - src2[0]; + dst[1] = src1[1] - src2[1]; + dst[2] = src1[2] - src2[2]; +} + +void +vcopy(const float *v1, float *v2) +{ + register int i; + for (i = 0 ; i < 3 ; i++) + v2[i] = v1[i]; +} + +void +vcross(const float *v1, const float *v2, float *cross) +{ + float temp[3]; + + temp[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]); + temp[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]); + temp[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]); + vcopy(temp, cross); +} + +float +vlength(const float *v) +{ + return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); +} + +void +vscale(float *v, float div) +{ + v[0] *= div; + v[1] *= div; + v[2] *= div; +} + +void +vnormal(float *v) +{ + vscale(v,1.0/vlength(v)); +} + +float +vdot(const float *v1, const float *v2) +{ + return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]; +} + +void +vadd(const float *src1, const float *src2, float *dst) +{ + dst[0] = src1[0] + src2[0]; + dst[1] = src1[1] + src2[1]; + dst[2] = src1[2] + src2[2]; +} + +/* + * Ok, simulate a track-ball. Project the points onto the virtual + * trackball, then figure out the axis of rotation, which is the cross + * product of P1 P2 and O P1 (O is the center of the ball, 0,0,0) + * Note: This is a deformed trackball-- is a trackball in the center, + * but is deformed into a hyperbolic sheet of rotation away from the + * center. This particular function was chosen after trying out + * several variations. + * + * It is assumed that the arguments to this routine are in the range + * (-1.0 ... 1.0) + */ +void +trackball(float q[4], float p1x, float p1y, float p2x, float p2y) +{ + float a[3]; /* Axis of rotation */ + float phi; /* how much to rotate about axis */ + float p1[3], p2[3], d[3]; + float t; + + if (p1x == p2x && p1y == p2y) { + /* Zero rotation */ + vzero(q); + q[3] = 1.0; + return; + } + + /* + * First, figure out z-coordinates for projection of P1 and P2 to + * deformed sphere + */ + vset(p1,p1x,p1y,tb_project_to_sphere(TRACKBALLSIZE,p1x,p1y)); + vset(p2,p2x,p2y,tb_project_to_sphere(TRACKBALLSIZE,p2x,p2y)); + + /* + * Now, we want the cross product of P1 and P2 + */ + vcross(p2,p1,a); + + /* + * Figure out how much to rotate around that axis. + */ + vsub(p1,p2,d); + t = vlength(d) / (2.0*TRACKBALLSIZE); + + /* + * Avoid problems with out-of-control values... + */ + if (t > 1.0) t = 1.0; + if (t < -1.0) t = -1.0; + phi = 2.0 * asin(t); + + axis_to_quat(a,phi,q); +} + +/* + * Given an axis and angle, compute quaternion. + */ +void +axis_to_quat(float a[3], float phi, float q[4]) +{ + vnormal(a); + vcopy(a,q); + vscale(q,sin(phi/2.0)); + q[3] = cos(phi/2.0); +} + +/* + * Project an x,y pair onto a sphere of radius r OR a hyperbolic sheet + * if we are away from the center of the sphere. + */ +static float +tb_project_to_sphere(float r, float x, float y) +{ + float d, t, z; + + d = sqrt(x*x + y*y); + if (d < r * 0.70710678118654752440) { /* Inside sphere */ + z = sqrt(r*r - d*d); + } else { /* On hyperbola */ + t = r / 1.41421356237309504880; + z = t*t / d; + } + return z; +} + +/* + * Given two rotations, e1 and e2, expressed as quaternion rotations, + * figure out the equivalent single rotation and stuff it into dest. + * + * This routine also normalizes the result every RENORMCOUNT times it is + * called, to keep error from creeping in. + * + * NOTE: This routine is written so that q1 or q2 may be the same + * as dest (or each other). + */ + +#define RENORMCOUNT 97 + +void +add_quats(float q1[4], float q2[4], float dest[4]) +{ + static int count=0; + float t1[4], t2[4], t3[4]; + float tf[4]; + + vcopy(q1,t1); + vscale(t1,q2[3]); + + vcopy(q2,t2); + vscale(t2,q1[3]); + + vcross(q2,q1,t3); + vadd(t1,t2,tf); + vadd(t3,tf,tf); + tf[3] = q1[3] * q2[3] - vdot(q1,q2); + + dest[0] = tf[0]; + dest[1] = tf[1]; + dest[2] = tf[2]; + dest[3] = tf[3]; + + if (++count > RENORMCOUNT) { + count = 0; + normalize_quat(dest); + } +} + +/* + * Quaternions always obey: a^2 + b^2 + c^2 + d^2 = 1.0 + * If they don't add up to 1.0, dividing by their magnitued will + * renormalize them. + * + * Note: See the following for more information on quaternions: + * + * - Shoemake, K., Animating rotation with quaternion curves, Computer + * Graphics 19, No 3 (Proc. SIGGRAPH'85), 245-254, 1985. + * - Pletinckx, D., Quaternion calculus as a basic tool in computer + * graphics, The Visual Computer 5, 2-13, 1989. + */ +static void +normalize_quat(float q[4]) +{ + int i; + float mag; + + mag = (q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); + for (i = 0; i < 4; i++) q[i] /= mag; +} + +/* + * Build a rotation matrix, given a quaternion rotation. + * + */ +void +build_rotmatrix(float m[4][4], float q[4]) +{ + m[0][0] = 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2]); + m[0][1] = 2.0 * (q[0] * q[1] - q[2] * q[3]); + m[0][2] = 2.0 * (q[2] * q[0] + q[1] * q[3]); + m[0][3] = 0.0; + + m[1][0] = 2.0 * (q[0] * q[1] + q[2] * q[3]); + m[1][1]= 1.0 - 2.0 * (q[2] * q[2] + q[0] * q[0]); + m[1][2] = 2.0 * (q[1] * q[2] - q[0] * q[3]); + m[1][3] = 0.0; + + m[2][0] = 2.0 * (q[2] * q[0] - q[1] * q[3]); + m[2][1] = 2.0 * (q[1] * q[2] + q[0] * q[3]); + m[2][2] = 1.0 - 2.0 * (q[1] * q[1] + q[0] * q[0]); + m[2][3] = 0.0; + + m[3][0] = 0.0; + m[3][1] = 0.0; + m[3][2] = 0.0; + m[3][3] = 1.0; +} + diff --git a/src/trackball.h b/src/trackball.h new file mode 100644 index 0000000..b676fb4 --- /dev/null +++ b/src/trackball.h @@ -0,0 +1,78 @@ +/* + * (c) Copyright 1993, 1994, Silicon Graphics, Inc. + * ALL RIGHTS RESERVED + * Permission to use, copy, modify, and distribute this software for + * any purpose and without fee is hereby granted, provided that the above + * copyright notice appear in all copies and that both the copyright notice + * and this permission notice appear in supporting documentation, and that + * the name of Silicon Graphics, Inc. not be used in advertising + * or publicity pertaining to distribution of the software without specific, + * written prior permission. + * + * THE MATERIAL EMBODIED ON THIS SOFTWARE IS PROVIDED TO YOU "AS-IS" + * AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE, + * INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR + * FITNESS FOR A PARTICULAR PURPOSE. IN NO EVENT SHALL SILICON + * GRAPHICS, INC. BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT, + * SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY + * KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION, + * LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF + * THIRD PARTIES, WHETHER OR NOT SILICON GRAPHICS, INC. HAS BEEN + * ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON + * ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE + * POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE. + * + * US Government Users Restricted Rights + * Use, duplication, or disclosure by the Government is subject to + * restrictions set forth in FAR 52.227.19(c)(2) or subparagraph + * (c)(1)(ii) of the Rights in Technical Data and Computer Software + * clause at DFARS 252.227-7013 and/or in similar or successor + * clauses in the FAR or the DOD or NASA FAR Supplement. + * Unpublished-- rights reserved under the copyright laws of the + * United States. Contractor/manufacturer is Silicon Graphics, + * Inc., 2011 N. Shoreline Blvd., Mountain View, CA 94039-7311. + * + * OpenGL(TM) is a trademark of Silicon Graphics, Inc. + */ +/* + * trackball.h + * A virtual trackball implementation + * Written by Gavin Bell for Silicon Graphics, November 1988. + */ + +/* + * Pass the x and y coordinates of the last and current positions of + * the mouse, scaled so they are from (-1.0 ... 1.0). + * + * The resulting rotation is returned as a quaternion rotation in the + * first paramater. + */ +void +trackball(float q[4], float p1x, float p1y, float p2x, float p2y); + +/* + * Given two quaternions, add them together to get a third quaternion. + * Adding quaternions to get a compound rotation is analagous to adding + * translations to get a compound translation. When incrementally + * adding rotations, the first argument here should be the new + * rotation, the second and third the total rotation (which will be + * over-written with the resulting new total rotation). + */ +void +add_quats(float *q1, float *q2, float *dest); + +/* + * A useful function, builds a rotation matrix in Matrix based on + * given quaternion. + */ +void +build_rotmatrix(float m[4][4], float q[4]); + +/* + * This function computes a quaternion based on an axis (defined by + * the given vector) and an angle about which to rotate. The angle is + * expressed in radians. The result is put into the third argument. + */ +void +axis_to_quat(float a[3], float phi, float q[4]); + diff --git a/src/utils.c b/src/utils.c new file mode 100644 index 0000000..ce7c7f3 --- /dev/null +++ b/src/utils.c @@ -0,0 +1,44 @@ +/* + * utils.c + * + * Utility stuff + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#include <math.h> + +/* Return the MOST POSITIVE of two numbers */ +unsigned int biggest(signed int a, signed int b) { + if ( a>b ) { + return a; + } + return b; +} + +/* Return the LEAST POSITIVE of two numbers */ +unsigned int smallest(signed int a, signed int b) { + if ( a<b ) { + return a; + } + return b; +} + +double distance(double x1, double y1, double x2, double y2) { + return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); +} + +double modulus(double x, double y, double z) { + return sqrt(x*x + y*y + z*z); +} + +/* Angle between two vectors. Answer in degrees */ +double angle_between(double x1, double y1, double z1, double x2, double y2, double z2) { + + double mod1 = modulus(x1, y1, z1); + double mod2 = modulus(x2, y2, z2); + return ((acos((x1*x2 + y1*y2 + z1*z2) / (mod1*mod2)))/M_PI) * 180; + +} diff --git a/src/utils.h b/src/utils.h new file mode 100644 index 0000000..10ef3eb --- /dev/null +++ b/src/utils.h @@ -0,0 +1,24 @@ +/* + * utils.h + * + * Utility stuff + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef UTILS_H +#define UTILS_H + +extern unsigned int biggest(signed int a, signed int b); +extern unsigned int smallest(signed int a, signed int b); +extern double distance(double x1, double y1, double x2, double y2); +extern double modulus(double x, double y, double z); +extern double angle_between(double x1, double y1, double z1, double x2, double y2, double z2); + +#endif /* UTILS_H */ |