Remove vestigial geoptimiser stuff

7 files changed, 2 insertions, 3110 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 27b73f5e..97494b93 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -73,20 +73,11 @@ if (GTK_FOUND)
   add_custom_target(crystfel-gresources
                     DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/crystfel-gresources.c)
 
-endif ()
-
-# If no version of GTK was found, try for GDK
-if (NOT GTK_FOUND)
+else()
   message(STATUS "GTK3 not found. GUI parts will not be compiled")
-  message(STATUS "Looking separately for GDK")
-  pkg_search_module(GDK gdk-3.0)
-  if (NOT GDK_FOUND)
-    message(STATUS "GDK not found.")
-  endif()
-endif()
+endif ()
 
 pkg_search_module(CAIRO cairo)
-pkg_search_module(GDKPIXBUF gdk-pixbuf-2.0)
 
 include(CheckCCompilerFlag)
 check_c_compiler_flag("-fdiagnostics-color=always" HAVE_DIAG_COLOR)
@@ -126,8 +117,6 @@ include(CheckLibraryExists)
 set(HAVE_CAIRO ${CAIRO_FOUND})
 set(HAVE_GTK ${GTK_FOUND})
 set(HAVE_OPENCL ${OpenCL_FOUND})
-set(HAVE_GDKPIXBUF ${GDKPIXBUF_FOUND})
-set(HAVE_GDK ${GDK_FOUND})
 set(HAVE_ZMQ ${ZMQ_FOUND})
 set(HAVE_HDF5 1)
 set(HAVE_ASAPO ${ASAPO_FOUND})
@@ -347,40 +336,6 @@ target_link_libraries(ambigator ${COMMON_LIBRARIES} ${HDF5_C_LIBRARIES})
 list(APPEND CRYSTFEL_EXECUTABLES ambigator)
 
 # ----------------------------------------------------------------------
-# geoptimiser
-
-      # FIXME!
-#if (GDKPIXBUF_FOUND)
-#  target_include_directories(${PROJECT_NAME} PRIVATE ${GDKPIXBUF_INCLUDE_DIRS})
-#  target_link_libraries(${PROJECT_NAME} PRIVATE ${GDKPIXBUF_LIBRARIES})
-#endif (GDKPIXBUF_FOUND)
-
-# FIXME: Restore!
-#set(GEOPTIMISER_SOURCES src/geoptimiser.c)
-#add_executable(geoptimiser ${GEOPTIMISER_SOURCES}
-#               ${CMAKE_CURRENT_BINARY_DIR}/version.c)
-#target_include_directories(geoptimiser PRIVATE ${COMMON_INCLUDES})
-#target_link_libraries(geoptimiser ${COMMON_LIBRARIES})
-#list(APPEND CRYSTFEL_EXECUTABLES geoptimiser)
-#
-## Add features one by one so that #ifdef HAVE_XX//#include XX.h always works
-## If Cairo, gdk-pixbuf and GDK are all found, geoptimiser will add PNG support
-#if (CAIRO_FOUND)
-#  target_include_directories(geoptimiser PRIVATE ${CAIRO_INCLUDE_DIRS})
-#  target_link_libraries(geoptimiser ${CAIRO_LIBRARIES})
-#endif (CAIRO_FOUND)
-#
-#if (GDKPIXBUF_FOUND)
-#  target_include_directories(geoptimiser PRIVATE ${GDKPIXBUF_INCLUDE_DIRS})
-#  target_link_libraries(geoptimiser ${GDKPIXBUF_LIBRARIES})
-#endif (GDKPIXBUF_FOUND)
-#
-#if (GDK_FOUND)
-#  target_include_directories(geoptimiser PRIVATE ${GDK_INCLUDE_DIRS})
-#  target_link_libraries(geoptimiser ${GDK_LIBRARIES})
-#endif (GDK_FOUND)
-
-# ----------------------------------------------------------------------
 # whirligig
 
 set(WHIRLIGIG_SOURCES src/whirligig.c)
@@ -453,7 +408,6 @@ install(FILES
 	doc/man/cell_explorer.1
 	doc/man/check_hkl.1
 	doc/man/compare_hkl.1
-	doc/man/geoptimiser.1
 	doc/man/get_hkl.1
 	doc/man/indexamajig.1
 	doc/man/list_events.1
diff --git a/config.h.cmake.in b/config.h.cmake.in
index c056cc4a..9e48d674 100644
--- a/config.h.cmake.in
+++ b/config.h.cmake.in
@@ -1,8 +1,6 @@
 /* config.h for CrystFEL main programs */
 
 #cmakedefine HAVE_GTK
-#cmakedefine HAVE_GDKPIXBUF
-#cmakedefine HAVE_GDK
 #cmakedefine HAVE_CAIRO
 #cmakedefine HAVE_OPENCL
 #cmakedefine HAVE_CL_CL_H
diff --git a/config.h.in b/config.h.in
index afa44161..d0b6e9c2 100644
--- a/config.h.in
+++ b/config.h.in
@@ -1,8 +1,6 @@
 /* config.h for CrystFEL main programs */
 
 #mesondefine HAVE_GTK
-#mesondefine HAVE_GDKPIXBUF
-#mesondefine HAVE_GDK
 #mesondefine HAVE_CAIRO
 #mesondefine HAVE_OPENCL
 #mesondefine HAVE_CL_CL_H
diff --git a/doc/man/geoptimiser.1 b/doc/man/geoptimiser.1
deleted file mode 100644
index b85bfe23..00000000
--- a/doc/man/geoptimiser.1
+++ /dev/null
@@ -1,153 +0,0 @@
-.\"
-.\" geoptimiser man page
-.\"
-.\" Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
-.\"                       a research centre of the Helmholtz Association.
-.\"
-.\" Part of CrystFEL - crystallography with a FEL
-.\"
-
-.TH GEOPTIMISER 1
-.SH NAME
-geoptimiser \- detector geometry refinement
-.SH SYNOPSIS
-.PP
-.BR geoptimiser
-\fB-i \fIinput.stream \fB-g \fIinput.geom \fB-o \fIoutput.geom \fB-c \fIconnected_rigidgroup_coll \fB-q \fI\quadrant_rigidgroup_coll\fR
-[\fBoptions\fR]
-.PP
-\fBgeoptimiser --help\fR
-
-.SH DESCRIPTION
-
-\fBgeoptimiser\fR refines and optimizes the detector geometry by comparing the locations of observed Bragg peaks in a set of indexed patterns with the spot locations predicted from the crystal indexing procedure. It can refine position, rotation and distance of each panel relative to the interaction region. It requires a stream file with indexed patterns, a geometry file with the detector geometry to refine, and some parameters to specify which panels are physically connected to each other and which are attached to the same physical support (for example, all panels in a quadrant of the CSPAD detector). The output is a geometry file with the optimized detector geometry and a set of diagnostic error maps in HDF5 format.  Several options are available to tweak the number of peaks included in the optimization procedure based on a range of criteria.
-For a complete description of the optimization algorithm, see the following paper:
-
-.IP
-O. Yefanov, V. Mariani, C. Gati, T. A. White, H. N. Chapman, and A. Barty. "Accurate determination of segmented X-ray detector geometry". Optics Express 23 (2015) 28459. doi:10.1364/OE.23.028459.
-
-.PP
-For minimal basic use, you need to provide a stream file with diffraction patterns, a filename for the output optimized geometry, and the name of two rigid group collections defined in the geometry file: one describing which panels in the detector are physically connected (and hence whose geometry should be optimized as if they were a single panel), and one to describe which panels are attached to the same underlying support (whose position and orientation are likely to be correlated).
-
-.PP
-See \fBman crystfel_geometry\fR for information on how to create a file describing the detector geometry, and guidelines to define the required rigid groups and rigid groups collections.
-
-.PP
-Geoptimizer saves error maps before and after the geometry optimization. These maps show an estimation of the geometry error as average displacement, for each pixel, between the predicted and the observed positions of the Bragg peaks. The maps are color-scaled PNG files and
-are named "error_map_before.png" and "error_map_after.png" respectively. In order to better visualize small local errors, the color range has been set to show all displacements bigger that 1 pixel as white (in other words, "displacement saturation" is set at 1 pixel).
-
-
-.SH OPTIONS
-.PD 0
-.IP "\fB-i\fR \fIfilename\fR"
-.IP \fB--input=\fR\fIfilename\fR
-.PD
-Give the filename of the stream from which to read the indexed patterns.
-
-.PD 0
-.IP "\fB-g\fR \fIfilename\fR"
-.IP \fB--geometry=\fR\fIfilename\fR
-.PD
-Read the detector geometry to optimize from \fIfilename\fR.  If this option is omitted, the geometry file from the stream's header will be used (see \fB--input\fR), if present.
-
-.PD 0
-.IP "\fB-o\fR \fIfilename\fR"
-.IP \fB--output=\fR\fIfilename\fR
-.PD
-Write the optimized detector geometry to \fIfilename\fR.
-
-.PD 0
-.IP "\fB-c \fIname\fR"
-.IP \fB--connected=\fIname\fR
-.PD
-Specifies the name of the rigid group collection for connected panels.  This rigid group collection describes how the panels are physically connected in the detector.
-A set of rigid groups must be defined in the geometry file, with each group containing only panels that are physically connected to each other (for example the pairs of ASICs sharing the same piece of detector silicon in the CSPAD detector).
-.sp
-If a given rigid group is a member of \fIname\fR, then panels which are members of that rigid group will be kept strictly in the same relative position and orientation relative to one another.
-
-.PD 0
-.IP "\fB-q\fR \fIname\fR"
-.IP \fB--quadrants=\fR\fIname\fR
-.PD
-Specifies the name of the rigid group collection for detector 'quadrants'.  This rigid group collection describes how panels are connected to the underlying support of the detector, for example, the panels belonging to the same quadrant of the CSPAD detector.
-.sp
-If a given rigid group is a member of \fIname\fR, then panels which are members of that rigid group and which do not contain enough peaks for positional refinement will be moved according to the average movement of the other panels in the group.
-
-.PD 0
-.IP \fB--no-error-maps\fR
-.PD
-Forces geoptimiser not to save error maps.
-
-.PD 0
-.IP "\fB-x\fR \fIn\fR"
-.IP \fB--min-num-peaks-per-pixel=\fR\fIn\fR
-.PD
-Sets the minimum number of peaks that should fall within a pixel, across all indexed patterns, to contribute to the geometry optimization. The default value is 3.
-
-.PD 0
-.IP "\fB-p\fR \fIn\fR"
-.IP \fB--min-num-pixels-per-conn-group=\fR\fIn\fR
-.PD
-Sets the minimum number of pixels that contribute to the geometry opimization (see the \fB--min-num-peaks-per-pixel\fR option above ) that a connected group should have for the group to be optimized independently. The default value is 100.
-
-.PD 0
-.IP \fB--min-num-peaks-per-panel=\fR\fIn\fR
-.PD
-This option has been renamd to  \fB--min-num-pixels-per-conn-group\fR. It has been deprecated and will soon be removed. It is currently mapped to the \fB--min-num-pixels-per-conn-group\fR option.
-
-.PD 0
-.IP "\fB-l\fR"
-.IP \fB--most-freq-clen\fR
-.PD
-Some stream files can contain patterns collected using different camera lengths. By default, detector distance is optimized using only patterns collected with the most frequent camera length in the stream file.
-However, all patterns are used to compute panel shifts and rotations. With this option, shifts are rotation are computed with the same subset of patterns used for the detector distance optimization.
-
-.PD 0
-.IP "\fB-s\fR"
-.IP \fB--individual-dist-offset\fR
-.PD
-By default, geoptimiser optimizes detector distance assuming that all panels have the same distance from the sample (i.e., a single distance is optimized for the whole detector). With this option, each panel's
-distance is instead optimized independently.
-
-.PD 0
-.IP "\fB-m\fR \fIdist\fR"
-.IP \fB--max-peak-dist=\fR\fIdist\fR
-.PD
-Geoptimiser refines detector geometry by comparing the predicted position of Bragg peaks with the location of detected peak in indexed patterns. This option sets the maximum distance in pixels between the predicted and the observed peaks for the pair
-to be included in the optimization process. The default maximum distance is half of the minimum distance between Bragg peaks derived from the data.
-
-.PD 0
-.IP \fB--no-stretch\fR
-.PD
-By default, geoptimiser refines the distance between the detector and the sample. This option turns off this optimization.
-
-.PD 0
-.IP \fB--no-cspad\fR
-.PD
-If a geometry file containing 64 panels (ASICs) is provided by the user, geoptimiser assumes that the detector described by the geometry file is a CSPAD, and performs some sanity-checks on the relative distance and orientation of connected panels. If the checks fail, the geometry optimization is stopped. This option turns off these checks.
-
-.PD 0
-.IP \fB--enforce-cspad-layout\fR
-.PD
-When this option is used, geoptimiser tries to fix the CSPAD layout problems detected by the checks described in the entry for the \fB--no-cspad\fR option. Immediately after performing this operation, geoptimser saves the new detector layout in the output refined geometry file and exits.
-
-.SH AUTHOR
-This page was written by Valerio Mariani, Oleksandr Yefanov and Thomas White.
-
-.SH REPORTING BUGS
-Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
-
-.SH COPYRIGHT AND DISCLAIMER
-Copyright © 2014-2021 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
-.P
-geoptimiser, and this manual, are part of CrystFEL.
-.P
-CrystFEL is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-.P
-CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-.P
-You should have received a copy of the GNU General Public License along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
-
-.SH SEE ALSO
-.BR crystfel (7),
-.BR crystfel_geometry (5)
diff --git a/libcrystfel/libcrystfel-config.h.cmake.in b/libcrystfel/libcrystfel-config.h.cmake.in
index 82d0611a..0b34950f 100644
--- a/libcrystfel/libcrystfel-config.h.cmake.in
+++ b/libcrystfel/libcrystfel-config.h.cmake.in
@@ -7,7 +7,6 @@
 #cmakedefine HAVE_FDIP
 #cmakedefine HAVE_ZLIB
 #cmakedefine HAVE_GZBUFFER
-#cmakedefine HAVE_GDKPIXBUF
 #cmakedefine HAVE_LIBCCP4
 #cmakedefine HAVE_MSGPACK
 #cmakedefine HAVE_CLOCK_GETTIME
diff --git a/meson.build b/meson.build
index 59bbcfb7..21f4534a 100644
--- a/meson.build
+++ b/meson.build
@@ -43,16 +43,6 @@ if gtkdep.found()
   conf_data.set10('HAVE_GTK', true)
 endif
 
-gdkdep = dependency('gdk-3.0', required: false)
-if gdkdep.found()
-  conf_data.set10('HAVE_GDK', true)
-endif
-
-gdkpixbufdep = dependency('gdk-pixbuf-2.0', required: false)
-if gdkpixbufdep.found()
-  conf_data.set10('HAVE_GDKPIXBUF', true)
-endif
-
 cairodep = dependency('cairo', required: false)
 if cairodep.found()
   conf_data.set10('HAVE_CAIRO', true)
@@ -193,13 +183,6 @@ if hdf5dep.found()
              install_rpath: '$ORIGIN/../lib64/:$ORIGIN/../lib')
 endif
 
