LogPolarFFTTemplate

LogPolarFFTTemplateMatcher/fftm.cpp

#include "opencv2/core.hpp"
#include "opencv2/opencv.hpp"

using namespace std;
using namespace cv;

//----------------------------------------------------------
// Recombinate image quaters
//----------------------------------------------------------
void Recomb(Mat &src, Mat &dst)
{
    int cx = src.cols >> 1;
    int cy = src.rows >> 1;
    Mat tmp;
    tmp.create(src.size(), src.type());
    src(Rect(0, 0, cx, cy)).copyTo(tmp(Rect(cx, cy, cx, cy)));
    src(Rect(cx, cy, cx, cy)).copyTo(tmp(Rect(0, 0, cx, cy)));
    src(Rect(cx, 0, cx, cy)).copyTo(tmp(Rect(0, cy, cx, cy)));
    src(Rect(0, cy, cx, cy)).copyTo(tmp(Rect(cx, 0, cx, cy)));
    dst = tmp;
}
//----------------------------------------------------------
// 2D Forward FFT
//----------------------------------------------------------
void ForwardFFT(Mat &Src, Mat *FImg, bool do_recomb = true)
{
    int M = getOptimalDFTSize(Src.rows);
    int N = getOptimalDFTSize(Src.cols);
    Mat padded;
    copyMakeBorder(Src, padded, 0, M - Src.rows, 0, N - Src.cols, BORDER_CONSTANT, Scalar::all(0));
    Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
    Mat complexImg;
    merge(planes, 2, complexImg);
    dft(complexImg, complexImg);
    split(complexImg, planes);
    planes[0] = planes[0](Rect(0, 0, planes[0].cols & -2, planes[0].rows & -2));
    planes[1] = planes[1](Rect(0, 0, planes[1].cols & -2, planes[1].rows & -2));
    if (do_recomb)
    {
        Recomb(planes[0], planes[0]);
        Recomb(planes[1], planes[1]);
    }
    planes[0] /= float(M*N);
    planes[1] /= float(M*N);
    FImg[0] = planes[0].clone();
    FImg[1] = planes[1].clone();
}
//----------------------------------------------------------
// 2D inverse FFT
//----------------------------------------------------------
void InverseFFT(Mat *FImg, Mat &Dst, bool do_recomb = true)
{
    if (do_recomb)
    {
        Recomb(FImg[0], FImg[0]);
        Recomb(FImg[1], FImg[1]);
    }
    Mat complexImg;
    merge(FImg, 2, complexImg);
    idft(complexImg, complexImg);
    split(complexImg, FImg);
    Dst = FImg[0].clone();
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
void highpass(Size sz, Mat& dst)
{
    Mat a = Mat(sz.height, 1, CV_32FC1);
    Mat b = Mat(1, sz.width, CV_32FC1);

    float step_y = CV_PI / sz.height;
    float val = -CV_PI*0.5;

    for (int i = 0; i < sz.height; ++i)
    {
        a.at<float>(i) = cos(val);
        val += step_y;
    }

    val = -CV_PI*0.5;
    float step_x = CV_PI / sz.width;
    for (int i = 0; i < sz.width; ++i)
    {
        b.at<float>(i) = cos(val);
        val += step_x;
    }

    Mat tmp = a*b;
    dst = (1.0 - tmp).mul(2.0 - tmp);
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
float logpolar(Mat& src, Mat& dst)
{
    float radii = src.cols;
    float angles = src.rows;
    Point2f center(src.cols / 2, src.rows / 2);
    float d = norm(Vec2f(src.cols - center.x, src.rows - center.y));
    float log_base = pow(10.0, log10(d) / radii);
    float d_theta = CV_PI / (float)angles;
    float theta = CV_PI / 2.0;
    float radius = 0;
    Mat map_x(src.size(), CV_32FC1);
    Mat map_y(src.size(), CV_32FC1);
    for (int i = 0; i < angles; ++i)
    {
        for (int j = 0; j < radii; ++j)
        {
            radius = pow(log_base, float(j));
            float x = radius * sin(theta) + center.x;
            float y = radius * cos(theta) + center.y;
            map_x.at<float>(i, j) = x;
            map_y.at<float>(i, j) = y;
        }
        theta += d_theta;
    }
    remap(src, dst, map_x, map_y, CV_INTER_LINEAR, BORDER_CONSTANT, Scalar(0, 0, 0));
    return log_base;
}
//-----------------------------------------------------------------------------------------------------
// As input we need equal sized images, with the same aspect ratio,
// scale difference should not exceed 1.8 times.
//-----------------------------------------------------------------------------------------------------
RotatedRect LogPolarFFTTemplateMatch(Mat& im0, Mat& im1, double canny_threshold1, double canny_threshold2)
{
    // Accept 1 or 3 channel CV_8U, CV_32F or CV_64F images.
    CV_Assert((im0.type() == CV_8UC1) || (im0.type() == CV_8UC3) ||
        (im0.type() == CV_32FC1) || (im0.type() == CV_32FC3) ||
        (im0.type() == CV_64FC1) || (im0.type() == CV_64FC3));

    CV_Assert(im0.rows == im1.rows && im0.cols == im1.cols);

    CV_Assert(im0.channels() == 1 || im0.channels() == 3 || im0.channels() == 4);

    CV_Assert(im1.channels() == 1 || im1.channels() == 3 || im1.channels() == 4);

    Mat im0_tmp = im0.clone();
    Mat im1_tmp = im1.clone();
    if (im0.channels() == 3)
    {
        cvtColor(im0, im0, cv::COLOR_BGR2GRAY);
    }

    if (im0.channels() == 4)
    {
        cvtColor(im0, im0, cv::COLOR_BGRA2GRAY);
    }

    if (im1.channels() == 3)
    {
        cvtColor(im1, im1, cv::COLOR_BGR2GRAY);
    }

    if (im1.channels() == 4)
    {
        cvtColor(im1, im1, cv::COLOR_BGRA2GRAY);
    }

    if (im0.type() == CV_32FC1)
    {
       im0.convertTo(im0, CV_8UC1, 255.0);
    }

    if (im1.type() == CV_32FC1)
    {
       im1.convertTo(im1, CV_8UC1, 255.0);
    }

    if (im0.type() == CV_64FC1)
    {
        im0.convertTo(im0, CV_8UC1, 255.0);
    }

    if (im1.type() == CV_64FC1)
    {
        im1.convertTo(im1, CV_8UC1, 255.0);
    }


    Canny(im0, im0, canny_threshold1, canny_threshold2); // you can change this
    Canny(im1, im1, canny_threshold1, canny_threshold2);
    
    // Ensure both images are of CV_32FC1 type
    im0.convertTo(im0, CV_32FC1, 1.0 / 255.0);
    im1.convertTo(im1, CV_32FC1, 1.0 / 255.0);

    Mat F0[2], F1[2];
    Mat f0, f1;
    ForwardFFT(im0, F0);
    ForwardFFT(im1, F1);
    magnitude(F0[0], F0[1], f0);
    magnitude(F1[0], F1[1], f1);

    // Create filter 
    Mat h;
    highpass(f0.size(), h);

    // Apply it in freq domain
    f0 = f0.mul(h);
    f1 = f1.mul(h);

    float log_base;
    Mat f0lp, f1lp;

    log_base = logpolar(f0, f0lp);
    log_base = logpolar(f1, f1lp);

    // Find rotation and scale
    Point2d rotation_and_scale = cv::phaseCorrelate(f1lp, f0lp);

    float angle = 180.0 * rotation_and_scale.y / f0lp.rows;
    float scale = pow(log_base, rotation_and_scale.x);
    // --------------
    if (scale > 1.8)
    {
        rotation_and_scale = cv::phaseCorrelate(f1lp, f0lp);
        angle = -180.0 * rotation_and_scale.y / f0lp.rows;
        scale = 1.0 / pow(log_base, rotation_and_scale.x);
        if (scale > 1.8)
        {
            cout << "Images are not compatible. Scale change > 1.8" << endl;
            return RotatedRect();
        }
    }
    // --------------
    if (angle < -90.0)
    {
        angle += 180.0;
    }
    else if (angle > 90.0)
    {
        angle -= 180.0;
    }

    // Now rotate and scale fragment back, then find translation
    Mat rot_mat = getRotationMatrix2D(Point(im1.cols / 2, im1.rows / 2), angle, 1.0 / scale);

    // rotate and scale
    Mat im1_rs;
    warpAffine(im1, im1_rs, rot_mat, im1.size());

    // find translation
    Point2d tr = cv::phaseCorrelate(im1_rs, im0);

    // compute rotated rectangle parameters
    RotatedRect rr;
    rr.center = tr + Point2d(im0.cols / 2, im0.rows / 2);
    rr.angle = -angle;
    rr.size.width = im1.cols / scale;
    rr.size.height = im1.rows / scale;

    im0 = im0_tmp.clone();
    im1 = im1_tmp.clone();

    return rr;
}

LogPolarFFTTemplateMatcher/fftm.hpp

#ifndef _OPENCV_fftm_HPP_
#define _OPENCV_fftm_HPP_
#ifdef __cplusplus

#include "opencv2/core.hpp"
#include "opencv2/opencv.hpp"
//-----------------------------------------------------------------------------------------------------
// As input we need equal sized images, with the same aspect ratio,
// scale difference should not exceed 1.8 times.
//-----------------------------------------------------------------------------------------------------
cv::RotatedRect LogPolarFFTTemplateMatch(cv::Mat& im0, cv::Mat& im1, double canny_threshold1=200, double canny_threshold2=100);
#endif
#endif

LogPolarFFTTemplateMatcher/test.cpp

#include "fftm.hpp"
#include "gtest/gtest.h"

using namespace std;
using namespace cv;

//----------------------------------------------
// Check if rotated box contained in rectangular region
//----------------------------------------------
bool boxInRange(cv::Rect r, cv::RotatedRect& rr)
{
    Point2f rect_points[4];
    rr.points(rect_points);
    bool result = true;
    for (int i = 0; i < 4; ++i)
    {
        if (!r.contains(rect_points[i]))
        {
            result = false;
            break;
        }
    }
    return result;
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
void generateRotRectROI(Mat& img, RotatedRect& rect)
{
    if (img.empty())
    {
        cout << "Empty image in generateRotRectROI. " << endl;
        rect = RotatedRect();
        return;
    }

    int w = img.cols;
    int h = img.rows;
    Rect roi(0, 0, w, h);
    RNG rng(cv::getTickCount());
    while (1)
    {
        int x = rng.uniform(w/2-w/10, w/2+w/10);
        int y = rng.uniform(h/2-h/10, h/2+h/10);
        float scale= rng.uniform(0.85f, 1.0f);
        int wr = float(w)*scale;
        int hr = float(h)*scale;
        float ang= rng.uniform(-5.0, 5.0);
        rect = cv::RotatedRect(Point2f(x, y), Size(wr, hr), ang);
        if (boxInRange(roi, rect))
        {
            break;
        }
    }
}
//----------------------------------------------------------
//
//----------------------------------------------------------
void getQuadrangleSubPix_8u32f_CnR(const uchar* src, size_t src_step, Size src_size,
    float* dst, size_t dst_step, Size win_size,
    const double *matrix, int cn)
{
    int x, y, k;
    double A11 = matrix[0], A12 = matrix[1], A13 = matrix[2];
    double A21 = matrix[3], A22 = matrix[4], A23 = matrix[5];

    src_step /= sizeof(src[0]);
    dst_step /= sizeof(dst[0]);

    for (y = 0; y < win_size.height; y++, dst += dst_step)
    {
        double xs = A12*y + A13;
        double ys = A22*y + A23;
        double xe = A11*(win_size.width - 1) + A12*y + A13;
        double ye = A21*(win_size.width - 1) + A22*y + A23;

        if ((unsigned)(cvFloor(xs) - 1) < (unsigned)(src_size.width - 3) &&
            (unsigned)(cvFloor(ys) - 1) < (unsigned)(src_size.height - 3) &&
            (unsigned)(cvFloor(xe) - 1) < (unsigned)(src_size.width - 3) &&
            (unsigned)(cvFloor(ye) - 1) < (unsigned)(src_size.height - 3))
        {
            for (x = 0; x < win_size.width; x++)
            {
                int ixs = cvFloor(xs);
                int iys = cvFloor(ys);
                const uchar *ptr = src + src_step*iys;
                float a = (float)(xs - ixs), b = (float)(ys - iys), a1 = 1.f - a, b1 = 1.f - b;
                float w00 = a1*b1, w01 = a*b1, w10 = a1*b, w11 = a*b;
                xs += A11;
                ys += A21;

                if (cn == 1)
                {
                    ptr += ixs;
                    dst[x] = ptr[0] * w00 + ptr[1] * w01 + ptr[src_step] * w10 + ptr[src_step + 1] * w11;
                }
                else if (cn == 3)
                {
                    ptr += ixs * 3;
                    float t0 = ptr[0] * w00 + ptr[3] * w01 + ptr[src_step] * w10 + ptr[src_step + 3] * w11;
                    float t1 = ptr[1] * w00 + ptr[4] * w01 + ptr[src_step + 1] * w10 + ptr[src_step + 4] * w11;
                    float t2 = ptr[2] * w00 + ptr[5] * w01 + ptr[src_step + 2] * w10 + ptr[src_step + 5] * w11;

                    dst[x * 3] = t0;
                    dst[x * 3 + 1] = t1;
                    dst[x * 3 + 2] = t2;
                }
                else
                {
                    ptr += ixs*cn;
                    for (k = 0; k < cn; k++)
                        dst[x*cn + k] = ptr[k] * w00 + ptr[k + cn] * w01 +
                        ptr[src_step + k] * w10 + ptr[src_step + k + cn] * w11;
                }
            }
        }
        else
        {
            for (x = 0; x < win_size.width; x++)
            {
                int ixs = cvFloor(xs), iys = cvFloor(ys);
                float a = (float)(xs - ixs), b = (float)(ys - iys), a1 = 1.f - a, b1 = 1.f - b;
                float w00 = a1*b1, w01 = a*b1, w10 = a1*b, w11 = a*b;
                const uchar *ptr0, *ptr1;
                xs += A11; ys += A21;

                if ((unsigned)iys < (unsigned)(src_size.height - 1))
                    ptr0 = src + src_step*iys, ptr1 = ptr0 + src_step;
                else
                    ptr0 = ptr1 = src + (iys < 0 ? 0 : src_size.height - 1)*src_step;

                if ((unsigned)ixs < (unsigned)(src_size.width - 1))
                {
                    ptr0 += ixs*cn; ptr1 += ixs*cn;
                    for (k = 0; k < cn; k++)
                        dst[x*cn + k] = ptr0[k] * w00 + ptr0[k + cn] * w01 + ptr1[k] * w10 + ptr1[k + cn] * w11;
                }
                else
                {
                    ixs = ixs < 0 ? 0 : src_size.width - 1;
                    ptr0 += ixs*cn; ptr1 += ixs*cn;
                    for (k = 0; k < cn; k++)
                        dst[x*cn + k] = ptr0[k] * b1 + ptr1[k] * b;
                }
            }
        }
    }
}

//----------------------------------------------------------
// 
//----------------------------------------------------------
void myGetQuadrangleSubPix(const Mat& src, Mat& dst, Mat& m)
{
    CV_Assert(src.channels() == dst.channels());

    cv::Size win_size = dst.size();
    double matrix[6];
    cv::Mat M(2, 3, CV_64F, matrix);
    m.convertTo(M, CV_64F);
    double dx = (win_size.width - 1)*0.5;
    double dy = (win_size.height - 1)*0.5;
    matrix[2] -= matrix[0] * dx + matrix[1] * dy;
    matrix[5] -= matrix[3] * dx + matrix[4] * dy;

    if (src.depth() == CV_8U && dst.depth() == CV_32F)
        getQuadrangleSubPix_8u32f_CnR(src.data, src.step, src.size(),
        (float*)dst.data, dst.step, dst.size(),
            matrix, src.channels());
    else
    {
        CV_Assert(src.depth() == dst.depth());
        cv::warpAffine(src, dst, M, dst.size(),
            cv::INTER_LINEAR + cv::WARP_INVERSE_MAP,
            cv::BORDER_REPLICATE);
    }
}
//----------------------------------------------------------
// 
//----------------------------------------------------------
void getRotRectImg(cv::RotatedRect rr, Mat &img, Mat& dst)
{
    Mat m(2, 3, CV_64FC1);
    float ang = rr.angle*CV_PI / 180.0;
    m.at<double>(0, 0) = cos(ang);
    m.at<double>(1, 0) = sin(ang);
    m.at<double>(0, 1) = -sin(ang);
    m.at<double>(1, 1) = cos(ang);
    m.at<double>(0, 2) = rr.center.x;
    m.at<double>(1, 2) = rr.center.y;
    myGetQuadrangleSubPix(img, dst, m);
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
float rrDist(RotatedRect r1, RotatedRect r2)
{
    return norm(r1.center - r2.center) + fabs(r1.size.width-r2.size.width) + fabs(r1.size.height - r2.size.height) + fabs(r1.angle-r2.angle);
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
TEST(imgProc_LogPolarFFTTemplateMatch, resultTest)
{
    //cvtest::TS::ptr()->get_data_path() + "denoising/lena_orig.png";
    Mat test_img1 = imread("lena_orig.png", 0);
    
    EXPECT_EQ(false, test_img1.empty());
    
    if (test_img1.empty())
    {
        cout << "Error loading imput image. " << endl;
        return;
    }
    
    Mat test_img2;
    RotatedRect rect;
    generateRotRectROI(test_img1, rect);

    EXPECT_NE(0, rect.size.area());

    test_img2 = Mat(rect.size, test_img1.type());
    getRotRectImg(rect, test_img1, test_img2);
    resize(test_img2, test_img2, test_img1.size());
    RotatedRect rr = LogPolarFFTTemplateMatch(test_img1, test_img2, 200,100);

    EXPECT_NE(0, rr.size.area());

    Point2f rect_points[4];
    rr.points(rect_points);
    for (int j = 0; j < 4; j++)
    {
        line(test_img1, rect_points[j], rect_points[(j + 1) % 4], Scalar(1, 0, 0), 2, CV_AA);
    }
    
    float dist = rrDist(rr, rect);
    EXPECT_LE(dist, 8);
}

//----------------------------------------------------------
// 
//----------------------------------------------------------
int main(int argc, char* argv[])
{
    testing::InitGoogleTest(&argc,argv);
    return RUN_ALL_TESTS();
}

LogPolarFFTTemplateMatcher/example.cpp

#include "fftm.hpp"
using namespace std;
using namespace cv;
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
int main(int argc, unsigned int** argv)
{
    Mat im0 = imread("cat.png", 1);
    Mat im1 = imread("cat_part.png", 1);

    imshow("im1", im1);
    imshow("im0", im0);

    // As input we need equal sized images, with the same aspect ratio,
    // scale difference should not exceed 1.8 times.

    RotatedRect rr = LogPolarFFTTemplateMatch(im0, im1,200,100);

    // Plot rotated rectangle, to check result correctness
    Point2f rect_points[4];
    rr.points(rect_points);
    for (int j = 0; j < 4; j++)
    {
        line(im0, rect_points[j], rect_points[(j + 1) % 4], Scalar(255, 0, 0), 2, CV_AA);
    }

    imshow("result", im0);

    waitKey();
}