aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Makefile.am3
-rw-r--r--src/Makefile.in8
-rw-r--r--src/cell.c15
-rw-r--r--src/cell.h4
-rw-r--r--src/facetron.c260
-rw-r--r--src/peaks.c6
-rw-r--r--src/peaks.h4
7 files changed, 239 insertions, 61 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index be423094..34f8a9a0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -55,7 +55,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
sfac.c
calibrate_detector_LDADD = @LIBS@
-facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c
+facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
+ image.c diffraction.c sfac.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 4fae1697..4716d60e 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -69,7 +69,9 @@ am_cubeit_OBJECTS = cubeit.$(OBJEXT) cell.$(OBJEXT) \
cubeit_OBJECTS = $(am_cubeit_OBJECTS)
cubeit_DEPENDENCIES =
am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \
- hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT)
+ hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \
+ peaks.$(OBJEXT) image.$(OBJEXT) diffraction.$(OBJEXT) \
+ sfac.$(OBJEXT)
facetron_OBJECTS = $(am_facetron_OBJECTS)
facetron_DEPENDENCIES =
am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
@@ -282,7 +284,9 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
sfac.c
calibrate_detector_LDADD = @LIBS@
-facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c
+facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
+ image.c diffraction.c sfac.c
+
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c
cubeit_LDADD = @LIBS@
diff --git a/src/cell.c b/src/cell.c
index 79261249..cc252b70 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -179,6 +179,21 @@ UnitCell *cell_new_from_cell(UnitCell *orig)
}
+void cell_set_reciprocal(UnitCell *cell,
+ double asx, double asy, double asz,
+ double bsx, double bsy, double bsz,
+ double csx, double csy, double csz)
+{
+ if ( cell == NULL ) return;
+
+ cell->axs = asx; cell->ays = asy; cell->azs = asz;
+ cell->bxs = bsx; cell->bys = bsy; cell->bzs = bsz;
+ cell->cxs = csx; cell->cys = csy; cell->czs = csz;
+
+ cell->rep = CELL_REP_RECIP;
+}
+
+
/************************* Getter helper functions ****************************/
static int cell_crystallographic_to_cartesian(UnitCell *cell,
diff --git a/src/cell.h b/src/cell.h
index 113e6512..42ef2eb6 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -62,6 +62,10 @@ extern int cell_get_reciprocal(UnitCell *cell,
double *bsx, double *bsy, double *bsz,
double *csx, double *csy, double *csz);
+extern void cell_set_reciprocal(UnitCell *cell,
+ double asx, double asy, double asz,
+ double bsx, double bsy, double bsz,
+ double csx, double csy, double csz);
extern double resolution(UnitCell *cell,
signed int h, signed int k, signed int l);
diff --git a/src/facetron.c b/src/facetron.c
index 7a1eec53..0fa6b643 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -26,6 +26,7 @@
#include "cell.h"
#include "hdf5-file.h"
#include "detector.h"
+#include "peaks.h"
#define MAX_HITS (1024)
@@ -46,9 +47,78 @@ static void show_help(const char *s)
}
-/* Predict peaks */
-static int find_intersections(struct image *image, UnitCell *cell,
- double divergence, double bandwidth)
+static UnitCell *read_orientation_matrix(FILE *fh)
+{
+ float u, v, w;
+ struct rvec as, bs, cs;
+ UnitCell *cell;
+ char line[1024];
+
+ if ( fgets(line, 1023, fh) == NULL ) return NULL;
+ if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
+ ERROR("Couldn't read a-star\n");
+ return NULL;
+ }
+ as.u = u*1e9; as.v = v*1e9; as.w = w*1e9;
+ if ( fgets(line, 1023, fh) == NULL ) return NULL;
+ if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
+ ERROR("Couldn't read b-star\n");
+ return NULL;
+ }
+ bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9;
+ if ( fgets(line, 1023, fh) == NULL ) return NULL;
+ if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
+ ERROR("Couldn't read c-star\n");
+ return NULL;
+ }
+ cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9;
+ cell = cell_new_from_axes(as, bs, cs);
+
+ return cell;
+}
+
+
+static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
+{
+ char line[1024];
+ char *rval = NULL;
+
+ do {
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+
+ chomp(line);
+
+ if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
+ continue;
+ }
+
+ *filename = strdup(line+29);
+
+ /* Skip two lines (while checking for errors) */
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+
+ *cell = read_orientation_matrix(fh);
+ if ( *cell == NULL ) {
+ STATUS("Got filename but no cell for %s\n", *filename);
+ continue;
+ }
+
+ return 0;
+
+ } while ( rval != NULL );
+
+ return 1;
+}
+
+
+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;
@@ -60,7 +130,10 @@ static int find_intersections(struct image *image, UnitCell *cell,
signed int h, k, l;
hits = malloc(sizeof(struct reflhit)*MAX_HITS);
- if ( hits == NULL ) return 0;
+ if ( hits == NULL ) {
+ *n = 0;
+ return NULL;
+ }
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
@@ -140,84 +213,155 @@ static int find_intersections(struct image *image, UnitCell *cell,
hits[np].k = k;
hits[np].l = l;
np++;
- printf("%i %i %i 0.0 (at %f,%f) %i\n", h, k, l, xda, yda, p);
+
+ if ( output ) {
+ printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda);
+ }
}
}
}
- free(hits); /* FIXME! */
-
- return np;
+ *n = np;
+ return hits;
}
-static UnitCell *read_orientation_matrix(FILE *fh)
+static double integrate_all(struct image *image, struct reflhit *hits, int n)
{
- float u, v, w;
- struct rvec as, bs, cs;
- UnitCell *cell;
- char line[1024];
+ double itot = 0.0;
+ int i;
- if ( fgets(line, 1023, fh) == NULL ) return NULL;
- if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
- ERROR("Couldn't read a-star\n");
- return NULL;
- }
- as.u = u*1e9; as.v = v*1e9; as.w = w*1e9;
- if ( fgets(line, 1023, fh) == NULL ) return NULL;
- if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
- ERROR("Couldn't read b-star\n");
- return NULL;
+ for ( i=0; i<n; i++ ) {
+
+ float x, y, intensity;
+
+ if ( integrate_peak(image, hits[i].x, hits[i].y, &x, &y,
+ &intensity, 0, 0) ) continue;
+
+ itot += intensity;
}
- bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9;
- if ( fgets(line, 1023, fh) == NULL ) return NULL;
- if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
- ERROR("Couldn't read c-star\n");
- return NULL;
+
+ return itot;
+}
+
+
+static void get_trial_cell(UnitCell *out, UnitCell *in, int i, double step)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ switch ( i ) {
+ case 0 : asx += step; break;
+ case 1 : asx -= step; break;
+ case 2 : asy += step; break;
+ case 3 : asy -= step; break;
+ case 4 : asz += step; break;
+ case 5 : asz -= step; break;
+ case 6 : bsx += step; break;
+ case 7 : bsx -= step; break;
+ case 8 : bsy += step; break;
+ case 9 : bsy -= step; break;
+ case 10 : bsz += step; break;
+ case 11 : bsz -= step; break;
+ case 12 : csx += step; break;
+ case 13 : csx -= step; break;
+ case 14 : csy += step; break;
+ case 15 : csy -= step; break;
+ case 16 : csz += step; break;
+ case 17 : csz -= step; break;
+ default : break;
}
- cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9;
- cell = cell_new_from_axes(as, bs, cs);
- return cell;
+ cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
}
-static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
+static void try_refine(struct image *image, UnitCell *cell,
+ double da, double dw, double step)
{
- char line[1024];
- char *rval = NULL;
+ struct reflhit *hits;
+ int np;
+ double itot;
+ UnitCell *trial_cell;
+ int did_something;
+
+ trial_cell = cell_new();
+
+ hits = find_intersections(image, cell, da, dw, &np, 0);
+ itot = integrate_all(image, hits, np);
+ STATUS("%f\n", itot);
do {
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
+ int i;
- chomp(line);
+ did_something = 0;
- if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
- continue;
- }
+ for ( i=0; i<18; i++ ) {
- *filename = strdup(line+29);
+ struct reflhit *trial_hits;
+ double trial_itot;
- /* Skip two lines (while checking for errors) */
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
+ get_trial_cell(trial_cell, cell, i, step);
- *cell = read_orientation_matrix(fh);
- if ( *cell == NULL ) {
- STATUS("Got filename but no cell for %s\n", *filename);
- continue;
+ trial_hits = find_intersections(image, trial_cell,
+ da, dw, &np, 0);
+ trial_itot = integrate_all(image, hits, np);
+
+ if ( trial_itot > itot ) {
+
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ cell_get_reciprocal(trial_cell, &asx, &asy,
+ &asz, &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ cell_set_reciprocal(cell, asx, asy, asz, bsx,
+ bsy, bsz, csx, csy, csz);
+
+ itot = trial_itot;
+ free(hits);
+ hits = trial_hits;
+
+ did_something = 1;
+
+ } else {
+
+ free(trial_hits);
+
+ }
}
- return 0;
+ } while ( did_something );
- } while ( rval != NULL );
+ free(trial_cell);
+}
- return 1;
+
+/* Predict peaks */
+static void pre_refine(struct image *image, UnitCell *cell,
+ double *da, double *dw)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double istep, step;
+
+ /* Start by changing parameters by 1% */
+ cell_get_reciprocal(cell, &asx, &asy,
+ &asz, &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ istep = (asx+asy+asz+bsx+bsy+bsz+csx+csy+csz)/900.0;
+
+ for ( step=istep; step>istep/100.0; step-=istep/100.0 ) {
+ try_refine(image, cell, *da, *dw, step);
+ }
}
@@ -288,6 +432,8 @@ int main(int argc, char *argv[])
struct image image;
struct hdfile *hdfile;
+ double da, dw;
+ int np;
STATUS("Integrating intensities from '%s'\n", filename);
@@ -302,7 +448,13 @@ int main(int argc, char *argv[])
}
hdf5_read(hdfile, &image, 0);
- find_intersections(&image, cell, 5.0e-3, 3.0/100.0);
+
+ da = 5.0e-3; /* Initial divergence */
+ dw = 3.0/100.0; /* Initial bandwidth */
+
+ pre_refine(&image, cell, &da, &dw);
+
+ find_intersections(&image, cell, da, dw, &np, 1);
hdfile_close(hdfile);
free(cell);
diff --git a/src/peaks.c b/src/peaks.c
index 5f0a8a27..383f4663 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -169,9 +169,9 @@ static int cull_peaks(struct image *image)
/* Returns non-zero if peak has been vetoed.
* i.e. don't use result if return value is not zero. */
-static int integrate_peak(struct image *image, int xp, int yp,
- float *xc, float *yc, float *intensity,
- int do_polar, int do_sa)
+int integrate_peak(struct image *image, int xp, int yp,
+ float *xc, float *yc, float *intensity,
+ int do_polar, int do_sa)
{
signed int x, y;
const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
diff --git a/src/peaks.h b/src/peaks.h
index 2a24c5fa..f262a54e 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -25,5 +25,7 @@ extern void output_intensities(struct image *image, UnitCell *cell,
pthread_mutex_t *mutex, int polar, int sa);
extern int peak_sanity_check(struct image *image, UnitCell *cell);
extern int find_projected_peaks(struct image *image, UnitCell *cell);
-
+extern int integrate_peak(struct image *image, int xp, int yp,
+ float *xc, float *yc, float *intensity,
+ int do_polar, int do_sa);
#endif /* PEAKS_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-25 11:32:01 +0000