diff options
-rw-r--r-- | Makefile.am | 6 | ||||
-rw-r--r-- | Makefile.in | 53 | ||||
-rw-r--r-- | src/index.c | 16 | ||||
-rw-r--r-- | src/index.h | 1 | ||||
-rw-r--r-- | src/templates.c | 326 | ||||
-rw-r--r-- | src/templates.h | 33 |
6 files changed, 27 insertions, 408 deletions
diff --git a/Makefile.am b/Makefile.am index 052b9e62..8719496a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -41,7 +41,7 @@ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \ src/cell.c src/image.c src/peaks.c src/index.c \ src/filters.c src/diffraction.c src/detector.c \ src/sfac.c src/dirax.c src/mosflm.c \ - src/reflections.c src/templates.c src/symmetry.c \ + src/reflections.c src/symmetry.c \ src/geometry.c src/thread-pool.c \ src/beam-parameters.c src/reflist.c if HAVE_OPENCL @@ -100,7 +100,7 @@ endif src_reintegrate_SOURCES = src/reintegrate.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ src/stream.c src/index.c src/dirax.c src/mosflm.c \ - src/templates.c src/geometry.c src/symmetry.c \ + src/geometry.c src/symmetry.c \ src/thread-pool.c src/reflist.c src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \ @@ -117,7 +117,7 @@ EXTRA_DIST += src/cell.h src/hdf5-file.h src/image.h src/utils.h \ src/render.h src/hdfsee.h src/dirax.h src/mosflm.h src/peaks.h \ src/index.h src/filters.h src/diffraction-gpu.h src/cl-utils.h \ src/symmetry.h src/povray.h src/index-priv.h src/geometry.h \ - src/templates.h src/render_hkl.h src/stream.h src/thread-pool.h \ + src/render_hkl.h src/stream.h src/thread-pool.h \ src/beam-parameters.h src/post-refinement.h src/hrs-scaling.h \ src/reflist.h diff --git a/Makefile.in b/Makefile.in index 3e8f1851..50e4837c 100644 --- a/Makefile.in +++ b/Makefile.in @@ -144,10 +144,9 @@ src_hdfsee_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am__src_indexamajig_SOURCES_DIST = src/indexamajig.c src/hdf5-file.c \ src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \ src/filters.c src/diffraction.c src/detector.c src/sfac.c \ - src/dirax.c src/mosflm.c src/reflections.c src/templates.c \ - src/symmetry.c src/geometry.c src/thread-pool.c \ - src/beam-parameters.c src/reflist.c src/diffraction-gpu.c \ - src/cl-utils.c + src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \ + src/geometry.c src/thread-pool.c src/beam-parameters.c \ + src/reflist.c src/diffraction-gpu.c src/cl-utils.c @HAVE_OPENCL_TRUE@am__objects_1 = src/diffraction-gpu.$(OBJEXT) \ @HAVE_OPENCL_TRUE@ src/cl-utils.$(OBJEXT) am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ @@ -156,10 +155,9 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ src/filters.$(OBJEXT) src/diffraction.$(OBJEXT) \ src/detector.$(OBJEXT) src/sfac.$(OBJEXT) src/dirax.$(OBJEXT) \ src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \ - src/templates.$(OBJEXT) src/symmetry.$(OBJEXT) \ - src/geometry.$(OBJEXT) src/thread-pool.$(OBJEXT) \ - src/beam-parameters.$(OBJEXT) src/reflist.$(OBJEXT) \ - $(am__objects_1) + src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \ + src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ + src/reflist.$(OBJEXT) $(am__objects_1) src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS) src_indexamajig_LDADD = $(LDADD) src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -208,9 +206,9 @@ am_src_reintegrate_OBJECTS = src/reintegrate.$(OBJEXT) \ src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) src/utils.$(OBJEXT) \ src/detector.$(OBJEXT) src/peaks.$(OBJEXT) src/image.$(OBJEXT) \ src/stream.$(OBJEXT) src/index.$(OBJEXT) src/dirax.$(OBJEXT) \ - src/mosflm.$(OBJEXT) src/templates.$(OBJEXT) \ - src/geometry.$(OBJEXT) src/symmetry.$(OBJEXT) \ - src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT) + src/mosflm.$(OBJEXT) src/geometry.$(OBJEXT) \ + src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT) \ + src/reflist.$(OBJEXT) src_reintegrate_OBJECTS = $(am_src_reintegrate_OBJECTS) src_reintegrate_LDADD = $(LDADD) src_reintegrate_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -555,17 +553,16 @@ EXTRA_DIST = configure m4/gnulib-cache.m4 src/cell.h src/hdf5-file.h \ src/dw-hdfsee.h src/render.h src/hdfsee.h src/dirax.h \ src/mosflm.h src/peaks.h src/index.h src/filters.h \ src/diffraction-gpu.h src/cl-utils.h src/symmetry.h \ - src/povray.h src/index-priv.h src/geometry.h src/templates.h \ - src/render_hkl.h src/stream.h src/thread-pool.h \ - src/beam-parameters.h src/post-refinement.h src/hrs-scaling.h \ - src/reflist.h data/sfac/Ca.nff data/sfac/C.nff \ - data/sfac/Fe.nff data/sfac/H.nff data/sfac/Mg.nff \ - data/sfac/N.nff data/sfac/O.nff data/sfac/P.nff \ - data/sfac/S.nff data/sfac/f0_WaasKirf.dat data/defs.h \ - data/diffraction.cl data/hdfsee.ui doc/geometry.txt \ - doc/indexamajig.txt doc/pattern_sim.txt doc/process_hkl.txt \ - doc/symmetry.txt doc/twin-calculator.pdf doc/0-INDEX \ - doc/examples/lcls-dec.geom \ + src/povray.h src/index-priv.h src/geometry.h src/render_hkl.h \ + src/stream.h src/thread-pool.h src/beam-parameters.h \ + src/post-refinement.h src/hrs-scaling.h src/reflist.h \ + data/sfac/Ca.nff data/sfac/C.nff data/sfac/Fe.nff \ + data/sfac/H.nff data/sfac/Mg.nff data/sfac/N.nff \ + data/sfac/O.nff data/sfac/P.nff data/sfac/S.nff \ + data/sfac/f0_WaasKirf.dat data/defs.h data/diffraction.cl \ + data/hdfsee.ui doc/geometry.txt doc/indexamajig.txt \ + doc/pattern_sim.txt doc/process_hkl.txt doc/symmetry.txt \ + doc/twin-calculator.pdf doc/0-INDEX doc/examples/lcls-dec.geom \ doc/examples/lcls-june-r0013-r0128.geom \ doc/examples/lcls-xpp-estimate.geom doc/examples/simple.geom \ doc/examples/lcls-dec.beam doc/examples/lcls-june.beam \ @@ -597,9 +594,9 @@ src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c \ src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \ src/filters.c src/diffraction.c src/detector.c src/sfac.c \ - src/dirax.c src/mosflm.c src/reflections.c src/templates.c \ - src/symmetry.c src/geometry.c src/thread-pool.c \ - src/beam-parameters.c src/reflist.c $(am__append_4) + src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \ + src/geometry.c src/thread-pool.c src/beam-parameters.c \ + src/reflist.c $(am__append_4) @HAVE_GTK_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \ @HAVE_GTK_TRUE@ src/hdf5-file.c src/utils.c src/image.c src/filters.c \ @HAVE_GTK_TRUE@ src/thread-pool.c src/detector.c @@ -647,7 +644,7 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src_reintegrate_SOURCES = src/reintegrate.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ src/stream.c src/index.c src/dirax.c src/mosflm.c \ - src/templates.c src/geometry.c src/symmetry.c \ + src/geometry.c src/symmetry.c \ src/thread-pool.c src/reflist.c src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \ @@ -860,8 +857,6 @@ src/diffraction.$(OBJEXT): src/$(am__dirstamp) \ src/dirax.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/mosflm.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) -src/templates.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) src/geometry.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/reflist.$(OBJEXT): src/$(am__dirstamp) \ @@ -968,7 +963,6 @@ mostlyclean-compile: -rm -f src/stream.$(OBJEXT) -rm -f src/sum_stack.$(OBJEXT) -rm -f src/symmetry.$(OBJEXT) - -rm -f src/templates.$(OBJEXT) -rm -f src/thread-pool.$(OBJEXT) -rm -f src/utils.$(OBJEXT) -rm -f tests/list_check.$(OBJEXT) @@ -1016,7 +1010,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/stream.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sum_stack.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/symmetry.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/templates.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/thread-pool.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/utils.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@ diff --git a/src/index.c b/src/index.c index 3d2246bc..d6375834 100644 --- a/src/index.c +++ b/src/index.c @@ -30,7 +30,6 @@ #include "detector.h" #include "index.h" #include "index-priv.h" -#include "templates.h" /* Base class constructor for unspecialised indexing private data */ @@ -67,10 +66,6 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, case INDEXING_MOSFLM : iprivs[n] = indexing_private(indm[n]); break; - case INDEXING_TEMPLATE : - iprivs[n] = generate_templates(cell, filename, det, - nominal_photon_energy); - break; } } @@ -98,8 +93,6 @@ void cleanup_indexing(IndexingPrivate **priv) case INDEXING_MOSFLM : free(priv[n]); break; - case INDEXING_TEMPLATE : - free_templates(priv[n]); } n++; @@ -153,18 +146,14 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, STATUS("Running MOSFLM...\n"); run_mosflm(image, cell); break; - case INDEXING_TEMPLATE : - match_templates(image, ipriv[n]); - break; } - if ( image->ncells == 0 ) { STATUS("No candidate cells found.\n"); n++; continue; } - if ( (cellr == CELLR_NONE) || (indm[n] == INDEXING_TEMPLATE) ) { + if ( cellr == CELLR_NONE ) { image->indexed_cell = cell_new_from_cell( image->candidate_cells[0]); if ( verbose ) { @@ -254,9 +243,6 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell) list[i] = INDEXING_DIRAX; } else if ( strcmp(methods[i], "mosflm") == 0) { list[i] = INDEXING_MOSFLM; - } else if ( strcmp(methods[i], "template") == 0) { - list[i] = INDEXING_TEMPLATE; - *need_cell = 1; } else { ERROR("Unrecognised indexing method '%s'\n", methods[i]); diff --git a/src/index.h b/src/index.h index 8f3e0ab7..005f0d93 100644 --- a/src/index.h +++ b/src/index.h @@ -28,7 +28,6 @@ typedef enum { INDEXING_NONE, INDEXING_DIRAX, INDEXING_MOSFLM, - INDEXING_TEMPLATE } IndexingMethod; diff --git a/src/templates.c b/src/templates.c deleted file mode 100644 index ed557afa..00000000 --- a/src/templates.c +++ /dev/null @@ -1,326 +0,0 @@ -/* - * templates.c - * - * Indexing by template matching - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - - -#include "index.h" -#include "index-priv.h" -#include "symmetry.h" -#include "utils.h" -#include "geometry.h" -#include "hdf5-file.h" -#include "peaks.h" -#include "reflist.h" - -#include <assert.h> - - -#define INTEGRATION_SQUARE_SIDE (10) -#define THRESHOLD (500) - - -/* Private data for template indexing */ -struct _indexingprivate_template -{ - struct _indexingprivate base; - UnitCell *cell; - int n_templates; - struct template *templates; -}; - - -struct template { - double omega; - double phi; - RefList *spots; -}; - - -/* Generate templates for the given cell using a representative image */ -IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, - struct detector *det, - double nominal_photon_energy) -{ - struct _indexingprivate_template *priv; - const char *sym; - double omega_max, phi_max; - int n_templates; - const double omega_step = deg2rad(0.5); - const double phi_step = deg2rad(0.5); - double omega, phi; - struct image image; - struct hdfile *hdfile; - int i; - - hdfile = hdfile_open(filename); - if ( hdfile == NULL ) { - return NULL; - } else if ( hdfile_set_image(hdfile, "/data/data0") ) { - ERROR("Couldn't select path\n"); - return NULL; - } - hdf5_read(hdfile, &image, 0); - hdfile_close(hdfile); - image.det = det; - if ( image.lambda < 0.0 ) image.lambda = nominal_photon_energy; - - priv = calloc(1, sizeof(struct _indexingprivate_template)); - priv->base.indm = INDEXING_TEMPLATE; - - sym = cell_get_pointgroup(cell); - - /* These define the orientation in space */ - if ( is_polyhedral(sym) ) { - ERROR("WARNING: Point group is polyhedral.\n"); - ERROR("This means I can't properly determine the orientation"); - ERROR(" ranges for template matching. Expect trouble.\n"); - } - omega_max = 2.0*M_PI / rotational_order(sym); - if ( has_bisecting_mirror_or_diad(sym) ) omega_max /= 2.0; - phi_max = M_PI; - if ( has_perpendicular_mirror(sym) ) phi_max /= 2.0; - - /* One more axis would define the rotation in the plane of the image */ - - STATUS("Orientation ranges in %s: %.0f-%.0f, %.0f-%.0f deg.\n", - sym, 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max)); - - n_templates = (omega_max * phi_max)/(omega_step * phi_step); - STATUS("%i templates to be calculated.\n", n_templates); - - priv->templates = malloc(n_templates * sizeof(struct template)); - - i = 0; - for ( omega = 0.0; omega < omega_max-omega_step; omega += omega_step ) { - for ( phi = 0.0; phi < phi_max-phi_step; phi += phi_step ) { - - RefList *reflections; - UnitCell *cell_rot; - - assert(i < n_templates); - - cell_rot = rotate_cell(cell, omega, phi, 0.0); - - reflections = find_intersections(&image, cell_rot, 0); - if ( reflections == NULL ) { - ERROR("Template calculation failed.\n"); - return NULL; - } - - priv->templates[i].omega = omega; - priv->templates[i].phi = phi; - priv->templates[i].spots = reflections; - i++; - - free(cell_rot); - - } - progress_bar(omega*1000.0, (omega_max-2.0*omega_step)*1000.0, - "Generating templates"); - } - - priv->n_templates = n_templates; - priv->cell = cell_new_from_cell(cell); - - free(image.data); - free(image.flags); - - return (struct _indexingprivate *)priv; -} - - -static int fast_integrate_peak(struct image *image, int xp, int yp) -{ - int x, y; - double total = 0; - int r = INTEGRATION_SQUARE_SIDE; - - for ( x=xp-r; x<=xp+r; x++ ) { - for ( y=yp-r; y<=yp+r; y++ ) { - - int val; - - if ( (x>=image->width) || (x<0) ) continue; - if ( (y>=image->height) || (y<0) ) continue; - - val = image->data[x+image->width*y]; - - if ( val < THRESHOLD ) continue; - total += val; - - } - } - - return total; -} - - -static double integrate_all_rot(struct image *image, RefList *reflections, - double rot) -{ - double itot = 0.0; - double cosr, sinr; - int num_int = 0; - Reflection *refl; - RefListIterator *iter; - - cosr = cos(rot); - sinr = sin(rot); - - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - float xp, yp; - double x, y; - - get_detector_pos(refl, &x, &y); - - xp = cosr*x + sinr*y; - yp = -sinr*x + cosr*y; - - itot += fast_integrate_peak(image, rint(xp), rint(yp)); - num_int++; - - } - - return itot / num_int; -} - - -/* Return the mean of the distances between peaks in the image and peaks from - * the given template. */ -static double mean_distance(struct image *image, RefList *reflections, - double rot) -{ - double dtot = 0.0; - double cosr, sinr; - int num_dist = 0; - Reflection *refl; - RefListIterator *iter; - - cosr = cos(rot); - sinr = sin(rot); - - /* For each template peak */ - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - float xp, yp; - int j; - double min_dsq; - double x, y; - - get_detector_pos(refl, &x, &y); - - xp = cosr*x + sinr*y; - yp = -sinr*x + cosr*y; - - /* Compare to every real peak */ - min_dsq = +INFINITY; - for ( j=0; j<image_feature_count(image->features); j++ ) { - - struct imagefeature *f; - double this_dsq; - - f = image_get_feature(image->features, j); - - this_dsq = pow(f->x - xp, 2.0) + pow(f->y - yp, 2.0); - if ( this_dsq < min_dsq ) min_dsq = this_dsq; - - } - - if ( min_dsq < pow(50.0, 2.0) ) { - dtot += sqrt(min_dsq); - num_dist++; - } - - } - - return dtot / num_dist; -} - - -void match_templates(struct image *image, IndexingPrivate *ipriv) -{ - struct _indexingprivate_template *priv - = (struct _indexingprivate_template *)ipriv; - int i, max_i; - double max, rot, rot_max, rot_step, rot_best; - const int peaks = 0; - - if ( peaks ) { - max = INFINITY; - } else { - max = 0.0; - } - max_i = 0; - rot_max = 2.0*M_PI; - rot_step = 2.0*M_PI / 360.0; /* 1 deg steps */ - rot_best = 0.0; - - for ( i=0; i<priv->n_templates; i++ ) { - for ( rot=0.0; rot<rot_max; rot+=rot_step ) { - - double val; - int best; - - if ( !peaks ) { - val = integrate_all_rot(image, priv->templates[i].spots, - rot); - best = val > max; - } else { - val = mean_distance(image, priv->templates[i].spots, - rot); - best = val < max; - } - - if ( best ) { - max = val; - max_i = i; - rot_best = rot; - } - - } - progress_bar(i, priv->n_templates-1, "Indexing"); - } - - STATUS("%i (%.2f, %.2f, %.2f): %.2f\n", max_i, - rad2deg(priv->templates[max_i].omega), - rad2deg(priv->templates[max_i].phi), - rad2deg(rot_best), max); - - image->ncells = 1; - image->candidate_cells[0] = rotate_cell(priv->cell, - priv->templates[max_i].omega, - priv->templates[max_i].phi, - rot_best); -} - - -void free_templates(IndexingPrivate *priv) -{ - int i; - struct _indexingprivate_template *tpriv - = (struct _indexingprivate_template *)priv; - - for ( i=0; i<tpriv->n_templates; i++ ) { - free(tpriv->templates[i].spots); - } - - free(tpriv->templates); - cell_free(tpriv->cell); - free(tpriv); -} diff --git a/src/templates.h b/src/templates.h deleted file mode 100644 index 33e6c75e..00000000 --- a/src/templates.h +++ /dev/null @@ -1,33 +0,0 @@ -/* - * templates.h - * - * Indexing by template matching - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifndef TEMPLATES_H -#define TEMPLATES_H - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - - -#include "index.h" -#include "image.h" -#include "cell.h" - -extern IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, - struct detector *det, - double nominal_photon_energy); - - -extern void match_templates(struct image *image, IndexingPrivate *ipriv); - -extern void free_templates(IndexingPrivate *priv); - -#endif /* TEMPLATES_H */ |