Tidy up (mostly powder_plot)

3 files changed, 549 insertions, 392 deletions

diff --git a/Makefile.am b/Makefile.am
index 00c9fd8b..db70f378 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -119,8 +119,6 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
                           src/symmetry.c src/post-refinement.c \
                           src/hrs-scaling.c src/reflist.c
 
-src_stream_SOURCES = src/reflist-utils.c
-
 if BUILD_CUBEIT
 src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \
                      src/detector.c src/render.c src/filters.c src/image.c \
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 4985ccbe..569a41bd 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -43,7 +43,6 @@
 #include "reflist-utils.h"
 
 
-
 /* Write statistics at APPROXIMATELY this interval */
 #define STATS_EVERY_N_SECONDS (5)
 
diff --git a/src/powder_plot.c b/src/powder_plot.c
index 515c8c4e..13875364 100644
--- a/src/powder_plot.c
+++ b/src/powder_plot.c
@@ -4,7 +4,7 @@
  * Plot powder patterns
  *
  * (c) 2011 Andrew Aquila <andrew.aquila@cfel.de>
- * (c) 2006-2010 Thomas White  <taw@physics.org>
+ * (c) 2006-2011 Thomas White  <taw@physics.org>
  *
  * Part of CrystFEL - crystallography with a FEL
  *
@@ -33,8 +33,9 @@
 #include "reflist-utils.h"
 #include "symmetry.h"
 
