Use unsigned types for image storage

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

12 files changed, 75 insertions, 60 deletions

diff --git a/src/control.c b/src/control.c
index 1b4e826..c79ca17 100644
--- a/src/control.c
+++ b/src/control.c
@@ -13,7 +13,7 @@
 
 #include "control.h"
 
-int control_add_image(ControlContext *ctx, int16_t *image, int width, int height, double tilt) {
+int control_add_image(ControlContext *ctx, uint16_t *image, int width, int height, double tilt) {
 
 	ctx->images[ctx->n_images].tilt = tilt;
 	ctx->images[ctx->n_images].omega = ctx->omega;
diff --git a/src/control.h b/src/control.h
index 5de0b86..2e16079 100644
--- a/src/control.h
+++ b/src/control.h
@@ -49,7 +49,7 @@ typedef enum {
 
 typedef struct imagerecord_struct {
 
-	int16_t		*image;
+	uint16_t		*image;
 	double		tilt;
 	double		omega;
 	
@@ -115,7 +115,7 @@ typedef struct cctx_struct {
 
 } ControlContext;
 
-extern int control_add_image(ControlContext *ctx, int16_t *image, int width, int height, double tilt);
+extern int control_add_image(ControlContext *ctx, uint16_t *image, int width, int height, double tilt);
 
 #endif	/* CONTROL_H */
 
diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 3a87aa5..d60b789 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -72,11 +72,11 @@ static gboolean imagedisplay_configure_event(GtkWidget *widget, GdkEventConfigur
 
 }
 
-static void imagedisplay_put_data(ImageDisplay *imagedisplay, int16_t *image16) {
+static void imagedisplay_put_data(ImageDisplay *imagedisplay, uint16_t *image16) {
 
 	unsigned int x, y;
 	unsigned int w, h;
-	int16_t min, max;
+	int min, max;
 	
 	h = imagedisplay->height;
 	w = imagedisplay->width;
@@ -84,19 +84,19 @@ static void imagedisplay_put_data(ImageDisplay *imagedisplay, int16_t *image16)
 	min = 0;  max = 0;
 	for ( y=0; y<h; y++ ) {
 		for ( x=0; x<w; x++ ) {
-			int16_t val;
+			uint16_t val;
 			val = image16[x+w*y];
 			if ( val > max ) max = val;
 			if ( val < min ) min = val;
 		}
 	}
-	
+	printf("min=%i, max=%i\n", min, max);
 	/* Turn 16-bit image data into 8-bit display data */
 	imagedisplay->data = malloc(3*w*h);
 	for ( y=0; y<h; y++ ) {
 		for ( x=0; x<w; x++ ) {
 		
-			int16_t val16, val8;
+			uint16_t val16, val8;
 			
 			val16 = image16[x+w*y];
 			val8 = (255*(val16-min)) / (max-min);
@@ -121,7 +121,7 @@ static void imagedisplay_close(GtkWidget *widget, ImageDisplay *imagedisplay) {
 }
 
 /* Display an image */
-ImageDisplay *imagedisplay_open(int16_t *image, unsigned int width, unsigned int height, const char *title) {
+ImageDisplay *imagedisplay_open(uint16_t *image, unsigned int width, unsigned int height, const char *title) {
 
 	ImageDisplay *imagedisplay;
 	GdkGeometry geom;
diff --git a/src/imagedisplay.h b/src/imagedisplay.h
index 61aa157..94e6c13 100644
--- a/src/imagedisplay.h
+++ b/src/imagedisplay.h
@@ -37,7 +37,7 @@ typedef struct struct_imagedisplay {
 	
 } ImageDisplay;
 
-extern ImageDisplay *imagedisplay_open(int16_t *image, unsigned int width, unsigned int height, const char *title);
+extern ImageDisplay *imagedisplay_open(uint16_t *image, unsigned int width, unsigned int height, const char *title);
 extern void imagedisplay_mark_point(ImageDisplay *imagedisplay, unsigned int x, unsigned int y);
 extern void imagedisplay_mark_circle(ImageDisplay *imagedisplay, unsigned int x, unsigned int y);
 #include "control.h"
diff --git a/src/itrans-lsq.c b/src/itrans-lsq.c
index 841a8b9..168a5b3 100644
--- a/src/itrans-lsq.c
+++ b/src/itrans-lsq.c
@@ -116,7 +116,7 @@ static int itrans_peaksearch_lsq_fdf(const gsl_vector *gaussian, void *params, g
 
 }
 
-static void itrans_interpolate(int16_t *image, int width, int x, int y) {
+static void itrans_interpolate(uint16_t *image, int width, int x, int y) {
 
 	int a, b, c, d;
 	double av_horiz, av_vert, av;
@@ -131,9 +131,9 @@ static void itrans_interpolate(int16_t *image, int width, int x, int y) {
 
 unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
 
-	int16_t max_val = 0;
+	uint16_t max_val = 0;
 	int width, height;
-	int16_t *image;
+	uint16_t *image;
 	unsigned int n_reflections = 0;
 	
 	width = imagerecord.width;
@@ -158,7 +158,7 @@ unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx,
 		gboolean sanity;
 		int dx, dy;
 		int bb_lh, bb_rh, bb_top, bb_bot;
-		int16_t ival;
+		uint16_t ival;
 		double sx, sy;
 	
 		/* Locate the highest point */
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index 6d64a08..f0677d4 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -29,7 +29,7 @@
  * can be done as matrices to save effort
  * Renormalises matrix to 0->1
  */
-static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(int16_t *image, int width, int height) {
+static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(uint16_t *image, int width, int height) {
 
 	gsl_matrix *raw;
 	int i, j;
@@ -496,7 +496,7 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx
 	int i;
 	double px,py;
 	double tilt_degrees = imagerecord.tilt;
-	int16_t *image = imagerecord.image;
+	uint16_t *image = imagerecord.image;
 
 	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord.width, imagerecord.height);
 	p = itrans_peaksearch_stat_iterate(m, &n_reflections);
diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c
index 6aadea2..1aee5ca 100644
--- a/src/itrans-threshold.c
+++ b/src/itrans-threshold.c
@@ -24,12 +24,12 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext
 	int width, height;
 	int x, y;
 	unsigned int n_reflections = 0;
-	int16_t *image = imagerecord.image;
+	uint16_t *image = imagerecord.image;
 	double tilt_degrees = imagerecord.tilt;
 	
 	width = imagerecord.width;
 	height = imagerecord.height;
-
+	
 	/* Simple Thresholding */
 	for ( y=0; y<height; y++ ) {
 		for ( x=0; x<width; x++ ) {
@@ -53,10 +53,10 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext
 
 unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, ControlContext *ctx) {
 
-	int16_t max_val = 0;
+	uint16_t max_val = 0;
 	int width, height;
 	unsigned int n_reflections = 0;
-	int16_t *image = imagerecord.image;
+	uint16_t *image = imagerecord.image;
 	double tilt_degrees = imagerecord.tilt;
 	
 	width = imagerecord.width;
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index a096e74..1454fd4 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -27,7 +27,7 @@ unsigned int itrans_peaksearch_zaefferer(ImageRecord imagerecord, ControlContext
 	int x, y;
 	unsigned int n_reflections;
 	int width, height;
-	int16_t *image;
+	uint16_t *image;
 	double tilt_degrees;
 	
 	tilt_degrees = imagerecord.tilt;
diff --git a/src/mapping.c b/src/mapping.c
index 24a6470..6d74ddd 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -21,20 +21,20 @@
 ReflectionContext *mapping_create(ControlContext *ctx) {
 
 	ReflectionContext *rctx;
-	ImageDisplay *sum_id;
-	int16_t max_val;
-	int max_x, max_y;
-	int twidth, theight;
-	int x, y, i;
-	int16_t *image_total;
-	double x_centre = 0;
-	double y_centre = 0;
+	int i;
 	
 	/* Create reflection context */
 	rctx = reflection_init();
 	
 	if ( !ctx->prealign ) {
 	
+		int max_x, max_y;
+		uint16_t max_val;
+		int twidth, theight;
+		int x, y;
+		uint16_t *image_total;
+		ImageDisplay *sum_id;
+		
 		/* Find centre of image stack
 		 * Determine maximum size of image to accommodate, and allocate memory */
 		twidth = 0;  theight = 0;
@@ -42,8 +42,8 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 			if ( ctx->images[i].width > twidth ) twidth = ctx->images[i].width;
 			if ( ctx->images[i].height > theight ) theight = ctx->images[i].height;
 		}
-		image_total = malloc(twidth * theight * sizeof(int16_t));
-		memset(image_total, 0, twidth * theight * sizeof(int16_t));
+		image_total = malloc(twidth * theight * sizeof(uint16_t));
+		memset(image_total, 0, twidth * theight * sizeof(uint16_t));
 		
 		/* Add the image stack together */
 		for ( i=0; i<ctx->n_images; i++ ) {
@@ -67,13 +67,17 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 				}
 			}
 		}
-		x_centre = max_x;
-		y_centre = max_y;
+		
+		/* Record this measurement on all images */
+		for ( i=0; i<ctx->n_images; i++ ) {
+			ctx->images[i].x_centre = max_x;
+			ctx->images[i].y_centre = max_y;
+		}
 		
 		/* Display */
 		sum_id = imagedisplay_open(image_total, twidth, theight, "Sum of All Images");
-		imagedisplay_mark_point(sum_id, x_centre, y_centre);
-		imagedisplay_add_tilt_axis(sum_id, x_centre, y_centre, ctx->omega);
+		imagedisplay_mark_point(sum_id, max_x, max_y);
+		imagedisplay_add_tilt_axis(sum_id, max_x, max_y, ctx->omega);
 		
 	}
 	
@@ -81,10 +85,6 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 	 *   (let itrans add the reflections to rctx for now) */
 	ctx->reflectionctx = rctx;
 	for ( i=0; i<ctx->n_images; i++ ) {
-		if ( !ctx->prealign ) {
-			ctx->images[i].x_centre = x_centre;
-			ctx->images[i].y_centre = y_centre;
-		}
 		itrans_process_image(ctx->images[i], ctx);
 	}
 	
diff --git a/src/mrc.c b/src/mrc.c
index 50596a5..1703b7b 100644
--- a/src/mrc.c
+++ b/src/mrc.c
@@ -33,6 +33,7 @@ int mrc_read(ControlContext *ctx) {
 	unsigned int i;
 	unsigned int extsize;
 	double pixel_size;
+	int x, y;
 	
 	fh = fopen(ctx->filename, "rb");
 
@@ -40,7 +41,7 @@ int mrc_read(ControlContext *ctx) {
 	fread(&mrc, sizeof(MRCHeader), 1, fh);
 	printf("%i images in series\n", mrc.nz);
 	if ( mrc.mode != 1 ) {
-		fprintf(stderr, "MR: Unknown MRC image mode\n");
+		fprintf(stderr, "MR: Unknown MRC image mode (must be 16-bit signed)\n");
 		fclose(fh);
 		return -1;
 	}
@@ -67,7 +68,8 @@ int mrc_read(ControlContext *ctx) {
 	
 	for ( i=0; i<mrc.nz; i++ ) {
 	
-		int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
+		int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(uint16_t));
+		uint16_t *uimage = malloc(mrc.ny * mrc.nx * sizeof(uint16_t));
 		
 		printf("Image #%3i: tilt=%f omega=%f L=%f\n", i, ext[i].a_tilt, ext[i].tilt_axis, ext[i].magnification);
 		ctx->camera_length = ext[i].magnification;
@@ -82,7 +84,15 @@ int mrc_read(ControlContext *ctx) {
 		fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET);
 		fread(image, mrc.nx*mrc.ny*2, 1, fh);
 		
-		control_add_image(ctx, image, mrc.nx, mrc.ny, ext[i].a_tilt);
+		for ( x=0; x<mrc.nx; x++ ) {
+			for ( y=0; y<mrc.ny; y++ ) {
+				uimage[x + mrc.nx*y] = image[x + mrc.nx*y] + 32767;
+				if ( uimage[x + mrc.nx*y] < 10 ) printf("%i\n", uimage[x + mrc.nx*y]);
+			}
+		}
+		free(image);
+		
+		control_add_image(ctx, uimage, mrc.nx, mrc.ny, ext[i].a_tilt);
 		
 	}
 
diff --git a/src/mrc.h b/src/mrc.h
index 16c7d9d..6414dfe 100644
--- a/src/mrc.h
+++ b/src/mrc.h
@@ -45,27 +45,27 @@ typedef struct struct_mrcheader {
 	float		amin;
 	float		amax;
 	float		amean;
-	int16_t		ispg;		/* Space group number */
-	int16_t		nsymbt;
+	uint16_t		ispg;		/* Space group number */
+	uint16_t		nsymbt;
 	int32_t		next;
-	int16_t		dvid;
+	uint16_t		dvid;
 	char		extra[30];
-	int16_t		numintegers;
-	int16_t		numfloats;
-	int16_t		sub;
-	int16_t		zfac;
+	uint16_t		numintegers;
+	uint16_t		numfloats;
+	uint16_t		sub;
+	uint16_t		zfac;
 	float		min2;
 	float		max2;
 	float		min3;	
 	float		max3;	
 	float		min4;	
 	float		max4;
-	int16_t		idtype;
-	int16_t		lens;
-	int16_t		nd1;
-	int16_t		nd2;
-	int16_t		vd1;
-	int16_t		vd2;
+	uint16_t		idtype;
+	uint16_t		lens;
+	uint16_t		nd1;
+	uint16_t		nd2;
+	uint16_t		vd1;
+	uint16_t		vd2;
 	float		tiltangles[9];
 	float		zorg;
 	float		xorg;
diff --git a/src/readpng.c b/src/readpng.c
index dc72e6e..0ea87fe 100644
--- a/src/readpng.c
+++ b/src/readpng.c
@@ -37,7 +37,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	png_bytep *row_pointers;
 	unsigned int x;
 	unsigned int y;
-	int16_t *image;
+	uint16_t *image;
 	
 	/* Open file */
 	fh = fopen(filename, "rb");
@@ -97,7 +97,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	height = png_get_image_height(png_ptr, info_ptr);
 	bit_depth = png_get_bit_depth(png_ptr, info_ptr);
 	channels = png_get_channels(png_ptr, info_ptr);
-	printf("RI: width=%i, height=%i, depth=%i, channels=%i\n", width, height, bit_depth, channels);
+	//printf("RI: width=%i, height=%i, depth=%i, channels=%i\n", width, height, bit_depth, channels);
 	if ( (bit_depth != 16) && (bit_depth != 8) ) {
 		printf("RI: Whoops! Can't handle images with other than 8 or 16 bpp yet...\n");
 		png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
@@ -108,7 +108,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	/* Get image data */
 	row_pointers = png_get_rows(png_ptr, info_ptr);
 	
-	image = malloc(height * width * sizeof(int16_t));
+	image = malloc(height * width * sizeof(uint16_t));
 	ctx->fmode = FORMULATION_CLEN;
 	
 	for ( y=0; y<height; y++ ) {
@@ -120,10 +120,12 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 				int i;
 				val = 0;
 				for ( i=0; i<channels; i++ ) {
-					/* PNG files are big-endian: */
+					/* PNG files are big-endian... */
 					val += row_pointers[y][(channels*x*2)+(2*i)] << 8;
 					val += row_pointers[y][(channels*x*2)+(2*i)+1];
 				}
+				val /= channels;
+				if ( val > 65535 ) printf("%i\n", val);
 			}
 			if ( bit_depth == 8 ) {
 				int i;
@@ -131,6 +133,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 				for ( i=0; i<channels; i++ ) {
 					val += row_pointers[y][(channels*x)+i];
 				}
+				val /= channels;
 			}
 			
 			image[x + width*(height-1-y)] = val;
@@ -141,7 +144,9 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
 	fclose(fh);
 
-	control_add_image(ctx, image, width, height, tilt);
+	if ( control_add_image(ctx, image, width, height, tilt) == 0 ) {
+		imagedisplay_open(image, width, height, "First image");
+	}
 	
 	return 0;