Quantification stuff

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

1 files changed, 109 insertions, 0 deletions

diff --git a/src/intensities.c b/src/intensities.c
new file mode 100644
index 0000000..96c7ad1
--- /dev/null
+++ b/src/intensities.c
@@ -0,0 +1,109 @@
+/*
+ * intensities.c
+ *
+ * Extract integrated intensities by relrod estimation
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "control.h"
+#include "reflections.h"
+#include "image.h"
+#include "reproject.h"
+
+/* Extract integrated reflection intensities by estimating the spike function
+ * based on the observed intensity and the calculated excitation error from
+ * the lattice refinement.  Easy.	*/
+void intensities_extract(ControlContext *ctx) {
+
+	int i, j;
+	int n_meas, n_dupl, n_notf;
+	double max;
+	Reflection *reflection;
+
+	/* Free previous analysis if required */
+	if ( ctx->integrated != NULL ) {
+		reflectionlist_free(ctx->integrated);
+	}
+	ctx->integrated = reflectionlist_new();
+	
+	n_meas = 0;
+	n_dupl = 0;
+	n_notf = 0;
+	max = 0;
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+	
+		ImageRecord *image;
+		
+		image = &ctx->images->images[i];
+		if ( image->rflist == NULL ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
+		
+		for ( j=0; j<image->rflist->n_features; j++ ) {
+		
+			ImageFeature *feature;
+			signed int h, k, l;
+			
+			feature = &image->rflist->features[j];
+			
+			h = feature->reflection->h;
+			k = feature->reflection->k;
+			l = feature->reflection->l;
+				
+			if ( feature->partner != NULL ) {
+			
+				if ( (h!=0) || (k!=0) || (l!=0) ) {
+				
+					double intensity;
+					Reflection *new;
+					
+					/* Perform relrod calculation of doom here.
+					 * TODO: Figure out if this is even possible. */
+					intensity = feature->partner->intensity;
+					
+					new = reflection_add(ctx->integrated,
+							feature->reflection->x, feature->reflection->y, feature->reflection->z,
+							intensity, REFLECTION_GENERATED);
+					
+					if ( new != NULL ) {
+						new->h = h;
+						new->k = k;
+						new->l = l;
+						//printf("IN: Adding %3i %3i %3i, intensity=%f\n", h, k, l, intensity);
+						if ( intensity > max ) max = intensity;
+						n_meas++;
+					} else {
+						printf("IN: Duplicate measurement for %3i %3i %3i\n", h, k, l);
+						n_dupl++;
+					}
+				
+				}
+				
+			} else {
+				//printf("IN: %3i %3i %3i not found\n", h, k, l);
+				n_notf++;
+			}
+		
+		}
+		
+	}
+	
+	/* Normalise all reflections to the most intense reflection */
+	reflection = ctx->integrated->reflections;
+	while ( reflection ) {
+		reflection->intensity /= max;
+		reflection = reflection->next;
+	}
+	
+	printf("IN: %5i intensities measured\n", n_meas);
+	printf("IN: %5i duplicated measurements\n", n_dupl);
+	printf("IN: %5i predicted reflections not found\n", n_notf);
+
+}
+