New cache format stores image features rather than 3D reflection coordinates

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@167 bf6ca9ba-c028-0410-8290-897cf20841d1

10 files changed, 168 insertions, 110 deletions

diff --git a/src/cache.c b/src/cache.c
index 8c328b3..8f05e05 100644
--- a/src/cache.c
+++ b/src/cache.c
@@ -1,7 +1,7 @@
 /*
  * cache.c
  *
- * Save the reflection datablock to save having to recalculate it
+ * Save a list of features from images to save recalculation
  *
  * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
  *	    Thomas White <taw27@cam.ac.uk>
@@ -18,103 +18,168 @@
 #include <stdio.h>
 #include <stdint.h>
 #include <string.h>
+#include <assert.h>
 
-#include "reflections.h"
-#include "cache.h"
+#include "image.h"
 
-typedef struct struct_cacheheader {
-	char top[16];
-	int count;
-	double scale;
+typedef struct {
+	char		top[16];
+	uint32_t	n_images;
 } CacheHeader;
 
-ReflectionList *cache_load(const char *filename) {
 
-	FILE *f;
-	CacheHeader ch;
-	ReflectionList *reflectionlist;
-	size_t cachedreflection_size;
-	int i;
-	
-	cachedreflection_size = sizeof(Reflection) - sizeof(Reflection *);
+int cache_save(ImageList *images, char *cache_filename) {
 
-	reflectionlist = reflectionlist_new();
-	f = fopen(filename, "rb");
-	if ( !f ) {
-		fprintf(stderr, "Couldn't open cache file\n");
-	}
-	if ( fread(&ch, sizeof(CacheHeader), 1, f) == 0 ) {
-		fprintf(stderr, "Couldn't read cache header\n");
-		fclose(f);
-		return NULL;
+	FILE		*fh;
+	CacheHeader	ch;
+	uint32_t	*index;
+	size_t		index_size;
+	int		i;
+	
+	fh = fopen(cache_filename, "wb");
+	if ( !fh ) {
+		fprintf(stderr, "Couldn't open feature cache file.\n");
+		return -1;
 	}
 	
-	for ( i=0; i<ch.count; i++ ) {
+	memcpy(&ch.top, "DTR-1.0.6\0\0\0\0\0\0\0", 16);
+	ch.n_images = images->n_images;
+	fwrite(&ch, sizeof(CacheHeader), 1, fh);
 	
-		Reflection *cr;
+	index_size = images->n_images * 2 * sizeof(uint32_t);
+	index = malloc(index_size);
+	/* Dummy write to reserve space in the file */
+	memset(index, 0, index_size);
+	if ( fwrite(index, index_size, 1, fh) != 1 ) {
+		fprintf(stderr, "Couldn't perform dummy index write.\n");
+		fclose(fh);
+		return -1;
+	}
+	
+	for ( i=0; i<images->n_images; i++ ) {
+		
+		uint32_t size;
+		
+		size = images->images[i].features->n_features * sizeof(ImageFeature);
+		index[2*i] = ftell(fh);
+		index[2*i + 1] = images->images[i].features->n_features;
 		
-		cr = malloc(sizeof(Reflection));
-		if ( fread(cr, cachedreflection_size, 1, f) == 0 ) {
-			fprintf(stderr, "Couldn't read reflections from cache\n");
-			fclose(f);
-			free(cr);
-			reflectionlist_clear(reflectionlist);
-			return NULL;
+		if ( fwrite(images->images[i].features->features, size, 1, fh) != 1 ) {
+			fprintf(stderr, "Couldn't write feature list for image %i.\n", i);
+			fclose(fh);
+			return -1;
 		}
-			
-		cr->next = NULL;	/* Guarantee swift failure in the event of a screw-up */
-		//printf("reading (%f,%f,%f) i=%f (%d,%d,%d) %d\n",cr->x,cr->y,cr->z,cr->intensity,cr->h,cr->k,cr->l,cr->type);
-		reflection_add_from_reflection(reflectionlist, cr);
 		
 	}
 	
-	fclose(f);
+	/* Actual write */
+	if ( fseek(fh, sizeof(CacheHeader), SEEK_SET) ) {
+		fprintf(stderr, "Couldn't seek for index write.\n");
+		fclose(fh);
+		return -1;
+	}
+	if ( fwrite(index, index_size, 1, fh) != 1 ) {
+		fprintf(stderr, "Couldn't write index write.\n");
+		fclose(fh);
+		return -1;
+	}
+	
+	fclose(fh);
+	
+	return 0;
 	
-	return reflectionlist;
 }
 
