From 0dcc9e2e1ea4fcbd37db5ab7ac74146098c2f4d7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 17 Aug 2010 18:14:44 +0200 Subject: Generation of templates (needs debugging) --- src/Makefile.am | 5 ++- src/Makefile.in | 12 ++--- src/cell.c | 3 ++ src/facetron.c | 135 +------------------------------------------------------- src/symmetry.c | 3 +- src/templates.c | 109 +++++++++++++++++++++++++++++++++++++++++++-- 6 files changed, 122 insertions(+), 145 deletions(-) (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index b7108918..7fdf6dc9 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -22,7 +22,8 @@ process_hkl_LDADD = @LIBS@ indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \ peaks.c index.c filters.c diffraction.c detector.c \ - sfac.c dirax.c reflections.c templates.c symmetry.c + sfac.c dirax.c reflections.c templates.c symmetry.c \ + geometry.c indexamajig_LDADD = @LIBS@ if HAVE_OPENCL indexamajig_SOURCES += diffraction-gpu.c cl-utils.c @@ -56,7 +57,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 diffraction.c sfac.c + image.c diffraction.c sfac.c geometry.c facetron_LDADD = @LIBS@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c diff --git a/src/Makefile.in b/src/Makefile.in index 9355ef39..4b4da1f7 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -71,7 +71,7 @@ cubeit_DEPENDENCIES = am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \ hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \ peaks.$(OBJEXT) image.$(OBJEXT) diffraction.$(OBJEXT) \ - sfac.$(OBJEXT) + sfac.$(OBJEXT) geometry.$(OBJEXT) facetron_OBJECTS = $(am_facetron_OBJECTS) facetron_DEPENDENCIES = am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \ @@ -89,7 +89,7 @@ hdfsee_DEPENDENCIES = am__indexamajig_SOURCES_DIST = indexamajig.c hdf5-file.c utils.c \ cell.c image.c peaks.c index.c filters.c diffraction.c \ detector.c sfac.c dirax.c reflections.c templates.c symmetry.c \ - diffraction-gpu.c cl-utils.c + geometry.c diffraction-gpu.c cl-utils.c @HAVE_OPENCL_TRUE@am__objects_1 = diffraction-gpu.$(OBJEXT) \ @HAVE_OPENCL_TRUE@ cl-utils.$(OBJEXT) am_indexamajig_OBJECTS = indexamajig.$(OBJEXT) hdf5-file.$(OBJEXT) \ @@ -97,7 +97,7 @@ am_indexamajig_OBJECTS = indexamajig.$(OBJEXT) hdf5-file.$(OBJEXT) \ index.$(OBJEXT) filters.$(OBJEXT) diffraction.$(OBJEXT) \ detector.$(OBJEXT) sfac.$(OBJEXT) dirax.$(OBJEXT) \ reflections.$(OBJEXT) templates.$(OBJEXT) symmetry.$(OBJEXT) \ - $(am__objects_1) + geometry.$(OBJEXT) $(am__objects_1) indexamajig_OBJECTS = $(am_indexamajig_OBJECTS) indexamajig_DEPENDENCIES = am__pattern_sim_SOURCES_DIST = pattern_sim.c diffraction.c utils.c \ @@ -260,7 +260,8 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \ process_hkl_LDADD = @LIBS@ indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \ peaks.c index.c filters.c diffraction.c detector.c sfac.c \ - dirax.c reflections.c templates.c symmetry.c $(am__append_3) + dirax.c reflections.c templates.c symmetry.c geometry.c \ + $(am__append_3) indexamajig_LDADD = @LIBS@ @HAVE_GTK_TRUE@hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c image.c \ @HAVE_GTK_TRUE@ filters.c @@ -286,7 +287,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 diffraction.c sfac.c + image.c diffraction.c sfac.c geometry.c facetron_LDADD = @LIBS@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c @@ -415,6 +416,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/displaywindow.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/facetron.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/filters.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geometry.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/get_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hdf5-file.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hdfsee.Po@am__quote@ diff --git a/src/cell.c b/src/cell.c index d6464a71..d6d101cc 100644 --- a/src/cell.c +++ b/src/cell.c @@ -207,6 +207,9 @@ static int cell_crystallographic_to_cartesian(UnitCell *cell, { double tmp, V, cosalphastar, cstar; + /* Note: Please consider and possibly change the ranges for template + * matching (in templates.c) if the calculations below are altered. */ + /* Firstly: Get a in terms of x, y and z * +a (cryst) is defined to lie along +x (cart) */ *ax = cell->a; diff --git a/src/facetron.c b/src/facetron.c index 0fa6b643..0c2b4d56 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -26,10 +26,7 @@ #include "cell.h" #include "hdf5-file.h" #include "detector.h" -#include "peaks.h" - - -#define MAX_HITS (1024) +#include "geometry.h" static void show_help(const char *s) @@ -116,136 +113,6 @@ static int find_chunk(FILE *fh, UnitCell **cell, char **filename) } -static struct reflhit *find_intersections(struct image *image, UnitCell *cell, - double divergence, double bandwidth, - int *n, int output) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - struct reflhit *hits; - int np = 0; - int hmax, kmax, lmax; - double mres; - signed int h, k, l; - - hits = malloc(sizeof(struct reflhit)*MAX_HITS); - if ( hits == NULL ) { - *n = 0; - return NULL; - } - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - - mres = 1.0 / 8.0e-10; /* 8 Angstroms */ - hmax = mres / modulus(asx, asy, asz); - kmax = mres / modulus(bsx, bsy, bsz); - lmax = mres / modulus(csx, csy, csz); - - for ( h=-hmax; hlambda - image->lambda*bandwidth/2.0; - lhigh = image->lambda + image->lambda*bandwidth/2.0; - - /* Get the coordinates of the reciprocal lattice point */ - zl = h*asz + k*bsz + l*csz; - if ( zl < 0.0 ) continue; /* Do this check very early */ - xl = h*asx + k*bsx + l*csx; - yl = h*asy + k*bsy + l*csy; - - ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */ - delta = divergence/image->lambda; - dps_sq = ds_sq + pow(delta, 2.0); /* d'*^2 */ - - /* In range? */ - divfact = 2.0 * delta * sqrt(xl*xl + yl*yl); - if ( ds_sq - 2.0*zl/llow > 0.0 ) continue; - if ( ds_sq - 2.0*zl/lhigh < 0.0 ) continue; - - /* Work out which panel this peak would fall on */ - for ( p=0; pdet->n_panels; p++ ) { - - /* Camera length for this panel */ - cl = image->det->panels[p].clen; - - /* Coordinates of peak relative to central beam, in m */ - xd = cl*xl / (ds_sq/(2.0*zl) - zl); - yd = cl*yl / (ds_sq/(2.0*zl) - zl); - - /* Convert to pixels */ - xd *= image->det->panels[p].res; - yd *= image->det->panels[p].res; - - /* Add the coordinates of the central beam */ - xda = xd + image->det->panels[p].cx; - yda = yd + image->det->panels[p].cy; - - /* Now, is this on this panel? */ - if ( xda < image->det->panels[p].min_x ) continue; - if ( xda > image->det->panels[p].max_x ) continue; - if ( yda < image->det->panels[p].min_y ) continue; - if ( yda > image->det->panels[p].max_y ) continue; - - /* Woohoo! */ - found = 1; - break; - - } - - if ( !found ) continue; - - hits[np].h = h; - hits[np].k = k; - hits[np].l = l; - np++; - - if ( output ) { - printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda); - } - - } - } - } - - *n = np; - return hits; -} - - -static double integrate_all(struct image *image, struct reflhit *hits, int n) -{ - double itot = 0.0; - int i; - - for ( i=0; ibase.indm = INDEXING_TEMPLATE; /* We can only distinguish orientations within the holohedral cell */ holo = get_holohedral(cell_get_pointgroup(cell)); - STATUS("%s\n", holo); /* These define the orientation in space */ if ( is_polyhedral(holo) ) { @@ -54,8 +127,38 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename) /* One more axis would define the rotation in the plane of the image */ - STATUS("Orientation ranges: %5.0f -> %5.0f, %5.0f -> %5.0f deg.\n", - 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max)); + STATUS("Orientation ranges in %s: %.0f-%.0f, %.0f-%.0f deg.\n", + holo, 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); + + for ( omega = 0.0; omega < omega_max; omega += omega_step ) { + for ( phi = 0.0; phi < phi_max; phi += phi_step ) { + + int n; + struct reflhit *hits; + UnitCell *cell_rot; + + cell_rot = rotate_cell(cell, omega, phi); + + hits = find_intersections(&image, cell_rot, 5.0e-3, + 3.0/100.0, &n, 0); + if ( hits == NULL ) { + ERROR("Template calculation failed.\n"); + return NULL; + } + + STATUS("%.2f, %.2f : %i features\n", + rad2deg(omega), rad2deg(phi), n); + + free(cell_rot); + + } + STATUS("Finished omega=%.2f\n", rad2deg(omega)); + } + + priv->n_templates = n_templates; return (struct _indexingprivate *)priv; } -- cgit v1.2.3