-# geoptimiser
-# FIXME: restore
-#executable('geoptimiser',
-#           ['src/geoptimiser.c', 'src/hdfsee-render.c', versionc],
-#           dependencies: [mdep, libcrystfeldep, gsldep, gdkpixbufdep, gdkdep],
-#           install: true,
-#           install_rpath: '$ORIGIN/../lib64/:$ORIGIN/../lib')
 
 # CrystFEL GUI
 if gtkdep.found()
@@ -300,7 +283,6 @@ install_man(['doc/man/ambigator.1',
              'doc/man/compare_hkl.1',
              'doc/man/crystfel.7',
              'doc/man/crystfel_geometry.5',
-             'doc/man/geoptimiser.1',
              'doc/man/get_hkl.1',
              'doc/man/indexamajig.1',
              'doc/man/list_events.1',
diff --git a/src/geoptimiser.c b/src/geoptimiser.c
deleted file mode 100644
index f28e9551..00000000
--- a/src/geoptimiser.c
+++ /dev/null
@@ -1,2886 +0,0 @@
-/*
- * geoptimiser.c
- *
- * Refine detector geometry
- *
- * Copyright © 2014-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2014-2015 Oleksandr Yefanov
- *   2014-2015 Valerio Mariani
- *   2014-2020 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <getopt.h>
-#include <math.h>
-#include <ctype.h>
-#include <time.h>
-#include <float.h>
-#include <assert.h>
-#include <gsl/gsl_statistics_double.h>
-#include <gsl/gsl_sort.h>
-
-#ifdef HAVE_CAIRO
-#ifdef HAVE_GDK
-#ifdef HAVE_GDKPIXBUF
-#include <cairo.h>
-#include <gdk-pixbuf/gdk-pixbuf.h>
-#include <gdk/gdk.h>
-#define CAN_SAVE_TO_PNG
-#endif /* HAVE_GDKPIXBUF */
-#endif /* HAVE_GDK */
-#endif /* HAVE_CAIRO */
-
-#include <detgeom.h>
-#include <datatemplate.h>
-#include <stream.h>
-#include <crystal.h>
-#include <image.h>
-#include <utils.h>
-#include <colscale.h>
-
-#include "version.h"
-
-
-struct imagefeature;
-
-static void show_help(const char *s)
-{
-	printf("Syntax: %s -i input.stream -g input.geom -o refined.geom "
-	       "-c connected_rgcollection -q quadrant_rgcollection [options]\n",
-	       s);
-	printf(
-"Refines detector geometry.\n"
-"\n"
-"  -h, --help                                   Display this help message.\n"
-"\n"
-"      --version                                Print CrystFEL version number and exit.\n"
-"  -i, --input=<filename>                       Input stream\n"
-"  -g. --geometry=<file>                        Input geometry file (if omitted: from stream).\n"
-"  -o, --output=<filename>                      Output geometry file.\n"
-"  -q, --quadrants=<rg_coll>                    Rigid group collection for quadrants.\n"
-"  -c, --connected=<rg_coll>                    Rigid group collection for connected ASICs.\n"
-"      --no-error-maps                          Do not generate error map PNGs.\n"
-"      --stretch-map                            Generate stretch map PNG (panels distance).\n"
-"  -x, --min-num-peaks-per-pixel=<num>          Minimum number of peaks per pixel.\n"
-"	                                         Default: 3. \n"
-"  -z, --max-num-peaks-per-pixel=<num>          Maximum number of peaks per pixel.\n"
-"	                                         Default is num_patterns / 10. \n"
-"  -p, --min-num-pixels-per-conn-group=<num>    Minimum number of useful pixels per\n"
-"                                                connected group.  Default: 100.\n"
-"  -l, --most-freq-clen                         Use only the most frequent camera length\n"
-"  -s, --individual-dist-offset                 Use a distance offset for each panel.\n"
-"                                                Default: whole-detector offset.\n"
-"      --no-stretch                             Do not optimize distance offset.\n"
-"                                                Default: distance offset is optimized\n"
-"      --no-rotation                            Do not optimize rotation of connectedd groups.\n"
-"                                                Default: rotation is optimized\n"
-"  -m  --max-peak-dist=<num>                    Maximum distance between predicted and\n"
-"                                                detected peaks (in pixels)\n"
-"                                                Default: half of minimal inter-Bragg distance\n"
-);
-}
-
-struct geoptimiser_params
-{
-	char *outfile;
-	int min_num_peaks_per_pix;
-	int max_num_peaks_per_pix;
-	int min_num_pix_per_conn_group;
-	int only_best_distance;
-	int nostretch;
-	int norotation;
-	int individual_coffset;
-	int error_maps;
-	int stretch_map;
-	int enforce_cspad_layout;
-	int no_cspad;
-	double max_peak_dist;
-	const char *command_line;
-};
-
-
-struct connected_data
-{
-	double sh_x;
-	double sh_y;
-	double cang;
-	double cstr;
-	int num_quad;
-	int num_peaks_per_pixel;
-	unsigned int n_peaks_in_conn;
-	char *name;
-};
-
-
-struct single_pixel_displ
-{
-	double dx;
-	double dy;
-	struct single_pixel_displ *ne;
-};
-
-
-struct gpanel
-{
-	struct detgeom_panel       *p;
-
-	/* Individual pixel displacements */
-	struct single_pixel_displ  *pix_displ_list;
-	struct single_pixel_displ **curr_pix_displ;
-	int                        *num_pix_displ;
-
-	/* Average displacements for each pixel */
-	double                     *avg_displ_x;
-	double                     *avg_displ_y;
-	double                     *avg_displ_abs;
-};
-
-
-static void free_pixbuf_data(guchar *data, gpointer p)
-{
-	free(data);
-}
-
-
-static GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale)
-{
-	guchar *data;
-	size_t x, y;
-	int max;
-
-	data = malloc(3*w*h);
-	if ( data == NULL ) return NULL;
-
-	max = h-(h/6);
-
-	for ( y=0; y<h; y++ ) {
-
-		double r, g, b;
-		int val;
-
-		val = y-(h/6);
-
-		colscale_lookup(val, max, scale, &r, &g, &b);
-
-		data[3*( 0+w*(h-1-y) )+0] = 0;
-		data[3*( 0+w*(h-1-y) )+1] = 0;
-		data[3*( 0+w*(h-1-y) )+2] = 0;
-		for ( x=1; x<w; x++ ) {
-			data[3*( x+w*(h-1-y) )+0] = 255*r;
-			data[3*( x+w*(h-1-y) )+1] = 255*g;
-			data[3*( x+w*(h-1-y) )+2] = 255*b;
-		}
-	}
-
-	y = h/6;
-	for ( x=1; x<w; x++ ) {
-		data[3*( x+w*(h-1-y) )+0] = 255;
-		data[3*( x+w*(h-1-y) )+1] = 255;
-		data[3*( x+w*(h-1-y) )+2] = 255;
-	}
-
-	return gdk_pixbuf_new_from_data(data, GDK_COLORSPACE_RGB,
-	                                FALSE, 8, w, h, w*3,
-	                                free_pixbuf_data, NULL);
-}
-
-
-static void compute_x_y(double fs, double ss, struct detgeom_panel *p,
-                        double *x, double *y)
-{
-	double xs, ys;
-
-	xs = fs*p->fsx + ss*p->ssx;
-	ys = fs*p->fsy + ss*p->ssy;
-
-	*x = xs + p->cnx;
-	*y = ys + p->cny;
-}
-
-
-static Reflection *find_closest_reflection(struct image *image,
-                                           double fx, double fy,
-                                           double *d, double *det_shift)
-{
-	double dmin = HUGE_VAL;
-	Reflection *closest = NULL;
-	int i;
-
-	for ( i=0; i<image->n_crystals; i++ ) {
-
-		Reflection *refl;
-		RefListIterator *iter;
-		RefList *rlist = crystal_get_reflections(image->crystals[i]);
-
-		for ( refl = first_refl(rlist, &iter);
-		      refl != NULL;
-		      refl = next_refl(refl, iter) )
-		{
-			double ds;
-			double rfs, rss;
-			double rx, ry;
-			int pn;
-
-			get_detector_pos(refl, &rfs, &rss);
-			pn = get_panel_number(refl);
-
-			compute_x_y(rfs, rss,
-			            &image->detgeom->panels[pn],
-			            &rx, &ry);
-
-			ds = distance(rx+det_shift[0], ry+det_shift[1], fx, fy);
-
-			if ( ds < dmin ) {
-				dmin = ds;
-				closest = refl;
-			}
-
-		}
-
-	}
-
-	if ( closest == NULL ) {
-		*d = +INFINITY;
-	} else {
-		*d = dmin;
-	}
-	return closest;
-}
-
-
-static double avg_clen(struct detgeom *det)
-{
-	int i;
-	double av = 0.0;
-
-	for ( i=0; i<det->n_panels; i++ ) {
-		av += det->panels[i].cnz * det->panels[i].pixel_pitch;
-	}
-	return av / i;
-}
-
-
-static struct image **read_patterns_from_stream(Stream *st,
-                                                const DataTemplate *dtempl,
-                                                int *n)
-{
-	struct image **images;
-	int n_chunks = 0;
-	int max_images = 1024;
-	int n_images = 0;
-
-	images = malloc(max_images * sizeof(struct image *));
-	if ( images == NULL ) {
-		ERROR("Failed to allocate memory for images.\n");
-		return NULL;
-	}
-
-	do {
-
-		images[n_images] = stream_read_chunk(st, dtempl,
-		                                     STREAM_REFLECTIONS
-		                                   | STREAM_PEAKS
-		                                   | STREAM_UNITCELL);
-		if ( images[n_images] == NULL ) break;
-
-		n_chunks++; /* Number of chunks processed */
-
-		/* Reject if there are no crystals (not indexed) */
-		if ( images[n_images]->n_crystals == 0 ) continue;
-
-		n_images++;  /* Number of images accepted */
-
-		if ( n_images == max_images ) {
-
-			struct image **images_new;
-
-			images_new = realloc(images,
-			      (max_images+1024)*sizeof(struct image *));
-			if ( images_new == NULL ) {
-				ERROR("Failed to allocate memory for "
-				      "patterns.\n");
-				free(images);
-				return NULL;
-			}
-
-			max_images += 1024;
-			images = images_new;
-		}
-
-		if ( n_images % 1000 == 0 ) {
-			STATUS("Loaded %i indexed patterns from %i total "
-			       "patterns.\n", n_images, n_chunks);
-		}
-
-
-	} while ( 1 );
-
-	*n = n_images;
-
-	STATUS("Found %i indexed patterns in stream (from a total of %i).\n",
-	       n_images, n_chunks);
-
-	return images;
-}
-
-
-static struct rvec get_q_from_xyz(double rx, double ry, double dist, double l)
-{
-
-	struct rvec q;
-	double r = sqrt(rx*rx + ry*ry);
-	double twotheta = atan2(r, dist);
-	double az = atan2(ry, rx);
-
-	q.u = 1.0/l * sin(twotheta)*cos(az);
-	q.v = 1.0/l * sin(twotheta)*sin(az);
-	q.w = 1.0/l * (cos(twotheta) - 1.0);
-
-	return q;
-}
-
-
-static UnitCell *compute_avg_cell_parameters(struct image **images,
-                                             int n)
-{
-	int numavc;
-	int j, i;
-	double minc[6];
-	double maxc[6];
-	double avg_cpar[6] = {0, 0, 0, 0, 0, 0};
-	UnitCell *avg;
-
-	for ( j=0; j<6; j++ ) {
-		minc[j] = 1e100;
-		maxc[j] = -1e100;
-	}
-
-	numavc = 0;
-	for ( i=0; i<n; i++ ) {
-
-		struct image *image;
-		double cpar[6];
-		int j, cri;
-
-		image = images[i];
-
-		for ( cri=0; cri<image->n_crystals; cri++ ) {
-
-			UnitCell *cell = crystal_get_cell(image->crystals[cri]);
-
-			cell_get_parameters(cell,
-			                    &cpar[0],  // a
-			                    &cpar[1],  // b
-			                    &cpar[2],  // c
-			                    &cpar[3],  // alpha
-			                    &cpar[4],  // beta
-			                    &cpar[5]); // gamma
-
-			for ( j=0; j<6; j++ ) {
-				avg_cpar[j] += cpar[j];
-				if ( cpar[j]<minc[j] ) minc[j] = cpar[j];
-				if ( cpar[j]>maxc[j] ) maxc[j] = cpar[j];
-			}
-			numavc++;
-
-		}
-
-	}
-
-	if ( numavc > 0 ) {
-		for ( j=0; j<6; j++ ) avg_cpar[j] /= numavc;
-	}
-
-	avg = cell_new();
-
-	cell_set_parameters(avg, avg_cpar[0], avg_cpar[1], avg_cpar[2],
-	                         avg_cpar[3], avg_cpar[4], avg_cpar[5]);
-
-	STATUS("Average cell parameters:\n");
-	STATUS("Average a, b, c (in A): %6.3f, %6.3f, %6.3f\n",
-	       avg_cpar[0]*1e10, avg_cpar[1]*1e10, avg_cpar[2]*1e10);
-	STATUS("Minimum -Maximum a, b, c:\n"
-	       "\t%6.3f - %6.3f A,\n"
-	       "\t%6.3f - %6.3f A,\n"
-	       "\t%6.3f - %6.3f A\n",
-	       minc[0]*1e10, maxc[0]*1e10, minc[1]*1e10,
-	       maxc[1]*1e10, minc[2]*1e10, maxc[2]*1e10);
-	STATUS("Average alpha,beta,gamma in degrees: %6.3f, %6.3f, %6.3f\n",
-	       rad2deg(avg_cpar[3]), rad2deg(avg_cpar[4]), rad2deg(avg_cpar[5]));
-	STATUS("Minimum - Maximum alpha,beta,gamma in degrees:\n"
-	       "\t%5.2f - %5.2f,\n"
-	       "\t%5.2f - %5.2f,\n"
-	       "\t%5.2f - %5.2f\n",
-	       rad2deg(minc[3]), rad2deg(maxc[3]),
-	       rad2deg(minc[4]), rad2deg(maxc[4]),
-	       rad2deg(minc[5]), rad2deg(maxc[5]));
-
-	return avg;
-}
-
-
-static double pick_clen_to_use(struct geoptimiser_params *gparams,
-                               struct image **images, int n,
-                               double avg_res, UnitCell *avg)
-{
-	int cp, i, u;
-	int num_clens;
-	int best_clen;
-	int *clens_population;
-	double *clens;
-	double *lambdas;
-	double min_braggp_dist;
-	double clen_to_use;
-	struct rvec cqu;
-	double a, b, c, al, be, ga;
-
-	/* These need to be big enough for the number of different camera
-	 * lengths in the data set.  There are probably only a few, but assume
-	 * the worst case here - a unique camera length for each frame */
-	clens = calloc(n, sizeof(double));
-	clens_population = calloc(n, sizeof(int));
-	lambdas = calloc(n, sizeof(double));
-	if ((lambdas == NULL) || (clens == NULL) || (clens_population == NULL))
-	{
-		ERROR("Failed to allocate memory for clen calculation.\n");
-		free(lambdas);
-		free(clens);
-		free(clens_population);
-		return -1.0;
-	}
-
-	num_clens = 0;
-
-	for ( cp=0; cp<n; cp++ ) {
-
-		int i;
-		int found = 0;
-
-		for ( i=0; i<num_clens; i++ ) {
-			if ( fabs(avg_clen(images[cp]->detgeom) - clens[i]) <0.0001 ) {
-				clens_population[i]++;
-				lambdas[i] += images[cp]->lambda;
-				found = 1;
-				break;
-			}
-		}
-
-		if ( found ) continue;
-
-		clens[num_clens] = avg_clen(images[cp]->detgeom);
-		clens_population[num_clens] = 1;
-		lambdas[num_clens] = images[cp]->lambda;
-		num_clens++;
-
-	}
-
-	for ( u=0; u<num_clens; u++ ) {
-		lambdas[u] /= clens_population[u];
-	}
-
-	if ( num_clens == 1 ) {
-		STATUS("All patterns have the same camera length: %f m.\n",
-		       clens[0]);
-	} else {
-		STATUS("%i different camera lengths were found for the input "
-		       "patterns:\n", num_clens);
-	}
-
-	best_clen = 0;
-	clen_to_use = clens[0];
-	cell_get_parameters(avg, &a, &b, &c, &al, &be, &ga);
-	for ( i=0; i<num_clens; i++ ) {
-
-		assert(clens_population[i] > 0);
-
-		cqu = get_q_from_xyz(1.0/avg_res, 0, clens[i], lambdas[i]);
-
-		min_braggp_dist = fmin(fmin((1.0/cqu.u)/a,
-		                            (1.0/cqu.u)/b),
-		                            (1.0/cqu.u)/c);
-
-		STATUS("Camera length %0.4f m was found %i times.\n"
-		       "Minimum inter-bragg peak distance (based on "
-		       "average cell parameters): %0.1f pixels.\n",
-		       clens[i], clens_population[i], min_braggp_dist);
-
-		if ( min_braggp_dist<1.2*gparams->max_peak_dist ) {
-			STATUS("WARNING: The distance between Bragg peaks is "
-			       "too small: %0.1f < 1.2*%0.1f pixels.\n",
-		               min_braggp_dist, gparams->max_peak_dist);
-		}
-
-		if ( gparams->max_peak_dist==0.0 ) {
-                        gparams->max_peak_dist = 0.5*min_braggp_dist;
-			STATUS("WARNING: Maximum distance between peaks is "
-			       "set to: %0.1f pixels.\n",
-			       gparams->max_peak_dist);
-		}
-
-		if ( clens_population[i] > clens_population[best_clen] ) {
-			best_clen = i;
-			clen_to_use = clens[best_clen];
-		}
-
-	}
-
-	if ( gparams->only_best_distance ) {
-		STATUS("Only %i patterns with camera length %0.4f m will be "
-		       "used.\n", clens_population[best_clen], clen_to_use);
-	}
-
-	free(clens);
-	free(lambdas);
-	free(clens_population);
-
-	return clen_to_use;
-}
-
-
-static double comp_median(double *arr, long n)
-{
-	gsl_sort(arr, 1, n);
-	return gsl_stats_median_from_sorted_data(arr, 1, n);
-}
-
-
-static int panel_in_rigid_group(struct rigid_group *rg,
-                                int panel_number)
-{
-	int i;
-	for ( i=0; i<rg->n_panels; i++ ) {
-		if ( panel_number == rg->panel_numbers[i] ) return 1;
-	}
-	return 0;
-}
-
-
-/* Go up a level in the detector hierarchy, e.g. from "dominoes" to
- * "quadrants", by finding the first panel in "rg" (a rigid group of
- * the finer grained type) in the collection of coarser-grained
- * groups */
-static int find_quad_for_connected(struct rigid_group *rg,
-                                   struct rg_collection *quadrants)
-{
-	int qi;
-
-	for ( qi=0; qi<quadrants->n_rigid_groups; qi++ ) {
-		if ( panel_in_rigid_group(quadrants->rigid_groups[qi],
-		                          rg->panel_numbers[0]) ) {
-			return qi;
-		}
-	}
-
-	/* Hopefully never reached */
-	ERROR("Couldn't find quadrant for connected group!\n");
-	abort();
-}
-
-
-/* Take all the (valid) displacements for pixel "i" in panel "gp", calculate
- * the median displacements in each direction and the modulus */
-static int fill_avg_pixel_displ(struct gpanel *gp, int i)
-{
-	double *list_dx;
-	double *list_dy;
-	int count = 0;
-	int ei;
-
-	list_dx = calloc(gp->num_pix_displ[i], sizeof(double));
-	list_dy = calloc(gp->num_pix_displ[i], sizeof(double));
-	if ( (list_dx == NULL) || (list_dy == NULL) ) {
-		ERROR("Failed to allocate memory for pixel statistics.\n");
-		free(list_dx);
-		free(list_dy);
-		return 1;
-	}
-
-	gp->curr_pix_displ[i] = &gp->pix_displ_list[i];
-
-	for ( ei=0; ei<gp->num_pix_displ[i]; ei++ ) {
-
-		struct single_pixel_displ *pix;
-
-		pix = gp->curr_pix_displ[i];
-
-		if ( pix->dx == -10000.0 ) break;
-		list_dx[count] = pix->dx;
-		list_dy[count] = pix->dy;
-		count++;
-		if ( pix->ne == NULL ) {
-			break;
-		} else {
-			gp->curr_pix_displ[i] = gp->curr_pix_displ[i]->ne;
-		}
-	}
-
-	if ( count < 1 ) {
-		free(list_dx);
-		free(list_dy);
-		return 0;
-	}
-
-	gp->avg_displ_x[i] = comp_median(list_dx, count);
-	gp->avg_displ_y[i] = comp_median(list_dy, count);
-	gp->avg_displ_abs[i] = modulus2d(gp->avg_displ_x[i],
-	                                 gp->avg_displ_y[i]);
-
-	free(list_dx);
-	free(list_dy);
-	return 0;
-}
-
-
-static int allocate_next_element(struct single_pixel_displ **curr_pix_displ,
-                                 int pix_index)
-{
-	curr_pix_displ[pix_index]->ne = malloc(sizeof(struct single_pixel_displ));
-	if ( curr_pix_displ[pix_index]->ne == NULL ) {
-		ERROR("Failed to allocate memory for pixel statistics.\n");
-		return 1;
-	}
-
-	curr_pix_displ[pix_index] = curr_pix_displ[pix_index]->ne;
-
-	return 0;
-}
-
-
-static int add_distance_to_list(struct gpanel *gp,
-				struct imagefeature *imfe,
-                                Reflection *refl, double fx, double fy,
-                                double *det_shift)
-{
-	int pix_index;
-	int ifs, iss;
-	double rfs, rss;
-	double crx, cry;
-
-	ifs = imfe->fs;
-	iss = imfe->ss;  /* Explicit rounding towards zero (truncation) */
-	pix_index = ifs + gp->p->w*iss;
-
-	if ( (ifs >= gp->p->w) || (iss >= gp->p->h) ) {
-		ERROR("Peak is outside panel!\n");
-		return 1;
-	}
-
-	if ( gp->num_pix_displ[pix_index] > 0 ) {
-
-		int ret;
-
-		ret = allocate_next_element(gp->curr_pix_displ, pix_index);
-
-		if ( ret != 0 ) return ret;
-
-	}
-
-	get_detector_pos(refl, &rfs, &rss);
-
-	compute_x_y(rfs, rss, gp->p, &crx, &cry);
-	gp->curr_pix_displ[pix_index]->dx = fx - crx - det_shift[0];
-	gp->curr_pix_displ[pix_index]->dy = fy - cry - det_shift[1];
-	gp->curr_pix_displ[pix_index]->ne = NULL;
-	gp->num_pix_displ[pix_index]++;
-
-	return 0;
-}
-
-
-static int count_pixels_with_min_peaks(struct gpanel *gp, int min_num_peaks,
-                                       int max_num_peaks)
-{
-	int pixel_count = 0;
-	int ifs, iss;
-
-	for ( iss=0; iss<gp->p->h; iss++ ) {
-	for ( ifs=0; ifs<gp->p->w; ifs++ ) {
-
-		int idx = ifs+gp->p->w*iss;
-		if ( gp->num_pix_displ[idx] >= min_num_peaks &&
-		     gp->num_pix_displ[idx] < max_num_peaks) {
-			pixel_count += 1;
-		}
-
-	}
-	}
-	return pixel_count;
-}
-
-
-static void adjust_min_peaks_per_conn(struct rg_collection *connected,
-                                      struct gpanel *gpanels,
-                                      struct geoptimiser_params *gparams,
-                                      struct connected_data *conn_data)
-{
-	int min_num_peaks, di, ip;
-
-	STATUS("Adjusting the minimum number of measurements per pixel in "
-	       "order to have enough measurements for each connected group.\n");
-
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-
-		for ( min_num_peaks=gparams->min_num_peaks_per_pix;
-		      min_num_peaks>0; min_num_peaks-- )
-		{
-			int di_count = 0;
-
-			for ( ip=0; ip<connected->rigid_groups[di]->n_panels;
-			      ip++ )
-			{
-				int pix_count;
-				int pn;
-				struct gpanel *gp;
-
-				pn = connected->rigid_groups[di]->panel_numbers[ip];
-				gp = &gpanels[pn];
-
-				pix_count = count_pixels_with_min_peaks(gp,
-				   min_num_peaks, gparams->max_num_peaks_per_pix);
-
-				di_count += pix_count;
-
-			}
-
-			conn_data[di].n_peaks_in_conn = di_count;
-			if ( di_count >= gparams->min_num_pix_per_conn_group ) {
-				conn_data[di].num_peaks_per_pixel = min_num_peaks;
-				break;
-			}
-		}
-
-		STATUS("Minimum number of measurements per pixel for connected "
-		       "group %s has been set to %i\n", conn_data[di].name,
-		       conn_data[di].num_peaks_per_pixel);
-	}
-}
-
-
-static int compute_pixel_displacements(struct image **images,
-                                       int n_images,
-                                       struct gpanel *gpanels,
-                                       struct rg_collection *connected,
-                                       struct geoptimiser_params *gparams,
-                                       double clen_to_use,
-                                       struct connected_data *conn_data)
-{
-	int cp;
-
-	STATUS("Computing pixel displacements.\n");
-
-	for ( cp=0; cp<n_images; cp++ ) {
-
-		int fi;
-		ImageFeatureList *flist = images[cp]->features;
-		double det_shift[2] = {0.0, 0.0};
-		double shift_x, shift_y;
-		struct detgeom *det = images[cp]->detgeom;
-
-		if ( gparams->only_best_distance ) {
-			if ( fabs(avg_clen(det) - clen_to_use) > 0.0001 ) {
-				continue;
-			}
-		}
-
-		crystal_get_det_shift(images[cp]->crystals[0], &shift_x,
-		                      &shift_y);
-		det_shift[0] = shift_x / det->panels[0].pixel_pitch;
-		det_shift[1] = shift_y / det->panels[0].pixel_pitch;
-
-		for ( fi=0; fi<image_feature_count(images[cp]->features); fi++ ) {
-
-			double min_dist;
-			double fx, fy;
-			Reflection *refl;
-			struct imagefeature *imfe;
-
-			imfe = image_get_feature(flist, fi);
-			if ( imfe == NULL ) continue;
-
-			compute_x_y(imfe->fs, imfe->ss,
-			            &det->panels[imfe->pn], &fx, &fy);
-
-			/* Find the closest reflection (from all crystals) */
-			refl = find_closest_reflection(images[cp], fx, fy,
-			                               &min_dist, det_shift);
-			if ( refl == NULL ) continue;
-
-			if ( min_dist < gparams->max_peak_dist ) {
-
-				struct gpanel *gp;
-				int r;
-				gp = &gpanels[imfe->pn];
-
-				r = add_distance_to_list(gp, imfe, refl, fx, fy,
-				                         det_shift);
-				if ( r ) {
-					ERROR("Error processing peak %f,%f "
-					      "(panel %s), image %s %s\n",
-					      imfe->fs, imfe->ss, gp->p->name,
-					      images[cp]->filename,
-					      images[cp]->ev);
-					return r;
-				}
-
-			}
-		}
-	}
-
-	return 0;
-}
-
-
-static int compute_avg_pix_displ(struct gpanel *gp, int idx,
-                                 int num_peaks_per_pixel,
-                                 int max_peaks_per_pixel)
-{
-	int ret;
-
-	if ( gp->num_pix_displ[idx] >= num_peaks_per_pixel &&
-	     gp->num_pix_displ[idx] < max_peaks_per_pixel ) {
-
-		ret = fill_avg_pixel_displ(gp, idx);
-		if ( ret != 0 ) return ret;
-
-	} else {
-
-		gp->num_pix_displ[idx] = -1;
-		gp->avg_displ_x[idx] = -10000.0;
-		gp->avg_displ_y[idx] = -10000.0;
-		gp->avg_displ_abs[idx] = -10000.0;
-
-	}
-
-	return 0;
-
-}
-
-
-static int compute_avg_displacements(struct rg_collection *connected,
-                                     struct connected_data *conn_data,
-                                     struct gpanel *gpanels,
-                                     struct geoptimiser_params *gparams)
-{
-	int di, ip, ifs, iss;
-	int pix_index, ret, max_peaks;
-	struct detgeom_panel *p;
-
-	max_peaks = 0;
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-
-		for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
-
-			int pn;
-			struct gpanel *gp;
-
-			pn = connected->rigid_groups[di]->panel_numbers[ip];
-			gp = &gpanels[pn];
-
-			for ( iss=0; iss<p->h; iss++ ) {
-			for ( ifs=0; ifs<p->w; ifs++ ) {
-
-				pix_index = ifs+p->w*iss;
-				if ( max_peaks < gp->num_pix_displ[pix_index] )
-					max_peaks = gp->num_pix_displ[pix_index];
-
-				ret = compute_avg_pix_displ(gp, pix_index,
-				            conn_data[di].num_peaks_per_pixel,
-				            gparams->max_num_peaks_per_pix);
-
-				if ( ret != 0 ) return ret;
-
-			}
-			}
-		}
-
-	}
-	STATUS("Maximum peaks per pixel is: %d\n", max_peaks);
-	return 0;
-}
-
-
-static double compute_error(struct rg_collection *connected,
-                            struct connected_data *conn_data,
-                            struct gpanel *gpanels)
-{
-	double total_error = 0;
-	int num_total_error = 0;
-	int di, ip;
-
-	for ( di=0;di<connected->n_rigid_groups;di++ ) {
-
-		double connected_error = 0;
-		int num_connected_error = 0;
-		int ifs, iss;
-
-		for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
-
-			struct gpanel *gp;
-			int pn = connected->rigid_groups[di]->panel_numbers[ip];
-
-			gp = &gpanels[pn];
-
-			for ( iss=0; iss<gp->p->h; iss++ ) {
-			for ( ifs=0; ifs<gp->p->w; ifs++ ) {
-
-				int pix_index;
-				pix_index = ifs+gp->p->w*iss;
-
-				if ( gp->num_pix_displ[pix_index]
-				       >= conn_data[di].num_peaks_per_pixel &&
-				       gp->avg_displ_abs[pix_index] > -9999.0 )
-				{
-					double cer;
-
-					cer = gp->avg_displ_abs[pix_index]
-					    * gp->avg_displ_abs[pix_index];
-					connected_error += cer;
-					num_connected_error++;
-					total_error += cer;
-					num_total_error++;
-				}
-			}
-			}
-
-		}
-
-		if ( num_connected_error > 0 ) {
-
-			connected_error /= (double)num_connected_error;
-			connected_error = sqrt(connected_error);
-
-			STATUS("Error for connected group %s: %d pixels with "
-			       "more than %d peaks: RMSD = %0.4f pixels.\n",
-			       conn_data[di].name, num_connected_error,
-			       conn_data[di].num_peaks_per_pixel,
-			       connected_error);
-		}
-	}
-
-	if ( num_total_error>0 ) {
-		total_error /= (double)num_total_error;
-		total_error = sqrt(total_error);
-	} else {
-		total_error = -1;
-	}
-
-	return total_error;
-}
-
-
-static int compute_rot_stretch_for_empty_panels(struct rg_collection *quads,
-                                                struct rg_collection *conn,
-                                                int min_pix,
-                                                struct connected_data *conn_data)
-{
-	int di,i;
-	double *aver_ang;
-	double *aver_str;
-	int *n;
-
-	STATUS("Computing rotation and elongation corrections for groups "
-	       "without the required number of measurements.\n");
-
-	aver_ang = malloc(quads->n_rigid_groups*sizeof(double));
-	aver_str = malloc(quads->n_rigid_groups*sizeof(double));
-	n = malloc(quads->n_rigid_groups*sizeof(double));
-
-	for ( i=0; i<quads->n_rigid_groups; i++ ) {
-		aver_ang[i] = 0.0;
-		aver_str[i] = 0.0;
-		n[i] = 0;
-	}
-
-	/* Calculate the mean values for the groups which DO have
-	 * enough measurements */
-	for ( di=0; di<conn->n_rigid_groups; di++ ) {
-		if ( conn_data[di].n_peaks_in_conn >= min_pix ) {
-			aver_ang[conn_data[di].num_quad] += conn_data[di].cang;
-			aver_str[conn_data[di].num_quad] += conn_data[di].cstr;
-			n[conn_data[di].num_quad]++;
-		}
-	}
-
-	/* Divide totals to get means */
-	for ( i=0; i<quads->n_rigid_groups; i++ ) {
-		aver_ang[i] /= (double)n[i];
-		aver_str[i] /= (double)n[i];
-	}
-
-	for ( di=0; di<conn->n_rigid_groups; di++ ) {
-
-		int qn = conn_data[di].num_quad;
-		assert(qn < quads->n_rigid_groups);
-
-		if ( conn_data[di].n_peaks_in_conn >= min_pix ) continue;
-
-		if ( n[qn] > 0 ) {
-
-			conn_data[di].cang = aver_ang[qn];
-			conn_data[di].cstr = aver_str[qn];
-			STATUS("Connected group %s has only %i useful "
-			       "pixels. Using average angle: %0.4f "
-			       "degrees\n", conn_data[di].name,
-			       conn_data[di].n_peaks_in_conn,
-			       rad2deg(conn_data[di].cang));
-
-		} else {
-
-			STATUS("Connected group %s does not have enough "
-			       " peaks (%i). It will not be moved.\n",
-			       conn_data[di].name,
-			       conn_data[di].n_peaks_in_conn);
-		}
-	}
-
-	free(aver_ang);
-	free(aver_str);
-	free(n);
-
-	return 0;
-}
-
-
-static DataTemplate *correct_rotation_and_stretch(struct rg_collection *connected,
-                                                  const DataTemplate *dtempl,
-                                                  struct gpanel *gpanels,
-                                                  struct connected_data *conn_data,
-                                                  double clen_to_use,
-                                                  double stretch_coeff,
-                                                  struct geoptimiser_params *gparams)
-{
-#if 0
-	/* FIXME (GitLab #29) */
-	int di, ip;
-
-	STATUS("Applying rotation and stretch corrections.\n");
-
-	if ( gparams->individual_coffset ) {
-
-		STATUS("Using individual distances for rigid panels.\n");
-		for ( di=0; di<connected->n_rigid_groups; di++ ) {
-			for ( ip=0; ip<connected->rigid_groups[di]->n_panels;
-			      ip++ ) {
-
-				struct detgeom_panel *p;
-
-				p = connected->rigid_groups[di]->panels[ip];
-				p->coffset -= (1.0-conn_data[di].cstr)*clen_to_use;
-			}
-		}
-
-	} else {
-
-		int pi;
-
-		for ( di=0; di<connected->n_rigid_groups; di++ ) {
-			conn_data[di].cstr = stretch_coeff;
-		}
-
-		for ( pi=0; pi<det->n_panels; pi++) {
-			det->panels[pi].coffset -= clen_to_use*(1.0-stretch_coeff);
-		}
-		STATUS("Using a single offset distance for the whole "
-		       "detector: %f m. Stretch ceofficient: %0.4f\n",
-		       det->panels[0].coffset, stretch_coeff);
-	}
-
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-
-		int npp = conn_data[di].num_peaks_per_pixel;
-		double cs = conn_data[di].cstr;
-
-		if ( fabs(cs)<FLT_EPSILON ) cs = stretch_coeff;
-
-		for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
-
-			struct detgeom_panel *p;
-			double new_fsx, new_fsy, new_ssx, new_ssy;
-			double new_cnx, new_cny;
-			int fs, ss;
-			struct gpanel *gp;
-
-			p = connected->rigid_groups[di]->panels[ip];
-			gp = &gpanels[panel_number(det, p)];
-
-			new_fsx = p->fsx*cos(conn_data[di].cang)-
-			          p->fsy*sin(conn_data[di].cang);
-			new_fsy = p->fsx*sin(conn_data[di].cang)+
-			          p->fsy*cos(conn_data[di].cang);
-			new_ssx = p->ssx*cos(conn_data[di].cang)-
-			          p->ssy*sin(conn_data[di].cang);
-			new_ssy = p->ssx*sin(conn_data[di].cang)+
-			          p->ssy*cos(conn_data[di].cang);
-
-			/* Calculate corner adjustment
-			 * NB All panels follow the first one */
-			if ( ip > 0 ) {
-
-				struct detgeom_panel *p0;
-				double dx_old, dy_old, dx_new, dy_new;
-
-				p0 = connected->rigid_groups[di]->panels[0];
-
-				dx_old = p->cnx - p0->cnx / cs;
-				dy_old = p->cny - p0->cny / cs;
-
-				dx_new = dx_old*cos(conn_data[di].cang)-
-				         dy_old*sin(conn_data[di].cang);
-				dy_new = dx_old*sin(conn_data[di].cang)+
-				         dy_old*cos(conn_data[di].cang);
-
-				new_cnx = p0->cnx + dx_new;
-				new_cny = p0->cny + dy_new;
-
-			} else {
-
-				new_cnx = p->cnx * cs;
-				new_cny = p->cny * cs;
-			}
-
-			/* The average displacements now need to be updated */
-			for ( ss=0; ss<p->h; ss++ ) {
-			for ( fs=0; fs<p->w; fs++ ) {
-
-				int i = fs + p->w*ss;
-
-				if ( gp->num_pix_displ[i] < npp ) continue;
-
-				gp->avg_displ_x[i] += fs*(new_fsx/cs - p->fsx)
-				                    + ss*(new_ssx/cs - p->ssx)
-				                    + new_cnx/cs - p->cnx;
-
-				gp->avg_displ_y[i] += fs*(new_fsy/cs - p->fsy)
-				                    + ss*(new_ssy/cs - p->ssy)
-				                    + new_cny/cs - p->cny;
-			}
-			}
-
-			p->fsx = new_fsx;
-			p->fsy = new_fsy;
-			p->ssx = new_ssx;
-			p->ssy = new_ssy;
-
-			p->cnx = new_cnx;
-			p->cny = new_cny;
-
-		}
-	}
-#endif
-	return NULL;
-}
-
-
-/* Collect together all the offsets for each group in "connected"
- * Only offsets which have enough peaks per pixel will be used. */
-static int collate_offsets_for_rg(struct rigid_group *group,
-                                  struct gpanel *gpanels,
-                                  int num_peaks_per_pixel,
-                                  double *list_dx, double *list_dy, int list_sz)
-{
-	int ip;
-	int counter = 0;
-
-	for ( ip=0; ip<group->n_panels; ip++ ) {
-
-		int ifs, iss;
-		int pn = group->panel_numbers[ip];
-		struct gpanel *gp = &gpanels[pn];
-
-		for ( iss=0; iss<gp->p->h; iss++ ) {
-		for ( ifs=0; ifs<gp->p->w; ifs++ ) {
-
-			int idx = ifs+gp->p->w*iss;
-
-			if ( gp->num_pix_displ[idx] < num_peaks_per_pixel ) {
-				continue;
-			}
-
-			if ( list_sz == counter ) {
-				ERROR("Too many offsets in group %s!\n",
-				      group->name);
-				return 0;
-			}
-
-			list_dx[counter] = gp->avg_displ_x[idx];
-			list_dy[counter] = gp->avg_displ_y[idx];
-			counter++;
-
-		}
-		}
-	}
-
-	return counter;
-}
-
-
-static void fill_conn_data_sh(struct connected_data *conn_data,
-                              double *list_dx, double *list_dy,
-                              int di, double mpd)
-{
-	conn_data[di].sh_x = comp_median(list_dx,
-	                                 conn_data[di].n_peaks_in_conn);
-	conn_data[di].sh_y = comp_median(list_dy,
-	                                 conn_data[di].n_peaks_in_conn);
-
-	STATUS("Group %s, num pixels: %i, shifts x,y: %0.8f, %0.8f px\n",
-	       conn_data[di].name, conn_data[di].n_peaks_in_conn,
-	       conn_data[di].sh_x, conn_data[di].sh_y);
-
-	if ( modulus2d(conn_data[di].sh_x, conn_data[di].sh_y ) > 0.8*mpd ) {
-		STATUS("WARNING: absolute shift is: %0.1f > 0.8*%0.1f px "
-		       "Increase the value of max_peak_distance!\n",
-		       modulus2d(conn_data[di].sh_x, conn_data[di].sh_y), mpd);
-	}
-}
-
-
-static int compute_shift(struct rg_collection *connected,
-                         struct connected_data *conn_data,
-                         struct geoptimiser_params *gparams,
-                         struct gpanel *gpanels)
-{
-	int di;
-
-        STATUS("Computing shift corrections.\n");
-
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-
-		double *list_dx;
-		double *list_dy;
-		int ct;
-
-		list_dx = malloc(conn_data[di].n_peaks_in_conn*sizeof(double));
-		list_dy = malloc(conn_data[di].n_peaks_in_conn*sizeof(double));
-		if  ( (list_dx == NULL) || (list_dy == NULL) ) {
-			ERROR("Failed to allocate memory for computing shifts\n");
-			free(list_dx);
-			free(list_dy);
-			return 1;
-		}
-
-		ct = collate_offsets_for_rg(connected->rigid_groups[di],
-		                            gpanels,
-		                            conn_data[di].num_peaks_per_pixel,
-		                            list_dx, list_dy,
-		                            conn_data[di].n_peaks_in_conn);
-
-		if ( ct != conn_data[di].n_peaks_in_conn ) {
-			ERROR("Wrong number of peaks for group %s!\n",
-			      connected->rigid_groups[di]->name);
-			ERROR("Counter: %i n_peaks_in_conn: %i\n",
-			       ct, conn_data[di].n_peaks_in_conn);
-			abort();
-		}
-
-		if ( conn_data[di].n_peaks_in_conn
-		          >= gparams->min_num_pix_per_conn_group )
-		{
-			fill_conn_data_sh(conn_data, list_dx, list_dy, di,
-			                  gparams->max_peak_dist);
-
-		} else {
-			conn_data[di].sh_x = -10000.0;
-			conn_data[di].sh_y = -10000.0;
-		}
-
-		free(list_dx);
-		free(list_dy);
-	}
-
-	return 0;
-}
-
-
-static int compute_shift_for_empty_panels(struct rg_collection *quadrants,
-                                            struct rg_collection *connected,
-                                            struct connected_data *conn_data,
-                                            int min_pix)
-{
-	int di, i;
-	double *aver_x;
-	double *aver_y;
-	int *n;
-
-	aver_x = malloc(quadrants->n_rigid_groups * sizeof(double));
-	aver_y = malloc(quadrants->n_rigid_groups * sizeof(double));
-	n = malloc(quadrants->n_rigid_groups * sizeof(int));
-
-	for ( i=0; i<quadrants->n_rigid_groups; i++ ) {
-		aver_x[i] = 0.0;
-		aver_y[i] = 0.0;
-		n[i] = 0;
-	}
-
-	STATUS("Computing shift corrections for groups without the required "
-	       "number of measurements.\n");
-
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-		if ( conn_data[di].n_peaks_in_conn >= min_pix &&
-		     conn_data[di].sh_x > -9999.0 ) {
-			aver_x[conn_data[di].num_quad] += conn_data[di].sh_x;
-			aver_y[conn_data[di].num_quad] += conn_data[di].sh_y;
-			n[conn_data[di].num_quad]++;
-		}
-	}
-
-	for ( i=0; i<quadrants->n_rigid_groups; i++ ) {
-		aver_x[i] /= (double)n[i];
-		aver_y[i] /= (double)n[i];
-	}
-
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-		if ( conn_data[di].n_peaks_in_conn >= min_pix ) continue;
-
-		int qn = conn_data[di].num_quad;
-
-		if ( n[qn] > 0 ) {
-
-			conn_data[di].sh_x = aver_x[qn];
-			conn_data[di].sh_y = aver_y[qn];
-			STATUS("Panel %s doesn't not have enough (%i) "
-			       "peaks. Using average shifts (in pixels) "
-			       "X,Y: %0.2f,%0.2f\n", conn_data[di].name,
-			       conn_data[di].n_peaks_in_conn,
-			       conn_data[di].sh_x, conn_data[di].sh_y);
-
-		} else {
-
-			STATUS("Panel %s has not enough (%i) peaks. "
-			       "It will not be moved.\n",
-			       conn_data[di].name,
-			       conn_data[di].n_peaks_in_conn);
-
-		}
-	}
-
-
-	free(aver_x);
-	free(aver_y);
-	free(n);
-
-	return 0;
-}
-
-
-static void correct_shift(struct rg_collection *connected,
-                          struct connected_data *conn_data,
-                          double clen_to_use)
-{
-#if 0
-	/* FIXME (GitLab #29) */
-	int di;
-	int ip;
-
-	STATUS("Applying shift corrections.\n");
-
-	for ( di=0;di<connected->n_rigid_groups;di++ ) {
-		for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) {
-
-			struct detgeom_panel *p;
-
-			p = connected->rigid_groups[di]->panels[ip];
-
-			if ( conn_data[di].sh_x > -9999.0 &&
-			     conn_data[di].sh_y > -9999.0 ) {
-
-				p->cnx -= conn_data[di].sh_x;
-				p->cny -= conn_data[di].sh_y;
-
-			} else {
-				STATUS("For some reason shifts for panel %s has "
-				       "not been computed!\n", p->name);
-			}
-		}
-	}
-#endif
-}
-
-
-static void scan_p1(int ip0, int ip1, struct rg_collection *connected,
-                    struct connected_data *conn_data,
-                    struct gpanel *gpanels,
-                    int di, double min_dist,
-                    long *num_ang, int ifs0, int iss0,
-                    double c_x0, double c_y0, double cd_x0, double cd_y0,
-                    int compute, double *angles, double *stretches)
-{
-
-	int iss1, ifs1;
-	int p1;
-	struct gpanel *gp1;
-
-	p1 = connected->rigid_groups[di]->panel_numbers[ip1];
-	gp1 = &gpanels[p1];
-
-	int min_ss_p1, min_fs_p1;
-
-	if ( ip0 == ip1 ) {
-		min_fs_p1 = ifs0;
-		min_ss_p1 = iss0;
-	} else {
-		min_fs_p1 = 0;
-		min_ss_p1 = 0;
-	}
-
-	for ( iss1=min_ss_p1; iss1<gp1->p->h; iss1++ ) {
-	for ( ifs1=min_fs_p1; ifs1<gp1->p->w; ifs1++ ) {
-
-		double dist;
-		double c_x1, c_y1, cd_x1, cd_y1;
-		double d_c_x, d_c_y, d_cd_x, d_cd_y;
-		double len1, len2;
-		int pix_index1 = ifs1+gp1->p->w*iss1;
-
-		if ( gp1->num_pix_displ[pix_index1]
-		           < conn_data[di].num_peaks_per_pixel ) continue;
-
-		compute_x_y(ifs1, iss1, gp1->p, &c_x1, &c_y1);
-		cd_x1 = c_x1 - gp1->avg_displ_x[pix_index1];
-		cd_y1 = c_y1 - gp1->avg_displ_y[pix_index1];
-		d_c_x = c_x1-c_x0;
-		d_c_y = c_y1-c_y0;
-		d_cd_x = cd_x1-cd_x0;
-		d_cd_y = cd_y1-cd_y0;
-
-		dist = modulus2d(d_c_x,d_c_y);
-		if ( dist < min_dist ) continue;
-
-		len1 = modulus2d(d_c_x, d_c_y);
-		len2 = modulus2d(d_cd_x, d_cd_y);
-		if ( len1<FLT_EPSILON || len2<FLT_EPSILON ) continue;
-
-		if ( compute ) {
-
-			double scal_m;
-			double multlen;
-
-			scal_m = d_c_x * d_cd_x+ d_c_y * d_cd_y - FLT_EPSILON;
-
-			multlen = len1*len2;
-			if ( fabs(scal_m)>=multlen ) {
-				angles[*num_ang] = 0.0;
-			} else {
-
-				angles[*num_ang] = acos(scal_m/multlen);
-
-				if (d_c_y * d_cd_x - d_c_x * d_cd_y < 0) {
-					angles[*num_ang] *= -1.;
-				}
-
-			}
-
-			stretches[*num_ang] = len1/len2;
-		}
-
-		*num_ang = *num_ang+1;
-	}
-	}
-}
-
-
-/* Executed for each panel in the connected group */
-static void scan_p0(int ip0, struct rg_collection *connected,
-                    struct connected_data *conn_data,
-                    struct gpanel *gpanels,
-                    int di, double min_dist,
-                    long *num_ang, int compute,
-                    double *angles, double *stretches)
-{
-	int iss0, ifs0, ip1;
-	struct gpanel *gp;
-	int p0;
-
-	p0 = connected->rigid_groups[di]->panel_numbers[ip0];
-	gp = &gpanels[p0];
-
-	for ( iss0=0; iss0<gp->p->h; iss0++ ) {
-	for ( ifs0=0; ifs0<gp->p->w; ifs0++ ) {
-
-		double c_x0, c_y0, cd_x0, cd_y0;
-		int pix_index0 = ifs0+gp->p->w*iss0;
-
-		if ( gp->num_pix_displ[pix_index0]
-		                < conn_data[di].num_peaks_per_pixel ) continue;
-
-		compute_x_y(ifs0, iss0, gp->p, &c_x0, &c_y0);
-		cd_x0 = c_x0 - gp->avg_displ_x[pix_index0];
-		cd_y0 = c_y0 - gp->avg_displ_y[pix_index0];
-
-		for ( ip1=ip0; ip1<connected->rigid_groups[di]->n_panels;
-		      ip1++ )
-		{
-			scan_p1(ip0, ip1, connected,
-			        conn_data, gpanels, di, min_dist,
-			        num_ang, ifs0, iss0, c_x0,
-			        c_y0, cd_x0, cd_y0, compute,
-			        angles, stretches);
-		}
-
-	}
-	}
-}
-
-
-static double compute_rotation_and_stretch(struct rg_collection *connected,
-                                           struct connected_data *conn_data,
-                                           struct gpanel *gpanels,
-                                           double dist_coeff_for_rot_str,
-                                           struct geoptimiser_params *gparams)
-{
-	int di;
-	double stretch_cf;
-	long num_coeff;
-	long *num_angles;
-	double *stretch_coeff;
-
-	STATUS("Computing rotation and stretch corrections.\n");
-
-	stretch_coeff = malloc(connected->n_rigid_groups*sizeof(double));
-	num_angles = malloc(connected->n_rigid_groups*sizeof(long));
-	num_coeff = 0;
-
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-
-		long max_num_ang = 0;
-
-		double min_dist;
-		double *angles;
-		double *stretches;
-
-		struct detgeom_panel *first_p;
-		long num_ang = 0;
-		int ip0;
-		int num_pix_first_p;
-
-		/* Enough peaks in this connected group? */
-		if ( conn_data[di].n_peaks_in_conn
-		     < gparams->min_num_pix_per_conn_group ) continue;
-
-		first_p = connected->rigid_groups[di]->panels[0];
-
-		num_pix_first_p = first_p->w * first_p->h;
-
-		/* FIXME: minrad here is not universal (GitLab #29) */
-		min_dist = dist_coeff_for_rot_str * sqrt(num_pix_first_p
-		                       * connected->rigid_groups[di]->n_panels);
-
-		for ( ip0=0; ip0<connected->rigid_groups[di]->n_panels; ip0++ ) {
-			scan_p0(ip0, connected, conn_data, gpanels,
-			        di, min_dist, &num_ang, 0, NULL, NULL);
-		}
-
-		max_num_ang = num_ang+1;
-
-		angles = malloc(max_num_ang*sizeof(double));
-		if ( angles == NULL ) {
-			ERROR("Error allocating memory for angle "
-			      "optimization\n");
-			return -1.0;
-		}
-		stretches = malloc(max_num_ang*sizeof(double));
-		if ( stretches == NULL ) {
-			ERROR("Error allocating memory for stretch "
-			      "optimization\n");
-			return -1.0;
-		}
-
-		num_ang = 0;
-
-		for ( ip0=0; ip0<connected->rigid_groups[di]->n_panels; ip0++ ) {
-			scan_p0(ip0, connected, conn_data, gpanels,
-			        di, min_dist, &num_ang, 1, angles, stretches);
-		}
-
-		if ( num_ang < 1 ) continue;
-
-		conn_data[di].cang = -comp_median(angles, num_ang);
-		conn_data[di].cstr = comp_median(stretches, num_ang);
-
-		STATUS("Panel %s, num: %li, angle: %0.4f deg, stretch coeff: "
-		       "%0.4f\n", conn_data[di].name, num_ang,
-		                  rad2deg(conn_data[di].cang),
-		                  conn_data[di].cstr);
-
-		stretch_coeff[num_coeff] = conn_data[di].cstr;
-		num_angles[num_coeff] = num_ang;
-		num_coeff++;
-
-		free(angles);
-		free(stretches);
-	}
-
-	if ( gparams->norotation ) {
-		STATUS("However, rotation will not be optimized, so the "
-		       "rotation angle has been set to 0\n");
-		for ( di=0; di<connected->n_rigid_groups; di++ ) {
-			conn_data[di].cang = 0.0;
-		}
-	}
-
-	stretch_cf = 1.0;
-
-	printf("Computing overall stretch coefficient.\n");
-
-	if ( num_coeff>0 ) {
-
-		int peaks_per_p;
-
-		peaks_per_p = gparams->min_num_peaks_per_pix;
-
-		while ( peaks_per_p>=0 ) {
-
-			double total_num;
-			long di;
-
-			stretch_cf = 0;
-			total_num = 0;
-			for ( di=0; di<num_coeff; di++ ) {
-				if ( conn_data[di].num_peaks_per_pixel >=
-				     peaks_per_p ) {
-					stretch_cf += stretch_coeff[di]*
-					         (double)num_angles[di];
-					total_num += num_angles[di];
-				}
-			}
-
-			if ( total_num > 0 ) {
-
-				stretch_cf /= total_num;
-
-				printf("(Using only connected groups for which "
-				       "the minimum number of measurements per "
-				       "pixel is %i).\n", peaks_per_p);
-				break;
-			}
-			peaks_per_p--;
-		}
-	}
-
-	if ( stretch_cf < FLT_EPSILON ) {
-		stretch_cf = 1.0;
-	}
-
-	STATUS("The global stretch coefficient for the patterns is %0.4f\n",
-	       stretch_cf);
-	if ( gparams->nostretch ) {
-		stretch_cf = 1.0;
-		STATUS("However, distance offset will not be optimized, so the "
-		       "stretching coefficient has been set to 1.0\n");
-		for ( di=0; di<connected->n_rigid_groups; di++ ) {
-			conn_data[di].cstr = stretch_cf;
-		}
-
-	}
-
-	free(stretch_coeff);
-	free(num_angles);
-
-	return stretch_cf;
-}
-
-
-#ifdef CAN_SAVE_TO_PNG
-
-static void draw_panel(struct image *image, cairo_t *cr,
-                       cairo_matrix_t *basic_m, GdkPixbuf **pixbufs, int i)
-{
-	struct detgeom_panel p = image->det->panels[i];
-	int w = gdk_pixbuf_get_width(pixbufs[i]);
-	int h = gdk_pixbuf_get_height(pixbufs[i]);
-	cairo_matrix_t m;
-
-	/* Start with the basic coordinate system */
-	cairo_set_matrix(cr, basic_m);
-
-	/* Move to the right location */
-	cairo_translate(cr, p.cnx, p.cny);
-
-	/* Twiddle directions according to matrix */
-	cairo_matrix_init(&m, p.fsx, p.fsy, p.ssx, p.ssy, 0.0, 0.0);
-	cairo_transform(cr, &m);
-
-	gdk_cairo_set_source_pixbuf(cr, pixbufs[i], 0.0, 0.0);
-	cairo_rectangle(cr, 0.0, 0.0, w, h);
-}
-
-
-struct rectangle
-{
-	int width, height;
-	double min_x, min_y, max_x, max_y;
-};
-
-
-static int draw_detector(cairo_surface_t *surf, struct image *image,
-                         struct rectangle rect)
-{
-	cairo_t *cr;
-	cairo_matrix_t basic_m;
-	cairo_matrix_t m;
-	GdkPixbuf **pixbufs;
-	int n_pixbufs;
-
-	cr = cairo_create(surf);
-
-	pixbufs = render_panels(image, 1, SCALE_GEOPTIMISER, 1, &n_pixbufs);
-
-	/* Blank grey background */
-	cairo_rectangle(cr, 0.0, 0.0, rect.width, rect.height);
-	cairo_set_source_rgb(cr, 0.5, 0.5, 0.5);
-	cairo_fill(cr);
-
-	/* Set up basic coordinate system
-	 *  - origin in the centre, y upwards. */
-	cairo_identity_matrix(cr);
-	cairo_matrix_init(&m, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0);
-	cairo_translate(cr, -rect.min_x , rect.max_y);
-	cairo_transform(cr, &m);
-	cairo_get_matrix(cr, &basic_m);
-
-	if (pixbufs != NULL) {
-
-		int i;
-
-		for (i = 0; i < image->det->n_panels; i++) {
-			draw_panel(image, cr, &basic_m, pixbufs, i);
-			cairo_fill(cr);
-		}
-
-	}
-
-	/* Free old pixbufs */
-	if (pixbufs != NULL) {
-		int i;
-		for (i = 0; i < n_pixbufs; i++) {
-			g_object_unref(pixbufs[i]);
-		}
-		free(pixbufs);
-	}
-
-	return 0;
-
-}
-
-
-static int save_data_to_png(char *filename, struct detector *det,
-                            struct gpanel *gpanels)
-{
-	struct image im;
-	int i;
-	struct rectangle rect;
-	GdkPixbuf *col_scale;
-	cairo_t *cr;
-
-	cairo_status_t r;
-	cairo_surface_t *surf;
-
-	im.det = det;
-	im.bad = NULL;
-	im.dp = malloc(det->n_panels*sizeof(float *));
-	if ( im.dp == NULL ) {
-		ERROR("Failed to allocate data\n");
-		return 1;
-	}
-	for ( i=0; i<det->n_panels; i++ ) {
-
-		int fs, ss;
-		struct detgeom_panel *p;
-
-		p = &det->panels[i];
-
-		im.dp[i] = calloc(p->w * p->h, sizeof(float));
-		if ( im.dp[i] == NULL ) {
-			ERROR("Failed to allocate data\n");
-			return 1;
-		}
-
-		for ( ss=0; ss<p->h; ss++ ) {
-		for ( fs=0; fs<p->w; fs++ ) {
-
-			int idx;
-			float val;
-
-			idx = fs + ss*p->w;
-
-			if ( gpanels[i].avg_displ_abs[idx] == -10000.0) {
-				val = 0.0;
-			} else if ( gpanels[i].avg_displ_abs[idx] > 1.0) {
-				val = 1.0;
-			} else {
-				val = (float)gpanels[i].avg_displ_abs[idx];
-			}
-			val *= 10.0; /* render_panels sets this as max */
-
-			im.dp[i][fs+p->w*ss] = val;
-
-		}
-		}
-	}
-
-	get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x,
-	                  &rect.max_y);
-
-	if (rect.min_x > 0.0) rect.min_x = 0.0;
-	if (rect.max_x < 0.0) rect.max_x = 0.0;
-	if (rect.min_y > 0.0) rect.min_y = 0.0;
-	if (rect.max_y < 0.0) rect.max_y = 0.0;
-
-	rect.width = rect.max_x - rect.min_x;
-	rect.height = rect.max_y - rect.min_y;
-
-	/* Add a thin border */
-	rect.width += 2.0;
-	rect.height += 2.0;
-	surf = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, rect.width + 20,
-	                                  rect.height);
-
-	draw_detector(surf, &im, rect);
-
-	col_scale = render_get_colour_scale(20, rect.height, SCALE_GEOPTIMISER);
-
-	cr = cairo_create(surf);
-	cairo_identity_matrix(cr);
-	cairo_translate(cr, rect.width, 0.0);
-	cairo_rectangle(cr, 0.0, 0.0, 20.0, rect.height);
-	gdk_cairo_set_source_pixbuf(cr, col_scale, 0.0, 0.0);
-	cairo_fill(cr);
-	cairo_destroy(cr);
-
-	for ( i=0; i<det->n_panels; i++ ) {
-		free(im.dp[i]);
-	}
-	free(im.dp);
-
-	r = cairo_surface_write_to_png(surf, filename);
-	if ( r != CAIRO_STATUS_SUCCESS ) return 1;
-
-	return 0;
-}
-
-
-static int save_stretch_to_png(char *filename, struct detector *det,
-                               struct rg_collection *connected,
-                               struct connected_data *conn_data)
-
-{
-	struct image im;
-	int i, di, ip;
-	float min_str, max_str;
-	struct rectangle rect;
-	GdkPixbuf *col_scale;
-	cairo_t *cr;
-
-	cairo_status_t r;
-	cairo_surface_t *surf;
-
-	im.det = det;
-	im.bad = NULL;
-	im.dp = malloc(det->n_panels*sizeof(float *));
-	if ( im.dp == NULL ) {
-		ERROR("Failed to allocate data\n");
-		return 1;
-	}
-
-	min_str = 10000;
-	max_str = 0;
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-		if (conn_data[di].cstr > max_str) max_str = conn_data[di].cstr;
-		if (conn_data[di].cstr < min_str) min_str = conn_data[di].cstr;
-	}
-
-        i = 0;
-	for ( di=0; di<connected->n_rigid_groups; di++ ) {
-	for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
-
-		int fs, ss;
-		float val;
-		struct detgeom_panel *p;
-
-		p = connected->rigid_groups[di]->panels[ip];
-
-		im.dp[i] = calloc(p->w * p->h, sizeof(float));
-		if ( im.dp[i] == NULL ) {
-			ERROR("Failed to allocate data\n");
-			return 1;
-		}
-
-		val = (conn_data[di].cstr - min_str) / (max_str - min_str);
-		val *= 10.0; /* render_panels sets this as max */
-
-		for ( ss=0; ss<p->h; ss++ ) {
-			for ( fs=0; fs<p->w; fs++ ) {
-
-				im.dp[i][fs+p->w*ss] = val;
-
-			}
-		}
-		i++;
-	}
-	}
-
-	get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x,
-	                  &rect.max_y);
-
-	if (rect.min_x > 0.0) rect.min_x = 0.0;
-	if (rect.max_x < 0.0) rect.max_x = 0.0;
-	if (rect.min_y > 0.0) rect.min_y = 0.0;
-	if (rect.max_y < 0.0) rect.max_y = 0.0;
-
-	rect.width = rect.max_x - rect.min_x;
-	rect.height = rect.max_y - rect.min_y;
-
-	/* Add a thin border */
-	rect.width += 2.0;
-	rect.height += 2.0;
-	surf = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, rect.width + 20,
-	                                  rect.height);
-
-	draw_detector(surf, &im, rect);
-
-	col_scale = render_get_colour_scale(20, rect.height, SCALE_GEOPTIMISER);
-
-	cr = cairo_create(surf);
-	cairo_identity_matrix(cr);
-	cairo_translate(cr, rect.width, 0.0);
-	cairo_rectangle(cr, 0.0, 0.0, 20.0, rect.height);
-	gdk_cairo_set_source_pixbuf(cr, col_scale, 0.0, 0.0);
-	cairo_fill(cr);
-	cairo_destroy(cr);
-
-	for ( i=0; i<det->n_panels; i++ ) {
-		free(im.dp[i]);
-	}
-	free(im.dp);
-
-	r = cairo_surface_write_to_png(surf, filename);
-	if ( r != CAIRO_STATUS_SUCCESS ) return 1;
-
-	return 0;
-}
-
-#else /* CAN_SAVE_TO_PNG */
-
-static int save_stretch_to_png(char *filename, struct detector *det,
-                               struct rg_collection *connected,
-                               struct connected_data *conn_data)
-{
-	ERROR("WARNING: geoptimiser was compiled without gdk-pixbuf and cairo "
-	      "support. Error maps will not be saved.\n");
-	return 0; /* Return "success" so that program continues */
-}
-
-
-static int save_data_to_png(char *filename, struct detector *det,
-                            struct gpanel *gpanels)
-{
-	ERROR("WARNING: geoptimiser was compiled without gdk-pixbuf and cairo "
-	      "support. Error maps will not be saved.\n");
-	return 0; /* Return "success" so that program continues */
-}
-
-#endif /* CAN_SAVE_TO_PNG */
-
-
-int check_and_enforce_cspad_dist(struct geoptimiser_params *gparams,
-				 struct detector *det)
-{
-	int np = 0;
-	int num_errors_found = 0;
-
-	double dist_to_check = 197.0;
-	double tol = 0.2;
-
-	for ( np=0; np<det->n_panels; np = np+2 ) {
-
-		double dist2;
-                double cnx2, cny2;
-
-		struct detgeom_panel *ep = &det->panels[np];
-		struct detgeom_panel *op = &det->panels[np+1];
-
-                cnx2 = ep->cnx + 197.0*ep->fsx;
-                cny2 = ep->cny + 197.0*ep->fsy;
-
-		dist2 = (( cnx2 - op->cnx )*( cnx2 - op->cnx ) +
-		         ( cny2 - op->cny )*( cny2 - op->cny ));
-
-		if ( dist2 > (tol*tol)) {
-
-			num_errors_found += 1;
-
-			STATUS("Warning: distance between panels %s and %s "
-			       "is outside acceptable margins (Corners are "
-                               "more than 0.2 pixels away: %3.2f).\n", ep->name,
-			       op->name, sqrt(dist2));
-
-			if ( gparams->enforce_cspad_layout ) {
-
-				double new_op_cx, new_op_cy;
-
-				new_op_cx = ep->cnx + ep->fsx*dist_to_check;
-				new_op_cy = ep->cny + ep->fsy*dist_to_check;
-
-				op->cnx = new_op_cx;
-				op->cny = new_op_cy;
-
-				STATUS("Enforcing distance....\n");
-			}
-
-		}
-
-		if ( ep->fsx != op->fsx || ep->ssx != op->ssx ||
-		     ep->fsy != op->fsy || ep->ssx != op->ssx ) {
-
-			num_errors_found += 1;
-
-			STATUS("Warning: relative orientation between panels "
-			       "%s and %s is incorrect.\n", ep->name, op->name);
-
-			if ( gparams->enforce_cspad_layout ) {
-
-				STATUS("Enforcing relative orientation....\n");
-
-				op->fsx = ep->fsx;
-				op->ssx = ep->ssx;
-				op->fsy = ep->fsy;
-				op->ssy = ep->ssy;
-
-				op->xfs = ep->xfs;
-				op->xss = ep->xss;
-				op->yfs = ep->yfs;
-				op->yss = ep->yss;
-			}
-
-		}
-
-	}
-	return num_errors_found;
-}
-
-
-int check_and_enforce_agipd_dist(struct geoptimiser_params *gparams,
-				 struct detector *det)
-{
-	int np = 0;
-	int npp = 0;
-	int num_errors_found = 0;
-
-	double dist_to_check = 66.0;
-	double tol = 0.1;
-
-	for ( np=0; np<det->n_panels; np = np+8 ) {
-
-		double dist2;
-                double cnx2, cny2;
-
-		struct detgeom_panel *ep = &det->panels[np];
-
-		for ( npp=0; npp<8; npp++ ) {
-
-			struct detgeom_panel *op = &det->panels[np+npp];
-
-	                cnx2 = ep->cnx + npp*dist_to_check*ep->ssx;
-	                cny2 = ep->cny + npp*dist_to_check*ep->ssy;
-
-			dist2 = (( cnx2 - op->cnx )*( cnx2 - op->cnx ) +
-			         ( cny2 - op->cny )*( cny2 - op->cny ));
-
-			if ( dist2 > (tol*tol)) {
-
-				num_errors_found += 1;
-
-				STATUS("Warning: distance between panels %s and %s "
-				       "is outside acceptable margins (Corners are "
-	                               "more than 0.2 pixels away: %3.2f).\n", ep->name,
-				       op->name, sqrt(dist2));
-
-				if ( gparams->enforce_cspad_layout ) {
-
-					double new_op_cx, new_op_cy;
-
-					new_op_cx = ep->cnx + ep->ssx*npp*dist_to_check;
-					new_op_cy = ep->cny + ep->ssy*npp*dist_to_check;
-
-					op->cnx = new_op_cx;
-					op->cny = new_op_cy;
-
-					STATUS("Enforcing distance....\n");
-				}
-
-			}
-
-			if ( ep->fsx != op->fsx || ep->ssx != op->ssx ||
-			     ep->fsy != op->fsy || ep->ssx != op->ssx ) {
-
-				num_errors_found += 1;
-
-				STATUS("Warning: relative orientation between panels "
-				       "%s and %s is incorrect.\n", ep->name, op->name);
-
-				if ( gparams->enforce_cspad_layout ) {
-
-					STATUS("Enforcing relative orientation....\n");
-
-					op->fsx = ep->fsx;
-					op->ssx = ep->ssx;
-					op->fsy = ep->fsy;
-					op->ssy = ep->ssy;
-
-					op->xfs = ep->xfs;
-					op->xss = ep->xss;
-					op->yfs = ep->yfs;
-					op->yss = ep->yss;
-				}
-
-			}
-		}
-	}
-        return num_errors_found;
-}
-
-
-static struct connected_data *initialize_conn_data(struct rg_collection *quadrants,
-						   struct rg_collection *connected)
-{
-
-	struct connected_data *conn_data;
-
-	conn_data = malloc(connected->n_rigid_groups*
-	                   sizeof(struct connected_data));
-
-	int di;
-
-	for (di=0; di<connected->n_rigid_groups; di++) {
-
-		conn_data[di].num_quad = find_quad_for_connected(
-		                         connected->rigid_groups[di],
-		                         quadrants);
-		conn_data[di].cang = 0.0;
-		conn_data[di].cstr = 1.0;
-		conn_data[di].sh_x = -10000.0;
-		conn_data[di].sh_y = -10000.0;
-		conn_data[di].num_peaks_per_pixel = 1;
-		conn_data[di].name = connected->rigid_groups[di]->name;
-		conn_data[di].n_peaks_in_conn = 0;
-	}
-
-	return conn_data;
-}
-
-
-static int initialize_pixel_displacement_list(struct gpanel *gp)
-{
-	int ipx;
-
-	gp->pix_displ_list = calloc(gp->p->w*gp->p->h,
-	                            sizeof(struct single_pixel_displ));
-	if ( gp->pix_displ_list == NULL ) {
-		ERROR("Error allocating memory for pixel displacement data.\n");
-		return 1;
-	}
-
-	gp->curr_pix_displ = calloc(gp->p->w*gp->p->h,
-	                            sizeof(struct single_pixel_displ *));
-	if ( gp->curr_pix_displ == NULL ) {
-		ERROR("Error allocating memory for pixel displacement data.\n");
-		free(gp->pix_displ_list);
-		return 1;
-	}
-	gp->num_pix_displ = calloc(gp->p->w*gp->p->h, sizeof(int));
-	if ( gp->num_pix_displ == NULL ) {
-		ERROR("Error allocating memory for pixel displacement data.\n");
-		free(gp->pix_displ_list);
-		free(gp->curr_pix_displ);
-		return 1;
-	}
-
-	for ( ipx=0; ipx<gp->p->w*gp->p->h; ipx++ ) {
-		gp->pix_displ_list[ipx].dx = -10000.0;
-		gp->pix_displ_list[ipx].dy = -10000.0;
-		gp->pix_displ_list[ipx].ne = NULL;
-		gp->curr_pix_displ[ipx] = &gp->pix_displ_list[ipx];
-		gp->num_pix_displ[ipx] = 0;
-	}
-
-	return 0;
-}
-
-
-static void free_displ_lists(struct gpanel *gpanels, int n)
-{
-	int j;
-	struct single_pixel_displ *curr = NULL;
-	struct single_pixel_displ *next = NULL;
-
-	for ( j=0; j<n; j++ ) {
-
-		int i;
-		struct gpanel *gp = &gpanels[j];
-
-		for ( i=0; i<gp->p->w*gp->p->h; i++ ) {
-
-			curr = &gp->pix_displ_list[i];
-
-			if ( curr->ne != NULL ) {
-				curr = curr->ne;
-				while ( curr != NULL ) {
-					next = curr->ne;
-					free(curr);
-					curr = next;
-				}
-			}
-		}
-
-		free(gp->curr_pix_displ);
-		free(gp->pix_displ_list);
-
-	}
-}
-
-
-static void recompute_panel_avg_displ(struct rg_collection *connected,
-                                      struct connected_data *conn_data,
-                                      struct detgeom_panel *p,
-                                      struct gpanel *gp,
-                                      int num_peaks_per_pixel,
-                                      double sh_x, double sh_y)
-{
-	int ifs, iss;
-
-	for ( iss=0; iss<p->h; iss++ ) {
-	for ( ifs=0; ifs<p->w; ifs++ ) {
-
-		int pix_index = ifs+p->w*iss;
-
-		if ( gp->num_pix_displ[pix_index] >= num_peaks_per_pixel ) {
-
-			gp->avg_displ_x[pix_index] -= sh_x;
-			gp->avg_displ_y[pix_index] -= sh_y;
-			gp->avg_displ_abs[pix_index] = modulus2d(
-			                            gp->avg_displ_x[pix_index],
-			                            gp->avg_displ_y[pix_index]);
-
-		} else {
-			gp->avg_displ_abs[pix_index] = -10000.0;
-		}
-
-	}
-	}
-}
-
-
-void recompute_avg_displ(struct rg_collection *connected,
-                         struct connected_data *conn_data,
-                         struct detector *det, struct gpanel *gpanels)
-{
-
-	int di, ip;
-
-	for ( di=0;di<connected->n_rigid_groups;di++ ) {
-
-		if (conn_data[di].sh_x < -9999.0) continue;
-
-		for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) {
-
-			struct detgeom_panel *p;
-			struct gpanel *gp;
-
-			p = connected->rigid_groups[di]->panels[ip];
-			gp = &gpanels[panel_number(det, p)];
-			recompute_panel_avg_displ(connected, conn_data, p, gp,
-			                     conn_data[di].num_peaks_per_pixel,
-			                     conn_data[di].sh_x,
-			                     conn_data[di].sh_y);
-
-		}
-	}
-}
-
-
-int optimize_geometry(struct geoptimiser_params *gparams, Stream *st,
-                      struct detector *det,
-                      struct rg_collection *quadrants,
-                      struct rg_collection *connected)
-{
-	int pi;
-	int ret;
-	int write_ret;
-
-	double res_sum;
-	double avg_res;
-	double clen_to_use;
-
-	/* for angles and stretch calculation use
-	 * only pixels which are distco*size_panel
-	 * away */
-	double dist_coeff_for_rot_str = 0.2;
-	double total_error;
-
-	double stretch_coeff = 1.0;
-
-	struct connected_data *conn_data = NULL;
-	struct image **images;
-	int n_images = 0;
-	UnitCell *avg_cell;
-	struct gpanel *gpanels;
-	const char *geometry_data = stream_geometry_file(st);
-
-	STATUS("Maximum distance between peaks: %0.1f pixels.\n",
-	       gparams->max_peak_dist);
-
-	STATUS("Minimum number of measurements for a pixel to be included in "
-	       "the refinement: %i\n", gparams->min_num_peaks_per_pix);
-	STATUS("Minimum number of measurements for connected group for "
-	       "accurate estimation of position/orientation: %i\n",
-	       gparams->min_num_pix_per_conn_group);
-
-	/* CS-PAD */
-	if ( (det->n_panels == 64) && !gparams->no_cspad ) {
-
-		int num_errors = 0;
-
-		STATUS("It looks like the detector is a CSPAD. "
-		       "Checking relative distance and orientation of "
-		       "connected ASICS.\n");
-		STATUS("If the detector is not a CSPAD, please rerun the "
-		       "program with the --no-cspad option.\n");
-
-		STATUS("Enforcing CSPAD layout...\n");
-		num_errors = check_and_enforce_cspad_dist(gparams, det);
-
-		if ( gparams->enforce_cspad_layout ) {
-
-			int geom_wr;
-
-			STATUS("Saving geometry with enforced CSPAD layout.\n"
-			       "Please restart geometry optimization using the "
-			       "optimized geometry from this run as input "
-			       "geometry file.\n");
-			geom_wr = write_detector_geometry_3(geometry_data,
-			                                    gparams->outfile, det,
-			                                    gparams->command_line, 1);
-			if ( geom_wr != 0 ) {
-				ERROR("Error in writing output geometry file.\n");
-				return 1;
-			}
-			STATUS("All done!\n");
-			return 0;
-		}
-
-		if ( !gparams->enforce_cspad_layout && num_errors > 0 ) {
-			ERROR("Relative distance and orientation of connected "
-			      "ASICS do not respect the CSPAD layout.\n"
-			      "Geometry optimization cannot continue.\n"
-			      "Please rerun the program with the "
-			      "--enforce-cspad-layout option.\n");
-			return 1;
-		}
-	}
-
-	/* AGIPD */
-	if ( (det->n_panels == 128) && det->panels[0].res > 4999 &&
-              det->panels[0].res < 5001 && !gparams->no_cspad ) {
-
-		int num_errors = 0;
-
-		STATUS("It looks like the detector is an AGIPD. "
-		       "Checking relative distance and orientation of "
-		       "connected ASICS.\n");
-		STATUS("If the detector is not an AGIPD , please rerun the "
-		       "program with the --no-agipd option.\n");
-
-		STATUS("Enforcing AGIPD layout...\n");
-		num_errors = check_and_enforce_agipd_dist(gparams, det);
-
-		if ( gparams->enforce_cspad_layout ) {
-
-			int geom_wr;
-
-			STATUS("Saving geometry with enforced AGIPD layout.\n"
-			       "Please restart geometry optimization using the "
-			       "optimized geometry from this run as input "
-			       "geometry file.\n");
-			geom_wr = write_detector_geometry_3(geometry_data,
-			                                    gparams->outfile, det,
-			                                    gparams->command_line, 1);
-			if ( geom_wr != 0 ) {
-				ERROR("Error in writing output geometry file.\n");
-				return 1;
-			}
-			STATUS("All done!\n");
-			return 0;
-		}
-
-		if ( !gparams->enforce_cspad_layout && num_errors > 0 ) {
-			ERROR("Relative distance and orientation of connected "
-			      "ASICS do not respect the AGIPD layout.\n"
-			      "Geometry optimization cannot continue.\n"
-			      "Please rerun the program with the "
-			      "--enforce-agipd-layout option.\n");
-			return 1;
-		}
-	}
-
-	images = read_patterns_from_stream(st, dtempl, &n_images);
-	if ( (n_images < 1) || (images == NULL) ) {
-		ERROR("Error reading stream file\n");
-		return 1;
-	}
-
-	avg_cell = compute_avg_cell_parameters(images, n_images);
-	if ( avg_cell == NULL ) {
-		free(images);
-		return 1;
-	}
-
-	res_sum = 0;
-	for ( pi=0; pi<det->n_panels; pi++ ) {
-		res_sum += det->panels[pi].res;
-	}
-	avg_res = res_sum/det->n_panels;
-
-	clen_to_use = pick_clen_to_use(gparams, images, n_images, avg_res,
-	                               avg_cell);
-	if ( clen_to_use < 0.0 ) return 1;
-
-	if ( gparams->max_num_peaks_per_pix == 0 && n_images > 100 ) {
-		gparams->max_num_peaks_per_pix = n_images / 10;
-	} else if ( gparams->max_num_peaks_per_pix == 0 ) {
-		gparams->max_num_peaks_per_pix = n_images;
-	}
-	STATUS("Maximum number of measurements for a pixel to be included in "
-	       "the refinement: %i\n", gparams->max_num_peaks_per_pix);
-
-	gpanels = calloc(det->n_panels, sizeof(struct gpanel));
-	if ( gpanels == NULL ) {
-		ERROR("Failed to allocate panels.\n");
-		return 1;
-	}
-	for ( pi=0; pi<det->n_panels; pi++ ) {
-
-		struct gpanel *gp = &gpanels[pi];
-		int npx = det->panels[pi].w * det->panels[pi].h;
-
-		gp->p = &det->panels[pi];
-
-		gp->avg_displ_x = malloc(npx*sizeof(double));
-		gp->avg_displ_y = malloc(npx*sizeof(double));
-		gp->avg_displ_abs = malloc(npx*sizeof(double));
-
-		if ( (gp->avg_displ_x == NULL) || (gp->avg_displ_y == NULL)
-		  || (gp->avg_displ_abs == NULL) )
-		{
-			ERROR("Failed to allocate displacements.\n");
-			return 1;
-		}
-
-		if ( initialize_pixel_displacement_list(gp) ) {
-			ERROR("Failed to allocate lists.\n");
-			return 1;
-		}
-
-	}
-
-	conn_data = initialize_conn_data(quadrants, connected);
-	if ( conn_data == NULL ) {
-		ERROR("Error allocating memory for connected structure data.\n");
-		return 1;
-	}
-
-	if ( compute_pixel_displacements(images, n_images, gpanels,
-	                                 connected, gparams, clen_to_use,
-	                                 conn_data) ) return 1;
-
-	adjust_min_peaks_per_conn(connected, gpanels, gparams, conn_data);
-
-	if ( compute_avg_displacements(connected, conn_data, gpanels, gparams) ) {
-		free(conn_data);
-		free(images);
-		return 1;
-	}
-
-	free_displ_lists(gpanels, det->n_panels);
-	/* gpanels[].num_pix_displ is still there */
-
-	STATUS("Computing error before correction.\n");
-	total_error = compute_error(connected, det, conn_data, gpanels);
-
-	STATUS("Detector-wide error before correction: RMSD = %0.4f pixels.\n",
-	       total_error);
-
-	if ( gparams->error_maps ) {
-
-		STATUS("Saving error map before correction.\n");
-
-		if ( save_data_to_png("error_map_before.png", det, gpanels) ) {
-			ERROR("Error while writing data to file.\n");
-			free(conn_data);
-			free(images);
-			return 1;
-		}
-
-
-	}
-
-	if ( !gparams->nostretch || !gparams->norotation ) {
-
-		stretch_coeff = compute_rotation_and_stretch(connected,
-	                                             conn_data, det, gpanels,
-                                                     dist_coeff_for_rot_str,
-                                                     gparams);
-		if ( stretch_coeff < 0.0 ) {
-			free(conn_data);
-			return 1;
-		}
-
-		ret = compute_rot_stretch_for_empty_panels(quadrants, connected,
-		                        gparams->min_num_pix_per_conn_group,
-		                        conn_data);
-		if ( ret ) {
-			free(conn_data);
-			return 1;
-		}
-
-		if ( gparams->stretch_map ) {
-
-#ifdef CAN_SAVE_TO_PNG
-
-			STATUS("Saving stretch map - for out-of-plane rotations.\n");
-
-			ret = save_stretch_to_png("stretch_co.png", det,
-			                          connected, conn_data);
-
-			if ( ret ) {
-				ERROR("Error while writing data to file.\n");
-				free(conn_data);
-				return 1;
-			}
-
-#else /* CAN_SAVE_TO_PNG */
-
-			STATUS("ERROR: geoptimiser was compiled without GTK and "
-			       "cairo support.\n Stretch map will not be saved.\n");
-
-#endif /* CAN_SAVE_TO_PNG */
-
-		}
-
-		correct_rotation_and_stretch(connected, det, gpanels,
-		                             conn_data, clen_to_use,
-		                             stretch_coeff, gparams);
-	}
-
-	ret = compute_shift(connected, conn_data, gparams, gpanels);
-	if ( ret != 0 ) {
-		free(conn_data);
-		return 1;
-	}
-
-	compute_shift_for_empty_panels(quadrants, connected, conn_data,
-	                               gparams->min_num_pix_per_conn_group);
-
-	correct_shift(connected, conn_data, clen_to_use);
-
-	recompute_avg_displ(connected, conn_data, det, gpanels);
-
-	if ( gparams->error_maps ) {
-
-#ifdef CAN_SAVE_TO_PNG
-
-		STATUS("Saving error map after correction.\n");
-
-		ret = save_data_to_png("error_map_after.png", det, gpanels);
-		if ( ret ) {
-			ERROR("Error while writing data to file.\n");
-			free(conn_data);
-			return 1;
-		}
-
-#else /* CAN_SAVE_TO_PNG */
-
-		STATUS("geoptimiser was compiled without GTK and Cairo "
-		       "support.\n Error maps will not be saved.\n");
-
-#endif /* CAN_SAVE_TO_PNG */
-
-	}
-
-	STATUS("Computing errors after correction.\n");
-	total_error = compute_error(connected, det, conn_data, gpanels);
-
-	STATUS("Detector-wide error after correction: RMSD = %0.4f pixels.\n",
-	       total_error);
-
-	write_ret = write_detector_geometry_3(geometry_data,
-	                                      gparams->outfile, det,
-	                                      gparams->command_line, 1);
-	if ( write_ret != 0 ) {
-		ERROR("Error in writing output geometry file.\n");
-		return 1;
-	}
-	STATUS("All done!\n");
-	if ( gparams->error_maps ) {
-
-#ifdef CAN_SAVE_TO_PNG
-
-		STATUS("Be sure to inspect error_map_before.png and "
-		       "error_map_after.png !!\n");
-
-#endif /* CAN_SAVE_TO_PNG */
-	}
-
-	free(conn_data);
-	return 0;
-}
-
-
-int main(int argc, char *argv[])
-{
-	int c, i;
-	int ret_val;
-	char buffer[256];
-	char command_line[1024];
-	char *infile = NULL;
-	char *quadrant_coll_name = NULL;
-	char *connected_coll_name = NULL;
-	Stream *st;
-	struct geoptimiser_params *gparams;
-	DataTemplate *dtempl;
-	struct rg_collection *quadrants;
-	struct rg_collection *connected;
-	char *geometry_filename = NULL;
-
-	gparams = malloc(sizeof(struct geoptimiser_params));
-	if ( gparams == NULL ) {
-		ERROR("Couldn't allocate params\n");
-		return 1;
-	}
-
-	gparams->outfile = NULL;
-	gparams->min_num_peaks_per_pix = 3;
-	gparams->max_num_peaks_per_pix = 0;
-	gparams->min_num_pix_per_conn_group = 100;
-	gparams->only_best_distance = 0;
-	gparams->enforce_cspad_layout = 0;
-	gparams->nostretch = 0;
-	gparams->norotation = 0;
-	gparams->individual_coffset = 0;
-	gparams->no_cspad = 0;
-	gparams->error_maps = 1;
-	gparams->stretch_map = 0;
-	gparams->max_peak_dist = 0.0;
-
-	const struct option longopts[] = {
-
-	/* Options with long and short versions */
-	{"help",                             0, NULL,               'h'},
-	{"version",                          0, NULL,               10 },
-	{"input",                            1, NULL,               'i'},
-	{"output",                           1, NULL,               'o'},
-	{"geometry",                         1, NULL,               'g'},
-	{"quadrants",                        1, NULL,               'q'},
-	{"connected",                        1, NULL,               'c'},
-	{"min-num-peaks-per-pixel",          1, NULL,               'x'},
-	{"max-num-peaks-per-pixel",          0, NULL,               'z'},
-	{"min-num-pixels-per-conn-group",    1, NULL,               'p'},
-	{"most-few-clen",                    0, NULL,               'l'},
-	{"max-peak-dist",                    1, NULL,               'm'},
-	{"individual-dist-offset",           0, NULL,               's'},
-
-	/* Long-only options with no arguments */
-	{"no-stretch",                       0, &gparams->nostretch,            1},
-	{"no-rotation",                      0, &gparams->norotation,           1},
-	{"no-error-maps",                    0, &gparams->error_maps,           0},
-	{"stretch-map",                      0, &gparams->stretch_map,          1},
-	{"enforce-cspad-layout",             0, &gparams->enforce_cspad_layout, 1},
-	{"enforce-agipd-layout",             0, &gparams->enforce_cspad_layout, 1},
-	{"no-cspad",                         0, &gparams->no_cspad,             1},
-	{"no-agipd",                         0, &gparams->no_cspad,             1},
-
-	/* Long-only options with arguments */
-	{"min-num-peaks-per-panel",          1, NULL,                          11},
-
-
-	{0, 0, NULL, 0}
-	};
-
-	/* Short options */
-	while ((c = getopt_long(argc, argv, "ho:i:g:q:c:o:x:z:p:lsm:",
-	                       longopts, NULL)) != -1) {
-
-		switch (c) {
-
-			case 'h' :
-			show_help(argv[0]);
-			return 0;
-
-			case 10 :
-			printf("CrystFEL: %s\n",
-			       crystfel_version_string());
-			printf("%s\n",
-			       crystfel_licence_string());
-			return 0;
-
-			case 'o' :
-			gparams->outfile = strdup(optarg);
-			break;
-
-			case 'i' :
-			infile = strdup(optarg);
-			break;
-
-			case 'g' :
-			geometry_filename = strdup(optarg);
-			break;
-
-			case 'q' :
-			quadrant_coll_name = strdup(optarg);
-			break;
-
-			case 'c' :
-			connected_coll_name = strdup(optarg);
-			break;
-
-			case 'x' :
-			gparams->min_num_peaks_per_pix = atoi(optarg);
-			break;
-
-			case 'z' :
-			gparams->max_num_peaks_per_pix = atoi(optarg);
-			break;
-
-			case 'p' :
-			gparams->min_num_pix_per_conn_group = atoi(optarg);
-			break;
-
-			case 'l' :
-			gparams->only_best_distance = 1;
-			break;
-
-			case 'm' :
-			gparams->max_peak_dist = strtof(optarg, NULL);
-			break;
-
-			case 's' :
-			gparams->individual_coffset = 1;
-			break;
-
-			case 11:
-			ERROR("WARNING: The --min-num-peaks-per-panel option has been "
-			      "renamed to --min-num-pixels-per-conn-group. The "
-			      "old option has been deprecated and will "
-			      "soon be removed. It is currently remapped to "
-			      "--min-num-pixels-per-conn-group\n");
-			gparams->min_num_pix_per_conn_group = atoi(optarg);
-			break;
-
-		}
-	}
-
-	if ( infile == NULL ) {
-		ERROR("You must provide an input stream file.\n");
-		return 1;
-	}
-
-	if ( gparams->outfile == NULL ) {
-		ERROR("You must provide an output filename.\n");
-		return 1;
-	}
-
-	if ( quadrant_coll_name == NULL ) {
-		ERROR("You must provide a rigid group collection for "
-		      "quadrants.\n");
-		return 1;
-	}
-
-	if ( connected_coll_name == NULL ) {
-		ERROR("You must provide a rigid group collection for connected "
-		      "panels.\n");
-		return 1;
-	}
-
-	command_line[0] = '\0';
-	for ( i=0; i<argc; i++ ) {
-		if ( i > 0 ) strcat(command_line, " ");
-		strcpy(buffer, argv[i]);
-		strcat(command_line, buffer);
-	}
-
-#ifdef CAN_SAVE_TO_PNG
-#if !GLIB_CHECK_VERSION(2,36,0)
-	g_type_init();
-#endif
-#endif
-
-	st = stream_open_for_read(infile);
-	if ( st == NULL ) {
-		ERROR("Failed to open input stream '%s'\n", infile);
-		return 1;
-	}
-
-	if ( geometry_filename == NULL ) {
-		const char *stgeom = stream_geometry_file(st);
-		if ( stgeom != NULL ) {
-			dtempl = data_template_new_from_string(stgeom);
-		} else {
-			ERROR("No geometry found in stream.  "
-			      "Try again with --geometry=<file>.\n");
-			return 1;
-		}
-	} else {
-		dtempl = data_template_new_from_file(geometry_filename);
-		free(geometry_filename);
-	}
-
-	if ( dtempl == NULL ) {
-		ERROR("Failed to read initial detector geometry.\n");
-		return 1;
-	}
-
-	quadrants = data_template_get_rigid_groups(dtempl,
-	                                           quadrant_coll_name);
-	if ( quadrants == NULL ) {
-		ERROR("Cannot find rigid group collection for quadrants: %s\n",
-		      quadrant_coll_name);
-		return 1;
-	}
-
-	connected = data_template_get_rigid_groups(dtempl,
-	                                           connected_coll_name);
-	if ( connected == NULL ) {
-		ERROR("Cannot find rigid group collection for connected "
-		      "asics: %s\n", connected_coll_name);
-		return 1;
-	}
-
-	ret_val = optimize_geometry(gparams, st, det, quadrants, connected);
-
-	stream_close(st);
-
-	return ret_val;
-}