diff options
author | Thomas White <taw@physics.org> | 2011-06-07 17:48:47 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:27 +0100 |
commit | b7050c2de9a6e990fac39554d65018b641397dee (patch) | |
tree | d94970d78b065a5ccc24afcdfff971ec1f3b89cc | |
parent | 6c06f74f4ab6d8089f55355abb4b137267f8ac6f (diff) |
Get rid of integrate_pixels()
-rw-r--r-- | src/peaks.c | 165 |
1 files changed, 0 insertions, 165 deletions
diff --git a/src/peaks.c b/src/peaks.c index e4a7cfb1..c542e44c 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -572,168 +572,3 @@ void integrate_reflections(struct image *image, int polar, int use_closer, } } - - -RefList *integrate_pixels(struct image *image, int circular_domain, - double domain_r, int do_polar) -{ - int i; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int fs, ss; - double aslen, bslen, cslen; - double *intensities; - double *xmom; - double *ymom; - ReflItemList *obs; - RefList *reflections; - - obs = new_items(); - intensities = new_list_intensity(); - xmom = new_list_intensity(); - ymom = new_list_intensity(); - reflections = reflist_new(); - - /* "Borrow" direction values to get reciprocal lengths */ - cell_get_reciprocal(image->indexed_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(image->indexed_cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz); - /* For each pixel */ - fesetround(1); /* Round towards nearest */ - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; ss++ ) { - - 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; - struct panel *p; - double twotheta; - - p = find_panel(image->det, fs, ss); - if ( p == NULL ) continue; - if ( p->no_index ) continue; - - q = get_q(image, fs, ss, &twotheta, 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; - double pix_area, Lsq, proj_area, dsq, sa; - double phi, pa, pb, pol; - double xs, ys, rx, ry; - - /* Veto if we want to integrate a bad region */ - if ( image->flags != NULL ) { - int flags; - flags = image->flags[fs+image->width*ss]; - if ( !(flags & 0x01) ) continue; - } - - val = image->data[fs+image->width*ss]; - - /* 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 */ - proj_area = pix_area * cos(twotheta); - - /* Calculate distance from crystal to pixel */ - xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx; - ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy; - rx = (xs + p->cnx) / p->res; - ry = (ys + p->cny) / p->res; - dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); - - /* Projected area of pixel / distance squared */ - sa = 1.0e7 * proj_area / (dsq + Lsq); - - /* Solid angle correction is needed in this case */ - val /= sa; - - if ( do_polar ) { - - phi = atan2(ry, rx); - pa = pow(sin(phi)*sin(twotheta), 2.0); - pb = pow(cos(twotheta), 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*fs); - integrate_intensity(ymom, h, k, l, val*ss); - - 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; - Reflection *refl; - - 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; - - refl = add_refl(reflections, it->h, it->k, it->l); - set_int(refl, intensity); - set_detector_pos(refl, 0.0, xp, yp); - - } - - free(xmom); - free(ymom); - free(intensities); - delete_items(obs); - - return reflections; -} |