Tidy up

3 files changed, 262 insertions, 244 deletions

diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
index 43013c40..abc92c82 100644
--- a/src/calibrate_detector.c
+++ b/src/calibrate_detector.c
@@ -34,7 +34,6 @@
 #include "peaks.h"
 
 
-
 static void show_help(const char *s)
 {
 	printf("Syntax: %s [options] -i <file.h5>\n\n", s);
@@ -47,21 +46,19 @@ static void show_help(const char *s)
 "  -m, --method=<method>      The calibration method.\n"
 "               xy            Determine panel shifts in plane of detector\n"
 "  -o, --output=<file>        Output results here"
-"  -n, --npeaks=<number>      Don't refine unless this many peaks observed\n"
-"                             (in the whole stream, not a single shot)\n"
+"  -n, --npeaks=<number>      Don't refine unless this many peaks are found\n"
+"                              in the whole stream.\n"
 "\n");
 }
 
-//FIXME: should this function be in detector.h?
-int calculate_projected_peak(struct panel *panel, struct rvec q,
-                             double kk, double *fs, double *ss)
-{
 
-	if (panel == NULL) return 1;
+static int calculate_projected_peak(struct panel *panel, struct rvec q,
+                                    double kk, double *fs, double *ss)
+{
+	if ( panel == NULL ) return 1;
 
 	double xd, yd, cl;
 	double plx, ply;
-
 	const double den = kk + q.w;
 
 	/* Camera length for this panel */
@@ -92,22 +89,23 @@ int calculate_projected_peak(struct panel *panel, struct rvec q,
 	if ( *ss > panel->max_ss ) return 3;
 
 	return 0;
-
 }
 
-//FIXME: should this function be in detector.h?
-struct rvec nearest_bragg(struct image image, struct rvec q)
-{
 
+static struct rvec nearest_bragg(struct image *image, struct rvec q)
+{
 	struct rvec g;
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
 	double hd, kd, ld;
 	double h, k, l;
+	int s;
+	double U[9];
+	double hvec[3];
 
-	/* miller indices of nearest Bragg reflection */
-	cell_get_cartesian(image.indexed_cell, &ax, &ay, &az,
+	/* Miller indices of nearest Bragg reflection */
+	cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
 	                                       &bx, &by, &bz,
 	                                       &cx, &cy, &cz);
 
@@ -119,56 +117,202 @@ struct rvec nearest_bragg(struct image image, struct rvec q)
 	k = lrint(kd);
 	l = lrint(ld);
 
-	/* now get scattering vector for reflectin [hkl]
-	 * this means solving the equation U*q = h */
-	double U[] = {ax, ay, az,
-	              bx, by, bz,
-	              cx, cy, cz};
-
-	double hvec[] = {h,k,l};
+	/* Now get scattering vector for reflection hkl
+	 * by solving the equation U*q = h */
+	U[0] = ax;  U[1] = ay;  U[2] = az;
+	U[3] = bx;  U[4] = by;  U[5] = bz;
+	U[6] = cx;  U[7] = cy;  U[8] = cz;
 
-	gsl_matrix_view m
-		= gsl_matrix_view_array (U, 3, 3);
+	hvec[0] = h;  hvec[1] = k;  hvec[2] = l;
 
-	gsl_vector_view b
-		= gsl_vector_view_array (hvec, 3);
+	gsl_matrix_view m = gsl_matrix_view_array(U, 3, 3);
+	gsl_vector_view b = gsl_vector_view_array(hvec, 3);
+	gsl_vector *x = gsl_vector_alloc(3);
 
-	gsl_vector *x = gsl_vector_alloc (3);
-
-	int s;
+	gsl_permutation * perm = gsl_permutation_alloc(3);
+	gsl_linalg_LU_decomp(&m.matrix, perm, &s);
+	gsl_linalg_LU_solve(&m.matrix, perm, &b.vector, x);
 
-	gsl_permutation * perm = gsl_permutation_alloc (3);
-	gsl_linalg_LU_decomp (&m.matrix, perm, &s);
-	gsl_linalg_LU_solve (&m.matrix, perm, &b.vector, x);
-
-	/* outgoing wavevector for [hkl] */
+	/* Outgoing wavevector for hkl */
 	g.u = x->data[0];
 	g.v = x->data[1];
 	g.w = x->data[2];
 
 	return g;
-
 }
 
 
-int main(int argc, char *argv[])
+static void refine_xy(FILE *fh, struct image *image, int minpeaks,
+                       const char *outfilename)
 {
+	double *weightedSumFS;
+	double *weightedSumSS;
+	double *summedWeights;
+	double *meanShiftFS;
+	double *meanShiftSS;
+	int *peaksFound;
+	double cnx, cny;
+	double xsh, ysh;
+	int pi;
+	int nChunks;
+
+	weightedSumFS = calloc(sizeof(double), image->det->n_panels);
+	weightedSumSS = calloc(sizeof(double), image->det->n_panels);
+	summedWeights = calloc(sizeof(double), image->det->n_panels);
+	peaksFound = calloc(sizeof(int), image->det->n_panels);
+	meanShiftFS = calloc(sizeof(double), image->det->n_panels);
+	meanShiftSS = calloc(sizeof(double), image->det->n_panels);
+
+	/* Initialize arrays  */
+	for (pi=0; pi<image->det->n_panels; pi++) {
+		weightedSumFS[pi] = 0;
+		weightedSumSS[pi] = 0;
+		summedWeights[pi] = 0;
+		peaksFound[pi] = 0;
+		meanShiftFS[pi] = 0;
+		meanShiftSS[pi] = 0;
+	}
+
+	fesetround(1);  /* Round towards nearest */
+
+	nChunks = 0;
+	while ( 1 )  {
+
+		int fail;
+		int nFeatures;
+		int i;
+
+		/* Get next chunk */
+		fail = read_chunk(fh, image);
+		if ( fail ) break;
+		nChunks += 1;
+
+		/* Skip if no peaks found */
+		if ( image->features == NULL ) {
+			continue;
+		}
+
+		/* Loop through peaks to determine mean panel shift */
+		nFeatures = image_feature_count(image->features);
+
+		for ( i=0; i<nFeatures; i++ ) {
+
+			struct panel *p;
+			struct imagefeature *thisFeature;
+			struct rvec q;
+			double twotheta;
+			struct rvec g;
+			double fs, ss;   /* Observed peaks */
+			double pfs, pss; /* Predicted peaks */
+			double dfs, dss; /* Observed - predicted */
+			int pi;
+			double thisWeight;
+
+			/* If we find a feature, determine peak
+			 * position */
+			thisFeature = image_get_feature(image->features, i);
+			if ( thisFeature == NULL ) {
+				continue;
+			}
+
+			fs = thisFeature->fs;
+			ss = thisFeature->ss;
+
+			p = find_panel(image->det, fs, ss);
+			if ( p == NULL ) {
+				continue;
+			}
+			if ( p->no_index ) continue;
+
+			/* Now determine the predicted peak position.
+			 * Scattering vector of this peak */
+			q = get_q(image, fs, ss, &twotheta, 1.0/image->lambda);
+
+			/* Scattering vector of nearest bragg peak */
+			g = nearest_bragg(image, q);
+
+			/* Coordinate of this predicted peak */
+			fail = calculate_projected_peak(p, g, 1.0/image->lambda,
+			                                &pfs, &pss);
+			if ( fail != 0 ) continue;
+
+			/* Finally, we have the shift for this peak */
+			dfs = pfs - fs;
+			dss = pss - ss;
+
+			/* Add this shift to the weighted sum over shifts */
+			pi = find_panel_number(image->det,fs,ss);
+			thisWeight = 1.0;
+			weightedSumFS[pi] += thisWeight*dfs;
+			weightedSumSS[pi] += thisWeight*dss;
+			summedWeights[pi] += thisWeight;
+			peaksFound[pi] += 1;
+
+		}
+
+	}
+
+	/* Calculate weighted average shift in peak positions */
+	for ( pi=0; pi<image->det->n_panels; pi++ ) {
+		meanShiftFS[pi] = weightedSumFS[pi]/summedWeights[pi];
+		meanShiftSS[pi] = weightedSumSS[pi]/summedWeights[pi];
+	}
+
+	/* Populate the image structure with new geometry info */
+	for ( pi=0; pi<image->det->n_panels; pi++ ) {
+
+		struct panel *p;
+
+		p = &image->det->panels[pi];
+
+		xsh = 0;
+		ysh = 0;
+
+		if ( peaksFound[pi] >= minpeaks ) {
+
+			/* Convert shifts from raw coords to lab frame */
+			xsh = meanShiftFS[pi]*p->fsx + meanShiftSS[pi]*p->ssx;
+			ysh = meanShiftFS[pi]*p->fsy + meanShiftSS[pi]*p->ssy;
+
+			/* Add shifts to original panel corner locations */
+			cnx = p->cnx + xsh;
+			cny = p->cny + ysh;
+
+		} else {
+
+			/* Not refined? use original coordinates */
+			cnx = p->cnx;
+			cny = p->cny;
+
+		}
+
+		image->det->panels[pi].cnx = cnx;
+		image->det->panels[pi].cny = cny;
+		if ( peaksFound[pi] < minpeaks) {
+			image->det->panels[pi].no_index = 1;
+		}
+
+		STATUS("Panel %s, # peaks: %10d, mean shifts: %f %f\n",
+		       p->name, peaksFound[pi], xsh, ysh);
+
+	}
+
+	/* Write the new geometry file */
+	write_detector_geometry(outfilename, image->det);
+}
 
 
+int main(int argc, char *argv[])
+{
 	char c;
 	struct image image;
-	struct panel p;
 	char *filename = NULL;
 	char *geometry = NULL;
 	char *method = NULL;
 	char *outputfile = NULL;
 	FILE *fh = NULL;
 	FILE *outfh = NULL;
-	int nFeatures;
-	int i;
-	int fail;
-	int nChunks;
-	int minpeaks=0;
+	int minpeaks = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -182,14 +326,13 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:g:m:o:n:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hi:g:m:o:n:",
+	                        longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' :
 			show_help(argv[0]);
 			return 0;
-		case 0 :
-			break;
 		case 'i' :
 			filename = strdup(optarg);
 			break;
@@ -205,6 +348,8 @@ int main(int argc, char *argv[])
 		case 'n' :
 			minpeaks = atoi(optarg);
 			break;
+		case 0 :
+			break;
 		default :
 			return 1;
 		}
@@ -215,12 +360,12 @@ int main(int argc, char *argv[])
 		ERROR("You must specify the input filename with -i\n");
 		return 1;
 	}
+
 	fh = fopen(filename,"r");
-   if ( fh == NULL ) {
-		ERROR("Problem opening file\n");
+	if ( fh == NULL ) {
+		ERROR("Couldn't open file '%s'\n", filename);
 		return 1;
 	}
-	printf("Read stream file: %s\n",filename);
 
 	if ( geometry == NULL ) {
 		ERROR("You need to specify a geometry file with --geometry\n");
@@ -232,193 +377,45 @@ int main(int argc, char *argv[])
 		ERROR("Failed to read detector geometry from %s\n", geometry);
 		return 1;
 	}
-	printf("Read geometry file: %s\n",geometry);
+	free(geometry);
+
 	image.width = image.det->max_fs;
 	image.height = image.det->max_ss;
-	free(geometry);
 
-	if ( !(outputfile == NULL) ) {
-		printf("Writing result to file: %s\n",outputfile);
-      outfh = fopen(outputfile,"w");
-   } else {
-		ERROR("You did not specify an output file.\n");
+	if ( outputfile != NULL ) {
+		STATUS("Writing result to '%s'\n", outputfile);
+		outfh = fopen(outputfile, "w");
+	} else {
+		ERROR("You need to specify an output file.\n");
 		return 1;
 	}
+
 	if ( minpeaks == 0 ) {
 		minpeaks = 1;
-		printf("You did not specify minimum number of peaks."
-             " Setting default value of %d\n",minpeaks);
+		STATUS("You did not specify minimum number of peaks.\n");
+		STATUS("Using default value of %d\n", minpeaks);
 	}
 
 	if ( method == NULL ) {
-		printf("You did not specify a refinement method-- setting default.\n");
+		STATUS("You did not specify a refinement method "
+		       "- using default.\n");
 		method = strdup("xy");
 	}
 
+	if ( strcmp(method,"xy") == 0 ) {
 
-	if ( !strcmp(method,"xy" ) ) {
-
-		printf("Using refinement method %s\n",method);
-
-		double * weightedSumFS, * weightedSumSS;
-		double * summedWeights;
-		double * meanShiftFS, * meanShiftSS;
-		int * peaksFound;
-		double cnx, cny;
-		double xsh, ysh;
-
-		weightedSumFS = (double *)calloc(sizeof(double), image.det->n_panels);
-		weightedSumSS = (double *)calloc(sizeof(double), image.det->n_panels);
-		summedWeights = (double *)calloc(sizeof(double), image.det->n_panels);
-		peaksFound = (int *)calloc(sizeof(int), image.det->n_panels);
-		meanShiftFS = (double *)calloc(sizeof(double), image.det->n_panels);
-		meanShiftSS = (double *)calloc(sizeof(double), image.det->n_panels);
-
-		/* initialize arrays (is there a standard function to do this?)  */
-		int pi;
-		for (pi=0; pi<image.det->n_panels; pi++) {
-			weightedSumFS[pi] = 0;
-	      weightedSumSS[pi] = 0;
-   	   summedWeights[pi] = 0;
-     		peaksFound[pi] = 0;
-      	meanShiftFS[pi] = 0;
-      	meanShiftSS[pi] = 0;
-		}
-
-
-		fesetround(1);  /* Round towards nearest */
-
-		nChunks = 0;
-		while (1) {
-
-		   // check for peaks in next file
-			fail = read_chunk(fh, &image);
-		   if ( fail == 1 ) {
-		      // FIXME:should check if this is EOF, or broken file handle
-		      break;
-		   }
-			nChunks += 1;
-
-
-			// move on to the next chunk if no peaks found
-			if ( image.features == NULL ) {
-				continue;
-			}
-
-			// now loop through peaks to determine mean panel shift
-			nFeatures = image_feature_count(image.features);
-
-			for (i=0; i<nFeatures; i++) {
-
-				struct panel * p;
-				struct imagefeature * thisFeature;
-				struct rvec q;
-				double twotheta;
-				struct rvec g;
-				double fs, ss;   /* observed peaks */
-				double pfs, pss; /* predicted peaks */
-				double dfs, dss; /* observed - predicted */
-				int pi;
-				double thisWeight;
-				//int fail;
+		STATUS("Using refinement method %s\n", method);
 
-				/* if we find a feature, determine peak position */
-				thisFeature = image_get_feature(image.features,i);
-				if ( thisFeature == NULL ) {
-					continue;
-				}
-
-				fs = thisFeature->fs;
-				ss = thisFeature->ss;
-
-				p = find_panel(image.det, fs, ss);
-				if ( p == NULL ) {
-					continue;
-				}
-				if ( p->no_index ) continue;
-
-				/* now determine the predicted peak position */
-
-				/* scattering vector of this peak */
-				q = get_q(&image, fs, ss, &twotheta, 1.0/image.lambda);
-
-				/* scattering vector of nearest bragg peak */
-				g = nearest_bragg(image, q);
-
-				/* coordinate of this predicted peak */
-				fail = calculate_projected_peak(p,g,1/image.lambda,&pfs,&pss);
-
-				/* check for error, e.g. out of panel */
-				if ( fail != 0 ) continue;
-
-				/* Finally, we have the shift in position of this peak */
-				dfs = pfs - fs;
-				dss = pss - ss;
-
-				/* Add this shift to the weighted sum over shifts */
-				pi = find_panel_number(image.det,fs,ss);
-				thisWeight = 1; // FIXME: use real weighting some day
-				weightedSumFS[pi] += thisWeight*dfs;
-				weightedSumSS[pi] += thisWeight*dss;
-				summedWeights[pi] += thisWeight;
-				peaksFound[pi] += 1;
-
-			} /* end loop through image features */
-
-		} /* end loop through stream chunks */
-
-		/* calculate weighted average shift in peak positions */
-		for (pi=0; pi < image.det->n_panels; pi++) {
-			meanShiftFS[pi] = weightedSumFS[pi]/summedWeights[pi];
-			meanShiftSS[pi] = weightedSumSS[pi]/summedWeights[pi];
-		}
-
-		/* first populate the image structure with new geometry info */
-		for (pi=0; pi < image.det->n_panels; pi++) {
-
-			p = image.det->panels[pi];
-
-			xsh = 0;
-			ysh = 0;
-
-			if ( peaksFound[pi] >= minpeaks ) {
-
-				/* convert shifts from raw coords to lab frame */
-				xsh = meanShiftFS[pi]*p.fsx + meanShiftSS[pi]*p.ssx;
-				ysh = meanShiftFS[pi]*p.fsy + meanShiftSS[pi]*p.ssy;
-
-				/* add shifts to original panel corner locations */
-				cnx = p.cnx + xsh;
-				cny = p.cny + ysh;
-
-			} else {
-				/* not refined? use original coordinates */
-				cnx = p.cnx;
-				cny = p.cny;
-			}
-
-			image.det->panels[pi].cnx = cnx;
-			image.det->panels[pi].cny = cny;
-			if ( peaksFound[pi] < minpeaks) image.det->panels[pi].no_index = 1;
-
-			printf("panel %s, # peaks: %10d, mean shifts: %f %f\n",p.name, peaksFound[pi],xsh,ysh);
-
-		}
-
-		/* write the new geometry file */
-		print_detector_geometry(outfh,&image);
+		refine_xy(fh, &image, minpeaks, outputfile);
 
 	} else {
 
 		printf("Refinement method %s not recognized\n",method);
-   	return 1;
+   		return 1;
 
 	}
 
-
 	fclose(outfh);
 
-	printf("Done!\n");
-
 	return 0;
 }
diff --git a/src/detector.c b/src/detector.c
index f705430d..da57c658 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -941,55 +941,76 @@ void get_pixel_extents(struct detector *det,
 }
 
 
-int print_detector_geometry(FILE * file, struct image * image) {
-
-	struct panel * p;
+int write_detector_geometry(const char *filename, struct detector *det)
+{
+	struct panel *p;
 	int pi;
+	FILE *fh;
 
-	if ( image == NULL ) return 1;
-	if ( file == NULL ) return 2;
-	if ( image->det->n_panels < 1 ) return 3;
-
-	for (pi=0; pi<image->det->n_panels; pi++) {
+	if ( filename == NULL ) return 2;
+	if ( det->n_panels < 1 ) return 3;
 
-		p = &(image->det->panels[pi]);
+	fh = fopen(filename, "w");
+	if ( fh == NULL ) return 1;
 
-		if ( p == NULL ) return 4;
+	for ( pi=0; pi<det->n_panels; pi++) {
 
-		fprintf(file,"%s/min_fs = %d\n",p->name,p->min_fs);
-		fprintf(file,"%s/min_ss = %d\n",p->name,p->min_ss);
-		fprintf(file,"%s/max_fs = %d\n",p->name,p->max_fs);
-		fprintf(file,"%s/max_ss = %d\n",p->name,p->max_ss);
-		fprintf(file,"%s/badrow_direction = %C\n",p->name,p->badrow);
-		fprintf(file,"%s/res = %g\n",p->name,p->res);
-		fprintf(file,"%s/peak_sep = %g\n",p->name,p->peak_sep);
-		fprintf(file,"%s/clen = %s\n",p->name,p->clen_from);
-		//FIXME: the following is sketchy, but it will work for now.  we need
-		//       to generalise the parser in detector.c
 		char coord;
 		char sign;
-		if (p->fsx != 0){
-			if (p->fsx>0){sign='+';}else{sign='-';}
+
+		p = &(det->panels[pi]);
+
+		if ( p == NULL ) return 4;
+
+		fprintf(fh, "%s/min_fs = %d\n", p->name, p->min_fs);
+		fprintf(fh, "%s/min_ss = %d\n", p->name, p->min_ss);
+		fprintf(fh, "%s/max_fs = %d\n", p->name, p->max_fs);
+		fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss);
+		fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow);
+		fprintf(fh, "%s/res = %g\n", p->name, p->res);
+		fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep);
+		fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from);
+
+		/* FIXME: The following is sketchy, but it will work for now.
+		 * We need to generalise the parser in detector.c */
+		if ( p->fsx != 0 ) {
+			if ( p->fsx > 0 ) {
+				sign='+';
+			} else {
+				sign='-';
+			}
 			coord = 'x';
 		} else {
-			if (p->fsy>0){sign='+';}else{sign='-';}
+			if ( p->fsy > 0 ) {
+				sign='+';
+			} else {
+				sign='-';
+			}
 			coord = 'y';
 		}
-		fprintf(file,"%s/fs = %C%C\n",p->name, sign, coord);
+		fprintf(fh, "%s/fs = %C%C\n", p->name, sign, coord);
 		if (p->ssx != 0){
-			if (p->ssx>0){sign='+';}else{sign='-';}
+			if ( p->ssx > 0 ) {
+				sign='+';
+			} else {
+				sign='-';
+			}
 			coord = 'x';
 		} else {
-			if (p->ssy>0){sign='+';}else{sign='-';}
+			if ( p->ssy > 0 ) {
+				sign='+';
+			} else {
+				sign='-';
+			}
 			coord = 'y';
 		}
-		fprintf(file,"%s/ss = %C%C\n",p->name, sign, coord);
-		fprintf(file,"%s/corner_x = %g\n",p->name,p->cnx);
-		fprintf(file,"%s/corner_y = %g\n",p->name,p->cny);
-		fprintf(file,"%s/no_index = %d\n",p->name,p->no_index);
+		fprintf(fh, "%s/ss = %C%C\n", p->name, sign, coord);
+		fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx);
+		fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny);
+		fprintf(fh, "%s/no_index = %d\n", p->name, p->no_index);
 
 	}
+	fclose(fh);
 
 	return 0;
-
 }
diff --git a/src/detector.h b/src/detector.h
index 4c5166a8..86e7eea3 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -115,7 +115,7 @@ extern int reverse_2d_mapping(double x, double y, double *pfs, double *pss,
 
 extern double largest_q(struct image *image);
 
-extern int print_detector_geometry(FILE * file, struct image * image);
+extern int write_detector_geometry(const char *filename, struct detector *det);