Add tests/prof2d_check

1 files changed, 192 insertions, 0 deletions

diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
new file mode 100644
index 00000000..24cf4be5
--- /dev/null
+++ b/tests/prof2d_check.c
@@ -0,0 +1,192 @@
+/*
+ * prof2d_check.c
+ *
+ * Check 2D profile fitting
+ *
+ * Copyright © 2013 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 <image.h>
+#include <utils.h>
+#include <beam-parameters.h>
+#include <histogram.h>
+
+#include "../libcrystfel/src/integration.c"
+
+
+#define ADD_PX(fs, ss, val) \
+	if ( ((fs)>0) && ((ss)>0) && ((fs)<w) && ((ss)<h) ) { \
+		image.dp[0][(fs)+w*(ss)] += (val); \
+	}
+
+int main(int argc, char *argv[])
+{
+	struct image image;
+	int fs, ss;
+	FILE *fh;
+	unsigned int seed;
+	int fail = 0;
+	const int w = 1024;
+	const int h = 1024;
+	RefList *list;
+	RefListIterator *iter;
+	Reflection *refl;
+	UnitCell *cell;
+	Crystal *cr;
+	const int ir_inn = 2;
+	const int ir_mid = 4;
+	const int ir_out = 6;
+	Histogram *hi;
+	double esd_sum = 0.0;
+	int n = 0;
+
+	fh = fopen("/dev/urandom", "r");
+	fread(&seed, sizeof(seed), 1, fh);
+	fclose(fh);
+	srand(seed);
+
+	image.flags = NULL;
+	image.beam = NULL;
+	image.lambda = ph_eV_to_lambda(9000.0);
+	image.bw = 0.000001;
+	image.div = 0.0;
+
+	image.det = calloc(1, sizeof(struct detector));
+	image.det->n_panels = 1;
+	image.det->panels = calloc(1, sizeof(struct panel));
+
+	image.dp = calloc(1, sizeof(float *));
+	image.bad = calloc(1, sizeof(int *));
+
+	image.width = w;
+	image.height = h;
+	image.det->panels[0].min_fs = 0;
+	image.det->panels[0].max_fs = w;
+	image.det->panels[0].min_ss = 0;
+	image.det->panels[0].max_ss = h;
+	image.det->panels[0].w = w;
+	image.det->panels[0].h = h;
+	image.det->panels[0].fsx = 1.0;
+	image.det->panels[0].fsy = 0.0;
+	image.det->panels[0].ssx = 0.0;
+	image.det->panels[0].ssy = 1.0;
+	image.det->panels[0].xfs = 1.0;
+	image.det->panels[0].yfs = 0.0;
+	image.det->panels[0].xss = 0.0;
+	image.det->panels[0].yss = 1.0;
+	image.det->panels[0].cnx = -w/2;
+	image.det->panels[0].cny = -h/2;
+	image.det->panels[0].clen = 60.0e-3;
+	image.det->panels[0].res = 100000;  /* 10 px per mm */
+	image.det->panels[0].adu_per_eV = 10.0/9000.0; /* 10 adu/ph */
+	image.det->panels[0].max_adu = +INFINITY;  /* No cutoff */
+
+	image.det->furthest_out_panel = &image.det->panels[0];
+	image.det->furthest_out_fs = 0;
+	image.det->furthest_out_ss = 0;
+
+	image.dp[0] = malloc(w*h*sizeof(float));
+	memset(image.dp[0], 0, w*h*sizeof(float));
+	image.bad[0] = malloc(w*h*sizeof(int));
+	memset(image.bad[0], 0, w*h*sizeof(int));
+
+	cell = cell_new();
+	cell_set_lattice_type(cell, L_CUBIC);
+	cell_set_centering(cell, 'P');
+	cell_set_parameters(cell, 300.0e-10, 300.0e-10, 300.0e-10,
+	                    deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
+	cell = cell_rotate(cell, random_quaternion());
+
+	cr = crystal_new();
+	crystal_set_profile_radius(cr, 0.001e9);
+	crystal_set_mosaicity(cr, 0.0);  /* radians */
+	crystal_set_image(cr, &image);
+	crystal_set_cell(cr, cell);
+
+	image.n_crystals = 1;
+	image.crystals = &cr;
+
+	list = find_intersections(&image, cr);
+
+	for ( fs=0; fs<w; fs++ ) {
+	for ( ss=0; ss<h; ss++ ) {
+		image.dp[0][fs+w*ss] = 10.0*poisson_noise(40);
+	}
+	}
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double pfs, pss;
+		int fs, ss;
+		const int pk_ph = 1000;
+
+		get_detector_pos(refl, &pfs, &pss);
+		fs = pfs;  ss = pss;
+
+		ADD_PX(fs, ss, 10.0*poisson_noise(pk_ph));
+		ADD_PX(fs-1, ss, 10.0*poisson_noise(pk_ph));
+		ADD_PX(fs+1, ss, 10.0*poisson_noise(pk_ph));
+		ADD_PX(fs, ss-1, 10.0*poisson_noise(pk_ph));
+		ADD_PX(fs, ss+1, 10.0*poisson_noise(pk_ph));
+
+	}
+
+	reflist_free(list);  /* integrate_prof2d() will predict again */
+	integrate_prof2d(INTEGRATION_PROF2D, cr, &image, INTDIAG_NONE, 0, 0, 0,
+	                 ir_inn, ir_mid, ir_out);
+
+	hi = histogram_init();
+
+	list = crystal_get_reflections(cr);
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		if ( get_redundancy(refl) == 0 ) continue;
+		histogram_add_value(hi, get_intensity(refl));
+		esd_sum += get_esd_intensity(refl);
+		n++;
+	}
+
+	printf("Mean calculated sigma(I) = %.2f (%i measurements)\n",
+	       esd_sum / n, n);
+
+	histogram_show(hi);
+
+	histogram_free(hi);
+	free(image.beam);
+	free(image.det->panels);
+	free(image.det);
+	free(image.dp[0]);
+	free(image.dp);
+
+	if ( fail ) return 1;
+
+	return 0;
+}