Move Mille stuff to separate file

5 files changed, 263 insertions, 119 deletions

diff --git a/libcrystfel/meson.build b/libcrystfel/meson.build
index cacf0e82..82937f48 100644
--- a/libcrystfel/meson.build
+++ b/libcrystfel/meson.build
@@ -161,6 +161,10 @@ libcrystfel_sources = ['src/image.c',
                        symop_lex_ch,
                        symop_parse_ch]
 
+if millepededep.found()
+  libcrystfel_sources += 'src/crystfel-mille.c'
+endif
+
 configure_file(input: 'libcrystfel-config.h.meson.in',
                output: 'libcrystfel-config.h',
                configuration: conf_data)
@@ -205,7 +209,8 @@ install_headers(['src/reflist.h',
                  'src/datatemplate.h',
                  'src/colscale.h',
                  'src/detgeom.h',
-                 'src/fom.h'],
+                 'src/fom.h',
+                 'src/crystfel-mille.h'],
                 subdir: 'crystfel')
 
 # API documentation (Doxygen)
diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c
new file mode 100644
index 00000000..c06dff01
--- /dev/null
+++ b/libcrystfel/src/crystfel-mille.c
@@ -0,0 +1,168 @@
+/*
+ * crystfel-mille.c
+ *
+ * Interface to Millepede geometry refinement
+ *
+ * Copyright © 2023 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2023 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/>.
+ *
+ */
+
+#include <libcrystfel-config.h>
+
+#include <stdlib.h>
+#include <assert.h>
+
+#include "image.h"
+#include "geometry.h"
+#include "cell-utils.h"
+#include "predict-refine.h"
+#include "profile.h"
+
+#include <mille_c_wrap.h>
+
+
+static const enum gparam rv[] =
+{
+	GPARAM_ASX,
+	GPARAM_ASY,
+	GPARAM_ASZ,
+	GPARAM_BSX,
+	GPARAM_BSY,
+	GPARAM_BSZ,
+	GPARAM_CSX,
+	GPARAM_CSY,
+	GPARAM_CSZ,
+};
+
+
+int mille_label(int hierarchy_level, int member_index, char param)
+{
+	int label;
+
+	assert(member_index < 1000);
+
+	label = 100000*hierarchy_level + 100*member_index;
+	switch ( param ) {
+		case 'x' : return label+1;
+		case 'y' : return label+2;
+		case 'z' : return label+3;
+		default : abort();
+	}
+}
+
+
+void write_mille(Mille *mille, int n, UnitCell *cell,
+                 struct reflpeak *rps, struct image *image,
+                 double dx, double dy,
+                 const struct detgeom_panel_group *group,
+                 int hierarchy_level, int member_index)
+{
+	int i;
+	float local_gradients[9];
+	float global_gradients[6];
+	int labels[6];
+
+	/* Spot x-position terms */
+	for ( i=0; i<n; i++ ) {
+
+		int j;
+		for ( j=0; j<9; j++ ) {
+			local_gradients[j] = x_gradient(rv[j], rps[i].refl,
+			                                cell, rps[i].panel);
+		}
+
+		global_gradients[0] = x_gradient(GPARAM_DETX, rps[i].refl,
+		                                 cell, rps[i].panel);
+		global_gradients[1] = x_gradient(GPARAM_DETY, rps[i].refl,
+		                                 cell, rps[i].panel);
+		global_gradients[2] = x_gradient(GPARAM_CLEN, rps[i].refl,
+		                                 cell, rps[i].panel);
+		labels[0] = mille_label(hierarchy_level, member_index, 'x');
+		labels[1] = mille_label(hierarchy_level, member_index, 'y');
+		labels[2] = mille_label(hierarchy_level, member_index, 'z');
+
+		mille_add_measurement(mille,
+		                      9, local_gradients,
+		                      3, global_gradients, labels,
+		                      x_dev(&rps[i], image->detgeom, dx, dy),
+		                      0.65*rps[i].panel->pixel_pitch);
+	}
+
+	/* Spot y-position terms */
+	for ( i=0; i<n; i++ ) {
+
+		int j;
+		for ( j=0; j<9; j++ ) {
+			local_gradients[j] = y_gradient(rv[j], rps[i].refl,
+			                                cell, rps[i].panel);
+		}
+
+		global_gradients[0] = y_gradient(GPARAM_DETX, rps[i].refl,
+		                                 cell, rps[i].panel);
+		global_gradients[1] = y_gradient(GPARAM_DETY, rps[i].refl,
+		                                 cell, rps[i].panel);
+		global_gradients[2] = y_gradient(GPARAM_CLEN, rps[i].refl,
+		                                 cell, rps[i].panel);
+		labels[0] = mille_label(hierarchy_level, member_index, 'x');
+		labels[1] = mille_label(hierarchy_level, member_index, 'y');
+		labels[2] = mille_label(hierarchy_level, member_index, 'z');
+
+		mille_add_measurement(mille,
+		                      9, local_gradients,
+		                      3, global_gradients, labels,
+		                      y_dev(&rps[i], image->detgeom, dx, dy),
+		                      0.65*rps[i].panel->pixel_pitch);
+	}
+
+	/* Next level of hierarchy */
+	for ( i=0; i<group->n_children; i++ ) {
+		write_mille(mille, n, cell, rps, image, dx, dy,
+		            group->children[i], hierarchy_level+1, i);
+	}
+
+}
+
+
+Mille *crystfel_mille_new(const char *outFileName,
+                          int asBinary,
+                          int writeZero)
+{
+	return mille_new(outFileName, asBinary, writeZero);
+}
+
+
+void crystfel_mille_free(Mille *m)
+{
+	mille_free(m);
+}
+
+
+void crystfel_mille_delete_last_record(Mille *m)
+{
+	mille_delete_last_record(m);
+}
+
+
+void crystfel_mille_write_record(Mille *m)
+{
+	mille_write_record(m);
+}
diff --git a/libcrystfel/src/crystfel-mille.h b/libcrystfel/src/crystfel-mille.h
new file mode 100644
index 00000000..7ba8de61
--- /dev/null
+++ b/libcrystfel/src/crystfel-mille.h
@@ -0,0 +1,61 @@
+/*
+ * crystfel-mille.h
+ *
+ * Interface to Millepede geometry refinement
+ *
+ * Copyright © 2023 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2023 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/>.
+ *
+ */
+
+#ifndef CRYSTFEL_MILLE_H
+#define CRYSTFEL_MILLE_H
+
+typedef void *Mille;
+
+#include "cell.h"
+#include "image.h"
+#include "predict-refine.h"
+
+/**
+ * \file crystfel-mille.h
+ * Detector geometry refinement using Millepede
+ */
+
+extern Mille *crystfel_mille_new(const char *outFileName,
+                                 int asBinary,
+                                 int writeZero);
+
+extern void crystfel_mille_free(Mille *m);
+
+extern int mille_label(int hierarchy_level, int member_index, char param);
+
+extern void write_mille(Mille *mille, int n, UnitCell *cell,
+                        struct reflpeak *rps, struct image *image,
+                        double dx, double dy,
+                        const struct detgeom_panel_group *group,
+                        int hierarchy_level, int member_index);
+
+extern void crystfel_mille_delete_last_record(Mille *m);
+
+extern void crystfel_mille_write_record(Mille *m);
+
+#endif	/* CRYSTFEL_MILLE_H */
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index 86a37014..585df4f8 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -39,10 +39,7 @@
 #include "cell-utils.h"
 #include "predict-refine.h"
 #include "profile.h"
