aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-11-20 22:42:57 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:06 +0100
commita4305d698c45fb68cd8e9ce2727b4cfdbda90f9e (patch)
treef4d763ba7fc8d2632a77b7b3e26fc065e90755a6
parented5e49e996cd4ee9d1118e19a27ee1ac7cc67c41 (diff)
Move post refinement stuff to a new file
-rw-r--r--src/Makefile.am4
-rw-r--r--src/Makefile.in8
-rw-r--r--src/facetron.c197
-rw-r--r--src/post-refinement.c213
-rw-r--r--src/post-refinement.h50
5 files changed, 272 insertions, 200 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index dbbdc9a2..5599da98 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -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
diff --git a/src/Makefile.in b/src/Makefile.in
index 00f9af02..fbdcd53c 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -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@
diff --git a/src/facetron.c b/src/facetron.c
index cf4f099e..9ea7bcea 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -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) );
diff --git a/src/post-refinement.c b/src/post-refinement.c
new file mode 100644
index 00000000..5ca0dc98
--- /dev/null
+++ b/src/post-refinement.c
@@ -0,0 +1,213 @@
+/*
+ * 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);
+}
diff --git a/src/post-refinement.h b/src/post-refinement.h
new file mode 100644
index 00000000..1ef7a1fd
--- /dev/null
+++ b/src/post-refinement.h
@@ -0,0 +1,50 @@
+/*
+ * 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 */