cross_corr_helpers.c

#include <math.h>

typedef struct {
    float mean;
    float std;
} MeanStd_f;

float compute_mean_f(
    const float *a,
    const int *x_dim,
    const int *y_dim
)
{
    float acc = 0.0;
    for(size_t i = 0; i < *x_dim; i++)
    {
        for(size_t j = 0; j < *y_dim; j++)
        {
            acc += a[i * *x_dim + j];
        }
    }
    return acc / (*x_dim * *y_dim);
}

void compute_mean_std_f(
    const float *a,
    const int *x_dim,
    const int *y_dim,
    MeanStd_f *mean_std
)
{
    // Get the mean
    mean_std->mean = compute_mean_f(a, x_dim, y_dim);
    
    // Compute the absolute square differences
    float val = 0.0;
    float acc = 0.0;
    for (int i = 0; i < *x_dim; i++)
    {
        for(int j = 0; j < *y_dim; j++)
        {
            val = fabsf(a[i * *x_dim + j] - mean_std->mean);
            acc += val * val;
        }
    }

    mean_std->std = sqrtf(acc / (*x_dim * *y_dim));
}

void zscore_normalize_f(
    const float *a,
    const int *x_dim,
    const int *y_dim,
    float *result
)
{
    MeanStd_f mean_std;
    compute_mean_std_f(a, x_dim, y_dim, &mean_std);

    for (int i = 0; i < *x_dim; i++)
    {
        for(int j = 0; j < *y_dim; j++)
        {
            result[i * *x_dim + j] = (a[i * *x_dim + j] - mean_std.mean) / mean_std.std;           
        }
    }
}

float elementwise_sum_f(
    const float *a,
    const int *x_dim,
    const int *y_dim
)
{
    float acc = 0.0;
    for (int i = 0; i < *x_dim; i++)
    {
        for(int j = 0; j < *y_dim; j++)
        {
            acc += a[i * *x_dim + j];
        }
    }
    return acc;
}

void elementwise_multiply_f(
    const float *a,
    const float *b,
    const int *x_dim,
    const int *y_dim,
    float *result
)
{
    for (int i = 0; i < *x_dim; i++)
    {
        for(int j = 0; j < *y_dim; j++)
        {
            result[i * *x_dim + j] = a[i * *x_dim + j] * b[i * *x_dim + j];
        }
    }
}

void zero_norm_cross_correlations_f(
    const float *image,
    const int *image_x_dim,
    const int *image_y_dim,
    const float *templ,
    const int *templ_x_dim,
    const int *templ_y_dim,
    float *result
)
{
    // Dimensions
    int o_dim = *image_x_dim - *templ_y_dim + 1;
    int i_dim = *image_y_dim - *templ_y_dim + 1;
    int image_len = *image_x_dim * *image_y_dim;
    int templ_len = *templ_x_dim * *templ_y_dim;
    int o_upper_bound = 0;
    int i_upper_bound = 0;

    // Local slices and z-score containers
    float *image_slice, *templ_z, *image_z, *to_sum;
    image_slice = (float*) malloc(templ_len * sizeof(float));
    templ_z = (float*) malloc(templ_len * sizeof(float));
    image_z = (float*) malloc(templ_len * sizeof(float));
    to_sum = (float*) malloc(templ_len * sizeof(float));
    
    // Z-score normalize template
    zscore_normalize_f(templ, templ_x_dim, templ_y_dim, templ_z);

    for (int o = 0; o < o_dim; o++)
    {
        for (int i = 0; i < i_dim; i++)
        {
            // Build the image slice
            for (int is = 0; is < *templ_x_dim; is++)
            {
                for (int js = 0; js < *templ_y_dim; js++)
                {
                    image_slice[is * *templ_x_dim + js] = image[(o + is) * *templ_x_dim + (i + js)];
                }
            }

            // Z-score normalize image_slice
            zscore_normalize_f(image_slice, templ_x_dim, templ_y_dim, image_z);

            // Compute the ZNCC
            elementwise_multiply_f(image_z, templ_z, templ_x_dim, templ_y_dim, to_sum);
            result[o * o_dim + i] = elementwise_sum_f(to_sum, templ_x_dim, templ_y_dim) / templ_len;
        }
    }

    // Free that memory
    free(image_slice);
    free(templ_z);
    free(image_z);
    free(to_sum);
}