More refinement stuff

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

1 files changed, 57 insertions, 30 deletions

diff --git a/src/refine.c b/src/refine.c
index 45aec8e..7871dc1 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -43,11 +43,38 @@ void refine_do_sequence(ControlContext *ctx) {
 	
 }
 
+static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) {
+
+	double fom = 0.0;
+	int f;
+	
+	for ( f=0; f<nf; f++ ) {
+	
+		double xdf, ydf, zdf;
+		
+		xdf = d[f].h*s[i].dax + d[f].k*s[i].dbx + d[f].l*s[i].dcx;
+		ydf = d[f].h*s[i].day + d[f].k*s[i].dby + d[f].l*s[i].dcy;
+		zdf = d[f].h*s[i].daz + d[f].k*s[i].dbz + d[f].l*s[i].dcz;
+		xdf -= d[f].dx;
+		ydf -= d[f].dy;
+		zdf -= d[f].dz;
+		
+		fom += sqrt(xdf*xdf + ydf*ydf + zdf*zdf);
+		
+	}
+	
+	return fom;
+
+}
+
+static void refine_simplex_reflect(SimplexVertex *s, int v_worst) {
+}
+
 void refine_do_cell(ControlContext *ctx) {
 
 	SimplexVertex s[10];
 	Deviation *d;
-	const double delta = 0.1e9;
+	double delta;
 	int i, nf, f, v_worst;
 	double fom_worst;
 	
@@ -61,20 +88,6 @@ void refine_do_cell(ControlContext *ctx) {
 		return;
 	}
 	
-	/* Initialise the simplex */
-	s[0].dax = 0.0;	s[0].dbx = 0.0;	s[0].dcx = 0.0;
-	s[0].day = 0.0;	s[0].dby = 0.0;	s[0].dcy = 0.0;
-	s[0].daz = 0.0;	s[0].dbz = 0.0;	s[0].dcz = 0.0;
-	memcpy(&s[1], &s[0], sizeof(SimplexVertex));	s[1].dax = delta;
-	memcpy(&s[2], &s[0], sizeof(SimplexVertex));	s[2].dbx = delta;
-	memcpy(&s[3], &s[0], sizeof(SimplexVertex));	s[3].dcx = delta;
-	memcpy(&s[4], &s[0], sizeof(SimplexVertex));	s[4].day = delta;
-	memcpy(&s[5], &s[0], sizeof(SimplexVertex));	s[5].dby = delta;
-	memcpy(&s[6], &s[0], sizeof(SimplexVertex));	s[6].dcy = delta;
-	memcpy(&s[7], &s[0], sizeof(SimplexVertex));	s[7].daz = delta;
-	memcpy(&s[8], &s[0], sizeof(SimplexVertex));	s[8].dbz = delta;
-	memcpy(&s[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
-	
 	/* Create the table of indicies and deviations */
 	nf = 0;
 	for ( i=0; i<ctx->images->n_images; i++ ) {
@@ -100,7 +113,7 @@ void refine_do_cell(ControlContext *ctx) {
 		
 		image = &ctx->images->images[i];
 		
-		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
+		for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) {
 		
 			ImageFeature *rf;
 			double dix, diy, dx, dy;
@@ -117,7 +130,7 @@ void refine_do_cell(ControlContext *ctx) {
 			/* Determine the difference vector */
 			dix = rf->partner->x - rf->x;
 			diy = rf->partner->y - rf->y;
-			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy);
+			//printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy);
 			
 			old_x = rf->partner->x;
 			old_y = rf->partner->y;
@@ -127,7 +140,7 @@ void refine_do_cell(ControlContext *ctx) {
 			mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
 			rf->partner->x = old_x;
 			rf->partner->y = old_y;
-			printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9);
+			//printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9);
 			
 			d[f].dx = dlx;
 			d[f].dy = dly;
@@ -140,27 +153,41 @@ void refine_do_cell(ControlContext *ctx) {
 	}
 	assert( f == nf );
 	
+	/* Initialise the simplex */
+	delta = 0.01e9;
+	s[0].dax = 0.0;	s[0].dbx = 0.0;	s[0].dcx = 0.0;
+	s[0].day = 0.0;	s[0].dby = 0.0;	s[0].dcy = 0.0;
+	s[0].daz = 0.0;	s[0].dbz = 0.0;	s[0].dcz = 0.0;
+	memcpy(&s[1], &s[0], sizeof(SimplexVertex));	s[1].dax = delta;
+	memcpy(&s[2], &s[0], sizeof(SimplexVertex));	s[2].dbx = delta;
+	memcpy(&s[3], &s[0], sizeof(SimplexVertex));	s[3].dcx = delta;
+	memcpy(&s[4], &s[0], sizeof(SimplexVertex));	s[4].day = delta;
+	memcpy(&s[5], &s[0], sizeof(SimplexVertex));	s[5].dby = delta;
+	memcpy(&s[6], &s[0], sizeof(SimplexVertex));	s[6].dcy = delta;
+	memcpy(&s[7], &s[0], sizeof(SimplexVertex));	s[7].daz = delta;
+	memcpy(&s[8], &s[0], sizeof(SimplexVertex));	s[8].dbz = delta;
+	memcpy(&s[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
+	
 	/* Find the least favourable vertex of the simplex */
 	v_worst = 0;
-	fom_worst = 0;
+	fom_worst = 0.0;
 	for ( i=0; i<10; i++ ) {
 	
-		double fom = 0;
-		int j;
+		double fom;
 		
-		for ( j=0; j<nf; j++ ) {
+		fom = refine_mean_dev(d, nf, s, i);
 		
-			double xdev, ydev, zdev;
-			
-			xdev = 0;
-			ydev = 0;
-			zdev = 0;
-			
-			fom += sqrt(xdev*xdev + ydev*ydev + zdev*zdev);
-			
+		printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, (fom/nf)/1e9);
+		if ( fom > fom_worst ) {
+			fom_worst = fom;
+			v_worst = i;
 		}
 		
 	}
+	printf("The worst vertex is number %i\n", v_worst);
+	
+	/* Reflect this vertex across the opposite face */
+	refine_simplex_reflect(s, v_worst);
 	
 	ctx->images->images[ctx->dw->cur_image].rflist = NULL;
 	reproject_lattice_changed(ctx);