diff --git a/modules/ximgproc/README.md b/modules/ximgproc/README.md index b566aa354..1afdd610b 100644 --- a/modules/ximgproc/README.md +++ b/modules/ximgproc/README.md @@ -9,3 +9,4 @@ Extended Image Processing 6. Superpixels 7. Graph segmentation 8. Selective search from segmentation +10. Paillou Filter diff --git a/modules/ximgproc/doc/ximgproc.bib b/modules/ximgproc/doc/ximgproc.bib index 3de1724c9..0fbbb374b 100644 --- a/modules/ximgproc/doc/ximgproc.bib +++ b/modules/ximgproc/doc/ximgproc.bib @@ -166,3 +166,13 @@ year={2014}, organization={IEEE} } + +@article{paillou1997detecting, + title={Detecting step edges in noisy SAR images: a new linear operator}, + author={Paillou, Philippe}, + journal={IEEE transactions on geoscience and remote sensing}, + volume={35}, + number={1}, + pages={191--196}, + year={1997} +} diff --git a/modules/ximgproc/include/opencv2/ximgproc.hpp b/modules/ximgproc/include/opencv2/ximgproc.hpp index 84015b910..b803ce3ef 100644 --- a/modules/ximgproc/include/opencv2/ximgproc.hpp +++ b/modules/ximgproc/include/opencv2/ximgproc.hpp @@ -48,6 +48,8 @@ #include "ximgproc/weighted_median_filter.hpp" #include "ximgproc/slic.hpp" #include "ximgproc/lsc.hpp" +#include "ximgproc/paillou_filter.hpp" + /** @defgroup ximgproc Extended Image Processing @{ diff --git a/modules/ximgproc/include/opencv2/ximgproc/paillou_filter.hpp b/modules/ximgproc/include/opencv2/ximgproc/paillou_filter.hpp new file mode 100644 index 000000000..e58b609d5 --- /dev/null +++ b/modules/ximgproc/include/opencv2/ximgproc/paillou_filter.hpp @@ -0,0 +1,67 @@ +/* + * By downloading, copying, installing or using the software you agree to this license. + * If you do not agree to this license, do not download, install, + * copy or use the software. + * + * + * License Agreement + * For Open Source Computer Vision Library + * (3 - clause BSD License) + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met : + * + * *Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and / or other materials provided with the distribution. + * + * * Neither the names of the copyright holders nor the names of the contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * This software is provided by the copyright holders and contributors "as is" and + * any express or implied warranties, including, but not limited to, the implied + * warranties of merchantability and fitness for a particular purpose are disclaimed. + * In no event shall copyright holders or contributors be liable for any direct, + * indirect, incidental, special, exemplary, or consequential damages + * (including, but not limited to, procurement of substitute goods or services; + * loss of use, data, or profits; or business interruption) however caused + * and on any theory of liability, whether in contract, strict liability, + * or tort(including negligence or otherwise) arising in any way out of + * the use of this software, even if advised of the possibility of such damage. + */ + +#ifndef __OPENCV_PAILLOUFILTER_HPP__ +#define __OPENCV_PAILLOUFILTER_HPP__ +#ifdef __cplusplus + +#include + +namespace cv { +namespace ximgproc { + +//! @addtogroup ximgproc_filters +//! @{ + +/** +* @brief Applies Paillou filter to an image. +* +* For more details about this implementation, please see @cite paillou1997detecting +* +* @param op Source 8-bit or 16bit image, 1-channel or 3-channel image. +* @param _dst result CV_32F image with same numeber of channel than op. +* @param omega double see paper +* @param alpha double see paper +* +* @sa GradientPaillouX, GradientPaillouY +*/ +CV_EXPORTS void GradientPaillouY(InputArray op, OutputArray _dst, double alpha, double omega); +CV_EXPORTS void GradientPaillouX(InputArray op, OutputArray _dst, double alpha, double omega); + +} +} +#endif +#endif diff --git a/modules/ximgproc/samples/paillou_demo.cpp b/modules/ximgproc/samples/paillou_demo.cpp new file mode 100644 index 000000000..7ac98b682 --- /dev/null +++ b/modules/ximgproc/samples/paillou_demo.cpp @@ -0,0 +1,107 @@ +/* + * By downloading, copying, installing or using the software you agree to this license. + * If you do not agree to this license, do not download, install, + * copy or use the software. + * + * + * License Agreement + * For Open Source Computer Vision Library + * (3 - clause BSD License) + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met : + * + * *Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and / or other materials provided with the distribution. + * + * * Neither the names of the copyright holders nor the names of the contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * This software is provided by the copyright holders and contributors "as is" and + * any express or implied warranties, including, but not limited to, the implied + * warranties of merchantability and fitness for a particular purpose are disclaimed. + * In no event shall copyright holders or contributors be liable for any direct, + * indirect, incidental, special, exemplary, or consequential damages + * (including, but not limited to, procurement of substitute goods or services; + * loss of use, data, or profits; or business interruption) however caused + * and on any theory of liability, whether in contract, strict liability, + * or tort(including negligence or otherwise) arising in any way out of + * the use of this software, even if advised of the possibility of such damage. + */ + +#include +#include +#include +#include +#include "opencv2/ximgproc/paillou_filter.hpp" + +using namespace cv; +using namespace cv::ximgproc; + +#include +using namespace std; + +int aa = 100, ww = 10; +Mat dx, dy; +UMat img; +const char* window_name = "Gradient Modulus"; + +static void DisplayImage(Mat x,string s) +{ + vector sx; + split(x, sx); + vector minVal(3), maxVal(3); + for (int i = 0; i < static_cast(sx.size()); i++) + { + minMaxLoc(sx[i], &minVal[i], &maxVal[i]); + } + maxVal[0] = *max_element(maxVal.begin(), maxVal.end()); + minVal[0] = *min_element(minVal.begin(), minVal.end()); + Mat uc; + x.convertTo(uc, CV_8U,255/(maxVal[0]-minVal[0]),-255*minVal[0]/(maxVal[0]-minVal[0])); + imshow(s, uc); +} + + +/** + * @function paillouFilter + * @brief Trackbar callback + */ +static void PaillouFilter(int, void*) +{ + Mat dst; + double a=aa/100.0,w=ww/100.0; + Mat rx,ry; + GradientPaillouX(img,rx,a,w); + GradientPaillouY(img,ry,a,w); + DisplayImage(rx, "Gx"); + DisplayImage(ry, "Gy"); + add(rx.mul(rx),ry.mul(ry),dst); + sqrt(dst,dst); + DisplayImage(dst, window_name ); +} + + +int main(int argc, char* argv[]) +{ + if (argc==2) + imread(argv[1]).copyTo(img); + if (img.empty()) + { + cout << "File not found or empty image\n"; + } + imshow("Original",img); + namedWindow( window_name, WINDOW_AUTOSIZE ); + + /// Create a Trackbar for user to enter threshold + createTrackbar( "a:",window_name, &aa, 400, PaillouFilter ); + createTrackbar( "w:", window_name, &ww, 400, PaillouFilter ); + PaillouFilter(0,NULL); + waitKey(); + return 0; +} \ No newline at end of file diff --git a/modules/ximgproc/src/paillou_filter.cpp b/modules/ximgproc/src/paillou_filter.cpp new file mode 100644 index 000000000..a82158051 --- /dev/null +++ b/modules/ximgproc/src/paillou_filter.cpp @@ -0,0 +1,486 @@ +#include "precomp.hpp" +#include "opencv2/highgui.hpp" +#include +#include +#include + +namespace cv { +namespace ximgproc { + +/* +If you use this code please cite this @cite paillou1997detecting +Detecting step edges in noisy SAR images: a new linear operator IEEE Transactions on Geoscience and Remote Sensing (Volume:35 , Issue: 1 ) 1997 +*/ + + +class ParallelGradientPaillouYCols: public ParallelLoopBody +{ +private: + Mat &img; + Mat &dst; + double a; + double w; + bool verbose; +public: + ParallelGradientPaillouYCols(Mat& imgSrc, Mat &d,double aa,double ww): + img(imgSrc), + dst(d), + a(aa), + w(ww), + verbose(false) + {} + void Verbose(bool b){verbose=b;} + virtual void operator()(const Range& range) const + { + CV_Assert(img.depth()==CV_8UC1 || img.depth()==CV_16SC1 || img.depth()==CV_16UC1); + CV_Assert(dst.depth()==CV_32FC1); + if (verbose) + std::cout << getThreadNum()<<"# :Start from row " << range.start << " to " << range.end-1<<" ("<img.cols)?img.rows:img.cols; + Mat matYp(1,tailleSequence,CV_64FC1), matYm(1,tailleSequence,CV_64FC1); + double *yp=matYp.ptr(0), *ym=matYm.ptr(0); + int rows=img.rows,cols=img.cols; + + // Equation 12 p193 + double b1=-2*exp(-a)*cosh(w); + double a1=2*exp(-a)*cosh(w)-exp(-2*a)-1; + double b2=exp(-2*a); + + switch(img.depth()){ + case CV_8U : + for (int j=range.start;j(0)+j; + double border=*c1; + yp[0] = *c1 ; + c1+=cols; + yp[1] = *c1 - b1*yp[0]-b2*border; + c1+=cols; + for (int i=2;i=0;i--,c1-=cols) + ym[i]=*c1-b1*ym[i+1]-b2*ym[i+2]; + // Equation 25 p193 + for (int i=0;i(0) + j; + f2 = dst.ptr(0) + j; + double border = *c1; + yp[0] = *c1; + c1 += cols; + yp[1] = *c1 - b1*yp[0] - b2*border; + c1 += cols; + for (int i = 2; i(rows - 1) + j; + border = *c1; + ym[rows - 1] = *c1; + c1 -= cols; + ym[rows - 2] = *c1 - b1*ym[rows - 1]; + c1 -= cols; + for (int i = rows - 3; i >= 0; i--, c1 -= cols) + ym[i] = *c1 - b1*ym[i + 1] - b2*ym[i + 2]; + // Equation 25 p193 + for (int i = 0; i(0) + j; + f2 = dst.ptr(0) + j; + double border = *c1; + yp[0] = *c1; + c1 += cols; + yp[1] = *c1 - b1*yp[0] - b2*border; + c1 += cols; + for (int i = 2; i(rows - 1) + j; + border = *c1; + ym[rows - 1] = *c1; + c1 -= cols; + ym[rows - 2] = *c1 - b1*ym[rows - 1]; + c1 -= cols; + for (int i = rows - 3; i >= 0; i--, c1 -= cols) + ym[i] = *c1 - b1*ym[i + 1] - b2*ym[i + 2]; + // Equation 25 p193 + for (int i = 0; iimg.cols)?img.rows:img.cols; + Mat matIym(1,tailleSequence,CV_64FC1), matIyp(1,tailleSequence,CV_64FC1); + double *iym=matIym.ptr(0), *iyp=matIyp.ptr(0); + int cols=img.cols; + + // Equation 13 p193 + double d=(1-2*exp(-a)*cosh(w)+exp(-2*a))/(2*a*exp(-a)*sinh(w)+w*(1-exp(-2*a))); + double c1=a*d; + double c2=w*d; + // Equation 12 p193 + double b1=-2*exp(-a)*cosh(w); + double b2=exp(-2*a); + // Equation 14 p193 + double a0p=c2; + double a1p=(c1*sinh(w)-c2*cosh(w))*exp(-a); + double a1m=a1p-c2*b1; + double a2m=-c2*b2; + + for (int i=range.start;i(i); + int j=0; + iyp[0] = a0p*iy0[0] ; + iyp[1] = a0p*iy0[1] + a1p*iy0[0] - b1*iyp[0]; + iy0 += 2; + for (j=2;j(i)+cols-1; + iym[cols-1] = 0; + iy0--; + iym[cols-2] = a1m*iy0[1] - b1*iym[cols-1]; + iy0--; + for (j=cols-3;j>=0;j--,iy0--) + iym[j] = a1m*iy0[1] + a2m*iy0[2] - b1*iym[j+1] - b2*iym[j+2]; + iy = dst.ptr(i); + for (j=0;jimg.cols) ? img.rows : img.cols; + Mat matIym(1,tailleSequence,CV_64FC1), matIyp(1,tailleSequence,CV_64FC1); + double *iym=matIym.ptr(0), *iyp=matIyp.ptr(0); + int rows = img.rows,cols=img.cols; + + // Equation 13 p193 + double d = (1 - 2 * exp(-a)*cosh(w) + exp(-2 * a)) / (2 * a*exp(-a)*sinh(w) + w*(1 - exp(-2 * a))); + double c1 = a*d; + double c2 = w*d; + // Equation 12 p193 + double b1 = -2 * exp(-a)*cosh(w); + double b2 = exp(-2 * a); + // Equation 14 p193 + double a0p = c2; + double a1p = (c1*sinh(w) - c2*cosh(w))*exp(-a); + double a1m = a1p - c2*b1; + double a2m = -c2*b2; + + for (int j = range.start; j(0)+j; + iyp[0] = a0p*iy0[0]; + iy0 +=cols; + iyp[1] = a0p*iy0[0] + a1p*iy0[-cols] - b1*iyp[0]; + iy0 +=cols; + for (int i = 2; i(rows-1) + j; + iym[rows - 1] = 0; + iy0 -=cols; + iym[rows - 2] = a1m*iy0[cols] - b1*iym[rows-1]; + iy0-=cols; + for (int i = rows - 3; i >= 0; i--, iy0-=cols) + iym[i] = a1m*iy0[cols] + a2m*iy0[2*cols] - b1*iym[i + 1] - b2*iym[i + 2]; + iy = dst.ptr(0)+j; + for (int i = 0; iimg.cols) ? img.rows : img.cols; + Mat matYp(1,tailleSequence,CV_64FC1), matYm(1,tailleSequence,CV_64FC1); + double *yp=matYp.ptr(0), *ym=matYm.ptr(0); + int cols = img.cols; + + // Equation 12 p193 + double b1 = -2 * exp(-a)*cosh(w); + double a1 = 2 * exp(-a)*cosh(w) - exp(-2 * a) - 1; + double b2 = exp(-2 * a); + + switch(img.depth()){ + case CV_8U : + for (int i = range.start; i= 0; j--, c1--) + ym[j] = *c1 - b1*ym[j + 1] - b2*ym[j + 2]; + // Equation 25 p193 + f2 = im1.ptr(i); + for (int j = 0; j(i); + double border = *c1; + yp[0] = *c1; + c1++; + yp[1] = *c1 - b1*yp[0] - b2*border; + c1++; + for (int j = 2; j(i)+cols-1; + border = *c1; + ym[cols - 1] = *c1; + c1--; + ym[cols - 2] = *c1 - b1*ym[cols - 1]; + c1--; + for (int j = cols - 3; j >= 0; j--, c1--) + ym[j] = *c1 - b1*ym[j + 1] - b2*ym[j + 2]; + // Equation 25 p193 + f2 = im1.ptr(i); + for (int j = 0; j(i); + f2 = im1.ptr(i); + double border = *c1; + yp[0] = *c1; + c1++; + yp[1] = *c1 - b1*yp[0] - b2*border; + c1++; + for (int j = 2; j(i) + cols - 1; + border = *c1; + ym[cols - 1] = *c1; + c1--; + ym[cols - 2] = *c1 - b1*ym[cols - 1]; + c1--; + for (int j = cols - 3; j >= 0; j--, c1--) + ym[j] = *c1 - b1*ym[j + 1] - b2*ym[j + 2]; + // Equation 25 p193 + for (int j = 0; j(i); + f2 = im1.ptr(i); + double border = *c1; + yp[0] = *c1; + c1++; + yp[1] = *c1 - b1*yp[0] - b2*border; + c1++; + for (int j = 2; j(i) + cols - 1; + border = *c1; + ym[cols - 1] = *c1; + c1--; + ym[cols - 2] = *c1 - b1*ym[cols - 1]; + c1--; + for (int j = cols - 3; j >= 0; j--, c1--) + ym[j] = *c1 - b1*ym[j + 1] - b2*ym[j + 2]; + // Equation 25 p193 + for (int j = 0; j planSrc; + split(opSrc,planSrc); + std::vector planTmp; + split(tmp,planTmp); + std::vector planDst; + split(_dst,planDst); + for (int i = 0; i < static_cast(planSrc.size()); i++) + { + if (planSrc[i].isContinuous() && planTmp[i].isContinuous() && planDst[i].isContinuous()) + { + ParallelGradientPaillouYCols x(planSrc[i],planTmp[i],alpha,omega); + parallel_for_(Range(0,opSrc.cols), x,getNumThreads()); + ParallelGradientPaillouYRows xr(planTmp[i],planDst[i],alpha,omega); + parallel_for_(Range(0,opSrc.rows), xr,getNumThreads()); + + } + else + std::cout << "PB"; + } + merge(planDst,_dst); +} + +void GradientPaillouX(InputArray _op, OutputArray _dst, double alpha, double omega) +{ + Mat tmp(_op.size(),CV_32FC(_op.channels())); + _dst.create( _op.size(),CV_32FC(tmp.channels()) ); + Mat opSrc = _op.getMat(); + std::vector planSrc; + split(opSrc,planSrc); + std::vector planTmp; + split(tmp,planTmp); + std::vector planDst; + split(_dst,planDst); + for (int i = 0; i < static_cast(planSrc.size()); i++) + { + if (planSrc[i].isContinuous() && planTmp[i].isContinuous() && planDst[i].isContinuous()) + { + ParallelGradientPaillouXRows x(planSrc[i],planTmp[i],alpha,omega); + parallel_for_(Range(0,opSrc.rows), x,getNumThreads()); + ParallelGradientPaillouXCols xr(planTmp[i],planDst[i],alpha,omega); + parallel_for_(Range(0,opSrc.cols), xr,getNumThreads()); + } + else + std::cout << "PB"; + } + merge(planDst,_dst); +} +} +}