Commits

N.Y. Lin committed 853a89f Merge

Merge: Warping and net programming examples

  • Participants
  • Parent commits f08bf92, 8634020

Comments (0)

Files changed (6)

image/npr/kuwahara_filter/Readme.txt

+
+附带的程序文件名为filter.exe,可以设置的参数有分割数N、G_(σ_s )和G_(σ_r ),
+默认为N=4,G_(σ_s )=1.0,G_(σ_r )=3.0。使用例子如下:
+    filter.exe  lion.png  8  1.0  3.0
+输出结果保存到当前目录,名字为filter_result.png。

image/npr/kuwahara_filter/images/lion.png

Added
New image

image/npr/kuwahara_filter/images/lion_proposed.PNG

Added
New image

image/npr/kuwahara_filter/src/filter.cpp

+// ======================================================================================
+// File         : filter.cpp
+// Author       : maxint <lnychina@gmail.com> 
+// Last Change  : 08/21/2010 | 21:52:45 PM | Saturday,August
+// Description  : CPU implementation of Anisotropic Kuwahara filtering
+//                Paper Home: http://www.kyprianidis.com/pg2009.html
+// ======================================================================================
+
+#include "filter.h"
+#include <highgui.h>
+#include <cassert>
+#include <sstream>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <ctime>
+
+using namespace cv;
+using namespace std;
+
+static clock_t t_cal_mean_and_std2 = 0;
+
+//CFilter::CFilter(void)
+//{
+//}
+
+CFilter::~CFilter(void)
+{
+    _Mdest.release();
+    _Msrc.release();
+    _Mest.release();
+}
+
+void CFilter::init(float s_r, float s_s, int N)
+{
+    clock_t tt = clock();
+
+    _segma_r = s_r;
+    _segma_s = s_s;
+    _N = N;
+    int h = int(_segma_r * 2);
+    _Ksize = h;
+    //_Ksize = 2*h+1;
+
+    Mat KK(2*_Ksize, 2*_Ksize, CV_64FC1, 0.0);
+    _meanK.resize(N);
+    _segma2K.resize(N);
+
+    vector<Point> pts;
+    pts.push_back(Point2d( _Ksize, _Ksize-1));
+    pts.push_back(Point2d(2*_Ksize-1, _Ksize-1-int(tan(CV_PI/N)*_Ksize) ) );
+    pts.push_back(Point2d(2*_Ksize-1, _Ksize+int(tan(CV_PI/N)*_Ksize) ) );
+    pts.push_back(Point2d(_Ksize, _Ksize));
+    fillConvexPoly(KK, &pts[0], 4, Scalar(1.0f));
+
+    Matd Mg = createGaussianMatrix(_Ksize, _segma_r);
+    GaussianBlur(KK, KK, Size(), _segma_s);
+    multiply(KK, Mg, KK);
+    normalizeMatrix(KK);
+    Mat K;
+    KK.convertTo(K, CV_32F);   // ��ΪwarpAffineֻ֧��float
+
+    _MKi.resize(_N);
+    int i;
+    for (i=0; i<_N; ++i)
+    {
+        Mat rotateM = getRotationMatrix2D(Point2f(_Ksize-0.5f, _Ksize-0.5f), i * 360/_N, 1);
+        Mat dest;
+        warpAffine(K, dest, rotateM, K.size());
+        _MKi[i] = dest;
+    }
+
+    cout << "CFilter::init()\n\tFinished in " << clock()-tt << " ms.\n";
+}
+
+bool CFilter::calculateStuctureTenser()
+{
+    clock_t tt = clock();
+    int x, y;
+
+    //////////////////////////////////////////////////////////////////////////
+    // Step 1: ������
+    // x, y derivative
+    Mat_<Vec3d> img_dx, img_dy;
+    Sobel(_Msrc, img_dx, CV_64F, 1, 0, CV_SCHARR, 1./16, 0, BORDER_REFLECT);
+    Sobel(_Msrc, img_dy, CV_64F, 0, 1, CV_SCHARR, 1./16, 0, BORDER_REFLECT);
+
+
+    //////////////////////////////////////////////////////////////////////////
+    // Step 2: ����ṹ����
+    // 1 -> dxx, 2 -> dxy, 3 -> dyy
+    Mat_<Vec3d> tensor(_Msrc.size());
+    for (y=0; y<_Msrc.rows; ++y)
+    {
+        for (x=0; x<_Msrc.cols; ++x)
+        {
+            tensor(y, x)[0] = img_dx(y, x).dot(img_dx(y, x));
+            tensor(y, x)[1] = img_dx(y, x).dot(img_dy(y, x));
+            tensor(y, x)[2] = img_dy(y, x).dot(img_dy(y, x));
+        }
+    }
+    // ���Խṹ����
+    Mat_<Vec3d> tensorBlur;
+    GaussianBlur(tensor, tensorBlur, Size(), 2.0, 0, BORDER_REFLECT);
+
+
+    //////////////////////////////////////////////////////////////////////////
+    // Step 3:��ͨ�����������������ط���͸������Գ̶�
+    // 1 -> local orientation, 2 -> anisotropy
+    _Mest.create(tensor.size());
+    Mat M(2, 2, CV_64FC1);
+    Mat_<double> D, V;
+    for (y=0; y<_Msrc.rows; ++y)
+    {
+        for (x=0; x<_Msrc.cols; ++x)
+        {
+            M = (Mat_<double>(2, 2) << tensorBlur(y,x)[0], tensorBlur(y,x)[1], tensorBlur(y,x)[1], tensorBlur(y,x)[2]);
+            if (!eigen(M, D, V)) return false;
+            _Mest(y, x)[0] = std::atan2(V(1,1), V(1,0));
+
+            if ( (_Mest(y, x)[1] = (D(0,0)-D(1,0)) / (std::abs(D(0,0)+D(1,0))+0.001)) >1 ) return false;
+            //double tmp4 = est(y, x)[1];
+        }
+    }
+
+    cout << "CFilter::calculateStuctureTenser()\n\tFinished in " << clock()-tt << " ms.\n";
+
+    return true;
+}
+
+template<typename T1, int n> static inline
+Vec<T1, n> mul(const Vec<T1, n>& a, const Vec<T1, n>& b)
+{
+    Vec<T1, n> c;
+    for (int i=0; i<n; ++i)
+    {
+        c.val[i] = saturate_cast<T1>(a.val[i]*b.val[i]);
+    }
+    return c;
+}
+
+
+void CFilter::apply(const Mat &src, Mat &dest, float alpha/* = 1.0*/, float q /*= 8.0*/ )
+{
+    clock_t tt = clock();
+
+    _Msrc = src;
+
+    if (!calculateStuctureTenser())
+        return;
+
+    dest.create(src.size(), src.type());
+    _Mdest = dest;
+    int i, j;
+    Mat transM(2, 3, CV_64FC1);
+    Mat roi(2*_Ksize, 2*_Ksize, src.type());        // �洢��ԭͼ��任����׼�ռ��ϵ�ͼ���
+    double alphaX, radian;
+    Vec2d estX;
+    for (i=0; i<_Msrc.rows; ++i)
+    {
+        for (j=0; j<_Msrc.cols; ++j)
+        {
+            estX = _Mest.at<Vec2d>(i, j);
+            alphaX = alpha / (alpha+estX[1]);    // anisotropic
+            radian = estX[0];     // orientation
+            transM = getTransformMatrix2D( Point2f(float(j)+0.5f, float(i)+0.5f),
+                Point2f(float(_Ksize), float(_Ksize)), radian, alphaX, 1./alphaX);
+
+            // �����׼�ռ��ϵ�ͼ���
+            warpAffine(src, roi, transM, roi.size(), INTER_LINEAR, BORDER_REFLECT);
+
+            // ����ƽ��ֵ�ͷ���
+            calMeanAndStd2(roi);
+
+            double sum_a_m, sum_a;
+            Vec3b pi;
+            for (int ii=0; ii<3; ++ii)
+            {
+                sum_a_m = sum_a = 0;
+                for (int n=0; n<_N; ++n)
+                {
+                    _segma2K[n][ii] = 1/(1+pow(_segma2K[n][ii], (double)q/2.0f));
+                    sum_a += _segma2K[n][ii];
+                    sum_a_m += _segma2K[n][ii] * _meanK[n][ii];
+                }
+                pi[ii] = saturate_cast<uchar>(sum_a_m / sum_a);
+                //cout << i*_Msrc.cols+j << ":" << (int)pi[0] << "," << (int)pi[1] << "," << (int)pi[2] << endl;
+            }
+            dest.at<Vec3b>(i, j) = pi;
+        }//forj
+    }//fori
+
+    cout << "CFilter::apply()\n\tFinished in " << clock()-tt << " ms.\n";
+    cout << "CFilter::calMeanAndStd2()\n\tFinished accumulatively in " << t_cal_mean_and_std2 << " ms.\n";
+    cout << "\n== Anisotropic Kuwahara Filtering parameters:\n"
+        << "   Image Size:\t" << _Msrc.rows << "(rows) X " << _Msrc.cols << "(cols)\n"
+        << "   Structure Tenser segma:\t0.2 (can't custom now)\n"
+        << "   Fliter segma_r:\t" << _segma_r << "\n"
+        << "   Filter segma_s:\t" << _segma_s << "\n"
+        << "   Number of sections N:\t" << _N << "\n";
+
+}
+
+void CFilter::calMeanAndStd2(const Mat &roi)
+{
+    clock_t tt = clock();
+
+    _meanK.assign(_N, Vec3d(0,0,0));
+    _segma2K.assign(_N, Vec3d(0,0,0));
+    int sqrf;
+    byte bf;
+    for (int y=0; y<2*_Ksize; ++y)
+    {
+        for (int x=0; x<2*_Ksize; ++x)
+        {
+            for (int i=0; i<3; ++i)
+            {
+                bf = roi.at<Vec3b>(y, x)[i];
+                sqrf = bf * bf;
+                for (int k=0; k<_N; ++k)
+                {
+                    _meanK[k][i] += bf * _MKi[k](y, x);
+                    _segma2K[k][i] += sqrf * _MKi[k](y, x);
+                }
+            }
+        }
+    }
+    for (int i=0; i<3; ++i)
+    {
+        for (int k=0; k<_N; ++k)
+        {
+            _segma2K[k][i] -= _meanK[k][i] * _meanK[k][i];
+        }
+    }
+
+    //s2 -= mul(m, m);
+    t_cal_mean_and_std2 += clock()-tt;
+}
+
+Mat_<double> CFilter::createGaussianMatrix(int radius, double sigma1)
+{
+    Matd kernel(2*radius, 2*radius);
+    Matd cd(radius, 1);
+    double radiusX = cvRound(sigma1*3);
+    double scale2X = -0.5/(sigma1*sigma1);
+    double sum = 0;
+
+    int i, j;
+    for (i=0; i<radius; ++i)
+    {
+        double x = (i+1) * radiusX / radius;
+        double t = std::exp(scale2X*x*x);
+        cd(i, 0) = t;
+        sum += t;
+    }
+    sum = 1./sum;
+    for (i=0; i<radius; ++i)
+    {
+        cd(i, 0) *= sum;
+    }
+    // copy to kernel
+    for (i=0; i<radius; ++i)
+    {
+        for (j=0; j<radius; ++j)
+        {
+            kernel(i, radius+j) = cd(radius-1-i, 0) * cd(j, 0);
+            // copy to other 3 points
+            kernel(i, radius-1-j) = kernel(i, radius+j);
+            kernel(2*radius-1-i, radius-1-j) = kernel(i, radius+j);
+            kernel(2*radius-1-i, radius+j) = kernel(i, radius+j);
+        }
+    }
+    return kernel;
+}
+
+void CFilter::printMatrix(const Mat &m, const char* file)
+{
+    std::ofstream of;
+    of.open(file);
+    of.precision(6);
+    of.width(3);
+    for (int i=0; i<m.rows; ++i)
+    {
+        for (int j=0; j<m.cols; ++j)
+        {
+            of << m.at<float>(i, j) << "\t";
+        }
+        of << std::endl;
+    }
+    of.close();
+}
+
+void CFilter::normalizeMatrix(const Mat &m)
+{
+    double sum = 0;
+    MatConstIterator_<double> it = m.begin<double>(), it_end = m.end<double>();
+    for (; it!=it_end; ++it)
+        sum += *it;
+    m /= sum;
+}
+
+//  [12/2/2009 maxint]
+// offsetTranMatrix * scaleMatrix * rotateMatrix * transformMatrix * X = X'
+// X = (x, y, 1)'
+Mat CFilter::getTransformMatrix2D( Point2f center, Point2f offset, double radian, double scale1, double scale2 )
+{
+    double alpha = cos(radian);
+    double beta = sin(radian);
+
+    Mat M(2, 3, CV_64F);
+    double* m = (double*)M.data;
+
+    m[0] = scale1*alpha;
+    m[1] = scale1*beta;
+    m[2] = (-alpha*center.x - beta*center.y)*scale1 + offset.x;
+    m[3] = -beta*scale2;
+    m[4] = alpha*scale2;
+    m[5] = (beta*center.x - alpha*center.y)*scale2 + offset.y;
+
+    return M;
+}

