From 636b71d41e8b2cd318e6df24159a6fb65301e5df Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 6 Sep 2011 15:39:02 +0200 Subject: Initial candidate ReAx stuff --- libcrystfel/src/reax.c | 478 ++++++++++++++++++++++++++++--------------------- 1 file changed, 278 insertions(+), 200 deletions(-) (limited to 'libcrystfel/src/reax.c') diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c index ff1ea582..49231438 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -45,6 +45,31 @@ struct dvec }; +struct reax_candidate +{ + struct dvec v; /* This is the vector for the candidate */ + double strength; /* How much intensity fell into the range */ +}; + + +struct reax_search_v +{ + unsigned int smin; + unsigned int smax; /* Search for vector in this range */ + + struct reax_candidate *cand; /* Candidate vectors go here */ + int n_cand; /* There are this many candidates */ +}; + + +struct reax_search +{ + struct reax_search_v *search; /* Search for these vectors */ + int n_search; /* There are this many vectors to find */ + double pmax; /* The maximum feature resolution */ +}; + + struct reax_private { IndexingPrivate base; @@ -65,14 +90,12 @@ struct reax_private }; -static double check_dir(struct dvec *dir, ImageFeatureList *flist, +static void fill_and_transform(struct dvec *dir, ImageFeatureList *flist, int nel, double pmax, double *fft_in, fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, const char *rg, struct detector *det) { int n, i; - double tot; for ( i=0; in_search; i++ ) { + + double tot = 0.0; + double peak = 0.0; + double peak_mod = 0.0; + int j; + + for ( j=s->search[i].smin; j<=s->search[i].smax; j++ ) { + + double re, im; + + re = fft_out[j][0]; + im = fft_out[j][1]; + peak += sqrt(re*re + im*im); + + } + + for ( j=0; j=s->search[i].smin) + && (j<=s->search[i].smax) ) { + if ( am > peak ) { + peak = am; + peak_mod = (double)j/(2.0*pmax); + } + } + + } + + /* If sufficiently strong, add to list of candidates */ + if ( peak > 2.0*(tot/nel) ) { + + size_t ns; + struct reax_candidate *cnew; + int cpos; + + cpos = s->search[i].n_cand; + + ns = (cpos+1) * sizeof(struct reax_candidate); + cnew = realloc(s->search[i].cand, ns); + if ( cnew == NULL ) { + ERROR("Failed to add candidate.\n"); + } else { + + s->search[i].cand = cnew; + s->search[i].cand[cpos].v.x = dir->x*peak_mod; + s->search[i].cand[cpos].v.y = dir->y*peak_mod; + s->search[i].cand[cpos].v.z = dir->z*peak_mod; + s->search[i].cand[cpos].v.th = dir->th; + s->search[i].cand[cpos].v.ph = dir->ph; + s->search[i].n_cand++; + + } + + } + } return tot; } +static void fine_search(struct reax_private *p, ImageFeatureList *flist, + double pmax, double *fft_in, fftw_complex *fft_out, + struct reax_search_v *sv, struct reax_candidate *c, + const char *rg, struct detector *det) +{ + double fom = 0.0; + double th, ph; + double inc; + struct dvec dir; + int i, s; + double max; + + inc = p->angular_inc / 100.0; + + for ( th=c->v.th-p->angular_inc; th<=c->v.th+p->angular_inc; th+=inc ) { + for ( ph=c->v.ph-p->angular_inc; ph<=c->v.ph+p->angular_inc; ph+=inc ) { + + double new_fom; + + dir.x = cos(ph) * sin(th); + dir.y = sin(ph) * sin(th); + dir.z = cos(th); + dir.th = th; + dir.ph = ph; + + fill_and_transform(&dir, flist, p->nel, pmax, fft_in, fft_out, + p->plan, rg, det); + + for ( i=sv->smin; i<=sv->smax; i++ ) { + + double re, im, m; + + re = fft_out[i][0]; + im = fft_out[i][1]; + m = sqrt(re*re + im*im); + if ( m > max ) { + max = m; + s = i; + } + } + + if ( new_fom > fom ) { + fom = new_fom; + c->v = dir; + } + + } + } +} + + +static double find_candidates(struct reax_private *p, + ImageFeatureList *flist, double pmax, + double *fft_in, fftw_complex *fft_out, + struct reax_search *s, + const char *rg, struct detector *det) +{ + double fom = 0.0; + int i; + double th, ph; + + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom; + + new_fom = check_dir(&p->directions[i], flist, + p->nel, pmax, fft_in, fft_out, p->plan, + s, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + + for ( i=0; in_search; i++ ) { + + struct reax_search_v *sv; + int j; + + sv = &s->search[i]; + for ( j=0; jn_cand; j++ ) { + fine_search(p, flist, pmax, fft_in, fft_out, sv, + &sv->cand[j], rg, det); + } + } + + return fom; +} + + +/* Set up search parameters to look for all three cell axes */ +static struct reax_search *search_all_axes(UnitCell *cell, double pmax) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + double mod_a, mod_b, mod_c; + double amin, amax; + double bmin, bmax; + double cmin, cmax; + unsigned int smin, smax; + const double ltol = 5.0; /* Direct space axis length tolerance in % */ + struct reax_search *s; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + mod_a = modulus(ax, ay, az); + amin = mod_a * (1.0-ltol/100.0); + amax = mod_a * (1.0+ltol/100.0); + + mod_b = modulus(bx, by, bz); + bmin = mod_b * (1.0-ltol/100.0); + bmax = mod_b * (1.0+ltol/100.0); + + mod_c = modulus(cx, cy, cz); + cmin = mod_c * (1.0-ltol/100.0); + cmax = mod_c * (1.0+ltol/100.0); + + s = malloc(3*sizeof(*s)); + s->pmax = pmax; + s->n_search = 1; + smin = 2.0*pmax * amin; smax = 2.0*pmax * amax; + s->search[0].smin = smin; s->search[0].smax = smax; + smin = 2.0*pmax * bmin; smax = 2.0*pmax * bmax; + s->search[1].smin = smin; s->search[1].smax = smax; + smin = 2.0*pmax * cmin; smax = 2.0*pmax * cmax; + s->search[2].smin = smin; s->search[2].smax = smax; + + return s; +} + + /* Refine a direct space vector. From Clegg (1984) */ static double iterate_refine_vector(double *x, double *y, double *z, ImageFeatureList *flist) @@ -239,63 +468,6 @@ static void refine_cell(struct image *image, UnitCell *cell, } -static void fine_search(struct reax_private *p, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, double th_cen, double ph_cen, - double *x, double *y, double *z) -{ - double fom = 0.0; - double th, ph; - double inc; - struct dvec dir; - int i, s; - double max, modv; - - inc = p->angular_inc / 100.0; - - for ( th=th_cen-p->angular_inc; th<=th_cen+p->angular_inc; th+=inc ) { - for ( ph=ph_cen-p->angular_inc; ph<=ph_cen+p->angular_inc; ph+=inc ) { - - double new_fom; - - dir.x = cos(ph) * sin(th); - dir.y = sin(ph) * sin(th); - dir.z = cos(th); - - new_fom = check_dir(&dir, flist, nel, pmax, fft_in, fft_out, - plan, smin, smax, NULL, NULL); - - if ( new_fom > fom ) { - fom = new_fom; - *x = dir.x; *y = dir.y; *z = dir.z; - } - - } - } - - s = -1; - max = 0.0; - for ( i=smin; i<=smax; i++ ) { - - double re, im, m; - - re = fft_out[i][0]; - im = fft_out[i][1]; - m = sqrt(re*re + im*im); - if ( m > max ) { - max = m; - s = i; - } - - } - assert(s>0); - - modv = (double)s / (2.0*pmax); - *x *= modv; *y *= modv; *z *= modv; -} - - static double get_model_phase(double x, double y, double z, ImageFeatureList *f, int nel, double pmax, double *fft_in, fftw_complex *fft_out, fftw_plan plan, @@ -309,8 +481,7 @@ static double get_model_phase(double x, double y, double z, ImageFeatureList *f, dir.x = x; dir.y = y; dir.z = z; - check_dir(&dir, f, nel, pmax, fft_in,fft_out, plan, smin, smax, - rg, det); + fill_and_transform(&dir, f, nel, pmax, fft_in, fft_out, plan, rg, det); s = -1; max = 0.0; @@ -336,7 +507,7 @@ static double get_model_phase(double x, double y, double z, ImageFeatureList *f, static void refine_rigid_group(struct image *image, UnitCell *cell, - const char *rg, int nel, double pmax, + const char *rg, double pmax, double *fft_in, fftw_complex *fft_out, fftw_plan plan, int smin, int smax, struct detector *det, struct reax_private *pr) @@ -353,8 +524,6 @@ static void refine_rigid_group(struct image *image, UnitCell *cell, signed int aix, aiy; signed int bix, biy; signed int cix, ciy; - double max; - int max_i, max_j; cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); @@ -362,12 +531,15 @@ static void refine_rigid_group(struct image *image, UnitCell *cell, mb = modulus(bx, by, bz); mc = modulus(cx, cy, cz); - pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, nel, pmax, - fft_in, fft_out, plan, smin, smax, rg, det); - phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, nel, pmax, - fft_in, fft_out, plan, smin, smax, rg, det); - phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, nel, pmax, - fft_in, fft_out, plan, smin, smax, rg, det); + pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, + pr->nel, pmax, fft_in, fft_out, plan, + smin, smax, rg, det); + phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, + pr->nel, pmax, fft_in, fft_out, plan, + smin, smax, rg, det); + phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, + pr->nel, pmax, fft_in, fft_out, plan, + smin, smax, rg, det); for ( i=0; in_panels; i++ ) { if ( det->panels[i].rigid_group == rg ) { @@ -461,7 +633,7 @@ static void refine_rigid_group(struct image *image, UnitCell *cell, static void refine_all_rigid_groups(struct image *image, UnitCell *cell, - int nel, double pmax, + double pmax, double *fft_in, fftw_complex *fft_out, fftw_plan plan, int smin, int smax, struct detector *det, @@ -471,67 +643,25 @@ static void refine_all_rigid_groups(struct image *image, UnitCell *cell, for ( i=0; idet->num_rigid_groups; i++ ) { refine_rigid_group(image, cell, image->det->rigid_groups[i], - nel, pmax, fft_in, fft_out, plan, smin, smax, + pmax, fft_in, fft_out, plan, smin, smax, det, p); } } -void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) +static double max_feature_resolution(ImageFeatureList *flist) { - int i; - struct reax_private *p; - double fom; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double mod_a, mod_b, mod_c; - double al, be, ga; - double th, ph; - double *fft_in; - fftw_complex *fft_out; - int smin, smax; - double amin, amax; - double bmin, bmax; - double cmin, cmax; double pmax; - int n; - const double ltol = 5.0; /* Direct space axis length - * tolerance in percent */ - const double angtol = deg2rad(1.5); /* Direct space angle tolerance - * in radians */ - - assert(pp->indm == INDEXING_REAX); - p = (struct reax_private *)pp; - - fft_in = fftw_malloc(p->nel*sizeof(double)); - fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - mod_a = modulus(ax, ay, az); - amin = mod_a * (1.0-ltol/100.0); - amax = mod_a * (1.0+ltol/100.0); - - mod_b = modulus(bx, by, bz); - bmin = mod_b * (1.0-ltol/100.0); - bmax = mod_b * (1.0+ltol/100.0); - - mod_c = modulus(cx, cy, cz); - cmin = mod_c * (1.0-ltol/100.0); - cmax = mod_c * (1.0+ltol/100.0); - - al = angle_between(bx, by, bz, cx, cy, cz); - be = angle_between(ax, ay, az, cx, cy, cz); - ga = angle_between(ax, ay, az, bx, by, bz); + int i, n; pmax = 0.0; - n = image_feature_count(image->features); + n = image_feature_count(flist); for ( i=0; ifeatures, i); + f = image_get_feature(flist, i); if ( f == NULL ) continue; val = modulus(f->rx, f->ry, f->rz); @@ -539,92 +669,40 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) } - /* Sanity check */ - if ( pmax < 1e4 ) return; - - /* Search for a */ - smin = 2.0*pmax * amin; - smax = 2.0*pmax * amax; - fom = 0.0; th = 0.0; ph = 0.0; - for ( i=0; in_dir; i++ ) { - - double new_fom; - - new_fom = check_dir(&p->directions[i], image->features, - p->nel, pmax, fft_in, fft_out, p->plan, - smin, smax, NULL, NULL); - if ( new_fom > fom ) { - fom = new_fom; - th = p->directions[i].th; - ph = p->directions[i].ph; - } - - } - fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, - p->plan, smin, smax, th, ph, &ax, &ay, &az); - - /* Search for b */ - smin = 2.0*pmax * bmin; - smax = 2.0*pmax * bmax; - fom = 0.0; th = 0.0; ph = 0.0; - for ( i=0; in_dir; i++ ) { - - double new_fom, ang; - - ang = angle_between(p->directions[i].x, p->directions[i].y, - p->directions[i].z, ax, ay, az); - if ( fabs(ang-ga) > angtol ) continue; - - new_fom = check_dir(&p->directions[i], image->features, - p->nel, pmax, fft_in, fft_out, p->plan, - smin, smax, NULL, NULL); - if ( new_fom > fom ) { - fom = new_fom; - th = p->directions[i].th; - ph = p->directions[i].ph; - } + return pmax; +} - } - fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, - p->plan, smin, smax, th, ph, &bx, &by, &bz); - /* Search for c */ - smin = 2.0*pmax * cmin; - smax = 2.0*pmax * cmax; - fom = 0.0; th = 0.0; ph = 0.0; - for ( i=0; in_dir; i++ ) { +void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) +{ + struct reax_private *p; + double *fft_in; + fftw_complex *fft_out; + double pmax; + struct reax_search *s; - double new_fom, ang; + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; - ang = angle_between(p->directions[i].x, p->directions[i].y, - p->directions[i].z, ax, ay, az); - if ( fabs(ang-be) > angtol ) continue; + fft_in = fftw_malloc(p->nel*sizeof(double)); + fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - ang = angle_between(p->directions[i].x, p->directions[i].y, - p->directions[i].z, bx, by, bz); - if ( fabs(ang-al) > angtol ) continue; + pmax = max_feature_resolution(image->features); + s = search_all_axes(cell, pmax); - new_fom = check_dir(&p->directions[i], image->features, - p->nel, pmax, fft_in, fft_out, p->plan, - smin, smax, NULL, NULL); - if ( new_fom > fom ) { - fom = new_fom; - th = p->directions[i].th; - ph = p->directions[i].ph; - } + /* Sanity check */ + if ( pmax < 1e4 ) return; - } - fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, - p->plan, smin, smax, th, ph, &cx, &cy, &cz); + find_candidates(p, image->features, pmax, fft_in, fft_out, s, + NULL, image->det); image->candidate_cells[0] = cell_new(); - cell_set_cartesian(image->candidate_cells[0], - ax, ay, az, bx, by, bz, cx, cy, cz); + /* FIXME: Assemble cell from candidates */ - refine_all_rigid_groups(image, image->candidate_cells[0], p->nel, pmax, - fft_in, fft_out, p->plan, smin, smax, - image->det, p); - map_all_peaks(image); +// refine_all_rigid_groups(image, image->candidate_cells[0], pmax, +// fft_in, fft_out, p->plan, smin, smax, +// image->det, p); +// map_all_peaks(image); refine_cell(image, image->candidate_cells[0], image->features); fftw_free(fft_in); -- cgit v1.2.3