Initial import

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@1 bf6ca9ba-c028-0410-8290-897cf20841d1

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 */