image/npr/kuwahara_filter/src/filter.h

+// ======================================================================================
+// File         : filter.h
+// Author       : maxint <lnychina@gmail.com> 
+// Last Change  : 08/21/2010 | 21:52:42 PM | Saturday,August
+// Description  : CPU implementation of Anisotropic Kuwahara filtering
+//                Paper Home: http://www.kyprianidis.com/pg2009.html
+// ======================================================================================
+
+#ifndef _FILTER_H_
+#define _FILTER_H_
+
+#include <cv.h>
+using namespace cv;
+
+class CFilter
+{
+private:
+    typedef Mat_<double> Matd;
+    typedef Mat_<float> Matf;
+    typedef vector<Matf> MatVec;
+
+public:
+    //CFilter(const Mat &src, Mat &dest, const Mat &est, double s_r, double s_s, int N);
+    ~CFilter(void);
+
+public:
+    // Ԥ��������˲��ֿ�
+    void init(float s_r, float s_s, int N);
+
+    // Ӧ��Kuwahara �˲�
+    void apply(const Mat &src, Mat &dest, float alpha = 1.0, float q = 8.0);
+
+private:
+    // ����ṹ����������ÿ�����ص�ı߽緽��͸������Գ̶�
+    bool calculateStuctureTenser();
+
+    // ����һ����˹���󣬰뾶Ϊradius����׼��Ϊsigmal
+    Matd createGaussianMatrix(int radius, double sigma1);
+
+    void printMatrix(const Mat &m, const char* file);
+
+    void normalizeMatrix(const Mat &m);
+
+    void calMeanAndStd2(const Mat &roi);
+
+    Mat getTransformMatrix2D( Point2f center, Point2f offset, double radian = 0, double scale1 = 1, double scale2 = 1 );
+
+private:
+    int _N;                 // �˲��˷ֿ�������ȡ4��8��16
+    double _segma_r;        // ��������Ԥ�����˲��ֿ��Gaussian�˲�����
+    double _segma_s;
+    Mat _Msrc;              // ԭͼ��
+    Mat _Mdest;             // Ŀ��ͼ��
+    Mat_<Vec2d> _Mest;      // ���������˲���������
+
+    vector<Vec3d> _meanK;   // �ֿ�ľ�ֵ
+    vector<Vec3d> _segma2K; // �ֿ�ķ���
+
+    MatVec _MKi;            // Ԥ�������˲��ֿ�
+    int _Ksize;             // �˲��ֿ��С
+};
+
+#endif//_FILTER_H_

