diff options
Diffstat (limited to 'src/silx/math/histogramnd/src/histogramnd_template.c')
| -rw-r--r-- | src/silx/math/histogramnd/src/histogramnd_template.c | 260 |
1 files changed, 260 insertions, 0 deletions
diff --git a/src/silx/math/histogramnd/src/histogramnd_template.c b/src/silx/math/histogramnd/src/histogramnd_template.c new file mode 100644 index 0000000..0276bb4 --- /dev/null +++ b/src/silx/math/histogramnd/src/histogramnd_template.c @@ -0,0 +1,260 @@ +/*########################################################################## +# Copyright (C) 2016 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ############################################################################*/ + +#include "templates.h" + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include <stdarg.h> + +#ifdef HISTO_SAMPLE_T +#ifdef HISTO_WEIGHT_T +#ifdef HISTO_CUMUL_T + +int TEMPLATE(histogramnd, HISTO_SAMPLE_T, HISTO_WEIGHT_T, HISTO_CUMUL_T) + (HISTO_SAMPLE_T *i_sample, + HISTO_WEIGHT_T *i_weights, + int i_n_dim, + int i_n_elem, + double *i_bin_ranges, + int *i_n_bins, + uint32_t *o_histo, + HISTO_CUMUL_T *o_cumul, + double *o_bin_edges, + int i_opt_flags, + HISTO_WEIGHT_T i_weight_min, + HISTO_WEIGHT_T i_weight_max) +{ + /* some counters */ + int i = 0, j = 0; + long elem_idx = 0; + + HISTO_WEIGHT_T * weight_ptr = 0; + HISTO_SAMPLE_T elem_coord = 0.; + + /* computed bin index (i_sample -> grid) */ + long bin_idx = 0; + + double * g_min = 0; + double * g_max = 0; + double * range = 0; + + /* ================================ + * Parsing options, if any. + * ================================ + */ + + int filt_min_weight = 0; + int filt_max_weight = 0; + int last_bin_closed = 0; + + /* Testing the option flags */ + if(i_opt_flags & HISTO_WEIGHT_MIN) + { + filt_min_weight = 1; + } + + if(i_opt_flags & HISTO_WEIGHT_MAX) + { + filt_max_weight = 1; + } + + if(i_opt_flags & HISTO_LAST_BIN_CLOSED) + { + last_bin_closed = 1; + } + + /* storing the min & max bin coordinates in their own arrays because + * i_bin_ranges = [[min0, max0], [min1, max1], ...] + * (mostly for the sake of clarity) + * (maybe faster access too?) + */ + g_min = (double *) malloc(i_n_dim *sizeof(double)); + g_max = (double *) malloc(i_n_dim * sizeof(double)); + /* range used to convert from i_coords to bin indices in the grid */ + range = (double *) malloc(i_n_dim * sizeof(double)); + + if(!g_min || !g_max || !range) + { + free(g_min); + free(g_max); + free(range); + return HISTO_ERR_ALLOC; + } + + j = 0; + for(i=0; i<i_n_dim; i++) + { + g_min[i] = i_bin_ranges[i*2]; + g_max[i] = i_bin_ranges[i*2+1]; + range[i] = g_max[i]-g_min[i]; + + for(bin_idx=0; bin_idx<i_n_bins[i]; j++, bin_idx++) + { + o_bin_edges[j] = g_min[i] + + bin_idx * (range[i] / i_n_bins[i]); + } + o_bin_edges[j++] = g_max[i]; + } + + weight_ptr = i_weights; + + if(!i_weights) + { + /* if weights are not provided there no point in trying to filter them + * (!! careful if you change this, some code below relies on it !!) + */ + filt_min_weight = 0; + filt_max_weight = 0; + + /* If the weights array is not provided then there is no point + * updating the weighted histogram, only the bin counts (o_histo) + * will be filled. + * (!! careful if you change this, some code below relies on it !!) + */ + o_cumul = 0; + } + + /* tried to use pointers instead of indices here, but it didn't + * seem any faster (probably because the compiler + * optimizes stuff anyway), + * so i'm keeping the "indices" version, for the sake of clarity + */ + for(elem_idx=0; + elem_idx<i_n_elem*i_n_dim; + elem_idx+=i_n_dim, weight_ptr++) + { + /* no testing the validity of weight_ptr here, because if it is NULL + * then filt_min_weight/filt_max_weight will be 0. + * (see code above) + */ + if(filt_min_weight && *weight_ptr<i_weight_min) + { + continue; + } + if(filt_max_weight && *weight_ptr>i_weight_max) + { + continue; + } + + bin_idx = 0; + + for(i=0; i<i_n_dim; i++) + { + elem_coord = i_sample[elem_idx+i]; + + /* ===================== + * Element is rejected if any of the following is NOT true : + * 1. coordinate is >= than the minimum value + * 2. coordinate is <= than the maximum value + * 3. coordinate==maximum value and last_bin_closed is True + * ===================== + */ + if(elem_coord<g_min[i]) + { + bin_idx = -1; + break; + } + + /* Here we make the assumption that most of the time + * there will be more coordinates inside the grid interval + * (one test) + * than coordinates higher or equal to the max + * (two tests) + */ + if(elem_coord<g_max[i]) + { + /* Warning : the following factorization seems to + * increase the effect of precision error. + * bin_idx = (long)floor( + * (bin_idx + + * (elem_coord-g_min[i])/range[i]) * + * i_n_bins[i] + * ); + */ + + /* Not using floor to speed up things. + * We don't (?) need all the error checking provided by + * the built-in floor(). + * Also the value is supposed to be always positive. + */ + bin_idx = bin_idx * i_n_bins[i] + + (long)( + ((elem_coord-g_min[i]) * i_n_bins[i]) / + range[i] + ); + } + else /* ===> elem_coord>=g_max[i] */ + { + /* if equal and the last bin is closed : + * put it in the last bin + * else : discard + */ + if(last_bin_closed && elem_coord==g_max[i]) + { + bin_idx = (bin_idx + 1) * i_n_bins[i] - 1; + } + else + { + bin_idx = -1; + break; + } + } /* if(elem_coord<g_max[i]) */ + + } /* for(i=0; i<i_n_dim; i++) */ + + /* element is out of the grid */ + if(bin_idx==-1) + { + continue; + } + + if(o_histo) + { + o_histo[bin_idx] += 1; + } + if(o_cumul) + { + /* not testing the pointer since o_cumul is null if + * i_weights is null. + */ + o_cumul[bin_idx] += (HISTO_CUMUL_T) *weight_ptr; + } + + } /* for(elem_idx=0; elem_idx<i_n_elem*i_n_dim; elem_idx+=i_n_dim) */ + + free(g_min); + free(g_max); + free(range); + + /* For now just returning 0 (OK) since all the checks are done in + * python. This might change later if people want to call this + * function directly from C (might have to implement error codes). + */ + return HISTO_OK; +} + +#endif +#endif +#endif |