Stuff for partialator

9 files changed, 163 insertions, 102 deletions

diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c
index 8ba955f6..6ff251e0 100644
--- a/libcrystfel/src/crystal.c
+++ b/libcrystfel/src/crystal.c
@@ -49,6 +49,9 @@
 
 struct _crystal
 {
+	/* The image containing the crystal */
+	struct image            *image;
+
 	/* Information about the crystal */
 	UnitCell                *cell;
 	double                  m;     /* Mosaicity in radians */
@@ -60,6 +63,9 @@ struct _crystal
 	/* Integrated (or about-to-be-integrated) reflections */
 	RefList                 *reflections;
 	long long int           n_saturated;  /* Number of overloads */
+
+	/* User flag, e.g. for "this is a bad crystal". */
+	int                     user_flag;
 };
 
 
@@ -139,6 +145,24 @@ long long int crystal_get_num_saturated_reflections(Crystal *cryst)
 }
 
 
+struct image *crystal_get_image(Crystal *cryst)
+{
+	return cryst->image;
+}
+
+
+double crystal_get_osf(Crystal *cryst)
+{
+	return cryst->osf;
+}
+
+
+int crystal_get_user_flag(Crystal *cryst)
+{
+	return cryst->user_flag;
+}
+
+
 /********************************** Setters ***********************************/
 
 
@@ -172,3 +196,21 @@ void crystal_set_num_saturated_reflections(Crystal *cryst, long long int n)
 {
 	cryst->n_saturated = n;
 }
+
+
+void crystal_set_image(Crystal *cryst, struct image *image)
+{
+	cryst->image = image;
+}
+
+
+void crystal_set_osf(Crystal *cryst, double osf)
+{
+	cryst->osf = osf;
+}
+
+
+void crystal_set_user_flag(Crystal *cryst, int user_flag)
+{
+	cryst->user_flag = user_flag;
+}
diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h
index f418f23a..f7c49a81 100644
--- a/libcrystfel/src/crystal.h
+++ b/libcrystfel/src/crystal.h
@@ -57,6 +57,7 @@ extern double crystal_get_resolution_limit(Crystal *cryst);
 extern long long int crystal_get_num_saturated_reflections(Crystal *cryst);
 extern int crystal_get_user_flag(Crystal *cryst);
 extern double crystal_get_osf(Crystal *cryst);
+extern struct image *crystal_get_image(Crystal *cryst);
 
 extern void crystal_set_cell(Crystal *cryst, UnitCell *cell);
 extern void crystal_set_profile_radius(Crystal *cryst, double r);
@@ -66,5 +67,6 @@ extern void crystal_set_num_saturated_reflections(Crystal *cryst,
                                                   long long int n);
 extern void crystal_set_user_flag(Crystal *cryst, int flag);
 extern void crystal_set_osf(Crystal *cryst, double osf);
+extern void crystal_set_image(Crystal *cryst, struct image *image);
 
 #endif	/* CRYSTAL_H */
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 9c494e2c..704baa51 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -346,7 +346,7 @@ RefList *select_intersections(struct image *image, Crystal *cryst)
 
 
 /* Calculate partialities and apply them to the image's reflections */
-void update_partialities(struct image *image, Crystal *cryst)
+void update_partialities(Crystal *cryst)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -354,6 +354,7 @@ void update_partialities(struct image *image, Crystal *cryst)
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
+	struct image *image = crystal_get_image(cryst);
 
 	cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz,
 	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 9d8fa475..aecdc28a 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -46,7 +46,7 @@ extern "C" {
 extern RefList *find_intersections(struct image *image, Crystal *cryst);
 extern RefList *select_intersections(struct image *image, Crystal *cryst);
 
-extern void update_partialities(struct image *image, Crystal *cryst);
+extern void update_partialities(Crystal *cryst);
 
 #ifdef __cplusplus
 }
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index f86ae0bb..0498a532 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -57,7 +57,7 @@
 struct scale_queue_args
 {
 	RefList *reference;
-	struct image *images;
+	Crystal **crystals;
 	int n_started;
 	double max_shift;
 };
