Sobel Edge detection

matrix.h

#include <initializer_list>
#include <vector>
#include <iostream>
#include <cmath>
#include <limits>

template <typename T>
class Matrix
{
public:
    std::vector<T> m_vec;
    const int m_rowNum;

    Matrix(const Matrix& m) : m_rowNum(m.m_rowNum), m_vec(m.m_vec) {};

    Matrix(const int m, const int n) : m_vec(m * n, 0), m_rowNum(n) {};

    Matrix(const std::initializer_list<std::initializer_list<T>>& il) : m_rowNum(il.size())
    {
        int row_s = 0;
        for (const auto& ix : il) {
            int col_s = 0;
            for (const auto& iy : ix) {
                m_vec.push_back(iy);
                col_s++;
            }
            row_s++;
        }
    }

    const int get_rows() const { return m_rowNum; };
    const int get_cols() const { return m_vec.size() / m_rowNum; };
    
    T& operator()(const int col, const int row)
    {
        return el(col, row);
    }

    const T& operator()(const int col, const int row) const
    {
        return el(col, row);
    }

    T& el(const int col, const int row)
    {
        return m_vec[row * get_cols() + col];
    }

    const T& el(const int col, const int row) const
    {
        return m_vec[row * get_cols() + col];
    }

    void el(const int col, const int row, const T val)
    {
        m_vec[row * get_cols() + col] = val;
    }

    const Matrix<T> operator+(const Matrix<T>& other) const
    {
        Matrix<T> result = *this;
        result += other;
        return result;
    }

    Matrix<T>& operator+=(const Matrix<T>& other)
    {
        for (int i = 0; i < m_vec.size(); i++) {
            m_vec[i] += other.m_vec[i];
        }

        return *this;
    }

    template <typename U>
    const Matrix<T> operator*(const U other) const
    {
        Matrix<T> result = *this;
        result *= other;
        return result;
    }

    Matrix<T>& operator*=(const Matrix<T>& other)
    {
        for (int i = 0; i < m_vec.size(); i++) {
            m_vec[i] *= other.m_vec[i];
        }

        return *this;
    }

    Matrix<T>& operator*=(const T other)
    {
        for (int i = 0; i < m_vec.size(); i++) {
            m_vec[i] *= other;
        }

        return *this;
    }

    const Matrix<T> operator/(const Matrix<T>& other) const
    {
        Matrix<T> result = *this;
        result /= other;
        return result;
    }

    Matrix<T>& operator/=(const Matrix<T>& other) const
    {
        for (int i = 0; i < m_vec.size(); i++) {
            m_vec[i] /= other.m_vec[i];
        }
        
        return *this;
    }

    const Matrix<T> operator^(const int val) const
    {
        Matrix<T> result = *this;
        result ^= val;
        return result;
    }

    Matrix& operator^=(const int val)
    {
        T temp;
        for (int i = 0; i < m_vec.size(); i++) {
            temp = m_vec[i];
            for (int j = 1; j < val; j++)
                m_vec[i] *= temp;
        }

        return *this;
    }

    Matrix& operator=(const T val)
    {
        for (int i = 0; i < m_vec.size(); i++) {
            m_vec[i] = val;
        }

        return *this;
    }

    const Matrix<double> sqrt() const
    {
        Matrix<double> result (get_cols(), get_rows());
        for (int i = 0; i < m_vec.size(); i++) {
            result.m_vec[i] = std::sqrt(m_vec[i]);
        }
        return result;
    }

    const Matrix<int> round() const
    {
        Matrix<int> result (get_cols(), get_rows());
        for (int i = 0; i < m_vec.size(); i++) {
            result.m_vec[i] = static_cast<int>(std::round(m_vec[i]));
        }
        return result;
    }

    Matrix& threshold(const T val)
    {
        for (int i = 0; i < m_vec.size(); i++) {
            if (m_vec[i] > val)
                m_vec[i] = val;
        }

        return *this;
    }

    void double_threshold(const T separator, const T val1, const T val2)
    {
        for (int i = 0; i < m_vec.size(); i++) {
            if (m_vec[i] < separator) {
                m_vec[i] = val1;
            }
            else {
                m_vec[i] = val2;
            }
        }
    }

    const T max()
    {
        T result = std::numeric_limits<T>::min();
        for (int i = 0; i < m_vec.size(); i++) {
            if (m_vec[i] > result)
                result = m_vec[i];
        }
        return result;
    }

    Matrix& normalize(const T val)
    {
        T max_val = max();
        for (int i = 0; i < m_vec.size(); i++) {
            m_vec[i] *= val / max_val;
        }

        return *this;
    }

