From 15bced1f814a8d1d50bc23c3a43621d8fd114bbf Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Oct 2015 12:01:44 +0200 Subject: Remove use of image->data in record_image() --- libcrystfel/src/detector.c | 115 ++++++++++++++++++++++++--------------------- 1 file changed, 61 insertions(+), 54 deletions(-) diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index ae70e406..ddb05f3a 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -412,33 +412,17 @@ int detector_has_clen_references(struct detector *det) } -void record_image(struct image *image, int do_poisson, double background, - gsl_rng *rng, double beam_radius, double nphotons) +static void record_panel(struct panel *p, float *dp, int do_poisson, + gsl_rng *rng, double ph_per_e, double background, + double lambda, + int *n_neg1, int *n_inf1, int *n_nan1, + int *n_neg2, int *n_inf2, int *n_nan2, + double *max_tt) { - int x, y; - double total_energy, energy_density; - double ph_per_e; - double area; - double max_tt = 0.0; - int n_inf1 = 0; - int n_neg1 = 0; - int n_nan1 = 0; - int n_inf2 = 0; - int n_neg2 = 0; - int n_nan2 = 0; - - /* How many photons are scattered per electron? */ - area = M_PI*pow(beam_radius, 2.0); - total_energy = nphotons * ph_lambda_to_en(image->lambda); - energy_density = total_energy / area; - ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0); - STATUS("Fluence = %8.2e photons, " - "Energy density = %5.3f kJ/cm^2, " - "Total energy = %5.3f microJ\n", - nphotons, energy_density/1e7, total_energy*1e6); + int fs, ss; - for ( x=0; xwidth; x++ ) { - for ( y=0; yheight; y++ ) { + for ( ss=0; ss>p->h; ss++ ) { + for ( fs=0; fs>p->w; fs++ ) { double counts; double cf; @@ -447,28 +431,27 @@ void record_image(struct image *image, int do_poisson, double background, double xs, ys, rx, ry; double dsq, proj_area; float dval; - struct panel *p; + double twotheta; - intensity = (double)image->data[x + image->width*y]; - if ( isinf(intensity) ) n_inf1++; - if ( intensity < 0.0 ) n_neg1++; - if ( isnan(intensity) ) n_nan1++; - - p = find_panel(image->det, x, y); + intensity = (double)dp[fs + p->w*ss]; + if ( isinf(intensity) ) (*n_inf1)++; + if ( intensity < 0.0 ) (*n_neg1)++; + if ( isnan(intensity) ) (*n_nan1)++; /* 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) */ - proj_area = pix_area * cos(image->twotheta[x + image->width*y]); - /* Calculate distance from crystal to pixel */ - xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx; - ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy; + xs = fs*p->fsx + ss*p->ssx; + ys = ss*p->fsy + ss*p->ssy; rx = (xs + p->cnx) / p->res; ry = (ys + p->cny) / p->res; dsq = pow(rx, 2.0) + pow(ry, 2.0); + twotheta = atan2(sqrt(dsq), p->clen); + + /* Area of pixel as seen from crystal (approximate) */ + proj_area = pix_area * cos(twotheta); /* Projected area of pixel divided by distance squared */ sa = proj_area / (dsq + Lsq); @@ -484,36 +467,60 @@ void record_image(struct image *image, int do_poisson, double background, dval = counts + poisson_noise(rng, background); /* Convert to ADU */ - dval *= p->adu_per_eV * ph_lambda_to_eV(image->lambda); + dval *= p->adu_per_eV * ph_lambda_to_eV(lambda); /* Saturation */ if ( dval > p->max_adu ) dval = p->max_adu; - image->data[x + image->width*y] = dval; + dp[fs + p->w*ss] = dval; /* Sanity checks */ - if ( isinf(image->data[x+image->width*y]) ) n_inf2++; - if ( isnan(image->data[x+image->width*y]) ) n_nan2++; - if ( image->data[x+image->width*y] < 0.0 ) n_neg2++; + if ( isinf(dp[fs + p->w*ss]) ) n_inf2++; + if ( isnan(dp[fs + p->w*ss]) ) n_nan2++; + if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++; - if ( image->twotheta[x + image->width*y] > max_tt ) { - max_tt = image->twotheta[x + image->width*y]; - } + if ( twotheta > *max_tt ) *max_tt = twotheta; } - progress_bar(x, image->width-1, "Post-processing"); } +} - STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n", - rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); - double tt_side = image->twotheta[(image->width/2)+image->width*0]; - STATUS("At middle of bottom edge: %.2f deg, min d = %.2f nm\n", - rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9); +void record_image(struct image *image, int do_poisson, double background, + gsl_rng *rng, double beam_radius, double nphotons) +{ + double total_energy, energy_density; + double ph_per_e; + double area; + double max_tt = 0.0; + int n_inf1 = 0; + int n_neg1 = 0; + int n_nan1 = 0; + int n_inf2 = 0; + int n_neg2 = 0; + int n_nan2 = 0; + int pn; - tt_side = image->twotheta[0+image->width*(image->height/2)]; - STATUS("At middle of left edge: %.2f deg, min d = %.2f nm\n", - rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9); + /* How many photons are scattered per electron? */ + area = M_PI*pow(beam_radius, 2.0); + total_energy = nphotons * ph_lambda_to_en(image->lambda); + energy_density = total_energy / area; + ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0); + STATUS("Fluence = %8.2e photons, " + "Energy density = %5.3f kJ/cm^2, " + "Total energy = %5.3f microJ\n", + nphotons, energy_density/1e7, total_energy*1e6); + + for ( pn=0; pndet->n_panels; pn++ ) { + record_panel(&image->det->panels[pn], image->dp[pn], + do_poisson, rng, ph_per_e, background, + image->lambda, + &n_neg1, &n_inf1, &n_nan1, + &n_neg2, &n_inf2, &n_nan2, &max_tt); + } + + STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n", + rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); STATUS("Halve the d values to get the voxel size for a synthesis.\n"); -- cgit v1.2.3