-
-#ifdef HAVE_MILLEPEDE
-#include <mille_c_wrap.h>
-#endif
+#include "crystfel-mille.h"
 
 
 /** \file predict-refine.h */
@@ -70,15 +67,6 @@ static const enum gparam rv[] =
 	GPARAM_DETY,
 };
 
-struct reflpeak {
-	Reflection *refl;
-	struct imagefeature *peak;
-	double Ih;   /* normalised */
-	struct detgeom_panel *panel;  /* panel the reflection appears on
-                                       * (we assume this never changes) */
-};
-
-
 static void twod_mapping(double fs, double ss, double *px, double *py,
                          struct detgeom_panel *p, double dx, double dy)
 {
@@ -99,8 +87,8 @@ static double r_dev(struct reflpeak *rp)
 }
 
 
-static double x_dev(struct reflpeak *rp, struct detgeom *det,
-                    double dx, double dy)
+double x_dev(struct reflpeak *rp, struct detgeom *det,
+             double dx, double dy)
 {
 	/* Peak position term */
 	double xpk, ypk, xh, yh;
@@ -112,8 +100,8 @@ static double x_dev(struct reflpeak *rp, struct detgeom *det,
 }
 
 
-static double y_dev(struct reflpeak *rp, struct detgeom *det,
-                    double dx, double dy)
+double y_dev(struct reflpeak *rp, struct detgeom *det,
+             double dx, double dy)
 {
 	/* Peak position term */
 	double xpk, ypk, xh, yh;
@@ -573,76 +561,6 @@ static void free_rps_noreflist(struct reflpeak *rps, int n)
 }
 
 
-static void write_mille(Mille *mille, int n, UnitCell *cell,
-                        struct reflpeak *rps, struct image *image,
-                        double dx, double dy)
-{
-#ifdef HAVE_MILLEPEDE
-	int i;
-	float local_gradients[9];
-	float global_gradients[6];
-	int labels[6];
-
-	profile_start("mille-calc");
-
-	/* Spot x-position terms */
-	for ( i=0; i<n; i++ ) {
-
-		int j;
-		for ( j=0; j<9; j++ ) {
-			local_gradients[j] = x_gradient(rv[j], rps[i].refl,
-			                                cell, rps[i].panel);
-		}
-
-		global_gradients[0] = x_gradient(GPARAM_DETX, rps[i].refl,
-		                                 cell, rps[i].panel);
-		global_gradients[1] = x_gradient(GPARAM_DETY, rps[i].refl,
-		                                 cell, rps[i].panel);
-		global_gradients[2] = x_gradient(GPARAM_CLEN, rps[i].refl,
-		                                 cell, rps[i].panel);
-		labels[0] = 1;
-		labels[1] = 2;
-		labels[2] = 3;
-
-		mille_add_measurement(mille,
-		                      9, local_gradients,
-		                      3, global_gradients, labels,
-		                      x_dev(&rps[i], image->detgeom, dx, dy),
-		                      0.65*rps[i].panel->pixel_pitch);
-	}
-
-	/* Spot y-position terms */
-	for ( i=0; i<n; i++ ) {
-
-		int j;
-		for ( j=0; j<9; j++ ) {
-			local_gradients[j] = y_gradient(rv[j], rps[i].refl,
-			                                cell, rps[i].panel);
-		}
-
-		global_gradients[0] = y_gradient(GPARAM_DETX, rps[i].refl,
-		                                 cell, rps[i].panel);
-		global_gradients[1] = y_gradient(GPARAM_DETY, rps[i].refl,
-		                                 cell, rps[i].panel);
-		global_gradients[2] = y_gradient(GPARAM_CLEN, rps[i].refl,
-		                                 cell, rps[i].panel);
-		labels[0] = 1;
-		labels[1] = 2;
-		labels[2] = 3;
-
-		mille_add_measurement(mille,
-		                      9, local_gradients,
-		                      3, global_gradients, labels,
-		                      y_dev(&rps[i], image->detgeom, dx, dy),
-		                      0.65*rps[i].panel->pixel_pitch);
-	}
-
-	profile_end("mille-calc");
-
-#endif /* HAVE_MILLEPEDE */
-}
-
-
 int refine_prediction(struct image *image, Crystal *cr, Mille *mille)
 {
 	int n;
@@ -725,10 +643,15 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille)
 	         pred_residual(rps, n, image->detgeom, total_x, total_y));
 	crystal_add_notes(cr, tmp);
 
