aboutsummaryrefslogtreecommitdiff
path: root/src/peaks.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-10-11 18:10:22 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:02 +0100
commit7dcd21fc5ebf18d39f145c92f5c2038ee6c0d0a5 (patch)
tree9844d03675f5f78e7ef7b497f2f0637de2de4928 /src/peaks.c
parent29f669dcbd11af84edaff1c2bb417a20bed2b1ba (diff)
Add a new way of doing the peak integration
Diffstat (limited to 'src/peaks.c')
-rw-r--r--src/peaks.c223
1 files changed, 203 insertions, 20 deletions
diff --git a/src/peaks.c b/src/peaks.c
index 313af7a9..99598fa1 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -630,40 +630,23 @@ int peak_sanity_check(struct image *image, UnitCell *cell,
}
-void output_intensities(struct image *image, UnitCell *cell,
- pthread_mutex_t *mutex, int polar, int sa,
- int use_closer, FILE *ofh,
- int circular_domain, double domain_r)
+static void output_header(FILE *ofh, UnitCell *cell, struct image *image)
{
- int i;
- int n_found;
- int n_indclose = 0;
- int n_foundclose = 0;
- int n_veto = 0;
- int n_veto_second = 0;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
double a, b, c, al, be, ga;
- if ( image->n_hits == 0 ) {
- find_projected_peaks(image, cell, circular_domain, domain_r);
- }
- if ( image->n_hits == 0 ) return;
-
- /* Get exclusive access to the output stream if necessary */
- if ( mutex != NULL ) pthread_mutex_lock(mutex);
-
- /* Explicit printf() used here (not normally allowed) because
- * we really want to output to stdout */
fprintf(ofh, "Reflections from indexing in %s\n", image->filename);
fprintf(ofh, "Orientation (wxyz): %7.5f %7.5f %7.5f %7.5f\n",
image->orientation.w, image->orientation.x,
image->orientation.y, image->orientation.z);
+
cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm, %7.5f %7.5f %7.5f deg\n",
a*1.0e9, b*1.0e9, c*1.0e9,
rad2deg(al), rad2deg(be), rad2deg(ga));
+
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
@@ -681,6 +664,38 @@ void output_intensities(struct image *image, UnitCell *cell,
fprintf(ofh, "f0 = invalid\n");
}
+}
+
+
+void output_intensities(struct image *image, UnitCell *cell,
+ pthread_mutex_t *mutex, int polar, int sa,
+ int use_closer, FILE *ofh,
+ int circular_domain, double domain_r)
+{
+ int i;
+ int n_found;
+ int n_indclose = 0;
+ int n_foundclose = 0;
+ int n_veto = 0;
+ int n_veto_second = 0;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ if ( image->n_hits == 0 ) {
+ find_projected_peaks(image, cell, circular_domain, domain_r);
+ }
+ if ( image->n_hits == 0 ) return;
+
+ /* Get exclusive access to the output stream if necessary */
+ if ( mutex != NULL ) pthread_mutex_lock(mutex);
+
+ output_header(ofh, cell, image);
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
for ( i=0; i<image->n_hits; i++ ) {
float x, y, intensity;
@@ -807,3 +822,171 @@ void output_intensities(struct image *image, UnitCell *cell,
if ( mutex != NULL ) pthread_mutex_unlock(mutex);
}
+
+
+void output_pixels(struct image *image, UnitCell *cell,
+ pthread_mutex_t *mutex, int do_polar, int do_sa,
+ FILE *ofh, int circular_domain, double domain_r)
+{
+ int i;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ int x, y;
+ double aslen, bslen, cslen;
+ double *intensities;
+ double *xmom;
+ double *ymom;
+ ReflItemList *obs;
+
+ /* Get exclusive access to the output stream if necessary */
+ if ( mutex != NULL ) pthread_mutex_lock(mutex);
+
+ output_header(ofh, cell, image);
+
+ obs = new_items();
+ intensities = new_list_intensity();
+ xmom = new_list_intensity();
+ ymom = new_list_intensity();
+
+ /* "Borrow" direction values to get reciprocal lengths */
+ cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ aslen = modulus(ax, ay, az);
+ bslen = modulus(bx, by, bz);
+ cslen = modulus(cx, cy, cz);
+
+ cell_get_cartesian(cell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
+ /* For each pixel */
+ fesetround(1); /* Round towards nearest */
+ for ( x=0; x<image->width; x++ ) {
+ for ( y=0; y<image->height; y++ ) {
+
+ double hd, kd, ld; /* Indices with decimal places */
+ double dh, dk, dl; /* Distances in h,k,l directions */
+ signed int h, k, l;
+ struct rvec q;
+ double dist;
+
+ q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
+
+ 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);
+
+ dh = hd - h; dk = kd - k; dl = ld - l;
+
+ if ( circular_domain ) {
+
+ /* Circular integration domain */
+ dist = sqrt(pow(dh*aslen, 2.0) + pow(dk*bslen, 2.0)
+ + pow(dl*cslen, 2.0));
+
+ } else {
+
+ /* "Crystallographic" integration domain */
+ dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
+ }
+
+ if ( dist < domain_r ) {
+
+ double val;
+ struct panel *p;
+ double pix_area, Lsq, proj_area, dsq, sa;
+ double phi, pa, pb, pol;
+ float tt = 0.0;
+
+ /* Veto if we want to integrate a bad region */
+ if ( image->flags != NULL ) {
+ int flags;
+ flags = image->flags[x+image->width*y];
+ if ( !(flags & 0x01) ) continue;
+ }
+
+ val = image->data[x+image->width*y];
+
+ p = find_panel(image->det, x, y);
+ if ( p == NULL ) continue;
+
+ if ( p->no_index ) continue;
+
+ if ( do_sa ) {
+
+ /* Area of one pixel */
+ pix_area = pow(1.0/p->res, 2.0);
+ Lsq = pow(p->clen, 2.0);
+
+ /* Area of pixel as seen from crystal */
+ tt = get_tt(image, x, y);
+ proj_area = pix_area * cos(tt);
+
+ /* Calculate distance from crystal to pixel */
+ dsq = pow(((double)x - p->cx) / p->res, 2.0);
+ dsq += pow(((double)y - p->cy) / p->res, 2.0);
+
+ /* Projected area of pixel / distance squared */
+ sa = 1.0e7 * proj_area / (dsq + Lsq);
+
+ val /= sa;
+
+ }
+
+ if ( do_polar ) {
+
+ if ( !do_sa ) tt = get_tt(image, x, y);
+
+ phi = atan2(y, x);
+ pa = pow(sin(phi)*sin(tt), 2.0);
+ pb = pow(cos(tt), 2.0);
+ pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
+
+ val /= pol;
+
+ }
+
+ /* Add value to sum */
+ integrate_intensity(intensities, h, k, l, val);
+
+ integrate_intensity(xmom, h, k, l, val*x);
+ integrate_intensity(ymom, h, k, l, val*y);
+
+ if ( !find_item(obs, h, k, l) ) {
+ add_item(obs, h, k, l);
+ }
+
+ }
+
+ }
+ }
+
+ for ( i=0; i<num_items(obs); i++ ) {
+
+ struct refl_item *it;
+ double intensity, xmomv, ymomv;
+ double xp, yp;
+
+ it = get_item(obs, i);
+ intensity = lookup_intensity(intensities, it->h, it->k, it->l);
+ xmomv = lookup_intensity(xmom, it->h, it->k, it->l);
+ ymomv = lookup_intensity(ymom, it->h, it->k, it->l);
+
+ xp = xmomv / (double)intensity;
+ yp = ymomv / (double)intensity;
+
+ fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f)\n",
+ image->hits[i].h, image->hits[i].k, image->hits[i].l,
+ intensity, xp, yp);
+
+ }
+
+ fprintf(ofh, "No peak statistics, because output_pixels() was used.");
+ /* Blank line at end */
+ fprintf(ofh, "\n");
+
+ if ( mutex != NULL ) pthread_mutex_unlock(mutex);
+}