Skip to content
GitLab
Commit 31851dfe authored by Thomas White's avatar Thomas White
Browse files

Get rid of "image.raw_reflections"

parent 0b87230e
Hide whitespace changes
Inline Side-by-side
src/geometry.c
View file @ 31851dfe
......@@ -24,6 +24,7 @@
#include "beam-parameters.h"
#include "reflist.h"
#include "reflist-utils.h"
#include "symmetry.h"
static signed int locate_peak(double x, double y, double z, double k,
......@@ -293,10 +294,8 @@ double integrate_all(struct image *image, RefList *reflections)
/* Decide which reflections can be scaled */
static int select_scalable_reflections(RefList *list)
static void select_scalable_reflections(RefList *list, ReflItemList *sc_l)
{
int n_scalable = 0;
Reflection *refl;
RefListIterator *iter;
......@@ -312,16 +311,27 @@ static int select_scalable_reflections(RefList *list)
if ( v < 0.1 ) scalable = 0;
set_scalable(refl, scalable);
if ( scalable ) n_scalable++;
if ( scalable && (sc_l != NULL) ) {
}
signed int h, k, l;
get_indices(refl, &h, &k, &l); /* Should already be
* asymmetric */
if ( !find_item(sc_l, h, k, l) ) {
add_item(sc_l, h, k, l);
}
}
return n_scalable;
}
}
/* Calculate partialities and apply them to the image's raw_reflections */
/* Calculate partialities and apply them to the image's raw_reflections,
* returning a ReflItemList of the currentl scalable (asymmetric) reflections.
*/
void update_partialities(struct image *image, const char *sym,
ReflItemList *scalable,
int *n_expected, int *n_found, int *n_notfound)
{
Reflection *refl;
......@@ -335,56 +345,54 @@ void update_partialities(struct image *image, const char *sym,
refl = next_refl(refl, iter) )
{
Reflection *peak_in_pattern;
Reflection *p_peak;
double r1, r2, p, x, y;
signed int h, k, l;
signed int ha, ka, la;
int clamp1, clamp2;
int found = 0;
/* Get predicted indices and location */
get_indices(refl, &h, &k, &l);
get_detector_pos(refl, &x, &y);
if ( n_expected != NULL ) (*n_expected)++;
/* Get the asymmetric indices, with which the reflection in the
* image's list will be indexed. */
get_asymm(h, k, l, &ha, &ka, &la, sym);
/* Look for this reflection in the pattern */
peak_in_pattern = find_refl(image->raw_reflections, h, k, l);
if ( peak_in_pattern == NULL ) {
if ( n_notfound != NULL ) (*n_notfound)++;
p_peak = find_refl(image->reflections, ha, ka, la);
do {
signed int hs, ks, ls;
if ( p_peak != NULL ) {
get_symmetric_indices(p_peak, &hs, &ks, &ls);
if ( (hs==h) && (ks==k) && (ls==l) ) found = 1;
}
if ( !found && (p_peak != NULL ) ) {
p_peak = next_found_refl(p_peak);
}
} while ( !found && (p_peak != NULL) );
if ( !found ) {
if (n_notfound != NULL) (*n_notfound)++;
continue;
}
if ( n_found != NULL ) (*n_found)++;
/* Transfer partiality stuff */
get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
set_partial(peak_in_pattern, r1, r2, p, clamp1, clamp2);
set_partial(p_peak, r1, r2, p, clamp1, clamp2);
/* Transfer detector location */
get_detector_pos(refl, &x, &y);
set_detector_pos(peak_in_pattern, 0.0, x, y);
set_detector_pos(p_peak, 0.0, x, y);
}
reflist_free(predicted);
}
void update_partialities_and_asymm(struct image *image, const char *sym,
ReflItemList *obs,
int *n_expected, int *n_found,
int *n_notfound)
{
/* Get rid of the old list, about to be replaced */
reflist_free(image->reflections);
image->reflections = NULL;
/* Fill in partialities */
update_partialities(image, sym, n_expected, n_found, n_notfound);
/* Rewrite the reflections with the asymmetric indices
* to get the list used for scaling and post refinement */
image->reflections = asymmetric_indices(image->raw_reflections,
sym, obs);
select_scalable_reflections(image->reflections);
/* Need these lists to work fast */
optimise_reflist(image->reflections);
select_scalable_reflections(image->reflections, scalable);
}
src/geometry.h
View file @ 31851dfe
......@@ -24,8 +24,8 @@ extern RefList *find_intersections(struct image *image, UnitCell *cell,
extern double integrate_all(struct image *image, RefList *reflections);
extern void update_partialities_and_asymm(struct image *image, const char *sym,
ReflItemList *obs,
int *n_expected, int *n_found,
int *n_notfound);
extern void update_partialities(struct image *image, const char *sym,
ReflItemList *scalable,
int *n_expected, int *n_found, int *n_notfound);
#endif /* GEOMETRY_H */
src/image.h
View file @ 31851dfe
......@@ -87,7 +87,6 @@ typedef struct _imagefeaturelist ImageFeatureList;
* int height;
*
* RefList *reflections;
* RefList *raw_reflections;
*
* ImageFeatureList *features;
* };
......@@ -151,9 +150,6 @@ struct image {
/* Integrated (or about-to-be-integrated) reflections */
RefList *reflections;
/* Raw version of "reflections", e.g. in point group 1 */
RefList *raw_reflections;
/* Detected peaks */
ImageFeatureList *features;
......
src/partialator.c
View file @ 31851dfe
......@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
FILE *fh;
int nthreads = 1;
struct detector *det;
ReflItemList *obs;
ReflItemList *scalable;
int i;
int n_total_patterns;
struct image *images;
......@@ -282,9 +282,11 @@ int main(int argc, char *argv[])
/* Fill in what we know about the images so far */
rewind(fh);
obs = new_items();
scalable = new_items();
for ( i=0; i<n_total_patterns; i++ ) {
RefList *as;
images[n_usable_patterns].det = NULL;
if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
......@@ -316,13 +318,14 @@ int main(int argc, char *argv[])
images[n_usable_patterns].beam = NULL;
/* This is the raw list of reflections */
images[n_usable_patterns].raw_reflections =
images[n_usable_patterns].reflections;
images[n_usable_patterns].reflections = NULL;
as = asymmetric_indices(images[n_usable_patterns].reflections,
sym);
optimise_reflist(as);
reflist_free(images[n_usable_patterns].reflections);
images[n_usable_patterns].reflections = as;
update_partialities_and_asymm(&images[n_usable_patterns], sym,
obs, &n_expected, &n_found,
&n_notfound);
update_partialities(&images[n_usable_patterns], sym, scalable,
&n_expected, &n_found, &n_notfound);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
n_usable_patterns++;
......@@ -331,14 +334,15 @@ int main(int argc, char *argv[])
fclose(fh);
STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
100.0 * (double)n_found / n_expected, n_notfound, n_expected);
STATUS("Mean measurements per unique reflection: %5.2f\n",
(double)n_found / num_items(obs));
STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
(double)n_found / num_items(scalable));
cref = find_common_reflections(images, n_usable_patterns);
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
full = scale_intensities(images, n_usable_patterns, sym, obs, cref);
full = scale_intensities(images, n_usable_patterns, sym,
scalable, cref);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
......@@ -364,13 +368,13 @@ int main(int argc, char *argv[])
}
/* Refine the geometry of all patterns to get the best fit */
refine_all(images, n_total_patterns, det, sym, obs, full,
refine_all(images, n_total_patterns, det, sym, scalable, full,
nthreads, fhg, fhp);
/* Re-estimate all the full intensities */
reflist_free(full);
full = scale_intensities(images, n_usable_patterns,
sym, obs, cref);
sym, scalable, cref);
fclose(fhg);
fclose(fhp);
......@@ -388,10 +392,9 @@ int main(int argc, char *argv[])
/* Clean up */
for ( i=0; i<n_usable_patterns; i++ ) {
reflist_free(images[i].reflections);
reflist_free(images[i].raw_reflections);
}
reflist_free(full);
delete_items(obs);
delete_items(scalable);
free(sym);
free(outfile);
free_detector_geometry(det);
......
src/post-refinement.c
View file @ 31851dfe
......@@ -95,21 +95,21 @@ static double gradient(struct image *image, int k, Reflection *refl, double r)
double bsx, bsy, bsz;
double csx, csy, csz;
double xl, yl, zl;
signed int hi, ki, li;
signed int hs, ks, ls;
double r1, r2, p;
int clamp_low, clamp_high;
get_indices(refl, &hi, &ki, &li);
get_symmetric_indices(refl, &hs, &ks, &ls);
cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
xl = hi*asx + ki*bsx + li*csx;
yl = hi*asy + ki*bsy + li*csy;
zl = hi*asz + ki*bsz + li*csz;
xl = hs*asx + ks*bsx + ls*csx;
yl = hs*asy + ks*bsy + ls*csy;
zl = hs*asz + ks*bsz + ls*csz;
ds = 2.0 * resolution(image->indexed_cell, hi, ki, li);
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+k);
ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls);
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/image->lambda);
azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
......@@ -142,23 +142,23 @@ static double gradient(struct image *image, int k, Reflection *refl, double r)
/* Cell parameters and orientation */
case REF_ASX :
return hi * sin(tt) * cos(azi) * g;
return hs * sin(tt) * cos(azi) * g;
case REF_BSX :
return ki * sin(tt) * g;
return ks * sin(tt) * g;
case REF_CSX :
return li * sin(tt) * g;
return ls * sin(tt) * g;
case REF_ASY :
return hi * sin(tt) * g;
return hs * sin(tt) * g;
case REF_BSY :
return ki * sin(tt) * g;
return ks * sin(tt) * g;
case REF_CSY :
return li * sin(tt) * g;
return ls * sin(tt) * g;
case REF_ASZ :
return hi * cos(tt) * g;
return hs * cos(tt) * g;
case REF_BSZ :
return ki * cos(tt) * g;
return ks * cos(tt) * g;
case REF_CSZ :
return li * cos(tt) * g;
return ls * cos(tt) * g;
}
......@@ -285,8 +285,7 @@ static double pr_iterate(struct image *image, const RefList *full,
/* Calculate all gradients for this reflection */
for ( k=0; k<NUM_PARAMS; k++ ) {
double gr;
gr = gradient(image, k, refl,
image->profile_radius);
gr = gradient(image, k, refl, image->profile_radius);
gradients[k] = gr;
}
......@@ -403,8 +402,7 @@ static void plot_curve(struct image *image, const RefList *full,
ax = origval + (double)i*shval;
cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
update_partialities_and_asymm(image, sym,
NULL, NULL, NULL, NULL);
update_partialities(image, sym, NULL, NULL, NULL, NULL);
dev = mean_partial_dev(image, full, sym);
STATUS("%i %e %e\n", i, ax, dev);
......@@ -431,8 +429,7 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
double dev;
max_shift = pr_iterate(image, full, sym);
update_partialities_and_asymm(image, sym,
NULL, NULL, NULL, NULL);
update_partialities(image, sym, NULL, NULL, NULL, NULL);
dev = mean_partial_dev(image, full, sym);
STATUS("PR Iteration %2i: max shift = %5.2f dev = %5.2f\n",
......
src/reflist-utils.c
View file @ 31851dfe
......@@ -362,7 +362,7 @@ RefList *read_reflections(const char *filename)
}
RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs)
RefList *asymmetric_indices(RefList *in, const char *sym)
{
Reflection *refl;
RefListIterator *iter;
......@@ -386,12 +386,7 @@ RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs)
assert(cr != NULL);
copy_data(cr, refl);
if ( obs != NULL ) {
if ( !find_item(obs, ha, ka, la) ) {
add_item(obs, ha, ka, la);
}
}
set_symmetric_indices(cr, h, k, l);
}
......
src/reflist-utils.h
View file @ 31851dfe
......@@ -41,7 +41,6 @@ extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
signed int l, const char *sym, signed int *hu,
signed int *ku, signed int *lu);
extern RefList *asymmetric_indices(RefList *in, const char *sym,
ReflItemList *obs);
extern RefList *asymmetric_indices(RefList *in, const char *sym);
#endif /* REFLIST_UTILS_H */
src/reflist.c
View file @ 31851dfe
......@@ -45,6 +45,11 @@
struct _refldata {
/* Symmetric indices (i.e. the "real" indices) */
signed int hs;
signed int ks;
signed int ls;
/* Partiality and related geometrical stuff */
double r1; /* First excitation error */
double r2; /* Second excitation error */
......@@ -287,6 +292,29 @@ void get_indices(const Reflection *refl,
}
/**
* get_symmetric_indices:
* @refl: A %Reflection
* @h: Location at which to store the 'h' index of the reflection
* @k: Location at which to store the 'k' index of the reflection
* @l: Location at which to store the 'l' index of the reflection
*
* This function gives the symmetric indices, that is, the "real" indices before
* squashing down to the asymmetric reciprocal unit. This may be useful if the
* list is indexed according to the asymmetric indices, but you still need
* access to the symmetric version. This happens during post-refinement.
*
**/
void get_symmetric_indices(const Reflection *refl,
signed int *hs, signed int *ks,
signed int *ls)
{
*hs = refl->data.hs;
*ks = refl->data.ks;
*ls = refl->data.ls;
}
/**
* get_partiality:
* @refl: A %Reflection
......@@ -489,6 +517,15 @@ void set_ph(Reflection *refl, double phase)
}
void set_symmetric_indices(Reflection *refl,
signed int hs, signed int ks, signed int ls)
{
refl->data.hs = hs;
refl->data.ks = ks;
refl->data.ls = ls;
}
/********************************* Insertion **********************************/
static void insert_node(Reflection *head, Reflection *new)
......
src/reflist.h
View file @ 31851dfe
......@@ -47,6 +47,9 @@ extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
extern double get_partiality(const Reflection *refl);
extern void get_indices(const Reflection *refl,
signed int *h, signed int *k, signed int *l);
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 *r1, double *r2,
double *p, int *clamp_low, int *clamp_high);
......@@ -68,6 +71,8 @@ extern void set_redundancy(Reflection *refl, int red);
extern void set_sum_squared_dev(Reflection *refl, double dev);
extern void set_esd_intensity(Reflection *refl, double esd);
extern void set_ph(Reflection *refl, double phase);
extern void set_symmetric_indices(Reflection *refl,
signed int hs, signed int ks, signed int ls);
/* Insertion */
extern Reflection *add_refl(RefList *list,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use