    template <typename U>
    const Matrix<double> atan2(const Matrix<U>& x) const
    {
        Matrix<double> result (get_cols(), get_rows());
        for (int i = 0; i < m_vec.size(); i++) {
            result.m_vec[i] = std::atan2(m_vec[i], x.m_vec[i]);
        }
        return result;
    }

    template <typename U>
    const Matrix<U> cast()
    {
        Matrix<U> result (get_cols(), get_rows());
        for (int i = 0; i < m_vec.size(); i++) {
            result.m_vec[i] = static_cast<U>(m_vec[i]);
        }
        return result;
    }

    const T sum() const
    {
        T result = 0;
        for (auto& el : m_vec) {
            result += el;
        }
        return result;
    }

    template <typename U>
    friend std::ostream& operator<<(std::ostream&, const Matrix<U>&);  
};

template <typename T>
std::ostream& operator<<(std::ostream& os, const Matrix<T>& m) {
    for (int i = 0; i < m.get_rows(); i++) {
        os << "[ ";
        for (int j = 0; j < m.get_cols(); j++) {
            os << m(j, i) << " ";
        }
        os << "]" << std::endl;
    }

    return os;
}

sobel.cpp

#include "matrix.h"
#include <iostream>
#include <fstream>
#include <string>
#include <cstring>
#include <cstdint>
#include <cmath>

const double PI = 3.1415926535;
const int HIGH = 45;
const int LOW = 10;

template <typename T, typename U>
Matrix<T> convolve2d(const Matrix<U>& modifier, const Matrix<T>& origin, bool normalize_modifier = false)
{
    Matrix<T> result (origin.get_cols(), origin.get_rows());
    result = 0;

    const int modifier_width = modifier.get_cols();
    const int modifier_height = modifier.get_rows();
    const int modifier_sum = modifier.sum();

    int origin_val, origin_x, origin_y;
    for (int x = 1; x < origin.get_cols() - 1; x++) {
        for (int y = 1; y < origin.get_rows() - 1; y++) {
            int cell_result = 0;
            for (int xi = 0; xi < modifier_width; xi++) {
                for (int yi = 0; yi < modifier_height; yi++) {
                    origin_x = x + xi - modifier_width / 2;
                    origin_y = y + yi - modifier_height / 2;
                    if (origin_x < 0 || origin_x >= origin.get_cols() || origin_y < 0 || origin_y >= origin.get_rows())
                        origin_val = 0;
                    else
                        origin_val = origin(origin_x, origin_y);
                    cell_result += modifier(xi, yi) * origin_val;
                }
            }

            if (normalize_modifier) {
                cell_result /= modifier_sum;
            }

            result(x, y) = cell_result;
        }
    }
    
    return result;
}

Matrix<int>& read_grayscale_ppm(const char* filename)
{
    std::string line;
    std::ifstream is (filename);
    int width = 0, height = 0;
    int rgb_ind = 0;
    int r, g, b;
    int x = 0, y = 0;
    int* dest;

    std::getline(is, line); // P3
    std::getline(is, line); // dimensions

    std::sscanf(line.c_str(), "%d %d", &width, &height);

    Matrix<int>* image = new Matrix<int>(width, height);

    std::getline(is, line); // max value

    while (std::getline(is, line)) {
        const char* cline = line.c_str();
        const char* token = std::strtok((char*) cline, " ");

        while (token != NULL && *token != '\r') {
            if (rgb_ind == 0) {
                dest = &r;
                rgb_ind++;
            }
            else if (rgb_ind == 1) {
                dest = &g;
                rgb_ind++;
            }
            else if (rgb_ind == 2) {
                dest = &b;
                rgb_ind = 0;
            }

            std::sscanf(token, "%d", dest);

            if (rgb_ind == 0) {
                image->el(x, y, static_cast<int>(std::round(0.3 * r + 0.59 * g + 0.11 * b)));
                x++;
                if (x == width) {
                    x = 0;
                    y++;
                }
            }
                
            token = std::strtok(NULL, " ");
        }
    }
    is.close();
    return *image;
}

int write_grayscale_ppm(const Matrix<int>& m, const char* filename)
{
    std::ofstream of (filename);
    of << "P2" << std::endl;
    of << m.get_cols() << " " << m.get_rows() << std::endl;
    of << "255" << std::endl;
    for (int y = 0; y < m.get_rows(); y++) {
        for (int x = 0; x < m.get_cols(); x++) {
            of << m(x, y) << " ";
        }
    }
    of << std::endl;
    of.close();
}

