Add simple convolution C++ code

Author: gjbex

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.18)
+project(algorithms LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 23)
+set(CMaKE_CXX_STANDARD_REQUIRED YES)
+set(CMAKE_CXX_EXTENSIONS NO)
+
+add_compile_options(-Wall -Wextra -Wpedantic)
+
+add_subdirectory(src)

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/CMakeLists.txt

@@ -0,0 +1,15 @@
+add_subdirectory(convolution)
+
+add_executable(test_convolution.exe
+        test_convolution.cpp)
+target_include_directories(test_convolution.exe
+    PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution)
+target_link_libraries(test_convolution.exe
+    PRIVATE convolution)
+
+add_executable(benchmark_convolution.exe
+        benchmark_convolution.cpp)
+target_include_directories(benchmark_convolution.exe
+    PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution)
+target_link_libraries(benchmark_convolution.exe
+    PRIVATE convolution)

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/benchmark_convolution.cpp

@@ -0,0 +1,57 @@
+#include <convolution.h>
+#include <chrono>
+#include <iostream>
+#include <numeric>
+#include <random>
+
+Matrix create_image(int rows, int cols) {
+    std::mt19937 gen(1234);
+    std::uniform_real_distribution<double> dis(0.0, 1.0);
+    Matrix image(rows, cols);
+    for (int i = 0; i < rows; ++i) {
+        for (int j = 0; j < cols; ++j) {
+            image(i, j) = dis(gen);
+        }
+    }
+    return image;
+}
+
+Matrix create_kernel(int rows, int cols) {
+    Matrix kernel(rows, cols);
+    for (int i = 0; i < rows; ++i) {
+        for (int j = 0; j < cols; ++j) {
+            kernel(i, j) = 1.0/(rows*cols);
+        }
+    }
+    return kernel;
+}
+
+double element_sum(const Matrix& matrix) {
+    return std::accumulate(matrix.data(), matrix.data() + matrix.rows()*matrix.cols(), 0.0);
+}
+
+int main(int argc, char** argv) {
+    int rows = 1000;
+    int cols = 1000;
+    int kernel_rows = 7;
+    int kernel_cols = 7;
+    if (argc > 1) {
+        rows = atoi(argv[1]);
+        cols = atoi(argv[1]);
+    }
+    if (argc > 2) {
+        kernel_rows = atoi(argv[2]);
+        kernel_cols = atoi(argv[2]);
+    }
+    std::cout << "Image size: " << rows << "x" << cols << "\n";
+    std::cout << "Kernel size: " << kernel_rows << "x" << kernel_cols << "\n";
+    Matrix image = create_image(rows, cols);
+    Matrix kernel = create_kernel(kernel_rows, kernel_cols);
+    auto start = std::chrono::high_resolution_clock::now();
+    auto result = convolve(image, kernel);
+    auto end = std::chrono::high_resolution_clock::now();
+    std::chrono::duration<double> diff = end - start;
+    std::cout << "Time: " << diff.count() << " s\n";
+    std::cout << "Sum: " << element_sum(result) << "\n";
+    return 0;
+}

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/convolution/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_library(convolution SHARED
+    convolution.cpp matrices.cpp)

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/convolution/convolution.cpp

@@ -0,0 +1,42 @@
+#include "convolution.h"
+#include <exception>
+
+Matrix convolve(const Matrix& image, const Matrix& kernel) {
+    if (kernel.rows() % 2 != 1 || kernel.cols() % 2 != 1) {
+        throw std::invalid_argument("Only odd dimensions on kernel supported");
+    }
+    /*
+      s_mid and t_mid are number of pixels between the center pixel
+      and the edge, ie for a 5x5 filter they will be 2.
+      
+      The output size is calculated by adding s_mid, t_mid to each
+      side of the dimensions of the input image.
+    */
+    auto s_mid {kernel.rows()/2};
+    auto t_mid {kernel.cols()/2};
+    auto x_max {image.rows() + 2*s_mid};
+    auto y_max {image.cols() + 2*t_mid};
+    // Allocate result image.
+    Matrix new_image(x_max, y_max);
+    // Do convolution
+    for (int x = 0; x < x_max; ++x) {
+        for (int y = 0; y < y_max; ++y) {
+            // Calculate pixel value for h at (x,y). Sum one component
+            // for each pixel (s, t) of the filter kernel.
+            auto s_from {std::max(s_mid - x, -s_mid)};
+            auto s_to {std::min((x_max - x) - s_mid, s_mid + 1)};
+            auto t_from {std::max(t_mid - y, -t_mid)};
+            auto t_to {std::min((y_max - y) - t_mid, t_mid + 1)};
+            double value {0.0};
+            for (int s = s_from; s < s_to; ++s) {
+                for (int t = t_from; t < t_to; ++t) {
+                    auto v {x - s_mid + s};
+                    auto w {y - t_mid + t};
+                    value += kernel(s_mid - s, t_mid - t)*image(v, w);
+                }
+            }
+            new_image(x, y) = value;
+        }
+    }
+    return new_image;
+}

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/convolution/convolution.h

