diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-04-16 15:29:40 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-04-16 15:29:40 +0000 |
commit | b705768d46b39d2f7d38ec92db5086c0b8b3671d (patch) | |
tree | 3cdf33e033b8a58dc71380380c62d9d78c795517 /src/refine.c | |
parent | 5a5e330a0477f9db2f5af85637e9190be0e47fcc (diff) |
Refinement progress
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@272 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src/refine.c')
-rw-r--r-- | src/refine.c | 112 |
1 files changed, 80 insertions, 32 deletions
diff --git a/src/refine.c b/src/refine.c index 7bb2585..ac2ba38 100644 --- a/src/refine.c +++ b/src/refine.c @@ -34,6 +34,7 @@ typedef struct { double dax; double dbx; double dcx; double day; double dby; double dcy; double daz; double dbz; double dcz; + int does_nothing; /* If non-zero, this vertex doesn't appear to affect the FoM */ } SimplexVertex; typedef struct { @@ -127,20 +128,14 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { } -#if 0 -static void refine_display_simplex(SimplexVertex *s) { +static void refine_display_simplex(SimplexVertex sv) { - int i; - - for ( i=0; i<10; i++ ) { - printf("Vertex %i: %8f %8f %8f %8f %8f %8f %8f %8f %8f\n", - i, s[i].dax/1e9, s[i].day/1e9, s[i].daz/1e9, - s[i].dbx/1e9, s[i].dby/1e9, s[i].dbz/1e9, - s[i].dcx/1e9, s[i].dcy/1e9, s[i].dcz/1e9); - } + printf("%8f %8f %8f %8f %8f %8f %8f %8f %8f\n", + sv.dax/1e9, sv.day/1e9, sv.daz/1e9, + sv.dbx/1e9, sv.dby/1e9, sv.dbz/1e9, + sv.dcx/1e9, sv.dcy/1e9, sv.dcz/1e9); } -#endif /* Expand the simplex across from vertex v_worst by factor 'fac'. * fac = -1 is a reflection @@ -148,21 +143,62 @@ static void refine_display_simplex(SimplexVertex *s) { */ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) { - int v_face; /* A simplex vertex which is on the face being reflected across */ + SimplexVertex centre; + int i, nv; + + /* Average the coordinates of the non-worst vertices to + * get the centre of the opposite face. */ + centre.dax = 0.0; centre.day = 0.0; centre.daz = 0.0; + centre.dbx = 0.0; centre.dby = 0.0; centre.dbz = 0.0; + centre.dcx = 0.0; centre.dcy = 0.0; centre.dcz = 0.0; + nv = 0; + for ( i=0; i<10; i++ ) { + if ( i != v_worst ) { + centre.dax += s[i].dax; centre.dax += s[i].day; centre.dax += s[i].daz; + centre.dax += s[i].dbx; centre.dax += s[i].dby; centre.dax += s[i].dbz; + centre.dax += s[i].dcx; centre.dax += s[i].dcy; centre.dax += s[i].dcz; + nv++; + } + } + centre.dax /= nv; centre.day /= nv; centre.daz /= nv; + centre.dbx /= nv; centre.dby /= nv; centre.dbz /= nv; + centre.dcx /= nv; centre.dcy /= nv; centre.dcz /= nv; + +// printf(" Midpoint: "); +// refine_display_simplex(centre); + +// printf("Before transformation: "); +// refine_display_simplex(s[v_worst]); - v_face = 0; - if ( v_worst == 0 ) v_face = 1; + /* Do the transformation */ + s[v_worst].dax = centre.dax + fac * (s[v_worst].dax - centre.dax); + s[v_worst].day = centre.day + fac * (s[v_worst].day - centre.day); + s[v_worst].daz = centre.daz + fac * (s[v_worst].daz - centre.daz); + s[v_worst].dbx = centre.dbx + fac * (s[v_worst].dbx - centre.dbx); + s[v_worst].dby = centre.dby + fac * (s[v_worst].dby - centre.dby); + s[v_worst].dbz = centre.dbz + fac * (s[v_worst].dbz - centre.dbz); + s[v_worst].dcx = centre.dcx + fac * (s[v_worst].dcx - centre.dcx); + s[v_worst].dcy = centre.dcy + fac * (s[v_worst].dcy - centre.dcy); + s[v_worst].dcz = centre.dcz + fac * (s[v_worst].dcz - centre.dcz); + +// printf(" After transformation: "); +// refine_display_simplex(s[v_worst]); + +} - s[v_worst].dax = s[v_face].dax + fac * (s[v_worst].dax - s[v_face].dax); - s[v_worst].day = s[v_face].day + fac * (s[v_worst].day - s[v_face].dax); - s[v_worst].daz = s[v_face].daz + fac * (s[v_worst].daz - s[v_face].dax); - s[v_worst].dbx = s[v_face].dbx + fac * (s[v_worst].dbx - s[v_face].dbx); - s[v_worst].dby = s[v_face].dby + fac * (s[v_worst].dby - s[v_face].dbx); - s[v_worst].dbz = s[v_face].dbz + fac * (s[v_worst].dbz - s[v_face].dbx); - s[v_worst].dcx = s[v_face].dcx + fac * (s[v_worst].dcx - s[v_face].dcx); - s[v_worst].dcy = s[v_face].dcy + fac * (s[v_worst].dcy - s[v_face].dcx); - s[v_worst].dcz = s[v_face].dcz + fac * (s[v_worst].dcz - s[v_face].dcx); +/* Contract vertex v of simplex s towards vertex v_best */ +static void refine_simplex_contract(SimplexVertex *s, int v, int v_best) { + s[v].dax = s[v_best].dax + 0.5 * (s[v].dax - s[v_best].dax); + s[v].day = s[v_best].day + 0.5 * (s[v].day - s[v_best].day); + s[v].daz = s[v_best].daz + 0.5 * (s[v].daz - s[v_best].daz); + s[v].dbx = s[v_best].dbx + 0.5 * (s[v].dbx - s[v_best].dbx); + s[v].dby = s[v_best].dby + 0.5 * (s[v].dby - s[v_best].dby); + s[v].dbz = s[v_best].dbz + 0.5 * (s[v].dbz - s[v_best].dbz); + s[v].dcx = s[v_best].dcx + 0.5 * (s[v].dcx - s[v_best].dcx); + s[v].dcy = s[v_best].dcy + 0.5 * (s[v].dcy - s[v_best].dcy); + s[v].dcz = s[v_best].dcz + 0.5 * (s[v].dcz - s[v_best].dcz); + } static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug) { @@ -177,14 +213,15 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug fom_best = 100e9; v_second_worst = 0; fom_second_worst = 0.0; + if ( debug ) printf("Vertex FoM/nm^-1 DoesNothing\n"); for ( i=0; i<10; i++ ) { double fom; fom = refine_mean_dev(d, nf, s, i); - if ( debug ) printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9); - if ( fom > fom_worst ) { + if ( debug ) printf("%6i %8f %s\n", i, fom/1e9, s[i].does_nothing?"*":" "); + if ( (s[i].does_nothing == 0) && (fom > fom_worst) ) { v_second_worst = v_worst; fom_second_worst = fom_worst; fom_worst = fom; @@ -210,7 +247,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug /* It's worse than before. Contract in 1D and see if that helps. */ if ( debug ) printf("Worse. Trying a 1D contraction...\n"); - refine_simplex_transform(s, v_worst, 0.5); + refine_simplex_transform(s, v_worst, -0.5); /* Minus puts it back on the original side of the 'good' face */ fom_new_new = refine_mean_dev(d, nf, s, v_worst); if ( debug ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9); if ( fom_new_new > fom_second_worst ) { @@ -219,14 +256,13 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug if ( debug ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); for ( i=0; i<10; i++ ) { - if ( i != v_best ) refine_simplex_transform(s, i, 0.5); + if ( i != v_best ) refine_simplex_contract(s, i, v_best); } } - } else { + } else if ( fom_new < fom_worst ) { - #if 0 /* It's better. Try to expand in this direction */ double fom_new_new; SimplexVertex save; @@ -234,7 +270,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug if ( debug ) printf("This is better. Trying to expand...\n"); save = s[v_worst]; - refine_simplex_transform(s, v_worst, 2); + refine_simplex_transform(s, v_worst, 2.0); /* +ve means stay on this side of the 'good' face */ /* Better? */ fom_new_new = refine_mean_dev(d, nf, s, v_worst); if ( debug ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9); @@ -242,9 +278,15 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug /* "Got too confident" */ s[v_worst] = save; if ( debug ) printf("Got too confident - reverting\n"); - } /* else yay. */ - #endif + } else { + if ( debug ) printf("Better still. Great.\n"); + } + + } else { + printf("No change!\n"); + s[v_worst].does_nothing = 1; + } return fom_worst - fom_best; @@ -341,6 +383,7 @@ double refine_do_cell(ControlContext *ctx) { 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; + s[0].does_nothing = 0; memcpy(&s[1], &s[0], sizeof(SimplexVertex)); s[1].dax = delta; memcpy(&s[2], &s[0], sizeof(SimplexVertex)); s[2].day = delta; memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].daz = delta; @@ -350,6 +393,11 @@ double refine_do_cell(ControlContext *ctx) { memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].dcx = delta; memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dcy = delta; memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; + s[3].does_nothing = 1; + s[6].does_nothing = 1; + s[7].does_nothing = 1; + s[8].does_nothing = 1; + s[9].does_nothing = 1; /* Iterate */ maxiter = 500; |