aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/hrs-scaling.c153
-rw-r--r--src/scaling-report.c33
2 files changed, 163 insertions, 23 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 5349e577..b28aae3d 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -38,17 +38,16 @@
#define SCALING_RESTRAINT (1.0)
-struct queue_args
+struct scale_queue_args
{
RefList *reference;
struct image *images;
int n_started;
- int n_to_do;
double max_shift;
};
-struct worker_args
+struct scale_worker_args
{
struct image *image;
double shift;
@@ -56,12 +55,12 @@ struct worker_args
};
-static void *create_job(void *vqargs)
+static void *create_scale_job(void *vqargs)
{
- struct worker_args *wargs;
- struct queue_args *qargs = vqargs;
+ struct scale_worker_args *wargs;
+ struct scale_queue_args *qargs = vqargs;
- wargs = malloc(sizeof(struct worker_args));
+ wargs = malloc(sizeof(struct scale_worker_args));
wargs->reference = qargs->reference;
wargs->image = &qargs->images[qargs->n_started++];
@@ -70,9 +69,9 @@ static void *create_job(void *vqargs)
}
-static void run_job(void *vwargs, int cookie)
+static void run_scale_job(void *vwargs, int cookie)
{
- struct worker_args *wargs = vwargs;
+ struct scale_worker_args *wargs = vwargs;
struct image *image = wargs->image;
RefList *reference = wargs->reference;
Reflection *refl;
@@ -129,10 +128,10 @@ static void run_job(void *vwargs, int cookie)
}
-static void finalise_job(void *vqargs, void *vwargs)
+static void finalise_scale_job(void *vqargs, void *vwargs)
{
- struct queue_args *qargs = vqargs;
- struct worker_args *wargs = vwargs;
+ struct scale_queue_args *qargs = vqargs;
+ struct scale_worker_args *wargs = vwargs;
if ( wargs->shift > qargs->max_shift ) qargs->max_shift = wargs->shift;
free(wargs);
@@ -142,18 +141,17 @@ static void finalise_job(void *vqargs, void *vwargs)
static double iterate_scale(struct image *images, int n, RefList *reference,
int n_threads)
{
- struct queue_args qargs;
+ struct scale_queue_args qargs;
assert(reference != NULL);
qargs.reference = reference;
qargs.n_started = 0;
- qargs.n_to_do = n;
qargs.images = images;
qargs.max_shift = 0.0;
- run_threads(n_threads, run_job, create_job, finalise_job,
- &qargs, n, 0, 0, 0);
+ run_threads(n_threads, run_scale_job, create_scale_job,
+ finalise_scale_job, &qargs, n, 0, 0, 0);
return qargs.max_shift;
}
@@ -165,7 +163,6 @@ struct merge_queue_args
pthread_mutex_t full_lock;
struct image *images;
int n_started;
- int n_to_do;
};
@@ -268,7 +265,6 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
qargs.full = full;
qargs.n_started = 0;
- qargs.n_to_do = n;
qargs.images = images;
pthread_mutex_init(&qargs.full_lock, NULL);
@@ -290,6 +286,124 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
}
+
+struct esd_queue_args
+{
+ RefList *full;
+ struct image *images;
+ int n_started;
+};
+
+
+struct esd_worker_args
+{
+ struct image *image;
+ RefList *full;
+};
+
+
+static void *create_esd_job(void *vqargs)
+{
+ struct esd_worker_args *wargs;
+ struct esd_queue_args *qargs = vqargs;
+
+ wargs = malloc(sizeof(struct esd_worker_args));
+ wargs->full = qargs->full;
+
+ wargs->image = &qargs->images[qargs->n_started++];
+
+ return wargs;
+}
+
+
+static void run_esd_job(void *vwargs, int cookie)
+{
+ struct esd_worker_args *wargs = vwargs;
+ struct image *image = wargs->image;
+ RefList *full = wargs->full;
+ Reflection *refl;
+ RefListIterator *iter;
+ double G;
+
+ if ( image->pr_dud ) return;
+
+ G = image->osf;
+
+ for ( refl = first_refl(image->reflections, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ Reflection *f;
+ signed int h, k, l;
+ double num;
+ double Ihl, Ih;
+
+ if ( !get_scalable(refl) ) continue;
+
+ get_indices(refl, &h, &k, &l);
+ f = find_refl(full, h, k, l);
+ assert(f != NULL);
+
+ lock_reflection(f);
+
+ num = get_temp1(f);
+
+ Ih = get_intensity(f);
+ Ihl = get_intensity(refl) / (get_partiality(refl) * G);
+
+ num += pow(Ihl - Ih, 2.0);
+
+ set_temp1(f, num);
+ unlock_reflection(f);
+ }
+}
+
+
+static void finalise_esd_job(void *vqargs, void *vwargs)
+{
+ free(vwargs);
+}
+
+
+static void calculate_esds(struct image *images, int n, RefList *full,
+ int n_threads, int min_red)
+{
+ struct esd_queue_args qargs;
+ Reflection *refl;
+ RefListIterator *iter;
+
+ qargs.full = full;
+ qargs.n_started = 0;
+ qargs.images = images;
+
+ for ( refl = first_refl(full, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ set_temp1(refl, 0.0);
+ set_temp2(refl, 0.0);
+ }
+
+ run_threads(n_threads, run_esd_job, create_esd_job,
+ finalise_esd_job, &qargs, n, 0, 0, 0);
+
+ for ( refl = first_refl(full, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double esd;
+ int red = get_redundancy(refl);
+ esd = sqrt(get_temp1(refl));
+ esd /= (double)red;
+ set_esd_intensity(refl, esd);
+
+ if ( red < min_red ) {
+ set_redundancy(refl, 0);
+ }
+ }
+}
+
+
/* Scale the stack of images */
RefList *scale_intensities(struct image *images, int n, RefList *gref,
int n_threads)
@@ -297,6 +411,7 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
int i;
double max_corr;
RefList *full = NULL;
+ const int min_redundancy = 3;
/* No reference -> create an initial list to refine against */
if ( gref == NULL ) {
@@ -336,5 +451,7 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
full = lsq_intensities(images, n, n_threads);
} /* else we already did it */
+ calculate_esds(images, n, full, n_threads, min_redundancy);
+
return full;
}
diff --git a/src/scaling-report.c b/src/scaling-report.c
index 719402bc..1aadad2c 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -506,7 +506,7 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
cairo_new_path(cr);
cairo_rectangle(cr, 0.0, 0.0, g_width, g_height);
cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
- cairo_set_line_width(cr, 1.5);
+ cairo_set_line_width(cr, 1.0);
cairo_stroke(cr);
}
@@ -526,15 +526,17 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
const double g_height = 55.0;
char tmp[64];
Reflection *f;
- double Ifull, pos;
+ double Ifull, esd_Ifull, pos, mI, bit;
int have_full;
f = find_refl(full, h, k, l);
if ( f != NULL ) {
Ifull = get_intensity(f);
+ esd_Ifull = get_esd_intensity(f);
have_full = 1;
} else {
Ifull = 0.0;
+ esd_Ifull = 0.0;
have_full = 0;
}
@@ -642,14 +644,35 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
cairo_set_line_width(cr, 1.5);
cairo_stroke(cr);
- pos = Ifull/int_high[nbins-1];
- cairo_arc(cr, g_width*pos, g_height, g_width/30.0, 0.0, 2.0*M_PI);
+ mI = int_high[nbins-1];
+ pos = Ifull/mI;
+ bit = g_height / 15.0;
+ cairo_arc(cr, g_width*pos, g_height+bit, bit, 0.0, 2.0*M_PI);
if ( have_full ) {
- cairo_set_source_rgb(cr, 0.0, 0.69, 1.0);
+ cairo_set_source_rgb(cr, 0.0, 0.67, 0.45);
} else {
cairo_set_source_rgb(cr, 0.86, 0.0, 0.0);
}
cairo_fill(cr);
+
+ if ( have_full ) {
+
+ double eW = g_width*esd_Ifull/mI;
+
+ cairo_new_path(cr);
+ cairo_rectangle(cr, 0.0, g_height+bit*2.0,
+ g_width, g_height+bit*2.0);
+ //cairo_clip(cr);
+
+ cairo_new_path(cr);
+ cairo_move_to(cr, g_width*pos - eW, g_height+bit);
+ cairo_line_to(cr, g_width*pos + eW, g_height+bit);
+ cairo_set_source_rgb(cr, 0.0, 0.67, 0.45);
+ cairo_set_line_width(cr, 2.0);
+ cairo_stroke(cr);
+
+ cairo_reset_clip(cr);
+ }
}