Commit 81d0cc92 authored by Thomas White's avatar Thomas White
Browse files

Merge branch 'master' into tom/speed

Conflicts:
	libcrystfel/src/peaks.c
	src/indexamajig.c
parents a4ffb65f d7cacb78
* Thomas White <taw@physics.org>
Lead author, architecture, "vision" and so on.
* Kenneth Beyerlein <kenneth.beyerlein@desy.de>
Peak integration improvements
* Richard Kirian <rkirian@asu.edu>
MOSFLM auto-indexing, bug fixing.
......
......@@ -7,6 +7,7 @@ Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
Authors:
Thomas White <taw@physics.org>
Richard Kirian <rkirian@asu.edu>
Kenneth Beyerlein <kenneth.beyerlein@desy.de>
Andrew Aquila <andrew.aquila@cfel.de>
Andrew Martin <andrew.martin@desy.de>
Lorenzo Galli <lorenzo.galli@desy.de>
......
......@@ -33,7 +33,7 @@ const sampler_t sampler_c = CLK_NORMALIZED_COORDS_TRUE
float4 get_q(float fs, float ss, float res, float clen, float k,
float *ttp, float corner_x, float corner_y,
float corner_x, float corner_y,
float fsx, float fsy, float ssx, float ssy)
{
float rx, ry, r;
......@@ -51,7 +51,6 @@ float4 get_q(float fs, float ss, float res, float clen, float k,
r = sqrt(pow(rx, 2.0f) + pow(ry, 2.0f));
tt = atan2(r, clen);
*ttp = tt;
az = atan2(ry, rx);
......@@ -211,7 +210,7 @@ float molecule_factor(global float *intensities, global float *flags,
}
kernel void diffraction(global float *diff, global float *ttp, float k,
kernel void diffraction(global float *diff, float k,
int w, float corner_x, float corner_y,
float res, float clen, float16 cell,
global float *intensities,
......@@ -219,7 +218,8 @@ kernel void diffraction(global float *diff, global float *ttp, float k,
read_only image2d_t func_b,
read_only image2d_t func_c,
global float *flags,
float fsx, float fsy, float ssx, float ssy)
float fsx, float fsy, float ssx, float ssy,
float xo, float yo)
{
float tt;
float fs, ss;
......@@ -230,11 +230,11 @@ kernel void diffraction(global float *diff, global float *ttp, float k,
int idx;
/* Calculate fractional coordinates in fs/ss */
fs = convert_float(get_global_id(0));
ss = convert_float(get_global_id(1));
fs = convert_float(get_global_id(0)) + xo;
ss = convert_float(get_global_id(1)) + yo;
/* Get the scattering vector */
q = get_q(fs, ss, res, clen, k, &tt,
q = get_q(fs, ss, res, clen, k,
corner_x, corner_y, fsx, fsy, ssx, ssy);
/* Calculate the diffraction */
......@@ -245,5 +245,4 @@ kernel void diffraction(global float *diff, global float *ttp, float k,
/* Write the value to memory */
idx = convert_int_rtz(fs) + w*convert_int_rtz(ss);
diff[idx] = I_molecule * I_lattice;
ttp[idx] = tt;
}
......@@ -47,7 +47,7 @@ of pixel values in the HDF5 file, defined in terms only of the "fast scan" and
+x completes the right-handed coordinate system.
.PP
Naively speaking, this means that CrystFEL at the images from the "into the
Naively speaking, this means that CrystFEL looks at the images from the "into the
beam" perspective, but please avoid thinking of things in this way. It's much
better to consider the precise way in which the coordinates are mapped.
......@@ -116,6 +116,8 @@ panel0/ss = -x
.br
; the HDF5 file, is now given a position in the lab coordinate system.
.br
; Units are pixels.
.br
; Note that "first point in the panel" is a conceptual simplification. We refer
.br
; to that corner, and to the very corner of the pixel - NOT, for example, to the
......
......@@ -608,8 +608,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key,
panel->res = atof(val);
} else if ( strcmp(key, "max_adu") == 0 ) {
panel->max_adu = atof(val);
} else if ( strcmp(key, "peak_sep") == 0 ) {
panel->peak_sep = atof(val);
} else if ( strcmp(key, "badrow_direction") == 0 ) {
panel->badrow = val[0]; /* First character only */
if ( (panel->badrow != 'x') && (panel->badrow != 'y')
......@@ -687,8 +685,6 @@ static void parse_toplevel(struct detector *det, const char *key,
det->mask_good = v;
}
} else if ( strcmp(key, "peak_sep") == 0 ) {
det->defaults.peak_sep = atof(val);
} else if ( strcmp(key, "coffset") == 0 ) {
det->defaults.coffset = atof(val);
} else if ( parse_field_for_panel(&det->defaults, key, val, det) ) {
......@@ -738,7 +734,6 @@ struct detector *get_detector_geometry(const char *filename)
det->defaults.res = -1.0;
det->defaults.badrow = '-';
det->defaults.no_index = 0;
det->defaults.peak_sep = 50.0;
det->defaults.fsx = 1.0;
det->defaults.fsy = 0.0;
det->defaults.ssx = 0.0;
......@@ -887,7 +882,6 @@ struct detector *get_detector_geometry(const char *filename)
}
/* It's OK if the badrow direction is '0' */
/* It's not a problem if "no_index" is still zero */
/* The default peak_sep is OK (maybe) */
/* The default transformation matrix is at least valid */
if ( det->panels[i].max_fs > max_fs ) {
......@@ -1250,7 +1244,6 @@ int write_detector_geometry(const char *filename, struct detector *det)
fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss);
fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow);
fprintf(fh, "%s/res = %g\n", p->name, p->res);
fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep);
fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from);
fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy);
fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy);
......
......@@ -62,7 +62,6 @@ struct panel
double res; /* Resolution in pixels per metre */
char badrow; /* 'x' or 'y' */
int no_index; /* Don't index peaks in this panel if non-zero */
double peak_sep; /* Characteristic peak separation */
char *rigid_group; /* Rigid group, or -1 for none */
double adu_per_eV; /* Number of ADU per eV */
double max_adu; /* Treat pixel as unreliable if higher than this */
......
......@@ -9,6 +9,7 @@
*
* Authors:
* 2010-2012 Thomas White <taw@physics.org>
* 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
* 2011 Andrew Martin <andrew.martin@desy.de>
* 2011 Richard Kirian
*
......@@ -98,7 +99,6 @@ static int cull_peaks_in_panel(struct image *image, struct panel *p)
if ( ncol <= 3 ) continue;
/* Yes? Delete them all... */
nelim = 0;
for ( j=0; j<n; j++ ) {
struct imagefeature *g;
g = image_get_feature(image->features, j);
......@@ -146,13 +146,76 @@ static int cull_peaks(struct image *image)
}
/* cfs, css relative to panel origin */
static int *make_BgMask(struct image *image, struct panel *p,
double ir_out, int cfs, int css, double ir_in)
{
Reflection *refl;
RefListIterator *iter;
int *mask;
int w, h;
w = p->max_fs - p->min_fs + 1;
h = p->max_ss - p->min_ss + 1;
mask = calloc(w*h, sizeof(int));
if ( mask == NULL ) return NULL;
/* Loop over all reflections */
for ( refl = first_refl(image->reflections, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
struct panel *p2;
double pk2_fs, pk2_ss;
signed int dfs, dss;
double pk2_cfs, pk2_css;
get_detector_pos(refl, &pk2_fs, &pk2_ss);
/* Determine if reflection is in the same panel */
p2 = find_panel(image->det, pk2_fs, pk2_ss);
if ( p2 != p ) continue;
pk2_cfs = pk2_fs - p->min_fs;
pk2_css = pk2_ss - p->min_ss;
for ( dfs=-ir_in; dfs<=ir_in; dfs++ ) {
for ( dss=-ir_in; dss<=ir_in; dss++ ) {
signed int fs, ss;
/* In peak region for this peak? */
if ( dfs*dfs + dss*dss > ir_in*ir_in ) continue;
fs = pk2_cfs + dfs;
ss = pk2_css + dss;
/* On panel? */
if ( fs >= w ) continue;
if ( ss >= h ) continue;
if ( fs < 0 ) continue;
if ( ss < 0 ) continue;
mask[fs + ss*w] = 1;
}
}
}
return mask;
}
/* Returns non-zero if peak has been vetoed.
* i.e. don't use result if return value is not zero. */
int integrate_peak(struct image *image, int cfs, int css,
double *pfs, double *pss, double *intensity, double *sigma,
double ir_inn, double ir_mid, double ir_out)
static int integrate_peak(struct image *image, int cfs, int css,
double *pfs, double *pss,
double *intensity, double *sigma,
double ir_inn, double ir_mid, double ir_out,
int use_max_adu)
{
signed int fs, ss;
signed int dfs, dss;
double lim_sq, out_lim_sq, mid_lim_sq;
double pk_total;
int pk_counts;
......@@ -164,11 +227,21 @@ int integrate_peak(struct image *image, int cfs, int css,
double bg_tot_sq = 0.0;
double var;
double aduph;
int *bgPkMask;
int p_cfs, p_css, p_w, p_h;
p = find_panel(image->det, cfs, css);
if ( p == NULL ) return 1;
if ( p->no_index ) return 1;
/* Determine regions where there is expected to be a peak */
p_cfs = cfs - p->min_fs;
p_css = css - p->min_ss; /* Panel-relative coordinates */
p_w = p->max_fs - p->min_fs + 1;
p_h = p->max_ss - p->min_ss + 1;
bgPkMask = make_BgMask(image, p, ir_out, p_cfs, p_css, ir_inn);
if ( bgPkMask == NULL ) return 1;
aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda);
lim_sq = pow(ir_inn, 2.0);
......@@ -176,23 +249,27 @@ int integrate_peak(struct image *image, int cfs, int css,
out_lim_sq = pow(ir_out, 2.0);
/* Estimate the background */
for ( fs=-ir_out; fs<=+ir_out; fs++ ) {
for ( ss=-ir_out; ss<=+ir_out; ss++ ) {
for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) {
for ( dss=-ir_out; dss<=+ir_out; dss++ ) {
double val;
uint16_t flags;
struct panel *p2;
int idx;
/* Restrict to annulus */
if ( fs*fs + ss*ss > out_lim_sq ) continue;
if ( fs*fs + ss*ss < mid_lim_sq ) continue;
if ( dfs*dfs + dss*dss > out_lim_sq ) continue;
if ( dfs*dfs + dss*dss < mid_lim_sq ) continue;
/* Strayed off one panel? */
p2 = find_panel(image->det, fs+cfs, ss+css);
if ( p2 != p ) return 1;
if ( p_cfs+dfs >= p_w ) continue;
if ( p_css+dss >= p_h ) continue;
if ( p_cfs+dfs < 0 ) continue;
if ( p_css+dss < 0 ) continue;
/* Check if there is a peak in the background region */
if ( bgPkMask[(p_cfs+dfs) + p_w*(p_css+dss)] ) continue;
idx = fs+cfs+image->width*(ss+css);
idx = dfs+cfs+image->width*(dss+css);
/* Veto this peak if we tried to integrate in a bad region */
if ( image->flags != NULL ) {
......@@ -211,7 +288,7 @@ int integrate_peak(struct image *image, int cfs, int css,
val = image->data[idx];
/* Veto peak if it contains saturation in bg region */
if ( val > p->max_adu ) return 1;
if ( use_max_adu && (val > p->max_adu) ) return 1;
bg_tot += val;
bg_tot_sq += pow(val, 2.0);
......@@ -228,22 +305,23 @@ int integrate_peak(struct image *image, int cfs, int css,
pk_total = 0.0;
pk_counts = 0;
fsct = 0.0; ssct = 0.0;
for ( fs=-ir_inn; fs<=+ir_inn; fs++ ) {
for ( ss=-ir_inn; ss<=+ir_inn; ss++ ) {
for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) {
for ( dss=-ir_inn; dss<=+ir_inn; dss++ ) {
double val;
uint16_t flags;
struct panel *p2;
int idx;
/* Inner mask radius */
if ( fs*fs + ss*ss > lim_sq ) continue;
if ( dfs*dfs + dss*dss > lim_sq ) continue;
/* Strayed off one panel? */
p2 = find_panel(image->det, fs+cfs, ss+css);
if ( p2 != p ) return 1;
if ( p_cfs+dfs >= p_w ) continue;
if ( p_css+dss >= p_h ) continue;
if ( p_cfs+dfs < 0 ) continue;
if ( p_css+dss < 0 ) continue;
idx = fs+cfs+image->width*(ss+css);
idx = dfs+cfs+image->width*(dss+css);
/* Veto this peak if we tried to integrate in a bad region */
if ( image->flags != NULL ) {
......@@ -262,13 +340,13 @@ int integrate_peak(struct image *image, int cfs, int css,
val = image->data[idx] - bg_mean;
/* Veto peak if it contains saturation */
if ( image->data[idx] > p->max_adu ) return 1;
if ( use_max_adu && (image->data[idx] > p->max_adu) ) return 1;
pk_counts++;
pk_total += val;
fsct += val*(cfs+fs);
ssct += val*(css+ss);
fsct += val*(cfs+dfs);
ssct += val*(css+dss);
}
}
......@@ -285,6 +363,8 @@ int integrate_peak(struct image *image, int cfs, int css,
if ( intensity != NULL ) *intensity = pk_total;
if ( sigma != NULL ) *sigma = sqrt(var);
free(bgPkMask);
return 0;
}
......@@ -309,7 +389,6 @@ static void search_peaks_in_panel(struct image *image, float threshold,
int nrej_snr = 0;
int nacc = 0;
int ncull;
const int pws = p->peak_sep/2;
data = image->data;
stride = image->width;
......@@ -352,16 +431,14 @@ static void search_peaks_in_panel(struct image *image, float threshold,
max = data[mask_fs+stride*mask_ss];
did_something = 0;
for ( s_ss=biggest(mask_ss-pws/2,
p->min_ss);
s_ss<=smallest(mask_ss+pws/2,
p->max_ss);
s_ss++ ) {
for ( s_fs=biggest(mask_fs-pws/2,
p->min_fs);
s_fs<=smallest(mask_fs+pws/2,
p->max_fs);
s_fs++ ) {
for ( s_ss=biggest(mask_ss-ir_inn, p->min_ss);
s_ss<=smallest(mask_ss+ir_inn, p->max_ss);
s_ss++ )
{
for ( s_fs=biggest(mask_fs-ir_inn, p->min_fs);
s_fs<=smallest(mask_fs+ir_inn, p->max_fs);
s_fs++ )
{
if ( data[s_fs+stride*s_ss] > max ) {
max = data[s_fs+stride*s_ss];
......@@ -374,8 +451,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
}
/* Abort if drifted too far from the foot point */
if ( distance(mask_fs, mask_ss, fs, ss) >
p->peak_sep/2.0 )
if ( distance(mask_fs, mask_ss, fs, ss) > ir_inn )
{
break;
}
......@@ -383,7 +459,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
} while ( did_something );
/* Too far from foot point? */
if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep/2.0 ) {
if ( distance(mask_fs, mask_ss, fs, ss) > ir_inn ) {
nrej_dis++;
continue;
}
......@@ -397,7 +473,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
/* Centroid peak and get better coordinates. */
r = integrate_peak(image, mask_fs, mask_ss,
&f_fs, &f_ss, &intensity, &sigma,
ir_inn, ir_mid, ir_out);
ir_inn, ir_mid, ir_out, 0);
if ( r ) {
/* Bad region - don't detect peak */
......@@ -419,7 +495,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
/* Check for a nearby feature */
image_feature_closest(image->features, f_fs, f_ss, &d, &idx);
if ( d < p->peak_sep/2.0 ) {
if ( d < 2.0*ir_inn ) {
nrej_pro++;
continue;
}
......@@ -445,6 +521,12 @@ static void search_peaks_in_panel(struct image *image, float threshold,
// "%i in bad regions, %i with SNR < %g, %i badrow culled.\n",
// nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
// nrej_snr, min_snr, ncull);
if ( ncull != 0 ) {
STATUS("WARNING: %i peaks were badrow culled. This feature"
" should not usually be used.\nConsider setting"
" badrow=- in the geometry file.\n", ncull);
}
}
......@@ -655,7 +737,11 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
}
r = integrate_peak(image, pfs, pss, &fs, &ss,
&intensity, &sigma, ir_inn, ir_mid, ir_out);
&intensity, &sigma, ir_inn, ir_mid, ir_out,
1);
/* I/sigma(I) cutoff */
if ( intensity/sigma < min_snr ) r = 1;
/* Record intensity and set redundancy to 1 on success */
if ( r == 0 ) {
......
......@@ -900,6 +900,7 @@ struct _reflistiterator {
**/
Reflection *first_refl(RefList *list, RefListIterator **piter)
{
Reflection *refl;
RefListIterator *iter;
iter = malloc(sizeof(struct _reflistiterator));
......@@ -908,7 +909,9 @@ Reflection *first_refl(RefList *list, RefListIterator **piter)
iter->stack_ptr = 0;
*piter = iter;
Reflection *refl = list->head;
if ( list == NULL ) return NULL;
refl = list->head;
do {
......
......@@ -474,6 +474,7 @@ int skip_some_files(FILE *fh, int n)
rval = fgets(line, 1023, fh);
if ( rval == NULL ) continue;
chomp(line);
if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
} while ( rval != NULL );
......
......@@ -167,57 +167,13 @@ static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val)
}
void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
int na, int nb, int nc, UnitCell *ucell)
static void do_panels(struct gpu_context *gctx, struct image *image,
int *n_inf, int *n_neg, int *n_nan, double k, double frac,
double xo, double yo)
{
cl_int err;
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
int i;
cl_float16 cell;
cl_int4 ncells;
int n_inf = 0;
int n_neg = 0;
int n_nan = 0;
if ( gctx == NULL ) {
ERROR("GPU setup failed.\n");
return;
}
cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az;
cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz;
cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz;
ncells.s[0] = na;
ncells.s[1] = nb;
ncells.s[2] = nc;
ncells.s[3] = 0; /* unused */
/* Ensure all required LUTs are available */
check_sinc_lut(gctx, na);
check_sinc_lut(gctx, nb);
check_sinc_lut(gctx, nc);
if ( set_arg_float(gctx, 2, 1.0/image->lambda) ) return;
if ( set_arg_mem(gctx, 9, gctx->intensities) ) return;
if ( set_arg_mem(gctx, 10, gctx->sinc_luts[na-1]) ) return;
if ( set_arg_mem(gctx, 11, gctx->sinc_luts[nb-1]) ) return;
if ( set_arg_mem(gctx, 12, gctx->sinc_luts[nc-1]) ) return;
if ( set_arg_mem(gctx, 13, gctx->flags) ) return;
/* Unit cell */
err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell);
if ( err != CL_SUCCESS ) {
ERROR("Couldn't set unit cell: %s\n", clError(err));
return;
}
/* Allocate memory for the result */
image->data = calloc(image->width * image->height, sizeof(float));
image->twotheta = calloc(image->width * image->height, sizeof(double));
if ( set_arg_float(gctx, 1, k) ) return;
/* Iterate over panels */
for ( i=0; i<image->det->n_panels; i++ ) {
......@@ -225,14 +181,12 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
size_t dims[2];
size_t ldims[2] = {1, 1};
struct panel *p;
cl_mem tt;
size_t tt_size;
cl_mem diff;
size_t diff_size;
float *diff_ptr;
float *tt_ptr;
int pan_width, pan_height;
int fs, ss;
cl_int err;
p = &image->det->panels[i];
......@@ -247,25 +201,19 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
ERROR("Couldn't allocate diffraction memory\n");
return;
}
tt_size = pan_width * pan_height * sizeof(cl_float);
tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size,
NULL, &err);
if ( err != CL_SUCCESS ) {
ERROR("Couldn't allocate twotheta memory\n");
return;
}
if ( set_arg_mem(gctx, 0, diff) ) return;
if ( set_arg_mem(gctx, 1, tt) ) return;
if ( set_arg_int(gctx, 3, pan_width) ) return;
if ( set_arg_float(gctx, 4, p->cnx) ) return;
if ( set_arg_float(gctx, 5, p->cny) ) return;
if ( set_arg_float(gctx, 6, p->res) ) return;
if ( set_arg_float(gctx, 7, p->clen) ) return;
if ( set_arg_float(gctx, 14, p->fsx) ) return;
if ( set_arg_float(gctx, 15, p->fsy) ) return;
if ( set_arg_float(gctx, 16, p->ssx) ) return;
if ( set_arg_float(gctx, 17, p->ssy) ) return;
if ( set_arg_int(gctx, 2, pan_width) ) return;
if ( set_arg_float(gctx, 3, p->cnx) ) return;
if ( set_arg_float(gctx, 4, p->cny) ) return;
if ( set_arg_float(gctx, 5, p->res) ) return;
if ( set_arg_float(gctx, 6, p->clen) ) return;
if ( set_arg_float(gctx, 13, p->fsx) ) return;
if ( set_arg_float(gctx, 14, p->fsy) ) return;
if ( set_arg_float(gctx, 15, p->ssx) ) return;
if ( set_arg_float(gctx, 16, p->ssy) ) return;
if ( set_arg_float(gctx, 17, xo) ) return;
if ( set_arg_float(gctx, 18, yo) ) return;
dims[0] = pan_width;
dims[1] = pan_height;
......@@ -288,43 +236,124 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
clError(err));
return;
}
tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ,
0, tt_size, 0, NULL, NULL, &err);
if ( err != CL_SUCCESS ) {
ERROR("Couldn't map tt buffer\n");
return;
}
for ( fs=0; fs<pan_width; fs++ ) {
for ( ss=0; ss<pan_height; ss++ ) {
float val, tt;