aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/displaywindow.c4
-rw-r--r--src/ipr.c208
-rw-r--r--src/ipr.h2
-rw-r--r--src/reflections.c57
-rw-r--r--src/reflections.h19
-rw-r--r--src/reproject.c9
-rw-r--r--src/utils.c4
-rw-r--r--src/utils.h1
8 files changed, 195 insertions, 109 deletions
diff --git a/src/displaywindow.c b/src/displaywindow.c
index eedad54..6a3adfe 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -484,7 +484,6 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
glEnd();
/* x, y, z pointers */
-#if 0
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
DRAW_SHORT_POINTER
@@ -499,8 +498,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
glRotatef(-90.0, 0.0, 1.0, 0.0);
DRAW_SHORT_POINTER
glPopMatrix();
-#endif
-
+
/* Tilt axis */
glPushMatrix();
/* Images rotate clockwise by omega to put tilt axis at +x,
diff --git a/src/ipr.c b/src/ipr.c
index 47b57e2..dd6ce59 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -27,64 +27,6 @@
#include "ipr.h"
#include "displaywindow.h"
-#if 0
-static int ipr_random_image(ControlContext *ctx) {
- double n;
- n = (double)random()/RAND_MAX;
- n *= ctx->n_images;
- return (int)n;
-}
-#endif
-
-static Basis *ipr_choose_random_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;
-
- /* 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);
-
- 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_MARKER);
-
- basis = malloc(sizeof(Basis));
-
- basis->a.x = 1.842e9; basis->a.y = 0.0; basis->a.z = 0.0;
- basis->b.x = 0.0; basis->b.y = 1.842e9; basis->b.z = 0.0;
- basis->c.x = 0.0; basis->c.y = 0.0; basis->c.z = 1.842e9;
-
- return basis;
-
-}
-
static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
int n_indexed, n_counted;
@@ -97,7 +39,7 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
if ( cur->type == REFLECTION_NORMAL ) {
- /* Can this basis approximately account for this reflection? */
+ /* Can this basis "approximately" account for this reflection? */
double det;
double a11, a12, a13, a21, a22, a23, a31, a32, a33;
double h, k, l;
@@ -135,40 +77,126 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
}
-static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
+static Basis *ipr_find_basis(ControlContext *ctx) {
- Basis *basis;
- Basis *best_basis;
- double best_efom;
- int i;
+ 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;
- /* Track the best basis throughout the whole procedure */
- best_basis = NULL; best_efom = 0;
+ /* 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);
- printf("IP: Finding initial basis using eFOM...\n");
- for ( i=1; i<=30; i++ ) {
-
- double efom;
+ /* 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;
- basis = ipr_choose_random_basis(ctx);
- efom = ipr_efom(ctx->reflectionctx, basis);
+ do {
+
+ Reflection *check;
+ int orders, j, lfom;
+ double tol;
+
+ accept = 1;
- if ( (efom > best_efom) || !best_basis ) {
- if ( best_basis ) free(best_basis);
- best_efom = efom;
- best_basis = basis;
- printf("IP: %3i: eFOM = %7.3f %%\n", i, efom*100);
+ /* 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: Initial basis:\n");
- printf("IP: a = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->a.x, basis->a.y, basis->a.z, basis->a.modulus);
- printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus);
- printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus);
+ printf("IP: eFOM = %7.3f %%\n", ipr_efom(ctx->reflectionctx, basis)*100);
return basis;
+
}
/* Generate list of reflections to analyse */
@@ -191,14 +219,14 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
}
max_res = max_res / 2;
} else {
- max_res = 2e10; //works until I put some more in the reflect.cache header
+ max_res = 1e10; //works until I put some more in the reflect.cache header
}
ordered = reflection_init();
- max_order_a = max_res/basis->a.modulus;
- max_order_b = max_res/basis->b.modulus;
- max_order_c = max_res/basis->c.modulus;
+ max_order_a = max_res/modulus(basis->a.x, basis->a.y, basis->a.z);
+ max_order_b = max_res/modulus(basis->b.x, basis->b.y, basis->b.z);
+ max_order_c = max_res/modulus(basis->c.x, basis->c.y, basis->c.z);
for ( h=-max_order_a; h<=max_order_a; h++ ) {
for ( k=-max_order_b; k<=max_order_b; k++ ) {
for ( l=-max_order_c; l<=max_order_c; l++ ) {
@@ -217,6 +245,7 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
}
if ( (h!=0) || (k!=0) || (l!=0) ) {
// reflection_add(ctx->reflectionctx, x, y, z, 1, REFLECTION_GENERATED);
+ reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
}
}
@@ -256,10 +285,11 @@ int ipr_refine(ControlContext *ctx) {
Basis *basis;
int finished;
- srand(time(NULL));
-
- basis = ipr_choose_initial_basis(ctx);
- if ( !basis ) return -1;
+ basis = ipr_find_basis(ctx);
+ if ( !basis ) {
+ printf("IP: Unable to find basis\n");
+ return -1;
+ }
ctx->ipr_basis = basis;
ctx->ipr_lat = ipr_generate(ctx, basis);
diff --git a/src/ipr.h b/src/ipr.h
index 3de3724..e511756 100644
--- a/src/ipr.h
+++ b/src/ipr.h
@@ -23,8 +23,6 @@ typedef struct {
double y;
double z;
- double modulus; /* Convenience */
-
} Vector;
typedef struct basis_struct {
diff --git a/src/reflections.c b/src/reflections.c
index 63eaa54..a1f6a2c 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -102,6 +102,7 @@ void reflection_free(ReflectionContext *reflectionctx) {
Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) {
Reflection *new_reflection;
+ Reflection *nearest;
if ( reflectionctx->list_capped ) return NULL;
@@ -110,6 +111,10 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
fprintf(stderr, "No further reflections will be stored. Go and fix the peak detection.\n");
reflectionctx->list_capped = 1;
}
+
+ nearest = reflection_find_nearest_type(reflectionctx, x, y, z, type);
+ if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL;
+
reflectionctx->n_reflections++;
new_reflection = malloc(sizeof(Reflection));
@@ -119,6 +124,7 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
new_reflection->z = z;
new_reflection->intensity = intensity;
new_reflection->type = type;
+ new_reflection->found = 0;
reflectionctx->last_reflection->next = new_reflection;
reflectionctx->last_reflection = new_reflection;
@@ -248,7 +254,56 @@ Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x,
reflection = reflectionctx->reflections;
while ( reflection ) {
- if ( reflection->type == REFLECTION_NORMAL ) {
+ if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) {
+ double mod;
+ mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
+ if ( mod < max ) {
+ max = mod;
+ best = reflection;
+ }
+ }
+ 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) {
+
+ double max = +INFINITY;
+ Reflection *reflection;
+ Reflection *best = NULL;
+
+ reflection = reflectionctx->reflections;
+ while ( reflection ) {
+ if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) {
+ double mod;
+ mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
+ if ( (mod < max) && (mod >= min_distance) ) {
+ max = mod;
+ best = reflection;
+ }
+ }
+ reflection = reflection->next;
+ };
+
+ if ( best ) best->found = 1;
+ return best;
+
+}
+
+Reflection *reflection_find_nearest_type(ReflectionContext *reflectionctx, double x, double y, double z, ReflectionType type) {
+
+ double max = +INFINITY;
+ Reflection *reflection;
+ Reflection *best = NULL;
+
+ reflection = reflectionctx->reflections;
+ while ( reflection ) {
+ if ( reflection->type == type ) {
double mod;
mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
if ( mod < max ) {
diff --git a/src/reflections.h b/src/reflections.h
index 824fcba..9dd281a 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -29,16 +29,17 @@ typedef enum {
typedef struct reflection_struct {
- double x;
- double y;
- double z;
- double intensity;
+ double x;
+ double y;
+ double z;
+ double intensity;
- signed int h;
- signed int k;
- signed int l;
+ signed int h;
+ signed int k;
+ signed int l;
- ReflectionType type;
+ ReflectionType type;
+ int found; /* This reflection has been used in the seed-finding process */
struct reflection_struct *next; /* MUST BE LAST in order for caching to work */
@@ -68,6 +69,8 @@ extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *
extern double reflection_largest_g(ReflectionContext *reflectionctx);
extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta);
extern Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z);
+extern Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, double x, double y, double z, double min_distance);
+extern Reflection *reflection_find_nearest_type(ReflectionContext *reflectionctx, double x, double y, double z, ReflectionType type);
#endif /* REFLECTION_H */
diff --git a/src/reproject.c b/src/reproject.c
index 0971b25..ed359fc 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -33,12 +33,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
tilt = deg2rad(image.tilt);
omega = deg2rad(image.omega);
- refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
-
- i = 0;
-
- reflection = rctx->reflections;
-
/* Calculate the (normalised) incident electron wavevector
* n is rotated "into" the reconstruction, so only one omega step. */
nxt = 0.0; nyt = 0.0; nzt = 1.0;
@@ -58,6 +52,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
ux = uxt*cos(-omega) + uyt*-sin(omega); uy = -uxt*sin(omega) + uyt*cos(omega); uz = uzt;
//reflection_add(ctx->reflectionctx, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
+ refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
+ reflection = rctx->reflections;
+ i = 0;
do {
double xl, yl, zl;
diff --git a/src/utils.c b/src/utils.c
index eaaa036..4618b44 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -39,6 +39,10 @@ double modulus(double x, double y, double z) {
return sqrt(x*x + y*y + z*z);
}
+double distance3d(double x1, double y1, double z1, double x2, double y2, double z2) {
+ return modulus(x1-x2, y1-y2, z1-z2);
+}
+
/* Angle between two vectors. Answer in radians */
double angle_between(double x1, double y1, double z1, double x2, double y2, double z2) {
diff --git a/src/utils.h b/src/utils.h
index 99fc748..1d268e6 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -27,6 +27,7 @@ extern double angle_between(double x1, double y1, double z1, double x2, double y
extern double angle_between_d(double x1, double y1, double z1, double x2, double y2, double z2);
extern double lambda(double voltage);
extern void matrix_renormalise(gsl_matrix *m);
+extern double distance3d(double x1, double y1, double z1, double x2, double y2, double z2);
#define rad2deg(a) ((a)*180/M_PI)
#define deg2rad(a) ((a)*M_PI/180)
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-04 21:20:04 +0000