Commit a4305d69 authored by Thomas White's avatar Thomas White Committed by Thomas White
Browse files

Move post refinement stuff to a new file

parent ed5e49e9
......@@ -64,7 +64,7 @@ calibrate_detector_LDADD = @LIBS@
facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
image.c geometry.c reflections.c stream.c thread-pool.c \
beam-parameters.c symmetry.c
beam-parameters.c symmetry.c post-refinement.c
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
......@@ -87,4 +87,4 @@ EXTRA_DIST = cell.h hdf5-file.h image.h utils.h diffraction.h detector.h \
render.h hdfsee.h dirax.h peaks.h index.h filters.h \
diffraction-gpu.h cl-utils.h symmetry.h \
povray.h index-priv.h geometry.h templates.h render_hkl.h \
stream.h thread-pool.h beam-parameters.h
stream.h thread-pool.h beam-parameters.h post-refinement.h
......@@ -84,7 +84,8 @@ am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \
hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \
peaks.$(OBJEXT) image.$(OBJEXT) geometry.$(OBJEXT) \
reflections.$(OBJEXT) stream.$(OBJEXT) thread-pool.$(OBJEXT) \
beam-parameters.$(OBJEXT) symmetry.$(OBJEXT)
beam-parameters.$(OBJEXT) symmetry.$(OBJEXT) \
post-refinement.$(OBJEXT)
facetron_OBJECTS = $(am_facetron_OBJECTS)
facetron_DEPENDENCIES =
am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
......@@ -324,7 +325,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
calibrate_detector_LDADD = @LIBS@
facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
image.c geometry.c reflections.c stream.c thread-pool.c \
beam-parameters.c symmetry.c
beam-parameters.c symmetry.c post-refinement.c
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
......@@ -345,7 +346,7 @@ EXTRA_DIST = cell.h hdf5-file.h image.h utils.h diffraction.h detector.h \
render.h hdfsee.h dirax.h peaks.h index.h filters.h \
diffraction-gpu.h cl-utils.h symmetry.h \
povray.h index-priv.h geometry.h templates.h render_hkl.h \
stream.h thread-pool.h beam-parameters.h
stream.h thread-pool.h beam-parameters.h post-refinement.h
all: all-am
......@@ -491,6 +492,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/indexamajig.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pattern_sim.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/peaks.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/post-refinement.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/povray.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/powder_plot.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/process_hkl.Po@am__quote@
......
......@@ -22,9 +22,6 @@
#include <getopt.h>
#include <assert.h>
#include <pthread.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include "utils.h"
#include "hdf5-file.h"
......@@ -35,18 +32,12 @@
#include "peaks.h"
#include "thread-pool.h"
#include "beam-parameters.h"
#include "post-refinement.h"
/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
#define MAX_CYCLES (100)
/* Refineable parameters */
enum {
REF_SCALE,
REF_DIV,
NUM_PARAMS
};
static void show_help(const char *s)
{
......@@ -83,190 +74,6 @@ struct refine_args
};
/* Return the gradient of parameter 'k' given the current status of 'image'. */
static double gradient(struct image *image, int k,
struct cpeak spot, double I_partial)
{
double ds;
double nom, den;
ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
switch ( k ) {
case REF_SCALE :
return I_partial;
case REF_DIV :
nom = sqrt(2.0) * ds * sin(image->div);
den = sqrt(1.0 - cos(image->div));
return nom/den;
}
ERROR("No gradient defined for parameter %i\n", k);
abort();
}
/* Apply the given shift to the 'k'th parameter of 'image'. */
static void apply_shift(struct image *image, int k, double shift)
{
switch ( k ) {
case REF_SCALE :
image->osf += shift;
break;
case REF_DIV :
STATUS("Shifting div by %e\n", shift);
image->div += shift;
break;
default :
ERROR("No shift defined for parameter %i\n", k);
abort();
}
}
static double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
const char *sym, double *i_full, FILE *graph)
{
int h;
double delta_I = 0.0;
for ( h=0; h<n; h++ ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full;
float I_partial;
float xc, yc;
hind = spots[h].h;
kind = spots[h].k;
lind = spots[h].l;
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( spots[h].p < 0.1 ) continue;
/* Actual measurement of this reflection from this
* pattern? */
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
continue;
}
I_partial *= image->osf;
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
delta_I += fabs(I_partial - spots[h].p * I_full);
if ( graph != NULL ) {
fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
" %5.2f %5.2f\n",
hind, kind, lind, I_partial/spots[h].p, xc, yc,
spots[h].p, I_partial / I_full);
}
}
return delta_I / (double)n;
}
static double iterate(struct image *image, double *i_full, const char *sym,
struct cpeak **pspots, int *n)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
int h, param;
struct cpeak *spots = *pspots;
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
/* Construct the equations, one per reflection in this image */
for ( h=0; h<*n; h++ ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full, delta_I;
float I_partial;
float xc, yc;
int k;
hind = spots[h].h;
kind = spots[h].k;
lind = spots[h].l;
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( spots[h].p < 0.1 ) continue;
/* Actual measurement of this reflection from this
* pattern? */
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
continue;
}
I_partial *= image->osf;
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
delta_I = I_partial - spots[h].p * I_full;
for ( k=0; k<NUM_PARAMS; k++ ) {
int g;
double v_c;
for ( g=0; g<NUM_PARAMS; g++ ) {
double M_curr, M_c;
M_curr = gsl_matrix_get(M, g, k);
M_c = gradient(image, g, spots[h], I_partial)
* gradient(image, k, spots[h], I_partial);
M_c *= pow(I_full, 2.0);
gsl_matrix_set(M, g, k, M_curr + M_c);
}
v_c = delta_I * I_full * gradient(image, k, spots[h],
I_partial);
gsl_vector_set(v, k, v_c);
}
}
shifts = gsl_vector_alloc(NUM_PARAMS);
gsl_linalg_HH_solve(M, v, shifts);
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
}
gsl_matrix_free(M);
gsl_vector_free(v);
gsl_vector_free(shifts);
free(spots);
spots = find_intersections(image, image->indexed_cell, n, 0);
*pspots = spots;
return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
}
static void refine_image(int mytask, void *tasks)
{
struct refine_args *all_args = tasks;
......@@ -299,7 +106,7 @@ static void refine_image(int mytask, void *tasks)
i = 0;
do {
last_dev = dev;
dev = iterate(image, pargs->i_full, pargs->sym, &spots, &n);
dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n);
STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
i++;
} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
......
/*
* post-refinement.c
*
* Post refinement
*
* (c) 2006-2010 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdlib.h>
#include <gsl/gsl_poly.h>
#include <assert.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include "image.h"
#include "post-refinement.h"
#include "peaks.h"
#include "symmetry.h"
#include "geometry.h"
/* Return the gradient of parameter 'k' given the current status of 'image'. */
double gradient(struct image *image, int k,
struct cpeak spot, double I_partial)
{
double ds;
double nom, den;
ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
switch ( k ) {
case REF_SCALE :
return I_partial;
case REF_DIV :
nom = sqrt(2.0) * ds * sin(image->div);
den = sqrt(1.0 - cos(image->div));
return nom/den;
}
ERROR("No gradient defined for parameter %i\n", k);
abort();
}
/* Apply the given shift to the 'k'th parameter of 'image'. */
void apply_shift(struct image *image, int k, double shift)
{
switch ( k ) {
case REF_SCALE :
image->osf += shift;
break;
case REF_DIV :
STATUS("Shifting div by %e\n", shift);
image->div += shift;
break;
default :
ERROR("No shift defined for parameter %i\n", k);
abort();
}
}
double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
const char *sym, double *i_full, FILE *graph)
{
int h;
double delta_I = 0.0;
for ( h=0; h<n; h++ ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full;
float I_partial;
float xc, yc;
hind = spots[h].h;
kind = spots[h].k;
lind = spots[h].l;
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( spots[h].p < 0.1 ) continue;
/* Actual measurement of this reflection from this
* pattern? */
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
continue;
}
I_partial *= image->osf;
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
delta_I += fabs(I_partial - spots[h].p * I_full);
if ( graph != NULL ) {
fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
" %5.2f %5.2f\n",
hind, kind, lind, I_partial/spots[h].p, xc, yc,
spots[h].p, I_partial / I_full);
}
}
return delta_I / (double)n;
}
/* Perform one cycle of post refinement on 'image' against 'i_full' */
double pr_iterate(struct image *image, double *i_full, const char *sym,
struct cpeak **pspots, int *n)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
int h, param;
struct cpeak *spots = *pspots;
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
/* Construct the equations, one per reflection in this image */
for ( h=0; h<*n; h++ ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full, delta_I;
float I_partial;
float xc, yc;
int k;
hind = spots[h].h;
kind = spots[h].k;
lind = spots[h].l;
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( spots[h].p < 0.1 ) continue;
/* Actual measurement of this reflection from this
* pattern? */
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
continue;
}
I_partial *= image->osf;
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
delta_I = I_partial - spots[h].p * I_full;
for ( k=0; k<NUM_PARAMS; k++ ) {
int g;
double v_c;
for ( g=0; g<NUM_PARAMS; g++ ) {
double M_curr, M_c;
M_curr = gsl_matrix_get(M, g, k);
M_c = gradient(image, g, spots[h], I_partial)
* gradient(image, k, spots[h], I_partial);
M_c *= pow(I_full, 2.0);
gsl_matrix_set(M, g, k, M_curr + M_c);
}
v_c = delta_I * I_full * gradient(image, k, spots[h],
I_partial);
gsl_vector_set(v, k, v_c);
}
}
shifts = gsl_vector_alloc(NUM_PARAMS);
gsl_linalg_HH_solve(M, v, shifts);
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
}
gsl_matrix_free(M);
gsl_vector_free(v);
gsl_vector_free(shifts);
free(spots);
spots = find_intersections(image, image->indexed_cell, n, 0);
*pspots = spots;
return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
}
/*
* post-refinement.h
*
* Post refinement
*
* (c) 2006-2010 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
*/
#ifndef POST_REFINEMENT_H
#define POST_REFINEMENT_H
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdio.h>
#include "image.h"
/* Refineable parameters */
enum {
REF_SCALE,
REF_DIV,
NUM_PARAMS
};
/* Return the gradient of parameter 'k' given the current status of 'image'. */
double gradient(struct image *image, int k, struct cpeak spot,
double I_partial);
/* Apply the given shift to the 'k'th parameter of 'image'. */
void apply_shift(struct image *image, int k, double shift);
double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
const char *sym, double *i_full, FILE *graph);
double pr_iterate(struct image *image, double *i_full, const char *sym,
struct cpeak **pspots, int *n);
#endif /* POST_REFINEMENT_H */
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