image/npr/kuwahara_filter/src/main.cpp

+// ======================================================================================
+// File         : main.cpp
+// Author       : maxint <lnychina@gmail.com> 
+// Last Change  : 08/21/2010 | 21:52:25 PM | Saturday,August
+// Description  : CPU implementation of Anisotropic Kuwahara filtering
+//                Paper Home: http://www.kyprianidis.com/pg2009.html
+// ======================================================================================
+
+#include <cv.h>
+#include <highgui.h>
+#include <vector>
+#include <iostream>
+#include "filter.h"
+
+using namespace cv;
+using namespace std;
+
+int main(int argc, char* argv[])
+{
+    vector<char*> varg(argv, argv+argc);
+
+    char imgfile[256] = "images/lion.png";
+    float segma_r = 3.0;
+    float segma_s = 1.0;
+    int N = 4;
+
+    if (varg.size()>1)
+        strcpy(imgfile, varg[1]);
+    else
+        cout << "Help:\n"
+            << "\tfilter <image> [sections [segma_r [segma_s] ] ]\n\n";
+
+    if (varg.size()>2) N = atoi(varg[2]);
+    if (varg.size()>3) segma_r = float(atof(varg[3]));
+    if (varg.size()>4) segma_s = float(atof(varg[4]));
+
+
+    Mat img = imread(imgfile, CV_LOAD_IMAGE_COLOR);
+    if (img.empty())
+        return -1;
+
+    //////////////////////////////////////////////////////////////////////////
+    // Kuwahara�˲�
+    CFilter filter;
+    filter.init(segma_r, segma_s, N);
+    Mat res;
+    filter.apply(img, res);
+    imwrite("filter_result.png", res);
+
+
+    //////////////////////////////////////////////////////////////////////////
+    // ��˹�˲������ڶԱ�
+    Mat img_g;
+    GaussianBlur(img, img_g, Size(), 2.0);
+
+
+    //////////////////////////////////////////////////////////////////////////
+    // ��ʾ���
+    namedWindow("Kuwahera", CV_WINDOW_AUTOSIZE);
+    imshow("Kuwahera", res);
+
+    namedWindow("original", CV_WINDOW_AUTOSIZE);
+    imshow("original", img);
+
+    namedWindow("Gaussian Blur", CV_WINDOW_AUTOSIZE);
+    imshow("Gaussian Blur", img_g);
+
+    waitKey();
+
+    return 0;
+}