Commit 329134c4 authored by Thomas White's avatar Thomas White Committed by Thomas White
Browse files

Remove solid angle correction

It's never correct when using "bucket" integration, and always correct
when using "pixel" integration, so don't give the option.
parent 5bc60000
......@@ -283,7 +283,7 @@ double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
float x, y, intensity;
if ( integrate_peak(image, cpeaks[i].x, cpeaks[i].y, &x, &y,
&intensity, NULL, NULL, 0, 0, 0) ) continue;
&intensity, NULL, NULL, 0, 0) ) continue;
itot += intensity;
}
......
......@@ -62,7 +62,6 @@ struct static_index_args
int config_polar;
int config_sanity;
int config_satcorr;
int config_sa;
int config_closer;
float threshold;
float min_gradient;
......@@ -182,8 +181,6 @@ static void show_help(const char *s)
" --unpolarized Don't correct for the polarisation of the X-rays.\n"
" --no-sat-corr Don't correct values of saturated peaks using a\n"
" table included in the HDF5 file.\n"
" --no-sa Don't correct for the differing solid angles of\n"
" the pixels.\n"
" --threshold=<n> Only accept peaks above <n> ADU. Default: 800.\n"
" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
" Default: 100,000.\n"
......@@ -415,7 +412,7 @@ static void process_image(void *pp, int cookie)
if ( config_nearbragg ) {
output_intensities(&image, image.indexed_cell,
pargs->static_args.output_mutex,
config_polar, pargs->static_args.config_sa,
config_polar,
pargs->static_args.config_closer,
pargs->static_args.ofh, 0, 0.1);
}
......@@ -525,7 +522,6 @@ int main(int argc, char *argv[])
int config_polar = 1;
int config_sanity = 0;
int config_satcorr = 1;
int config_sa = 1;
int config_checkprefix = 1;
int config_closer = 1;
float threshold = 800.0;
......@@ -586,7 +582,6 @@ int main(int argc, char *argv[])
{"check-sanity", 0, &config_sanity, 1},
{"no-sat-corr", 0, &config_satcorr, 0},
{"sat-corr", 0, &config_satcorr, 1}, /* Compat */
{"no-sa", 0, &config_sa, 0},
{"threshold", 1, NULL, 't'},
{"min-gradient", 1, NULL, 4},
{"no-check-prefix", 0, &config_checkprefix, 0},
......@@ -882,7 +877,6 @@ int main(int argc, char *argv[])
qargs.static_args.config_polar = config_polar;
qargs.static_args.config_sanity = config_sanity;
qargs.static_args.config_satcorr = config_satcorr;
qargs.static_args.config_sa = config_sa;
qargs.static_args.config_closer = config_closer;
qargs.static_args.cellr = cellr;
qargs.static_args.threshold = threshold;
......
......@@ -208,8 +208,7 @@ static void integrate_image(struct image *image, ReflItemList *obs)
* pattern? */
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[j].x, spots[j].y,
&xc, &yc, &i_partial, NULL, NULL,
1, 1, 0) ) {
&xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
spots[j].valid = 0;
continue;
}
......
......@@ -565,7 +565,7 @@ int main(int argc, char *argv[])
if ( config_nearbragg ) {
find_projected_peaks(&image, cell, 0, 0.1);
output_intensities(&image, cell, NULL, 0, 1, 0, stdout,
output_intensities(&image, cell, NULL, 0, 0, stdout,
0, 0.1);
free(image.cpeaks);
}
......
......@@ -182,7 +182,7 @@ static int cull_peaks(struct image *image)
int integrate_peak(struct image *image, int xp, int yp,
float *xc, float *yc, float *intensity,
double *pbg, double *pmax,
int do_polar, int do_sa, int centroid)
int do_polar, int centroid)
{
signed int x, y;
int lim, out_lim;
......@@ -204,7 +204,7 @@ int integrate_peak(struct image *image, int xp, int yp,
for ( x=-out_lim; x<+out_lim; x++ ) {
for ( y=-out_lim; y<+out_lim; y++ ) {
double val, sa, pix_area, Lsq, dsq, proj_area;
double val;
float tt = 0.0;
double phi, pa, pb, pol;
uint16_t flags;
......@@ -238,30 +238,9 @@ int integrate_peak(struct image *image, int xp, int yp,
if ( val > max ) max = val;
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 (approximate) */
tt = get_tt(image, x+xp, y+yp);
proj_area = pix_area * cos(tt);
/* Calculate distance from crystal to pixel */
dsq = pow(((double)(x+xp) - p->cx) / p->res, 2.0);
dsq += pow(((double)(y+yp) - p->cy) / p->res, 2.0);
/* Projected area of pixel divided by distance squared */
sa = 1.0e7 * proj_area / (dsq + Lsq);
val /= sa;
}
if ( do_polar ) {
if ( !do_sa ) tt = get_tt(image, x+xp, y+yp);
tt = get_tt(image, x+xp, y+yp);
phi = atan2(y+yp, x+xp);
pa = pow(sin(phi)*sin(tt), 2.0);
......@@ -420,7 +399,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
* Don't bother doing polarisation/SA correction, because the
* intensity of this peak is only an estimate at this stage. */
r = integrate_peak(image, mask_x, mask_y,
&fx, &fy, &intensity, NULL, NULL, 0, 0, 1);
&fx, &fy, &intensity, NULL, NULL, 0, 1);
if ( r ) {
/* Bad region - don't detect peak */
nrej_bad++;
......@@ -702,7 +681,7 @@ static void output_header(FILE *ofh, UnitCell *cell, struct image *image)
void output_intensities(struct image *image, UnitCell *cell,
pthread_mutex_t *mutex, int polar, int sa,
pthread_mutex_t *mutex, int polar,
int use_closer, FILE *ofh,
int circular_domain, double domain_r)
{
......@@ -766,7 +745,7 @@ void output_intensities(struct image *image, UnitCell *cell,
* revised coordinates. */
r = integrate_peak(image, f->x, f->y, &x, &y,
&intensity, &bg, &max,
polar, sa, 1);
polar, 1);
if ( r ) {
/* The original peak (which also went
* through integrate_peak(), but with
......@@ -788,7 +767,7 @@ void output_intensities(struct image *image, UnitCell *cell,
image->cpeaks[i].x,
image->cpeaks[i].y,
&x, &y, &intensity, &bg, &max,
polar, sa, 1);
polar, 1);
if ( r ) {
/* Plain old ordinary peak veto */
n_veto++;
......@@ -809,7 +788,7 @@ void output_intensities(struct image *image, UnitCell *cell,
image->cpeaks[i].x,
image->cpeaks[i].y,
&x, &y, &intensity, &bg, &max, polar,
sa, 0);
0);
if ( r ) {
/* Plain old ordinary peak veto */
n_veto++;
......@@ -867,7 +846,7 @@ void output_intensities(struct image *image, UnitCell *cell,
void output_pixels(struct image *image, UnitCell *cell,
pthread_mutex_t *mutex, int do_polar, int do_sa,
pthread_mutex_t *mutex, int do_polar,
FILE *ofh, int circular_domain, double domain_r)
{
int i;
......@@ -962,30 +941,27 @@ void output_pixels(struct image *image, UnitCell *cell,
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);
/* Area of one pixel */
pix_area = pow(1.0/p->res, 2.0);
Lsq = pow(p->clen, 2.0);
/* 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);
/* Area of pixel as seen from crystal */
tt = get_tt(image, x, y);
proj_area = pix_area * cos(tt);
/* Projected area of pixel / distance squared */
sa = 1.0e7 * proj_area / (dsq + Lsq);
/* 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);
val /= sa;
/* 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 ) {
if ( !do_sa ) tt = get_tt(image, x, y);
tt = get_tt(image, x, y);
phi = atan2(y, x);
pa = pow(sin(phi)*sin(tt), 2.0);
......
......@@ -24,12 +24,12 @@ extern void search_peaks(struct image *image, float threshold,
extern void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex);
extern void output_intensities(struct image *image, UnitCell *cell,
pthread_mutex_t *mutex, int polar, int sa,
pthread_mutex_t *mutex, int polar,
int use_closer, FILE *ofh, int circular_domain,
double domain_r);
extern void output_pixels(struct image *image, UnitCell *cell,
pthread_mutex_t *mutex, int do_polar, int do_sa,
pthread_mutex_t *mutex, int do_polar,
FILE *ofh, int circular_domain, double domain_r);
extern int peak_sanity_check(struct image *image, UnitCell *cell,
......@@ -39,6 +39,6 @@ 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,
double *pbg, double *pmax,
int do_polar, int do_sa, int centroid);
int do_polar, int centroid);
#endif /* PEAKS_H */
......@@ -246,7 +246,7 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL,
1, 1, 0) ) {
1, 0) ) {
continue;
}
......@@ -318,7 +318,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
/* FIXME: Coordinates aren't whole numbers */
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL,
1, 1, 0) ) {
1, 0) ) {
continue;
}
......
......@@ -140,8 +140,7 @@ static void process_image(void *pg, int cookie)
output_intensities(&image, apargs->cell,
pargs->output_mutex, pargs->config_polar,
pargs->config_sa, pargs->config_closer,
pargs->ofh, 0, 0.1);
pargs->config_closer, pargs->ofh, 0, 0.1);
}
free(image.data);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment