Compute Zero-norm Cross-Correlations from Numpy Arrays

zero_norm_cross_correlation.cpp

#include <iostream>
#include <valarray>
#include <sstream>
#include <fstream>
#include <string>
#include <numeric>
#include <cmath>
#include <chrono>

using namespace std;

valarray<float> read_numpy_file(string file_name)
{
    ifstream ifs(file_name);
    size_t x_dim;
    size_t y_dim;
    ifs >> x_dim >> y_dim;
    valarray<float> v(x_dim * y_dim);
    for(size_t i = 0; i < x_dim; i++)
    {
        for(size_t j = 0; j < y_dim; j++)
        {
            ifs >> v[i * x_dim + j];
        }
        
    }
    return v;
}

float compute_mean(
    const valarray<float> &a,
    const size_t &x_dim,
    const size_t &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 / a.size();
}

valarray<float> compute_mean_and_standard_deviation(
        const valarray<float> &a,
        const size_t &x_dim,
        const size_t &y_dim
)
{
    // Compute the mean, store in first entry of mean_and_std
    valarray<float> mean_and_std;
    mean_and_std.resize(2);
    mean_and_std[0] = compute_mean(a, x_dim, y_dim);

    // Compute the absolute square differences
    float acc = 0.0;
    float val = 0.0;
    for(size_t i = 0; i < x_dim; i++)
    {
        for(size_t j = 0; j < y_dim; j++)
        {
            val = abs(a[i * x_dim + j] - mean_and_std[0]);
            acc += val * val;
        }
    }

    // First item in mean_and_std becomes the standard deviation
    mean_and_std[1] = sqrt(acc / a.size());
    return mean_and_std;
}

valarray<float> zscore_normalize(
    const valarray<float> &a,
    const size_t &x_dim,
    const size_t &y_dim
)
{
    valarray<float> output;
    output.resize(x_dim * y_dim);
    valarray<float> mean_and_std = compute_mean_and_standard_deviation(a, x_dim, y_dim);
    for(size_t i = 0; i < x_dim; i++)
    {
        for(size_t j = 0; j < y_dim; j++)
        {
            output[i * x_dim + j] = (a[i * x_dim + j] - mean_and_std[0]) / mean_and_std[1];
        }
    }
    return output;
}

float elementwise_sum(
    const valarray<float> &a,
    const size_t &x_dim,
    const size_t &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;
}

valarray<float> elementwise_multiply(
    const valarray<float> &a,
    const valarray<float> &b,
    const size_t x_dim,
    const size_t y_dim
)
{
    valarray<float> output;
    output.resize(x_dim * y_dim);
    for(size_t i = 0; i < x_dim; i++)
    {
        for(size_t j = 0; j < y_dim; j++)
        {
            output[i * x_dim + j] = a[i * x_dim + j] * b[i * x_dim + j];
        }
    }
    return output;
}

valarray<float> zero_norm_cross_correlations(
    const valarray<float> &image,
    const size_t &image_x_dim,
    const size_t &image_y_dim,
    const valarray<float> &templ,
    const size_t &templ_x_dim,
    const size_t &templ_y_dim
)
{
    size_t o_dim = image_x_dim - templ_x_dim + 1;
    size_t i_dim = image_y_dim - templ_y_dim + 1;
    size_t o_upper_bound, i_upper_bound;

    valarray<float> output;
    output.resize(o_dim * i_dim);
    valarray<float> image_slice;
    image_slice.resize(templ_x_dim * templ_y_dim);

    valarray<float> templ_z = zscore_normalize(templ, templ_x_dim, templ_y_dim);
    valarray<float> image_z;
    image_z.resize(templ_x_dim * templ_y_dim);

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

            image_z = zscore_normalize(image_slice, templ_x_dim, templ_y_dim);

            output[o * o_dim + i] = elementwise_sum(
                elementwise_multiply(
                    image_z,
                    templ_z,
                    templ_x_dim,
                    templ_y_dim
                ),
                templ_x_dim,
                templ_y_dim
            ) / templ_z.size();
        }
    }
    return output;
}

int main(int argc, char const *argv[])
{
    // Image and dimensions
    auto image = read_numpy_file("image.txt");
    size_t image_x_dim = 224;
    size_t image_y_dim = 1294;

    // Template and dimensions
    auto templ = read_numpy_file("template.txt"); 
    size_t templ_o_min = 130;
    size_t templ_o_max = 214;
    size_t templ_i_min = 1005;
    size_t templ_i_max = 1093;
    size_t templ_o_dim = templ_o_max - templ_o_min;
    size_t templ_i_dim = templ_i_max - templ_i_min;
    
    // Additional Image Dimensions
    size_t image_o_start = templ_o_min - 5;
    size_t image_o_finis = templ_o_max + 5;
    size_t image_o_dim = image_o_finis - image_o_start;
    size_t image_i_dim = image_y_dim;

    // Slices
    valarray<float> templ_slice;
    templ_slice.resize(templ_o_dim * templ_i_dim);
    valarray<float> image_slice;
    image_slice.resize(image_o_dim * image_i_dim);

    // Fill templ_slice
    for (size_t o = 0; o < templ_o_dim; o++)
    {
        for (size_t i = 0; i < templ_i_dim; i++)
        {
            templ_slice[o * templ_o_dim + i] = templ[(o + templ_o_min) * templ_o_dim + (i + templ_i_min)];
        }
    }

    // Fill image_slice
    for (size_t o = 0; o < image_o_dim; o++)
    {
        for (size_t i = 0; i < image_i_dim; i++)
        {
            image_slice[o * image_o_dim + i] = image[(o + image_o_start) * image_o_dim + i];
        }
    }

    auto start = chrono::high_resolution_clock::now();
    auto output_ccs = zero_norm_cross_correlations(
        image_slice,
        image_o_dim,
        image_i_dim,
        templ_slice,
        templ_o_dim,
        templ_i_dim
    );
    auto finish = chrono::high_resolution_clock::now();
    cout << chrono::duration_cast<chrono::milliseconds>(finish - start).count() << "\n";
    return 0;
}