diff options
-rw-r--r-- | libcrystfel/src/detector.c | 136 | ||||
-rw-r--r-- | libcrystfel/src/detector.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/hdf5-file.c | 136 | ||||
-rw-r--r-- | libcrystfel/src/image.h | 9 | ||||
-rw-r--r-- | libcrystfel/src/peaks.c | 60 | ||||
-rw-r--r-- | src/diffraction-gpu.c | 36 | ||||
-rw-r--r-- | src/diffraction.c | 62 | ||||
-rw-r--r-- | src/dw-hdfsee.c | 10 | ||||
-rw-r--r-- | src/geoptimiser.c | 110 | ||||
-rw-r--r-- | src/hdfsee-render.c | 11 | ||||
-rw-r--r-- | src/partial_sim.c | 59 | ||||
-rw-r--r-- | src/pattern_sim.c | 84 | ||||
-rw-r--r-- | src/process_image.c | 4 | ||||
-rw-r--r-- | tests/gpu_sim_check.c | 80 | ||||
-rw-r--r-- | tests/integration_check.c | 1 | ||||
-rw-r--r-- | tests/prof2d_check.c | 1 | ||||
-rw-r--r-- | tests/ring_check.c | 33 |
17 files changed, 473 insertions, 361 deletions
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 99526ff9..53510d69 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; x<image->width; x++ ) { - for ( y=0; y<image->height; 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; - - intensity = (double)image->data[x + image->width*y]; - if ( isinf(intensity) ) n_inf1++; - if ( intensity < 0.0 ) n_neg1++; - if ( isnan(intensity) ) n_nan1++; + double twotheta; - 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; pn<image->det->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"); @@ -554,9 +561,7 @@ struct panel *find_panel(struct detector *det, double fs, double ss) } -/* Like find_panel(), but uses the original panel bounds, i.e. referring to - * what's in the HDF5 file */ -struct panel *find_orig_panel(struct detector *det, double fs, double ss) +signed int find_orig_panel_number(struct detector *det, double fs, double ss) { int p; @@ -564,13 +569,20 @@ struct panel *find_orig_panel(struct detector *det, double fs, double ss) if ( (fs >= det->panels[p].orig_min_fs) && (fs < det->panels[p].orig_max_fs+1) && (ss >= det->panels[p].orig_min_ss) - && (ss < det->panels[p].orig_max_ss+1) ) - { - return &det->panels[p]; - } + && (ss < det->panels[p].orig_max_ss+1) ) return p; } - return NULL; + return -1; +} + + +/* Like find_panel(), but uses the original panel bounds, i.e. referring to + * what's in the HDF5 file */ +struct panel *find_orig_panel(struct detector *det, double fs, double ss) +{ + signed int pn = find_orig_panel_number(det, fs, ss); + if ( pn == -1 ) return NULL; + return &det->panels[pn]; } diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index b3474b7a..5e544b75 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -205,6 +205,8 @@ extern struct panel *find_panel(struct detector *det, double fs, double ss); extern signed int find_panel_number(struct detector *det, double fs, double ss); extern struct panel *find_orig_panel(struct detector *det, double fs, double ss); +extern signed int find_orig_panel_number(struct detector *det, + double fs, double ss); extern struct detector *get_detector_geometry(const char *filename, struct beam_params *beam); diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index e7b6eed4..b159327e 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -748,7 +748,7 @@ static struct hdf5_write_location *make_location_list(struct detector *det, } -static void write_location(hid_t fh, const struct image *image, +static void write_location(hid_t fh, struct detector *det, float *data, struct hdf5_write_location *loc) { hid_t sh, dh, ph; @@ -784,7 +784,7 @@ static void write_location(hid_t fh, const struct image *image, struct panel p; int r; - p = image->det->panels[loc->panel_idxs[pi]]; + p = det->panels[loc->panel_idxs[pi]]; f_offset[0] = p.orig_min_ss; f_offset[1] = p.orig_min_fs; @@ -810,14 +810,14 @@ static void write_location(hid_t fh, const struct image *image, m_count[0] = p.max_ss - p.min_ss +1; m_count[1] = p.max_fs - p.min_fs +1; - dimsm[0] = image->height; - dimsm[1] = image->width; + dimsm[0] = det->max_fs + 1; + dimsm[1] = det->max_ss + 1; memspace = H5Screate_simple(2, dimsm, NULL); r = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, m_offset, NULL, m_count, NULL); r = H5Dwrite(dh, H5T_NATIVE_FLOAT, memspace, - dh_dataspace, H5P_DEFAULT, image->data); + dh_dataspace, H5P_DEFAULT, data); if ( r < 0 ) { ERROR("Couldn't write data\n"); H5Pclose(ph); @@ -948,6 +948,34 @@ static void write_spectrum(hid_t fh, struct sample *spectrum, int spectrum_size, } +static float *make_array_from_dp(const struct image *image) +{ + int i; + float *data; + + data = malloc(image->width * image->height * sizeof(float)); + if ( data == NULL ) { + ERROR("Failed to allocate data\n"); + return NULL; + } + + for ( i=0; i<image->det->n_panels; i++ ) { + + int fs, ss; + struct panel *p = &image->det->panels[i]; + + for ( ss=0; ss<p->h; ss++ ) { + for ( fs=0; fs<p->w; fs++ ) { + int idx = p->min_fs+fs + image->width*(p->min_ss+ss); + data[idx] = image->dp[i][fs + p->w*ss]; + } + } + } + + return data; +} + + int hdf5_write_image(const char *filename, const struct image *image, char *element) { @@ -957,12 +985,16 @@ int hdf5_write_image(const char *filename, const struct image *image, struct hdf5_write_location *locations; int num_locations; const char *ph_en_loc; + float *data; if ( image->det == NULL ) { ERROR("Geometry not available\n"); return 1; } + data = make_array_from_dp(image); + if ( data == NULL ) return 1; + fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); if ( fh < 0 ) { ERROR("Couldn't create file: %s\n", filename); @@ -979,7 +1011,7 @@ int hdf5_write_image(const char *filename, const struct image *image, &num_locations); for ( li=0; li<num_locations; li++ ) { - write_location(fh, image, &locations[li]); + write_location(fh, image->det, data, &locations[li]); } if ( image->beam == NULL || image->beam->photon_energy_from == NULL ) { @@ -1001,7 +1033,7 @@ int hdf5_write_image(const char *filename, const struct image *image, for ( li=0; li<num_locations; li ++ ) { free(locations[li].panel_idxs); } - + free(data); free(locations); return 0; } @@ -1062,18 +1094,28 @@ static void debodge_saturation(struct hdfile *f, struct image *image) for ( i=0; i<size[0]; i++ ) { - unsigned int x, y; + unsigned int fs, ss; float val; + struct panel *p; + signed int pn; - x = buf[3*i+0]; - y = buf[3*i+1]; + fs = buf[3*i+0]; + ss = buf[3*i+1]; val = buf[3*i+2]; - image->data[x+image->width*y] = val / 5.0; - image->data[x+1+image->width*y] = val / 5.0; - image->data[x-1+image->width*y] = val / 5.0; - image->data[x+image->width*(y+1)] = val / 5.0; - image->data[x+image->width*(y-1)] = val / 5.0; + /* Turn "original" position into "panel" position */ + pn = find_orig_panel_number(image->det, fs, ss); + if ( pn == -1 ) { + ERROR("Failed to find panel!\n"); + continue; + } + p = &image->det->panels[pn]; + + image->dp[pn][fs+p->w*ss] = val/5.0; + image->dp[pn][fs+1+p->w*ss] = val/5.0; + image->dp[pn][fs-1+p->w*ss] = val/5.0; + image->dp[pn][fs+p->w*(ss-1)] = val/5.0; + image->dp[pn][fs+p->w*(ss+1)] = val/5.0; } @@ -1083,7 +1125,8 @@ static void debodge_saturation(struct hdfile *f, struct image *image) } -static int unpack_panels(struct image *image, struct detector *det) +static int unpack_panels(struct image *image, struct detector *det, + float *data, uint16_t *flags) { int pi; @@ -1118,7 +1161,7 @@ static int unpack_panels(struct image *image, struct detector *det) css = ss+p->min_ss; idx = cfs + css*image->width; - image->dp[pi][fs+p->w*ss] = image->data[idx]; + image->dp[pi][fs+p->w*ss] = data[idx]; if ( p->no_index ) bad = 1; @@ -1126,18 +1169,18 @@ static int unpack_panels(struct image *image, struct detector *det) bad = 1; } - if ( image->flags != NULL ) { + if ( flags != NULL ) { - int flags; + int f; - flags = image->flags[idx]; + f = flags[idx]; /* Bad if it's missing any of the "good" bits */ - if ( (flags & image->det->mask_good) + if ( (f & image->det->mask_good) != image->det->mask_good ) bad = 1; /* Bad if it has any of the "bad" bits. */ - if ( flags & image->det->mask_bad ) bad = 1; + if ( f & image->det->mask_bad ) bad = 1; } image->bad[pi][fs+p->w*ss] = bad; @@ -1403,17 +1446,15 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, free(buf); return 1; } - image->data = buf; if ( image->det != NULL ) { ERROR("WARNING: hdf5_read() called with geometry structure.\n"); } image->det = simple_geometry(image); + unpack_panels(image, image->det, buf, NULL); if ( satcorr ) debodge_saturation(f, image); - unpack_panels(image, image->det); - if ( image->beam != NULL ) { fill_in_beam_parameters(image->beam, f, NULL, image); @@ -1432,12 +1473,12 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, } -static void load_mask(struct hdfile *f, struct event *ev, char *mask, - const char *mask_file, - const char *pname, struct image *image, - size_t p_w, size_t sum_p_h, - hsize_t *f_offset, hsize_t *f_count, - hsize_t *m_offset, hsize_t *m_count) +static int load_mask(struct hdfile *f, struct event *ev, char *mask, + const char *mask_file, + const char *pname, struct image *image, + size_t p_w, size_t sum_p_h, uint16_t *flags, + hsize_t *f_offset, hsize_t *f_count, + hsize_t *m_offset, hsize_t *m_count) { hid_t mask_dataspace, mask_dh; int exists; @@ -1450,7 +1491,7 @@ static void load_mask(struct hdfile *f, struct event *ev, char *mask, fh = H5Fopen(mask_file, H5F_ACC_RDONLY, H5P_DEFAULT); if ( fh < 0 ) { ERROR("Couldn't open mask file '%s'\n", mask_file); - return; + return 1; } } else { fh = f->fh; @@ -1491,7 +1532,7 @@ static void load_mask(struct hdfile *f, struct event *ev, char *mask, } r = H5Dread(mask_dh, H5T_NATIVE_UINT16, memspace, - mask_dataspace, H5P_DEFAULT, image->flags); + mask_dataspace, H5P_DEFAULT, flags); if ( r < 0 ) { ERROR("Couldn't read flags for panel %s\n", pname); goto err; @@ -1501,14 +1542,12 @@ static void load_mask(struct hdfile *f, struct event *ev, char *mask, H5Dclose(mask_dh); if ( ev != NULL ) free(mask); - return; + return 0; err: if ( mask_file != NULL ) H5Fclose(fh); if ( ev != NULL ) free(mask); - free(image->flags); - image->flags = NULL; - return; + return 1; } @@ -1517,6 +1556,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, { herr_t r; float *buf; + uint16_t *flags; int sum_p_h; int p_w; int pi; @@ -1548,8 +1588,8 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, image->width = p_w; image->height = sum_p_h; - image->flags = calloc(p_w*sum_p_h,sizeof(uint16_t)); - if ( image->flags == NULL ) { + flags = calloc(p_w*sum_p_h,sizeof(uint16_t)); + if ( flags == NULL ) { ERROR("Failed to allocate memory for flags\n"); return 1; } @@ -1695,9 +1735,11 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, H5Sclose(memspace); if ( p->mask != NULL ) { - load_mask(f, ev, p->mask, p->mask_file, p->name, - image, p_w, sum_p_h, - f_offset, f_count, m_offset, m_count); + if ( load_mask(f, ev, p->mask, p->mask_file, p->name, + image, p_w, sum_p_h, flags, + f_offset, f_count, m_offset, m_count) ) { + ERROR("Error loading bad pixel mask!\n"); + } } free(f_offset); @@ -1705,13 +1747,10 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, } - image->data = buf; - - if ( satcorr ) debodge_saturation(f, image); - fill_in_values(image->det, f, ev); - unpack_panels(image, image->det); + unpack_panels(image, image->det, buf, flags); + if ( satcorr ) debodge_saturation(f, image); if ( image->beam != NULL ) { @@ -1727,6 +1766,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, } + free(buf); + free(flags); + return 0; } diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 10e83905..52614d06 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -113,10 +113,6 @@ struct beam_params * <programlisting> * struct image * { - * float *data; - * uint16_t *flags; - * double *twotheta; - * * Crystal **crystals; * int n_crystals; * IndexingMethod indexed_by; @@ -163,11 +159,6 @@ struct image; struct image { - /* The following three fields will be going away in the future */ - float *data; - uint16_t *flags; - double *twotheta; - float **dp; /* Data in panel */ int **bad; /* Bad pixels by panel */ diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 9e3c5e77..69aad4c0 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -242,7 +242,7 @@ static int integrate_peak(struct image *image, int cfs, int css, double bg_tot_sq = 0.0; double var; double aduph; - int p_cfs, p_css, p_w, p_h; + int p_cfs, p_css; signed int pn; pn = find_panel_number(image->det, cfs, css); @@ -254,8 +254,6 @@ static int integrate_peak(struct image *image, int cfs, int css, /* 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; aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda); @@ -264,8 +262,8 @@ static int integrate_peak(struct image *image, int cfs, int css, out_lim_sq = pow(ir_out, 2.0); /* Estimate the background */ - for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) { for ( dss=-ir_out; dss<=+ir_out; dss++ ) { + for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) { double val; int idx; @@ -275,7 +273,7 @@ static int integrate_peak(struct image *image, int cfs, int css, if ( dfs*dfs + dss*dss < mid_lim_sq ) continue; /* Strayed off one panel? */ - if ( (p_cfs+dfs >= p_w) || (p_css+dss >= p_h) + if ( (p_cfs+dfs >= p->w) || (p_css+dss >= p->h) || (p_cfs+dfs < 0 ) || (p_css+dss < 0) ) return 4; /* Wandered into a bad region? */ @@ -283,8 +281,8 @@ static int integrate_peak(struct image *image, int cfs, int css, return 14; } - idx = dfs+cfs+image->width*(dss+css); - val = image->data[idx]; + idx = dfs+p_cfs+p->w*(dss+p_css); + val = image->dp[pn][idx]; /* Check if peak contains saturation in bg region */ if ( (saturated != NULL) && (val > p->max_adu) ) *saturated = 1; @@ -304,8 +302,8 @@ static int integrate_peak(struct image *image, int cfs, int css, pk_total = 0.0; pk_counts = 0; fsct = 0.0; ssct = 0.0; - for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) { for ( dss=-ir_inn; dss<=+ir_inn; dss++ ) { + for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) { double val; int idx; @@ -314,7 +312,7 @@ static int integrate_peak(struct image *image, int cfs, int css, if ( dfs*dfs + dss*dss > lim_sq ) continue; /* Strayed off one panel? */ - if ( (p_cfs+dfs >= p_w) || (p_css+dss >= p_h) + if ( (p_cfs+dfs >= p->w) || (p_css+dss >= p->h) || (p_cfs+dfs < 0 ) || (p_css+dss < 0) ) return 8; /* Wandered into a bad region? */ @@ -322,8 +320,8 @@ static int integrate_peak(struct image *image, int cfs, int css, return 15; } - idx = dfs+cfs+image->width*(dss+css); - val = image->data[idx]; + idx = dfs+p_cfs+p->w*(dss+p_css); + val = image->dp[pn][idx]; /* Check if peak contains saturation */ if ( (saturated != NULL) && (val > p->max_adu) ) *saturated = 1; @@ -356,13 +354,13 @@ static int integrate_peak(struct image *image, int cfs, int css, static void search_peaks_in_panel(struct image *image, float threshold, - float min_gradient, float min_snr, - struct panel *p, + float min_gradient, float min_snr, int pn, double ir_inn, double ir_mid, double ir_out, int use_saturated) { int fs, ss, stride; float *data; + struct panel *p; double d; int idx; double f_fs = 0.0; @@ -378,11 +376,12 @@ static void search_peaks_in_panel(struct image *image, float threshold, int nacc = 0; int ncull; - data = image->data; - stride = image->width; + p = &image->det->panels[pn]; + data = image->dp[pn]; + stride = p->w; - for ( fs = p->min_fs+1; fs <= p->max_fs-1; fs++ ) { - for ( ss = p->min_ss+1; ss <= p->max_ss-1; ss++ ) { + for ( ss=0; ss<=p->h; ss++ ) { + for ( fs=0; fs<=p->w; fs++ ) { double dx1, dx2, dy1, dy2; double dxs, dys; @@ -425,12 +424,12 @@ 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-ir_inn, p->min_ss); - s_ss<=smallest(mask_ss+ir_inn, p->max_ss); + for ( s_ss=biggest(mask_ss-ir_inn, 0); + s_ss<=smallest(mask_ss+ir_inn, p->h-1); s_ss++ ) { - for ( s_fs=biggest(mask_fs-ir_inn, p->min_fs); - s_fs<=smallest(mask_fs+ir_inn, p->max_fs); + for ( s_fs=biggest(mask_fs-ir_inn, 0); + s_fs<=smallest(mask_fs+ir_inn, p->w-1); s_fs++ ) { @@ -459,13 +458,13 @@ static void search_peaks_in_panel(struct image *image, float threshold, } /* Should be enforced by bounds used above. Muppet check. */ - assert(mask_fs <= p->max_fs); - assert(mask_ss <= p->max_ss); - assert(mask_fs >= p->min_fs); - assert(mask_ss >= p->min_ss); + assert(mask_fs <= p->w); + assert(mask_ss <= p->h); + assert(mask_fs >= 0); + assert(mask_ss >= 0); /* Centroid peak and get better coordinates. */ - r = integrate_peak(image, mask_fs, mask_ss, + r = integrate_peak(image, mask_fs+p->min_fs, mask_ss+p->min_ss, &f_fs, &f_ss, &intensity, &sigma, ir_inn, ir_mid, ir_out, &saturated); @@ -504,8 +503,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, } /* Add using "better" coordinates */ - image_add_feature(image->features, f_fs, f_ss, image, intensity, - NULL); + image_add_feature(image->features, f_fs, f_ss, + image, intensity, NULL); nacc++; } @@ -551,11 +550,10 @@ void search_peaks(struct image *image, float threshold, float min_gradient, for ( i=0; i<image->det->n_panels; i++ ) { - struct panel *p = &image->det->panels[i]; + if ( image->det->panels[i].no_index ) continue; - if ( p->no_index ) continue; search_peaks_in_panel(image, threshold, min_gradient, - min_snr, p, ir_inn, ir_mid, ir_out, + min_snr, i, ir_inn, ir_mid, ir_out, use_saturated); } diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c index bdfdb13b..1f937bf8 100644 --- a/src/diffraction-gpu.c +++ b/src/diffraction-gpu.c @@ -3,11 +3,11 @@ * * Calculate diffraction patterns by Fourier methods (GPU version) * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2009-2014 Thomas White <taw@physics.org> + * 2009-2015 Thomas White <taw@physics.org> * 2013 Alexandra Tolstikova * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de> * @@ -253,20 +253,17 @@ static int do_panels(struct gpu_context *gctx, struct image *image, return 1; } - for ( fs=0; fs<pan_width; fs++ ) { for ( ss=0; ss<pan_height; ss++ ) { + for ( fs=0; fs<pan_width; fs++ ) { float val; - int tfs, tss; val = diff_ptr[fs + pan_width*ss]; if ( isinf(val) ) (*n_inf)++; if ( val < 0.0 ) (*n_neg)++; if ( isnan(val) ) (*n_nan)++; - tfs = p->min_fs + fs; - tss = p->min_ss + ss; - image->data[tfs + image->width*tss] += val; + image->dp[i][fs + pan_width*ss] += val; } } @@ -294,7 +291,6 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int n_inf = 0; int n_neg = 0; int n_nan = 0; - int fs, ss; int i; if ( gctx == NULL ) { @@ -326,11 +322,19 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return 1; /* Allocate memory for the result */ - image->data = calloc(image->width * image->height, sizeof(float)); - if ( image->data == NULL ) { + image->dp = malloc(image->det->n_panels * sizeof(float *)); + if ( image->dp == NULL ) { ERROR("Couldn't allocate memory for result.\n"); return 1; } + for ( i=0; i<image->det->n_panels; i++ ) { + struct panel *p = &image->det->panels[i]; + image->dp[i] = calloc(p->w * p->h, sizeof(float)); + if ( image->dp[i] == NULL ) { + ERROR("Couldn't allocate memory for panel %i\n", i); + return 1; + } + } double tot = 0.0; for ( i=0; i<image->nsamples; i++ ) { @@ -356,18 +360,6 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, n_neg, n_inf, n_nan); } - /* Calculate "2theta" values for detector geometry */ - image->twotheta = calloc(image->width * image->height, sizeof(double)); - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; ss++ ) { - - double twotheta; - get_q(image, fs, ss, &twotheta, 1.0/image->lambda); - image->twotheta[fs + image->width*ss] = twotheta; - - } - } - return 0; } diff --git a/src/diffraction.c b/src/diffraction.c index 5a936809..93674f52 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -361,31 +361,30 @@ static double molecule_factor(const double *intensities, const double *phases, } - -static void diffraction_at_k(struct image *image, const double *intensities, - const double *phases, const unsigned char *flags, - UnitCell *cell, GradientMethod m, - const SymOpList *sym, double k, - double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - double *lut_a, double *lut_b, double *lut_c, - double weight) +static void diffraction_panel(struct image *image, const double *intensities, + const double *phases, const unsigned char *flags, + UnitCell *cell, GradientMethod m, + const SymOpList *sym, double k, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz, + double *lut_a, double *lut_b, double *lut_c, + int pn, double weight) { - unsigned int fs, ss; + int fs, ss; const int nxs = 4; const int nys = 4; + struct panel *p = &image->det->panels[pn]; weight /= nxs*nys; - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; ss++ ) { + for ( ss=0; ss<p->h; ss++ ) { + for ( fs=0; fs<p->w; fs++ ) { int idx; double f_lattice, I_lattice; double I_molecule; struct rvec q; - double twotheta; int xs, ys; float xo, yo; @@ -395,7 +394,7 @@ static void diffraction_at_k(struct image *image, const double *intensities, xo = (1.0/nxs) * xs; yo = (1.0/nys) * ys; - q = get_q(image, fs+xo, ss+yo, &twotheta, k); + q = get_q_for_panel(p, fs+xo, ss+yo, NULL, k); f_lattice = lattice_factor(q, ax, ay, az, bx, by, bz, @@ -411,14 +410,33 @@ static void diffraction_at_k(struct image *image, const double *intensities, I_lattice = pow(f_lattice, 2.0); - idx = fs + image->width*ss; - image->data[idx] += I_lattice * I_molecule * weight; - image->twotheta[idx] = twotheta; + idx = fs + p->w*ss; + image->dp[pn][idx] += I_lattice * I_molecule * weight; } } } - progress_bar(fs, image->width-1, "Calculating diffraction"); + progress_bar(ss, p->h-1, "Calculating diffraction"); + } +} + + +static void diffraction_at_k(struct image *image, const double *intensities, + const double *phases, const unsigned char *flags, + UnitCell *cell, GradientMethod m, + const SymOpList *sym, double k, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz, + double *lut_a, double *lut_b, double *lut_c, + double weight) +{ + int i; + + for ( i=0; i<image->det->n_panels; i++ ) { + diffraction_panel(image, intensities, phases, flags, cell, m, + sym, k, ax, ay, az, bx, by, bz, cx, cy, cz, + lut_a, lut_b, lut_c, i, weight); } } @@ -710,12 +728,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc, cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - /* Allocate (and zero) the "diffraction array" */ - image->data = calloc(image->width * image->height, sizeof(float)); - - /* Needed later for Lorentz calculation */ - image->twotheta = malloc(image->width * image->height * sizeof(double)); - lut_a = get_sinc_lut(na, no_fringes); lut_b = get_sinc_lut(nb, no_fringes); lut_c = get_sinc_lut(nc, no_fringes); diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c index 004c92ae..e260e4a0 100644 --- a/src/dw-hdfsee.c +++ b/src/dw-hdfsee.c @@ -94,8 +94,6 @@ static gint displaywindow_closed(GtkWidget *window, DisplayWindow *dw) if ( dw->image != NULL ) { free(dw->image->filename); - free(dw->image->data); - free(dw->image->flags); free(dw->image); } @@ -207,7 +205,7 @@ static int render_adsc_uint16(DisplayWindow *dw, const char *filename) fs = dfs; ss = dss; - val = image->data[fs + image->width * ss]; + val = 0.0;//image->data[fs + image->width * ss]; FIXME! if ( val < 0 ) { out = 0; } else if ( val > 65535 ) { @@ -925,8 +923,6 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event) if ( dw->not_ready_yet ) return 0; - float *old_data = dw->image->data; - uint16_t *old_flags = dw->image->flags; float **old_dp = dw->image->dp; int **old_bad = dw->image->bad; @@ -934,8 +930,6 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event) dw->ev_list->events[new_event], 0); if ( fail ) { ERROR("Couldn't load image"); - dw->image->data = old_data; - dw->image->flags = old_flags; dw->image->dp = old_dp; dw->image->bad = old_bad; return 1; @@ -956,8 +950,6 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event) free(old_dp[i]); free(old_bad[i]); } - free(old_data); - free(old_flags); free(old_dp); free(old_bad); return 0; diff --git a/src/geoptimiser.c b/src/geoptimiser.c index 0ac59ddf..4144bba2 100644 --- a/src/geoptimiser.c +++ b/src/geoptimiser.c @@ -2193,51 +2193,6 @@ struct rectangle }; -static int unpack_slab(struct image *image) -{ - struct detector *det = image->det; - int pi; - - image->dp = malloc(det->n_panels * sizeof(float *)); - image->bad = malloc(det->n_panels * sizeof(int *)); - if ( (image->dp == NULL) || (image->bad == NULL) ) { - ERROR("Failed to allocate panels.\n"); - return 1; - } - - for ( pi=0; pi<det->n_panels; pi++ ) { - - struct panel *p; - int fs, ss; - - p = &det->panels[pi]; - image->dp[pi] = malloc(p->w*p->h*sizeof(float)); - image->bad[pi] = calloc(p->w*p->h, sizeof(int)); - if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) ) { - ERROR("Failed to allocate panel\n"); - return 1; - } - - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - int idx; - int cfs, css; - - cfs = fs+p->min_fs; - css = ss+p->min_ss; - idx = cfs + css*image->width; - - image->dp[pi][fs+p->w*ss] = image->data[idx]; - image->bad[pi][fs+p->w*ss] = 0; - } - } - } - - return 0; -} - - static int draw_detector(cairo_surface_t *surf, struct image *image, struct rectangle rect) { @@ -2249,7 +2204,6 @@ static int draw_detector(cairo_surface_t *surf, struct image *image, cr = cairo_create(surf); - unpack_slab(image); pixbufs = render_panels(image, 1, SCALE_GEOPTIMISER, 1, &n_pixbufs); /* Blank grey background */ @@ -2304,25 +2258,51 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet, cairo_status_t r; cairo_surface_t *surf; - im.data = malloc((edet->width)*(edet->height)*sizeof(float)); - if ( im.data == NULL ) { - ERROR("Failed to allocate memory to save data.\n"); - return 1; - } im.det = edet->det; im.width = edet->width; im.height = edet->height; - im.flags = NULL; + im.dp = malloc(edet->det->n_panels*sizeof(float *)); + if ( im.dp == NULL ) { + ERROR("Failed to allocate data\n"); + return 1; + } + for ( i=0; i<edet->det->n_panels; i++ ) { - for ( i=0; i<(edet->width)*(edet->height); i++) { - if ( data[i] == -10000.0) { - im.data[i] = 0.0; - } else if ( data[i] > 1.0) { - im.data[i] = 1.0; - } else { - im.data[i] = (float)data[i]; + int fs, ss; + struct panel *p; + + p = &edet->det->panels[i]; + + im.dp[i] = calloc(p->w * p->h, sizeof(float)); + if ( im.dp[i] == NULL ) { + ERROR("Failed to allocate data\n"); + return 1; + } + + for ( ss=0; ss<p->h; ss++ ) { + for ( fs=0; fs<p->w; fs++ ) { + + int idx; + int cfs, css; + float val; + + cfs = fs+p->min_fs; + css = ss+p->min_ss; + idx = cfs + css*edet->width; + + if ( data[idx] == -10000.0) { + val = 0.0; + } else if ( data[idx] > 1.0) { + val = 1.0; + } else { + val = (float)data[idx]; + } + val *= 10.0; /* render_panels sets this as max */ + + im.dp[i][fs+p->w*ss] = val; + + } } - im.data[i] *= 10.0; /* render_panels sets this as max */ } get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x, @@ -2354,13 +2334,13 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet, cairo_fill(cr); cairo_destroy(cr); - r = cairo_surface_write_to_png(surf, filename); - if (r != CAIRO_STATUS_SUCCESS) { - free(im.data); - return 1; + for ( i=0; i<edet->det->n_panels; i++ ) { + free(im.dp[i]); } + free(im.dp); - free(im.data); + r = cairo_surface_write_to_png(surf, filename); + if ( r != CAIRO_STATUS_SUCCESS ) return 1; return 0; } diff --git a/src/hdfsee-render.c b/src/hdfsee-render.c index 71e2e1e8..4c6c3114 100644 --- a/src/hdfsee-render.c +++ b/src/hdfsee-render.c @@ -283,8 +283,9 @@ int render_tiff_fp(struct image *image, const char *filename) line = _TIFFmalloc(TIFFScanlineSize(th)); for ( y=0; y<image->height; y++ ) { - memcpy(line, &image->data[(image->height-1-y)*image->width], - image->width*4); + //memcpy(line, &image->data[(image->height-1-y)*image->width], + // image->width*4); + // FIXME! TIFFWriteScanline(th, line, y, 0); } _TIFFfree(line); @@ -325,7 +326,8 @@ int render_tiff_int16(struct image *image, const char *filename, double boost) for ( y=0; y<image->height; y++ ) { for ( x=0;x<image->width; x++ ) { double val; - val = image->data[x+image->height*y]; + //val = image->data[x+image->height*y]; + val = 0.0; // FIXME! if ( val > max ) max = val; } } @@ -336,7 +338,8 @@ int render_tiff_int16(struct image *image, const char *filename, double boost) double val; - val = image->data[x+(image->height-1-y)*image->width]; + //val = image->data[x+(image->height-1-y)*image->width]; + val = 0.0; // FIXME! val *= ((double)boost/max); /* Clamp to 16-bit range, diff --git a/src/partial_sim.c b/src/partial_sim.c index 867183da..16e9cb88 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -170,6 +170,18 @@ static void calculate_partials(Crystal *cr, } +static signed int panel_number(struct detector *det, struct panel *p) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + if ( &det->panels[i] == p ) return i; + } + + return -1; +} + + static void draw_and_write_image(struct image *image, RefList *reflections, gsl_rng *rng, double background) { @@ -177,7 +189,26 @@ static void draw_and_write_image(struct image *image, RefList *reflections, RefListIterator *iter; int i; - image->data = calloc(image->width*image->height, sizeof(float)); + image->dp = malloc(image->det->n_panels*sizeof(float *)); + if ( image->dp == NULL ) { + ERROR("Failed to allocate data\n"); + return; + } + for ( i=0; i<image->det->n_panels; i++ ) { + + int j; + struct panel *p = &image->det->panels[i]; + + image->dp[i] = calloc(p->w * p->h, sizeof(float)); + if ( image->dp[i] == NULL ) { + ERROR("Failed to allocate data\n"); + return; + } + for ( j=0; j<p->w*p->h; j++ ) { + image->dp[i][j] = poisson_noise(rng, background); + } + + } for ( refl = first_refl(reflections, &iter); refl != NULL; @@ -186,26 +217,32 @@ static void draw_and_write_image(struct image *image, RefList *reflections, double Ip; double dfs, dss; int fs, ss; + struct panel *p; + signed int pn; Ip = get_intensity(refl); get_detector_pos(refl, &dfs, &dss); - fs = nearbyint(dfs); - ss = nearbyint(dss); + p = get_panel(refl); + pn = panel_number(image->det, p); /* Yeurgh */ + assert(pn != -1); + + fs = nearbyint(dfs) - p->min_fs; + ss = nearbyint(dss) - p->min_ss; assert(fs >= 0); assert(ss >= 0); - assert(fs < image->width); - assert(ss < image->height); - image->data[fs + image->width*ss] += Ip; - - } + assert(fs < p->w); + assert(ss < p->h); + image->dp[pn][fs + p->w*ss] += Ip; - for ( i=0; i<image->width*image->height; i++ ) { - image->data[i] += poisson_noise(rng, background); } hdf5_write_image(image->filename, image, NULL); - free(image->data); + + for ( i=0; i<image->det->n_panels; i++ ) { + free(image->dp[i]); + } + free(image->dp); } diff --git a/src/pattern_sim.c b/src/pattern_sim.c index bd97ead7..f724057e 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -249,8 +249,7 @@ int main(int argc, char *argv[]) int c; struct image image; struct gpu_context *gctx = NULL; - struct image *powder; - float *powder_data; + struct image powder; char *intfile = NULL; double *intensities; char *rval; @@ -293,6 +292,7 @@ int main(int argc, char *argv[]) double bandwidth = 0.01; double photon_energy = 9000.0; struct beam_params beam; + int i; /* Long options */ const struct option longopts[] = { @@ -667,7 +667,21 @@ int main(int argc, char *argv[]) /* Initialise stuff */ image.filename = NULL; image.features = NULL; - image.flags = NULL; + image.bad = NULL; + image.dp = malloc(image.det->n_panels*sizeof(float *)); + if ( image.dp == NULL ) { + ERROR("Failed to allocate data\n"); + return 1; + } + for ( i=0; i<image.det->n_panels; i++ ) { + image.dp[i] = calloc(image.det->panels[i].w * + image.det->panels[i].h, + sizeof(float)); + if ( image.dp[i] == NULL ) { + ERROR("Failed to allocate data\n"); + return 1; + } + } rng = gsl_rng_alloc(gsl_rng_mt19937); if ( config_random ) { @@ -679,12 +693,23 @@ int main(int argc, char *argv[]) gsl_rng_set(rng, seed); } - powder = calloc(1, sizeof(struct image)); - powder->width = image.width; - powder->height = image.height; - powder->det = image.det; - powder_data = calloc(image.width*image.height, sizeof(float)); - powder->data = powder_data; + powder.width = image.width; + powder.height = image.height; + powder.det = image.det; + powder.dp = malloc(image.det->n_panels*sizeof(float *)); + if ( powder.dp == NULL ) { + ERROR("Failed to allocate powder data\n"); + return 1; + } + for ( i=0; i<image.det->n_panels; i++ ) { + powder.dp[i] = calloc(image.det->panels[i].w * + image.det->panels[i].h, + sizeof(float)); + if ( powder.dp[i] == NULL ) { + ERROR("Failed to allocate powder data\n"); + return 1; + } + } /* Splurge a few useful numbers */ STATUS("Simulation parameters:\n"); @@ -732,6 +757,7 @@ int main(int argc, char *argv[]) int na, nb, nc; double a, b, c, d; UnitCell *cell; + int err = 0; if ( random_size ) { @@ -825,10 +851,6 @@ int main(int argc, char *argv[]) } - /* Ensure no residual information */ - image.data = NULL; - image.twotheta = NULL; - cell_get_parameters(cell, &a, &b, &c, &d, &d, &d); STATUS("Particle size = %i x %i x %i" " ( = %5.2f x %5.2f x %5.2f nm)\n", @@ -836,8 +858,6 @@ int main(int argc, char *argv[]) if ( config_gpu ) { - int err; - if ( gctx == NULL ) { gctx = setup_gpu(config_nosfac, intensities, flags, sym_str, @@ -845,13 +865,12 @@ int main(int argc, char *argv[]) } err = get_diffraction_gpu(gctx, &image, na, nb, nc, cell, no_fringes); - if ( err ) image.data = NULL; } else { get_diffraction(&image, na, nb, nc, intensities, phases, flags, cell, grad, sym, no_fringes); } - if ( image.data == NULL ) { + if ( err ) { ERROR("Diffraction calculation failed.\n"); done = 1; goto skip; @@ -862,18 +881,23 @@ int main(int argc, char *argv[]) if ( powder_fn != NULL ) { - int x, y, w; + int pn; - w = image.width; + for ( pn=0; pn<image.det->n_panels; pn++ ) { + + size_t w, i; + + w = image.det->panels[pn].w + * image.det->panels[pn].h; + + for ( i=0; i<w; i++ ) { + powder.dp[pn][i] += (double)image.dp[pn][i]; + } - for ( x=0; x<image.width; x++ ) { - for ( y=0; y<image.height; y++ ) { - powder->data[x+w*y] += (double)image.data[x+w*y]; - } } if ( !(ndone % 10) ) { - hdf5_write_image(powder_fn, powder, NULL); + hdf5_write_image(powder_fn, &powder, NULL); } } @@ -895,10 +919,6 @@ int main(int argc, char *argv[]) } - /* Clean up */ - free(image.data); - free(image.twotheta); - cell_free(cell); skip: @@ -909,18 +929,18 @@ skip: } while ( !done ); if ( powder_fn != NULL ) { - hdf5_write_image(powder_fn, powder, NULL); + hdf5_write_image(powder_fn, &powder, NULL); } if ( gctx != NULL ) { cleanup_gpu(gctx); } - free(image.det->panels); + for ( i=0; i<image.det->n_panels; i++ ) { + free(powder.dp[i]); + } free(intfile); free(image.det); - free(powder->data); - free(powder); cell_free(input_cell); free(intensities); free(outfile); diff --git a/src/process_image.c b/src/process_image.c index e7e3aa78..4fd19176 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -135,8 +135,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, int any_crystals; image.features = NULL; - image.data = NULL; - image.flags = NULL; image.copyme = iargs->copyme; image.id = cookie; image.filename = pargs->filename_p_e->filename; @@ -343,8 +341,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, free(image.dp); free(image.bad); - free(image.data); - if ( image.flags != NULL ) free(image.flags); image_feature_list_free(image.features); free_detector_geometry(image.det); hdfile_close(hdfile); diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c index dbe59ce2..7c7c8e0c 100644 --- a/tests/gpu_sim_check.c +++ b/tests/gpu_sim_check.c @@ -3,11 +3,11 @@ * * Check that GPU simulation agrees with CPU version * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2012-2014 Thomas White <taw@physics.org> + * 2012-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -110,11 +110,19 @@ int main(int argc, char *argv[]) det = calloc(1, sizeof(struct detector)); det->n_panels = 2; det->panels = calloc(2, sizeof(struct panel)); + det->max_fs = 1024; + det->max_ss = 1024; det->panels[0].min_fs = 0; det->panels[0].max_fs = 1023; det->panels[0].min_ss = 0; det->panels[0].max_ss = 511; + det->panels[0].orig_min_fs = 0; + det->panels[0].orig_max_fs = 1023; + det->panels[0].orig_min_ss = 0; + det->panels[0].orig_max_ss = 511; + det->panels[0].w = 1024; + det->panels[0].h = 512; det->panels[0].fsx = 1; det->panels[0].fsy = 0; det->panels[0].ssx = 0; @@ -128,11 +136,18 @@ int main(int argc, char *argv[]) det->panels[0].clen = 100.0e-3; det->panels[0].res = 9090.91; det->panels[0].adu_per_eV = 1.0; + det->panels[0].data = NULL; det->panels[1].min_fs = 0; det->panels[1].max_fs = 1023; det->panels[1].min_ss = 512; det->panels[1].max_ss = 1023; + det->panels[1].orig_min_fs = 0; + det->panels[1].orig_max_fs = 1023; + det->panels[1].orig_min_ss = 512; + det->panels[1].orig_max_ss = 1023; + det->panels[1].w = 1024; + det->panels[1].h = 512; det->panels[1].fsx = 1; det->panels[1].fsy = 0; det->panels[1].ssx = 0; @@ -145,12 +160,15 @@ int main(int argc, char *argv[]) det->panels[1].cny = sep; det->panels[1].clen = 100.0e-3; det->panels[1].res = 9090.91; - det->panels[0].adu_per_eV = 1.0; + det->panels[1].adu_per_eV = 1.0; + det->panels[1].data = NULL; cpu_image.det = det; gpu_image.det = det; cpu_image.beam = NULL; gpu_image.beam = NULL; + cpu_image.spectrum = NULL; + gpu_image.spectrum = NULL; cpu_image.lambda = ph_en_to_lambda(eV_to_J(6000)); gpu_image.lambda = ph_en_to_lambda(eV_to_J(6000)); @@ -171,6 +189,20 @@ int main(int argc, char *argv[]) sym = get_pointgroup("1"); + cpu_image.dp = malloc(det->n_panels * sizeof(float *)); + if ( cpu_image.dp == NULL ) { + ERROR("Couldn't allocate memory for result.\n"); + return 1; + } + for ( i=0; i<det->n_panels; i++ ) { + struct panel *p = &det->panels[i]; + cpu_image.dp[i] = calloc(p->w * p->h, sizeof(float)); + if ( cpu_image.dp[i] == NULL ) { + ERROR("Couldn't allocate memory for panel %i\n", i); + return 1; + } + } + start = get_hires_seconds(); get_diffraction(&cpu_image, 8, 8, 8, NULL, NULL, NULL, cell, GRADIENT_MOSAIC, sym, 1); @@ -184,18 +216,25 @@ int main(int argc, char *argv[]) gpu_min = +INFINITY; gpu_max = -INFINITY; gpu_tot = 0.0; cpu_min = +INFINITY; cpu_max = -INFINITY; cpu_tot = 0.0; dev = 0.0; - for ( i=0; i<1024*1024; i++ ) { + for ( i=0; i<det->n_panels; i++ ) { + + int j; + struct panel *p = &det->panels[i]; + + for ( j=0; j<p->w*p->h; j++ ) { + + const double cpu = cpu_image.dp[i][j]; + const double gpu = gpu_image.dp[i][j]; - const double cpu = cpu_image.data[i]; - const double gpu = gpu_image.data[i]; + if ( cpu > cpu_max ) cpu_max = cpu; + if ( cpu < cpu_min ) cpu_min = cpu; + if ( gpu > gpu_max ) gpu_max = gpu; + if ( gpu < gpu_min ) gpu_min = gpu; + gpu_tot += gpu; + cpu_tot += cpu; + dev += fabs(gpu - cpu); - if ( cpu > cpu_max ) cpu_max = cpu; - if ( cpu < cpu_min ) cpu_min = cpu; - if ( gpu > gpu_max ) gpu_max = gpu; - if ( gpu < gpu_min ) gpu_min = gpu; - gpu_tot += gpu; - cpu_tot += cpu; - dev += fabs(gpu - cpu); + } } perc = 100.0*dev/cpu_tot; @@ -204,25 +243,22 @@ int main(int argc, char *argv[]) STATUS("CPU: min=%8e, max=%8e, total=%8e\n", cpu_min, cpu_max, cpu_tot); STATUS("dev = %8e (%5.2f%% of CPU total)\n", dev, perc); - cell_free(cell); - free_detector_geometry(det); - - if ( perc > 1.0 ) { + if ( 1 ) { STATUS("Test failed! I'm writing cpu-sim.h5 and gpu-sim.h5" " for you to inspect.\n"); - hdf5_write("cpu-sim.h5", cpu_image.data, cpu_image.width, - cpu_image.height, H5T_NATIVE_FLOAT); - - hdf5_write("gpu-sim.h5", gpu_image.data, gpu_image.width, - gpu_image.height, H5T_NATIVE_FLOAT); + hdf5_write_image("cpu-sim.h5", &cpu_image, NULL); + hdf5_write_image("gpu-sim.h5", &gpu_image, NULL); return 1; } gsl_rng_free(rng); + cell_free(cell); + free_detector_geometry(det); + return 0; } diff --git a/tests/integration_check.c b/tests/integration_check.c index ccb613c3..2acd0f34 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -65,7 +65,6 @@ int main(int argc, char *argv[]) fclose(fh); gsl_rng_set(rng, seed); - image.flags = NULL; image.beam = NULL; image.lambda = ph_eV_to_lambda(9000.0); diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c index c7ebd8b2..f66b6513 100644 --- a/tests/prof2d_check.c +++ b/tests/prof2d_check.c @@ -77,7 +77,6 @@ int main(int argc, char *argv[]) fclose(fh); gsl_rng_set(rng, seed); - image.flags = NULL; image.beam = NULL; image.lambda = ph_eV_to_lambda(9000.0); image.bw = 0.000001; diff --git a/tests/ring_check.c b/tests/ring_check.c index 8ced12d5..99c433c2 100644 --- a/tests/ring_check.c +++ b/tests/ring_check.c @@ -62,7 +62,7 @@ static void third_integration_check(struct image *image, int n_trials, for ( fs=0; fs<image->width; fs++ ) { for ( ss=0; ss<image->height; ss++ ) { - image->data[fs+image->width*ss] + image->dp[0][fs+image->width*ss] = poisson_noise(rng, 1000.0); } } @@ -121,9 +121,9 @@ static void fourth_integration_check(struct image *image, int n_trials, for ( fs=0; fs<image->width; fs++ ) { for ( ss=0; ss<image->height; ss++ ) { int idx = fs+image->width*ss; - image->data[idx] = poisson_noise(rng, 1000.0); + image->dp[0][idx] = poisson_noise(rng, 1000.0); if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue; - image->data[idx] += 1000.0; + image->dp[0][idx] += 1000.0; pcount++; } } @@ -170,7 +170,6 @@ int main(int argc, char *argv[]) int r, npx; double ex; gsl_rng *rng; - int *bad; rng = gsl_rng_alloc(gsl_rng_mt19937); @@ -179,10 +178,10 @@ int main(int argc, char *argv[]) fclose(fh); gsl_rng_set(rng, seed); - image.data = malloc(128*128*sizeof(float)); - bad = calloc(128*128, sizeof(int)); - image.bad = &bad; - image.flags = NULL; + image.dp = malloc(sizeof(float *)); + image.dp[0] = malloc(128*128*sizeof(float)); + image.bad = malloc(sizeof(uint16_t *)); + image.bad[0] = calloc(128*128, sizeof(uint16_t)); image.beam = NULL; image.lambda = ph_eV_to_lambda(1000.0); @@ -191,9 +190,9 @@ int main(int argc, char *argv[]) image.det->panels = calloc(1, sizeof(struct panel)); image.det->panels[0].min_fs = 0; - image.det->panels[0].max_fs = 128; + image.det->panels[0].max_fs = 127; image.det->panels[0].min_ss = 0; - image.det->panels[0].max_ss = 128; + image.det->panels[0].max_ss = 127; image.det->panels[0].fsx = 1.0; image.det->panels[0].fsy = 0.0; image.det->panels[0].ssx = 0.0; @@ -206,12 +205,14 @@ int main(int argc, char *argv[]) image.det->panels[0].cny = -64.0; image.det->panels[0].clen = 1.0; image.det->panels[0].res = 1.0; + image.det->panels[0].w = 128; + image.det->panels[0].h = 128; image.det->panels[0].adu_per_eV = 1.0/1000.0; /* -> 1 adu per photon */ image.det->panels[0].max_adu = +INFINITY; /* No cutoff */ image.width = 128; image.height = 128; - memset(image.data, 0, 128*128*sizeof(float)); + memset(image.dp[0], 0, 128*128*sizeof(float)); image.n_crystals = 0; image.crystals = NULL; @@ -238,7 +239,7 @@ int main(int argc, char *argv[]) for ( fs=0; fs<image.width; fs++ ) { for ( ss=0; ss<image.height; ss++ ) { if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue; - image.data[fs+image.width*ss] = 1000.0; + image.dp[0][fs+image.width*ss] = 1000.0; npx++; } } @@ -278,9 +279,9 @@ int main(int argc, char *argv[]) npx = 0; for ( fs=0; fs<image.width; fs++ ) { for ( ss=0; ss<image.height; ss++ ) { - image.data[fs+image.width*ss] = 1000.0; + image.dp[0][fs+image.width*ss] = 1000.0; if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue; - image.data[fs+image.width*ss] += 1000.0; + image.dp[0][fs+image.width*ss] += 1000.0; npx++; } } @@ -310,11 +311,11 @@ int main(int argc, char *argv[]) } - free(image.beam); free(image.det->panels); free(image.det); - free(image.data); + free(image.dp[0]); + free(image.dp); gsl_rng_free(rng); if ( fail ) return 1; |