Orientation detection filter

..Orinent2D..

Calculate edge orientation for digital images.

.gitignore

*.o
orient2d

lena.png

Makefile

CC			:= gcc
CXX			:= g++

RM			:= rm
SH			:= bash

SRC			:= orient2d.cpp
OBJS		:= $(patsubst %.cpp, %.o, $(SRC))

CFLAGS		:= -Wall -g -O2 -I/usr/include -I/usr/local/include -I/usr/local/include/opencv4
CXXFLAGS	:= -std=c++11 -Wall -g -O2 -I/usr/include -I/usr/local/include -I/usr/local/include/opencv4
LDFLAGS		:= -L/usr/local/lib/opencv4 -lopencv_core -lopencv_highgui -lopencv_imgproc -lopencv_imgcodecs

PROGRAM		:= orient2d

.PHONY: all
all: $(PROGRAM)

$(OBJS): $(SRC)
	$(CXX) $(CXXFLAGS) -c $<

$(PROGRAM): $(OBJS)
	$(CXX) $(LDFLAGS) -o $(PROGRAM) $^

.PHONY: clean
clean:
	@$(RM) -f *.o $(PROGRAM)

orient2d.cpp

#include <cmath>
#include <string>
#include <opencv2/opencv.hpp>

const double Pi = 4.0 * std::atan(1.0);

inline static double convolute(const cv::Mat &img, int x, int y, const cv::Mat &kernel) {
    const int width = img.cols;
    const int height = img.rows;
    const int rx = (kernel.cols - 1) / 2;
    const int ry = (kernel.rows - 1) / 2;

    double accum = 0.0;
    double sumwgt = 0.0;
    int count = 0;
    for (int v = -ry; v <= ry; v++) {
        for (int u = -rx; u <= rx; u++) {
            int nx = x + u;
            int ny = y + v;
            if (nx < 0 || ny < 0 || nx >= width || ny >= height) continue;

            int fx = u + rx;
            int fy = v + ry;
            double w = kernel.at<float>(fx, fy);
            accum += w * img.at<float>(ny, nx);
            sumwgt += w;
            count += 1;
        }
    }

    return accum;
}

