Feature quantification

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

2 files changed, 76 insertions, 3 deletions

diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 9774a33..db1f3db 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -219,9 +219,16 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 	cur = imagedisplay->marks;
 	max = 0.0;
 	while ( cur ) {
-		if ( cur->weight > max ) max = cur->weight;
-		if ( cur->weight < 0.0 ) printf("ID: Warning: ImageDisplayMark with negative weight\n");
+	
+		if ( cur->weight < 0.0 ) {
+			printf("ID: Warning: ImageDisplayMark with negative weight\n");
+			cur = cur->next;
+			continue;
+		}
+		
+		if ( log(1+0.1*cur->weight) > max ) max = log(1+0.1*cur->weight);
 		cur = cur->next;
+		
 	}
 	
 	cur = imagedisplay->marks;
@@ -243,7 +250,12 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 		
 			double r;
 		
-			r = 20 * (cur->weight/max);
+			if ( cur->weight < 0.0 ) {
+				cur = cur->next;
+				continue;
+			}
+			
+			r = 20.0 * (log(1+0.1*cur->weight)/max);
 			
 			gdk_draw_arc(drawingarea->window, gc, FALSE,
 						xoffs + cur->x*scale - r,
diff --git a/src/itrans.c b/src/itrans.c
index ee932a9..1018a66 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -14,6 +14,8 @@
 #include <config.h>
 #endif
 
+#include <math.h>
+
 #include "image.h"
 #include "itrans-threshold.h"
 #include "itrans-zaefferer.h"
@@ -36,7 +38,66 @@ ImageFeatureList *itrans_process_image(ImageRecord *image, PeakSearchMode psmode
 	
 }
 
+/* Quantification radius */
+#define R 20
+
 /* Quantify each feature in the image's feature list */
 void itrans_quantify_features(ImageRecord *image) {
+
+	int i;
+	
+	for ( i=0; i<image->features->n_features; i++ ) {
+
+		double theta;
+		double background, total;
+		int n;
+		int xi, yi;
+		
+		background = 0.0;
+		n = 0;
+		for ( theta=0; theta<2*M_PI; theta+=2*M_PI/10 ) {
+		
+			double x, y;
+			long int xd, yd;
+			
+			x = image->features->features[i].x + R*sin(theta);
+			y = image->features->features[i].y + R*cos(theta);
+			xd = rint(x);
+			yd = rint(y);
+			
+			if ( (xd>=0) && (yd>=0) && (xd<image->width) && (yd<image->height) ) {
+				background += image->image[xd + yd*image->width];
+				n++;
+			}
+
+		}
+		background /= n;
+		
+		total = 0.0;
+		for ( xi=-20; xi<=20; xi++ ) {
+			for ( yi=-20; yi<=20; yi++ ) {
+			
+				double x, y;
+				long int xd, yd;
+				
+				if ( xi*xi + yi*yi > 20 ) continue;
+				
+				x = image->features->features[i].x + xi;
+				y = image->features->features[i].y + yi;
+				
+				xd = rint(x);
+				yd = rint(y);
+				
+				if ( (xd>=0) && (yd>=0) && (xd<image->width) && (yd<image->height) ) {
+					total += image->image[xd + yd*image->width] - background;
+				}
+			
+			}
+		}
+		
+		image->features->features[i].intensity = total;
+
+	}
+
 }