const Matrix<int> get_quadrants(const Matrix<double>& theta, const Matrix<double>& Gy)
{
    double temp;
    Matrix<int> result (theta.get_cols(), theta.get_rows());
    result = 0;

    for (int x = 1; x < theta.get_cols() - 1; x++) {
        for (int y = 1; y < theta.get_rows() - 1; y++) {
            temp = theta(x, y);
            if (Gy(x, y) < 0)
                temp += PI;

            result(x, y) = static_cast<int>(std::round(4 * temp / PI));
            if (result(x, y) == 4)
                result(x, y) = 0;
        }
    }

    return result;
}

const Matrix<uint8_t> preprocess_angles(const Matrix<double>& G, const Matrix<int>& D)
{
    double cell, cell1, cell2;
    Matrix<uint8_t> result (G.get_cols(), G.get_rows());
    result = 0;

    for (int x = 1; x < G.get_cols() - 1; x++) {
        for (int y = 1; y < G.get_rows() - 1; y++) {
            cell = G(x, y);

            switch (D(x, y)) {
                case 0:
                    cell1 = G(x - 1, y);
                    cell2 = G(x + 1, y);
                    break;
                case 3:
                    cell1 = G(x - 1, y - 1);
                    cell2 = G(x + 1, y + 1);
                    break;
                case 2:
                    cell1 = G(x, y - 1);
                    cell2 = G(x, y + 1);
                    break;
                case 1:
                    cell1 = G(x + 1, y - 1);
                    cell2 = G(x - 1, y + 1);
                    break;
                default:
                    cell1 = 0;
                    cell2 = 0;
                    break;
            }

            if (cell > cell1 && cell > cell2)
                result(x, y) = 1;
        }
    }

    return result;
}

void fix_cells_at(const Matrix<double>& G, const Matrix<uint8_t>& processed_quadrants, Matrix<uint8_t>& visited, Matrix<int>& result, int x, int y)
{
    if (visited(x, y))
        return;
    visited(x, y) = 1;

    if (y > 0) {
        std::vector<std::pair<int, int>> check = {{x - 1, y}, {x + 1, y}, {x, y - 1}, {x, y + 1}};
        for (auto& p : check) {
            int xi = p.first;
            int yi = p.second;

            if (processed_quadrants(xi, yi) == 1 && G(xi, yi) > LOW) {
                result(xi, yi) = 255;
                fix_cells_at(G, processed_quadrants, visited, result, xi, yi);
            }
        }
    }
}


const Matrix<int> fix_cells(const Matrix<double>& G, const Matrix<uint8_t>& processed_quadrants)
{
    Matrix<uint8_t> visited (G.get_cols(), G.get_rows());
    visited = 0;
    Matrix<int> result (G.get_cols(), G.get_rows());
    result = 0;

    for (int x = 0; x < G.get_cols(); x++) {
        for (int y = 0; y < G.get_rows(); y++) {
            if (processed_quadrants(x, y) == 0)
                continue;

            if (G(x, y) > HIGH) {
                result(x, y) = 255;
                fix_cells_at(G, processed_quadrants, visited, result, x, y);
            }
        }
    }

    return result;
}

    
int main(const int argc, const char** argv)
{
    Matrix<int> Gx {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}};
    Matrix<int> Gy {{1, 2, 1}, {0, 0, 0}, {-1, -2, -1}};
    Matrix<int> gaussian {{1, 2, 1}, {2, 4, 2}, {1, 2, 1}};

    Matrix<int>& in = read_grayscale_ppm(argv[1]);

    Matrix<double> blurred = convolve2d(gaussian, in.cast<double>(), true);
    delete &in;

    Matrix<double> convolved_gx = convolve2d(Gx, blurred);
    Matrix<double> convolved_gy = convolve2d(Gy, blurred);

    Matrix<double> G = ((convolved_gx ^ 2) + (convolved_gy ^ 2)).sqrt();
    G.normalize(255);

    Matrix<double> theta = convolved_gy.atan2(convolved_gx);

    Matrix<int> D = get_quadrants(theta, convolved_gy);
    
    write_grayscale_ppm(G.cast<int>(), "G.ppm");
    write_grayscale_ppm(D.cast<int>() * 63, "quadrants.ppm");

    Matrix<uint8_t> processed_quadrants = preprocess_angles(G, D);
    write_grayscale_ppm(processed_quadrants.cast<int>() * 255, "processed.ppm");

    Matrix<int> final_image = fix_cells(G, processed_quadrants);

    write_grayscale_ppm(final_image, "out.ppm");
}