aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile.am2
-rw-r--r--src/Makefile.am2
-rw-r--r--src/basis.c252
-rw-r--r--src/basis.h38
-rw-r--r--src/displaywindow.c1
-rw-r--r--src/ipr.c175
-rw-r--r--src/ipr.h14
-rw-r--r--src/mapping.c2
-rw-r--r--src/reflections.c5
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,
diff --git a/src/ipr.c b/src/ipr.c
index dd6ce59..56f871c 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -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;
diff --git a/src/ipr.h b/src/ipr.h
index e511756..8fb4862 100644
--- a/src/ipr.h
+++ b/src/ipr.h
@@ -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;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-21 06:46:56 +0000