From cd6f7d1291898cf9a1185e5974a6f20b4127b3d0 Mon Sep 17 00:00:00 2001 From: LaurentBerger Date: Thu, 8 Nov 2018 22:07:42 +0100 Subject: [PATCH] Rewrite deriche filter- add test - add python wrapper --- .../opencv2/ximgproc/deriche_filter.hpp | 20 +- modules/ximgproc/samples/dericheSample.py | 60 ++++++ modules/ximgproc/src/deriche_filter.cpp | 194 ++++++++++-------- modules/ximgproc/test/test_deriche_filter.cpp | 41 ++++ 4 files changed, 214 insertions(+), 101 deletions(-) create mode 100644 modules/ximgproc/samples/dericheSample.py create mode 100644 modules/ximgproc/test/test_deriche_filter.cpp diff --git a/modules/ximgproc/include/opencv2/ximgproc/deriche_filter.hpp b/modules/ximgproc/include/opencv2/ximgproc/deriche_filter.hpp index 2371feb6d..26d3b6759 100644 --- a/modules/ximgproc/include/opencv2/ximgproc/deriche_filter.hpp +++ b/modules/ximgproc/include/opencv2/ximgproc/deriche_filter.hpp @@ -51,25 +51,25 @@ namespace ximgproc { * * For more details about this implementation, please see http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.476.5736&rep=rep1&type=pdf * -* @param _op Source 8-bit or 16bit image, 1-channel or 3-channel image. -* @param _dst result CV_32FC image with same number of channel than _op. -* @param alphaDerive double see paper -* @param alphaMean double see paper +* @param op Source 8-bit or 16bit image, 1-channel or 3-channel image. +* @param dst result CV_32FC image with same number of channel than _op. +* @param alpha double see paper +* @param omega double see paper * */ -CV_EXPORTS void GradientDericheY(InputArray _op, OutputArray _dst, double alphaDerive,double alphaMean); +CV_EXPORTS_W void GradientDericheY(InputArray op, OutputArray dst, double alpha,double omega); /** * @brief Applies X Deriche filter to an image. * * For more details about this implementation, please see http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.476.5736&rep=rep1&type=pdf * -* @param _op Source 8-bit or 16bit image, 1-channel or 3-channel image. -* @param _dst result CV_32FC image with same number of channel than _op. -* @param alphaDerive double see paper -* @param alphaMean double see paper +* @param op Source 8-bit or 16bit image, 1-channel or 3-channel image. +* @param dst result CV_32FC image with same number of channel than _op. +* @param alpha double see paper +* @param omega double see paper * */ -CV_EXPORTS void GradientDericheX(InputArray _op, OutputArray _dst, double alphaDerive,double alphaMean); +CV_EXPORTS_W void GradientDericheX(InputArray op, OutputArray dst, double alpha,double omega); } } diff --git a/modules/ximgproc/samples/dericheSample.py b/modules/ximgproc/samples/dericheSample.py new file mode 100644 index 000000000..6468b07c8 --- /dev/null +++ b/modules/ximgproc/samples/dericheSample.py @@ -0,0 +1,60 @@ +import sys +import numpy as np +import cv2 as cv + +def AddSlider(sliderName,windowName,minSlider,maxSlider,valDefault, update=[]): + if update is None: + cv.createTrackbar(sliderName, windowName, valDefault,maxSlider-minSlider+1) + else: + cv.createTrackbar(sliderName, windowName, valDefault,maxSlider-minSlider+1, update) + cv.setTrackbarMin(sliderName, windowName, minSlider) + cv.setTrackbarMax(sliderName, windowName, maxSlider) + cv.setTrackbarPos(sliderName, windowName, valDefault) +class Filtrage: + def __init__(self): + self.s =0 + self.alpha = 100 + self.omega = 100 + self.updateFiltre=True + self.img=[] + self.dximg=[] + self.dyimg=[] + self.module=[] + def DericheFilter(self): + self.dximg = cv.ximgproc.GradientDericheX( self.img, self.alpha/100., self.omega/1000. ) + self.dyimg = cv.ximgproc.GradientDericheY( self.img, self.alpha/100., self.omega/1000. ) + dx2=self.dximg*self.dximg + dy2=self.dyimg*self.dyimg + self.module = np.sqrt(dx2+dy2) + cv.normalize(src=self.module,dst=self.module,norm_type=cv.NORM_MINMAX) + def SlideBarDeriche(self): + cv.destroyWindow(self.filename) + cv.namedWindow(self.filename) + AddSlider("alpha",self.filename,1,400,self.alpha,self.UpdateAlpha) + AddSlider("omega",self.filename,1,1000,self.omega,self.UpdateOmega) + + def UpdateOmega(self,x ): + self.updateFiltre=True + self.omega=x + def UpdateAlpha(self,x ): + self.updateFiltre=True + self.alpha=x + def run(self,argv): + # Load the source image + self.filename = argv[0] if len(argv) > 0 else "../doc/pics/corridor_fld.jpg" + self.img=cv.imread(self.filename,cv.IMREAD_GRAYSCALE) + if self.img is None: + print ('cannot read file') + return + self.SlideBarDeriche() + while True: + cv.imshow(self.filename,self.img) + if self.updateFiltre: + self.DericheFilter() + cv.imshow("module",self.module) + self.updateFiltre =False + code = cv.waitKey(10) + if code==27: + break +if __name__ == '__main__': + Filtrage().run(sys.argv[1:]) diff --git a/modules/ximgproc/src/deriche_filter.cpp b/modules/ximgproc/src/deriche_filter.cpp index 3a20b2cf1..635004972 100644 --- a/modules/ximgproc/src/deriche_filter.cpp +++ b/modules/ximgproc/src/deriche_filter.cpp @@ -46,23 +46,23 @@ Using Canny's criteria to derive a recursively implemented optimal edge detector namespace cv { namespace ximgproc { template static void -VerticalIIRFilter(Mat &img,Mat &dst,const Range &r,double alphaDerive) +VerticalIIRFilter(Mat &img,Mat &dst,const Range &r,double alpha,double omega) { float *f2; int tailleSequence = (img.rows>img.cols) ? img.rows : img.cols; Mat matG1(1, tailleSequence, CV_64FC1), matG2(1, tailleSequence, CV_64FC1); double *g1 = matG1.ptr(0), *g2 = (double*)matG2.ptr(0); - double kp = pow(1 - exp(-alphaDerive), 2.0) / exp(-alphaDerive); - double a1, a2, a3, a4; + double a2, a3; double b1, b2; int rows = img.rows, cols = img.cols; + double c = (1 - 2 * exp(-alpha)*cos(omega) + exp(-2 * alpha)) / (exp(-alpha)*sin(omega)); + + double a = -c * exp(-alpha)*sin(omega); + a2 = 1;// kp*exp(-alpha); + a3 = 1;//-kp*exp(-alpha); + b1 = -2 * exp(-alpha)*cos(omega); + b2 = exp(-2 * alpha); - kp = pow(1 - exp(-alphaDerive), 2.0) / exp(-alphaDerive); - a1 = 0; - a2 = kp*exp(-alphaDerive), a3 = -kp*exp(-alphaDerive); - a4 = 0; - b1 = 2 * exp(-alphaDerive); - b2 = -exp(-2 * alphaDerive); for (int j = r.start; j(0); c1 += (rows - 1)*cols + j; i = rows - 1; - g2[i] = (a3 + a4)* *c1; + g2[i] = a3 * *c1; i--; c1 -= cols; - g2[i] = a3* c1[cols] + a4 * c1[cols] + (b1)*g2[i + 1]; + g2[i] = a3* c1[cols] + (b1)*g2[i + 1]; i--; c1 -= cols; for (i = rows - 3; i >= 0; i--, c1 -= cols) - g2[i] = a3*c1[cols] + a4* c1[2 * cols] + - b1*g2[i + 1] + b2*g2[i + 2]; + g2[i] = a3*c1[cols] - + b1*g2[i + 1] - b2*g2[i + 2]; for (i = 0; i(a*(g1[i] - g2[i])); } } template static void -HorizontalIIRFilter(Mat &img, Mat &dst, const Range &r, double alphaDerive) +HorizontalIIRFilter(Mat &img, Mat &dst, const Range &r, double alpha, double omega) { float *f1; int rows = img.rows, cols = img.cols; int tailleSequence = (rows>cols) ? rows : cols; Mat matG1(1, tailleSequence, CV_64FC1), matG2(1, tailleSequence, CV_64FC1); double *g1 = (double*)matG1.ptr(0), *g2 = (double*)matG2.ptr(0); - double kp;; - double a1, a2, a3, a4; + double a,a2, a3; double b1, b2; + double c = (1 - 2 * exp(-alpha)*cos(omega) + exp(-2 * alpha)) / (exp(-alpha)*sin(omega)); - kp = pow(1 - exp(-alphaDerive), 2.0) / exp(-alphaDerive); - a1 = 0; - a2 = kp*exp(-alphaDerive); - a3 = -kp*exp(-alphaDerive); - a4 = 0; - b1 = 2 * exp(-alphaDerive); - b2 = -exp(-2 * alphaDerive); + a = -c*exp(-alpha)*sin(omega); + a2 = 1;// kp*exp(-alpha); + a3 = 1;//-kp*exp(-alpha); + b1 = -2 * exp(-alpha)*cos(omega); + b2 = exp(-2 * alpha); for (int i = r.start; i(i); T *c1 = img.ptr(i); int j = 0; - g1[j] = (a1 + a2)* *c1; + g1[j] = a2* *c1; j++; c1++; - g1[j] = a1 * c1[0] + a2*c1[j - 1] + (b1)* g1[j - 1]; + g1[j] = a2*c1[j - 1] - (b1)* g1[j - 1]; j++; c1++; for (j = 2; j(0); c1 += i*cols + cols - 1; j = cols - 1; - g2[j] = (a3 + a4)* *c1; + g2[j] = a3* *c1; j--; c1--; - g2[j] = (a3 + a4) * c1[1] + b1 * g2[j + 1]; + g2[j] = a3 * c1[1] - b1 * g2[j + 1]; j--; c1--; for (j = cols - 3; j >= 0; j--, c1--) - g2[j] = a3*c1[1] + a4*c1[2] + b1*g2[j + 1] + b2*g2[j + 2]; + g2[j] = a3*c1[1] - b1*g2[j + 1] - b2*g2[j + 2]; for (j = 0; j(a*(g1[j] - g2[j])); } } @@ -151,15 +149,17 @@ class ParallelGradientDericheYCols : public ParallelLoopBody private: Mat &img; Mat &dst; - double alphaDerive; + double alpha; + double omega; bool verbose; public: - ParallelGradientDericheYCols(Mat &imgSrc, Mat &d, double ald) : + ParallelGradientDericheYCols(Mat &imgSrc, Mat &d, double ald,double o) : img(imgSrc), dst(d), - alphaDerive(ald), + alpha(ald), + omega(o), verbose(false) { int type = img.depth(); @@ -175,19 +175,19 @@ public: switch (img.depth()) { case CV_8U: - VerticalIIRFilter(img,dst,range, alphaDerive); + VerticalIIRFilter(img,dst,range, alpha,omega); break; case CV_8S: - VerticalIIRFilter(img, dst, range, alphaDerive); + VerticalIIRFilter(img, dst, range, alpha, omega); break; case CV_16U: - VerticalIIRFilter(img, dst, range, alphaDerive); + VerticalIIRFilter(img, dst, range, alpha, omega); break; case CV_16S: - VerticalIIRFilter(img, dst, range, alphaDerive); + VerticalIIRFilter(img, dst, range, alpha, omega); break; case CV_32F: - VerticalIIRFilter(img, dst, range, alphaDerive); + VerticalIIRFilter(img, dst, range, alpha, omega); break; default: return; @@ -204,14 +204,16 @@ class ParallelGradientDericheYRows : public ParallelLoopBody private: Mat &img; Mat &dst; - double alphaMoyenne; + double alpha; + double omega; bool verbose; public: - ParallelGradientDericheYRows(Mat& imgSrc, Mat &d, double alm) : + ParallelGradientDericheYRows(Mat& imgSrc, Mat &d, double ald,double o) : img(imgSrc), dst(d), - alphaMoyenne(alm), + alpha(ald), + omega(o), verbose(false) { int type = img.depth(); @@ -228,42 +230,44 @@ public: int tailleSequence = (img.rows>img.cols) ? img.rows : img.cols; Mat matG1(1,tailleSequence,CV_64FC1), matG2(1,tailleSequence,CV_64FC1); double *g1 = matG1.ptr(0), *g2 = matG2.ptr(0); - double k, a5, a6, a7, a8; - double b3, b4; int cols = img.cols; - k = pow(1 - exp(-alphaMoyenne), 2.0) / (1 + 2 * alphaMoyenne*exp(-alphaMoyenne) - exp(-2 * alphaMoyenne)); - a5 = k; - a6 = k*exp(-alphaMoyenne)*(alphaMoyenne - 1); - a7 = k*exp(-alphaMoyenne)*(alphaMoyenne + 1); - a8 = -k*exp(-2 * alphaMoyenne); - b3 = 2 * exp(-alphaMoyenne); - b4 = -exp(-2 * alphaMoyenne); + + double a2po2 = (alpha*alpha + omega * omega); + double k = (1 - 2 * exp(-alpha)*cos(omega) + exp(-2 * alpha))*a2po2; + k = k / (2 * alpha*exp(-alpha)*sin(omega) + omega - omega * exp(-2 * alpha)); + double c1 = k * alpha / a2po2; + double c2 = k * omega / a2po2; + double a0 = c2; + double a1 = (-c2 * cos(omega) + c1 * sin(omega))*exp(-alpha); + double b1 = -2 * exp(-alpha)*cos(omega); + double b2 = exp(-2 * alpha); + double a2 = a1 - c2 * b1, a3 = -c2 * b2; for (int i = range.start; i(i); f1 = img.ptr(i); int j = 0; - g1[j] = (a5 + a6)* *f1; + g1[j] = (a0 + a1)* *f1; j++; f1++; - g1[j] = a5 * f1[0] + a6*f1[j - 1] + (b3)* g1[j - 1]; + g1[j] = a0 * f1[0] + a1*f1[j - 1] - b1* g1[j - 1]; j++; f1++; for (j = 2; j= 0; j--, f1--) - g2[j] = a7*f1[1] + a8*f1[2] + b3*g2[j + 1] + b4*g2[j + 2]; + g2[j] = a2*f1[1] + a3*f1[2] - b1*g2[j + 1] - b2*g2[j + 2]; for (j = 0; jcols) ? rows : cols; Mat matG1(1,tailleSequence,CV_64FC1), matG2(1,tailleSequence,CV_64FC1); double *g1 = (double*)matG1.ptr(0), *g2 = (double*)matG2.ptr(0); - double k, a5, a6, a7, a8 = 0; - double b3, b4; + double a2po2 = (alpha*alpha + omega * omega); + double k = (1 - 2 * exp(-alpha)*cos(omega) + exp(-2 * alpha))*a2po2; + k = k / (2 * alpha*exp(-alpha)*sin(omega) + omega - omega * exp(-2 * alpha)); + double c1 = k * alpha / a2po2; + double c2 = k * omega / a2po2; + double a0 = c2; + double a1 = (-c2 * cos(omega) + c1 * sin(omega))*exp(-alpha); + double b1 = -2 * exp(-alpha)*cos(omega); + double b2=exp(-2*alpha); + double a2=a1-c2*b1, a3=-c2*b2; - k = pow(1 - exp(-alphaMoyenne), 2.0) / (1 + 2 * alphaMoyenne*exp(-alphaMoyenne) - exp(-2 * alphaMoyenne)); - a5 = k, a6 = k*exp(-alphaMoyenne)*(alphaMoyenne - 1); - a7 = k*exp(-alphaMoyenne)*(alphaMoyenne + 1), a8 = -k*exp(-2 * alphaMoyenne); - b3 = 2 * exp(-alphaMoyenne); - b4 = -exp(-2 * alphaMoyenne); for (int j = range.start; j(0); f1 += j; int i = 0; - g1[i] = (a5 + a6)* *f1; + g1[i] = (a0 + a1)* *f1; i++; f1 += cols; - g1[i] = a5 * *f1 + a6 * f1[-cols] + (b3)* g1[i - 1]; + g1[i] = a0 * *f1 + a1 * f1[-cols] - (b1)* g1[i - 1]; i++; f1 += cols; for (i = 2; i(0); f1 += (rows - 1)*cols + j; i = rows - 1; - g2[i] = (a7 + a8)* *f1; + g2[i] = (a2 + a3)* *f1; i--; f1 -= cols; - g2[i] = (a7 + a8)* f1[cols] + (b3)*g2[i + 1]; + g2[i] = (a2 + a3)* f1[cols] - b2*g2[i + 1]; i--; f1 -= cols; for (i = rows - 3; i >= 0; i--, f1 -= cols) - g2[i] = a7*f1[cols] + a8* f1[2 * cols] + - b3*g2[i + 1] + b4*g2[i + 2]; + g2[i] = a2*f1[cols] + a3* f1[2 * cols] - + b1*g2[i + 1] - b2*g2[i + 2]; for (i = 0; i(i)) + (j*img.channels()); @@ -358,14 +367,16 @@ class ParallelGradientDericheXRows : public ParallelLoopBody private: Mat &img; Mat &dst; - double alphaDerive; + double alpha; + double omega; bool verbose; public: - ParallelGradientDericheXRows(Mat& imgSrc, Mat &d, double ald) : + ParallelGradientDericheXRows(Mat& imgSrc, Mat &d, double ald, double o) : img(imgSrc), dst(d), - alphaDerive(ald), + alpha(ald), + omega(o), verbose(false) { int type = img.depth(); @@ -381,19 +392,19 @@ public: switch (img.depth()) { case CV_8U: - HorizontalIIRFilter(img,dst,range,alphaDerive); + HorizontalIIRFilter(img,dst,range, alpha,omega); break; case CV_8S: - HorizontalIIRFilter(img, dst, range, alphaDerive); + HorizontalIIRFilter(img, dst, range, alpha, omega); break; case CV_16U: - HorizontalIIRFilter(img, dst, range, alphaDerive); + HorizontalIIRFilter(img, dst, range, alpha, omega); break; case CV_16S: - HorizontalIIRFilter(img, dst, range, alphaDerive); + HorizontalIIRFilter(img, dst, range, alpha, omega); break; case CV_32F: - HorizontalIIRFilter(img, dst, range, alphaDerive); + HorizontalIIRFilter(img, dst, range, alpha, omega); break; default: return; @@ -404,10 +415,11 @@ public: }; }; -void GradientDericheY(InputArray _op, OutputArray _dst,double alphaDerive, double alphaMean) +void GradientDericheY(InputArray _op, OutputArray _dst,double alphaDerive, double omega) { std::vector planSrc; split(_op, planSrc); + std::vector planTmp; std::vector planDst; for (size_t i = 0; i < planSrc.size(); i++) @@ -415,15 +427,15 @@ void GradientDericheY(InputArray _op, OutputArray _dst,double alphaDerive, doubl planTmp.push_back(Mat(_op.size(), CV_32FC1)); planDst.push_back(Mat(_op.size(), CV_32FC1)); CV_Assert(planSrc[i].isContinuous() && planTmp[i].isContinuous() && planDst[i].isContinuous()); - ParallelGradientDericheYCols x(planSrc[i], planTmp[i], alphaDerive); + ParallelGradientDericheYCols x(planSrc[i], planTmp[i], alphaDerive,omega); parallel_for_(Range(0, planSrc[i].cols), x, getNumThreads()); - ParallelGradientDericheYRows xr(planTmp[i], planDst[i], alphaMean); + ParallelGradientDericheYRows xr(planTmp[i], planDst[i], alphaDerive, omega); parallel_for_(Range(0, planTmp[i].rows), xr, getNumThreads()); } merge(planDst, _dst); } -void GradientDericheX(InputArray _op, OutputArray _dst, double alphaDerive, double alphaMean) +void GradientDericheX(InputArray _op, OutputArray _dst, double alpha, double omega) { std::vector planSrc; split(_op, planSrc); @@ -435,9 +447,9 @@ void GradientDericheX(InputArray _op, OutputArray _dst, double alphaDerive, doub planDst.push_back(Mat(_op.size(), CV_32FC1)); CV_Assert(planSrc[i].isContinuous() && planTmp[i].isContinuous() && planDst[i].isContinuous()); - ParallelGradientDericheXRows x(planSrc[i], planTmp[i], alphaDerive); + ParallelGradientDericheXRows x(planSrc[i], planTmp[i], alpha, omega); parallel_for_(Range(0, planSrc[i].rows), x, getNumThreads()); - ParallelGradientDericheXCols xr(planTmp[i], planDst[i], alphaMean); + ParallelGradientDericheXCols xr(planTmp[i], planDst[i], alpha, omega); parallel_for_(Range(0, planTmp[i].cols), xr, getNumThreads()); } merge(planDst, _dst); diff --git a/modules/ximgproc/test/test_deriche_filter.cpp b/modules/ximgproc/test/test_deriche_filter.cpp new file mode 100644 index 000000000..65aabae84 --- /dev/null +++ b/modules/ximgproc/test/test_deriche_filter.cpp @@ -0,0 +1,41 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. +#include "test_precomp.hpp" + +namespace opencv_test { namespace { + +TEST(ximgproc_DericheFilter, regression) +{ + Mat img = Mat::zeros(64, 64, CV_8UC3); + Mat res = Mat::zeros(64, 64, CV_32FC3); + img.at(31, 31) = Vec3b(1, 2, 4); + double a = 0.5; + double w = 0.0005; + Mat dst; + + ximgproc::GradientDericheX(img, dst, a, w); + double c = pow(1 - exp(-a), 2.0) * exp(a); + double k = pow(a*(1 - exp(-a)), 2.0) / (1 + 2 * a*exp(-a) - exp(-2 * a)); + for (int i = 0; i < img.rows; i++) + { + double n = -31 + i; + for (int j = 0; j < img.cols; j++) + { + double m = -31 + j; + double x = -c * exp(-a * fabs(m))*sin(w*m); + x = x * (k*(a*sin(w*fabs(n)) + w * cos(w*fabs(n)))*exp(-a * fabs(n))) / (a*a + w * w); + x = x / (w*w); + float xx=static_cast(x); + res.at(i, j) = Vec3f(xx, 2 * xx, 4 * xx); + } + } + EXPECT_LE(cv::norm(res, dst, NORM_INF), 1e-5); + + Mat dst2; + ximgproc::GradientDericheY(img, dst2, a, w); + cv::transpose(dst2, dst2); + EXPECT_LE(cv::norm(dst2, dst, NORM_INF), 1e-5); +} +} +} // namespace