-int cache_save(ReflectionList *reflectionlist, char *cache_filename) {
+int cache_load(ImageList *images, const char *filename) {
 
-	FILE *f;
-	CacheHeader ch;
-	Reflection *r;
-	int count;
-	const char top[16] = "DTRCACHE\0\0\0\0\0\0\0\0";
-	size_t cachedreflection_size;
-			
-	printf("Caching reflections to %s\n", cache_filename);
-	
-	cachedreflection_size = sizeof(Reflection) - sizeof(Reflection *);
-	
-	count = 0;
-	r = reflectionlist->reflections;
-	while ( r != NULL ) {
-		count++;
-		r = r->next;
-	};
-	
-	f = fopen(cache_filename, "wb");
-	
-	if ( f == NULL ) {
-		printf("Couldn't save reflection cache\n");
+	FILE		*fh;
+	CacheHeader	ch;
+	uint32_t	*index;
+	size_t		index_size;
+	int		i;
+	
+	fh = fopen(filename, "rb");
+	if ( !fh ) {
+		fprintf(stderr, "Couldn't open cache file.\n");
 		return -1;
 	}
-	memcpy(&ch.top, &top, sizeof(top));
-	ch.count = count;
-	ch.scale = 0.; //temp, currently doesn't do anything
-	fwrite(&ch, sizeof(CacheHeader), 1, f);
 	
-	r = reflectionlist->reflections;
-	while ( r != NULL ) {
+	if ( fread(&ch, sizeof(CacheHeader), 1, fh) != 1 ) {
+		fprintf(stderr, "Couldn't read cache header.\n");
+		fclose(fh);
+		return -1;
+	}
+	
+	/* Format check */
+	if ( strncmp(ch.top, "DTR-1.0.6\0", 8) != 0 ) {
+		fprintf(stderr, "Can't read this cache format.\n");
+		fclose(fh);
+		return -1;
+	}
+	
+	/* Muppet check */
+	if ( ch.n_images != images->n_images ) {
+		fprintf(stderr, "Number of images in cache doesn't match.\n");
+		fclose(fh);
+		return -1;
+	}
+	
+	index_size = images->n_images * 2 * sizeof(uint32_t);
+	index = malloc(index_size);
+	if ( fread(index, index_size, 1, fh) != 1 ) {
+		fprintf(stderr, "Couldn't read index.\n");
+		fclose(fh);
+		return -1;
+	}
+	
+	for ( i=0; i<images->n_images; i++ ) {
+	
+		uint32_t size;
+		uint32_t offs;
+		int j;
 		
-		fwrite(r, cachedreflection_size, 1, f);	/* Write the reflection block, stopping just short of the "next" pointer */
-		r = r->next;
+		offs = index[2*i];
+		size = index[2*i + 1];
 		
-	};
+		if ( fseek(fh, offs, SEEK_SET) ) {
+			fprintf(stderr, "Couldn't seek to feature list for image %i.\n", i);
+			fclose(fh);
+			return -1;
+		}
+		
+		if ( images->images[i].features ) image_feature_list_free(images->images[i].features);
+		images->images[i].features = image_feature_list_new();
+		if ( !images->images[i].features ) {
+			fprintf(stderr, "Couldn't allocate feature list for image %i.\n", i);
+			fclose(fh);
+			return -1;
+		}
+		
+		assert(images->images[i].features->features == NULL);
+		images->images[i].features->features = malloc(size*sizeof(ImageFeature));
+		if ( !images->images[i].features->features ) {
+			fprintf(stderr, "Couldn't allocate feature list block for image %i.\n", i);
+			fclose(fh);
+			return -1;
+		}
+		
+		if ( fread(images->images[i].features->features, size*sizeof(ImageFeature), 1, fh) != 1 ) {
+			fprintf(stderr, "Couldn't read feature list for image %i.\n", i);
+			fclose(fh);
+			return -1;
+		}
+		
+		images->images[i].features->n_features = size;
+		
+		/* Set "parent" fields for all the features */
+		for ( j=0; j<images->images[i].features->n_features; j++ ) {
+			images->images[i].features->features[j].parent = &images->images[i];
+		}
+		
+	}
 	
-	fclose(f);
+	fclose(fh);
 	
 	return 0;
-	
 }
 
diff --git a/src/cache.h b/src/cache.h
index 04dbdb9..2f1908a 100644
--- a/src/cache.h
+++ b/src/cache.h
@@ -1,7 +1,7 @@
 /*
  * cache.c
  *
- * Save the reflection datablock to save having to recalculate it
+ * Save a list of features from images to save recalculation
  *
  * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
  *	    Thomas White <taw27@cam.ac.uk>
@@ -17,8 +17,10 @@
 #include <config.h>
 #endif
 
-extern ReflectionList *cache_load(const char *filename);
-extern int cache_save(ReflectionList *reflectionlist, char *cache_filename);
+#include "image.h"
 
-#endif /*CACHE_H_*/
+extern int cache_load(ImageList *images, const char *filename);
+extern int cache_save(ImageList *images, const char *filename);
+
+#endif /*CACHE_H */
 
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 5a5f25a..fa3e5ba 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -884,7 +884,9 @@ static gint displaywindow_savecache_response(GtkWidget *widget, gint response, C
 	if ( response == GTK_RESPONSE_ACCEPT ) {
 		char *cache_filename;
 		cache_filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(widget));
-		cache_save(ctx->reflectionlist, cache_filename);
+		if ( cache_save(ctx->images, cache_filename) ) {
+			displaywindow_error("Failed to save cache file.", ctx->dw);
+		}
 		g_free(cache_filename);
 	}
 	
@@ -979,7 +981,7 @@ static void displaywindow_addmenubar(DisplayWindow *dw) {
 		{ "DirAxAction", GTK_STOCK_EXECUTE, "Start _DirAx", "<Ctrl>D", NULL, G_CALLBACK(displaywindow_dirax) },
 		{ "StopDirAxAction", GTK_STOCK_CLOSE, "Stop DirAx", NULL, NULL, G_CALLBACK(displaywindow_dirax_stop) },
 		{ "ReprojectAction", NULL, "_Reproject Diffraction Patterns", NULL, NULL, G_CALLBACK(displaywindow_reproject) },
-		{ "SaveCacheAction", NULL, "Save Reflections to _Cache", NULL, NULL, G_CALLBACK(displaywindow_savecache) },
+		{ "SaveCacheAction", NULL, "Save Image Analysis to _Cache", NULL, NULL, G_CALLBACK(displaywindow_savecache) },
 		{ "SetAxisAction", NULL, "Set Tilt Axis Position", NULL, NULL, G_CALLBACK(displaywindow_setaxis) },
 		
 		{ "HelpAction", NULL, "_Help", NULL, NULL, NULL },
diff --git a/src/image.c b/src/image.c
index 7205e0b..83ae271 100644
--- a/src/image.c
+++ b/src/image.c
@@ -77,7 +77,6 @@ void image_add_feature(ImageFeatureList *flist, double x, double y, ImageRecord
 	flist->features[flist->n_features].x = x;
 	flist->features[flist->n_features].y = y;
 	flist->features[flist->n_features].intensity = intensity;
-	flist->features[flist->n_features].x = x;
 	flist->features[flist->n_features].parent = parent;
 	
 	flist->n_features++;
@@ -104,3 +103,4 @@ void image_feature_list_free(ImageFeatureList *flist) {
 	free(flist);
 
 }
+
diff --git a/src/image.h b/src/image.h
index 642cc3a..5061893 100644
--- a/src/image.h
+++ b/src/image.h
@@ -21,8 +21,8 @@
 typedef struct imagefeature_struct {
 
 	struct imagerecord_struct	*parent;
-	int				x;
-	int				y;
+	double				x;
+	double				y;
 	double				intensity;
 	
 } ImageFeature;
diff --git a/src/itrans.h b/src/itrans.h
index ec2e852..838dd3d 100644
--- a/src/itrans.h
+++ b/src/itrans.h
@@ -17,6 +17,7 @@
 #endif
 
 #include "control.h"
+#include "image.h"
 
 extern ImageFeatureList *itrans_process_image(ImageRecord *image, PeakSearchMode psmode);
 
diff --git a/src/main.c b/src/main.c
index 0a8362c..2d8b8ed 100644
--- a/src/main.c
+++ b/src/main.c
@@ -36,14 +36,10 @@
 
 void main_do_reconstruction(ControlContext *ctx) {
 
-	if ( ctx->cache_filename ) {
-		prealign_sum_stack(ctx->images, ctx->have_centres);
-		ctx->reflectionlist = cache_load(ctx->cache_filename);
-		printf("Loading cached reflections from '%s'\n", ctx->cache_filename);
-	} else if ( ctx->inputfiletype != INPUT_DRX ) {
+	if ( ctx->inputfiletype != INPUT_DRX ) {
 		prealign_sum_stack(ctx->images, ctx->have_centres);
 		mapping_create(ctx);
-	}
+	} /* else has already been created by dirax_load() */
 	
 	if ( ctx->reflectionlist ) {
 		ctx->dw = displaywindow_open(ctx);
@@ -151,7 +147,7 @@ void main_method_dialog_open(ControlContext *ctx) {
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Zaefferer Gradient Search");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Iterative Statistical Analysis");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Get From Cache File");
-	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 3);
+	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 2);
 	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 1, 2);
 	g_signal_connect(G_OBJECT(ctx->combo_peaksearch), "changed", G_CALLBACK(main_peaksearch_changed), ctx);
 	
diff --git a/src/mapping.c b/src/mapping.c
index 544a9ac..def29c4 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -19,6 +19,7 @@
 #include "itrans.h"
 #include "image.h"
 #include "displaywindow.h"
+#include "cache.h"
 
 static int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
 
@@ -110,7 +111,8 @@ static void mapping_map_features(ControlContext *ctx) {
 			double nx, ny, nz, twotheta;
 			
 			if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j], &nx, &ny, &nz, &twotheta) ) {
-				reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
+				reflection_add(ctx->reflectionlist, nx, ny, nz,
+					ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
 			}
 
 		}
@@ -124,12 +126,17 @@ ReflectionList *mapping_create(ControlContext *ctx) {
 
 	int i;
 	
-	/* Find all the features */
-	printf("MP: Analysing images..."); fflush(stdout);
-	for ( i=0; i<ctx->images->n_images; i++ ) {
-		ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
+	if ( !ctx->cache_filename ) {
+		/* Find all the features */
+		printf("MP: Analysing images..."); fflush(stdout);
+		for ( i=0; i<ctx->images->n_images; i++ ) {
+			ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
+		}
+		printf("done.\n");
+	} else {
+		printf("MP: Loading previous image analysis from '%s'\n", ctx->cache_filename);
+		if ( cache_load(ctx->images, ctx->cache_filename) ) return NULL;
 	}
-	printf("done.\n");
 	
 	mapping_map_features(ctx);
 	
diff --git a/src/reflections.c b/src/reflections.c
index 3e37a14..9f3f1c8 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -126,18 +126,6 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d
 	
 }
 
-void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r) {
-	if ( reflectionlist->last_reflection ) {
-		reflectionlist->last_reflection->next = r;
-		reflectionlist->last_reflection = r;
-	} else {
-		reflectionlist->reflections = r;
-		reflectionlist->last_reflection = r;
-	}
-	r->next = NULL;
-	reflectionlist->n_reflections++;
-}
-
 double reflectionlist_largest_g(ReflectionList *reflectionlist) {
 
 	double max = 0.0;
diff --git a/src/reflections.h b/src/reflections.h
index 2cc7f10..7f287a1 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -61,10 +61,7 @@ extern void reflectionlist_clear(ReflectionList *reflectionlist);
 extern void reflectionlist_clear_markers(ReflectionList *reflectionlist);
 extern void reflectionlist_free(ReflectionList *reflectionlist);
 
-#include "control.h"
-#include "image.h"
 extern Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, double z, double intensity, ReflectionType type);
-extern void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r);
 
 extern double reflectionlist_largest_g(ReflectionList *reflectionlist);
 extern Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x, double y, double z);