aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore1
-rw-r--r--Makefile.am4
-rw-r--r--Makefile.in40
-rw-r--r--src/calibrate_detector.c460
4 files changed, 13 insertions, 492 deletions
diff --git a/.gitignore b/.gitignore
index 8cdbd120..ce0c0350 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,7 +30,6 @@ src/indexamajig
src/compare_hkl
src/powder_plot
src/render_hkl
-src/calibrate_detector
src/partialator
src/check_hkl
src/sum_stack
diff --git a/Makefile.am b/Makefile.am
index 48dd7633..857a7962 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -3,7 +3,7 @@ SUBDIRS = lib doc/reference libcrystfel
ACLOCAL_AMFLAGS = -I m4
bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
- src/compare_hkl src/powder_plot src/calibrate_detector \
+ src/compare_hkl src/powder_plot \
src/partialator src/check_hkl src/partial_sim
noinst_PROGRAMS = tests/list_check tests/integration_check \
@@ -65,8 +65,6 @@ if HAVE_CAIRO
src_render_hkl_SOURCES = src/render_hkl.c
endif
-src_calibrate_detector_SOURCES = src/calibrate_detector.c
-
src_partialator_SOURCES = src/partialator.c src/post-refinement.c \
src/hrs-scaling.c
diff --git a/Makefile.in b/Makefile.in
index 7b2f6f84..0e2a6d2b 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -38,9 +38,8 @@ host_triplet = @host@
bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
src/get_hkl$(EXEEXT) src/indexamajig$(EXEEXT) \
src/compare_hkl$(EXEEXT) src/powder_plot$(EXEEXT) \
- src/calibrate_detector$(EXEEXT) src/partialator$(EXEEXT) \
- src/check_hkl$(EXEEXT) src/partial_sim$(EXEEXT) \
- $(am__EXEEXT_1) $(am__EXEEXT_2)
+ src/partialator$(EXEEXT) src/check_hkl$(EXEEXT) \
+ src/partial_sim$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
noinst_PROGRAMS = tests/list_check$(EXEEXT) \
tests/integration_check$(EXEEXT) \
tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
@@ -93,19 +92,14 @@ am__installdirs = "$(DESTDIR)$(bindir)" "$(DESTDIR)$(man1dir)" \
@HAVE_OPENCL_TRUE@am__EXEEXT_3 = tests/gpu_sim_check$(EXEEXT)
PROGRAMS = $(bin_PROGRAMS) $(noinst_PROGRAMS)
am__dirstamp = $(am__leading_dot)dirstamp
-am_src_calibrate_detector_OBJECTS = src/calibrate_detector.$(OBJEXT)
-src_calibrate_detector_OBJECTS = $(am_src_calibrate_detector_OBJECTS)
-src_calibrate_detector_LDADD = $(LDADD)
-src_calibrate_detector_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \
- $(top_builddir)/libcrystfel/libcrystfel.la
-AM_V_lt = $(am__v_lt_$(V))
-am__v_lt_ = $(am__v_lt_$(AM_DEFAULT_VERBOSITY))
-am__v_lt_0 = --silent
am_src_check_hkl_OBJECTS = src/check_hkl.$(OBJEXT)
src_check_hkl_OBJECTS = $(am_src_check_hkl_OBJECTS)
src_check_hkl_LDADD = $(LDADD)
src_check_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \
$(top_builddir)/libcrystfel/libcrystfel.la
+AM_V_lt = $(am__v_lt_$(V))
+am__v_lt_ = $(am__v_lt_$(AM_DEFAULT_VERBOSITY))
+am__v_lt_0 = --silent
am_src_compare_hkl_OBJECTS = src/compare_hkl.$(OBJEXT)
src_compare_hkl_OBJECTS = $(am_src_compare_hkl_OBJECTS)
src_compare_hkl_LDADD = $(LDADD)
@@ -233,18 +227,16 @@ am__v_CCLD_0 = @echo " CCLD " $@;
AM_V_GEN = $(am__v_GEN_$(V))
am__v_GEN_ = $(am__v_GEN_$(AM_DEFAULT_VERBOSITY))
am__v_GEN_0 = @echo " GEN " $@;
-SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
- $(src_compare_hkl_SOURCES) $(src_get_hkl_SOURCES) \
- $(src_hdfsee_SOURCES) $(src_indexamajig_SOURCES) \
- $(src_partial_sim_SOURCES) $(src_partialator_SOURCES) \
- $(src_pattern_sim_SOURCES) $(src_powder_plot_SOURCES) \
- $(src_process_hkl_SOURCES) $(src_render_hkl_SOURCES) \
- $(tests_gpu_sim_check_SOURCES) \
+SOURCES = $(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \
+ $(src_get_hkl_SOURCES) $(src_hdfsee_SOURCES) \
+ $(src_indexamajig_SOURCES) $(src_partial_sim_SOURCES) \
+ $(src_partialator_SOURCES) $(src_pattern_sim_SOURCES) \
+ $(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \
+ $(src_render_hkl_SOURCES) $(tests_gpu_sim_check_SOURCES) \
$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
$(tests_pr_gradient_check_SOURCES) \
$(tests_symmetry_check_SOURCES)
-DIST_SOURCES = $(src_calibrate_detector_SOURCES) \
- $(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \
+DIST_SOURCES = $(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \
$(src_get_hkl_SOURCES) $(am__src_hdfsee_SOURCES_DIST) \
$(src_indexamajig_SOURCES) $(src_partial_sim_SOURCES) \
$(am__src_partialator_SOURCES_DIST) \
@@ -602,7 +594,6 @@ src_compare_hkl_SOURCES = src/compare_hkl.c
src_check_hkl_SOURCES = src/check_hkl.c
src_powder_plot_SOURCES = src/powder_plot.c
@HAVE_CAIRO_TRUE@src_render_hkl_SOURCES = src/render_hkl.c
-src_calibrate_detector_SOURCES = src/calibrate_detector.c
src_partialator_SOURCES = src/partialator.c src/post-refinement.c \
src/hrs-scaling.c $(am__append_6)
tests_list_check_SOURCES = tests/list_check.c
@@ -752,11 +743,6 @@ src/$(am__dirstamp):
src/$(DEPDIR)/$(am__dirstamp):
@$(MKDIR_P) src/$(DEPDIR)
@: > src/$(DEPDIR)/$(am__dirstamp)
-src/calibrate_detector.$(OBJEXT): src/$(am__dirstamp) \
- src/$(DEPDIR)/$(am__dirstamp)
-src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibrate_detector_DEPENDENCIES) src/$(am__dirstamp)
- @rm -f src/calibrate_detector$(EXEEXT)
- $(AM_V_CCLD)$(LINK) $(src_calibrate_detector_OBJECTS) $(src_calibrate_detector_LDADD) $(LIBS)
src/check_hkl.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/check_hkl$(EXEEXT): $(src_check_hkl_OBJECTS) $(src_check_hkl_DEPENDENCIES) src/$(am__dirstamp)
@@ -862,7 +848,6 @@ tests/symmetry_check$(EXEEXT): $(tests_symmetry_check_OBJECTS) $(tests_symmetry_
mostlyclean-compile:
-rm -f *.$(OBJEXT)
- -rm -f src/calibrate_detector.$(OBJEXT)
-rm -f src/check_hkl.$(OBJEXT)
-rm -f src/cl-utils.$(OBJEXT)
-rm -f src/compare_hkl.$(OBJEXT)
@@ -891,7 +876,6 @@ mostlyclean-compile:
distclean-compile:
-rm -f *.tab.c
-@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/calibrate_detector.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/check_hkl.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/cl-utils.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/compare_hkl.Po@am__quote@
diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
deleted file mode 100644
index f7316da4..00000000
--- a/src/calibrate_detector.c
+++ /dev/null
@@ -1,460 +0,0 @@
-/*
- * calibrate_detector.c
- *
- * Attempt to refine detector geometry
- *
- * (c) 2011 Rick Kirian <rkirian@asu.edu>
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <getopt.h>
-#include <fenv.h>
-#include <math.h>
-#include <gsl/gsl_linalg.h>
-
-#include "utils.h"
-#include "image.h"
-#include "detector.h"
-#include "index.h"
-#include "hdf5-file.h"
-#include "stream.h"
-#include "peaks.h"
-
-
-static void show_help(const char *s)
-{
- printf("Syntax: %s [options] -i <file.h5>\n\n", s);
- printf(
-"Stream-based optimisation of detector geometry.\n"
-"\n"
-" -h, --help Display this help message.\n"
-" -g. --geometry=<file> Get detector geometry from file.\n"
-" -i, --input=<file> Input filename.\n"
-" -m, --method=<method> The calibration method. Choose from:\n"
-" xy Determine panel shifts in plane of detector\n"
-" -o, --output=<file> Name of output geometry file.\n"
-" -n, --npeaks=<number> Don't refine unless this many peaks are found\n"
-" in the whole stream.\n"
-"\n");
-}
-
-
-static int write_detector_geometry(const char *filename, struct detector *det)
-{
- struct panel *p;
- int pi;
- FILE *fh;
-
- if ( filename == NULL ) return 2;
- if ( det->n_panels < 1 ) return 3;
-
- fh = fopen(filename, "w");
- if ( fh == NULL ) return 1;
-
- for ( pi=0; pi<det->n_panels; pi++) {
-
- p = &(det->panels[pi]);
-
- if ( p == NULL ) return 4;
-
- fprintf(fh, "%s/min_fs = %d\n", p->name, p->min_fs);
- fprintf(fh, "%s/min_ss = %d\n", p->name, p->min_ss);
- fprintf(fh, "%s/max_fs = %d\n", p->name, p->max_fs);
- fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss);
- fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow);
- fprintf(fh, "%s/res = %g\n", p->name, p->res);
- fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep);
- fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from);
- fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy);
- fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy);
- fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx);
- fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny);
- fprintf(fh, "%s/no_index = %d\n", p->name, p->no_index);
-
- }
- fclose(fh);
-
- return 0;
-}
-
-
-static int calculate_projected_peak(struct panel *panel, struct rvec q,
- double kk, double *fs, double *ss)
-{
- if ( panel == NULL ) return 1;
-
- double xd, yd, cl;
- double plx, ply;
- const double den = kk + q.w;
-
- /* Camera length for this panel */
- cl = panel->clen;
-
- /* Coordinates of peak relative to central beam, in m */
- xd = cl * q.u / den;
- yd = cl * q.v / den;
-
- /* Convert to pixels */
- xd *= panel->res;
- yd *= panel->res;
-
- /* Convert to relative to the panel corner */
- plx = xd - panel->cnx;
- ply = yd - panel->cny;
-
- *fs = panel->xfs*plx + panel->yfs*ply;
- *ss = panel->xss*plx + panel->yss*ply;
-
- *fs += panel->min_fs;
- *ss += panel->min_ss;
-
- /* Now, is this on this panel? */
- if ( *fs < panel->min_fs ) return 3;
- if ( *fs > panel->max_fs ) return 3;
- if ( *ss < panel->min_ss ) return 3;
- if ( *ss > panel->max_ss ) return 3;
-
- return 0;
-}
-
-
-static struct rvec nearest_bragg(struct image *image, struct rvec q)
-{
- struct rvec g;
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- double hd, kd, ld;
- double h, k, l;
- int s;
- double U[9];
- double hvec[3];
-
- /* Miller indices of nearest Bragg reflection */
- cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
-
- hd = q.u * ax + q.v * ay + q.w * az;
- kd = q.u * bx + q.v * by + q.w * bz;
- ld = q.u * cx + q.v * cy + q.w * cz;
-
- h = lrint(hd);
- k = lrint(kd);
- l = lrint(ld);
-
- /* Now get scattering vector for reflection hkl
- * by solving the equation U*q = h */
- U[0] = ax; U[1] = ay; U[2] = az;
- U[3] = bx; U[4] = by; U[5] = bz;
- U[6] = cx; U[7] = cy; U[8] = cz;
-
- hvec[0] = h; hvec[1] = k; hvec[2] = l;
-
- gsl_matrix_view m = gsl_matrix_view_array(U, 3, 3);
- gsl_vector_view b = gsl_vector_view_array(hvec, 3);
- gsl_vector *x = gsl_vector_alloc(3);
-
- gsl_permutation *perm = gsl_permutation_alloc(3);
- gsl_linalg_LU_decomp(&m.matrix, perm, &s);
- gsl_linalg_LU_solve(&m.matrix, perm, &b.vector, x);
-
- /* Outgoing wavevector for hkl */
- g.u = x->data[0];
- g.v = x->data[1];
- g.w = x->data[2];
-
- return g;
-}
-
-
-static void refine_xy(FILE *fh, struct image *image, int minpeaks,
- const char *outfilename)
-{
- double *weightedSumFS;
- double *weightedSumSS;
- double *summedWeights;
- double *meanShiftFS;
- double *meanShiftSS;
- int *peaksFound;
- double cnx, cny;
- double xsh, ysh;
- int pi;
- int nChunks;
-
- weightedSumFS = calloc(image->det->n_panels, sizeof(double));
- weightedSumSS = calloc(image->det->n_panels, sizeof(double));
- summedWeights = calloc(image->det->n_panels, sizeof(double));
- peaksFound = calloc(image->det->n_panels, sizeof(double));
- meanShiftFS = calloc(image->det->n_panels, sizeof(double));
- meanShiftSS = calloc(image->det->n_panels, sizeof(double));
-
- /* Initialize arrays */
- for (pi=0; pi<image->det->n_panels; pi++) {
- weightedSumFS[pi] = 0;
- weightedSumSS[pi] = 0;
- summedWeights[pi] = 0;
- peaksFound[pi] = 0;
- meanShiftFS[pi] = 0;
- meanShiftSS[pi] = 0;
- }
-
- fesetround(1); /* Round towards nearest */
-
- nChunks = 0;
- while ( 1 ) {
-
- int fail;
- int nFeatures;
- int i;
-
- /* Get next chunk */
- fail = read_chunk(fh, image);
- if ( fail ) break;
- nChunks += 1;
-
- /* Skip if no peaks found */
- if ( image->features == NULL ) {
- continue;
- }
-
- /* Loop through peaks to determine mean panel shift */
- nFeatures = image_feature_count(image->features);
-
- for ( i=0; i<nFeatures; i++ ) {
-
- struct panel *p;
- struct imagefeature *thisFeature;
- struct rvec q;
- double twotheta;
- struct rvec g;
- double fs, ss; /* Observed peaks */
- double pfs, pss; /* Predicted peaks */
- double dfs, dss; /* Observed - predicted */
- int pi;
- double thisWeight;
-
- /* If we find a feature, determine peak
- * position */
- thisFeature = image_get_feature(image->features, i);
- if ( thisFeature == NULL ) {
- continue;
- }
-
- fs = thisFeature->fs;
- ss = thisFeature->ss;
-
- p = find_panel(image->det, fs, ss);
- if ( p == NULL ) {
- continue;
- }
- if ( p->no_index ) continue;
-
- /* Now determine the predicted peak position.
- * Scattering vector of this peak */
- q = get_q(image, fs, ss, &twotheta, 1.0/image->lambda);
-
- /* Scattering vector of nearest bragg peak */
- g = nearest_bragg(image, q);
-
- /* Coordinate of this predicted peak */
- fail = calculate_projected_peak(p, g, 1.0/image->lambda,
- &pfs, &pss);
- if ( fail != 0 ) continue;
-
- /* Finally, we have the shift for this peak */
- dfs = pfs - fs;
- dss = pss - ss;
-
- /* Add this shift to the weighted sum over shifts */
- pi = find_panel_number(image->det,fs,ss);
- thisWeight = 1.0;
- weightedSumFS[pi] += thisWeight*dfs;
- weightedSumSS[pi] += thisWeight*dss;
- summedWeights[pi] += thisWeight;
- peaksFound[pi] += 1;
-
- }
-
- }
-
- /* Calculate weighted average shift in peak positions */
- for ( pi=0; pi<image->det->n_panels; pi++ ) {
- meanShiftFS[pi] = weightedSumFS[pi]/summedWeights[pi];
- meanShiftSS[pi] = weightedSumSS[pi]/summedWeights[pi];
- }
-
- /* Populate the image structure with new geometry info */
- for ( pi=0; pi<image->det->n_panels; pi++ ) {
-
- struct panel *p;
-
- p = &image->det->panels[pi];
-
- xsh = 0;
- ysh = 0;
-
- if ( peaksFound[pi] >= minpeaks ) {
-
- /* Convert shifts from raw coords to lab frame */
- xsh = meanShiftFS[pi]*p->fsx + meanShiftSS[pi]*p->ssx;
- ysh = meanShiftFS[pi]*p->fsy + meanShiftSS[pi]*p->ssy;
-
- /* Add shifts to original panel corner locations */
- cnx = p->cnx + xsh;
- cny = p->cny + ysh;
-
- } else {
-
- /* Not refined? use original coordinates */
- cnx = p->cnx;
- cny = p->cny;
-
- }
-
- image->det->panels[pi].cnx = cnx;
- image->det->panels[pi].cny = cny;
- if ( peaksFound[pi] < minpeaks) {
- image->det->panels[pi].no_index = 1;
- }
-
- STATUS("Panel %s, # peaks: %10d, mean shifts: %f %f\n",
- p->name, peaksFound[pi], xsh, ysh);
-
- }
-
- /* Write the new geometry file */
- write_detector_geometry(outfilename, image->det);
-}
-
-
-int main(int argc, char *argv[])
-{
- char c;
- struct image image;
- char *filename = NULL;
- char *geometry = NULL;
- char *method = NULL;
- char *outputfile = NULL;
- FILE *fh = NULL;
- FILE *outfh = NULL;
- int minpeaks = 0;
-
- /* Long options */
- const struct option longopts[] = {
- {"help", 0, NULL, 'h'},
- {"input", 1, NULL, 'i'},
- {"geometry", 1, NULL, 'g'},
- {"method", 1, NULL, 'm'},
- {"output", 1, NULL, 'o'},
- {"npeaks", 1, NULL, 'n'},
- {0, 0, NULL, 0}
- };
-
- /* Short options */
- while ((c = getopt_long(argc, argv, "hi:g:m:o:n:",
- longopts, NULL)) != -1) {
-
- switch (c) {
- case 'h' :
- show_help(argv[0]);
- return 0;
- case 'i' :
- filename = strdup(optarg);
- break;
- case 'g' :
- geometry = strdup(optarg);
- break;
- case 'm' :
- method = strdup(optarg);
- break;
- case 'o' :
- outputfile = strdup(optarg);
- break;
- case 'n' :
- minpeaks = atoi(optarg);
- break;
- case 0 :
- break;
- default :
- return 1;
- }
-
- }
-
- if ( filename == NULL ) {
- ERROR("You must specify the input filename with -i\n");
- return 1;
- }
-
- fh = fopen(filename,"r");
- if ( fh == NULL ) {
- ERROR("Couldn't open file '%s'\n", filename);
- return 1;
- }
-
- if ( geometry == NULL ) {
- ERROR("You need to specify a geometry file with --geometry\n");
- return 1;
- }
-
- image.det = get_detector_geometry(geometry);
- if ( image.det == NULL ) {
- ERROR("Failed to read detector geometry from %s\n", geometry);
- return 1;
- }
- free(geometry);
-
- image.width = image.det->max_fs;
- image.height = image.det->max_ss;
-
- if ( outputfile != NULL ) {
- STATUS("Writing result to '%s'\n", outputfile);
- outfh = fopen(outputfile, "w");
- } else {
- ERROR("You need to specify an output file.\n");
- return 1;
- }
-
- if ( minpeaks == 0 ) {
- minpeaks = 1;
- STATUS("You did not specify minimum number of peaks.\n");
- STATUS("Using default value of %d\n", minpeaks);
- }
-
- if ( method == NULL ) {
- STATUS("You did not specify a refinement method "
- "- using default.\n");
- method = strdup("xy");
- }
-
- if ( strcmp(method,"xy") == 0 ) {
-
- STATUS("Using refinement method %s\n", method);
-
- refine_xy(fh, &image, minpeaks, outputfile);
-
- } else {
-
- printf("Refinement method %s not recognized\n",method);
- return 1;
-
- }
-
- fclose(outfh);
-
- return 0;
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-20 16:51:49 +0000