@@ -0,0 +1,8 @@
+#ifndef CONVOLUTION_HDR
+#define CONVOLUTION_HDR
+
+#include "matrices.h"
+
+Matrix convolve(const Matrix& image, const Matrix& kernel);
+
+#endif

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/convolution/matrices.cpp

@@ -0,0 +1,49 @@
+#include "matrices.h"
+
+Matrix::Matrix(const Matrix& other) :
+        rows_(other.rows_), cols_(other.cols_),
+        data_(new double[rows_ * cols_]) {
+    for (int i = 0; i < rows_ * cols_; ++i) {
+        data_[i] = other.data_[i];
+    }
+}
+
+Matrix& Matrix::operator=(const Matrix& other) {
+    if (this != &other) {
+        rows_ = other.rows_;
+        cols_ = other.cols_;
+        data_.reset(new double[rows_ * cols_]);
+        for (int i = 0; i < rows_ * cols_; ++i) {
+            data_[i] = other.data_[i];
+        }
+    }
+    return *this;
+}
+Matrix::Matrix(Matrix&& other) noexcept :
+        rows_{other.rows_}, cols_{other.cols_},
+        data_{std::move(other.data_)} {
+    other.rows_ = 0;
+    other.cols_ = 0;
+}
+
+Matrix& Matrix::operator=(Matrix&& other) noexcept {
+    if (&other != this) {
+        rows_ = other.rows_;
+        cols_ = other.cols_;
+        data_ = std::move(other.data_);
+
+        other.rows_ = 0;
+        other.cols_ = 0;
+    }
+    return *this;
+}
+
+std::ostream& operator<<(std::ostream& os, const Matrix& m) {
+    for (int i = 0; i < m.rows_; ++i) {
+        for (int j = 0; j < m.cols_; ++j) {
+            os << m(i, j) << " ";
+        }
+        os << std::endl;
+    }
+    return os;
+}

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/convolution/matrices.h

@@ -0,0 +1,37 @@
+#ifndef MATRICES_HDR
+#define MATRICES_HDR
+
+#include <iostream>
+#include <memory>
+
+
+struct Matrix {
+    private:
+        int rows_;
+        int cols_;
+        std::unique_ptr<double[]> data_;
+    public:
+        Matrix(int rows, int cols) :
+            rows_(rows), cols_(cols), data_(new double[rows * cols]) {}
+        // copy constructor & assignment operator
+        Matrix(const Matrix& other);
+        Matrix& operator=(const Matrix& other);
+        // move constructor & assignment operator
+        Matrix(Matrix&& other) noexcept;
+        Matrix& operator=(Matrix&& other) noexcept;
+        // matrix indexing by row and column
+        double& operator()(int i, int j) { return data_[i * cols_ + j]; }
+        double operator()(int i, int j) const { return data_[i * cols_ + j]; }
+        // getters for number of rows and columns
+        int rows() const { return rows_; }
+        int cols() const { return cols_; }
+        // accessors for the data
+        double* data() { return data_.get(); }
+        const double* data() const { return data_.get(); }
+        // destructor
+        ~Matrix() = default;
+        // textual representation of the matrix
+        friend std::ostream& operator<<(std::ostream& os, const Matrix& m);
+};
+
+#endif

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/test_convolution.cpp

@@ -0,0 +1,34 @@
+#include <convolution.h>
+#include <iostream>
+
+int main() {
+    // Create a 10x10 image
+    Matrix image(10, 10);
+    for (int i = 0; i < image.rows(); ++i) {
+        for (int j = 0; j < image.cols(); ++j) {
+            image(i, j) = i*image.cols() + j;
+        }
+    }
+
+    // Print the image
+    std::cout << image << std::endl;
+
+    // Create a 3x3 kernel
+    Matrix kernel(3, 3);
+    for (int i = 0; i < kernel.rows(); ++i) {
+        for (int j = 0; j < kernel.cols(); ++j) {
+            kernel(i, j) = 1.0/(kernel.rows()*kernel.cols());
+        }
+    }
+
+    // Print the kernel
+    std::cout << kernel << std::endl;
+
+    // Create a convolution object
+    auto new_image = convolve(image, kernel);
+
+    // Print the result
+    std::cout << new_image << std::endl;
+
+    return 0;
+}