@@ -65,7 +65,7 @@ struct scale_queue_args
 
 struct scale_worker_args
 {
-	struct image *image;
+	Crystal *crystal;
 	double shift;
 	RefList *reference;
 };
@@ -79,7 +79,7 @@ static void *create_scale_job(void *vqargs)
 	wargs = malloc(sizeof(struct scale_worker_args));
 	wargs->reference = qargs->reference;
 
-	wargs->image = &qargs->images[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started++];
 
 	return wargs;
 }
@@ -88,20 +88,21 @@ static void *create_scale_job(void *vqargs)
 static void run_scale_job(void *vwargs, int cookie)
 {
 	struct scale_worker_args *wargs = vwargs;
-	struct image *image = wargs->image;
+	Crystal *cr = wargs->crystal;
 	RefList *reference = wargs->reference;
 	Reflection *refl;
 	RefListIterator *iter;
 	double num = 0.0;
 	double den = 0.0;
 	double corr;
+	const double osf = crystal_get_osf(cr);
 
-	if ( image->pr_dud ) {
+	if ( crystal_get_user_flag(cr) ) {
 		wargs->shift = 0.0;
 		return;
 	}
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -128,8 +129,8 @@ static void run_scale_job(void *vwargs, int cookie)
 		Ihl = get_intensity(refl) / get_partiality(refl);
 		esd = get_esd_intensity(refl) / get_partiality(refl);
 
-		num += Ih * (Ihl/image->osf) / pow(esd/image->osf, 2.0);
-		den += pow(Ih, 2.0)/pow(esd/image->osf, 2.0);
+		num += Ih * (Ihl/osf) / pow(esd/osf, 2.0);
+		den += pow(Ih, 2.0)/pow(esd/osf, 2.0);
 
 	}
 
@@ -138,7 +139,7 @@ static void run_scale_job(void *vwargs, int cookie)
 
 	corr = num / den;
 	if ( !isnan(corr) && !isinf(corr) ) {
-		image->osf *= corr;
+		crystal_set_osf(cr, osf*corr);
 	}
 	wargs->shift = fabs(corr-1.0);
 
@@ -155,7 +156,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs)
 }
 
 
