diff options
Diffstat (limited to 'src')
-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 |
8 files changed, 204 insertions, 172 deletions
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); |