Merge branch 'julia'

17 files changed, 373 insertions, 349 deletions

diff --git a/src/ambigator.c b/src/ambigator.c
index 3777cf0a..91e621ac 100644
--- a/src/ambigator.c
+++ b/src/ambigator.c
@@ -1287,7 +1287,8 @@ int main(int argc, char *argv[])
 			RefList *list;
 			UnitCell *cell;
 
-			cr = image->crystals[i];
+			cr = image->crystals[i].cr;
+			list = image->crystals[i].refls;
 			cell = crystal_get_cell(cr);
 
 			if ( n_crystals == max_crystals ) {
@@ -1308,7 +1309,6 @@ int main(int argc, char *argv[])
 
 			}
 
-			list = crystal_get_reflections(cr);
 			crystals[n_crystals] = asymm_and_merge(list, s_sym,
 			                                       cell,
 			                                       rmin, rmax,
diff --git a/src/cell_explorer.c b/src/cell_explorer.c
index 99cb284f..a92de347 100644
--- a/src/cell_explorer.c
+++ b/src/cell_explorer.c
@@ -2002,9 +2002,7 @@ static int add_stream(CellWindow *w, const char *stream_filename,
 
 		for ( i=0; i<image->n_crystals; i++ ) {
 
-			Crystal *cr;
-
-			cr = image->crystals[i];
+			Crystal *cr = image->crystals[i].cr;
 
 			if ( w->n_cells == max_cells ) {
 
diff --git a/src/crystfelimageview.c b/src/crystfelimageview.c
index 03350204..1659ae4f 100644
--- a/src/crystfelimageview.c
+++ b/src/crystfelimageview.c
@@ -707,9 +707,8 @@ static gint draw_sig(GtkWidget *window, cairo_t *cr, CrystFELImageView *iv)
 	if ( iv->show_refls ) {
 		int i;
 		for ( i=0; i<iv->image->n_crystals; i++ ) {
-			Crystal *cry = iv->image->crystals[i];
 			draw_refls(cr, iv,
-			           crystal_get_reflections(cry),
+			           iv->image->crystals[i].refls,
 			           iv->label_refls,
 			           crystal_cols[i % n_crystal_cols]);
 		}
diff --git a/src/gui_index.c b/src/gui_index.c
index 2cc8e8db..218331eb 100644
--- a/src/gui_index.c
+++ b/src/gui_index.c
@@ -636,19 +636,19 @@ static void run_indexing_once(struct crystfelproject *proj)
 	index_pattern(proj->cur_image, ipriv);
 
 	for ( i=0; i<proj->cur_image->n_crystals; i++ ) {
-		crystal_set_mosaicity(proj->cur_image->crystals[i], 0.0);
+		crystal_set_mosaicity(proj->cur_image->crystals[i].cr, 0.0);
 		if ( proj->indexing_params.use_fix_profile_radius ) {
 
 			/* Manual radius */
-			crystal_set_profile_radius(proj->cur_image->crystals[i],
+			crystal_set_profile_radius(proj->cur_image->crystals[i].cr,
 			                           proj->indexing_params.fix_profile_radius);
 
 		} else {
 
 			/* Auto radius determination */
-			crystal_set_profile_radius(proj->cur_image->crystals[i],
+			crystal_set_profile_radius(proj->cur_image->crystals[i].cr,
 			                           0.02e9);
-			if ( refine_radius(proj->cur_image->crystals[i],
+			if ( refine_radius(proj->cur_image->crystals[i].cr,
 		                           proj->cur_image) )
 			{
 				ERROR("WARNING: Radius determination failed\n");
@@ -682,7 +682,7 @@ static void run_indexing_once(struct crystfelproject *proj)
 	STATUS("Number of crystals: %i\n",
 	       proj->cur_image->n_crystals);
 	for ( i=0; i<proj->cur_image->n_crystals; i++ ) {
-		cell_print(crystal_get_cell(proj->cur_image->crystals[i]));
+		cell_print(crystal_get_cell(proj->cur_image->crystals[i].cr));
 	}
 
 }
diff --git a/src/gui_peaksearch.c b/src/gui_peaksearch.c
index 976e08ec..314fd5d4 100644
--- a/src/gui_peaksearch.c
+++ b/src/gui_peaksearch.c
@@ -41,6 +41,7 @@
 
 #include <datatemplate.h>
 #include <peaks.h>
+#include <peakfinder8.h>
 
 #include "crystfelimageview.h"
 #include "gui_project.h"
@@ -60,52 +61,55 @@ void update_peaks(struct crystfelproject *proj)
 
 	switch ( proj->peak_search_params.method ) {
 
+		ImageFeatureList *peaks;
+
 		case PEAK_ZAEF:
-		search_peaks(proj->cur_image,
-		             proj->peak_search_params.threshold,
-		             proj->peak_search_params.min_sq_gradient,
-		             proj->peak_search_params.min_snr,
-		             proj->peak_search_params.pk_inn,
-		             proj->peak_search_params.pk_mid,
-		             proj->peak_search_params.pk_out,
-		             1);
+		proj->cur_image->features = search_peaks(proj->cur_image,
+		                                         proj->peak_search_params.threshold,
+		                                         proj->peak_search_params.min_sq_gradient,
+		                                         proj->peak_search_params.min_snr,
+		                                         proj->peak_search_params.pk_inn,
+		                                         proj->peak_search_params.pk_mid,
+		                                         proj->peak_search_params.pk_out,
+		                                         1);
 		break;
 
 		case PEAK_PEAKFINDER8:
-		search_peaks_peakfinder8(proj->cur_image, 2048,
-		                         proj->peak_search_params.threshold,
-		                         proj->peak_search_params.min_snr,
-		                         proj->peak_search_params.min_pix_count,
-		                         proj->peak_search_params.max_pix_count,
-		                         proj->peak_search_params.local_bg_radius,
-		                         proj->peak_search_params.min_res,
-		                         proj->peak_search_params.max_res,
-		                         1, 0, NULL);
+		proj->cur_image->features = peakfinder8(proj->cur_image, 2048,
+		                                        proj->peak_search_params.threshold,
+		                                        proj->peak_search_params.min_snr,
+		                                        proj->peak_search_params.min_pix_count,
+		                                        proj->peak_search_params.max_pix_count,
+		                                        proj->peak_search_params.local_bg_radius,
+		                                        proj->peak_search_params.min_res,
+		                                        proj->peak_search_params.max_res,
+		                                        1, 0, NULL);
 		break;
 
 		case PEAK_PEAKFINDER9:
-		search_peaks_peakfinder9(proj->cur_image,
-		                         proj->peak_search_params.min_snr_biggest_pix,
-		                         proj->peak_search_params.min_snr_peak_pix,
-		                         proj->peak_search_params.min_snr,
-		                         proj->peak_search_params.min_sig,
-		                         proj->peak_search_params.min_peak_over_neighbour,
-		                         proj->peak_search_params.local_bg_radius);
+		proj->cur_image->features = search_peaks_peakfinder9(proj->cur_image,
+		                                                     proj->peak_search_params.min_snr_biggest_pix,
+		                                                     proj->peak_search_params.min_snr_peak_pix,
+		                                                     proj->peak_search_params.min_snr,
+		                                                     proj->peak_search_params.min_sig,
+		                                                     proj->peak_search_params.min_peak_over_neighbour,
+		                                                     proj->peak_search_params.local_bg_radius);
 		break;
 
 		case PEAK_HDF5:
 		case PEAK_CXI:
-		proj->cur_image->features = image_read_peaks(proj->dtempl,
-		                                             proj->cur_image->filename,
-		                                             proj->cur_image->ev,
-		                                             proj->peak_search_params.half_pixel_shift);
+		peaks = image_read_peaks(proj->dtempl,
+		                         proj->cur_image->filename,
+		                         proj->cur_image->ev,
+		                         proj->peak_search_params.half_pixel_shift);
 		if ( proj->peak_search_params.revalidate ) {
-			validate_peaks(proj->cur_image,
-			               proj->peak_search_params.min_snr,
-			               proj->peak_search_params.pk_inn,
-			               proj->peak_search_params.pk_mid,
-			               proj->peak_search_params.pk_out,
-			               1, 0);
+			proj->cur_image->features = validate_peaks(proj->cur_image, peaks,
+			                                           proj->peak_search_params.min_snr,
+			                                           proj->peak_search_params.pk_inn,
+			                                           proj->peak_search_params.pk_mid,
+			                                           proj->peak_search_params.pk_out,
+			                                           1, 0);
+			image_feature_list_free(peaks);
 		}
 		break;
 
diff --git a/src/merge.c b/src/merge.c
index 0a5c5f14..4a2dae66 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -55,7 +55,7 @@ struct merge_queue_args
 {
 	RefList *full;
 	pthread_rwlock_t full_lock;
-	Crystal **crystals;
+	struct crystal_refls *crystals;
 	int n_started;
 	double push_res;
 	int use_weak;
@@ -68,6 +68,7 @@ struct merge_worker_args
 {
 	struct merge_queue_args *qargs;
 	Crystal *crystal;
+	RefList *refls;
 	int crystal_number;
 	int n_reflections;
 };
@@ -81,7 +82,9 @@ static void *create_merge_job(void *vqargs)
 	wargs = malloc(sizeof(struct merge_worker_args));
 	wargs->qargs = qargs;
 	wargs->crystal_number = qargs->n_started;
-	wargs->crystal = qargs->crystals[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started].cr;
+	wargs->refls = qargs->crystals[qargs->n_started].refls;
+	qargs->n_started++;
 
 	return wargs;
 }
@@ -180,7 +183,7 @@ static void run_merge_job(void *vwargs, int cookie)
 	G = crystal_get_osf(cr);
 	B = crystal_get_Bfac(cr);
 
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	for ( refl = first_refl(wargs->refls, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -262,8 +265,8 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 }
 
 
-RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
-                           int min_meas,
+RefList *merge_intensities(struct crystal_refls *crystals, int n,
+                           int n_threads, int min_meas,
                            double push_res, int use_weak, int ln_merge)
 {
 	RefList *full;
@@ -364,7 +367,7 @@ double correct_reflection(double val, Reflection *refl, double osf, double Bfac,
 }
 
 
-double residual(Crystal *cr, const RefList *full, int free,
+double residual(RefList *list, Crystal *cr, const RefList *full, int free,
                 int *pn_used, const char *filename)
 {
 	Reflection *refl;
@@ -376,7 +379,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 	double B = crystal_get_Bfac(cr);
 	UnitCell *cell = crystal_get_cell(cr);
 
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -420,7 +423,8 @@ double residual(Crystal *cr, const RefList *full, int free,
 }
 
 
-double log_residual(Crystal *cr, const RefList *full, int free,
+double log_residual(RefList *list, Crystal *cr,
+                    const RefList *full, int free,
                     int *pn_used, const char *filename)
 {
 	double dev = 0.0;
@@ -439,7 +443,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 		}
 	}
 
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -487,7 +491,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 
 
 /* Has to match run_merge_job to be useful */
-void write_unmerged(const char *fn, Crystal **crystals, int n_crystals)
+void write_unmerged(const char *fn, struct crystal_refls *crystals, int n_crystals)
 {
 	FILE *fh;
 	int i;
@@ -506,17 +510,17 @@ void write_unmerged(const char *fn, Crystal **crystals, int n_crystals)
 		UnitCell *cell;
 
 		fprintf(fh, "Crystal %i\n", i);
-		if ( crystal_get_user_flag(crystals[i]) ) {
+		if ( crystal_get_user_flag(crystals[i].cr) ) {
 			fprintf(fh, "Flagged: yes\n");
 		} else {
 			fprintf(fh, "Flagged: no\n");
 		}
 
-		G = crystal_get_osf(crystals[i]);
-		B = crystal_get_Bfac(crystals[i]);
-		cell = crystal_get_cell(crystals[i]);
+		G = crystal_get_osf(crystals[i].cr);
+		B = crystal_get_Bfac(crystals[i].cr);
+		cell = crystal_get_cell(crystals[i].cr);
 
-		for ( refl = first_refl(crystal_get_reflections(crystals[i]), &iter);
+		for ( refl = first_refl(crystals[i].refls, &iter);
 		      refl != NULL;
 		      refl = next_refl(refl, iter) )
 		{
diff --git a/src/merge.h b/src/merge.h
index 5f7e0e64..e304db78 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -43,7 +43,7 @@
 #define MIN_PART_MERGE (0.3)
 
 
-extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
+extern RefList *merge_intensities(struct crystal_refls *crystals, int n, int n_threads,
                                   int min_meas, double push_res, int use_weak,
                                   int ln_merge);
 
@@ -53,12 +53,12 @@ extern double correct_reflection_nopart(double val, Reflection *refl,
 extern double correct_reflection(double val, Reflection *refl, double osf,
                                  double Bfac, double res);
 
-extern double residual(Crystal *cr, const RefList *full, int free,
+extern double residual(RefList *list, Crystal *cr, const RefList *full, int free,
                        int *pn_used, const char *filename);
 
-extern double log_residual(Crystal *cr, const RefList *full, int free,
+extern double log_residual(RefList *list, Crystal *cr, const RefList *full, int free,
                            int *pn_used, const char *filename);
 
-extern void write_unmerged(const char *fn, Crystal **crystals, int n_crystals);
+extern void write_unmerged(const char *fn, struct crystal_refls *crystals, int n_crystals);
 
 #endif	/* MERGE */
diff --git a/src/partialator.c b/src/partialator.c
index 9b1803b0..7d3f19f6 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -169,14 +169,14 @@ static void add_to_csplit(struct custom_split *csplit, const char *id,
 
 
 /* Write two-way split results (i.e. for CC1/2 etc) for this list of crystals */
-static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
-                        int nthreads, PartialityModel pmodel,
+static void write_split(struct crystal_refls *crystals, int n_crystals,
+                        const char *outfile, int nthreads, PartialityModel pmodel,
                         int min_measurements, SymOpList *sym, double push_res)
 {
 	char tmp[1024];
 	RefList *split;
-	Crystal **crystals1;
-	Crystal **crystals2;
+	struct crystal_refls *crystals1;
+	struct crystal_refls *crystals2;
 	int n_crystals1 = 0;
 	int n_crystals2 = 0;
 	int i;
@@ -186,13 +186,12 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 		return;
 	}
 
-	crystals1 = malloc(n_crystals * sizeof(Crystal *));
+	crystals1 = malloc(n_crystals * sizeof(struct crystal_refls));
 	if ( crystals1 == NULL ) return;
 
-	crystals2 = malloc(n_crystals * sizeof(Crystal *));
+	crystals2 = malloc(n_crystals * sizeof(struct crystal_refls));
 	if ( crystals2 == NULL ) return;
 
-
 	for ( i=0; i<n_crystals; i++ ) {
 		if ( i % 2 ) {
 			crystals1[n_crystals1] = crystals[i];
@@ -259,14 +258,15 @@ static char *insert_into_filename(const char *fn, const char *add)
 
 /* Write custom split results (including a two-way split) */
 static void write_custom_split(struct custom_split *csplit, int dsn,
-                               Crystal **crystals, int n_crystals,
+                               struct crystal_refls *crystals,
+                               struct image **images, int n_crystals,
                                PartialityModel pmodel, int min_measurements,
                                double push_res, SymOpList *sym, int nthreads,
                                const char *outfile)
 {
 	char *tmp;
 	RefList *split;
-	Crystal *crystalsn[n_crystals];
+	struct crystal_refls crystalsn[n_crystals];
 	int n_crystalsn = 0;
 	int i;
 
@@ -277,8 +277,8 @@ static void write_custom_split(struct custom_split *csplit, int dsn,
 		char *id;
 		int dsn_crystal;
 
-		fn = crystal_get_image(crystals[i])->filename;
-		evs = crystal_get_image(crystals[i])->ev;
+		fn = images[i]->filename;
+		evs = images[i]->ev;
 		if ( evs == NULL ) evs = "//";
 
 		id = malloc(strlen(evs)+strlen(fn)+2);
@@ -362,7 +362,7 @@ static void show_help(const char *s)
 }
 
 
-static signed int find_first_crystal(Crystal **crystals, int n_crystals,
+static signed int find_first_crystal(struct image **images, int n_crystals,
                                      struct custom_split *csplit, int dsn)
 {
 	int i;
@@ -373,8 +373,8 @@ static signed int find_first_crystal(Crystal **crystals, int n_crystals,
 		char *id;
 		int dsn_crystal;
 
-		fn = crystal_get_image(crystals[i])->filename;
-		evs = crystal_get_image(crystals[i])->ev;
+		fn = images[i]->filename;
+		evs = images[i]->ev;
 		if ( evs == NULL ) evs = "//";
 
 		id = malloc(strlen(evs)+strlen(fn)+2);
@@ -394,8 +394,8 @@ static signed int find_first_crystal(Crystal **crystals, int n_crystals,
 }
 
 
-static void check_csplit(Crystal **crystals, int n_crystals,
-                         struct custom_split *csplit)
+static void check_csplit(struct crystal_refls *crystals, struct image **images,
+                         int n_crystals, struct custom_split *csplit)
 {
 	int i;
 	int n_nosplit = 0;
@@ -412,8 +412,8 @@ static void check_csplit(Crystal **crystals, int n_crystals,
 		char *id;
 		int dsn_crystal;
 
-		fn = crystal_get_image(crystals[i])->filename;
-		evs = crystal_get_image(crystals[i])->ev;
+		fn = images[i]->filename;
+		evs = images[i]->ev;
 		if ( evs == NULL ) evs = "//";
 
 		id = malloc(strlen(evs)+strlen(fn)+2);
@@ -437,7 +437,7 @@ static void check_csplit(Crystal **crystals, int n_crystals,
 	for ( i=0; i<csplit->n_datasets; i++ ) {
 
 		/* Try to find a crystal with dsn = i */
-		if ( find_first_crystal(crystals, n_crystals, csplit, i) != -1 )
+		if ( find_first_crystal(images, n_crystals, csplit, i) != -1 )
 		{
 			n_cry++;
 		} else {
@@ -621,11 +621,10 @@ static void skip_to_end(FILE *fh)
 }
 
 
-static int set_initial_params(Crystal *cr, FILE *fh, double force_bandwidth,
-                              double force_radius, double force_lambda)
+static int set_initial_params(Crystal *cr, struct image *image, FILE *fh,
+                              double force_bandwidth, double force_radius,
+                              double force_lambda)
 {
-	struct image *image = crystal_get_image(cr);
-
 	if ( fh != NULL ) {
 
 		int err;
@@ -736,7 +735,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 }
 
 
-static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
+static void write_pgraph(RefList *full, struct crystal_refls *crystals, int n_crystals,
                          signed int iter, const char *suff,
                          const char *log_folder)
 {
@@ -762,16 +761,16 @@ static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
 	}
 
 	for ( i=0; i<n_crystals; i++ ) {
-		if ( crystal_get_user_flag(crystals[i]) != 0 ) continue;
-		write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full,
-		                crystals[i], i, iter);
+		if ( crystal_get_user_flag(crystals[i].cr) != 0 ) continue;
+		write_to_pgraph(fh, crystals[i].refls, full,
+		                crystals[i].cr, i, iter);
 	}
 
 	fclose(fh);
 }
 
 
-static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
+static void all_residuals(struct crystal_refls *crystals, int n_crystals, RefList *full,
                           int no_free,
                           double *presidual, double *pfree_residual,
                           double *plog_residual, double *pfree_log_residual,
@@ -798,28 +797,28 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 		double r, free_r, log_r, free_log_r;
 		int n;
 
-		if ( crystal_get_user_flag(crystals[i]) ) continue;
+		if ( crystal_get_user_flag(crystals[i].cr) ) continue;
 
 		/* Scaling should have been done right before calling this */
-		r = residual(crystals[i], full, 0, &n, NULL);
+		r = residual(crystals[i].refls, crystals[i].cr, full, 0, &n, NULL);
 		if ( n == 0 ) {
 			n_non_linear++;
 		} else if ( isnan(r) ) {
 			n_nan_linear++;
 		}
-		free_r = residual(crystals[i], full, 1, &n, NULL);
+		free_r = residual(crystals[i].refls, crystals[i].cr, full, 1, &n, NULL);
 		if ( n == 0 ) {
 			n_non_linear_free++;
 		} else if ( isnan(free_r) ) {
 			n_nan_linear_free++;
 		}
-		log_r = log_residual(crystals[i], full, 0, &n, NULL);
+		log_r = log_residual(crystals[i].refls, crystals[i].cr, full, 0, &n, NULL);
 		if ( n == 0 ) {
 			n_non_log++;
 		} else if ( isnan(log_r) ) {
 			n_nan_log++;
 		}
-		free_log_r = log_residual(crystals[i], full, 1, &n, NULL);
+		free_log_r = log_residual(crystals[i].refls, crystals[i].cr, full, 1, &n, NULL);
 		if ( n == 0 ) {
 			n_non_log_free++;
 		} else if ( isnan(free_log_r) ) {
@@ -876,7 +875,7 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 }
 
 
-static void show_all_residuals(Crystal **crystals, int n_crystals,
+static void show_all_residuals(struct crystal_refls *crystals, int n_crystals,
                                RefList *full, int no_free)
 {
 	double dev, free_dev, log_dev, free_log_dev;
@@ -896,7 +895,8 @@ struct log_qargs
 {
 	int iter;
 	int next;
-	Crystal **crystals;
+	struct crystal_refls *crystals;
+	struct image **images;
 	int n_crystals;
 	RefList *full;
 	int scaleflags;
@@ -909,6 +909,8 @@ struct log_qargs
 struct log_args
 {
 	Crystal *cr;
+	RefList *refls;
+	struct image *image;
 	RefList *full;
 	int scaleflags;
 	PartialityModel pmodel;
@@ -928,7 +930,9 @@ static void *get_log_task(void *vp)
 	task = malloc(sizeof(struct log_args));
 	if ( task == NULL ) return NULL;
 
-	task->cr = qargs->crystals[qargs->next];
+	task->cr = qargs->crystals[qargs->next].cr;
+	task->refls = qargs->crystals[qargs->next].refls;
+	task->image = qargs->images[qargs->next];
 	task->full = qargs->full;
 	task->iter = qargs->iter;
 	task->cnum = qargs->next;
@@ -944,12 +948,13 @@ static void *get_log_task(void *vp)
 static void write_logs(void *vp, int cookie)
 {
 	struct log_args *args = vp;
-	write_specgraph(args->cr, args->full, args->iter, args->cnum,
-	                args->log_folder);
-	write_gridscan(args->cr, args->full, args->iter, args->cnum,
-	               args->scaleflags, args->pmodel, args->log_folder);
-	write_test_logs(args->cr, args->full, args->iter, args->cnum,
-	                args->log_folder);
+	write_specgraph(args->refls, args->cr, args->image, args->full,
+	                args->iter, args->cnum, args->log_folder);
+	write_gridscan(args->refls, args->cr, args->image, args->full,
+	               args->iter, args->cnum, args->scaleflags, args->pmodel,
+	               args->log_folder);
+	write_test_logs(args->cr, args->image, args->full, args->iter,
+	                args->cnum, args->log_folder);
 }
 
 
@@ -963,7 +968,8 @@ static void done_log(void *vqargs, void *vp)
 }
 
 
-static void write_logs_parallel(Crystal **crystals, int n_crystals,
+static void write_logs_parallel(struct crystal_refls *crystals, struct image **images,
+                                int n_crystals,
                                 RefList *full, int iter, int n_threads,
                                 int scaleflags, PartialityModel pmodel,
                                 const char *log_folder)
@@ -974,6 +980,7 @@ static void write_logs_parallel(Crystal **crystals, int n_crystals,
 	qargs.next = 0;
 	qargs.full = full;
 	qargs.crystals = crystals;
+	qargs.images = images;
 	qargs.n_done = 0;
 	qargs.n_crystals = n_crystals;
 	qargs.scaleflags = scaleflags;
@@ -1084,7 +1091,8 @@ int main(int argc, char *argv[])
 	int no_scale = 0;
 	int no_Bscale = 0;
 	int no_pr = 0;
-	Crystal **crystals;
+	struct crystal_refls *crystals;
+	struct image **images;
 	char *pmodel_str = NULL;
 	PartialityModel pmodel = PMODEL_XSPHERE;
 	int min_measurements = 2;
@@ -1526,6 +1534,7 @@ int main(int argc, char *argv[])
 	n_crystals = 0;
 	n_crystals_seen = 0;
 	crystals = NULL;
+	images = NULL;
 	if ( sparams_fn != NULL ) {
 		char line[1024];
 		sparams_fh = fopen(sparams_fn, "r");
@@ -1560,7 +1569,6 @@ int main(int argc, char *argv[])
 		do {
 
 			struct image *image;
-			RefList *as;
 			int i;
 
 			image = stream_read_chunk(st, STREAM_REFLECTIONS);
@@ -1574,37 +1582,44 @@ int main(int argc, char *argv[])
 			for ( i=0; i<image->n_crystals; i++ ) {
 
 				Crystal *cr;
-				Crystal **crystals_new;
+				struct crystal_refls *crystals_new;
+				struct image **images_new;
 				RefList *cr_refl;
-				RefList *cr_refl_raw;
 				struct image *image_for_crystal;
 				double lowest_r;
 
 				n_crystals_seen++;
 				if ( n_crystals_seen <= start_after ) continue;
 
-				if ( crystal_get_resolution_limit(image->crystals[i]) < min_res ) continue;
+				if ( crystal_get_resolution_limit(image->crystals[i].cr) < min_res ) continue;
 
-				lowest_r = lowest_reflection(crystal_get_cell(image->crystals[i]));
-				if ( crystal_get_profile_radius(image->crystals[i]) > 0.5*lowest_r ) {
+				lowest_r = lowest_reflection(crystal_get_cell(image->crystals[i].cr));
+				if ( crystal_get_profile_radius(image->crystals[i].cr) > 0.5*lowest_r ) {
 					ERROR("Rejecting %s %s crystal %i because "
 					      "profile radius is obviously too big (%e %e).\n",
 					      image->filename, image->ev, i,
-					      crystal_get_profile_radius(image->crystals[i]),
+					      crystal_get_profile_radius(image->crystals[i].cr),
 					      lowest_r);
 					continue;
 				}
 
 				crystals_new = realloc(crystals,
-				                       (n_crystals+1)*sizeof(Crystal *));
+				                       (n_crystals+1)*sizeof(struct crystal_refls));
 				if ( crystals_new == NULL ) {
 					ERROR("Failed to allocate memory for crystal "
 					      "list.\n");
 					return 1;
 				}
 				crystals = crystals_new;
-				crystals[n_crystals] = crystal_copy_deep(image->crystals[i]);
-				cr = crystals[n_crystals];
+				cr = crystal_copy(image->crystals[i].cr);
+				crystals[n_crystals].cr = cr;
+
+				images_new = realloc(images, (n_crystals+1)*sizeof(struct image *));
+				if ( images_new == NULL ) {
+					ERROR("Failed to allocate memory for image list\n");
+					return 1;
+				}
+				images = images_new;
 
 				/* Create a completely new, separate image
 				 * structure for this crystal. */
@@ -1614,11 +1629,10 @@ int main(int argc, char *argv[])
 					return 1;
 				}
 
-				crystal_set_image(cr, image_for_crystal);
 				*image_for_crystal = *image;
-				image_for_crystal->n_crystals = 1;
-				image_for_crystal->crystals = malloc(sizeof(Crystal *));
-				image_for_crystal->crystals[0] = cr;
+				images[n_crystals] = image_for_crystal;
+				image_for_crystal->n_crystals = 0;
+				image_for_crystal->crystals = NULL;
 				image_for_crystal->filename = strdup(image->filename);
 				image_for_crystal->ev = safe_strdup(image->ev);
 				image_for_crystal->detgeom = NULL;
@@ -1629,20 +1643,18 @@ int main(int argc, char *argv[])
 				image_for_crystal->bad = NULL;
 				image_for_crystal->sat = NULL;
 
-				/* This is the raw list of reflections */
-				cr_refl_raw = crystal_get_reflections(cr);
-
-				cr_refl = apply_max_adu(cr_refl_raw, max_adu);
-				reflist_free(cr_refl_raw);
-
+				cr_refl = apply_max_adu(image->crystals[i].refls, max_adu);
 				if ( !no_free ) select_free_reflections(cr_refl, rng);
-
-				as = asymmetric_indices(cr_refl, sym);
-				crystal_set_reflections(cr, as);
-				crystal_set_user_flag(cr, PRFLAG_OK);
+				image_add_crystal_refls(image_for_crystal, cr,
+				                        asymmetric_indices(cr_refl, sym));
 				reflist_free(cr_refl);
 
-				if ( set_initial_params(cr, sparams_fh, force_bandwidth,
+				crystals[n_crystals].cr = image_for_crystal->crystals[0].cr;
+				crystals[n_crystals].refls = image_for_crystal->crystals[0].refls;
+
+				crystal_set_user_flag(cr, PRFLAG_OK);
+				if ( set_initial_params(cr, image_for_crystal,
+				                        sparams_fh, force_bandwidth,
 				                        force_radius, force_lambda) )
 				{
 					ERROR("Failed to set initial parameters\n");
@@ -1681,17 +1693,16 @@ int main(int argc, char *argv[])
 	STATUS("Initial partiality calculation...\n");
 	for ( icryst=0; icryst<n_crystals; icryst++ ) {
 
-		Crystal *cr = crystals[icryst];
-		update_predictions(cr);
+		Crystal *cr = crystals[icryst].cr;
+		RefList *refls = crystals[icryst].refls;
+		update_predictions(refls, cr, images[icryst]);
 
 		/* Polarisation correction requires kpred values */
-		polarisation_correction(crystal_get_reflections(cr),
-		                        crystal_get_cell(cr), polarisation);
-
-		calculate_partialities(cr, pmodel);
+		polarisation_correction(refls, crystal_get_cell(cr), polarisation);
+		calculate_partialities(refls, cr, images[icryst], pmodel);
 	}
 
-	if (csplit != NULL) check_csplit(crystals, n_crystals, csplit);
+	if (csplit != NULL) check_csplit(crystals, images, n_crystals, csplit);
 
 	/* Make a first pass at cutting out crap */
 	//STATUS("Early rejection...\n");
@@ -1716,7 +1727,7 @@ int main(int argc, char *argv[])
 
 	if ( do_write_logs ) {
 		write_pgraph(full, crystals, n_crystals, 0, "", log_folder);
-		write_logs_parallel(crystals, n_crystals, full, 0, nthreads,
+		write_logs_parallel(crystals, images, n_crystals, full, 0, nthreads,
 		                    scaleflags, pmodel, log_folder);
 	}
 
@@ -1726,7 +1737,7 @@ int main(int argc, char *argv[])
 		STATUS("Scaling and refinement cycle %i of %i\n", itn+1, n_iter);
 
 		if ( !no_pr ) {
-			refine_all(crystals, n_crystals, full, nthreads, pmodel,
+			refine_all(crystals, images, n_crystals, full, nthreads, pmodel,
 			           itn+1, no_logs, sym, amb, scaleflags,
 			           log_folder);
 		}
@@ -1771,7 +1782,7 @@ int main(int argc, char *argv[])
 				int j;
 				for ( j=0; j<csplit->n_datasets; j++ ) {
 					write_custom_split(csplit, j, crystals,
-					                   n_crystals, pmodel,
+					                   images, n_crystals, pmodel,
 					                   min_measurements,
 							   push_res, sym,
 					                   nthreads, tmp);
@@ -1805,7 +1816,7 @@ int main(int argc, char *argv[])
 	show_all_residuals(crystals, n_crystals, full, no_free);
 	if ( do_write_logs ) {
 		write_pgraph(full, crystals, n_crystals, -1, "", log_folder);
-		write_logs_parallel(crystals, n_crystals, full, -1, nthreads,
+		write_logs_parallel(crystals, images, n_crystals, full, -1, nthreads,
 		                    scaleflags, pmodel, log_folder);
 	}
 
@@ -1827,7 +1838,7 @@ int main(int argc, char *argv[])
 	if ( csplit != NULL ) {
 		int i;
 		for ( i=0; i<csplit->n_datasets; i++ ) {
-			write_custom_split(csplit, i, crystals, n_crystals,
+			write_custom_split(csplit, i, crystals, images, n_crystals,
 			                   pmodel, min_measurements, push_res,
 			                   sym, nthreads, outfile);
 		}
@@ -1835,10 +1846,6 @@ int main(int argc, char *argv[])
 
 	/* Clean up */
 	gsl_rng_free(rng);
-	for ( icryst=0; icryst<n_crystals; icryst++ ) {
-		struct image *image = crystal_get_image(crystals[icryst]);
-		image_free(image);
-	}
 	free_contribs(full);
 	reflist_free(full);
 	free_symoplist(sym);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 70b14a0f..4a02ed08 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -56,14 +56,16 @@ struct rf_alteration
 
 struct rf_priv
 {
-	Crystal *cr;      /**< Original crystal (before any refinement) */
+	Crystal *cr;         /**< Original crystal (before any refinement) */
+	struct image *image; /**< Original image (before any refinement) */
 	const RefList *full;
 	int serial;
 	int scaleflags;
 	PartialityModel pmodel;
 
-	Crystal *cr_tgt;        /**< Crystal to use for testing modifications */
-	struct image image_tgt; /**< Image structure to go with cr_tgt */
+	Crystal *cr_tgt;         /**< Crystal to use for testing modifications */
+	struct image image_tgt;  /**< Image structure to go with cr_tgt */
+	RefList *refls;          /**< The reflections to use */
 };
 
 
@@ -151,32 +153,30 @@ static void rotate_cell_xy(UnitCell *source, UnitCell *tgt,
 
 
 static void apply_parameters(Crystal *cr_orig, Crystal *cr_tgt,
+                             struct image *im_orig, struct image *im_tgt,
                              struct rf_alteration alter)
 {
-	double R, lambda;
+	double R;
 	struct gaussian g;
-	struct image *image;
 
-	image = crystal_get_image(cr_tgt);
 	R = crystal_get_profile_radius(cr_orig);
-	lambda = crystal_get_image_const(cr_orig)->lambda;
 
 	rotate_cell_xy(crystal_get_cell(cr_orig), crystal_get_cell(cr_tgt),
 	               alter.rot_x, alter.rot_y);
 	crystal_set_profile_radius(cr_tgt, R+alter.delta_R);
-	image->lambda = lambda+alter.delta_wave;
+	im_tgt->lambda = im_orig->lambda+alter.delta_wave;
 
-	g.kcen = 1.0/image->lambda;
-	g.sigma = image->bw/image->lambda;
+	g.kcen = 1.0/im_tgt->lambda;
+	g.sigma = im_tgt->bw/im_tgt->lambda;
 	g.area = 1;
-	spectrum_set_gaussians(image->spectrum, &g, 1);
+	spectrum_set_gaussians(im_tgt->spectrum, &g, 1);
 }
 
 
 static double calc_residual(struct rf_priv *pv, struct rf_alteration alter,
                             int free)
 {
-	apply_parameters(pv->cr, pv->cr_tgt, alter);
+	apply_parameters(pv->cr, pv->cr_tgt, pv->image, &pv->image_tgt, alter);
 
 	if ( fabs(crystal_get_profile_radius(pv->cr_tgt)) > 5e9 ) {
 		STATUS("radius > 5e9\n");
@@ -189,16 +189,16 @@ static double calc_residual(struct rf_priv *pv, struct rf_alteration alter,
 		return NAN;
 	}
 
-	if ( crystal_get_image(pv->cr_tgt)->lambda <= 0.0 ) {
+	if ( pv->image_tgt.lambda <= 0.0 ) {
 		STATUS("lambda < 0\n");
 		return NAN;
 	}
 
-	update_predictions(pv->cr_tgt);
-	calculate_partialities(pv->cr_tgt, pv->pmodel);
+	update_predictions(pv->refls, pv->cr_tgt, &pv->image_tgt);
+	calculate_partialities(pv->refls, pv->cr_tgt, &pv->image_tgt, pv->pmodel);
 
 
-	return residual(pv->cr_tgt, pv->full, free, NULL, NULL);
+	return residual(pv->refls, pv->cr_tgt, pv->full, free, NULL, NULL);
 }
 
 
@@ -264,64 +264,58 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
 }
 
 
-static void try_reindex(Crystal *crin, const RefList *full,
-                        SymOpList *sym, SymOpList *amb, int scaleflags,
-                        PartialityModel pmodel)
+static void try_reindex(RefList **prefls, Crystal *crin, struct image *image,
+                        const RefList *full, SymOpList *sym, SymOpList *amb,
+                        int scaleflags, PartialityModel pmodel)
 {
-	Crystal *cr;
 	double residual_original;
 	int idx, n;
 
 	if ( sym == NULL || amb == NULL ) return;
 
-	if ( scale_one_crystal(crin, full, scaleflags) ) return;
-	residual_original = residual(crin, full, 0, NULL, NULL);
-
-	cr = crystal_copy(crin);
+	if ( scale_one_crystal(*prefls, crin, full, scaleflags) ) return;
+	residual_original = residual(*prefls, crin, full, 0, NULL, NULL);
 
 	n = num_equivs(amb, NULL);
 
 	for ( idx=0; idx<n; idx++ ) {
 
 		RefList *list;
-		UnitCell *cell;
-		double residual_flipped;
-
-		cell = cell_new_from_cell(crystal_get_cell(crin));
-		if ( cell == NULL ) return;
-		reindex_cell(cell, amb, idx);
-		crystal_set_cell(cr, cell);
-
-		list = reindex_reflections(crystal_get_reflections(crin),
-		                           amb, sym, idx);
-		crystal_set_reflections(cr, list);
-
-		update_predictions(cr);
-		calculate_partialities(cr, pmodel);
-
-		if ( scale_one_crystal(cr, full, scaleflags) ) return;
-		residual_flipped = residual(cr, full, 0, NULL, NULL);
-
-		if ( residual_flipped < residual_original ) {
-			crystal_set_cell(crin, cell);
-			crystal_set_reflections(crin, list);
-			residual_original = residual_flipped;
-		} else {
-			cell_free(crystal_get_cell(cr));
-			reflist_free(crystal_get_reflections(cr));
+		Crystal *cr;
+
+		cr = crystal_copy(crin);
+
+		reindex_cell(crystal_get_cell(cr), amb, idx);
+		list = reindex_reflections(*prefls, amb, sym, idx);
+
+		update_predictions(list, cr, image);
+		calculate_partialities(list, cr, image, pmodel);
+		if ( scale_one_crystal(list, cr, full, scaleflags) == 0 ) {
+
+			double residual_flipped;
+			residual_flipped = residual(list, cr, full, 0, NULL, NULL);
+
+			if ( residual_flipped < residual_original ) {
+				crystal_set_cell(crin, crystal_relinquish_cell(cr));
+				reflist_free(*prefls);
+				*prefls = list;
+				residual_original = residual_flipped;
+			} else {
+				reflist_free(list);
+			}
 		}
+
+		crystal_free(cr);
 	}
 
-	crystal_free(cr);
 }
 
 
-void write_test_logs(Crystal *crystal, const RefList *full,
+void write_test_logs(Crystal *crystal, struct image *image, const RefList *full,
                      signed int cycle, int serial,
                      const char *log_folder)
 {
 	FILE *fh;
-	struct image *image = crystal_get_image(crystal);
 	char tmp[256];
 	char ins[16];
 
@@ -367,7 +361,7 @@ void write_test_logs(Crystal *crystal, const RefList *full,
 }
 
 
-void write_specgraph(Crystal *crystal, const RefList *full,
+void write_specgraph(RefList *list, Crystal *crystal, struct image *image, const RefList *full,
                      signed int cycle, int serial,
                      const char *log_folder)
 {
@@ -378,7 +372,6 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 	double G = crystal_get_osf(crystal);
 	double B = crystal_get_Bfac(crystal);
 	UnitCell *cell;
-	struct image *image = crystal_get_image(crystal);
 	char ins[16];
 
 	snprintf(tmp, 256, "%s/specgraph-crystal%i.dat", log_folder, serial);
@@ -408,7 +401,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 		ins[1] = '\0';
 	}
 
-	for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
+	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -443,7 +436,8 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 }
 
 
-static void write_angle_grid(Crystal *cr, const RefList *full,
+static void write_angle_grid(RefList *list_in, Crystal *cr, struct image *image,
+                             const RefList *full,
                              signed int cycle, int serial, int scaleflags,
                              PartialityModel pmodel,
                              const char *log_folder)
@@ -452,7 +446,6 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
 	char fn[64];
 	char ins[16];
 	struct rf_priv priv;
-	RefList *list;
 	UnitCell *cell;
 	Spectrum *spectrum;
 
@@ -462,12 +455,11 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
 	priv.scaleflags = scaleflags;
 	priv.pmodel = pmodel;
 	priv.cr_tgt = crystal_copy(cr);
-	priv.image_tgt = *crystal_get_image(cr);
+	priv.image = image;
+	priv.image_tgt = *image;
 	spectrum = spectrum_new();
 	priv.image_tgt.spectrum = spectrum;
-	crystal_set_image(priv.cr_tgt, &priv.image_tgt);
-	list = copy_reflist(crystal_get_reflections(cr));
-	crystal_set_reflections(priv.cr_tgt, list);
+	priv.refls = copy_reflist(list_in);
 	cell = cell_new_from_cell(crystal_get_cell(cr));
 	crystal_set_cell(priv.cr_tgt, cell);
 
@@ -504,13 +496,14 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
 		fclose(fh);
 	}
 
-	reflist_free(crystal_get_reflections(priv.cr_tgt));
+	reflist_free(priv.refls);
 	crystal_free(priv.cr_tgt);
 	spectrum_free(spectrum);
 }
 
 
-static void write_radius_grid(Crystal *cr, const RefList *full,
+static void write_radius_grid(RefList *list_in, Crystal *cr, struct image *image,
+                              const RefList *full,
                               signed int cycle, int serial, int scaleflags,
                               PartialityModel pmodel,
                               const char *log_folder)
@@ -519,7 +512,6 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
 	char fn[64];
 	char ins[16];
 	struct rf_priv priv;
-	RefList *list;
 	UnitCell *cell;
 	Spectrum *spectrum;
 
@@ -529,12 +521,11 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
 	priv.scaleflags = scaleflags;
 	priv.pmodel = pmodel;
 	priv.cr_tgt = crystal_copy(cr);
-	priv.image_tgt = *crystal_get_image(cr);
+	priv.image = image;
+	priv.image_tgt = *image;
 	spectrum = spectrum_new();
 	priv.image_tgt.spectrum = spectrum;
-	crystal_set_image(priv.cr_tgt, &priv.image_tgt);
-	list = copy_reflist(crystal_get_reflections(cr));
-	crystal_set_reflections(priv.cr_tgt, list);
+	priv.refls = copy_reflist(list_in);
 	cell = cell_new_from_cell(crystal_get_cell(cr));
 	crystal_set_cell(priv.cr_tgt, cell);
 
@@ -571,19 +562,20 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
 		fclose(fh);
 	}
 
-	reflist_free(crystal_get_reflections(priv.cr_tgt));
+	reflist_free(priv.refls);
 	crystal_free(priv.cr_tgt);
 	spectrum_free(spectrum);
 }
 
 
-void write_gridscan(Crystal *cr, const RefList *full,
+void write_gridscan(RefList *list,
+                    Crystal *cr, struct image *image, const RefList *full,
                     signed int cycle, int serial, int scaleflags,
                     PartialityModel pmodel,
                     const char *log_folder)
 {
-	write_angle_grid(cr, full, cycle, serial, scaleflags, pmodel, log_folder);
-	write_radius_grid(cr, full, cycle, serial, scaleflags, pmodel, log_folder);
+	write_angle_grid(list, cr, image, full, cycle, serial, scaleflags, pmodel, log_folder);
+	write_radius_grid(list, cr, image, full, cycle, serial, scaleflags, pmodel, log_folder);
 }
 
 
@@ -631,7 +623,7 @@ static int refine_loop(struct rf_alteration *cur, struct rf_alteration *dirns,
 			        *total_iter, lowest_fom, freefom,
 			        cur->rot_x, cur->rot_y,
 			        crystal_get_profile_radius(priv->cr)+cur->delta_R,
-			        crystal_get_image_const(priv->cr)->lambda+cur->delta_wave);
+			        priv->image->lambda+cur->delta_wave);
 		}
 
 	} while ( (best != 0) && (n_iter < 10) );
@@ -649,7 +641,8 @@ static void zero_alter(struct rf_alteration *alter)
 }
 
 
-static void do_pr_refine(Crystal *cr, const RefList *full,
+static void do_pr_refine(RefList **plist_in, Crystal *cr, struct image *image,
+                         const RefList *full,
                          PartialityModel pmodel, int serial,
                          int cycle, int write_logs,
                          SymOpList *sym, SymOpList *amb, int scaleflags,
@@ -659,12 +652,11 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	struct rf_alteration alter;
 	int n_iter = 0;
 	double fom, freefom;
-	RefList *list;
 	FILE *fh = NULL;
 	UnitCell *cell;
 	Spectrum *spectrum;
 
-	try_reindex(cr, full, sym, amb, scaleflags, pmodel);
+	try_reindex(plist_in, cr, image, full, sym, amb, scaleflags, pmodel);
 
 	zero_alter(&alter);
 
@@ -674,12 +666,11 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	priv.scaleflags = scaleflags;
 	priv.pmodel = pmodel;
 	priv.cr_tgt = crystal_copy(cr);
-	priv.image_tgt = *crystal_get_image(cr);
+	priv.image = image;
+	priv.image_tgt = *image;
 	spectrum = spectrum_new();
 	priv.image_tgt.spectrum = spectrum;
-	crystal_set_image(priv.cr_tgt, &priv.image_tgt);
-	list = copy_reflist(crystal_get_reflections(cr));
-	crystal_set_reflections(priv.cr_tgt, list);
+	priv.refls = copy_reflist(*plist_in);
 	cell = cell_new_from_cell(crystal_get_cell(cr));
 	crystal_set_cell(priv.cr_tgt, cell);
 
@@ -700,7 +691,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 			        n_iter, fom, freefom,
 			        alter.rot_x, alter.rot_y,
 			        crystal_get_profile_radius(cr)+alter.delta_R,
-			        crystal_get_image(cr)->lambda+alter.delta_wave);
+			        image->lambda+alter.delta_wave);
 		}
 
 	}
@@ -724,15 +715,15 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	refine_loop(&alter, dirns, 2, &priv, &n_iter, fh);
 
 	/* Apply the final shifts */
-	apply_parameters(cr, cr, alter);
-	update_predictions(cr);
-	calculate_partialities(cr, pmodel);
+	apply_parameters(cr, cr, image, image, alter);
+	update_predictions(*plist_in, cr, image);
+	calculate_partialities(*plist_in, cr, image, pmodel);
 
 	if ( write_logs ) {
-		write_gridscan(cr, full, cycle, serial, scaleflags,
+		write_gridscan(*plist_in, cr, image, full, cycle, serial, scaleflags,
 		               pmodel, log_folder);
-		write_specgraph(cr, full, cycle, serial, log_folder);
-		write_test_logs(cr, full, cycle, serial, log_folder);
+		write_specgraph(*plist_in, cr, image, full, cycle, serial, log_folder);
+		write_test_logs(cr, image, full, cycle, serial, log_folder);
 	}
 
 	if ( crystal_get_profile_radius(cr) > 5e9 ) {
@@ -743,7 +734,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 		fclose(fh);
 	}
 
-	reflist_free(crystal_get_reflections(priv.cr_tgt));
+	reflist_free(priv.refls);
 	crystal_free(priv.cr_tgt);
 	spectrum_free(spectrum);
 }
@@ -753,6 +744,8 @@ struct refine_args
 {
 	RefList *full;
 	Crystal *crystal;
+	RefList **prefls;
+	struct image *image;
 	PartialityModel pmodel;
 	int serial;
 	int cycle;
@@ -768,7 +761,8 @@ struct pr_queue_args
 {
 	int n_started;
 	int n_done;
-	Crystal **crystals;
+	struct crystal_refls *crystals;
+	struct image **images;
 	int n_crystals;
 	struct refine_args task_defaults;
 };
@@ -777,12 +771,12 @@ struct pr_queue_args
 static void refine_image(void *task, int id)
 {
 	struct refine_args *pargs = task;
-	Crystal *cr = pargs->crystal;
 	int write_logs = 0;
 
 	write_logs = !pargs->no_logs && (pargs->serial % 20 == 0);
 
-	do_pr_refine(cr, pargs->full, pargs->pmodel,
+	do_pr_refine(pargs->prefls, pargs->crystal, pargs->image,
+	             pargs->full, pargs->pmodel,
 	             pargs->serial, pargs->cycle, write_logs,
 	             pargs->sym, pargs->amb, pargs->scaleflags,
 	             pargs->log_folder);
@@ -797,7 +791,9 @@ static void *get_image(void *vqargs)
 	task = malloc(sizeof(struct refine_args));
 	memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
 
-	task->crystal = qargs->crystals[qargs->n_started];
+	task->crystal = qargs->crystals[qargs->n_started].cr;
+	task->prefls = &qargs->crystals[qargs->n_started].refls;
+	task->image = qargs->images[qargs->n_started];
 	task->serial = qargs->n_started;
 
 	qargs->n_started++;
@@ -817,7 +813,7 @@ static void done_image(void *vqargs, void *task)
 }
 
 
-void refine_all(Crystal **crystals, int n_crystals,
+void refine_all(struct crystal_refls *crystals, struct image **images, int n_crystals,
                 RefList *full, int nthreads, PartialityModel pmodel,
                 int cycle, int no_logs,
                 SymOpList *sym, SymOpList *amb, int scaleflags,
@@ -828,6 +824,8 @@ void refine_all(Crystal **crystals, int n_crystals,
 
 	task_defaults.full = full;
 	task_defaults.crystal = NULL;
+	task_defaults.prefls = NULL;
+	task_defaults.image = NULL;
 	task_defaults.pmodel = pmodel;
 	task_defaults.cycle = cycle;
 	task_defaults.no_logs = no_logs;
@@ -842,6 +840,7 @@ void refine_all(Crystal **crystals, int n_crystals,
 	qargs.n_done = 0;
 	qargs.n_crystals = n_crystals;
 	qargs.crystals = crystals;
+	qargs.images = images;
 
 	/* Don't have threads which are doing nothing */
 	if ( n_crystals < nthreads ) nthreads = n_crystals;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 546050a4..91c0f4f8 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -58,18 +58,21 @@ enum prflag
 
 extern const char *str_prflag(enum prflag flag);
 
-extern void refine_all(Crystal **crystals, int n_crystals,
+extern void refine_all(struct crystal_refls *crystals, struct image **images,
+                       int n_crystals,
                        RefList *full, int nthreads, PartialityModel pmodel,
                        int cycle, int no_logs,
                        SymOpList *sym, SymOpList *amb, int scaleflags,
                        const char *log_folder);
 
-extern void write_gridscan(Crystal *cr, const RefList *full,
+extern void write_gridscan(RefList *list, Crystal *cr, struct image *image,
+                           const RefList *full,
                            int cycle, int serial, int scaleflags,
                            PartialityModel model,
                            const char *log_folder);
 
-extern void write_specgraph(Crystal *crystal, const RefList *full,
+extern void write_specgraph(RefList *list, Crystal *crystal, struct image *image,
+                            const RefList *full,
                             signed int cycle, int serial,
                             const char *log_folder);
 
@@ -77,7 +80,8 @@ extern void write_specgraph(Crystal *crystal, const RefList *full,
 extern double gradient(Crystal *cr, int k, Reflection *refl,
                        PartialityModel pmodel);
 
-extern void write_test_logs(Crystal *crystal, const RefList *full,
+extern void write_test_logs(Crystal *crystal, struct image *image,
+                            const RefList *full,
                             signed int cycle, int serial,
                             const char *log_folder);
 
diff --git a/src/process_hkl.c b/src/process_hkl.c
index be1c95d4..31e4b1fc 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -271,7 +271,7 @@ static void apply_kpred(double k, RefList *list)
 
 
 static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
-                         RefList *reference, const SymOpList *sym,
+                         RefList *new_refl, RefList *reference, const SymOpList *sym,
                          double **hist_vals, signed int hist_h,
                          signed int hist_k, signed int hist_l, int *hist_n,
                          struct polarisation p, double min_snr, double max_adu,
@@ -280,11 +280,8 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
 {
 	Reflection *refl;
 	RefListIterator *iter;
-	RefList *new_refl;
 	double scale;
 
-	new_refl = crystal_get_reflections(cr);
-
 	/* First, correct for polarisation */
 	apply_kpred(1.0/image->lambda, new_refl);
 	polarisation_correction(new_refl, crystal_get_cell(cr), p);
@@ -309,7 +306,7 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
 		scale = 1.0;
 	}
 
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	for ( refl = first_refl(new_refl, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -431,7 +428,8 @@ static int merge_stream(Stream *st,
 
 		for ( i=0; i<image->n_crystals; i++ ) {
 
-			Crystal *cr = image->crystals[i];
+			Crystal *cr = image->crystals[i].cr;
+			RefList *refls = image->crystals[i].refls;
 
 			n_crystals_seen++;
 			if ( (n_crystals_seen > start_after)
@@ -440,8 +438,8 @@ static int merge_stream(Stream *st,
 			{
 				int r;
 				n_crystals++;
-				r = merge_crystal(model, image, cr, reference,
-				                  sym, hist_vals,
+				r = merge_crystal(model, image, cr, refls,
+				                  reference, sym, hist_vals,
 						  hist_h, hist_k, hist_l,
 				                  hist_i, p,
 						  min_snr, max_adu, push_res,
diff --git a/src/process_image.c b/src/process_image.c
index 2f856594..c26842c3 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -57,6 +57,8 @@
 #include "im-sandbox.h"
 #include "im-zmq.h"
 #include "im-asapo.h"
+#include "peaks.h"
+#include "peakfinder8.h"
 
 static float **backup_image_data(float **dp, struct detgeom *det)
 {
@@ -292,35 +294,42 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	profile_start("peak-search");
 	switch ( iargs->peak_search.method ) {
 
+		ImageFeatureList *peaks;
+
 		case PEAK_HDF5:
 		case PEAK_CXI:
 		set_last_task(last_task, "peaksearch:hdf5orcxi");
-		image->features = image_read_peaks(iargs->dtempl,
-		                                   pargs->filename,
-		                                   pargs->event,
-		                                   iargs->peak_search.half_pixel_shift);
-		if ( image->features == NULL ) {
-			ERROR("Failed to get peaks from HDF5 file.\n");
-		}
+		peaks = image_read_peaks(iargs->dtempl,
+		                         pargs->filename,
+		                         pargs->event,
+		                         iargs->peak_search.half_pixel_shift);
 		if ( iargs->peak_search.revalidate ) {
-			validate_peaks(image, iargs->peak_search.min_snr,
-				       iargs->peak_search.pk_inn, iargs->peak_search.pk_mid,
-		                       iargs->peak_search.pk_out, iargs->peak_search.use_saturated,
-				       iargs->peak_search.check_hdf5_snr);
+			ImageFeatureList *npeaks = validate_peaks(image, peaks,
+			                                          iargs->peak_search.min_snr,
+			                                          iargs->peak_search.pk_inn,
+			                                          iargs->peak_search.pk_mid,
+			                                          iargs->peak_search.pk_out,
+			                                          iargs->peak_search.use_saturated,
+			                                          iargs->peak_search.check_hdf5_snr);
+			image_feature_list_free(peaks);
+			image->features = npeaks;
 		}
 		break;
 
 		case PEAK_ZAEF:
 		set_last_task(last_task, "peaksearch:zaef");
-		search_peaks(image, iargs->peak_search.threshold,
-		             iargs->peak_search.min_sq_gradient, iargs->peak_search.min_snr,
-		             iargs->peak_search.pk_inn, iargs->peak_search.pk_mid, iargs->peak_search.pk_out,
-		             iargs->peak_search.use_saturated);
+		image->features = search_peaks(image, iargs->peak_search.threshold,
+		                               iargs->peak_search.min_sq_gradient,
+		                               iargs->peak_search.min_snr,
+		                               iargs->peak_search.pk_inn,
+		                               iargs->peak_search.pk_mid,
+		                               iargs->peak_search.pk_out,
+		                               iargs->peak_search.use_saturated);
 		break;
 
 		case PEAK_PEAKFINDER8:
 		set_last_task(last_task, "peaksearch:pf8");
-		if ( search_peaks_peakfinder8(image, 2048,
+		image->features = peakfinder8(image, 2048,
 		                              iargs->peak_search.threshold,
 		                              iargs->peak_search.min_snr,
 		                              iargs->peak_search.min_pix_count,
@@ -329,40 +338,36 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		                              iargs->peak_search.min_res,
 		                              iargs->peak_search.max_res,
 		                              iargs->peak_search.use_saturated,
-					      iargs->peak_search.peakfinder8_fast,
-					      iargs->pf_private) ) {
-			ERROR("Failed to find peaks in image %s"
-			      "(event %s).\n",
-			      image->filename, image->ev);
-		}
+		                              iargs->peak_search.peakfinder8_fast,
+		                              iargs->pf_private);
 		break;
 
 		case PEAK_PEAKFINDER9:
 		set_last_task(last_task, "peaksearch:pf9");
-		if ( search_peaks_peakfinder9(image,
-		                              iargs->peak_search.min_snr_biggest_pix,
-		                              iargs->peak_search.min_snr_peak_pix,
-		                              iargs->peak_search.min_snr,
-		                              iargs->peak_search.min_sig,
-		                              iargs->peak_search.min_peak_over_neighbour,
-		                              iargs->peak_search.local_bg_radius) )
-		{
-			ERROR("Failed to find peaks in image %s"
-			      "(event %s).\n",
-			      image->filename, image->ev);
-		}
+		image->features = search_peaks_peakfinder9(image,
+		                                           iargs->peak_search.min_snr_biggest_pix,
+		                                           iargs->peak_search.min_snr_peak_pix,
+		                                           iargs->peak_search.min_snr,
+		                                           iargs->peak_search.min_sig,
+		                                           iargs->peak_search.min_peak_over_neighbour,
+		                                           iargs->peak_search.local_bg_radius);
 		break;
 
 		case PEAK_MSGPACK:
-		image->features = image_msgpack_read_peaks(iargs->dtempl,
-		                                           pargs->zmq_data,
-		                                           pargs->zmq_data_size,
-		                                           iargs->peak_search.half_pixel_shift);
+		peaks = image_msgpack_read_peaks(iargs->dtempl,
+		                                 pargs->zmq_data,
+		                                 pargs->zmq_data_size,
+		                                 iargs->peak_search.half_pixel_shift);
 		if ( iargs->peak_search.revalidate ) {
-			validate_peaks(image, iargs->peak_search.min_snr,
-				       iargs->peak_search.pk_inn, iargs->peak_search.pk_mid,
-		                       iargs->peak_search.pk_out, iargs->peak_search.use_saturated,
-				       iargs->peak_search.check_hdf5_snr);
+			ImageFeatureList *npeaks = validate_peaks(image, peaks,
+			                                          iargs->peak_search.min_snr,
+			                                          iargs->peak_search.pk_inn,
+			                                          iargs->peak_search.pk_mid,
+			                                          iargs->peak_search.pk_out,
+			                                          iargs->peak_search.use_saturated,
+			                                          iargs->peak_search.check_hdf5_snr);
+			image_feature_list_free(peaks);
+			image->features = npeaks;
 		}
 		break;
 
@@ -371,6 +376,10 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		break;
 
 	}
+	if ( image->features == NULL ) {
+		ERROR("Peaksearch failed for image %s" "(event %s).\n",
+		      image->filename, image->ev);
+	}
 	profile_end("peak-search");
 
 	image->peak_resolution = estimate_peak_resolution(image->features,
@@ -434,20 +443,20 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	set_last_task(last_task, "prediction params");
 	if ( iargs->fix_profile_r >= 0.0 ) {
 		for ( i=0; i<image->n_crystals; i++ ) {
-			crystal_set_profile_radius(image->crystals[i],
+			crystal_set_profile_radius(image->crystals[i].cr,
 			                           iargs->fix_profile_r);
-			crystal_set_mosaicity(image->crystals[i], 0.0);
+			crystal_set_mosaicity(image->crystals[i].cr, 0.0);
 		}
 	} else {
 		for ( i=0; i<image->n_crystals; i++ ) {
-			crystal_set_profile_radius(image->crystals[i], 0.02e9);
-			crystal_set_mosaicity(image->crystals[i], 0.0);
+			crystal_set_profile_radius(image->crystals[i].cr, 0.02e9);
+			crystal_set_mosaicity(image->crystals[i].cr, 0.0);
 		}
 	}
 
 	if ( iargs->fix_profile_r < 0.0 ) {
 		for ( i=0; i<image->n_crystals; i++ ) {
-			if ( refine_radius(image->crystals[i], image) ) {
+			if ( refine_radius(image->crystals[i].cr, image) ) {
 				ERROR("WARNING: Radius determination failed\n");
 			}
 		}
@@ -479,7 +488,7 @@ streamwrite:
 
 	int n = 0;
 	for ( i=0; i<image->n_crystals; i++ ) {
-		n += crystal_get_num_implausible_reflections(image->crystals[i]);
+		n += crystal_get_num_implausible_reflections(image->crystals[i].cr);
 	}
 	if ( n > 0 ) {
 		STATUS("WARNING: %i implausibly negative reflection%s in %s "
@@ -495,7 +504,7 @@ out:
 	pthread_mutex_lock(&sb_shared->totals_lock);
 	any_crystals = 0;
 	for ( i=0; i<image->n_crystals; i++ ) {
-		if ( crystal_get_user_flag(image->crystals[i]) == 0 ) {
+		if ( crystal_get_user_flag(image->crystals[i].cr) == 0 ) {
 			sb_shared->n_crystals++;
 			any_crystals = 1;
 		}
diff --git a/src/rejection.c b/src/rejection.c
index 4e3ff6ba..d3aa9b1f 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -63,7 +63,7 @@ static double mean_intensity(RefList *list)
 
 
 /* Reject really obvious outliers */
-void early_rejection(Crystal **crystals, int n)
+void early_rejection(struct crystal_refls *crystals, int n)
 {
 	int i;
 	double m = 0.0;
@@ -73,18 +73,16 @@ void early_rejection(Crystal **crystals, int n)
 
 	for ( i=0; i<n; i++ ) {
 		double u;
-		RefList *list = crystal_get_reflections(crystals[i]);
-		u = mean_intensity(list);
+		u = mean_intensity(crystals[i].refls);
 		m += u;
 		fprintf(fh, "%i %f\n", i, u);
 	}
 	mean_m = m/n;
 	for ( i=0; i<n; i++ ) {
 		double u;
-		RefList *list = crystal_get_reflections(crystals[i]);
-		u = mean_intensity(list);
+		u = mean_intensity(crystals[i].refls);
 		if ( u/mean_m < 0.2 ) {
-			crystal_set_user_flag(crystals[i], 5);
+			crystal_set_user_flag(crystals[i].cr, 5);
 			n_flag++;
 		}
 	}
@@ -219,7 +217,7 @@ static double calculate_cchalf(RefList *template, RefList *full,
 		 * If the crystal is marked as bad, we should not remove it
 		 * because it did not contribute in the first place.  */
 		if ( exclude != NULL && !crystal_get_user_flag(exclude) ) {
-			exrefl = find_refl(crystal_get_reflections(exclude), h, k, l);
+			exrefl = find_refl(template, h, k, l);
 		} else {
 			exrefl = NULL;
 		}
@@ -275,7 +273,7 @@ static double calculate_cchalf(RefList *template, RefList *full,
 struct deltacchalf_queue_args
 {
 	RefList *full;
-	Crystal **crystals;
+	struct crystal_refls *crystals;
 	int n_crystals;
 	int n_done;
 	int n_started;
@@ -289,6 +287,7 @@ struct deltacchalf_worker_args
 {
 	RefList *full;
 	Crystal *crystal;
+	RefList *refls;
 	int crystal_number;
 	int non;
 	int nan;
@@ -304,7 +303,8 @@ static void *create_deltacchalf_job(void *vqargs)
 	wargs = malloc(sizeof(struct deltacchalf_worker_args));
 
 	wargs->full = qargs->full;
-	wargs->crystal = qargs->crystals[qargs->n_started];
+	wargs->crystal = qargs->crystals[qargs->n_started].cr;
+	wargs->refls = qargs->crystals[qargs->n_started].refls;
 	wargs->crystal_number = qargs->n_started;
 	wargs->non = 0;
 	wargs->nan = 0;
@@ -320,9 +320,8 @@ static void run_deltacchalf_job(void *vwargs, int cookie)
 	double cchalf, cchalfi;
 	struct deltacchalf_worker_args *wargs = vwargs;
 	int nref = 0;
-	RefList *template = crystal_get_reflections(wargs->crystal);
-	cchalf = calculate_cchalf(template, wargs->full, NULL, &nref);
-	cchalfi = calculate_cchalf(template, wargs->full, wargs->crystal, &nref);
+	cchalf = calculate_cchalf(wargs->refls, wargs->full, NULL, &nref);
+	cchalfi = calculate_cchalf(wargs->refls, wargs->full, wargs->crystal, &nref);
 	//STATUS("Frame %i:", i);
 	//STATUS("   With = %f  ", cchalf*100.0);
 	//STATUS("Without = %f", cchalfi*100.0);
@@ -355,8 +354,8 @@ static void finalise_deltacchalf_job(void *vqargs, void *vwargs)
 }
 
 
-static void check_deltacchalf(Crystal **crystals, int n, RefList *full,
-                              int n_threads)
+static void check_deltacchalf(struct crystal_refls *crystals, int n,
+                              RefList *full, int n_threads)
 {
 	double cchalf;
 	int i;
@@ -409,7 +408,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full,
 		if ( isnan(vals[i]) || isinf(vals[i])
 		  || ((vals[i]<0.0) && (vals[i] < mean-2.0*sd)) )
 		{
-			crystal_set_user_flag(crystals[i], PRFLAG_DELTACCHALF);
+			crystal_set_user_flag(crystals[i].cr, PRFLAG_DELTACCHALF);
 		}
 	}
 
@@ -417,7 +416,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full,
 }
 
 
-static void show_duds(Crystal **crystals, int n_crystals)
+static void show_duds(struct crystal_refls *crystals, int n_crystals)
 {
 	int j;
 	int bads[32];
@@ -427,7 +426,7 @@ static void show_duds(Crystal **crystals, int n_crystals)
 
 	for ( j=0; j<n_crystals; j++ ) {
 		int flag;
-		flag = crystal_get_user_flag(crystals[j]);
+		flag = crystal_get_user_flag(crystals[j].cr);
 		assert(flag < 32);
 		bads[flag]++;
 		if ( flag != PRFLAG_OK ) any_bad++;
@@ -444,8 +443,8 @@ static void show_duds(Crystal **crystals, int n_crystals)
 }
 
 
-void check_rejection(Crystal **crystals, int n, RefList *full, double max_B,
-                     int no_deltacchalf, int n_threads)
+void check_rejection(struct crystal_refls *crystals, int n, RefList *full,
+                     double max_B, int no_deltacchalf, int n_threads)
 {
 	int i;
 	int n_acc = 0;
@@ -456,10 +455,10 @@ void check_rejection(Crystal **crystals, int n, RefList *full, double max_B,
 	}
 
 	for ( i=0; i<n; i++ ) {
-		if ( fabs(crystal_get_Bfac(crystals[i])) > max_B ) {
-			crystal_set_user_flag(crystals[i], PRFLAG_BIGB);
+		if ( fabs(crystal_get_Bfac(crystals[i].cr)) > max_B ) {
+			crystal_set_user_flag(crystals[i].cr, PRFLAG_BIGB);
 		}
-		if ( crystal_get_user_flag(crystals[i]) == 0 ) n_acc++;
+		if ( crystal_get_user_flag(crystals[i].cr) == 0 ) n_acc++;
 	}
 
 	show_duds(crystals, n);
diff --git a/src/rejection.h b/src/rejection.h
index 7dc2c92b..152c004e 100644
--- a/src/rejection.h
+++ b/src/rejection.h
@@ -35,10 +35,10 @@
 #endif
 
 
-#include "crystal.h"
+#include "image.h"
 
-extern void early_rejection(Crystal **crystals, int n);
-extern void check_rejection(Crystal **crystals, int n, RefList *full,
+extern void early_rejection(struct crystal_refls *crystals, int n);
+extern void check_rejection(struct crystal_refls *crystals, int n, RefList *full,
                             double max_B, int no_deltacchalf, int n_threads);
 
 #endif	/* REJECTION_H */
diff --git a/src/scaling.c b/src/scaling.c
index f007b757..52954f9a 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -53,6 +53,7 @@ struct scale_args
 {
 	RefList *full;
 	Crystal *crystal;
+	RefList *refls;
 	int flags;
 };
 
@@ -61,7 +62,7 @@ struct scale_queue_args
 {
 	int n_started;
 	int n_done;
-	Crystal **crystals;
+	struct crystal_refls *crystals;
 	int n_crystals;
 	struct scale_args task_defaults;
 };
@@ -70,7 +71,7 @@ struct scale_queue_args
 static void scale_crystal(void *task, int id)
 {
 	struct scale_args *pargs = task;
-	scale_one_crystal(pargs->crystal, pargs->full, pargs->flags);
+	scale_one_crystal(pargs->refls, pargs->crystal, pargs->full, pargs->flags);
 }
 
 
@@ -82,7 +83,8 @@ static void *get_crystal(void *vqargs)
 	task = malloc(sizeof(struct scale_args));
 	memcpy(task, &qargs->task_defaults, sizeof(struct scale_args));
 
-	task->crystal = qargs->crystals[qargs->n_started];
+	task->crystal = qargs->crystals[qargs->n_started].cr;
+	task->refls = qargs->crystals[qargs->n_started].refls;
 
 	qargs->n_started++;
 
@@ -99,8 +101,8 @@ static void done_crystal(void *vqargs, void *task)
 }
 
 
-static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
-                          int *ninc)
+static double total_log_r(struct crystal_refls *crystals, int n_crystals,
+                          RefList *full, int *ninc)
 {
 	int i;
 	double total = 0.0;
@@ -108,8 +110,9 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
 
 	for ( i=0; i<n_crystals; i++ ) {
 		double r;
-		if ( crystal_get_user_flag(crystals[i]) ) continue;
-		r = log_residual(crystals[i], full, 0, NULL, NULL);
+		if ( crystal_get_user_flag(crystals[i].cr) ) continue;
+		r = log_residual(crystals[i].refls, crystals[i].cr,
+		                 full, 0, NULL, NULL);
 		if ( isnan(r) ) continue;
 		total += r;
 		n++;
@@ -121,7 +124,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
 
 
 /* Perform iterative scaling, all the way to convergence */
-void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
+void scale_all(struct crystal_refls *crystals, int n_crystals, int nthreads, int scaleflags)
 {
 	struct scale_args task_defaults;
 	struct scale_queue_args qargs;
@@ -163,7 +166,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
 		int i;
 		double meanB = 0.0;
 		for ( i=0; i<n_crystals; i++ ) {
-			meanB += crystal_get_Bfac(crystals[i]);
+			meanB += crystal_get_Bfac(crystals[i].cr);
 		}
 		meanB /= n_crystals;
 		STATUS("Mean B = %e\n", meanB);
@@ -181,7 +184,8 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
 
 
 /* Calculates G and B, by which cr's reflections should be multiplied to fit reference */
-int scale_one_crystal(Crystal *cr, const RefList *listR, int flags)
+int scale_one_crystal(const RefList *listS, Crystal *cr,
+                      const RefList *listR, int flags)
 {
 	const Reflection *reflS;
 	RefListIterator *iter;
@@ -202,7 +206,6 @@ int scale_one_crystal(Crystal *cr, const RefList *listR, int flags)
 	int n_part = 0;
 	int n_nom = 0;
 	int n_red = 0;
-	RefList *listS = crystal_get_reflections(cr);
 	UnitCell *cell = crystal_get_cell(cr);
 	double G, B;
 
diff --git a/src/scaling.h b/src/scaling.h
index e47a845d..27ace6b7 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -45,9 +45,10 @@ enum ScaleFlags
 	SCALE_VERBOSE_ERRORS = 1<<1,
 };
 
-extern int scale_one_crystal(Crystal *cr, const RefList *reference, int flags);
+extern int scale_one_crystal(const RefList *listS, Crystal *cr,
+                             const RefList *reference, int flags);
 
-extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
-                      int flags);
+extern void scale_all(struct crystal_refls *crystals, int n_crystals,
+                      int nthreads, int flags);
 
 #endif	/* SCALING_H */
diff --git a/src/whirligig.c b/src/whirligig.c
index 20fcfa40..aa2cf090 100644
--- a/src/whirligig.c
+++ b/src/whirligig.c
@@ -175,12 +175,11 @@ static void process_series(struct image *images, signed int *ser,
 	fprintf(fh, "%i frames in series\n\n", len);
 	fprintf(fh, "   # Serial Filename   EventID   Crystal\n");
 	for ( i=0; i<len; i++ ) {
-		Crystal *cr = images[i].crystals[ser[i]];
 		fprintf(fh, "%4i %5i %s %s %i\n", i, images[i].serial,
 		                             images[i].filename,
 		                             images[i].ev,
 		                             ser[i]);
-		p[i] = transform_reflections(crystal_get_reflections(cr),
+		p[i] = transform_reflections(images[i].crystals[ser[i]].refls,
 		                             mat[i]);
 	}
 
@@ -316,8 +315,8 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2,
 	for ( i=0; i<i1->n_crystals; i++ ) {
 	for ( j=0; j<i2->n_crystals; j++ ) {
 
-		if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]),
-		                                                      crystal_get_cell(i2->crystals[j]),
+		if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i].cr),
+		                                                      crystal_get_cell(i2->crystals[j].cr),
 		                                                      tols, &m) )
 		{
 			if ( !crystal_used(win, n1, i)
@@ -377,12 +376,12 @@ static int try_join(struct window *win, int sn)
 
 	/* Get the appropriately transformed cell from the last crystal in this
 	 * series */
-	cr = win->img[sp].crystals[win->ser[sn][sp]];
+	cr = win->img[sp].crystals[win->ser[sn][sp]].cr;
 	ref = cell_transform_intmat(crystal_get_cell(cr), win->mat[sn][sp]);
 
 	for ( j=0; j<win->img[win->join_ptr].n_crystals; j++ ) {
 		Crystal *cr2;
-		cr2 = win->img[win->join_ptr].crystals[j];
+		cr2 = win->img[win->join_ptr].crystals[j].cr;
 		if ( compare_permuted_cell_parameters_and_orientation(ref, crystal_get_cell(cr2),
 		                                                      tols,
 		                                                      &win->mat[sn][win->join_ptr]) )