-static double iterate_scale(struct image *images, int n, RefList *reference,
+static double iterate_scale(Crystal **crystals, int n, RefList *reference,
                             int n_threads)
 {
 	struct scale_queue_args qargs;
@@ -164,7 +165,7 @@ static double iterate_scale(struct image *images, int n, RefList *reference,
 
 	qargs.reference = reference;
 	qargs.n_started = 0;
-	qargs.images = images;
+	qargs.crystals = crystals;
 	qargs.max_shift = 0.0;
 
 	run_threads(n_threads, run_scale_job, create_scale_job,
@@ -178,14 +179,14 @@ struct merge_queue_args
 {
 	RefList *full;
 	pthread_mutex_t full_lock;
-	struct image *images;
+	Crystal **crystals;
 	int n_started;
 };
 
 
 struct merge_worker_args
 {
-	struct image *image;
+	Crystal *crystal;
 	RefList *full;
 	pthread_mutex_t *full_lock;
 };
@@ -200,7 +201,7 @@ static void *create_merge_job(void *vqargs)
 	wargs->full = qargs->full;
 	wargs->full_lock = &qargs->full_lock;
 
-	wargs->image = &qargs->images[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started++];
 
 	return wargs;
 }
@@ -209,17 +210,17 @@ static void *create_merge_job(void *vqargs)
 static void run_merge_job(void *vwargs, int cookie)
 {
 	struct merge_worker_args *wargs = vwargs;
-	struct image *image = wargs->image;
+	Crystal *cr = wargs->crystal;
 	RefList *full = wargs->full;
 	Reflection *refl;
 	RefListIterator *iter;
 	double G;
 
-	if ( image->pr_dud ) return;
+	if ( crystal_get_user_flag(cr)) return;
 
-	G = image->osf;
+	G = crystal_get_osf(cr);
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -272,7 +273,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 }
 
 
-static RefList *lsq_intensities(struct image *images, int n, int n_threads)
+static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads)
 {
 	RefList *full;
 	struct merge_queue_args qargs;
@@ -283,7 +284,7 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
 
 	qargs.full = full;
 	qargs.n_started = 0;
-	qargs.images = images;
+	qargs.crystals = crystals;
 	pthread_mutex_init(&qargs.full_lock, NULL);
 
 	run_threads(n_threads, run_merge_job, create_merge_job,
@@ -308,14 +309,14 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
 struct esd_queue_args
 {
 	RefList *full;
-	struct image *images;
+	Crystal **crystals;
 	int n_started;
 };
 
 
 struct esd_worker_args
 {
-	struct image *image;
+	Crystal *crystal;
 	RefList *full;
 };
 
@@ -328,7 +329,7 @@ static void *create_esd_job(void *vqargs)
 	wargs = malloc(sizeof(struct esd_worker_args));
 	wargs->full = qargs->full;
 
-	wargs->image = &qargs->images[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started++];
 
 	return wargs;
 }
@@ -337,17 +338,17 @@ static void *create_esd_job(void *vqargs)
 static void run_esd_job(void *vwargs, int cookie)
 {
 	struct esd_worker_args *wargs = vwargs;
-	struct image *image = wargs->image;
+	Crystal *cr = wargs->crystal;
 	RefList *full = wargs->full;
 	Reflection *refl;
 	RefListIterator *iter;
 	double G;
 
-	if ( image->pr_dud ) return;
+	if ( crystal_get_user_flag(cr) ) return;
 
-	G = image->osf;
+	G = crystal_get_osf(cr);
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -383,7 +384,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs)
 }
 
 
-static void calculate_esds(struct image *images, int n, RefList *full,
+static void calculate_esds(Crystal **crystals, int n, RefList *full,
                            int n_threads, int min_red)
 {
 	struct esd_queue_args qargs;
@@ -392,7 +393,7 @@ static void calculate_esds(struct image *images, int n, RefList *full,
 
 	qargs.full = full;
 	qargs.n_started = 0;
-	qargs.images = images;
+	qargs.crystals = crystals;
 
 	for ( refl = first_refl(full, &iter);
 	      refl != NULL;
@@ -423,7 +424,7 @@ static void calculate_esds(struct image *images, int n, RefList *full,
 
 
 /* Scale the stack of images */
-RefList *scale_intensities(struct image *images, int n, RefList *gref,
+RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
                            int n_threads, int noscale)
 {
 	int i;
@@ -431,17 +432,17 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
 	RefList *full = NULL;
 	const int min_redundancy = 3;
 
-	for ( i=0; i<n; i++ ) images[i].osf = 1.0;
+	for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0);
 
 	if ( noscale ) {
-		full = lsq_intensities(images, n, n_threads);
-		calculate_esds(images, n, full, n_threads, min_redundancy);
+		full = lsq_intensities(crystals, n, n_threads);
+		calculate_esds(crystals, n, full, n_threads, min_redundancy);
 		return full;
 	}
 
 	/* No reference -> create an initial list to refine against */
 	if ( gref == NULL ) {
-		full = lsq_intensities(images, n, n_threads);
+		full = lsq_intensities(crystals, n, n_threads);
 	}
 
 	/* Iterate */
@@ -457,14 +458,14 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
 			reference = full;
 		}
 
-		max_corr = iterate_scale(images, n, reference, n_threads);
+		max_corr = iterate_scale(crystals, n, reference, n_threads);
 		//STATUS("Scaling iteration %2i: max correction = %5.2f\n",
 		//       i+1, max_corr);
 
 		/* No reference -> generate list for next iteration */
 		if ( gref == NULL ) {
 			reflist_free(full);
-			full = lsq_intensities(images, n, n_threads);
+			full = lsq_intensities(crystals, n, n_threads);
 		}
 
 		//show_scale_factors(images, n);
@@ -474,10 +475,10 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
 	} while ( (max_corr > 0.01) && (i < MAX_CYCLES) );
 
 	if ( gref != NULL ) {
-		full = lsq_intensities(images, n, n_threads);
+		full = lsq_intensities(crystals, n, n_threads);
 	} /* else we already did it */
 
-	calculate_esds(images, n, full, n_threads, min_redundancy);
+	calculate_esds(crystals, n, full, n_threads, min_redundancy);
 
 	return full;
 }
diff --git a/src/partialator.c b/src/partialator.c
index 35e25320..30f9bfe4 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -506,6 +506,7 @@ int main(int argc, char *argv[])
 			crystals = crystals_new;
 			crystals[n_crystals] = crystal_new();
 			cr = crystals[n_crystals];
+			crystal_set_image(cr, cur);
 
 			/* Fill in initial estimates of stuff */
 			crystal_set_osf(cr, 1.0);
@@ -516,7 +517,7 @@ int main(int argc, char *argv[])
 			as = asymmetric_indices(crystal_get_reflections(cr),
 			                        sym);
 			crystal_set_reflections(cr, as);
-			update_partialities(cur, cr);
+			update_partialities(cr);
 
 			nobs += select_scalable_reflections(as, reference);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 7410d931..31558ec7 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius)
 
 
 /* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, Reflection *refl, double r)
+double gradient(Crystal *cr, int k, Reflection *refl, double r)
 {
 	double ds, azix, aziy;
 	double ttlow, tthigh, tt;
@@ -103,17 +103,18 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
 	int clamp_low, clamp_high;
 	double klow, khigh;
 	double gr;
+	struct image *image = crystal_get_image(cr);
 
 	get_symmetric_indices(refl, &hs, &ks, &ls);
 
-	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
-	                                         &bsx, &bsy, &bsz,
-	                                         &csx, &csy, &csz);
+	cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
+	                                          &bsx, &bsy, &bsz,
+	                                          &csx, &csy, &csz);
 	xl = hs*asx + ks*bsx + ls*csx;
 	yl = hs*asy + ks*bsy + ls*csy;
 	zl = hs*asz + ks*bsz + ls*csz;
 
-	ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls);
+	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
 	get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
 
 	klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
@@ -237,8 +238,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
 
 
 /* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(struct image *image, int k, double shift)
+static void apply_shift(Crystal *cr, int k, double shift)
 {
+	double t;
+	struct image *image = crystal_get_image(cr);
+
 	switch ( k ) {
 
 		case REF_DIV :
@@ -251,7 +255,9 @@ static void apply_shift(struct image *image, int k, double shift)
 		break;
 
 		case REF_R :
-		image->profile_radius += shift;
+		t = crystal_get_profile_radius(cr);
+		t += shift;
+		crystal_set_profile_radius(cr, t);
 		break;
 
 		case REF_ASX :
@@ -263,7 +269,7 @@ static void apply_shift(struct image *image, int k, double shift)
 		case REF_CSX :
 		case REF_CSY :
 		case REF_CSZ :
-		apply_cell_shift(image->indexed_cell, k, shift);
+		apply_cell_shift(crystal_get_cell(cr), k, shift);
 		break;
 
 		default :
@@ -357,7 +363,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
 
 
 /* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(struct image *image, const RefList *full)
+static double pr_iterate(Crystal *cr, const RefList *full)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -369,7 +375,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 	double max_shift;
 	int nref = 0;
 
-	reflections = image->reflections;
+	reflections = crystal_get_reflections(cr);
 
 	M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
 	v = gsl_vector_calloc(NUM_PARAMS);
@@ -387,6 +393,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 		double p;
 		Reflection *match;
 		double gradients[NUM_PARAMS];
+		const double osf = crystal_get_osf(cr);
 
 		if ( !get_refinable(refl) ) continue;
 
@@ -398,7 +405,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 			      " is it marked as refinable?\n", ha, ka, la);
 			continue;
 		}
-		I_full = image->osf * get_intensity(match);
+		I_full = osf * get_intensity(match);
 
 		/* Actual measurement of this reflection from this pattern? */
 		I_partial = get_intensity(refl);
@@ -412,7 +419,8 @@ static double pr_iterate(struct image *image, const RefList *full)
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 			double gr;
-			gr = gradient(image, k, refl, image->profile_radius);
+			gr = gradient(cr, k, refl,
+			              crystal_get_profile_radius(cr));
 			gradients[k] = gr;
 		}
 
@@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 				if ( g > k ) continue;
 
 				M_c = gradients[g] * gradients[k];
-				M_c *= w * pow(image->osf * I_full, 2.0);
+				M_c *= w * pow(osf * I_full, 2.0);
 
 				M_curr = gsl_matrix_get(M, k, g);
 				gsl_matrix_set(M, k, g, M_curr + M_c);
@@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 
 	//STATUS("%i reflections went into the equations.\n", nref);
 	if ( nref == 0 ) {
-		image->pr_dud = 1;
+		crystal_set_user_flag(cr, 1);
 		gsl_matrix_free(M);
 		gsl_vector_free(v);
 		return 0.0;
@@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 
 		for ( param=0; param<NUM_PARAMS; param++ ) {
 			double shift = gsl_vector_get(shifts, param);
-			apply_shift(image, param, shift);
+			apply_shift(cr, param, shift);
 			//STATUS("Shift %i: %e\n", param, shift);
 			if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
 		}
@@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 }
 
 
-static double guide_dev(struct image *image, const RefList *full)
+static double guide_dev(Crystal *cr, const RefList *full)
 {
 	double dev = 0.0;
 
@@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full)
 	Reflection *refl;
 	RefListIterator *iter;
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) ) {
 
@@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full)
 		 * scale_intensities() might not yet have been called, so the
 		 * full version may not have been calculated yet. */
 
