diff options
author | Thomas White <taw@physics.org> | 2015-11-03 10:54:19 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2015-11-03 10:54:19 +0100 |
commit | 2cfd03ca948b53dc5e57622f8c3eb50156c03231 (patch) | |
tree | e094d48f39ed33859f1da45a404d70ccadf8355b /libcrystfel | |
parent | 832374fb4578ed5252bd3bcbdf699833115088f9 (diff) | |
parent | 478c76bc45ae1a7cbb7be1f10306a0ae985a41b6 (diff) |
Merge branch 'tom/imagedata'
Diffstat (limited to 'libcrystfel')
-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 |
5 files changed, 194 insertions, 149 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); } |