aboutsummaryrefslogtreecommitdiff
path: root/src/ipr.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/ipr.c')
-rw-r--r--src/ipr.c202
1 files changed, 40 insertions, 162 deletions
diff --git a/src/ipr.c b/src/ipr.c
index 71444af..03824de 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -27,18 +27,6 @@
#include "ipr.h"
#include "displaywindow.h"
-static int ipr_choose_max(Basis *basis) {
-
- if ( (basis->a.modulus >= basis->b.modulus) && (basis->a.modulus >= basis->c.modulus) ) return 1;
- if ( (basis->b.modulus >= basis->a.modulus) && (basis->b.modulus >= basis->c.modulus) ) return 2;
- if ( (basis->c.modulus >= basis->a.modulus) && (basis->c.modulus >= basis->b.modulus) ) return 3;
-
- /* Never reaches here... */
- fprintf(stderr, "IP: Couldn't determine maximum basis vector modulus\n");
- abort();
-
-}
-
#if 0
static int ipr_random_image(ControlContext *ctx) {
double n;
@@ -48,160 +36,47 @@ static int ipr_random_image(ControlContext *ctx) {
}
#endif
-static Reflection *ipr_random_reflection(ReflectionContext *rctx) {
-
- double n;
- int i;
- Reflection *cur;
-
- n = (double)random()/RAND_MAX;
- n *= rctx->n_reflections/2;
- cur = rctx->reflections; /* Jump a maximum of half way in */
- for ( i=1; i<n; i++ ) { /* Skip initial placeholder */
- cur = cur->next;
- }
-
- return cur;
-
-}
-
-/* Try to turn the basis into a more sensible one */
-static void ipr_try_compact(Basis *basis) {
-
- int i;
-
- for ( i=1; i<5; i++ ) {
-
- if ( modulus(basis->a.x - basis->b.x, basis->a.y - basis->b.y, basis->a.z - basis->b.z) < basis->a.modulus ) {
- basis->a.x = basis->a.x - basis->b.x; basis->a.y = basis->a.y - basis->b.y; basis->a.y = basis->a.y - basis->b.y;
- basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
- }
- if ( modulus(basis->a.x - basis->c.x, basis->a.z - basis->c.z, basis->a.z - basis->c.z) < basis->a.modulus ) {
- basis->a.x = basis->a.x - basis->c.x; basis->a.y = basis->a.y - basis->c.y; basis->a.y = basis->a.y - basis->c.y;
- basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
- }
- if ( modulus(basis->a.x + basis->b.x, basis->a.y + basis->b.y, basis->a.z + basis->b.z) < basis->a.modulus ) {
- basis->a.x = basis->a.x + basis->b.x; basis->a.y = basis->a.y + basis->b.y; basis->a.y = basis->a.y + basis->b.y;
- basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
- }
- if ( modulus(basis->a.x + basis->c.x, basis->a.z + basis->c.z, basis->a.z + basis->c.z) < basis->a.modulus ) {
- basis->a.x = basis->a.x + basis->c.x; basis->a.y = basis->a.y + basis->c.y; basis->a.y = basis->a.y + basis->c.y;
- basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
- }
-
- if ( modulus(basis->b.x - basis->a.x, basis->b.y - basis->a.y, basis->b.z - basis->a.z) < basis->b.modulus ) {
- basis->b.x = basis->b.x - basis->a.x; basis->b.y = basis->b.y - basis->a.y; basis->b.y = basis->b.y - basis->a.y;
- basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
- }
- if ( modulus(basis->b.x - basis->c.x, basis->b.y - basis->c.y, basis->b.z - basis->c.z) < basis->b.modulus ) {
- basis->b.x = basis->b.x - basis->c.x; basis->b.y = basis->b.y - basis->c.y; basis->b.y = basis->b.y - basis->c.y;
- basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
- }
- if ( modulus(basis->b.x + basis->a.x, basis->b.y + basis->a.y, basis->b.z + basis->a.z) < basis->b.modulus ) {
- basis->b.x = basis->b.x + basis->a.x; basis->b.y = basis->b.y + basis->a.y; basis->b.y = basis->b.y + basis->a.y;
- basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
- }
- if ( modulus(basis->b.x + basis->c.x, basis->b.y + basis->c.y, basis->b.z + basis->c.z) < basis->b.modulus ) {
- basis->b.x = basis->b.x + basis->c.x; basis->b.y = basis->b.y + basis->c.y; basis->b.y = basis->b.y + basis->c.y;
- basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
- }
-
- if ( modulus(basis->c.x - basis->a.x, basis->c.y - basis->a.z, basis->c.z - basis->a.z) < basis->c.modulus ) {
- basis->c.x = basis->c.x - basis->a.x; basis->c.y = basis->c.y - basis->a.y; basis->c.y = basis->c.y - basis->a.y;
- basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
- }
- if ( modulus(basis->c.x - basis->b.x, basis->c.y - basis->b.y, basis->c.z - basis->b.z) < basis->c.modulus ) {
- basis->c.x = basis->c.x - basis->b.x; basis->c.y = basis->c.y - basis->b.y; basis->c.y = basis->c.y - basis->b.y;
- basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
- }
- if ( modulus(basis->c.x + basis->a.x, basis->c.y + basis->a.z, basis->c.z + basis->a.z) < basis->c.modulus ) {
- basis->c.x = basis->c.x + basis->a.x; basis->c.y = basis->c.y + basis->a.y; basis->c.y = basis->c.y + basis->a.y;
- basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
- }
- if ( modulus(basis->c.x + basis->b.x, basis->c.y + basis->b.y, basis->c.z + basis->b.z) < basis->c.modulus ) {
- basis->c.x = basis->c.x + basis->b.x; basis->c.y = basis->c.y + basis->b.y; basis->c.y = basis->c.y + basis->b.y;
- basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
- }
-
- }
-
-}
-
static Basis *ipr_choose_random_basis(ControlContext *ctx) {
-
- Basis *basis;
- Reflection *reflection;
- int cur_max;
+
+ Basis *basis;
+ double tilt_min;
+ double tilt_max;
+ double tilt_mid;
+ ImageRecord *imagerecord;
+ double nx, ny;
+ double scale;
+ ImageReflection refl;
+ double x, y, z, twotheta;
+ Reflection *centre;
+
+ /* Locate the 'plane' in the middle of the "wedge" */
+ 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_min);
+
+ /* Find the point in the middle of the "wedge" */
+ scale = reflection_largest_g(ctx->reflectionctx);
+ nx = imagerecord->x_centre - scale*sin(imagerecord->omega);
+ ny = imagerecord->y_centre + scale*cos(imagerecord->omega);
+
+ /* Find where this point is in 3D space */
+ refl.x = nx;
+ refl.y = ny;
+ refl.parent = imagerecord;
+ if ( reflection_map_to_space(&refl, &x, &y, &z, &twotheta) ) return NULL;
+
+ /* Find the nearest 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_MARKER);
+ printf("Nearest %f,%f,%f is %f,%f,%f\n", x, y, z, centre->x, centre->y, centre->z);
basis = malloc(sizeof(Basis));
-restart:
- basis->a.x = 0; basis->a.y = 0; basis->a.z = 0; basis->a.modulus = 1e30;
- basis->b.x = 0; basis->b.y = 0; basis->b.z = 0; basis->b.modulus = 1e30;
- basis->c.x = 0; basis->c.y = 0; basis->c.z = 0; basis->c.modulus = 1e30;
- cur_max = 1;
-
- reflection = ipr_random_reflection(ctx->reflectionctx);
-
- do {
- double mod;
- unsigned int angle_works;
-
- mod = modulus(reflection->x, reflection->y, reflection->z);
-
- angle_works = 1;
- if ( basis->a.modulus != 1e30 ) {
- double angle = angle_between_d(reflection->x, reflection->y, reflection->z, basis->a.x, basis->a.y, basis->a.z);
- if ( angle < 20 ) angle_works = 0;
- if ( angle > 160 ) angle_works = 0;
- }
- if ( basis->b.modulus != 1e30 ) {
- double angle = angle_between_d(reflection->x, reflection->y, reflection->z, basis->b.x, basis->b.y, basis->b.z);
- if ( angle < 20 ) angle_works = 0;
- if ( angle > 160 ) angle_works = 0;
- }
- if ( basis->c.modulus != 1e30 ) {
- double angle = angle_between_d(reflection->x, reflection->y, reflection->z, basis->c.x, basis->c.y, basis->c.z);
- if ( angle < 20 ) angle_works = 0;
- if ( angle > 160 ) angle_works = 0;
- }
- if ( (mod > 1e9) && angle_works ) {
- cur_max = ipr_choose_max(basis);
- if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 1) ) {
- basis->a.x = reflection->x;
- basis->a.y = reflection->y;
- basis->a.z = reflection->z;
- basis->a.modulus = mod;
- }
- if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 2) ) {
- basis->b.x = reflection->x;
- basis->b.y = reflection->y;
- basis->b.z = reflection->z;
- basis->b.modulus = mod;
- }
- if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 3) ) {
- basis->c.x = reflection->x;
- basis->c.y = reflection->y;
- basis->c.z = reflection->z;
- basis->c.modulus = mod;
- }
- }
- reflection = reflection->next;
-
- /* End of list? */
- if ( !reflection ) {
- if ( (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) ) {
- goto restart; /* Probably the one circumstance under which this is allowed... */
- }
- }
-
- } while ( reflection || (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) );
-
- ipr_try_compact(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;
+ 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;
@@ -335,6 +210,9 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
if ( ref ) { /* Check it's not NULL */
ref->h = h; ref->k = k; ref->l = l;
}
+ if ( (h!=0) || (k!=0) || (l!=0) ) {
+ // reflection_add(ctx->reflectionctx, x, y, z, 1, REFLECTION_GENERATED);
+ }
}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-12 10:02:15 +0000