diff options
| -rw-r--r-- | src/hrs-scaling.c | 153 | ||||
| -rw-r--r-- | src/scaling-report.c | 33 |
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); + } } |