Peak search algorithms should not mutate image structure

This changes all the peak search procedures into pure functions that
return a new ImageFeatureList.  This takes the management of
image->features out of the hands of the peak search routines, and into
the calling code's responsibility.

In turn, this allows a load of stuff to become const.

2 files changed, 89 insertions, 76 deletions

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/process_image.c b/src/process_image.c
index 782c97b9..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,