diff options
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/src/cell.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/events.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/geometry.c | 465 | ||||
-rw-r--r-- | libcrystfel/src/geometry.h | 21 | ||||
-rw-r--r-- | libcrystfel/src/image.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 8 | ||||
-rw-r--r-- | libcrystfel/src/integration.c | 2 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 22 | ||||
-rw-r--r-- | libcrystfel/src/peakfinder8.c | 14 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 8 | ||||
-rw-r--r-- | libcrystfel/src/reflist.c | 92 | ||||
-rw-r--r-- | libcrystfel/src/reflist.h | 13 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 86 | ||||
-rw-r--r-- | libcrystfel/src/stream.h | 8 | ||||
-rw-r--r-- | libcrystfel/src/taketwo.c | 1116 | ||||
-rw-r--r-- | libcrystfel/src/xds.c | 273 |
16 files changed, 1400 insertions, 738 deletions
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 2a4e45be..cc18b49d 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -382,7 +382,7 @@ static int cell_invert(double ax, double ay, double az, return 1; } if ( gsl_linalg_LU_invert(m, perm, inv) ) { - ERROR("Couldn't invert matrix\n"); + ERROR("Couldn't invert cell matrix:\n"); gsl_matrix_free(m); gsl_permutation_free(perm); return 1; @@ -655,7 +655,7 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) return NULL; } if ( gsl_linalg_LU_invert(m, perm, inv) ) { - ERROR("Couldn't invert matrix\n"); + ERROR("Couldn't invert transformation matrix\n"); gsl_permutation_free(perm); return NULL; } diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c index 8e4eb861..d0e521d3 100644 --- a/libcrystfel/src/events.c +++ b/libcrystfel/src/events.c @@ -150,6 +150,8 @@ struct event *copy_event(struct event *ev) struct event *new_ev; int pi, di; + if ( ev == NULL ) return NULL; + if ( ev->dim_length == 0 && ev->path_length == 0) { new_ev = initialize_event(); @@ -252,6 +254,8 @@ void free_event(struct event *ev) { int pi; + if ( ev == NULL ) return; + if ( ev->path_length != 0 ) { for ( pi=0; pi<ev->path_length; pi++ ) { free(ev->path_entries[pi]); diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 576655de..9ae6d96f 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -238,76 +238,38 @@ static double random_partiality(signed int h, signed int k, signed int l, } -static double partiality(PartialityModel pmodel, - signed int h, signed int k, signed int l, - int serial, - double rlow, double rhigh, double pr) -{ - double D = rlow - rhigh; - - /* Convert to partiality */ - switch ( pmodel ) { - - default: - case PMODEL_UNITY: - return 1.0; - - case PMODEL_SCSPHERE: - return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D); - - case PMODEL_SCGAUSSIAN: - return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D); - - case PMODEL_RANDOM: - return random_partiality(h, k, l, serial); - - } -} - - static Reflection *check_reflection(struct image *image, Crystal *cryst, signed int h, signed int k, signed int l, double xl, double yl, double zl, Reflection *updateme) { - double tl; - double rlow, rhigh; /* "Excitation error" */ - double klow, khigh; /* Wavenumber */ Reflection *refl; - double cet, cez; /* Centre of Ewald sphere */ - double pr; - double del; + double R, top; + double kmin, kmax, k0, knom, k1, khalf; + double dcs, exerr; /* Don't predict 000 */ if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL; - pr = crystal_get_profile_radius(cryst); - del = image->div + crystal_get_mosaicity(cryst); - - /* "low" gives the largest Ewald sphere (wavelength short => k large) - * "high" gives the smallest Ewald sphere (wavelength long => k small) - */ - klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); - khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); - - /* If the point is looking "backscattery", reject it straight away */ - if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL; - - tl = sqrt(xl*xl + yl*yl); - - cet = -sin(del/2.0) * khigh; - cez = -cos(del/2.0) * khigh; - rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */ - - cet = sin(del/2.0) * klow; - cez = -cos(del/2.0) * klow; - rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */ - - /* Condition for reflection to be excited at all */ - if ( (updateme == NULL) - && (signbit(rlow) == signbit(rhigh)) - && (fabs(rlow) > pr) - && (fabs(rhigh) > pr) ) return NULL; + /* Calculate the limiting wavelengths, lambda0 and lambda1 + * = 1/k0 and 1/k1 respectively */ + R = fabs(crystal_get_profile_radius(cryst)); + top = R*R - xl*xl - yl*yl - zl*zl; + k0 = top/(2.0*(zl+R)); + khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl); + k1 = top/(2.0*(zl-R)); + + /* The reflection is excited if any of the reflection is within 2sigma + * of the nominal * wavelength of the X-ray beam + * (NB image->bw is full width) */ + kmin = 1.0/(image->lambda + image->lambda*image->bw); + knom = 1.0/image->lambda; + kmax = 1.0/(image->lambda - image->lambda*image->bw); + if ( (updateme == NULL) && ((k1>kmax) || (k0<kmin)) ) return NULL; + + /* Calculate excitation error */ + dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl); + exerr = 1.0/image->lambda - dcs; /* Loss of precision */ if ( updateme == NULL ) { refl = reflection_new(h, k, l); @@ -321,7 +283,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, if ( (image->det != NULL) && (updateme != NULL) ) { double fs, ss; - locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda, + locate_peak_on_panel(xl, yl, zl, knom, get_panel(updateme), &fs, &ss); set_detector_pos(refl, fs, ss); @@ -333,7 +295,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, double fs, ss; /* Position on detector */ signed int p; /* Panel number */ - p = locate_peak(xl, yl, zl, 1.0/image->lambda, + p = locate_peak(xl, yl, zl, knom, image->det, &fs, &ss); if ( p == -1 ) { reflection_free(refl); @@ -344,20 +306,9 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } - if ( unlikely(rlow < rhigh) ) { - ERROR("Reflection with rlow < rhigh!\n"); - ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", - h, k, l, rlow, rhigh); - ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw); - /* If we are updating, this is (kind of) OK */ - if ( updateme == NULL ) { - reflection_free(refl); - return NULL; - } - } - - set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by - * calculate_partialities() */ + set_kpred(refl, knom); + set_khalf(refl, khalf); + set_exerr(refl, exerr); set_lorentz(refl, 1.0); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); @@ -368,18 +319,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image) { - double azi; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double xl, yl, zl; signed int hs, ks, ls; - double rlow, rhigh, p; - double philow, phihigh, phi; - double khigh, klow; - double tl, cet, cez; - - get_partial(refl, &rlow, &rhigh, &p); + double tl, phi, azi; get_symmetric_indices(refl, &hs, &ks, &ls); @@ -390,26 +335,9 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image) yl = hs*asy + ks*bsy + ls*csy; zl = hs*asz + ks*bsz + ls*csz; - /* "low" gives the largest Ewald sphere (wavelength short => k large) - * "high" gives the smallest Ewald sphere (wavelength long => k small) - */ - klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); - khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); - tl = sqrt(xl*xl + yl*yl); - - cet = -sin(image->div/2.0) * klow; - cez = -cos(image->div/2.0) * klow; - philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0); - - cet = -sin(image->div/2.0) * khigh; - cez = -cos(image->div/2.0) * khigh; - phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0); - - /* Approximation: philow and phihigh are very similar */ - phi = (philow + phihigh) / 2.0; - - azi = atan2(yl, xl); + phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */ + azi = atan2(yl, xl); /* azimuth */ switch ( k ) { @@ -531,11 +459,17 @@ RefList *predict_to_res(Crystal *cryst, double max_res) } -static void set_unity_partialities(RefList *list) +static void set_unity_partialities(Crystal *cryst) { + RefList *list; Reflection *refl; RefListIterator *iter; + list = crystal_get_reflections(cryst); + if ( list == NULL ) { + ERROR("No reflections for partiality calculation!\n"); + return; + } for ( refl = first_refl(list, &iter); refl != NULL; refl = next_refl(refl, iter) ) @@ -546,25 +480,11 @@ static void set_unity_partialities(RefList *list) } -/** - * calculate_partialities: - * @cryst: A %Crystal - * @pmodel: A %PartialityModel - * - * Calculates the partialities for the reflections in @cryst, given the current - * crystal and image parameters. The crystal's image and reflection lists - * must be set. The specified %PartialityModel will be used. - * - * You must not have changed the crystal or image parameters since you last - * called predict_to_res() or update_predictions(), because this function - * relies on the limiting wavelength values calculated by those functions. - */ -void calculate_partialities(Crystal *cryst, PartialityModel pmodel) +static void set_random_partialities(Crystal *cryst) { RefList *list; Reflection *refl; RefListIterator *iter; - double pr; struct image *image; list = crystal_get_reflections(cryst); @@ -573,8 +493,166 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel) return; } - if ( pmodel == PMODEL_UNITY ) { - set_unity_partialities(list); + image = crystal_get_image(cryst); + if ( image == NULL ) { + ERROR("No image structure for partiality calculation!\n"); + return; + } + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + get_symmetric_indices(refl, &h, &k, &l); + set_partiality(refl, random_partiality(h, k, l, image->serial)); + set_lorentz(refl, 1.0); + } +} + + +static double do_integral(double q2, double zl, double R, + double lambda, double sig, char *verbose) +{ + int i; + double kmin, kmax, kstart, kfinis; + double inc; + double total = 0.0; + double k0, k1; + const int SAMPLES = 50; /* Number of samples for integration */ + const double N = 1.5; /* Pointiness of spectrum */ + FILE *fh = NULL; + + assert(R*R < q2); + assert(R > 0.0); + + /* Range over which P is different from zero */ + k0 = (R*R - q2)/(2.0*(zl+R)); + k1 = (R*R - q2)/(2.0*(zl-R)); + + /* Check for reflections partially "behind the beam" */ + if ( k0 < 0.0 ) k0 = +INFINITY; + if ( k1 < 0.0 ) k1 = +INFINITY; + + /* Range over which E is significantly different from zero */ + kmin = 1.0 / (lambda + 3.0*sig); + kmax = 1.0 / (lambda - 3.0*sig); + + /* Calculate range over which E*P is different from zero */ + if ( k0 < k1 ) { + STATUS("%e %e\n", k0, k1); + STATUS("%e %e %e\n", q2, zl, R); + STATUS("%e %e\n", kmin, kmax); + } + assert((isinf(k0) && isinf(k1)) || (k0 > k1)); + assert(kmax > kmin); + if ( kmin < k1 ) { + if ( kmax < k1 ) { + /* Case 1 */ + kstart = k1; kfinis = k0; + return 0.0; + } else if ( kmax < k0 ) { + /* Case 2 (kmax > k1)*/ + kstart = k1; kfinis = kmax; + } else { + /* Case 3 (kmax > k1 and kmax > k0) */ + kstart = k1; kfinis = k0; + } + } else if ( kmin < k0 ) { /* kmin > k1 */ + if ( kmax < k0 ) { + /* Case 4 */ + kstart = kmin; kfinis = kmax; + } else { + /* Case 5 (kmax > k0) */ + kstart = kmin; kfinis = k0; + } + } else { + /* Case 6 (kmin > k1 and (kmax>)kmin > k0) */ + kstart = k1; kfinis = k0; + return 0.0; + } + + if ( kstart < 0.0 ) kstart = kmin; + if ( kfinis < 0.0 ) kfinis = kmax; + + inc = (kfinis - kstart) / SAMPLES; + + if ( verbose ) { + char fn[64]; + snprintf(fn, 63, "partial%s.graph", verbose); + fh = fopen(fn, "w"); + fprintf(fh, " n p wavelength E P\n"); + STATUS("Nominal k = %e m^-1\n", 1.0/lambda); + STATUS(" (wavelength %e m)\n", lambda); + STATUS("Bandwidth %e m\n", sig); + STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl)); + STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl))); + STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0); + STATUS(" (wavelengths from %e to %e m\n", 1.0/k1, 1.0/k0); + STATUS("Beam goes from %e to %e m^-1\n", kmin, kmax); + STATUS(" (wavelengths from %e to %e m\n", 1.0/kmin, 1.0/kmax); + STATUS("Integration goes from %e to %e m^-1\n", kstart, kfinis); + STATUS(" (wavelengths from %e to %e m\n", 1.0/kstart, 1.0/kfinis); + } + + for ( i=0; i<SAMPLES; i++ ) { + + double p, kp, lrel; + double E, P; + + kp = kstart + i*inc; + double pref = sqrt(q2 + kp*kp + 2.0*zl*kp)/(2.0*R); + p = pref + 0.5 - kp/(2.0*R); + + /* Spectral energy term */ + lrel = fabs(1.0/kp - lambda); + E = exp(-0.5 * pow(lrel / sig, N)); + E /= sqrt(2.0 * M_PI * sig); + + /* RLP profile term */ + P = 4.0*p * (1.0 - p); + + total += E*P*inc; + + if ( fh != NULL ) { + fprintf(fh, "%3i %f %e %e %e\n", i, p, 1.0/kp, E, P); + } + } + + if ( fh != NULL ) fclose(fh); + + if ( isnan(total) ) { + ERROR("NaN total!\n"); + STATUS("Nominal k = %e m^-1\n", 1.0/lambda); + STATUS(" (wavelength %e m)\n", lambda); + STATUS("Bandwidth %e m\n", sig); + STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl)); + STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl))); + STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0); + STATUS(" (wavelengths from %e to %e m\n", 1.0/k1, 1.0/k0); + STATUS("Beam goes from %e to %e m^-1\n", kmin, kmax); + STATUS(" (wavelengths from %e to %e m\n", 1.0/kmin, 1.0/kmax); + STATUS("Integration goes from %e to %e m^-1\n", kstart, kfinis); + STATUS(" (wavelengths from %e to %e m\n", 1.0/kstart, 1.0/kfinis); + } + + return total; +} + + +static void ginn_spectrum_partialities(Crystal *cryst) +{ + RefList *list; + Reflection *refl; + RefListIterator *iter; + double r0, m, lambda, sig; + struct image *image; + UnitCell *cell; + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + + list = crystal_get_reflections(cryst); + if ( list == NULL ) { + ERROR("No reflections for partiality calculation!\n"); return; } @@ -584,31 +662,88 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel) return; } - pr = crystal_get_profile_radius(cryst); + cell = crystal_get_cell(cryst); + if ( cell == NULL ) { + ERROR("No unit cell for partiality calculation!\n"); + return; + } + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + r0 = fabs(crystal_get_profile_radius(cryst)); + m = crystal_get_mosaicity(cryst); + lambda = image->lambda; + sig = image->bw * lambda / 2.0; for ( refl = first_refl(crystal_get_reflections(cryst), &iter); refl != NULL; refl = next_refl(refl, iter) ) { - double rlow, rhigh, part; + double R; signed int h, k, l; + double xl, yl, zl; + double q2; + double total, norm; get_symmetric_indices(refl, &h, &k, &l); - get_partial(refl, &rlow, &rhigh, &part); + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; + + /* Radius of rlp profile */ + q2 = xl*xl + yl*yl + zl*zl; + + R = r0 + m * sqrt(q2); + + //char tmp[256]; + //snprintf(tmp, 255, "-%i,%i,%i", h, k, l); + char *tmp = NULL; - part = partiality(pmodel, h, k, l, image->serial, - rlow, rhigh, pr); + total = do_integral(q2, zl, R, lambda, sig, tmp); + norm = do_integral(q2, -0.5*q2*lambda, R, lambda, sig, NULL); - if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) { - ERROR("Assigning NAN partiality!\n"); - ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", - h, k, l, rlow, rhigh); - ERROR("R = %e, bw = %e\n", pr, image->bw); - ERROR("D = %e\n", rlow - rhigh); - abort(); - } + set_partiality(refl, total/norm); + set_lorentz(refl, 1.0); - set_partiality(refl, part); + assert(total <= 2.0*norm); + + } +} + + +/** + * calculate_partialities: + * @cryst: A %Crystal + * @pmodel: A %PartialityModel + * + * Calculates the partialities for the reflections in @cryst, given the current + * crystal and image parameters. The crystal's image and reflection lists + * must be set. The specified %PartialityModel will be used. + * + * You must not have changed the crystal or image parameters since you last + * called predict_to_res() or update_predictions(), because this function + * relies on the limiting wavelength values calculated by those functions. + */ +void calculate_partialities(Crystal *cryst, PartialityModel pmodel) +{ + switch ( pmodel ) { + + case PMODEL_UNITY : + set_unity_partialities(cryst); + break; + + case PMODEL_XSPHERE : + ginn_spectrum_partialities(cryst); + break; + + case PMODEL_RANDOM : + set_random_partialities(cryst); + break; + + default : + ERROR("Unknown partiality model %i\n", pmodel); + break; } } @@ -694,29 +829,30 @@ void polarisation_correction(RefList *list, UnitCell *cell, struct image *image) /* Returns dx_h/dP, where P = any parameter */ -double x_gradient(int param, Reflection *refl, UnitCell *cell, - struct panel *p, double lambda) +double x_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p) { signed int h, k, l; - double x, z, wn; - double ax, ay, az, bx, by, bz, cx, cy, cz; + double xl, zl, kpred; + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; get_indices(refl, &h, &k, &l); - wn = 1.0 / lambda; - cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - x = h*ax + k*bx + l*cx; - z = h*az + k*bz + l*cz; + kpred = get_kpred(refl); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + xl = h*asx + k*bsx + l*csx; + zl = h*asz + k*bsz + l*csz; switch ( param ) { case GPARAM_ASX : - return h * p->clen / (wn+z); + return h * p->clen / (kpred + zl); case GPARAM_BSX : - return k * p->clen / (wn+z); + return k * p->clen / (kpred + zl); case GPARAM_CSX : - return l * p->clen / (wn+z); + return l * p->clen / (kpred + zl); case GPARAM_ASY : return 0.0; @@ -728,13 +864,13 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell, return 0.0; case GPARAM_ASZ : - return -h * x * p->clen / (wn*wn + 2*wn*z + z*z); + return -h * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl); case GPARAM_BSZ : - return -k * x * p->clen / (wn*wn + 2*wn*z + z*z); + return -k * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl); case GPARAM_CSZ : - return -l * x * p->clen / (wn*wn + 2*wn*z + z*z); + return -l * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl); case GPARAM_DETX : return -1; @@ -743,7 +879,7 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell, return 0; case GPARAM_CLEN : - return x / (wn+z); + return xl / (kpred+zl); } @@ -753,18 +889,19 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell, /* Returns dy_h/dP, where P = any parameter */ -double y_gradient(int param, Reflection *refl, UnitCell *cell, - struct panel *p, double lambda) +double y_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p) { signed int h, k, l; - double y, z, wn; - double ax, ay, az, bx, by, bz, cx, cy, cz; + double yl, zl, kpred; + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; get_indices(refl, &h, &k, &l); - wn = 1.0 / lambda; - cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - y = h*ay + k*by + l*cy; - z = h*az + k*bz + l*cz; + kpred = get_kpred(refl); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; switch ( param ) { @@ -778,22 +915,22 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell, return 0.0; case GPARAM_ASY : - return h * p->clen / (wn+z); + return h * p->clen / (kpred + zl); case GPARAM_BSY : - return k * p->clen / (wn+z); + return k * p->clen / (kpred + zl); case GPARAM_CSY : - return l * p->clen / (wn+z); + return l * p->clen / (kpred + zl); case GPARAM_ASZ : - return -h * y * p->clen / (wn*wn + 2*wn*z + z*z); + return -h * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl); case GPARAM_BSZ : - return -k * y * p->clen / (wn*wn + 2*wn*z + z*z); + return -k * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl); case GPARAM_CSZ : - return -l * y * p->clen / (wn*wn + 2*wn*z + z*z); + return -l * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl); case GPARAM_DETX : return 0; @@ -802,7 +939,7 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell, return -1; case GPARAM_CLEN : - return y / (wn+z); + return yl / (kpred+zl); } diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 4fbd0bca..529f112e 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -3,7 +3,7 @@ * * Geometry of diffraction * - * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * @@ -47,8 +47,7 @@ extern "C" { /** * PartialityModel: * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. - * @PMODEL_SCSPHERE : Sphere model with source coverage factor included - * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included + * @PMODEL_XSPHERE : Flat sphere model with super-Gaussian spectrum * @PMODEL_RANDOM : Randomly assigned partialities * * A %PartialityModel describes a geometrical model which can be used to @@ -57,8 +56,7 @@ extern "C" { typedef enum { PMODEL_UNITY, - PMODEL_SCSPHERE, - PMODEL_SCGAUSSIAN, + PMODEL_XSPHERE, PMODEL_RANDOM, } PartialityModel; @@ -80,8 +78,13 @@ enum gparam { GPARAM_DETX, GPARAM_DETY, GPARAM_CLEN, - GPARAM_OSF, - GPARAM_BFAC + GPARAM_OSF, /* Linear scale factor */ + GPARAM_BFAC, /* D-W scale factor */ + GPARAM_ANG1, /* Out of plane rotation angles of crystal */ + GPARAM_ANG2, + GPARAM_WAVELENGTH, + + GPARAM_EOL /* End of list */ }; @@ -99,9 +102,9 @@ extern double sphere_fraction(double rlow, double rhigh, double pr); extern double gaussian_fraction(double rlow, double rhigh, double pr); extern double x_gradient(int param, Reflection *refl, UnitCell *cell, - struct panel *p, double lambda); + struct panel *p); extern double y_gradient(int param, Reflection *refl, UnitCell *cell, - struct panel *p, double lambda); + struct panel *p); #ifdef __cplusplus } diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 99c3d57e..9769c66a 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -216,7 +216,7 @@ struct image { /* Per-shot radiation values */ double lambda; /* Wavelength in m */ double div; /* Divergence in radians */ - double bw; /* Bandwidth as a fraction */ + double bw; /* FWHM bandwidth as a fraction */ /* Detected peaks */ long long num_peaks; diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 20595604..b20c689a 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -316,9 +316,9 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, ERROR("To disable prediction refinement ('norefine'), use --no-refine.\n"); ERROR("To check cell axes only ('axes'), use --no-cell-combinations.\n"); ERROR("To disable all unit cell checks ('raw'), use --no-check-cell.\n"); + ERROR("To disable peak alignment check ('bad'), use --no-check-peaks.\n"); ERROR("To disable indexing retry ('noretry'), use --no-retry.\n"); - ERROR("Multi-lattice indexing ('multi') is now the default: " - "use --no-multi to disable it.\n"); + ERROR("To enable multi-lattice indexing by 'delete and retry', use --multi\n"); ERROR("------------------\n"); free(methods); return NULL; @@ -1044,8 +1044,8 @@ char *detect_indexing_methods(UnitCell *cell) do_probe(dirax_probe, cell, methods); do_probe(asdf_probe, cell, methods); do_probe(xds_probe, cell, methods); - do_probe(taketwo_probe, cell, methods); - /* Don't automatically use Felix (yet) */ + /* Don't automatically use TakeTwo or Felix (yet) */ + //do_probe(taketwo_probe, cell, methods); //do_probe(felix_probe, cell, methods); if ( strlen(methods) == 0 ) { diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index d7a34954..880ae87e 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -1760,7 +1760,7 @@ void integrate_all_2(struct image *image, IntegrationMethod meth, IntDiag int_diag, signed int idh, signed int idk, signed int idl) { - integrate_all_3(image, meth, PMODEL_SCSPHERE, 0.0, + integrate_all_3(image, meth, PMODEL_XSPHERE, 0.0, ir_inn, ir_mid, ir_out, int_diag, idh, idk, idl); } diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index e59cf4e2..d1b7ffed 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -500,12 +500,18 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) switch ( mosflm->step ) { case 1 : + /* Backwards-compatible workaround for different Mosflm behaviour + * in version 7.2.2 */ + mosflm_sendline("DETECTOR ADSC\n", mosflm); + break; + + case 2 : mosflm_sendline("DETECTOR ROTATION HORIZONTAL" " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" " RECTANGULAR\n", mosflm); break; - case 2 : + case 3 : if ( (mosflm->mp->indm & INDEXING_USE_LATTICE_TYPE) && (mosflm->mp->template != NULL) ) { @@ -527,32 +533,32 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) } break; - case 3 : + case 4 : snprintf(tmp, 255, "DISTANCE %f\n", FAKE_CLEN*1e3); mosflm_sendline(tmp, mosflm); break; - case 4 : + case 5 : mosflm_sendline("BEAM 0.0 0.0\n", mosflm); break; - case 5 : + case 6 : wavelength = image->lambda*1e10; snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength); mosflm_sendline(tmp, mosflm); break; - case 6 : + case 7 : snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile); mosflm_sendline(tmp, mosflm); break; - case 7 : + case 8 : snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile); mosflm_sendline(tmp, mosflm); break; - case 8 : + case 9 : if ( mosflm->mp->indm & INDEXING_USE_CELL_PARAMETERS ) { cell_get_parameters(mosflm->mp->template, &a, &b, &c, &alpha, &beta, &gamma); @@ -572,7 +578,7 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) mosflm_sendline(tmp, mosflm); break; - case 9 : + case 10 : mosflm_sendline("GO\n", mosflm); mosflm->finished_ok = 1; break; diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 417465df..f96750de 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -32,6 +32,7 @@ #include <config.h> #endif +#include <float.h> #include <math.h> #include <stdlib.h> @@ -363,13 +364,15 @@ static void fill_radial_bins(float *data, int curr_r; float value; - for ( iss = 0; iss<h ; iss++ ) { - for ( ifs = 0; ifs<w ; ifs++ ) { + for ( iss=0; iss<h; iss++ ) { + for ( ifs=0; ifs<w; ifs++ ) { pidx = iss * w + ifs; if ( mask[pidx] != 0 ) { curr_r = (int)rint(r_map[pidx]); value = data[pidx]; - if ( value < rthreshold[curr_r ] && value>lthreshold[curr_r]) { + if ( value < rthreshold[curr_r] + && value > lthreshold[curr_r] ) + { roffset[curr_r] += value; rsigma[curr_r] += (value * value); rcount[curr_r] += 1; @@ -397,8 +400,8 @@ static void compute_radial_stats(float *rthreshold, if ( rcount[ri] == 0 ) { roffset[ri] = 0; rsigma[ri] = 0; - rthreshold[ri] = 1e9; - lthreshold[ri] = -1e9; + rthreshold[ri] = FLT_MAX; + lthreshold[ri] = FLT_MIN; } else { this_offset = roffset[ri] / rcount[ri]; this_sigma = rsigma[ri] / rcount[ri] - (this_offset * this_offset); @@ -1064,6 +1067,7 @@ int peakfinder8(struct image *img, int max_n_peaks, for ( i=0 ; i<rstats->n_rad_bins ; i++) { rstats->rthreshold[i] = 1e9; + rstats->lthreshold[i] = -1e9; } for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) { diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 3246dbec..14a60bc7 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -91,9 +91,7 @@ static void twod_mapping(double fs, double ss, double *px, double *py, static double r_dev(struct reflpeak *rp) { /* Excitation error term */ - double rlow, rhigh, p; - get_partial(rp->refl, &rlow, &rhigh, &p); - return (rlow+rhigh)/2.0; + return get_exerr(rp->refl); } @@ -425,7 +423,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, /* Positional x terms */ for ( k=0; k<num_params; k++ ) { gradients[k] = x_gradient(rv[k], rps[i].refl, cell, - rps[i].panel, image->lambda); + rps[i].panel); } for ( k=0; k<num_params; k++ ) { @@ -457,7 +455,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, /* Positional y terms */ for ( k=0; k<num_params; k++ ) { gradients[k] = y_gradient(rv[k], rps[i].refl, cell, - rps[i].panel, image->lambda); + rps[i].panel); } for ( k=0; k<num_params; k++ ) { diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index 707c802e..a35aa575 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -67,10 +67,11 @@ struct _refldata { signed int ls; /* Partiality and related geometrical stuff */ - double rlow; /* Low excitation error */ - double rhigh; /* High excitation error */ - double p; /* Partiality */ - double L; /* Lorentz factor */ + double khalf; /* Wavenumber of middle of reflection */ + double kpred; /* Wavenumber for prediction */ + double exerr; /* Excitation error */ + double p; /* Partiality */ + double L; /* Lorentz factor */ /* Location in image */ double fs; @@ -422,22 +423,43 @@ double get_intensity(const Reflection *refl) /** - * get_partial: + * get_khalf * @refl: A %Reflection - * @rlow: Location at which to store the "low" excitation error - * @rhigh: Location at which to store the "high" excitation error - * @p: Location at which to store the partiality * - * This function is used during post refinement (in conjunction with - * set_partial()) to get access to the details of the partiality calculation. + * Returns: the wavenumber at the centre of the reflection * **/ -void get_partial(const Reflection *refl, double *rlow, double *rhigh, - double *p) +double get_khalf(const Reflection *refl) { - *rlow = refl->data.rlow; - *rhigh = refl->data.rhigh; - *p = get_partiality(refl); + return refl->data.khalf; +} + + + + +/** + * get_kpred: + * @refl: A %Reflection + * + * Returns: the wavenumber which should be used for prediction of this reflection + * + **/ +double get_kpred(const Reflection *refl) +{ + return refl->data.kpred; +} + + +/** + * get_exerr: + * @refl: A %Reflection + * + * Returns: the excitation error (in m^-1) for this reflection + * + **/ +double get_exerr(const Reflection *refl) +{ + return refl->data.exerr; } @@ -612,23 +634,43 @@ void set_panel(Reflection *refl, struct panel *p) refl->data.panel = p; } +/** + * set_khalf: + * @refl: A %Reflection + * @khalf: The wavenumber at which the reflection should be predicted + * + * Sets the wavenumber at the centre of the reflection. + **/ +void set_khalf(Reflection *refl, double khalf) +{ + refl->data.khalf = khalf; +} + + /** - * set_partial: + * set_kpred: * @refl: A %Reflection - * @rlow: The "low" excitation error - * @rhigh: The "high" excitation error - * @p: The partiality + * @kpred: The wavenumber at which the reflection should be predicted * - * This function is used during post refinement (in conjunction with - * get_partial()) to get access to the details of the partiality calculation. + * Sets the wavenumber at which the reflection should be predicted. + * Used by predict_to_res() and update_predictions() + **/ +void set_kpred(Reflection *refl, double kpred) +{ + refl->data.kpred = kpred; +} + + +/** + * set_exerr: + * @refl: A %Reflection + * @exerr: The excitation error for the reflection * **/ -void set_partial(Reflection *refl, double rlow, double rhigh, double p) +void set_exerr(Reflection *refl, double exerr) { - refl->data.rlow = rlow; - refl->data.rhigh = rhigh; - refl->data.p = p; + refl->data.exerr = exerr; } diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index f2126ac9..9f494474 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2016 Thomas White <taw@physics.org> + * 2011-2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -87,6 +87,9 @@ extern Reflection *next_found_refl(Reflection *refl); extern void get_detector_pos(const Reflection *refl, double *fs, double *ss); extern struct panel *get_panel(const Reflection *refl); extern double get_partiality(const Reflection *refl); +extern double get_khalf(const Reflection *refl); +extern double get_kpred(const Reflection *refl); +extern double get_exerr(const Reflection *refl); extern double get_lorentz(const Reflection *refl); extern void get_indices(const Reflection *refl, signed int *h, signed int *k, signed int *l); @@ -94,8 +97,6 @@ extern void get_symmetric_indices(const Reflection *refl, signed int *hs, signed int *ks, signed int *ls); extern double get_intensity(const Reflection *refl); -extern void get_partial(const Reflection *refl, double *rlow, double *rhigh, - double *p); extern int get_redundancy(const Reflection *refl); extern double get_temp1(const Reflection *refl); extern double get_temp2(const Reflection *refl); @@ -109,7 +110,9 @@ extern int get_flag(const Reflection *refl); extern void copy_data(Reflection *to, const Reflection *from); extern void set_detector_pos(Reflection *refl, double fs, double ss); extern void set_panel(Reflection *refl, struct panel *p); -extern void set_partial(Reflection *refl, double rlow, double rhigh, double p); +extern void set_kpred(Reflection *refl, double kpred); +extern void set_khalf(Reflection *refl, double khalf); +extern void set_exerr(Reflection *refl, double exerr); extern void set_partiality(Reflection *refl, double p); extern void set_lorentz(Reflection *refl, double L); extern void set_intensity(Reflection *refl, double intensity); diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 0f4dcd6b..9d97ff98 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -64,6 +64,7 @@ struct _stream int major_version; int minor_version; + char *audit_info; long long int ln; @@ -1334,6 +1335,64 @@ void write_stream_header(FILE *ofh, int argc, char *argv[]) } +char *stream_audit_info(Stream *st) +{ + if ( st->audit_info == NULL ) return NULL; + return strdup(st->audit_info); +} + + +static void read_audit_lines(Stream *st) +{ + int done = 0; + size_t len = 0; + int first = 1; + + st->audit_info = malloc(4096); + if ( st->audit_info == NULL ) { + ERROR("Failed to allocate memory for audit information\n"); + return; + } + + /* Read lines from stream until one of them starts with "-----", + * then rewind to the start of that line */ + do { + + char line[1024]; + char *rval; + long pos; + + pos = ftell(st->fh); + + rval = fgets(line, 1023, st->fh); + if ( rval == NULL ) { + ERROR("Failed to read stream audit info.\n"); + close_stream(st); + return; + } + + if ( strncmp(line, "-----", 5) == 0 ) { + fseek(st->fh, pos, SEEK_SET); + done = 1; + } else { + chomp(line); + len += strlen(line); + if ( len > 4090 ) { + ERROR("Too much audit information.\n"); + return; + } else { + if ( !first ) { + strcat(st->audit_info, "\n"); + } + first = 0; + strcat(st->audit_info, line); + } + } + + } while ( !done ); +} + + Stream *open_stream_for_read(const char *filename) { Stream *st; @@ -1341,6 +1400,7 @@ Stream *open_stream_for_read(const char *filename) st = malloc(sizeof(struct _stream)); if ( st == NULL ) return NULL; st->old_indexers = 0; + st->audit_info = NULL; if ( strcmp(filename, "-") == 0 ) { st->fh = stdin; @@ -1383,6 +1443,8 @@ Stream *open_stream_for_read(const char *filename) st->ln = 1; + read_audit_lines(st); + return st; } @@ -1408,6 +1470,7 @@ Stream *open_stream_fd_for_write(int fd) st = malloc(sizeof(struct _stream)); if ( st == NULL ) return NULL; st->old_indexers = 0; + st->audit_info = NULL; st->fh = fdopen(fd, "w"); if ( st->fh == NULL ) { @@ -1435,12 +1498,13 @@ static void write_cell_to_stream(Stream *st, UnitCell *cell) /** - * open_stream_for_write_3 + * open_stream_for_write_4 * @filename: Filename of new stream * @geom_filename: The geometry filename to copy * @cell: A %UnitCell to write into the stream * @argc: The number of arguments to the program * @argv: The arguments to the program + * @indm_str: The list of indexing methods * * Creates a new stream with name @filename, and adds the stream format * and version header, plus a verbatim copy of the geometry file and the unit @@ -1448,9 +1512,9 @@ static void write_cell_to_stream(Stream *st, UnitCell *cell) * * Returns: a %Stream, or NULL on failure. */ -Stream *open_stream_for_write_3(const char *filename, +Stream *open_stream_for_write_4(const char *filename, const char *geom_filename, UnitCell *cell, - int argc, char *argv[]) + int argc, char *argv[], const char *indm_str) { Stream *st; @@ -1458,6 +1522,7 @@ Stream *open_stream_for_write_3(const char *filename, st = malloc(sizeof(struct _stream)); if ( st == NULL ) return NULL; st->old_indexers = 0; + st->audit_info = NULL; st->fh = fopen(filename, "w"); if ( st->fh == NULL ) { @@ -1477,6 +1542,10 @@ Stream *open_stream_for_write_3(const char *filename, if ( (argc > 0) && (argv != NULL) ) { write_command(st, argc, argv); } + + if ( indm_str != NULL ) { + fprintf(st->fh, "Indexing methods selected: %s\n", indm_str); + } if ( geom_filename != NULL ) { write_geometry_file(st, geom_filename); } @@ -1488,6 +1557,15 @@ Stream *open_stream_for_write_3(const char *filename, } +Stream *open_stream_for_write_3(const char *filename, + const char *geom_filename, UnitCell *cell, + int argc, char *argv[]) +{ + return open_stream_for_write_4(filename, geom_filename, cell, + argc, argv, NULL); +} + + /** * open_stream_for_write_2 * @filename: Filename of new stream @@ -1510,6 +1588,7 @@ Stream *open_stream_for_write_2(const char *filename, st = malloc(sizeof(struct _stream)); if ( st == NULL ) return NULL; st->old_indexers = 0; + st->audit_info = NULL; st->fh = fopen(filename, "w"); if ( st->fh == NULL ) { @@ -1576,6 +1655,7 @@ int get_stream_fd(Stream *st) void close_stream(Stream *st) { + free(st->audit_info); fclose(st->fh); free(st); } diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index 2ed7b050..a8e3e2ee 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -3,11 +3,11 @@ * * Stream tools * - * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2017 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * 2014 Valerio Mariani * 2011 Andrew Aquila * @@ -99,6 +99,9 @@ extern Stream *open_stream_for_write_2(const char *filename, extern Stream *open_stream_for_write_3(const char *filename, const char* geom_filename, UnitCell *cell, int argc, char *argv[]); +extern Stream *open_stream_for_write_4(const char *filename, + const char *geom_filename, UnitCell *cell, + int argc, char *argv[], const char *indm_str); extern Stream *open_stream_fd_for_write(int fd); extern int get_stream_fd(Stream *st); extern void close_stream(Stream *st); @@ -121,6 +124,7 @@ extern void write_command(Stream *st, int argc, char *argv[]); extern void write_geometry_file(Stream *st, const char *geom_filename); extern int rewind_stream(Stream *st); extern int is_stream(const char *filename); +extern char *stream_audit_info(Stream *st); #ifdef __cplusplus } diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 0857f813..234b00c3 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -28,6 +28,69 @@ * */ +/** + * \class TakeTwo + * Code outline. + * --- Get ready for calculation --- + * Pre-calculate symmetry operations (generate_rotation_symops()) + * Pre-calculate theoretical vectors from unit cell dimensions + * (gen_theoretical_vecs()) + * Generate observed vectors from data (gen_observed_vecs()) + * Match observed vectors to theoretical vectors (match_obs_to_cell_vecs()) + * + * --- Business bit --- + * ... n.b. rearranging to find all seeds in advance. + * + * Find starting seeds (find_seeds()): + * - Loop through pairs of observed vectors + * - If they share a spot, find matching pairs of theoretical vectors + * - Remove all duplicate matches due to symmetry operations + * - For the remainder, loop through the matches and extend the seeds + * (start_seed()). + * - If it returns a membership greater than the highest member threshold, + * return the matrix to CrystFEL. + * + * Extending a seed (start_seed()): + * - Generate a rotation matrix which matches the chosen start seed. + * - Loop through all observed vectors starting from 0. + * - Find another vector to add to the network, if available + * (find_next_index()). + * - If the index is not available, then give up if there were too many dead + * ends. Otherwise, remove the last member and pretend like that didn't + * happen, so it won't happen again. + * - Add the vector to increment the membership list. + * - If the membership number exceeds the maximum, tidy up the solution and + * return a success. + * - If the membership does not, then resume the loop and search for the + * next vector. + * + * Finding the next member (find_next_index()): + * - Go through the observed vectors, starting from the last index + 1 to + * explore only the "new" vectors. + * - If the vector does not share a spot with the current array of vectors, + * then skip it. + * - We must loop through all the current vectors in the network, and see if + * they match the newcomer for a given matching theoretical vector. + * - We only accept a match if the rotation matrix matches the seed matrix + * for a single matching theoretical vector. + * - If it does match, we can return a success. + * + * Tidying the solution (finish_solution()): + * - This chooses the most common rotation matrix of the bunch to choose to + * send to CrystFEL. But this should probably take the average instead, + * which is very possible. + * + * Clean up the mess (cleanup_taketwo_obs_vecs()) + */ + +/** + * Helen's to-do list + * - + * + * + * - Improve the final solution + */ + #include <gsl/gsl_matrix.h> #include <gsl/gsl_blas.h> #include <float.h> @@ -57,40 +120,90 @@ struct SpotVec { struct rvec obsvec; - struct rvec *matches; + struct TheoryVec *matches; int match_num; - struct rvec *asym_matches; - int asym_match_num; + int assignment; + int in_network; double distance; struct rvec *her_rlp; struct rvec *his_rlp; }; +/** + * theoryvec + */ + +struct TheoryVec +{ + struct rvec vec; + int asym; +}; + + +/** + * seed + */ + +struct Seed +{ + int obs1; + int obs2; + int idx1; + int idx2; + double score; +}; struct taketwo_private { IndexingMethod indm; UnitCell *cell; + int serial_num; /**< Serial of last image, -1 if unassigned */ + unsigned int xtal_num; /**< last number of crystals recorded */ + + struct TheoryVec *theory_vecs; /**< Theoretical vectors for given unit cell */ + unsigned int vec_count; /**< Number of theoretical vectors */ + + gsl_matrix **prevSols; /**< Previous solutions to be ignored */ + unsigned int numPrevs; /**< Previous solution count */ + double *prevScores; /**< previous solution scores */ + unsigned int *membership; /**< previous solution was success or failure */ + }; -// These rotation symmetry operators +/** + * Internal structure which gets passed the various functions looking after + * the indexing bits and bobs. */ struct TakeTwoCell { - UnitCell *cell; + UnitCell *cell; /**< Contains unit cell dimensions */ gsl_matrix **rotSymOps; unsigned int numOps; - struct SpotVec **obs_vecs; // Pointer to an array + + struct SpotVec *obs_vecs; + struct Seed *seeds; + int seed_count; int obs_vec_count; /* Options */ int member_thresh; - double len_tol; /* In reciprocal metres */ - double angle_tol; /* In radians */ - double trace_tol; /* Contains sqrt(4*(1-cos(angle))) */ + double len_tol; /**< In reciprocal metres */ + double angle_tol; /**< In radians */ + double trace_tol; /**< Contains sqrt(4*(1-cos(angle))) */ + + /** A given solution to refine */ + gsl_matrix *solution; + double x_ang; /**< Rotations in radians to apply to x axis of solution */ + double y_ang; /**< Rotations in radians to apply to y axis of solution */ + double z_ang; /**< Rotations in radians to apply to z axis of solution */ + + /**< Temporary memory always allocated for calculations */ gsl_matrix *twiz1Tmp; + /**< Temporary memory always allocated for calculations */ gsl_matrix *twiz2Tmp; + /**< Temporary memory always allocated for calculations */ gsl_vector *vec1Tmp; + /**< Temporary memory always allocated for calculations */ gsl_vector *vec2Tmp; }; @@ -104,6 +217,9 @@ struct TakeTwoCell /* Threshold for network members to consider a potential solution */ #define NETWORK_MEMBER_THRESHOLD (20) +/* Minimum for network members to consider a potential solution */ +#define MINIMUM_MEMBER_THRESHOLD (5) + /* Maximum dead ends for a single branch extension during indexing */ #define MAX_DEAD_ENDS (10) @@ -117,6 +233,12 @@ struct TakeTwoCell /* Tolerance for rot_mats_are_similar */ #define TRACE_TOLERANCE (deg2rad(3.0)) +/* Initial step size for refinement of solutions */ +#define ANGLE_STEP_SIZE (deg2rad(0.5)) + +/* Final required step size for refinement of solutions */ +#define ANGLE_CONVERGE_SIZE (deg2rad(0.01)) + /* TODO: Multiple lattices */ @@ -145,6 +267,15 @@ static struct rvec new_rvec(double new_u, double new_v, double new_w) return new_rvector; } +static struct rvec rvec_add_rvec(struct rvec first, struct rvec second) +{ + struct rvec diff = new_rvec(second.u + first.u, + second.v + first.v, + second.w + first.w); + + return diff; +} + static struct rvec diff_vec(struct rvec from, struct rvec to) { @@ -240,6 +371,41 @@ static void rotation_around_axis(struct rvec c, double th, gsl_matrix_set(res, 2, 2, cos(th) + c.w*c.w*omc); } +/** Rotate GSL matrix by three angles along x, y and z axes */ +static void rotate_gsl_by_angles(gsl_matrix *sol, double x, double y, + double z, gsl_matrix *result) +{ + gsl_matrix *x_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *y_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *z_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *xy_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *xyz_rot = gsl_matrix_alloc(3, 3); + + struct rvec x_axis = new_rvec(1, 0, 0); + struct rvec y_axis = new_rvec(1, 0, 0); + struct rvec z_axis = new_rvec(1, 0, 0); + + rotation_around_axis(x_axis, x, x_rot); + rotation_around_axis(y_axis, y, y_rot); + rotation_around_axis(z_axis, z, z_rot); + + /* Collapse the rotations in x and y onto z */ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, x_rot, + y_rot, 0.0, xy_rot); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xy_rot, + z_rot, 0.0, xyz_rot); + + /* Apply the whole rotation offset to the solution */ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xyz_rot, + sol, 0.0, result); + + gsl_matrix_free(x_rot); + gsl_matrix_free(y_rot); + gsl_matrix_free(z_rot); + gsl_matrix_free(xy_rot); + gsl_matrix_free(xyz_rot); +} + /* Rotate vector (vec1) around axis (axis) by angle theta. Find value of * theta for which the angle between (vec1) and (vec2) is minimised. */ @@ -297,6 +463,34 @@ static void closest_rot_mat(struct rvec vec1, struct rvec vec2, rotation_around_axis(axis, bestAngle, twizzle); } +/* +static double matrix_angle(gsl_matrix *m) +{ + double a = gsl_matrix_get(m, 0, 0); + double b = gsl_matrix_get(m, 1, 1); + double c = gsl_matrix_get(m, 2, 2); + + double cos_t = (a + b + c - 1) / 2; + double theta = acos(cos_t); + + return theta; +} + +static struct rvec matrix_axis(gsl_matrix *a) +{ + double ang = matrix_angle(a); + double cos_t = cos(ang); + double p = gsl_matrix_get(a, 0, 0); + double q = gsl_matrix_get(a, 1, 1); + double r = gsl_matrix_get(a, 2, 2); + double x = sqrt((p - cos_t) / (1 - cos_t)); + double y = sqrt((q - cos_t) / (1 - cos_t)); + double z = sqrt((r - cos_t) / (1 - cos_t)); + + struct rvec v = new_rvec(x, y, z); + return v; +} +*/ static double matrix_trace(gsl_matrix *a) { @@ -487,9 +681,6 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2, gsl_matrix *fullMat; rvec_to_gsl(cell->vec1Tmp, cell2); -// gsl_vector *cell2v = rvec_to_gsl(cell2); -// gsl_vector *cell2vr = gsl_vector_calloc(3); - normalise_rvec(&obs1); normalise_rvec(&obs2); normalise_rvec(&cell1); @@ -553,96 +744,118 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, } -static int obs_vecs_match_angles(struct SpotVec *her_obs, - struct SpotVec *his_obs, - int **her_match_idxs, int **his_match_idxs, - int *match_count, struct TakeTwoCell *cell) +static int obs_vecs_match_angles(int her, int his, + struct Seed **seeds, int *match_count, + struct TakeTwoCell *cell) { - int i, j; - *match_count = 0; + struct SpotVec *obs_vecs = cell->obs_vecs; + struct SpotVec *her_obs = &obs_vecs[her]; + struct SpotVec *his_obs = &obs_vecs[his]; + + *match_count = 0; + + double min_angle = deg2rad(2.5); + double max_angle = deg2rad(187.5); - double min_angle = deg2rad(2.5); - double max_angle = deg2rad(187.5); + /* calculate angle between observed vectors */ + double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec); - /* calculate angle between observed vectors */ - double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec); + /* calculate angle between all potential theoretical vectors */ + + int i, j; + for ( i=0; i<her_obs->match_num; i++ ) { + for ( j=0; j<his_obs->match_num; j++ ) { + double score = 0; - /* calculate angle between all potential theoretical vectors */ + struct rvec *her_match = &her_obs->matches[i].vec; + struct rvec *his_match = &his_obs->matches[j].vec; - for ( i=0; i<her_obs->match_num; i++ ) { - for ( j=0; j<his_obs->match_num; j++ ) { + double her_dist = rvec_length(*her_match); + double his_dist = rvec_length(*his_match); - struct rvec *her_match = &her_obs->matches[i]; - struct rvec *his_match = &his_obs->matches[j]; + double theory_angle = rvec_angle(*her_match, + *his_match); - double theory_angle = rvec_angle(*her_match, - *his_match); + /* is this angle a match? */ - /* is this angle a match? */ + double angle_diff = fabs(theory_angle - obs_angle); - double angle_diff = fabs(theory_angle - obs_angle); + /* within basic tolerance? */ + if ( angle_diff != angle_diff || + angle_diff > cell->angle_tol ) { + continue; + } + + double add = angle_diff; + if (add == add) { + score += add * her_dist * his_dist; + } - if ( angle_diff < cell->angle_tol ) { + /* If the angles are too close to 0 + or 180, one axis ill-determined */ + if (theory_angle < min_angle || + theory_angle > max_angle) { + continue; + } - // in the case of a brief check only - if (!her_match_idxs || !his_match_idxs) { - return 1; - } + /* check that third vector adequately completes + * triangle */ - /* If the angles are too close to 0 - or 180, one axis ill-determined */ - if (theory_angle < min_angle || - theory_angle > max_angle) { - continue; - } + struct rvec theory_diff = diff_vec(*his_match, *her_match); + struct rvec obs_diff = diff_vec(his_obs->obsvec, + her_obs->obsvec); - // check the third vector + theory_angle = rvec_angle(*her_match, + theory_diff); + obs_angle = rvec_angle(her_obs->obsvec, obs_diff); + angle_diff = fabs(obs_angle - theory_angle); - struct rvec theory_diff = diff_vec(*his_match, *her_match); - struct rvec obs_diff = diff_vec(his_obs->obsvec, - her_obs->obsvec); + double diff_dist = rvec_length(obs_diff); - theory_angle = rvec_angle(*her_match, - theory_diff); - obs_angle = rvec_angle(her_obs->obsvec, obs_diff); + if (angle_diff > ANGLE_TOLERANCE) { + continue; + } - if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { - continue; - } + add = angle_diff; + if (add == add) { + score += add * her_dist * diff_dist; + } - theory_angle = rvec_angle(*his_match, - theory_diff); - obs_angle = rvec_angle(his_obs->obsvec, obs_diff); + theory_angle = rvec_angle(*his_match, + theory_diff); + obs_angle = rvec_angle(his_obs->obsvec, obs_diff); - if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { - continue; - } + if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { + continue; + } - size_t new_size = (*match_count + 1) * - sizeof(int); - if (her_match_idxs && his_match_idxs) - { - /* Reallocate the array to fit in another match */ - int *temp_hers; - int *temp_his; - temp_hers = realloc(*her_match_idxs, new_size); - temp_his = realloc(*his_match_idxs, new_size); + add = angle_diff; + if (add == add) { + score += add * his_dist * diff_dist; + } - if ( temp_hers == NULL || temp_his == NULL ) { - apologise(); - } + /* we add a new seed to the array */ + size_t new_size = (*match_count + 1); + new_size *= sizeof(struct Seed); - (*her_match_idxs) = temp_hers; - (*his_match_idxs) = temp_his; + /* Reallocate the array to fit in another match */ + struct Seed *tmp_seeds = realloc(*seeds, new_size); - (*her_match_idxs)[*match_count] = i; - (*his_match_idxs)[*match_count] = j; - } + if ( tmp_seeds == NULL ) { + apologise(); + } + + (*seeds) = tmp_seeds; + + (*seeds)[*match_count].obs1 = her; + (*seeds)[*match_count].obs2 = his; + (*seeds)[*match_count].idx1 = i; + (*seeds)[*match_count].idx2 = j; + (*seeds)[*match_count].score = score * 1000; (*match_count)++; } } - } return (*match_count > 0); } @@ -672,8 +885,8 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, struct rvec i_obsvec = i_vec.obsvec; struct rvec j_obsvec = j_vec.obsvec; - struct rvec i_cellvec = i_vec.matches[match_members[i]]; - struct rvec j_cellvec = j_vec.matches[match_members[j]]; + struct rvec i_cellvec = i_vec.matches[match_members[i]].vec; + struct rvec j_cellvec = j_vec.matches[match_members[j]].vec; rotations[count] = generate_rot_mat(i_obsvec, j_obsvec, i_cellvec, j_cellvec, @@ -716,10 +929,26 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, return 1; } -static int weed_duplicate_matches(struct SpotVec *her_obs, - struct SpotVec *his_obs, - int **her_match_idxs, int **his_match_idxs, - int *match_count, struct TakeTwoCell *cell) +gsl_matrix *rot_mat_from_indices(int her, int his, + int her_match, int his_match, + struct TakeTwoCell *cell) +{ + struct SpotVec *obs_vecs = cell->obs_vecs; + struct SpotVec *her_obs = &obs_vecs[her]; + struct SpotVec *his_obs = &obs_vecs[his]; + struct rvec i_obsvec = her_obs->obsvec; + struct rvec j_obsvec = his_obs->obsvec; + struct rvec i_cellvec = her_obs->matches[her_match].vec; + struct rvec j_cellvec = his_obs->matches[his_match].vec; + + gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, + i_cellvec, j_cellvec, cell); + + return mat; +} + +static int weed_duplicate_matches(struct Seed **seeds, + int *match_count, struct TakeTwoCell *cell) { int num_occupied = 0; gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *)); @@ -731,19 +960,18 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, } signed int i, j; + int duplicates = 0; - for (i = *match_count - 1; i >= 0; i--) { - int her_match = (*her_match_idxs)[i]; - int his_match = (*his_match_idxs)[i]; + /* Now we weed out the self-duplicates from the remaining batch */ - struct rvec i_obsvec = her_obs->obsvec; - struct rvec j_obsvec = his_obs->obsvec; - struct rvec i_cellvec = her_obs->matches[her_match]; - struct rvec j_cellvec = his_obs->matches[his_match]; + for (i = *match_count - 1; i >= 0; i--) { + int her_match = (*seeds)[i].idx1; + int his_match = (*seeds)[i].idx2; - gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, - i_cellvec, j_cellvec, cell); + gsl_matrix *mat; + mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2, + her_match, his_match, cell); int found = 0; @@ -752,8 +980,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, symm_rot_mats_are_similar(old_mats[j], mat, cell)) { // we have found a duplicate, so flag as bad. - (*her_match_idxs)[i] = -1; - (*his_match_idxs)[i] = -1; + (*seeds)[i].idx1 = -1; + (*seeds)[i].idx2 = -1; found = 1; duplicates++; @@ -785,7 +1013,7 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, int *match_members, int start, int member_num, int *match_found, struct TakeTwoCell *cell) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); + struct SpotVec *obs_vecs = cell->obs_vecs; int obs_vec_count = cell->obs_vec_count; gsl_matrix *sub = gsl_matrix_calloc(3, 3); gsl_matrix *mul = gsl_matrix_calloc(3, 3); @@ -794,42 +1022,42 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, for ( i=start; i<obs_vec_count; i++ ) { + /* If we've considered this vector before, ignore it */ + if (obs_vecs[i].in_network == 1) + { + continue; + } + /* first we check for a shared spot - harshest condition */ int shared = obs_shares_spot_w_array(obs_vecs, i, obs_members, member_num); if ( !shared ) continue; - int skip = 0; - for ( j=0; j<member_num && skip == 0; j++ ) { - if (i == obs_members[j]) { - skip = 1; - } - } - - if (skip) { - continue; - } - int all_ok = 1; int matched = -1; + /* Check all existing members are happy to let in the newcomer */ for ( j=0; j<member_num && all_ok; j++ ) { struct SpotVec *me = &obs_vecs[i]; struct SpotVec *you = &obs_vecs[obs_members[j]]; - struct rvec you_cell = you->matches[match_members[j]]; + struct rvec you_cell; + you_cell = you->matches[match_members[j]].vec; struct rvec me_obs = me->obsvec; struct rvec you_obs = you->obsvec; int one_is_okay = 0; + /* Loop through all possible theoretical vector + * matches for the newcomer.. */ + for ( k=0; k<me->match_num; k++ ) { gsl_matrix *test_rot; - struct rvec me_cell = me->matches[k]; + struct rvec me_cell = me->matches[k].vec; test_rot = generate_rot_mat(me_obs, you_obs, me_cell, you_cell, @@ -844,6 +1072,11 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, if (ok) { one_is_okay = 1; + /* We are only happy if the vector + * matches for only one kind of + * theoretical vector. We don't want to + * accept mixtures of theoretical vector + * matches. */ if (matched >= 0 && k == matched) { *match_found = k; } else if (matched < 0) { @@ -863,12 +1096,6 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, if (all_ok) { - - for ( j=0; j<member_num; j++ ) { - // STATUS("%i ", obs_members[j]); - } - //STATUS("%i\n", i); - return i; } } @@ -878,16 +1105,212 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, return -1; } +/** + * Reward target function for refining solution to all vectors in a + * given image. Sum of exponentials of the negative distances, which + * means that the reward decays as the distance from the nearest + * theoretical vector reduces. */ +static double obs_to_sol_score(struct TakeTwoCell *ttCell) +{ + double total = 0; + int count = 0; + int i; + gsl_matrix *solution = ttCell->solution; + gsl_matrix *full_rot = gsl_matrix_alloc(3, 3); + rotate_gsl_by_angles(solution, ttCell->x_ang, ttCell->y_ang, + ttCell->z_ang, full_rot); + + for (i = 0; i < ttCell->obs_vec_count; i++) + { + struct rvec *obs = &ttCell->obs_vecs[i].obsvec; + rvec_to_gsl(ttCell->vec1Tmp, *obs); + + /* Rotate all the observed vectors by the modified soln */ + /* ttCell->vec2Tmp = 1.0 * full_rot * ttCell->vec1Tmp */ + gsl_blas_dgemv(CblasTrans, 1.0, full_rot, ttCell->vec1Tmp, + 0.0, ttCell->vec2Tmp); + struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); + + int j = ttCell->obs_vecs[i].assignment; + + if (j < 0) continue; + + struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; + struct rvec diff = diff_vec(rotated, *match); + + double length = rvec_length(diff); + + double addition = exp(-(1 / RECIP_TOLERANCE) * length); + total += addition; + count++; + } + + total /= (double)-count; + + gsl_matrix_free(full_rot); + + return total; +} + +/** + * Matches every observed vector in the image to its closest theoretical + * neighbour after applying the rotation matrix, in preparation for + * refining the rotation matrix to match these. */ +static void match_all_obs_to_sol(struct TakeTwoCell *ttCell) +{ + int i, j; + double total = 0; + int count = 0; + gsl_matrix *solution = ttCell->solution; + + for (i = 0; i < ttCell->obs_vec_count; i++) + { + struct rvec *obs = &ttCell->obs_vecs[i].obsvec; + rvec_to_gsl(ttCell->vec1Tmp, *obs); + + /* ttCell->vec2Tmp = 1.0 * solution * ttCell->vec1Tmp */ + gsl_blas_dgemv(CblasTrans, 1.0, solution, ttCell->vec1Tmp, + 0.0, ttCell->vec2Tmp); + struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); + + double smallest = FLT_MAX; + int assigned = -1; + + for (j = 0; j < ttCell->obs_vecs[i].match_num; j++) + { + struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; + struct rvec diff = diff_vec(rotated, *match); + + double length = rvec_length(diff); + if (length < smallest) + { + smallest = length; + assigned = j; + } + } + + ttCell->obs_vecs[i].assignment = assigned; + + if (smallest != FLT_MAX) + { + double addition = exp(-(1 / RECIP_TOLERANCE) * smallest); + total += addition; + count++; -static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, - int match_idx1, int match_idx2, int *max_members, - struct TakeTwoCell *cell) + } + } + + total /= (double)count; +} + +/** + * Refines a matrix against all of the observed vectors against their + * closest theoretical neighbour, by perturbing the matrix along the principle + * axes until it maximises a reward function consisting of the sum of + * the distances of individual observed vectors to their closest + * theoretical neighbour. Reward function means that noise and alternative + * lattices do not dominate the equation! + **/ +static void refine_solution(struct TakeTwoCell *ttCell) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); + match_all_obs_to_sol(ttCell); + + int i, j, k; + const int total = 3 * 3 * 3; + const int middle = (total - 1) / 2; + + struct rvec steps[total]; + double start = obs_to_sol_score(ttCell); + const int max_tries = 100; + + int count = 0; + double size = ANGLE_STEP_SIZE; + + /* First we create our combinations of steps */ + for (i = -1; i <= 1; i++) { + for (j = -1; j <= 1; j++) { + for (k = -1; k <= 1; k++) { + struct rvec vec = new_rvec(i, j, k); + steps[count] = vec; + count++; + } + } + } + + struct rvec current = new_rvec(ttCell->x_ang, ttCell->y_ang, + ttCell->z_ang); + + double best = start; + count = 0; + + while (size > ANGLE_CONVERGE_SIZE && count < max_tries) + { + struct rvec sized[total]; + + int best_num = middle; + for (i = 0; i < total; i++) + { + struct rvec sized_step = steps[i]; + sized_step.u *= size; + sized_step.v *= size; + sized_step.w *= size; + + struct rvec new_angles = rvec_add_rvec(current, + sized_step); + + sized[i] = new_angles; + + ttCell->x_ang = sized[i].u; + ttCell->y_ang = sized[i].v; + ttCell->z_ang = sized[i].w; + + double score = obs_to_sol_score(ttCell); + + if (score < best) + { + best = score; + best_num = i; + } + } + + if (best_num == middle) + { + size /= 2; + } + + current = sized[best_num]; + count++; + } + + ttCell->x_ang = 0; + ttCell->y_ang = 0; + ttCell->z_ang = 0; + + gsl_matrix *tmp = gsl_matrix_alloc(3, 3); + rotate_gsl_by_angles(ttCell->solution, current.u, + current.v, current.w, tmp); + gsl_matrix_free(ttCell->solution); + ttCell->solution = tmp; +} + + +static unsigned int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, + int match_idx1, int match_idx2, + struct TakeTwoCell *cell) +{ + + struct SpotVec *obs_vecs = cell->obs_vecs; int obs_vec_count = cell->obs_vec_count; int *obs_members; int *match_members; + /* Clear the in_network status of all vectors to start */ + int i; + for (i = 0; i < obs_vec_count; i++) + { + obs_vecs[i].in_network = 0; + } + /* indices of members of the self-consistent network of vectors */ obs_members = malloc((cell->member_thresh+3)*sizeof(int)); match_members = malloc((cell->member_thresh+3)*sizeof(int)); @@ -902,7 +1325,6 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, match_members[0] = match_idx1; match_members[1] = match_idx2; int member_num = 2; - *max_members = 2; /* counter for dead ends which must not exceed MAX_DEAD_ENDS * before it is reset in an additional branch */ @@ -923,7 +1345,7 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, signed int next_index = find_next_index(rot, obs_members, match_members, - start, member_num, + 0, member_num, &match_found, cell); if ( member_num < 2 ) { @@ -942,25 +1364,19 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, /* We have not had too many dead ends. Try removing the last member and continue. */ - start = obs_members[member_num - 1] + 1; member_num--; dead_ends++; continue; } - /* we have elongated membership - so reset dead_ends counter */ - // dead_ends = 0; - + /* Elongation of the network was successful */ + obs_vecs[next_index].in_network = 1; obs_members[member_num] = next_index; match_members[member_num] = match_found; member_num++; - if (member_num > *max_members) { - *max_members = member_num; - } - /* If member_num is high enough, we want to return a yes */ if ( member_num > cell->member_thresh ) break; @@ -972,126 +1388,210 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, free(obs_members); free(match_members); - return ( member_num > cell->member_thresh ); + return ( member_num ); } -static int start_seed(int i, int j, int i_match, int j_match, - gsl_matrix **rotation, int *max_members, - struct TakeTwoCell *cell) +static unsigned int start_seed(int i, int j, int i_match, int j_match, + gsl_matrix **rotation, struct TakeTwoCell *cell) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); + struct SpotVec *obs_vecs = cell->obs_vecs; gsl_matrix *rot_mat; rot_mat = generate_rot_mat(obs_vecs[i].obsvec, obs_vecs[j].obsvec, - obs_vecs[i].matches[i_match], - obs_vecs[j].matches[j_match], + obs_vecs[i].matches[i_match].vec, + obs_vecs[j].matches[j_match].vec, cell); /* Try to expand this rotation matrix to a larger network */ - int success = grow_network(rot_mat, i, j, i_match, j_match, max_members, - cell); + int member_num = grow_network(rot_mat, i, j, i_match, j_match, + cell); /* return this matrix and if it was immediately successful */ *rotation = rot_mat; - return success; + return member_num; } -static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell) +static int sort_seed_by_score(const void *av, const void *bv) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); - int obs_vec_count = cell->obs_vec_count; + struct Seed *a = (struct Seed *)av; + struct Seed *b = (struct Seed *)bv; + return a->score > b->score; +} + +static void remove_old_solutions(struct TakeTwoCell *cell, + struct taketwo_private *tp) +{ + int duplicates = 0; + struct Seed *seeds = cell->seeds; + unsigned int total = cell->seed_count; + + /* First we remove duplicates with previous solutions */ + + int i, j; + for (i = total - 1; i >= 0; i--) { + int her_match = seeds[i].idx1; + int his_match = seeds[i].idx2; + + gsl_matrix *mat; + mat = rot_mat_from_indices(seeds[i].obs1, seeds[i].obs2, + her_match, his_match, cell); + + if (mat == NULL) + { + continue; + } + + for (j = 0; j < tp->numPrevs; j++) + { + int sim = symm_rot_mats_are_similar(tp->prevSols[j], + mat, cell); + + /* Found a duplicate with a previous solution */ + if (sim) + { + seeds[i].idx1 = -1; + seeds[i].idx2 = -1; + duplicates++; + break; + } + } - /* META: Maximum achieved maximum network membership */ - int max_max_members = 0; - gsl_matrix *best_rotation = NULL; + gsl_matrix_free(mat); + } + +// STATUS("Removing %i duplicates due to prev solutions.\n", duplicates); +} -// unsigned long start_time = time(NULL); +static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) +{ + struct SpotVec *obs_vecs = cell->obs_vecs; + int obs_vec_count = cell->obs_vec_count; /* loop round pairs of vectors to try and find a suitable * seed to start building a self-consistent network of vectors */ int i, j; - for ( i=0; i<obs_vec_count; i++ ) { + for ( i=1; i<obs_vec_count; i++ ) { + for ( j=0; j<i; j++ ) { + /** Only check distances which are accumulatively less + * than the limit if we can easily generate seeds */ + if (obs_vecs[j].distance + obs_vecs[i].distance > + MAX_RECIP_DISTANCE && cell->seed_count > 100) { + continue; + } + /** Check to see if there is a shared spot - opportunity * for optimisation by generating a look-up table * by spot instead of by vector. */ - int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]); - + int shared = obs_vecs_share_spot(&obs_vecs[i], + &obs_vecs[j]); if ( !shared ) continue; /* cell vector index matches stored in i, j and total * number stored in int matches. */ - int *i_idx = 0; - int *j_idx = 0; - int matches; - - /* Check to see if any angles match from the cell vectors */ - obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j], - &i_idx, &j_idx, &matches, cell); - if ( matches == 0 ) { - free(i_idx); - free(j_idx); + int seed_num = 0; + struct Seed *seeds = NULL; + + /* Check to see if any angles match from the cell + * vectors */ + obs_vecs_match_angles(i, j, &seeds, &seed_num, cell); + + if (seed_num == 0) + { + /* Nothing to clean up here */ continue; } /* Weed out the duplicate seeds (from symmetric - * reflection pairs) - */ + * reflection pairs) */ + weed_duplicate_matches(&seeds, &seed_num, cell); + + /* Add all the new seeds to the full list */ + + size_t new_size = cell->seed_count + seed_num; + new_size *= sizeof(struct Seed); + + struct Seed *tmp = realloc(cell->seeds, new_size); - weed_duplicate_matches(&obs_vecs[i], &obs_vecs[j], - &i_idx, &j_idx, &matches, cell); + if (tmp == NULL) { + apologise(); + } + + cell->seeds = tmp; - /* We have seeds! Pass each of them through the seed-starter */ - /* If a seed has the highest achieved membership, make note...*/ - int k; - for ( k=0; k<matches; k++ ) { - if (i_idx[k] < 0 || j_idx[k] < 0) { + for (int i = 0; i < seed_num; i++) + { + if (seeds[i].idx1 < 0 || seeds[i].idx2 < 0) + { continue; } - int max_members = 0; - - int success = start_seed(i, j, - i_idx[k], j_idx[k], - rotation, &max_members, - cell); - - if (success) { - free(i_idx); free(j_idx); - gsl_matrix_free(best_rotation); - return success; - } else { - if (max_members > max_max_members) { - max_max_members = max_members; - gsl_matrix_free(best_rotation); - best_rotation = *rotation; - *rotation = NULL; - } else { - gsl_matrix_free(*rotation); - *rotation = NULL; - } - } + cell->seeds[cell->seed_count] = seeds[i]; + cell->seed_count++; } + } + } + + qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), + sort_seed_by_score); + + + return 1; +} + +static unsigned int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) +{ + struct Seed *seeds = cell->seeds; + int seed_num = cell->seed_count; + int member_num = 0; + int max_members = 0; + gsl_matrix *rot = NULL; + + /* We have seeds! Pass each of them through the seed-starter */ + /* If a seed has the highest achieved membership, make note...*/ + int k; + for ( k=0; k<seed_num; k++ ) { + int seed_idx1 = seeds[k].idx1; + int seed_idx2 = seeds[k].idx2; + + if (seed_idx1 < 0 || seed_idx2 < 0) { + continue; + } + + int seed_obs1 = seeds[k].obs1; + int seed_obs2 = seeds[k].obs2; - free(i_idx); - free(j_idx); + member_num = start_seed(seed_obs1, seed_obs2, seed_idx1, + seed_idx2, &rot, cell); + + if (member_num > max_members) + { + *rotation = rot; + max_members = member_num; + } + + if (member_num >= NETWORK_MEMBER_THRESHOLD) { + free(seeds); + return max_members; } - } /* yes this } is meant to be here */ + } + + free(seeds); - *rotation = best_rotation; - return (best_rotation != NULL); + return max_members; } + static void set_gsl_matrix(gsl_matrix *mat, double asx, double asy, double asz, double bsx, double bsy, double bsz, double csx, double csy, double csz) @@ -1184,25 +1684,24 @@ static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell) return 1; } - struct sortme { - struct rvec v; + struct TheoryVec v; double dist; }; -static int sort_func(const void *av, const void *bv) +static int sort_theory_distances(const void *av, const void *bv) { struct sortme *a = (struct sortme *)av; struct sortme *b = (struct sortme *)bv; return a->dist > b->dist; } - -static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, - struct SpotVec *obs_vecs, int obs_vec_count, - int is_asymmetric, struct TakeTwoCell *cell) +static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count, + struct TakeTwoCell *cell) { + struct SpotVec *obs_vecs = cell->obs_vecs; + int obs_vec_count = cell->obs_vec_count; int i, j; for ( i=0; i<obs_vec_count; i++ ) { @@ -1212,7 +1711,7 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, for ( j=0; j<cell_vec_count; j++ ) { /* get distance for unit cell vector */ - double cell_length = rvec_length(cell_vecs[j]); + double cell_length = rvec_length(cell_vecs[j].vec); double obs_length = obs_vecs[i].distance; /* check if this matches the observed length */ @@ -1235,20 +1734,15 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, /* Pointers to relevant things */ - struct rvec **match_array; + struct TheoryVec **match_array; int *match_count; - if (!is_asymmetric) { - match_array = &(obs_vecs[i].matches); - match_count = &(obs_vecs[i].match_num); - } else { - match_array = &(obs_vecs[i].asym_matches); - match_count = &(obs_vecs[i].asym_match_num); - } + match_array = &(obs_vecs[i].matches); + match_count = &(obs_vecs[i].match_num); /* Sort in order to get most agreeable matches first */ - qsort(for_sort, count, sizeof(struct sortme), sort_func); - *match_array = malloc(count*sizeof(struct rvec)); + qsort(for_sort, count, sizeof(struct sortme), sort_theory_distances); + *match_array = malloc(count*sizeof(struct TheoryVec)); *match_count = count; for ( j=0; j<count; j++ ) { (*match_array)[j] = for_sort[j].v; @@ -1268,7 +1762,7 @@ static int compare_spot_vecs(const void *av, const void *bv) } static int gen_observed_vecs(struct rvec *rlps, int rlp_count, - struct SpotVec **obs_vecs, int *obs_vec_count) + struct TakeTwoCell *cell) { int i, j; int count = 0; @@ -1279,7 +1773,6 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, for ( i=0; i<rlp_count-1 && count < MAX_OBS_VECTORS; i++ ) { for ( j=i+1; j<rlp_count; j++ ) { - /* calculate difference vector between rlps */ struct rvec diff = diff_vec(rlps[i], rlps[j]); @@ -1291,13 +1784,13 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, count++; struct SpotVec *temp_obs_vecs; - temp_obs_vecs = realloc(*obs_vecs, + temp_obs_vecs = realloc(cell->obs_vecs, count*sizeof(struct SpotVec)); if ( temp_obs_vecs == NULL ) { return 0; } else { - *obs_vecs = temp_obs_vecs; + cell->obs_vecs = temp_obs_vecs; /* initialise all SpotVec struct members */ @@ -1305,29 +1798,27 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, spot_vec.obsvec = diff; spot_vec.distance = sqrt(sqlength); spot_vec.matches = NULL; + spot_vec.assignment = -1; spot_vec.match_num = 0; spot_vec.her_rlp = &rlps[i]; spot_vec.his_rlp = &rlps[j]; - (*obs_vecs)[count - 1] = spot_vec; + cell->obs_vecs[count - 1] = spot_vec; } } } - /* Sort such that the shortest and least error-prone distances - are searched first. - */ - qsort(*obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs); + /* Sort such that the shortest distances are searched first. */ + qsort(cell->obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs); - *obs_vec_count = count; + cell->obs_vec_count = count; return 1; } -static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs, - struct rvec **asym_vecs, int *vec_count, - int *asym_vec_count) +static int gen_theoretical_vecs(UnitCell *cell, struct TheoryVec **cell_vecs, + unsigned int *vec_count) { double a, b, c, alpha, beta, gamma; int h_max, k_max, l_max; @@ -1351,7 +1842,6 @@ static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs, int h, k, l; int _h, _k, _l; int count = 0; - int asym_count = 0; for ( h=-h_max; h<=+h_max; h++ ) { for ( k=-k_max; k<=+k_max; k++ ) { @@ -1367,50 +1857,37 @@ static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs, if (h == _h && k == _k && l == _l) { asymmetric = 1; - asym_count++; } cell_vec.u = h*asx + k*bsx + l*csx; cell_vec.v = h*asy + k*bsy + l*csy; cell_vec.w = h*asz + k*bsz + l*csz; + struct TheoryVec theory; + theory.vec = cell_vec; + theory.asym = asymmetric; + /* add this to our array - which may require expanding */ count++; - struct rvec *temp_cell_vecs; - temp_cell_vecs = realloc(*cell_vecs, count*sizeof(struct rvec)); - struct rvec *temp_asym_vecs = NULL; - - if (asymmetric) { - temp_asym_vecs = realloc(*asym_vecs, - count*sizeof(struct rvec)); - } + struct TheoryVec *temp_cell_vecs; + temp_cell_vecs = realloc(*cell_vecs, + count*sizeof(struct TheoryVec)); if ( temp_cell_vecs == NULL ) { return 0; - } else if (asymmetric && temp_asym_vecs == NULL) { - return 0; } else { *cell_vecs = temp_cell_vecs; - (*cell_vecs)[count - 1] = cell_vec; - - if (!asymmetric) { - continue; - } - - *asym_vecs = temp_asym_vecs; - (*asym_vecs)[asym_count - 1] = cell_vec; + (*cell_vecs)[count - 1] = theory; } } } } *vec_count = count; - *asym_vec_count = asym_count; free_symoplist(rawList); - return 1; } @@ -1425,7 +1902,6 @@ static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs, int i; for ( i=0; i<obs_vec_count; i++ ) { free(obs_vecs[i].matches); - free(obs_vecs[i].asym_matches); } free(obs_vecs); @@ -1433,11 +1909,16 @@ static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs, static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) { + /* n.b. solutions in ttCell are taken care of in the + * partial taketwo cleanup. */ int i; for ( i=0; i<ttCell->numOps; i++ ) { gsl_matrix_free(ttCell->rotSymOps[i]); } + cleanup_taketwo_obs_vecs(ttCell->obs_vecs, + ttCell->obs_vec_count); + free(ttCell->vec1Tmp); free(ttCell->vec2Tmp); free(ttCell->twiz1Tmp); @@ -1459,39 +1940,38 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) * successful. **/ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, - struct rvec *rlps, int rlp_count) + struct rvec *rlps, int rlp_count, + struct taketwo_private *tp) { - int cell_vec_count = 0; - int asym_vec_count = 0; - struct rvec *cell_vecs = NULL; - struct rvec *asym_vecs = NULL; UnitCell *result; - int obs_vec_count = 0; - struct SpotVec *obs_vecs = NULL; int success = 0; gsl_matrix *solution = NULL; + /* Initialise TakeTwoCell */ struct TakeTwoCell ttCell; ttCell.cell = cell; + ttCell.seeds = NULL; + ttCell.seed_count = 0; ttCell.rotSymOps = NULL; + ttCell.obs_vecs = NULL; ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3); ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3); ttCell.vec1Tmp = gsl_vector_calloc(3); ttCell.vec2Tmp = gsl_vector_calloc(3); ttCell.numOps = 0; + ttCell.obs_vec_count = 0; + ttCell.solution = NULL; + ttCell.x_ang = 0; + ttCell.y_ang = 0; + ttCell.z_ang = 0; success = generate_rotation_sym_ops(&ttCell); - success = gen_theoretical_vecs(cell, &cell_vecs, &asym_vecs, - &cell_vec_count, &asym_vec_count); if ( !success ) return NULL; - success = gen_observed_vecs(rlps, rlp_count, &obs_vecs, &obs_vec_count); + success = gen_observed_vecs(rlps, rlp_count, &ttCell); if ( !success ) return NULL; - ttCell.obs_vecs = &obs_vecs; - ttCell.obs_vec_count = obs_vec_count; - if ( opts->member_thresh < 0 ) { ttCell.member_thresh = NETWORK_MEMBER_THRESHOLD; } else { @@ -1516,32 +1996,79 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, ttCell.trace_tol = sqrt(4.0*(1.0-cos(opts->trace_tol))); } - success = match_obs_to_cell_vecs(asym_vecs, asym_vec_count, - obs_vecs, obs_vec_count, 1, &ttCell); - - success = match_obs_to_cell_vecs(cell_vecs, cell_vec_count, - obs_vecs, obs_vec_count, 0, &ttCell); - - free(cell_vecs); - free(asym_vecs); + success = match_obs_to_cell_vecs(tp->theory_vecs, tp->vec_count, + &ttCell); if ( !success ) return NULL; - find_seed(&solution, &ttCell); + /* Find all the seeds, then take each one and extend them, returning + * a solution if it exceeds the NETWORK_MEMBER_THRESHOLD. */ + find_seeds(&ttCell, tp); + remove_old_solutions(&ttCell, tp); + unsigned int members = start_seeds(&solution, &ttCell); if ( solution == NULL ) { - cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); return NULL; } + /* If we have a solution, refine against vectors in the entire image */ + ttCell.solution = solution; + refine_solution(&ttCell); + solution = ttCell.solution; + double score = obs_to_sol_score(&ttCell); + + /* Add the current solution to the previous solutions list */ + int new_size = (tp->numPrevs + 1) * sizeof(gsl_matrix *); + gsl_matrix **tmp = realloc(tp->prevSols, new_size); + double *tmpScores = realloc(tp->prevScores, + (tp->numPrevs + 1) * sizeof(double)); + unsigned int *tmpSuccesses; + tmpSuccesses = realloc(tp->membership, + (tp->numPrevs + 1) * sizeof(unsigned int)); + + if (!tmp) { + apologise(); + } + + tp->prevSols = tmp; + tp->prevScores = tmpScores; + tp->membership = tmpSuccesses; + + tp->prevSols[tp->numPrevs] = solution; + tp->prevScores[tp->numPrevs] = score; + tp->membership[tp->numPrevs] = members; + tp->numPrevs++; + + /* Prepare the solution for CrystFEL friendliness */ result = transform_cell_gsl(cell, solution); - gsl_matrix_free(solution); - cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); cleanup_taketwo_cell(&ttCell); return result; } +/* Cleans up the per-image information for taketwo */ + +static void partial_taketwo_cleanup(struct taketwo_private *tp) +{ + if (tp->prevSols != NULL) + { + for (int i = 0; i < tp->numPrevs; i++) + { + gsl_matrix_free(tp->prevSols[i]); + } + + free(tp->prevSols); + } + + free(tp->prevScores); + free(tp->membership); + tp->prevScores = NULL; + tp->membership = NULL; + tp->xtal_num = 0; + tp->numPrevs = 0; + tp->prevSols = NULL; + +} /* CrystFEL interface hooks */ @@ -1555,6 +2082,34 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts, int i; struct taketwo_private *tp = (struct taketwo_private *)priv; + /* Check serial number against previous for solution tracking */ + int this_serial = image->serial; + + if (tp->serial_num == this_serial) + { + tp->xtal_num = image->n_crystals; + } + else + { + /* + for (i = 0; i < tp->numPrevs; i++) + { + STATUS("score, %i, %.5f, %i\n", + this_serial, tp->prevScores[i], + tp->membership[i]); + } + */ + + partial_taketwo_cleanup(tp); + tp->serial_num = this_serial; + tp->xtal_num = image->n_crystals; + } + + /* + STATUS("Indexing %i with %i attempts, %i crystals\n", this_serial, tp->attempts, + image->n_crystals); + */ + rlps = malloc((image_feature_count(image->features)+1)*sizeof(struct rvec)); for ( i=0; i<image_feature_count(image->features); i++ ) { struct imagefeature *pk = image_get_feature(image->features, i); @@ -1568,7 +2123,7 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts, rlps[n_rlps].v = 0.0; rlps[n_rlps++].w = 0.0; - cell = run_taketwo(tp->cell, opts, rlps, n_rlps); + cell = run_taketwo(tp->cell, opts, rlps, n_rlps, tp); free(rlps); if ( cell == NULL ) return 0; @@ -1631,14 +2186,27 @@ void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell) tp->cell = cell; tp->indm = *indm; + tp->serial_num = -1; + tp->xtal_num = 0; + tp->prevSols = NULL; + tp->numPrevs = 0; + tp->prevScores = NULL; + tp->membership = NULL; + tp->vec_count = 0; + tp->theory_vecs = NULL; + + gen_theoretical_vecs(cell, &tp->theory_vecs, &tp->vec_count); return tp; } - void taketwo_cleanup(IndexingPrivate *pp) { struct taketwo_private *tp = (struct taketwo_private *)pp; + + partial_taketwo_cleanup(tp); + free(tp->theory_vecs); + free(tp); } diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 41580e41..b28dc93a 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -3,12 +3,12 @@ * * Invoke xds for crystal autoindexing * - * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2013 Cornelius Gati * * Authors: - * 2010-2017 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * 2013 Cornelius Gati <cornelius.gati@cfel.de> * * This file is part of CrystFEL. @@ -54,7 +54,8 @@ #include "cell-utils.h" -#define XDS_VERBOSE 0 +/* Fake pixel size and camera length, both in metres */ +#define FAKE_PIXEL_SIZE (70e-6) #define FAKE_CLEN (0.1) @@ -66,142 +67,7 @@ struct xds_private }; -struct xds_data { - - /* Low-level stuff */ - int pty; - pid_t pid; - char *rbuffer; - int rbufpos; - int rbuflen; - - /* High-level stuff */ - int step; - int finished_ok; - UnitCell *target_cell; -}; - -static void xds_parseline(const char *line, struct image *image, - struct xds_data *xds) -{ -#if XDS_VERBOSE - char *copy; - int i; - - copy = strdup(line); - for ( i=0; i<strlen(copy); i++ ) { - if ( copy[i] == '\r' ) copy[i]='r'; - if ( copy[i] == '\n' ) copy[i]='\0'; - } - STATUS("XDS: %s\n", copy); - free(copy); -#endif -} - - -static int xds_readable(struct image *image, struct xds_data *xds) -{ - int rval; - int no_string = 0; - - rval = read(xds->pty, xds->rbuffer+xds->rbufpos, - xds->rbuflen-xds->rbufpos); - if ( (rval == -1) || (rval == 0) ) return 1; - - xds->rbufpos += rval; - assert(xds->rbufpos <= xds->rbuflen); - - while ( (!no_string) && (xds->rbufpos > 0) ) { - - int i; - int block_ready = 0; - - /* See if there's a full line in the buffer yet */ - for ( i=0; i<xds->rbufpos-1; i++ ) { - /* Means the last value looked at is rbufpos-2 */ - - if ( (xds->rbuffer[i] == '\r') - && (xds->rbuffer[i+1] == '\n') ) { - block_ready = 1; - break; - } - - } - - if ( block_ready ) { - - unsigned int new_rbuflen; - unsigned int endbit_length; - char *block_buffer = NULL; - - block_buffer = malloc(i+1); - memcpy(block_buffer, xds->rbuffer, i); - block_buffer[i] = '\0'; - - if ( block_buffer[0] == '\r' ) { - memmove(block_buffer, block_buffer+1, i); - } - - xds_parseline(block_buffer, image, xds); - free(block_buffer); - endbit_length = i+2; - - /* Now the block's been parsed, it should be - * forgotten about */ - memmove(xds->rbuffer, - xds->rbuffer + endbit_length, - xds->rbuflen - endbit_length); - - /* Subtract the number of bytes removed */ - xds->rbufpos = xds->rbufpos - - endbit_length; - new_rbuflen = xds->rbuflen - endbit_length; - if ( new_rbuflen == 0 ) new_rbuflen = 256; - xds->rbuffer = realloc(xds->rbuffer, - new_rbuflen); - xds->rbuflen = new_rbuflen; - - } else { - - if ( xds->rbufpos==xds->rbuflen ) { - - /* More buffer space is needed */ - xds->rbuffer = realloc( - xds->rbuffer, - xds->rbuflen + 256); - xds->rbuflen = xds->rbuflen + 256; - /* The new space gets used at the next - * read, shortly... */ - - } - no_string = 1; - - } - - } - - return 0; -} - - -static int check_cell(struct xds_private *xp, struct image *image, - UnitCell *cell) -{ - Crystal *cr; - - cr = crystal_new(); - if ( cr == NULL ) { - ERROR("Failed to allocate crystal.\n"); - return 0; - } - crystal_set_cell(cr, cell); - image_add_crystal(image, cr); - - return 1; -} - - -static int read_cell(struct image *image, struct xds_private *xp) +static int read_cell(struct image *image) { FILE * fh; float axstar, aystar, azstar; @@ -213,12 +79,10 @@ static int read_cell(struct image *image, struct xds_private *xp) char *rval, line[1024]; int r; UnitCell *cell; + Crystal *cr; fh = fopen("IDXREF.LP", "r"); - if ( fh == NULL ) { - ERROR("Couldn't open 'IDXREF.LP'\n"); - return 0; - } + if ( fh == NULL ) return 0; /* Not indexable */ do { rval = fgets(line, 1023, fh); @@ -236,30 +100,41 @@ static int read_cell(struct image *image, struct xds_private *xp) fclose(fh); return 0; } + + /* Get first vector */ rval = fgets(line, 1023, fh); if ( rval == NULL ) { fclose(fh); return 0; } - + if ( line[4] != '1' ) { + ERROR("No first vector from XDS.\n"); + return 0; + } memcpy(asx, line+7, 10); asx[10] = '\0'; memcpy(asy, line+17, 10); asy[10] = '\0'; memcpy(asz, line+27, 10); asz[10] = '\0'; + /* Get second vector */ rval = fgets(line, 1023, fh); if ( rval == NULL ) { fclose(fh); return 0; } - + if ( line[4] != '2' ) { + ERROR("No second vector from XDS.\n"); + return 0; + } memcpy(bsx, line+7, 10); bsx[10] = '\0'; memcpy(bsy, line+17, 10); bsy[10] = '\0'; memcpy(bsz, line+27, 10); bsz[10] = '\0'; + /* Get third vector */ rval = fgets(line, 1023, fh); fclose(fh); if ( rval == NULL ) return 0; - + if ( line[4] != '3' ) return 0; /* No error message this time + * - happens a lot */ memcpy(csx, line+7, 10); csx[10] = '\0'; memcpy(csy, line+17, 10); csy[10] = '\0'; memcpy(csz, line+27, 10); csz[10] = '\0'; @@ -284,9 +159,16 @@ static int read_cell(struct image *image, struct xds_private *xp) axstar*10e9, aystar*10e9, azstar*10e9, bxstar*10e9, bystar*10e9, bzstar*10e9, -cxstar*10e9, -cystar*10e9, -czstar*10e9); - r = check_cell(xp, image, cell); - return r; + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + crystal_set_cell(cr, cell); + image_add_crystal(image, cr); + + return 1; } @@ -322,8 +204,8 @@ static void write_spot(struct image *image) x = tan(ttx)*FAKE_CLEN; y = tan(tty)*FAKE_CLEN; - x = (x / 70e-6) + 1500; - y = (y / 70e-6) + 1500; + x = (x / FAKE_PIXEL_SIZE) + 1500; + y = (y / FAKE_PIXEL_SIZE) + 1500; fprintf(fh, "%10.2f %10.2f %10.2f %10.0f.\n", x, y, 0.5, f->intensity); @@ -450,8 +332,8 @@ static int write_inp(struct image *image, struct xds_private *xp) fprintf(fh, "NX= 3000\n"); fprintf(fh, "NY= 3000\n"); - fprintf(fh, "QX= 0.07\n"); - fprintf(fh, "QY= 0.07\n"); + fprintf(fh, "QX= %f\n", FAKE_PIXEL_SIZE*1e3); /* Pixel size in mm */ + fprintf(fh, "QY= %f\n", FAKE_PIXEL_SIZE*1e3); /* Pixel size in mm */ fprintf(fh, "INDEX_ORIGIN=0 0 0\n"); fprintf(fh, "DIRECTION_OF_DETECTOR_X-AXIS=1 0 0\n"); fprintf(fh, "DIRECTION_OF_DETECTOR_Y-AXIS=0 1 0\n"); @@ -474,46 +356,33 @@ static int write_inp(struct image *image, struct xds_private *xp) int run_xds(struct image *image, void *priv) { - unsigned int opts; int status; int rval; int n; - struct xds_data *xds; + pid_t pid; + int pty; struct xds_private *xp = (struct xds_private *)priv; - xds = malloc(sizeof(struct xds_data)); - if ( xds == NULL ) { - ERROR("Couldn't allocate memory for xds data.\n"); - return 0; - } - - xds->target_cell = xp->cell; - if ( write_inp(image, xp) ) { ERROR("Failed to write XDS.INP file for XDS.\n"); - free(xds); return 0; } n = image_feature_count(image->features); - if ( n < 25 ) { - free(xds); - return 0; - } + if ( n < 25 ) return 0; write_spot(image); /* Delete any old indexing result which may exist */ remove("IDXREF.LP"); - xds->pid = forkpty(&xds->pty, NULL, NULL, NULL); + pid = forkpty(&pty, NULL, NULL, NULL); - if ( xds->pid == -1 ) { + if ( pid == -1 ) { ERROR("Failed to fork for XDS\n"); - free(xds); return 0; } - if ( xds->pid == 0 ) { + if ( pid == 0 ) { /* Child process: invoke XDS */ struct termios t; @@ -528,66 +397,10 @@ int run_xds(struct image *image, void *priv) _exit(0); } + waitpid(pid, &status, 0); - xds->rbuffer = malloc(256); - xds->rbuflen = 256; - xds->rbufpos = 0; - - /* Set non-blocking */ - opts = fcntl(xds->pty, F_GETFL); - fcntl(xds->pty, F_SETFL, opts | O_NONBLOCK); - - //xds->step = 1; /* This starts the "initialisation" procedure */ - xds->finished_ok = 0; - - do { - - fd_set fds; - struct timeval tv; - int sval; - - FD_ZERO(&fds); - FD_SET(xds->pty, &fds); - - tv.tv_sec = 30; - tv.tv_usec = 0; - - sval = select(xds->pty+1, &fds, NULL, NULL, &tv); - - if ( sval == -1 ) { - - const int err = errno; - - switch ( err ) { - - case EINTR: - STATUS("Restarting select()\n"); - rval = 0; - break; - - default: - ERROR("select() failed: %s\n", strerror(err)); - rval = 1; - break; - - } - - } else if ( sval != 0 ) { - rval = xds_readable(image, xds); - } else { - ERROR("No response from XDS..\n"); - rval = 1; - } - - } while ( !rval ); - - close(xds->pty); - free(xds->rbuffer); - waitpid(xds->pid, &status, 0); - - rval = read_cell(image, xp); - - free(xds); + close(pty); + rval = read_cell(image); return rval; } |