Skip to content

Commit

Permalink
27.09.2022 - 25.10.2022: add forgotten changes from Skin MBP: histogr…
Browse files Browse the repository at this point in the history 
…am on multiple blocks of data
  • Loading branch information
@diemoschwarz
diemoschwarz committed Jul 17, 2024
1 parent e506823 commit ef79c9d
Show file tree
Hide file tree
Showing 2 changed files with 150 additions and 66 deletions.
171 changes: 107 additions & 64 deletions src/statistics/rta_histogram.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,139 +36,182 @@
* POSSIBILITY OF SUCH DAMAGE.
*/

#include <float.h>
#include "rta_histogram.h"


// init parameter struct to default values
void rta_histogram_init (rta_histogram_params_t *self)
void rta_histogram_init (rta_histogram_params_t *params)
{
self->nhist = 100;
self->lo_given = 0;
self->hi_given = 0;
self->lo = 0;
self->hi = 0;
self->norm = 0;
params->nhist = 100;
params->lo_given = 0;
params->hi_given = 0;
params->lo = 0;
params->hi = 0;
params->norm = 0;
}

/* Calculate histogram */
void rta_histogram_stride (rta_histogram_params_t *self,
void rta_histogram_stride (rta_histogram_params_t *params,
rta_real_t *input, const unsigned int i_stride, const unsigned int i_size,
rta_real_t *output, const unsigned int out_stride,
rta_real_t *bpfout, const unsigned int bpf_stride)
{
rta_real_t one = 1.0;
rta_real_t one = 1;

rta_histogram_weighted_stride(self, input, i_stride, i_size, &one, 0, // unweighted: all weights == 1
output, out_stride, bpfout, bpf_stride);
rta_histogram_weighted_stride(params,
input, i_stride, i_size,
&one, 0, // unweighted: all weights == 1
output, out_stride,
bpfout, bpf_stride);
}

void rta_histogram_stride_multi (rta_histogram_params_t *params, int num_input,
rta_real_t *input[], const unsigned int i_stride, const unsigned int i_size[],
rta_real_t *output, const unsigned int out_stride,
rta_real_t *bpfout, const unsigned int bpf_stride)
{
rta_real_t ones[num_input];

for (int i = 0; i < num_input; i++)
ones[i] = 1.0;

rta_histogram_weighted_stride_multi(params, num_input,
input, i_stride, i_size,
ones, 0, // unweighted: all weights == 1
output, out_stride,
bpfout, bpf_stride);
}

/* Calculate weighted histogram */
void rta_histogram_weighted_stride (rta_histogram_params_t *self,
void rta_histogram_weighted_stride (rta_histogram_params_t *params,
rta_real_t *input, const unsigned int i_stride, const unsigned int i_size,
rta_real_t *weights, const unsigned int w_stride,
rta_real_t *output, const unsigned int out_stride,
rta_real_t *bpfout, const unsigned int bpf_stride)
{
rta_histogram_weighted_stride_multi(params, 1,
&input, i_stride, &i_size,
&weights, w_stride, // unweighted: all weights == 1
output, out_stride,
bpfout, bpf_stride);
}

void rta_histogram_weighted_stride_multi (rta_histogram_params_t *params, int num_input,
rta_real_t *input[], const unsigned int i_stride, const unsigned int i_size[],
rta_real_t *weights[], const unsigned int w_stride,
rta_real_t *output, const unsigned int out_stride,
rta_real_t *bpfout, const unsigned int bpf_stride)
{
rta_real_t *ptr;
int i;
int i, j, k;

/* clear result matrix */
ptr = output;
for (i = 0; i < self->nhist; i++, ptr += out_stride)
for (i = 0; i < params->nhist; i++, ptr += out_stride)
*ptr = 0;

if ((ptr = bpfout)) // clear bpf domain output if given
for (i = 0; i < self->nhist; i++, ptr += bpf_stride)
for (i = 0; i < params->nhist; i++, ptr += bpf_stride)
*ptr = 0;

if (i_size <= 0)
// calculate and check input data
unsigned int numdata = 0;
for (k = 0; k < num_input; k++)
numdata += i_size[k];

if (numdata == 0)
return; // no data

/* find min/max if not given by attributes */
if (!self->lo_given && !self->hi_given)
if (!params->lo_given && !params->hi_given)
{
self->lo = *input; /* lower histogram limit */
self->hi = *input; /* upper histogram limit */

for (i = i_stride; i < i_size * i_stride; i += i_stride)
{
float x = input[i];
if (x < self->lo)
self->lo = x;
if (x > self->hi)
self->hi = x;
}
params->lo = FLT_MAX; /* lower histogram limit */
params->hi = -FLT_MAX; /* upper histogram limit */

for (k = 0; k < num_input; k++)
for (i = 0; i < i_size[k] * i_stride; i += i_stride)
{
float x = input[k][i];
if (x < params->lo)
params->lo = x;
if (x > params->hi)
params->hi = x;
}
}
else
{
if (!self->lo_given)
if (!params->lo_given)
{
self->lo = *input; /* lower histogram limit */
params->lo = FLT_MAX; /* lower histogram limit */

for (i = i_stride; i < i_size * i_stride; i += i_stride)
{
float x = input[i];
if (x < self->lo)
self->lo = x;
}
for (k = 0; k < num_input; k++)
for (i = 0; i < i_size[k] * i_stride; i += i_stride)
{
float x = input[k][i];
if (x < params->lo)
params->lo = x;
}
}

if (!self->hi_given)
if (!params->hi_given)
{
self->hi = *input; /* upper histogram limit */
params->hi = -FLT_MAX; /* upper histogram limit */

for (i = i_stride; i < i_size * i_stride; i += i_stride)
{
float x = input[i];
if (x > self->hi)
self->hi = x;
}
for (k = 0; k < num_input; k++)
for (i = 0; i < i_size[k] * i_stride; i += i_stride)
{
float x = input[k][i];
if (x > params->hi)
params->hi = x;
}
}
}

float xfact = self->nhist / (self->hi - self->lo + 1); // bin step
float xfact = params->nhist / (params->hi - params->lo + 1); // bin step

if (bpfout)
{ /* create bin values */
float bfact = 1 / xfact;

for (i = 0; i < self->nhist; i++)
for (i = 0; i < params->nhist; i++)
{
*bpfout = i * bfact + self->lo;
*bpfout = i * bfact + params->lo;
bpfout += bpf_stride;
}
}

/* calculate histogram */
for (i = 0; i < i_size * i_stride; i += i_stride, weights += w_stride)
{
/* find bin index */
int ind = (input[i] - self->lo) * xfact;
for (k = 0; k < num_input; k++)
for (i = 0, j = 0; i < i_size[k] * i_stride; i += i_stride, j += w_stride)
{
/* find bin index */
int ind = (input[k][i] - params->lo) * xfact;

/* clip in case bounds were given as attributes */
if (ind < 0)
ind = 0;
if (ind >= self->nhist)
ind = self->nhist - 1;
/* clip in case bounds were given as attributes */
if (ind < 0)
ind = 0;
if (ind >= params->nhist)
ind = params->nhist - 1;

output[ind * out_stride] += *weights;
}
output[ind * out_stride] += weights[k][j];
}

// normalise histogram
if (self->norm > 0)
if (params->norm > 0)
{
float normfact = 0;
rta_real_t *histcol = output;

if (self->norm == 1)
if (params->norm == 1)
{ // find max
for (i = 0; i < self->nhist; i++, histcol += out_stride)
for (i = 0; i < params->nhist; i++, histcol += out_stride)
if (*histcol > normfact)
normfact = *histcol;
}
else if (self->norm == 2)
else if (params->norm == 2)
{ // calc sum
for (i = 0; i < self->nhist; i++, histcol += out_stride)
for (i = 0; i < params->nhist; i++, histcol += out_stride)
normfact += *histcol;
}

Expand All @@ -178,7 +221,7 @@ void rta_histogram_weighted_stride (rta_histogram_params_t *self,
normfact = 1 / normfact;
histcol = output;

for (i = 0; i < self->nhist; i++, histcol += out_stride)
for (i = 0; i < params->nhist; i++, histcol += out_stride)
*histcol *= normfact;
}
}
Expand Down
45 changes: 43 additions & 2 deletions src/statistics/rta_histogram.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,28 @@ void rta_histogram_init (rta_histogram_params_t *self);
* @param bin_stride stride for bin index data
*/
void rta_histogram_stride (rta_histogram_params_t *params,
rta_real_t *input, const unsigned int i_stride, const unsigned int i_size,
rta_real_t *input, const unsigned int i_stride, const unsigned int i_size,
rta_real_t *output, const unsigned int out_stride,
rta_real_t *binout, const unsigned int bin_stride);


/**
* Calculate histogram over multiple blocks of data
*
* @param params pointer to histogram parameter struct
* @param num_input number of blocks of input data
* @param input array[num_input] of pointers to input data (at least i_size * i_stride elements)
* @param i_stride stride for input data
* @param i_size array[num_input] of number of elements per block of input data
* @param output pointer to output data (at least params->nhist * out_stride elements)
* @param out_stride stride for output data
* @param binout NULL or pointer to bin index output data (at least params->nhist * bpf_stride elements)
* @param bin_stride stride for bin index data
*/
void rta_histogram_stride_multi (rta_histogram_params_t *params, int num_input,
rta_real_t *input[], const unsigned int i_stride, const unsigned int i_size[],
rta_real_t *output, const unsigned int out_stride,
rta_real_t *binout, const unsigned int bin_stride);

/**
* Calculate weighted histogram
* (bins don't count occurrences, but sum weights given with each data value)
Expand All @@ -99,6 +117,29 @@ void rta_histogram_weighted_stride (rta_histogram_params_t *params,
rta_real_t *weights, const unsigned int w_stride,
rta_real_t *output, const unsigned int out_stride,
rta_real_t *binout, const unsigned int bin_stride);


/**
* Calculate weighted histogram on multiple input blocks
* (bins don't count occurrences, but sum weights given with each data value)
*
* @param params pointer to histogram parameter struct
* @param num_input number of blocks of input data
* @param input array[num_input] of pointers to input data (at least i_size * i_stride elements)
* @param i_stride stride for input data
* @param i_size array[num_input] of number of input elements
* @param weights array[num_input] of pointers to weights data (at least i_size * w_stride elements)
* @param w_stride stride for weights data
* @param output pointer to output data (at least params->nhist * out_stride elements)
* @param out_stride stride for output data
* @param binout NULL or pointer to bin index output data (at least params->nhist * bpf_stride elements)
* @param bin_stride stride for bin index data
*/
void rta_histogram_weighted_stride_multi (rta_histogram_params_t *params, int num_input,
rta_real_t *input[], const unsigned int i_stride, const unsigned int i_size[],
rta_real_t *weights[], const unsigned int w_stride,
rta_real_t *output, const unsigned int out_stride,
rta_real_t *binout, const unsigned int bin_stride);

#ifdef __cplusplus
}
Expand Down

0 comments on commit ef79c9d

Please sign in to comment.