+
 struct bin_stats {
-	unsigned int N;	
+	unsigned int N;
 	double total;
 	double mean;
 	double std_dev;
@@ -72,10 +73,11 @@ enum {
 	VOLUME
 };
 
-static int find_q_bin_index(double q, struct histogram_info *info, 
-                               struct bin_stats *hist)
+
+static int find_q_bin_index(double q, struct histogram_info *info,
+                            struct bin_stats *hist)
 {
-	/* bisection search alg. find q_bin index of order Log(n) time */	
+	/* bisection search alg. find q_bin index of order Log(n) time */
 	int mid;
 	int min = 0;
 	int max = (*info).histsize-1;
@@ -94,211 +96,263 @@ static int find_q_bin_index(double q, struct histogram_info *info,
 	return mid;
 }
 
-static int add_peak_to_histogram(double fs, double ss, double intensity, 
-                           struct image *image, struct histogram_info *info, 
-                           struct bin_stats *hist)
+
+/* Used for HDF5 files, peak list and stream positions */
+static int add_peak_to_histogram(double fs, double ss, double intensity,
+                                 struct image *image,
+                                 struct histogram_info *info,
+                                 struct bin_stats *hist)
 {
-	/* used for h5 files, peak/feature list and stream positions*/	
-	/* See Knuth TAOCP vol 2, 3rd ed, pg 232 for running variance*/	
 	struct rvec r;
 	double q, delta;
 	int i;
 
 	r = get_q(image, fs, ss, NULL, 1.0/ (*image).lambda);
 	q = modulus(r.u, r.v, r.w);
-	
-	/* ignore q value if outside of range */
-	if ((q<(*info).q_min) || (q>(*info).q_max)) {
+
+	/* Ignore q value if outside of range */
+	if ( (q<(*info).q_min) || (q>(*info).q_max) ) {
 		return 1;
-	}	
+	}
 	i = find_q_bin_index(q, info, hist);
-	
+
+	/* See Knuth TAOCP vol 2, 3rd ed, pg 232 for running variance */
 	delta = intensity - hist[i].mean;
 	hist[i].N++;
-	hist[i].total += intensity;	
+	hist[i].total += intensity;
 	hist[i].mean = hist[i].mean + delta /hist[i].N;
-	hist[i].std_dev = hist[i].std_dev + 
-	                      (delta *(intensity - hist[i].mean));
-	return 0;	
+	hist[i].std_dev = hist[i].std_dev + (delta *(intensity - hist[i].mean));
+
+	return 0;
 }
 
+
+/* Used for d and hkl of stream files where redundancy = 1 */
 static int add_d_to_histogram(double q, double intensity,
-                        struct histogram_info *info, struct bin_stats *hist)
+                              struct histogram_info *info,
+                              struct bin_stats *hist)
 {
-	/* used for d and hkl of stream files where redundancy = 1 */
 	double delta;
 	int i;
 
-	/* ignore q value if outside of range */
-	if ((q<(*info).q_min) || (q>(*info).q_max)) {
+	/* Ignore q value if outside of range */
+	if ( (q<(*info).q_min) || (q>(*info).q_max) ) {
 		return 1;
 	}
 	i = find_q_bin_index(q, info, hist);
-	
+
 	delta = intensity - hist[i].mean;
 	hist[i].N++;
-	hist[i].total += intensity;	
+	hist[i].total += intensity;
 	hist[i].mean = hist[i].mean + delta /hist[i].N;
-	hist[i].std_dev = hist[i].std_dev + 
-	                      (delta *(intensity - hist[i].mean));
-	return 0;		
+	hist[i].std_dev = hist[i].std_dev + (delta *(intensity - hist[i].mean));
+
+	return 0;
 }
 
-static int add_hkl_to_histogram(double q, double intensity,
-   int redundancy, int q_scaling, struct histogram_info *info, 
-   struct bin_stats *hist)
+
+static int add_hkl_to_histogram(double q, double intensity, int redundancy,
+                                int q_scaling, struct histogram_info *info,
+                                struct bin_stats *hist)
 {
-	int i = 0; 
+	int i = 0;
 
-	/* ignore q value if outside of range */
-	if ((q<(*info).q_min) || (q>(*info).q_max)) {
+	/* Ignore q value if outside of range */
+	if ( (q<(*info).q_min) || (q>(*info).q_max) ) {
 		return 1;
 	}
-	
-	/* The accounting is the intensity of the reflection times the 
-	   The number of occurance of that reflection smeered out over the
-	   surface area which is 4 Pi q^2  the 4 pi is left out as it is a 
-	   common consant and the total is in arbiturary units.
-	*/
-	for (i=0; i<redundancy; i++) {
-		if (q_scaling) {
+
+	/* The accounting is the intensity of the reflection times the
+	 * number of occurance of that reflection smeared out over the
+	 * surface area which is 4*pi*q^2  the 4*pi is left out since it is a
+	 * common constant and the total is in arbitrary units.
+	 */
+	for ( i=0; i<redundancy; i++ ) {
+		if ( q_scaling ) {
 			add_d_to_histogram(q, intensity/(q*q), info, hist);
 		} else {
 			add_d_to_histogram(q, intensity, info, hist);
 		}
 	}
-	return 0;	
+
+	return 0;
 }
 
 
-static int histogram_setup(struct histogram_info *info, 
+static int histogram_setup(struct histogram_info *info,
                            struct bin_stats *histdata)
 {
 	int i;
 	double x;
-	/* setup histogram*/
-	if ((*info).spacing == LINEAR) {x = 1.0;}
-	else if ((*info).spacing == q2) {x = 2.0;}
-	else {x = 3.0;}
-	for ( i=0; i<(*info).histsize; i++) {
-		histdata[i].N             =   0;
-		histdata[i].total         = 0.0;
-		histdata[i].mean          = 0.0;
-		histdata[i].std_dev       = 0.0;
+
+	if ( info->spacing == LINEAR ) {
+		x = 1.0;
+	} else if ( info->spacing == q2 ) {
+		x = 2.0;
+	} else {
+		x = 3.0;
+	}
+
+	for ( i=0; i<info->histsize; i++ ) {
+
+		double qd, qm;
+
+		histdata[i].N  = 0;
+		histdata[i].total = 0.0;
+		histdata[i].mean = 0.0;
+		histdata[i].std_dev = 0.0;
 		histdata[i].fract = 0.0;
-		histdata[i].q_min         = pow((  i    *(*info).q_delta) +
-		                             pow((*info).q_min, x), 1.0/x);
-		histdata[i].q_max         = pow(((i+1.0)*(*info).q_delta) + 
-		                             pow((*info).q_min, x), 1.0/x);
-		histdata[i].q_value       = pow(((i+0.5)*(*info).q_delta) +
-		                             pow((*info).q_min, x), 1.0/x);		
+
+		qd = info->q_delta;
+		qm = info->q_min;
+
+		histdata[i].q_min = pow( (i*qd) + pow(qm, x), 1.0/x);
+
+		histdata[i].q_max = pow( ((i+1.0)*qd) + pow(qm, x), 1.0/x);
+
+		histdata[i].q_value= pow( ((i+0.5)*qd) + pow(qm, x), 1.0/x);
+
 	}
 	return 0;
 }
 
+
 static int ring_fraction_calc(struct histogram_info *info,
-              struct bin_stats *hist, struct image *image)
+                              struct bin_stats *hist, struct image *image)
 {
-	struct panel *p;
-	struct rvec r;
-	int i,fs,ss;
-	double q, q_fs, q_ss;
-	double pi = 3.14159265;
+	int fs,ss;
+	int bin;
 
-	/* check if detector geometry is loaded */
-	if ( ((*image).det == NULL) || ((*image).lambda < 0.0)) {
-		return 1;
-	}
+	/* Check that detector geometry is present and wavelength is valid */
+	if ( (image->det == NULL) || (image->lambda < 0.0) ) return 1;
 
-	/* loop over all pixels */	
+	/* Loop over all pixels */
 	for ( ss=0; ss<(*image).height; ss++ ) {
 	for ( fs=0; fs<(*image).width;  fs++ ) {
-		r = get_q(image, fs, ss, NULL, 1.0/ (*image).lambda);
+
+		struct panel *p;
+		struct rvec r;
+		int i;
+		double q, q_fs, q_ss;
+
+		r = get_q(image, fs, ss, NULL, 1.0/image->lambda);
 		q = modulus(r.u, r.v, r.w);
-		
-		/* if pixel is valid (not a bad pixel and not out of range) */		
-		if ( (q>(*info).q_min) && (q<(*info).q_max) && 
-		     (in_bad_region((*image).det,fs,ss) == 0) ) {
-		
-		/* using the panel to determine the pixel size 
-		   takes care of edge pixels correctly */
-			p = find_panel((*image).det, fs, ss);
+
+		/* If pixel is valid (not a bad pixel and not out of range) */
+		if ( (q>info->q_min) && (q<info->q_max) &&
+		     (in_bad_region(image->det,fs,ss) == 0) ) {
+
+			/* Select the panel, then (sometimes) ask for the q
+			 * of the (corner of the) pixel one step beyond the
+			 * edge, to get the exact size of the required pixel.
+			 */
+			p = find_panel(image->det, fs, ss);
+
 			r = get_q_for_panel(p, (fs+1)-(double)p->min_fs,
-			       ss-(double)p->min_ss, 
-			       NULL, 1.0/(*image).lambda);
+			                       ss-(double)p->min_ss,
+			                       NULL, 1.0/image->lambda);
 			q_fs =  modulus(r.u, r.v, r.w);
+
 			r = get_q_for_panel(p, fs-(double)p->min_fs,
-			       (ss+1) -(double)p->min_ss, 
-			        NULL, 1.0/(*image).lambda);
+			                       (ss+1) -(double)p->min_ss,
+			                       NULL, 1.0/image->lambda);
 			q_ss =  modulus(r.u, r.v, r.w);
+
 			i = find_q_bin_index(q, info, hist);
-			hist[i].fract = hist[i].fract + 
-			   fabs((q_fs - q) * (q_ss - q)); //approx pixel size
-		}  
+
+			hist[i].fract = hist[i].fract + fabs((q_fs-q)*(q_ss-q));
+
+		}
+
 	}
 	}
-	/* divide measured area by ring area */
-	for ( i=0; i<(*info).histsize; i++) {
-		hist[i].fract = hist[i].fract /
-		   (pi * (pow(hist[i].q_max,2.0) - pow(hist[i].q_min,2.0)));
+
+	/* Divide measured area by ring area */
+	for ( bin=0; bin<info->histsize; bin++ ) {
+
+		double inner_area, outer_area, ring_area;
+
+		outer_area = pow(hist[bin].q_max, 2.0);
+		inner_area = pow(hist[bin].q_min, 2.0);
+		ring_area = M_PI*(outer_area - inner_area);
+
+		hist[bin].fract = hist[bin].fract / ring_area;
+
 	}
+
 	return 0;
 }
 
+
 static unsigned int process_h5(struct image *image, struct histogram_info *info,
-                     struct bin_stats *histdata, unsigned int processing_total)
+                               struct bin_stats *histdata)
 {
 	int fs, ss;
 	double intensity;
-	for ( ss=0; ss<(*image).height; ss++ ) {
-	for ( fs=0; fs<(*image).width; fs++ ) {
-		intensity = (*image).data[fs + (*image).width*ss];
-		progress_bar(fs + (*image).width*ss, processing_total,
-                             "Processing");	
-		if (in_bad_region((*image).det,fs,ss) == 0) {
-			add_peak_to_histogram(fs, ss, intensity, 
+
+	for ( ss=0; ss<image->height; ss++ ) {
+	for ( fs=0; fs<image->width; fs++ ) {
+
+		intensity = image->data[fs + image->width*ss];
+		if ( !in_bad_region(image->det,fs,ss) ) {
+			add_peak_to_histogram(fs, ss, intensity,
 		                              image, info, histdata);
-			}
 		}
-		}
-	return 0;	
-}	
+	}
+	progress_bar(ss, image->height, "Processing");
+	}
+	return 0;
+}
+
 
 static unsigned int process_hkl(struct image *image, char *sym, UnitCell *cell,
-                       struct histogram_info *info, struct bin_stats *histdata,
-                       int q_scaling,unsigned int processing_total, int use_redundency)
+                                struct histogram_info *info,
+                                struct bin_stats *histdata,
+                                int q_scaling, int use_redundancy)
 {
-		Reflection *refl;
-		RefListIterator *iter;
-		unsigned int i = 0;
-		unsigned int n_peaks = 0;
-		int h,k,l,redundancy; 
-		double q, intensity;
-		for ( refl = first_refl((*image).reflections, &iter);
-		      refl != NULL;
-		      refl = next_refl(refl, iter) ) {
-			i++;	
-			progress_bar(i, processing_total, "Processing");
-			get_indices(refl, &h, &k, &l);
-			intensity = get_intensity(refl);
-			if (use_redundency == 1) {
-				redundancy = get_redundancy(refl);
-			} else {			
-				redundancy = num_equivs(h, k, l, sym);
-			}
-			 /*note: resolution returns 1/d/2 in m */
-			q = 2.0 * resolution(cell, h, k, l);
-			//sigma = get_esd_intensity(refl);
-			add_hkl_to_histogram(q, intensity, redundancy,
-			   q_scaling, info, histdata);
-			n_peaks+=redundancy;
+	Reflection *refl;
+	RefListIterator *iter;
+	unsigned int i = 0;
+	unsigned int n_peaks = 0;
+	int h, k, l, redundancy;
+	double q, intensity;
+	unsigned int nref;
+
+	nref = num_reflections(image->reflections);
+
+	for ( refl = first_refl(image->reflections, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		get_indices(refl, &h, &k, &l);
+		intensity = get_intensity(refl);
+		if ( use_redundancy ) {
+			redundancy = get_redundancy(refl);
+		} else {
+			redundancy = num_equivs(h, k, l, sym);
 		}
+
+		/* Multiply by 2 to get 1/d (in m^-1) */
+		q = 2.0 * resolution(cell, h, k, l);
+
+		add_hkl_to_histogram(q, intensity, redundancy, q_scaling,
+		                     info, histdata);
+
+		n_peaks += redundancy;
+
+		i++;
+		progress_bar(i, nref, "Processing");
+
+	}
+
 	return n_peaks;
 }
 
-static unsigned int process_stream_reflection(FILE *fh,  struct image *image,
-      struct histogram_info *info, struct bin_stats *histdata, 
-      unsigned int processing_total, unsigned int *n_patterns)
+
+static unsigned int process_stream_reflection(FILE *fh, struct image *image,
+                                              struct histogram_info *info,
+                                              struct bin_stats *histdata,
+                                              unsigned int *n_patterns)
 {
 	int rval;
 	unsigned int i = 0;
@@ -306,175 +360,263 @@ static unsigned int process_stream_reflection(FILE *fh,  struct image *image,
 	Reflection *refl;
 	RefListIterator *iter;
 	double intensity, fs_double, ss_double;
+	unsigned int processing_total;
+
+	processing_total = count_patterns(fh);
+	rewind(fh);
+
 	do {
-	/* Get data from next chunk */
-		rval = read_chunk(fh, image);				
+		/* Get data from next chunk */
+		rval = read_chunk(fh, image);
 		if ( rval ) continue;
-		i++;
-		progress_bar(i, processing_total, "Processing");			
-		// check if the pattern indexed, if so use those peaks	
-		if ((*image).reflections != NULL) { 
-			(*n_patterns)++; //inc number of patterns used
-			for ( refl = first_refl((*image).reflections, &iter);
-		      	      refl != NULL;	
+
+		/* Check if the pattern indexed, if so use those peaks */
+		if ( image->reflections != NULL ) {
+
+			(*n_patterns)++;
+			for ( refl = first_refl(image->reflections, &iter);
+			      refl != NULL;
 		      	      refl = next_refl(refl, iter) ) {
-				/* note added fs_double as fs is an int */										
+				/* note added fs_double as fs is an int */
 				intensity = get_intensity(refl);
 				get_detector_pos(refl, &fs_double, &ss_double);
-				if(add_peak_to_histogram(fs_double, ss_double, 
-				    intensity, image, info, histdata) == 0) {
+
+				if ( !add_peak_to_histogram(fs_double,
+				                            ss_double,
+				                            intensity,
+				                            image, info,
+				                            histdata) )
+				{
 					n_peaks++;
 				}
+
 			}
+
 		}
-		free((*image).filename);
-		reflist_free((*image).reflections);
-		image_feature_list_free((*image).features);
-		cell_free((*image).indexed_cell);	
-	} while (rval == 0);	
+
+		free(image->filename);
+		reflist_free(image->reflections);
+		image_feature_list_free(image->features);
+		cell_free(image->indexed_cell);
+
+		i++;
+		progress_bar(i, processing_total, "Processing");
+
+	} while ( rval == 0 );
+
 	return n_peaks;
 }
 
-static unsigned int process_stream_d(FILE *fh,  struct image *image,
-      struct histogram_info *info, struct bin_stats *histdata, 
-      unsigned int processing_total, unsigned int *n_patterns)
+
+static unsigned int process_stream_d(FILE *fh, struct image *image,
+                                     struct histogram_info *info,
+                                     struct bin_stats *histdata,
+                                     unsigned int *n_patterns)
 {
-	int h, k, l, rval; 
+	int h, k, l, rval;
 	unsigned int i = 0;
-	unsigned int n_peaks = 0;	
+	unsigned int n_peaks = 0;
 	Reflection *refl;
 	RefListIterator *iter;
 	double intensity, q;
+	unsigned int processing_total;
+
+	processing_total = count_patterns(fh);
+	rewind(fh);
+
 	do {
-	/* Get data from next chunk */
-		rval = read_chunk(fh, image);				
+
+		/* Get data from next chunk */
+		rval = read_chunk(fh, image);
 		if ( rval ) continue;
-		i++;
-		progress_bar(i, processing_total, "Processing");			
-		// check if the pattern indexed, if so use those peaks				
-		if ((*image).reflections != NULL) { 
+
+		if ( image->reflections != NULL ) {
+
 			(*n_patterns)++; //inc number of patterns used
-			for ( refl = first_refl((*image).reflections, &iter);
-	      		      refl != NULL;	
-      			      refl = next_refl(refl, iter) ) {
+
+			for ( refl = first_refl(image->reflections, &iter);
+			      refl != NULL;
+			      refl = next_refl(refl, iter) )
+			{
 				get_indices(refl, &h, &k, &l);
 				intensity = get_intensity(refl);
-				/*resolution returns 1/d/2 in nm */
-				q = 2.0 * resolution((*image).indexed_cell, h, k, l); 
-				if (add_d_to_histogram(q, intensity, info, 
-				                       histdata) == 0) {
-					n_peaks++;
-				}
+				q = 2.0 * resolution(image->indexed_cell,
+				                     h, k, l);
+				if ( !add_d_to_histogram(q, intensity, info,
+				                        histdata)) n_peaks++;
 			}
 		}
-		free((*image).filename);
-		reflist_free((*image).reflections);
-		image_feature_list_free((*image).features);
-		cell_free((*image).indexed_cell);	
-	} while (rval == 0);	
+
+		free(image->filename);
+		reflist_free(image->reflections);
+		image_feature_list_free(image->features);
+		cell_free(image->indexed_cell);
+
+		i++;
+		progress_bar(i, processing_total, "Processing");
+
+	} while ( rval == 0 );
+
 	return n_peaks;
 }
 
 
-static unsigned int process_stream_hkl(FILE *fh,  struct image *image,
-      struct histogram_info *info, struct bin_stats *histdata,  UnitCell *cell,
-      unsigned int processing_total, unsigned int *n_patterns)
+static unsigned int process_stream_hkl(FILE *fh, struct image *image,
+                                       struct histogram_info *info,
+                                       struct bin_stats *histdata,
+                                       UnitCell *cell, unsigned int *n_patterns)
 {
-	int h, k, l, rval; 
+	int rval;
 	unsigned int i = 0;
-	unsigned int n_peaks = 0;	
+	unsigned int n_peaks = 0;
 	Reflection *refl;
 	RefListIterator *iter;
 	double intensity, q;
+	unsigned int processing_total;
+
+	processing_total = count_patterns(fh);
+	rewind(fh);
+
 	do {
-	/* Get data from next chunk */
-		rval = read_chunk(fh, image);				
+
+		/* Get data from next chunk */
+		rval = read_chunk(fh, image);
 		if ( rval ) continue;
-		i++;
-		progress_bar(i, processing_total, "Processing");
-		if ((*image).reflections != NULL) {
-			n_patterns++; //inc number of patterns used
-			for ( refl = first_refl((*image).reflections, &iter);
+		if ( image->reflections != NULL ) {
+
+			(*n_patterns)++;
+
+			for ( refl = first_refl(image->reflections, &iter);
 		      	      refl != NULL;
-		      	      refl = next_refl(refl, iter) ) {
+			      refl = next_refl(refl, iter) )
+			{
+				int h, k, l;
+
 				get_indices(refl, &h, &k, &l);
 				intensity = get_intensity(refl);
 				q = 2.0 * resolution(cell, h, k, l);
-				if(add_d_to_histogram(q, intensity, info,
-				                      histdata) == 0) {
-					n_peaks++;
-				}
+
+				if ( !add_d_to_histogram(q, intensity, info,
+				                         histdata) ) n_peaks++;
 			}
-		}	
-		free((*image).filename);
-		reflist_free((*image).reflections);
-		image_feature_list_free((*image).features);
-		cell_free((*image).indexed_cell);	
-	} while (rval == 0);	
+		}
+
+		free(image->filename);
+		reflist_free(image->reflections);
+		image_feature_list_free(image->features);
+		cell_free(image->indexed_cell);
+
+		i++;
+		progress_bar(i, processing_total, "Processing");
+
+	} while ( rval == 0 );
+
 	return n_peaks;
 }
 
-static unsigned int process_stream_peaks(FILE *fh,  struct image *image,
-      struct histogram_info *info, struct bin_stats *histdata, 
-      unsigned int processing_total, unsigned int *n_patterns, int only_indexed)
+
+static int add_features_to_histogram(struct image *image,
+                                     struct histogram_info *info,
+                                     struct bin_stats *histdata)
+{
+	int j;
+	int n_peaks;
+
+	n_peaks = 0;
+	for ( j=0; j<image_feature_count(image->features); j++) {
+
+		struct imagefeature *f;
+
+		f = image_get_feature(image->features, j);
+		if ( !f->valid ) continue;
+
+		if ( !in_bad_region(image->det, f->fs,f->ss) ) {
+
+			int r;
+
+			r = add_peak_to_histogram(f->fs, f->ss, f->intensity,
+			                          image, info, histdata);
+
+			if ( !r ) n_peaks++;
+
+		}
+
+	}
+
+	return n_peaks;
+}
+
+
+static unsigned int process_stream_peaks(FILE *fh, struct image *image,
+                                         struct histogram_info *info,
+                                         struct bin_stats *histdata,
+                                         unsigned int *n_patterns,
+                                         int only_indexed)
 {
-	struct imagefeature *f;	
-	int rval; 
+	int rval;
 	unsigned int i = 0;
 	unsigned int n_peaks = 0;
-	unsigned int j;		
+	unsigned int processing_total;
+
+	processing_total = count_patterns(fh);
+	rewind(fh);
+
 	do {
-	/* Get data from next chunk */
-		rval = read_chunk(fh, image);				
+
+		/* Get data from next chunk */
+		rval = read_chunk(fh, image);
 		if ( rval ) continue;
-		i++;
-		progress_bar(i, processing_total, "Processing");
-		if ((*image).features != NULL) {
-			if ((only_indexed == 0) || 
-			   ((only_indexed == 1) && 
-			    ((*image).reflections != NULL))) {											
-				(*n_patterns)++;	
-				for (j = 0;
-				     j<image_feature_count((*image).features);
-				     j++) {			
-					f = image_get_feature((*image).features, j);
-					if (in_bad_region((*image).det,
-					   (*f).fs,(*f).ss) == 0) {
-						if (add_peak_to_histogram((*f).fs,
-						   (*f).ss, (*f).intensity, image,
-						   info, histdata) == 0) {
-							n_peaks++;
-						}
-					}
-				}
+
+		if ( image->features != NULL ) {
+
+			if ( (!only_indexed)
+			  || ( only_indexed && (image->reflections != NULL)) )
+			{
+				(*n_patterns)++;
+				n_peaks += add_features_to_histogram(image,
+				                                     info,
+				                                     histdata);
 			}
 		}
-		free((*image).filename);
-		reflist_free((*image).reflections);
-		image_feature_list_free((*image).features);
-		cell_free((*image).indexed_cell);	
-	} while (rval == 0);	
+
+		free(image->filename);
+		reflist_free(image->reflections);
+		image_feature_list_free(image->features);
+		cell_free(image->indexed_cell);
+
+		i++;
+		progress_bar(i, processing_total, "Processing");
+
+	} while ( rval == 0 );
+
 	return n_peaks;
 }
 
-static unsigned int process_stream_h5(FILE *fh,  struct image *image,
-      struct histogram_info *info, struct bin_stats *histdata, int config_satcorr, 
-       int only_indexed, unsigned int processing_total, unsigned int *n_patterns)
+
+static unsigned int process_stream_h5(FILE *fh, struct image *image,
+                                      struct histogram_info *info,
+                                      struct bin_stats *histdata,
+                                      int config_satcorr, int only_indexed,
+                                      unsigned int *n_patterns)
 {
 	int fs, ss, rval;
 	double intensity;
 	unsigned int i = 0;
 	unsigned int n_peaks = 0;
-	struct hdfile *hdfile = NULL;	
+	struct hdfile *hdfile = NULL;
+	unsigned int processing_total;
+
+	processing_total = count_patterns(fh);
+	rewind(fh);
+
 	do {
-	/* Get data from next chunk */
-		rval = read_chunk(fh, image);				
+
+		/* Get data from next chunk */
+		rval = read_chunk(fh, image);
 		if ( rval ) continue;
-		i++;
-		progress_bar(i, processing_total, "Processing");
-		if ((only_indexed == 0) || 
-		   ((only_indexed == 1) && 
-		    ((*image).reflections != NULL))) {	
+		if ((only_indexed == 0) ||
+		   ((only_indexed == 1) &&
+		    ((*image).reflections != NULL))) {
 			hdfile = hdfile_open((*image).filename);
 			hdfile_set_image(hdfile, "/data/data");
 			hdf5_read(hdfile, image, config_satcorr);
@@ -482,23 +624,35 @@ static unsigned int process_stream_h5(FILE *fh,  struct image *image,
 			n_patterns++;
 			for ( ss=0; ss<(*image).height; ss++ ) {
 			for ( fs=0; fs<(*image).width;  fs++ ) {
-				intensity = (*image).data[fs + (*image).width*ss];
-				if (in_bad_region((*image).det,fs,ss) == 0) {
-					add_peak_to_histogram(fs, ss, intensity, 
-					   image, info, histdata);
+
+				intensity = image->data[fs + image->width*ss];
+				if ( !in_bad_region(image->det,fs,ss) ) {
+
+					add_peak_to_histogram(fs, ss, intensity,
+					                      image, info,
+					                      histdata);
+
 				}
+
 			}
 			}
 		}
-		free((*image).data);
-		free((*image).filename);
-		reflist_free((*image).reflections);
-		image_feature_list_free((*image).features);
-		cell_free((*image).indexed_cell);		
-	} while (rval == 0);	
+
+		free(image->data);
+		free(image->filename);
+		reflist_free(image->reflections);
+		image_feature_list_free(image->features);
+		cell_free(image->indexed_cell);
+
+		i++;
+		progress_bar(i, processing_total, "Processing");
+
+	} while ( rval == 0 );
+
 	return n_peaks;
 }
 
+
 static void show_help(const char *s)
 {
 	printf("Syntax: %s [options]\n\n", s);
@@ -507,19 +661,19 @@ static void show_help(const char *s)
 "\n"
 "  -h, --help              Display this help message.\n"
 "  -i, --input=<file>      Input filename. (*.stream, *.hkl, or *.h5)\n"
-"  -o, --output=<file>     Output filename.  (default stdout)\n"
-"  -g. --geometry=<file>   Get detector geometry from file.\n"
-"  -b, --beam=<file>       Get beam parameters (wavelength) from file.\n"
+"  -o, --output=<file>     Output filename. Default: stdout.\n"
+"  -g. --geometry=<file>   Get detector geometry from <file>.\n"
+"  -b, --beam=<file>       Get beam parameters (wavelength) from <file>.\n"
 "  -p, --pdb=<file>        Get unit cell from PDB file. (.hkl files only)\n"
 "  -y, --symmetry=<sym>    The symmetry of crystal (.hkl files only)\n"
-"  -s, --bins=n            Makes histogram wiht n bins (default is 100).\n"
+"  -s, --bins=n            Makes histogram with n bins (default is 100).\n"
 "      --spacing=<type>    Use 'type' to select the q spacing.\n"
 "                          Choose from:\n"
 "                            linear      : linear (default)\n"
 "                            q2          : even spacing in Wilson plots\n"
-"                            volume      : constant volume spacing\n"
-"      --q-max=n           The max q to be considered in plot.\n"
-"      --q-min=n           The min q to be considered in plot.\n"
+"                            volume      : constant volume\n"
+"      --q-max=n           The maximum q to be considered in plot.\n"
+"      --q-min=n           The minimum q to be considered in plot.\n"
 "  -d, --data=<type>       Use to select the kind of stream data in histogram.\n"
 "                          Choose from:\n"
 "                            reflection  : uses peak positons from indexed\n"
@@ -534,51 +688,53 @@ static void show_help(const char *s)
 "                                          regions\n"
 "     --no-sat-corr        Don't correct values of saturated peaks using a\n"
 "                            table included in the HDF5 file.\n"
-"     --only-indexed       Use wiht -data=peaks or h5 if you want to use the\n"
+"     --only-indexed       Use with -data=peaks or h5 if you want to use the\n"
 "                            peak list of only indexed patterns\n"
 "     --no-q-scaling       Use with .hkl files if you want to not scale the\n"
 "                            powder by 1/q^2\n"
 "     --ring-corr          Use if you want to scale the powder plot to\n"
 "                            correct for the fractional area sampled of the\n"
 "                            powder ring\n"
-"    --use-redundency      Use with .hkl files if you want to use the number\n"
+"     --use-redundancy     Use with .hkl files if you want to use the number\n"
 "                            of measurements and not the number of symetrical\n"
 "                            equivelent reflections as the number of time a\n"
 "                            reflection occurs in the powder\n"
 "\n");
 }
 
+
 int main(int argc, char *argv[])
 {
 	FILE *fh = NULL;
-	UnitCell *cell = NULL;	
-	struct image image;	
+	UnitCell *cell = NULL;
+	struct image image;
 	struct hdfile *hdfile = NULL;
 	struct bin_stats *histdata = NULL;
 	struct histogram_info hist_info;
 
-	hist_info.histsize =  100;	//default settings
+	/* Default settings */
+	hist_info.histsize =  100;
 	hist_info.q_min    = -1.0;
 	hist_info.q_max    = -1.0;
 	hist_info.spacing  = LINEAR;
 	image.lambda = -1.0;
 	image.beam = NULL;
 	image.det  = NULL;
-	
+
 	unsigned int n_patterns = 0;
-	unsigned int n_peaks    = 0;	
-	
+	unsigned int n_peaks    = 0;
+
 	int c, rval, file_type, data_type;
-	int config_satcorr = 1; //true by default
-	int need_geometry  = 0; //false
-	int need_beam      = 0; //false
-	int need_pdb       = 0; //false
-	int only_indexed   = 0; //false
-	int q_scaling      = 1; //true by default
-	int ring_corr      = 0; //false
-	int use_redundency = 0; //false
-	unsigned int i, processing_total;	
-	
+	int config_satcorr = 1;
+	int need_geometry  = 0;
+	int need_beam      = 0;
+	int need_pdb       = 0;
+	int only_indexed   = 0;
+	int q_scaling      = 1;
+	int ring_corr      = 0;
+	int use_redundancy = 0;
+	unsigned int i;
+
 	char *filename = NULL;
 	char *geometry = NULL;
 	char *beamf    = NULL;
@@ -605,16 +761,18 @@ int main(int argc, char *argv[])
 		{"only-indexed",       0, &only_indexed,       1 },
 		{"no-q-scaling",       0, &q_scaling,          0 },
 		{"ring-corr",          0, &ring_corr,          1 },
-		{"use-redundency",     0, &use_redundency,     1 },
+		{"use-redundancy",     0, &use_redundancy,     1 },
 		{"data",               1, NULL,               'd'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:o:g:b:p:s:d:y:", 
-	        longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hi:o:g:b:p:s:d:y:",
+	        longopts, NULL)) != -1)
+	{
 
 		switch (c) {
+
 		case 'h' :
 			show_help(argv[0]);
 			return 0;
@@ -647,15 +805,15 @@ int main(int argc, char *argv[])
 			hist_info.histsize = atoi(optarg);
 			break;
 
-		case 1 :			
+		case 1 :
 			hist_info.q_max = atof(optarg);
 			break;
 
-		case 2 :			
+		case 2 :
 			hist_info.q_min = atof(optarg);
 			break;
 
-		case 3 :			
+		case 3 :
 			if (strcmp(optarg, "linear") == 0 ) {
 				hist_info.spacing = LINEAR;
 			} else if (strcmp(optarg, "q2") == 0 ) {
@@ -663,8 +821,7 @@ int main(int argc, char *argv[])
 			} else if (strcmp(optarg, "volume") == 0) {
 				hist_info.spacing = VOLUME;
 			} else {
-				ERROR("Failed to read spacing plot type: '%s'\n",
-				      optarg);
+				ERROR("Invalid spacing type: '%s'\n", optarg);
 				return 1;
 			}
 			break;
@@ -679,79 +836,93 @@ int main(int argc, char *argv[])
 		default :
 			return 1;
 		}
+
 	}
 
 	/* Process input file type */
 	if ( filename == NULL ) {
+
 		ERROR("You must specify the input filename with -i\n");
 		return 1;
-	} else if ( is_stream(filename) == 1 ) { //open .stream file
+
+	}
+
+	if ( is_stream(filename) ) {
+
 		file_type = FILE_STREAM;
 		fh = fopen(filename, "r");
 		if ( fh == NULL ) {
 			ERROR("Failed to open input file\n");
 			return 1;
 		}
-		rval = read_chunk(fh, &image); //read chunk to get wavelength
+		/* Use wavelength from first chunk */
+		rval = read_chunk(fh, &image);
 		rewind(fh);
-		processing_total = count_patterns(fh);
-		rewind(fh);
-		
-	} else if ( H5Fis_hdf5(filename) > 0) { //open .h5 file
+
+	} else if ( H5Fis_hdf5(filename) > 0 ) {
+
 		file_type = FILE_H5;
-		need_geometry = 1; //true
+		need_geometry = 1;
 		hdfile = hdfile_open(filename);
 		hdfile_set_image(hdfile, "/data/data");
 		hdf5_read(hdfile, &image, config_satcorr);
 		hdfile_close(hdfile);
-		processing_total = image.width * image.height;
+
 	} else {
+
 		image.reflections = read_reflections(filename);
-		if (image.reflections != NULL) { //open .hkl file
+
+		if ( image.reflections != NULL ) {
 			file_type = FILE_HKL;
 			need_pdb = 1;
-			processing_total = num_reflections(image.reflections);
 		} else {
-			ERROR("Input file: %s is an invalid type.\n",filename);
+			ERROR("Couldn't recognise %s as reflection list,"
+			      " stream or image.\n", filename);
 			return 1;
 		}
+
 	}
 	free(filename);
 
-	/*Process data read type */
 	if ( datatype == NULL ) {
 		data_type = PLOT_D;
 		if ((hist_info.q_min < 0.0) || (hist_info.q_max < 0.0)) {
 			need_geometry = 1;
 		}
+
 	} else if ( strcmp(datatype, "reflection") == 0 ) {
 		data_type = PLOT_REFL;
 		need_geometry = 1;
+
 	} else if ( strcmp(datatype, "hkl") == 0 ) {
 		data_type = PLOT_HKL;
 		need_pdb = 1;
 		if ((hist_info.q_min <= 0.0) || (hist_info.q_max <= 0.0)) {
 			need_geometry = 1;
 		}
+
 	} else if ( strcmp(datatype, "d") == 0 ) {
 		data_type = PLOT_D;
 		if ((hist_info.q_min <= 0.0) || (hist_info.q_max <= 0.0)) {
 			need_geometry = 1;
 		}
+
 	} else if ( strcmp(datatype, "peaks") == 0 ) {
 		data_type = PLOT_PEAKS;
 		need_geometry = 1;
-	}
-	else if ( strcmp(datatype, "h5") == 0 ) {
+
+	} else if ( strcmp(datatype, "h5") == 0 ) {
 		data_type = PLOT_H5;
 		need_geometry = 1;
+
 	} else {
+
 		ERROR("Failed to read data plot type: '%s'\n", datatype);
 		return 1;
 	}
 
-	/* logic checks */	
-	if ((need_geometry) && (image.lambda < 0.0)) {
+	/* Logic checks */
+	if ( need_geometry && (image.lambda < 0.0) ) {
 		need_beam = 1;
 	}
 	if ( hist_info.histsize <= 0 ) {
@@ -759,16 +930,16 @@ int main(int argc, char *argv[])
                       "bins\n");
 		return 1;
 	}
-	if (hist_info.q_min > hist_info.q_max) {
+	if ( hist_info.q_min > hist_info.q_max ) {
 		ERROR("the minimum q value of: %e "
-	              "is greator then your max q value of: %e\n", 
+	              "is greator then your max q value of: %e\n",
                       hist_info.q_min, hist_info.q_max);
 		return 1;
 	}
 
-	/*get geometry, beam and pdb files and parameters as needed */
-	if (need_geometry) {
-		if (geometry == NULL) {
+	/* Get geometry, beam and pdb files and parameters as needed */
+	if ( need_geometry ) {
+		if ( geometry == NULL ) {
 			ERROR("You need to specify a geometry file with "
 		              "--geometry\n");
 			return 1;
@@ -782,10 +953,11 @@ int main(int argc, char *argv[])
 			image.width  = image.det->max_fs;
 			image.height = image.det->max_ss;
 		}
-	}	
+	}
 	free(geometry);
-	
-	if (need_beam) {
+
+	if ( need_beam ) {
+
 		if ( beamf == NULL ) {
 			ERROR("No wavelength in file, so you need to specify "
 			      "a beam parameters file with --beam\n");
@@ -800,10 +972,12 @@ int main(int argc, char *argv[])
 			image.lambda = ph_en_to_lambda(eV_to_J(
 			                 image.beam->photon_energy));
 		}
+
 	}
 	free(beamf);
 
-	if (need_pdb) {
+	if ( need_pdb ) {
+
 		if (pdb == NULL) {
 			ERROR("You need to specify a pdb file with --pdb.\n");
 			return 1;
@@ -815,6 +989,7 @@ int main(int argc, char *argv[])
 				return 1;
 			}
 		}
+
 	}
 	free(pdb);
 
@@ -822,82 +997,80 @@ int main(int argc, char *argv[])
 		sym = strdup("1");
 	}
 
-	/* setup histogram info*/
-	if (hist_info.q_min <= 0 ) {
+	/* Set up histogram info*/
+	if (hist_info.q_min <= 0.0 ) {
 		hist_info.q_min = smallest_q(&image);
 	}
-	if (hist_info.q_max <= 0 ) {
+	if (hist_info.q_max <= 0.0 ) {
 		hist_info.q_max = largest_q(&image);
 	}
 	if (hist_info.spacing == LINEAR) {
 		hist_info.q_delta = (hist_info.q_max - hist_info.q_min)/
 		                    hist_info.histsize;
 	} else if (hist_info.spacing == q2) {
-		hist_info.q_delta = (pow(hist_info.q_max, 2.0) - 
-		                     pow(hist_info.q_min, 2.0)) / 
+		hist_info.q_delta = (pow(hist_info.q_max, 2.0) -
+		                     pow(hist_info.q_min, 2.0)) /
 		                     hist_info.histsize;
 	} else {  //by default must be in VOLUME
 		hist_info.q_delta = (pow(hist_info.q_max, 3.0) -
-		                     pow(hist_info.q_min, 3.0)) / 
+		                     pow(hist_info.q_min, 3.0)) /
 		                     hist_info.histsize;
 	}
-	/* setup histogram data */
+	/* Set up histogram data */
 	histdata = malloc((hist_info.histsize) * sizeof(struct bin_stats));
 	histogram_setup(&hist_info, histdata);
 
-	/* setup ring scaling */
-	if (ring_corr == 1) {
-		if (ring_fraction_calc(&hist_info, histdata, &image) == 1) {
-			ERROR("Detector is not loaded could not correct for"
-			" finiate ring measurement. Do not use --ring-corr\n");
+	/* Set up ring scaling */
+	if ( ring_corr ) {
+
+		if ( ring_fraction_calc(&hist_info, histdata, &image) ) {
+
+			ERROR("Detector geometry is required to correct for"
+			      " finite ring area.\n");
 			return 1;
 		}
 	}
 
-	/* process reflections based on file type and data type*/
+	/* Process reflections based on file type and data type */
 	switch (file_type) {
 	case FILE_H5 :
 		n_patterns++;
-		n_peaks = process_h5(&image, &hist_info, histdata, processing_total);
+		n_peaks = process_h5(&image, &hist_info, histdata);
 		free(image.data);
 		break;
 
 	case FILE_HKL :
 		n_patterns++; //inc number of patterns used
 		n_peaks = process_hkl(&image, sym, cell, &hist_info, histdata,
-                                      q_scaling, processing_total, use_redundency);
+                                      q_scaling, use_redundancy);
 		break;
 
 	case FILE_STREAM :
-			switch (data_type) {
-			case PLOT_REFL :
-				n_peaks = process_stream_reflection(fh, &image,
-				             &hist_info, histdata, processing_total,
-				             &n_patterns);
-				break;
-			case PLOT_D :
-				n_peaks = process_stream_d(fh, &image, &hist_info, 
-					     histdata, processing_total,
-				             &n_patterns);
-				break;
-			case PLOT_HKL :
-				n_peaks = process_stream_hkl(fh, &image, &hist_info, 
-					     histdata, cell, processing_total,
-				             &n_patterns);				
-				break;
-			case PLOT_PEAKS :
-				n_peaks = process_stream_peaks(fh, &image, &hist_info,
-				             histdata, processing_total, 
-					     &n_patterns, only_indexed);
-				break;
-			case PLOT_H5 : 			
-				n_peaks =  process_stream_h5(fh, &image, &hist_info,
-				             histdata, config_satcorr, only_indexed,
-				             processing_total, &n_patterns);
-				break;
-			default :
-				break;
-			}
+		switch (data_type) {
+		case PLOT_REFL :
+			n_peaks = process_stream_reflection(fh, &image,
+			             &hist_info, histdata, &n_patterns);
+			break;
+		case PLOT_D :
+			n_peaks = process_stream_d(fh, &image, &hist_info,
+				     histdata, &n_patterns);
+			break;
+		case PLOT_HKL :
+			n_peaks = process_stream_hkl(fh, &image, &hist_info,
+				     histdata, cell, &n_patterns);
+			break;
+		case PLOT_PEAKS :
+			n_peaks = process_stream_peaks(fh, &image, &hist_info,
+			             histdata, &n_patterns, only_indexed);
+			break;
+		case PLOT_H5 :
+			n_peaks =  process_stream_h5(fh, &image, &hist_info,
+			             histdata, config_satcorr, only_indexed,
+			             &n_patterns);
+			break;
+		default :
+			break;
+		}
 		fclose(fh);
 		break;
 
@@ -906,65 +1079,52 @@ int main(int argc, char *argv[])
 
 	}
 
-	/* sqrt the variance to get the std_dev */	
-
-	for(i=0; i<hist_info.histsize; i++) {		
+	/* Sqrt the variance to get the std_dev */
+	for( i=0; i<hist_info.histsize; i++ ) {
 		if (histdata[i].N > 1) {
 			histdata[i].std_dev = sqrt(histdata[i].std_dev/
 			                          (histdata[i].N-1));
 		}
 	}
-	if (ring_corr == 1) {
-		for(i=0; i<hist_info.histsize; i++) {
-			histdata[i].N       = histdata[i].N      / 
+	if ( ring_corr ) {
+		for( i=0; i<hist_info.histsize; i++ ) {
+			histdata[i].N       = histdata[i].N      /
 			                      histdata[i].fract;
-			histdata[i].total   = histdata[i].total  / 
+			histdata[i].total   = histdata[i].total  /
 			                      histdata[i].fract;
-			histdata[i].mean    = histdata[i].mean   / 
+			histdata[i].mean    = histdata[i].mean   /
 			                      histdata[i].fract;
-			histdata[i].std_dev = histdata[i].total  / 
+			histdata[i].std_dev = histdata[i].total  /
 			                      histdata[i].fract;
 		}
 	}
 
-	/* print out the results */
-	if (output != NULL) {
+	/* Print out the results */
+	if ( output != NULL ) {
 		fh = fopen(output, "w");
 		if ( fh == NULL ) {
 			ERROR("Failed to open output file\n");
 			return 1;
 		}
-		fprintf(fh,"I read %i patterns with %i peaks\n", 
-		        n_patterns, n_peaks);
-		fprintf(fh,"q\tN\ttotal\tmean\tstd dev\t std dev of mean\n");
-	}
-	else {	
-		printf("I read %i patterns with %i peaks\n",
-		       n_patterns, n_peaks);
-		printf("q\tN\ttotal\tmean\tstd dev\t std dev of mean\n");
-	}	
-	for(i=0; i<hist_info.histsize; i++) {
-		if (output != NULL) {
-			fprintf(fh,"%5e\t%i\t%5e\t%5e\t%5e\t%5e\n", 
-			   histdata[i].q_min, histdata[i].N, histdata[i].total,
-		           histdata[i].mean, histdata[i].std_dev,
-			   histdata[i].std_dev/sqrt(histdata[i].N));
-		}
-		else {
-			printf("%5e\t%i\t%5e\t%5e\t%5e\t%5e\n",
-		          histdata[i].q_min, histdata[i].N, histdata[i].total,
-			  histdata[i].mean, histdata[i].std_dev,
-		          histdata[i].std_dev/sqrt(histdata[i].N));
-		}
+	} else {
+		fh = stdout;
+	}
+
+	fprintf(fh, "I read %i patterns with %i peaks\n", n_patterns, n_peaks);
+	fprintf(fh, "q\tN\ttotal\tmean\tstd dev\t std dev of mean\n");
+
+	for( i=0; i<hist_info.histsize; i++ ) {
+		fprintf(fh, "%5e\t%i\t%5e\t%5e\t%5e\t%5e\n",
+		        histdata[i].q_min, histdata[i].N,
+		        histdata[i].total, histdata[i].mean,
+		        histdata[i].std_dev,
+		        histdata[i].std_dev/sqrt(histdata[i].N));
 	}
 
 	if ( cell != NULL ) cell_free(cell);
 	if ( image.det  != NULL ) free(image.det);
 	if ( image.beam != NULL ) free(image.beam);
-	if (output != NULL) {
-		fclose(fh);
-		free(output);
-	}
+	fclose(fh);
 	free(histdata);
 	free(sym);