diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-28 23:52:58 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-28 23:52:58 +0000 |
commit | 555bc4f60c845bf47aee94e4b7963382838c9f57 (patch) | |
tree | aa012672259995addc1d01420aee2df608bd83d8 | |
parent | 12271165c0948536f9b34603432ade4953c97b4e (diff) |
Move basis stuff to separate file
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@136 bf6ca9ba-c028-0410-8290-897cf20841d1
-rw-r--r-- | Makefile.am | 2 | ||||
-rw-r--r-- | src/Makefile.am | 2 | ||||
-rw-r--r-- | src/basis.c | 252 | ||||
-rw-r--r-- | src/basis.h | 38 | ||||
-rw-r--r-- | src/displaywindow.c | 1 | ||||
-rw-r--r-- | src/ipr.c | 175 | ||||
-rw-r--r-- | src/ipr.h | 14 | ||||
-rw-r--r-- | src/mapping.c | 2 | ||||
-rw-r--r-- | src/reflections.c | 5 |
9 files changed, 297 insertions, 194 deletions
diff --git a/Makefile.am b/Makefile.am index bd1c637..4c3182e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,4 +1,4 @@ EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src/control.h src/readpng.h src/mrc.h src/imagedisplay.h src/main.h \ data/displaywindow.ui src/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h src/itrans-threshold.h \ - src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h + src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h src/basis.h SUBDIRS = src data diff --git a/src/Makefile.am b/src/Makefile.am index 59c995c..bdfe197 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ bin_PROGRAMS = dtr dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c cache.c \ - itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.c + itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" diff --git a/src/basis.c b/src/basis.c new file mode 100644 index 0000000..a7ff978 --- /dev/null +++ b/src/basis.c @@ -0,0 +1,252 @@ +/* + * basis.c + * + * Find approximate lattices to feed various procedures + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <math.h> +#include <stdio.h> +#include <stdlib.h> + +#include "utils.h" +#include "reflections.h" +#include "basis.h" + +static double basis_efom(ReflectionContext *rctx, Basis *basis) { + + int n_indexed, n_counted; + Reflection *cur; + + cur = rctx->reflections; + n_indexed = 0; + n_counted = 0; + while ( cur ) { + + if ( cur->type == REFLECTION_NORMAL ) { + + /* Can this basis "approximately" account for this reflection? */ + double det; + double a11, a12, a13, a21, a22, a23, a31, a32, a33; + double h, k, l; + + /* Set up the coordinate transform from hkl to xyz */ + a11 = basis->a.x; a12 = basis->a.y; a13 = basis->a.z; + a21 = basis->b.x; a22 = basis->b.y; a23 = basis->b.z; + a31 = basis->c.x; a32 = basis->c.y; a33 = basis->c.z; + + /* Invert the matrix to get hkl from xyz */ + det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31); + h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z) / det; + k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det; + l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det; + + /* Calculate the deviations in terms of |a|, |b| and |c| */ + h = fabs(h); k = fabs(k); l = fabs(l); + h -= floor(h); k -= floor(k); l -= floor(l); + if ( h == 1.0 ) h = 0.0; + if ( k == 1.0 ) k = 0.0; + if ( l == 1.0 ) l = 0.0; + + /* Define "approximately" here. Circle in basis space becomes an ellipsoid in reciprocal space */ + if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++; + + n_counted++; + + } + + cur = cur->next; + + } + + return (double)n_indexed / n_counted; + +} + +static int basis_lfom(ControlContext *ctx, double vx, double vy, double vz) { + + Reflection *tcentre; + int lfom; + double tol; + int j; + + lfom = 0; + tol = modulus(vx, vy, vz)/10.0; + tcentre = ctx->reflectionctx->reflections; + do { + + for ( j=-20; j<=20; j++ ) { + + Reflection *check; + + check = reflection_find_nearest(ctx->reflectionctx, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j); + if ( check && (distance3d(check->x, check->y, check->z, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j) < tol) ) { + lfom++; + } + + } + + tcentre = tcentre->next; + + } while ( tcentre ); + + return lfom; + +} + +static ReflectionContext *basis_find_seeds(ControlContext *ctx) { + + double tilt_min; + double tilt_max; + double tilt_mid; + ImageRecord *imagerecord; + double x_temp, y_temp, z_temp; + double scale; + double x, y, z; + Reflection *centre; + int i; + ReflectionContext *seeds; + + seeds = reflection_init(); + + /* Locate the 'plane' in the middle of the "wedge". + * This whole procedure assumes there is just one tilt axis. */ + tilt_min = control_min_tilt(ctx); + tilt_max = control_max_tilt(ctx); + tilt_mid = tilt_min + (tilt_max-tilt_min)/2; + imagerecord = control_image_nearest_tilt(ctx, tilt_mid); + + /* Apply the last two steps of the mapping transform to get the direction from the origin + * towards the middle of the wedge */ + x_temp = 0.0; + y_temp = cos(deg2rad(imagerecord->tilt)); + z_temp = -sin(deg2rad(imagerecord->tilt)); + x = x_temp*cos(-deg2rad(imagerecord->omega)) + y_temp*sin(-deg2rad(imagerecord->omega)); + y = -x_temp*sin(-deg2rad(imagerecord->omega)) + y_temp*cos(-deg2rad(imagerecord->omega)); + z = z_temp; + + /* Find the point in the middle of the "wedge" */ + scale = reflection_largest_g(ctx->reflectionctx)/4; + x *= scale; + y *= scale; + z *= scale; + reflection_add(ctx->reflectionctx, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2); + + /* Find an "origin" reflection */ + centre = reflection_find_nearest(ctx->reflectionctx, x, y, z); + if ( !centre ) return NULL; + centre->found = 1; + reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED); + + for ( i=1; i<=10; i++ ) { + + Reflection *vector; + int accept; + double vx, vy, vz; + + do { + + Reflection *check; + int lfom; + + accept = 1; + + /* Find a "candidate vector" reflection */ + vector = reflection_find_nearest_longer(ctx->reflectionctx, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */ + if ( !vector ) { + printf("BS: Couldn't find enough seeds\n"); + return NULL; + } + vector->found = 1; + + /* Get vector components (not the coordinates the vector was calculated from!) */ + vx = vector->x - centre->x; + vy = vector->y - centre->y; + vz = vector->z - centre->z; + + /* Proximity test */ + check = reflection_find_nearest_type(ctx->reflectionctx, vx, vy, vz, REFLECTION_NORMAL); + if ( check ) { + if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) { + /* Too close to another seed */ + accept = 0; + continue; + } + } + + /* lFOM test */ + lfom = basis_lfom(ctx, vx, vy, vz); + if ( lfom < 1 ) { + accept = 0; + continue; + } + printf("lfom=%i\n", lfom); + + } while ( !accept ); + + reflection_add(seeds, vx, vy, vz, 1.0, REFLECTION_NORMAL); + reflection_add(ctx->reflectionctx, vx, vy, vz, 1.0, REFLECTION_MARKER); + + } + + return seeds; + +} + +Basis *basis_find(ControlContext *ctx) { + + Basis *basis; + ReflectionContext *seeds; + Reflection *ref; + int i; + + /* Get the shortlist of seeds */ + seeds = basis_find_seeds(ctx); + + /* Assemble the seeds into a basis */ + basis = malloc(sizeof(Basis)); + ref = seeds->reflections; + for ( i=1; i<=3; i++ ) { + + double vx, vy, vz; + + vx = ref->x; + vy = ref->y; + vz = ref->z; + + switch ( i ) { + case 1 : { + basis->a.x = vx; + basis->a.y = vy; + basis->a.z = vz; + } + case 2 : { + basis->b.x = vx; + basis->b.y = vy; + basis->b.z = vz; + } + case 3 : { + basis->c.x = vx; + basis->c.y = vy; + basis->c.z = vz; + } + } + + ref = ref->next; + + } + + printf("BS: eFOM = %7.3f %%\n", basis_efom(ctx->reflectionctx, basis)*100); + + return basis; + +} + diff --git a/src/basis.h b/src/basis.h new file mode 100644 index 0000000..1d12be1 --- /dev/null +++ b/src/basis.h @@ -0,0 +1,38 @@ +/* + * basis.h + * + * Find approximate lattices to feed various procedures + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifndef BASIS_H +#define BASIS_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "control.h" + +typedef struct { + + double x; + double y; + double z; + +} Vector; + +typedef struct basis_struct { + Vector a; + Vector b; + Vector c; +} Basis; + +extern Basis *basis_find(ControlContext *ctx); + +#endif /* BASIS_H */ + diff --git a/src/displaywindow.c b/src/displaywindow.c index 6a3adfe..867b2cb 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -30,6 +30,7 @@ #include "main.h" #include "ipr.h" #include "displaywindow.h" +#include "basis.h" enum { DW_ORTHO, @@ -26,178 +26,7 @@ #include "reproject.h" #include "ipr.h" #include "displaywindow.h" - -static double ipr_efom(ReflectionContext *rctx, Basis *basis) { - - int n_indexed, n_counted; - Reflection *cur; - - cur = rctx->reflections; - n_indexed = 0; - n_counted = 0; - while ( cur ) { - - if ( cur->type == REFLECTION_NORMAL ) { - - /* Can this basis "approximately" account for this reflection? */ - double det; - double a11, a12, a13, a21, a22, a23, a31, a32, a33; - double h, k, l; - - /* Set up the coordinate transform from hkl to xyz */ - a11 = basis->a.x; a12 = basis->a.y; a13 = basis->a.z; - a21 = basis->b.x; a22 = basis->b.y; a23 = basis->b.z; - a31 = basis->c.x; a32 = basis->c.y; a33 = basis->c.z; - - /* Invert the matrix to get hkl from xyz */ - det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31); - h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z) / det; - k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det; - l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det; - - /* Calculate the deviations in terms of |a|, |b| and |c| */ - h = fabs(h); k = fabs(k); l = fabs(l); - h -= floor(h); k -= floor(k); l -= floor(l); - if ( h == 1.0 ) h = 0.0; - if ( k == 1.0 ) k = 0.0; - if ( l == 1.0 ) l = 0.0; - - /* Define "approximately" here. Circle in basis space becomes an ellipsoid in reciprocal space */ - if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++; - - n_counted++; - - } - - cur = cur->next; - - } - - return (double)n_indexed / n_counted; - -} - -static Basis *ipr_find_basis(ControlContext *ctx) { - - Basis *basis; - double tilt_min; - double tilt_max; - double tilt_mid; - ImageRecord *imagerecord; - double x_temp, y_temp, z_temp; - double scale; - double x, y, z; - Reflection *centre; - int i; - - /* Locate the 'plane' in the middle of the "wedge". - * This whole procedure assumes there is just one tilt axis. */ - tilt_min = control_min_tilt(ctx); - tilt_max = control_max_tilt(ctx); - tilt_mid = tilt_min + (tilt_max-tilt_min)/2; - imagerecord = control_image_nearest_tilt(ctx, tilt_mid); - - /* Apply the last two steps of the mapping transform to get the direction from the origin - * towards the middle of the wedge */ - x_temp = 0.0; - y_temp = cos(deg2rad(imagerecord->tilt)); - z_temp = -sin(deg2rad(imagerecord->tilt)); - x = x_temp*cos(-deg2rad(imagerecord->omega)) + y_temp*sin(-deg2rad(imagerecord->omega)); - y = -x_temp*sin(-deg2rad(imagerecord->omega)) + y_temp*cos(-deg2rad(imagerecord->omega)); - z = z_temp; - - /* Find the point in the middle of the "wedge" */ - scale = reflection_largest_g(ctx->reflectionctx)/4; - x *= scale; - y *= scale; - z *= scale; - reflection_add(ctx->reflectionctx, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2); - - /* Find an "origin" reflection */ - centre = reflection_find_nearest(ctx->reflectionctx, x, y, z); - if ( !centre ) return NULL; - reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED); - - basis = malloc(sizeof(Basis)); - for ( i=1; i<=3; i++ ) { - - Reflection *vector; - int accept; - double vx, vy, vz; - - do { - - Reflection *check; - int orders, j, lfom; - double tol; - - accept = 1; - - /* Find a "candidate vector" reflection */ - vector = reflection_find_nearest_longer(ctx->reflectionctx, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */ - if ( !vector ) { - printf("IP: Couldn't find enough seeds\n"); - return NULL; - } - - /* Get vector components (not the coordinates the vector was calculated from!) */ - vx = vector->x - centre->x; - vy = vector->y - centre->y; - vz = vector->z - centre->z; - - /* Proximity test */ - check = reflection_find_nearest_type(ctx->reflectionctx, vx, vy, vz, REFLECTION_MARKER); - if ( check ) { - if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) { - /* Too close to another seed */ - accept = 0; - continue; - } - } - - /* Line FoM test */ - orders = 100e9 / modulus(vx, vy, vz); /* Go out by +/- this number of iterations. Rounding error ignored */ - lfom = 0; - tol = modulus(vx, vy, vz)/10.0; - for ( j=-orders; j<=orders; j++ ) { - check = reflection_find_nearest(ctx->reflectionctx, centre->x+vx*j, centre->y+vy*j, centre->z+vz*j); - if ( check && (distance3d(check->x, check->y, check->z, centre->x+vx*j, centre->y+vy*j, centre->z+vz*j) < tol) ) { - lfom++; - } - } - if ( lfom < 1 ) { - accept = 0; - continue; - } - - } while ( !accept ); - reflection_add(ctx->reflectionctx, vx, vy, vz, 1.0, REFLECTION_MARKER); - - switch ( i ) { - case 1 : { - basis->a.x = vx; - basis->a.y = vy; - basis->a.z = vz; - } - case 2 : { - basis->b.x = vx; - basis->b.y = vy; - basis->b.z = vz; - } - case 3 : { - basis->c.x = vx; - basis->c.y = vy; - basis->c.z = vz; - } - } - - } - - printf("IP: eFOM = %7.3f %%\n", ipr_efom(ctx->reflectionctx, basis)*100); - - return basis; - -} +#include "basis.h" /* Generate list of reflections to analyse */ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) { @@ -285,7 +114,7 @@ int ipr_refine(ControlContext *ctx) { Basis *basis; int finished; - basis = ipr_find_basis(ctx); + basis = basis_find(ctx); if ( !basis ) { printf("IP: Unable to find basis\n"); return -1; @@ -17,20 +17,6 @@ #ifndef IPR_H #define IPR_H -typedef struct { - - double x; - double y; - double z; - -} Vector; - -typedef struct basis_struct { - Vector a; - Vector b; - Vector c; -} Basis; - #include "control.h" extern int ipr_refine(ControlContext *ctx); diff --git a/src/mapping.c b/src/mapping.c index cac1ce5..fdce9ea 100644 --- a/src/mapping.c +++ b/src/mapping.c @@ -32,7 +32,7 @@ ReflectionContext *mapping_create(ControlContext *ctx) { for ( i=0; i<ctx->n_images; i++ ) { itrans_process_image(&ctx->images[i], ctx); } - printf("done\n"); + printf("done: %i 'measured' reflections\n", ctx->reflectionctx->n_reflections); return rctx; diff --git a/src/reflections.c b/src/reflections.c index a1f6a2c..3fa1b6d 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -254,7 +254,7 @@ Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, reflection = reflectionctx->reflections; while ( reflection ) { - if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) { + if ( reflection->type == REFLECTION_NORMAL ) { double mod; mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z); if ( mod < max ) { @@ -265,10 +265,8 @@ Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, reflection = reflection->next; }; - if ( best ) best->found = 1; return best; - } Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, double x, double y, double z, double min_distance) { @@ -290,7 +288,6 @@ Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, dou reflection = reflection->next; }; - if ( best ) best->found = 1; return best; } |