From 3ffa0ce286bac2c3c4d9e7d41bd2ca8972d6288a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 30 Mar 2009 12:14:17 +0100 Subject: Tidy up and remove debug --- src/dirax.c | 10 +++---- src/itrans-stat.c | 6 +--- src/reproject.c | 83 +++++++++++++++++++++++++++---------------------------- 3 files changed, 46 insertions(+), 53 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index c4afbce..79b20ea 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -387,11 +387,9 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) ra = ((long long int)random() * (long long int)n); ra /= RAND_MAX; if ( used[ra] == 'U' ) { - printf("%lli - already used\n", ra); continue; } - printf("Selected reflection %lli\n", ra); /* Dig out the correct reflection. A little faffy * because of the linked list */ @@ -401,10 +399,10 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) } /* Limit resolution of reflections */ - if ( ( (ref->x/1e9)*(ref->x/1e9) - + (ref->y/1e9)*(ref->y/1e9) - + (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) ) - continue; +// if ( ( (ref->x/1e9)*(ref->x/1e9) +// + (ref->y/1e9)*(ref->y/1e9) +// + (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) ) +// continue; fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, diff --git a/src/itrans-stat.c b/src/itrans-stat.c index ee41817..0e03ead 100644 --- a/src/itrans-stat.c +++ b/src/itrans-stat.c @@ -441,7 +441,7 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, gsl_matrix_set(p, 2, n, val); n++; if ( n == size ) { - printf("expanding %i->%i\n", size, size*2); + p = itrans_peaksearch_stat_matrix_expand(p, size, size*2); size *= 2; } @@ -450,14 +450,10 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, } } } - //printf("ff: ending loop, found %d\n",n); *count = n; if ( n > 0 ) { - //printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); - printf("expandingg %i->%i\n", size, n); p = itrans_peaksearch_stat_matrix_expand(p, size, n); - //printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); } return p; diff --git a/src/reproject.c b/src/reproject.c index e1ce13e..c41aa3f 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -23,24 +23,24 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { int i; - + for ( i=0; in_features; i++ ) { - + //if ( !reflection_is_easy(rflist->features[i].reflection) ) continue; - + double d = 0.0; ImageFeature *partner; int idx; - + partner = image_feature_closest(image->features, rflist->features[i].x, rflist->features[i].y, &d, &idx); - + if ( (d <= 20.0) && partner ) { rflist->features[i].partner = partner; rflist->features[i].partner_d = d; } else { rflist->features[i].partner = NULL; } - + } } @@ -55,12 +55,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */ double ux, uy, uz, uxt, uyt, uzt; /* "up" vector (and calculation intermediates) */ //int done_debug = 0; - + flist = image_feature_list_new(); - + tilt = image->tilt; omega = image->omega; - + /* 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; @@ -71,7 +71,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * ky = ny / image->lambda; kz = nz / image->lambda; /* This is the centre of the Ewald sphere */ //reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1); - + /* Determine where "up" is * See above. */ uxt = 0.0; uyt = 1.0; uzt = 0.0; @@ -79,19 +79,19 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * uxt = ux; uyt = uy; uzt = uz; ux = uxt*cos(-omega) + uyt*-sin(omega); uy = -uxt*sin(omega) + uyt*cos(omega); uz = uzt; //reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); - + reflection = reflectionlist->reflections; do { - + double xl, yl, zl; double a, b, c; double A1, A2, s1, s2, s; - + /* Get the coordinates of the reciprocal lattice point */ xl = reflection->x; yl = reflection->y; zl = reflection->z; - + /* Next, solve the relrod equation to calculate the excitation error */ a = 1.0; b = -2.0*(xl*nx + yl*ny + zl*nz); @@ -101,37 +101,37 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * s1 = 1.0/image->lambda - A1; s2 = 1.0/image->lambda - A2; if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2; - + /* Skip this reflection if s is large */ if ( fabs(s) <= smax ) { - + double xi, yi, zi; double gx, gy, gz; double cx, cy, cz; double theta; double x, y; double rx, ry, rz; - + /* Determine the intersection point */ xi = xl + s*nx; yi = yl + s*ny; zi = zl + s*nz; - + /* Calculate Bragg angle */ gx = xi - kx; gy = yi - ky; gz = zi - kz; /* This is the vector from the centre of the sphere to the intersection */ theta = angle_between(-kx, -ky, -kz, gx, gy, gz); - + if ( theta > 0.0 ) { - + double psi, disc; - + //reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED); //reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER); - + /* Calculate azimuth of point in image (clockwise from "up", will be changed later) */ cx = yi*nz-zi*ny; cy = nx*zi-nz*xi; cz = ny*xi-nx*yi; /* c = i x n */ psi = angle_between(cx, cy, cz, ux, uy, uz); - + /* Finally, resolve the +/- Pi ambiguity from the previous step */ rx = cy*nz-cz*ny; ry = nx*cz-nz*cx; rz = ny*cx-nx*cy; /* r = [i x n] x n */ disc = angle_between(rx, ry, rz, ux, uy, uz); @@ -154,15 +154,15 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * psi = 3*M_PI_2 + psi; /* Case 4 */ // if ( (i==20) && !done_debug ) printf("case 4\n"); } - + // if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi)); psi = M_PI_2 - psi; /* Anticlockwise from "+x" instead of clockwise from "up" */ // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(psi)); - + psi += omega; // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi)); // if ( (i==20) && !done_debug ) done_debug = 1; - + /* Calculate image coordinates from polar representation */ if ( image->fmode == FORMULATION_CLEN ) { x = image->camera_len*tan(theta)*cos(psi); @@ -178,36 +178,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n"); return NULL; } - + x += image->x_centre; y += image->y_centre; - + /* Sanity check */ if ( (x>=0) && (xwidth) && (y>=0) && (yheight) ) { - + /* Record the reflection */ image_add_feature_reflection(flist, x, y, image, reflection->intensity, reflection); /* Intensity should be multiplied by relrod spike function except * reprojected reflections aren't used quantitatively for anything */ - + //printf("Reflection at %f, %f\n", x, y); - + } /* else it's outside the picture somewhere */ - + } /* else this is the central beam so don't worry about it */ - + } - + reflection = reflection->next; - + } while ( reflection ); - + /* Partner features only if the image has a feature list. This allows the test * program to use this function to generate simulated data. */ if ( image->features != NULL ) { reproject_partner_features(flist, image); } - + return flist; } @@ -216,21 +216,21 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * void reproject_cell_to_lattice(ControlContext *ctx) { int i; - + if ( ctx->cell_lattice ) { reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell); } else { ctx->cell_lattice = reflection_list_from_cell(ctx->cell); } - + displaywindow_enable_cell_functions(ctx->dw, TRUE); - + /* Invalidate all the reprojected feature lists */ for ( i=0; iimages->n_images; i++ ) { image_feature_list_free(ctx->images->images[i].rflist); ctx->images->images[i].rflist = NULL; } - + } /* Notify that ctx->cell has changed (also need to call displaywindow_update) */ @@ -240,4 +240,3 @@ void reproject_lattice_changed(ControlContext *ctx) { displaywindow_update_imagestack(ctx->dw); } - -- cgit v1.2.3