-		G = image->osf;
+		G = crystal_get_osf(cr);
 		p = get_partiality(refl);
 		I_partial = get_intensity(refl);
 		I_full = get_intensity(full_version);
@@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full)
 }
 
 
-void pr_refine(struct image *image, const RefList *full)
+void pr_refine(Crystal *cr, const RefList *full)
 {
 	double max_shift, dev;
 	int i;
 	const int verbose = 0;
 
 	if ( verbose ) {
-		dev = guide_dev(image, full);
+		dev = guide_dev(cr, full);
 		STATUS("\n");  /* Deal with progress bar */
 		STATUS("Before iteration:                       dev = %10.5e\n",
 		       dev);
 	}
 
 	i = 0;
-	image->pr_dud = 0;
+	crystal_set_user_flag(cr, 0);
 	do {
 
 		double asx, asy, asz;
@@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full)
 		double csx, csy, csz;
 		double dev;
 
-		cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+		cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
 			               &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
-		max_shift = pr_iterate(image, full);
+		max_shift = pr_iterate(cr, full);
 
-		update_partialities(image);
+		update_partialities(cr);
 
 		if ( verbose ) {
-			dev = guide_dev(image, full);
+			dev = guide_dev(cr, full);
 			STATUS("PR Iteration %2i: max shift = %10.2f"
 			       " dev = %10.5e\n", i+1, max_shift, dev);
 		}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 637c0bb0..52e70b2b 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -62,7 +62,7 @@ enum {
 extern void pr_refine(Crystal *cr, const RefList *full);
 
 /* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(struct image *image, int k, Reflection *refl, double r);
+extern double gradient(Crystal *cr, int k, Reflection *refl, double r);
 
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/scaling-report.c b/src/scaling-report.c
index ba425b96..3de82517 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -183,7 +183,7 @@ static void plot_point(cairo_t *cr, double g_width, double g_height,
 }
 
 
-static void partiality_graph(cairo_t *cr, const struct image *images, int n,
+static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
                              RefList *full)
 {
 	const double g_width = 200.0;
@@ -219,7 +219,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 
 	num_nondud = 0;
 	for ( i=0; i<n; i++ ) {
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
 		num_nondud++;
 	}
 
@@ -229,10 +229,13 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 
 		Reflection *refl;
 		RefListIterator *iter;
+		Crystal *cryst;
 
-		if ( images[i].pr_dud ) continue;
+		cryst = crystals[i];
 
-		for ( refl = first_refl(images[i].reflections, &iter);
+		if ( crystal_get_user_flag(cryst) ) continue;
+
+		for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
 		      refl != NULL;
 		      refl = next_refl(refl, iter) )
 		{
@@ -249,7 +252,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 			if ( get_redundancy(f) < 2 ) continue;
 
 			Ipart = get_intensity(refl);
-			Ifull = images[i].osf * get_intensity(f);
+			Ifull = crystal_get_osf(cryst) * get_intensity(f);
 
 			//if ( Ifull < 10 ) continue;  /* FIXME: Ugh */
 
@@ -267,7 +270,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 					double esd_pobs, esd_Ip, esd_If;
 					esd_Ip = get_esd_intensity(refl);
 					esd_If = get_esd_intensity(f);
-					esd_If *= images[i].osf;
+					esd_If *= crystal_get_osf(cryst);
 					esd_pobs  = pow(esd_Ip/Ipart, 2.0);
 					esd_pobs += pow(esd_If/Ifull, 2.0);
 					esd_pobs = sqrt(esd_pobs);
@@ -313,8 +316,8 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 }
 
 
-static void partiality_histogram(cairo_t *cr, const struct image *images,
-                                 int n, RefList *full, int calc, int backwards)
+static void partiality_histogram(cairo_t *cr, Crystal **crystals, int n,
+                                 RefList *full, int calc, int backwards)
 {
 	int f_max;
 	int i, b;
@@ -344,10 +347,11 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
 
 		Reflection *refl;
 		RefListIterator *iter;
+		Crystal *cryst = crystals[i];
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(cryst) ) continue;
 
-		for ( refl = first_refl(images[i].reflections, &iter);
+		for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
 		      refl != NULL;
 		      refl = next_refl(refl, iter) )
 		{
@@ -364,7 +368,7 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
 			Ipart = get_intensity(refl);
 			Ifull = get_intensity(f);
 
-			pobs = Ipart/(images[i].osf*Ifull);
+			pobs = Ipart/(crystal_get_osf(cryst)*Ifull);
 			pcalc = get_partiality(refl);
 
 			if ( calc ) {
@@ -426,8 +430,8 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
 }
 
 
-static void scale_factor_histogram(cairo_t *cr, const struct image *images,
-                                   int n, const char *title)
+static void scale_factor_histogram(cairo_t *cr, Crystal **crystals, int n,
+                                   const char *title)
 {
 	int f_max;
 	int i, b;
@@ -451,8 +455,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
 
 	osf_max = 0.0;
 	for ( i=0; i<n; i++ ) {
-		double osf = images[i].osf;
-		if ( images[i].pr_dud ) continue;
+		double osf = crystal_get_osf(crystals[i]);
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
 		if ( osf > osf_max ) osf_max = osf;
 	}
 	osf_max = ceil(osf_max+osf_max/10000.0);
@@ -473,9 +477,9 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
 
 		for ( i=0; i<n; i++ ) {
 
-			double osf = images[i].osf;
+			double osf = crystal_get_osf(crystals[i]);
 
-			if ( images[i].pr_dud ) continue;
+			if ( crystal_get_user_flag(crystals[i]) ) continue;
 
 			for ( b=0; b<nbins; b++ ) {
 				if ( (osf >= osf_low[b])
@@ -544,8 +548,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
 }
 
 
-static void intensity_histogram(cairo_t *cr, const struct image *images,
-                                int n, RefList *full,
+static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n,
+                                RefList *full,
                                 signed int h, signed int k, signed int l)
 {
 	int f_max;
@@ -582,12 +586,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
 	for ( i=0; i<n; i++ ) {
 
 		double osf;
+		Crystal *cryst = crystals[i];
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(cryst) ) continue;
 
-		osf = images[i].osf;
+		osf = crystal_get_osf(cryst);
 
-		for ( f = find_refl(images[i].reflections, h, k, l);
+		for ( f = find_refl(crystal_get_reflections(cryst), h, k, l);
 		      f != NULL;
 		      f = next_found_refl(f) )
 		{
@@ -616,12 +621,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
 	for ( i=0; i<n; i++ ) {
 
 		double osf;
+		Crystal *cryst = crystals[i];
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(cryst) ) continue;
 
-		osf = images[i].osf;
+		osf = crystal_get_osf(cryst);
 
-		for ( f = find_refl(images[i].reflections, h, k, l);
+		for ( f = find_refl(crystal_get_reflections(cryst), h, k, l);
 		      f != NULL;
 		      f = next_found_refl(f) )
 		{
@@ -771,7 +777,7 @@ static void find_most_sampled_reflections(RefList *list, int n, signed int *h,
 
 
 
-SRContext *sr_titlepage(struct image *images, int n,
+SRContext *sr_titlepage(Crystal **crystals, int n,
                         const char *filename, const char *stream_filename,
                         const char *cmdline)
 {
@@ -805,7 +811,7 @@ SRContext *sr_titlepage(struct image *images, int n,
 }
 
 
-void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
+void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n,
                   RefList *full)
 {
 	int i;
@@ -822,7 +828,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 
 	cairo_save(sr->cr);
 	cairo_translate(sr->cr, 480.0, 350.0);
-	scale_factor_histogram(sr->cr, images, n,
+	scale_factor_histogram(sr->cr, crystals, n,
 	                       "Distribution of overall scale factors");
 	cairo_restore(sr->cr);
 
@@ -830,7 +836,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 	cairo_save(sr->cr);
 
 	cairo_translate(sr->cr, 70.0, 330.0);
-	partiality_graph(sr->cr, images, n, full);
+	partiality_graph(sr->cr, crystals, n, full);
 
 	cairo_save(sr->cr);
 	cairo_move_to(sr->cr, 0.0, 0.0);
@@ -841,7 +847,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 	cairo_stroke(sr->cr);
 	cairo_set_dash(sr->cr, NULL, 0, 0.0);
 	cairo_translate(sr->cr, 0.0, -150.0);
-	partiality_histogram(sr->cr, images, n, full, 1, 0);
+	partiality_histogram(sr->cr, crystals, n, full, 1, 0);
 	cairo_restore(sr->cr);
 
 	cairo_save(sr->cr);
@@ -854,7 +860,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 	cairo_set_dash(sr->cr, NULL, 0, 0.0);
 	cairo_translate(sr->cr, 230.0, 200.0);
 	cairo_rotate(sr->cr, -M_PI_2);
-	partiality_histogram(sr->cr, images, n, full, 0, 1);
+	partiality_histogram(sr->cr, crystals, n, full, 0, 1);
 	cairo_restore(sr->cr);
 
 	cairo_restore(sr->cr);
@@ -873,7 +879,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 
 		cairo_save(sr->cr);
 		cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y);
-		intensity_histogram(sr->cr, images, n, full,
+		intensity_histogram(sr->cr, crystals, n, full,
 		                    sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]);
 		cairo_restore(sr->cr);