+	#ifdef HAVE_MILLEPEDE
 	if ( mille != NULL ) {
+		profile_start("mille-calc");
 		write_mille(mille, n, crystal_get_cell(cr), rps, image,
-		            total_x, total_y);
+		            total_x, total_y,
+		            image->detgeom->top_group, 0, 0);
+		profile_end("mille-calc");
 	}
+	#endif
 
 	crystal_set_det_shift(cr, total_x, total_y);
 
@@ -741,7 +664,7 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille)
 		crystal_set_det_shift(cr, orig_shift_x, orig_shift_y);
 		#ifdef HAVE_MILLEPEDE
 		if ( mille != NULL ) {
-			mille_delete_last_record(mille);
+			crystfel_mille_delete_last_record(mille);
 		}
 		#endif
 		return 1;
@@ -750,30 +673,10 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille)
 	#ifdef HAVE_MILLEPEDE
 	if ( mille != NULL ) {
 		profile_start("mille-write");
-		mille_write_record(mille);
+		crystfel_mille_write_record(mille);
 		profile_end("mille-write");
 	}
 	#endif
 
 	return 0;
 }
-
-
-Mille *crystfel_mille_new(const char *outFileName,
-                          int asBinary,
-                          int writeZero)
-{
-	#ifdef HAVE_MILLEPEDE
-	return mille_new(outFileName, asBinary, writeZero);
-	#else
-	return NULL;
-	#endif
-}
-
-
-void crystfel_mille_free(Mille *m)
-{
-	#ifdef HAVE_MILLEPEDE
-	mille_free(m);
-	#endif
-}
diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h
index 91b44e07..5e2552e6 100644
--- a/libcrystfel/src/predict-refine.h
+++ b/libcrystfel/src/predict-refine.h
@@ -29,24 +29,31 @@
 #ifndef PREDICT_REFINE_H
 #define PREDICT_REFINE_H
 
+struct reflpeak;
 #include "crystal.h"
+#include "crystfel-mille.h"
 
-struct image;
-
-typedef void *Mille;
+struct reflpeak {
+	Reflection *refl;
+	struct imagefeature *peak;
+	double Ih;   /* normalised */
+	struct detgeom_panel *panel;  /* panel the reflection appears on
+                                       * (we assume this never changes) */
+};
 
 /**
  * \file predict-refine.h
  * Prediction refinement: refinement of indexing solutions before integration.
  */
 
-extern Mille *crystfel_mille_new(const char *outFileName,
-                                 int asBinary,
-                                 int writeZero);
-extern void crystfel_mille_free(Mille *m);
-
 extern int refine_prediction(struct image *image, Crystal *cr, Mille *mille);
+
 extern int refine_radius(Crystal *cr, struct image *image);
 
+extern double x_dev(struct reflpeak *rp, struct detgeom *det,
+                    double dx, double dy);
+
+extern double y_dev(struct reflpeak *rp, struct detgeom *det,
+                    double dx, double dy);
 
 #endif	/* PREDICT_REFINE_H */