Calculate reflection ESDs, apply minimum redundancy and show ESD on histograms

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);
+	}
 }