inline static void orientFilter(cv::InputArray src, cv::OutputArray dst,
                                int nFilters = 72, int ksize = 13, double sigma_u = 1.8,
                                double sigma_v = 2.4, double lambda = 4.0) {
    cv::Mat org = src.getMat();
    cv::Mat &output = dst.getMatRef();
        
    if (org.channels() != 3) {
        fprintf(stderr, "Input image must have colors!\n");
        std::exit(1);
    }
        
    if (org.depth() != CV_32F) {
        fprintf(stderr, "Input image must be 32-bit floating point type!\n");
        std::exit(1);
    }

    cv::Mat gray;
    cv::cvtColor(org, gray, cv::COLOR_BGR2GRAY);
    
    const int width  = gray.cols;
    const int height = gray.rows;
    const double stepSize = Pi / (double)nFilters;
    
    // Compute filter bank
    std::vector<cv::Mat> filters;
    const int rx = (ksize - 1) / 2;
    const int ry = (ksize - 1) / 2;
    for (int i = 0; i < nFilters; i++) {
        cv::Mat filter(cv::Size(ksize, ksize), CV_32FC1, 0.0);
        for (int v = -ry; v <= ry; v++) {
            for (int u = -rx; u <= rx; u++) {
                if (u * u + v * v > rx * ry) continue;

                const double theta = i * stepSize;
                double u_tilde = u * std::cos(theta) - v * std::sin(theta);
                double v_tilde = u * std::sin(theta) + v * std::cos(theta);
                double ax = u_tilde / sigma_u;
                double ay = v_tilde / sigma_v;
                double w = std::exp(-0.5 * (ax * ax + ay * ay)) * 
                           std::cos(2.0 * Pi * u_tilde / lambda);

                int fx = u + rx;
                int fy = v + ry;
                filter.at<float>(fx, fy) = w;
            }
        }
        filters.push_back(filter);
    }

    cv::Mat orient = cv::Mat(height, width, CV_32FC1, 0.0);
    cv::Mat confid = cv::Mat(height, width, CV_32FC1, 0.0);
    for (int iter = 0; iter < 3; iter++) {
        // Initialize with zeros
        orient.setTo(cv::Scalar(0.0));
        confid.setTo(cv::Scalar(0.0));

        // Convolution
        orient.forEach<float>([&](float &value, const int p[2]) -> void {
            const int x = p[1];
            const int y = p[0];

            double maxval = -1.0e20;
            int maxidx = 0;
            std::vector<std::pair<double,double>> response;
            for (int i = 0; i < nFilters; i++) {            
                double res = convolute(gray, x, y, filters[i]);
                response.emplace_back(res, i * stepSize);

                if (maxval < res) {
                    maxval = res;
                    maxidx = i;
                }
            }

            double confidence = 0.0;
            for (int i = 0; i < nFilters; i++) {
                double t1 = response[i].second;
                double t2 = response[maxidx].second;
                double dist = std::min(std::abs(t1 - t2), std::min(std::abs(t1 - t2 + Pi), std::abs(t1 - t2 - Pi)));
                double diff = std::min(std::abs(response[i].first - response[maxidx].first), nFilters - std::abs(response[i].first - response[maxidx].first));
                confidence += std::sqrt(dist * diff * diff);
            }

            orient.at<float>(y, x) = response[maxidx].second;
            confid.at<float>(y, x) = confidence;
        });

        // Update confidence and copy to image
        double maxval, minval;
        cv::minMaxIdx(confid, &minval, &maxval);
        confid.forEach<float>([&](float &v, const int p[2]) -> void {
            v = std::pow(v - minval, 1.0) / std::pow(maxval - minval, 1.0);
            v = std::min(v * 2.0, 1.0);
        });
        confid.copyTo(gray);
    }

    // Bilateral filtering
    const double sigma_d = 3.0;
    const double sigma_p = 1.0;
    const double sigma_g = 0.1;
    output = cv::Mat(height, width, CV_32FC3);
    output.forEach<cv::Vec3f>([&](cv::Vec3f &value, const int p[2]) -> void {
        const int x = p[1];
        const int y = p[0];
        double accum = 0.0;
        double sumWgt = 0.0;
        for (int dy = -ry; dy <= ry; dy++) {
            for (int dx = -rx; dx <= rx; dx++) {
                int nx = x + dx;
                int ny = y + dy;
                if (nx >= 0 && ny >= 0 && nx < width && ny < height) {
                    double ds2 = (dx * dx + dy * dy) / (sigma_d * sigma_d);

                    double vp = confid.at<float>(y, x);
                    double vq = confid.at<float>(ny, nx);
                    double dv = vp / (vq + 1.0e-6) / (sigma_p * sigma_p);

                    cv::Vec3f cp = org.at<cv::Vec3f>(y, x);
                    cv::Vec3f cq = org.at<cv::Vec3f>(ny, nx);
                    float dR = cp[0] - cq[0];
                    float dG = cp[1] - cq[1];
                    float dB = cp[2] - cq[2];
                    double dc2 = (dR * dR + dG * dG + dB * dB) / (sigma_g * sigma_g);

                    double w = std::exp(-(ds2 + dv + dc2));
                    accum += w * orient.at<float>(ny, nx);
                    sumWgt += w;
                }
            }
        }
        const double confidence = confid.at<float>(y, x);
        const double theta = sumWgt != 0.0 ? 2.0 * accum / sumWgt : 0.0f;
        value = cv::Vec3f(0.0f, std::cos(theta),  std::sin(theta)) * confidence * 0.5f + cv::Vec3f(0.5f, 0.5f, 0.5f);
    });
}

int main(int argc, char** argv) {
    if (argc <= 1) {
        printf("usage: ./orient2d [image file]\n");
        std::exit(1);
    }

    cv::Mat img = cv::imread(argv[1], cv::IMREAD_COLOR);
    if (img.empty()) {
        fprintf(stderr, "Failed to load file: %s\n", argv[0]);
        std::exit(1);
    }

    img.convertTo(img, CV_32F, 1.0 / 255.0);

    // Filter
    cv::Mat orient;
    orientFilter(img, orient);

    // Save
    cv::Mat output;
    orient.convertTo(output, CV_8U, 255.0);
    cv::imwrite("output.png", output);
}

output.png