aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/cell.c4
-rw-r--r--libcrystfel/src/events.c4
-rw-r--r--libcrystfel/src/geometry.c465
-rw-r--r--libcrystfel/src/geometry.h21
-rw-r--r--libcrystfel/src/image.h2
-rw-r--r--libcrystfel/src/index.c8
-rw-r--r--libcrystfel/src/integration.c2
-rw-r--r--libcrystfel/src/mosflm.c22
-rw-r--r--libcrystfel/src/peakfinder8.c14
-rw-r--r--libcrystfel/src/predict-refine.c8
-rw-r--r--libcrystfel/src/reflist.c92
-rw-r--r--libcrystfel/src/reflist.h13
-rw-r--r--libcrystfel/src/stream.c86
-rw-r--r--libcrystfel/src/stream.h8
-rw-r--r--libcrystfel/src/taketwo.c1116
-rw-r--r--libcrystfel/src/xds.c273
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;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-28 09:23:49 +0000