Skip to content

O
opencv_contrib
Commit 7360642d authored by LaurentBerger's avatar LaurentBerger
Browse files
Options

Improve paillou filter

parent 86342522
Show whitespace changes
Inline Side-by-side
Showing with 152 additions and 233 deletions
modules/ximgproc/include/opencv2/ximgproc/paillou_filter.hpp
View file @ 7360642d
......@@ -51,8 +51,8 @@ namespace ximgproc {
*
* 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 op Source CV_8U(S) or CV_16U(S), 1-channel or 3-channels image.
* @param _dst result CV_32F image with same number of channel than op.
* @param omega double see paper
* @param alpha double see paper
*
......
modules/ximgproc/samples/paillou_demo.cpp
View file @ 7360642d
......@@ -47,8 +47,8 @@ using namespace cv::ximgproc;
using namespace std;
int aa = 100, ww = 10;
Mat dx, dy;
UMat img;
Ptr<Mat> img;
const char* window_name = "Gradient Modulus";
static void DisplayImage(Mat x,string s)
......@@ -77,8 +77,8 @@ 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);
GradientPaillouX(*img.get(),rx,a,w);
GradientPaillouY(*img.get(),ry,a,w);
DisplayImage(rx, "Gx");
DisplayImage(ry, "Gy");
add(rx.mul(rx),ry.mul(ry),dst);
......@@ -89,13 +89,15 @@ static void PaillouFilter(int, void*)
int main(int argc, char* argv[])
{
if (argc==2)
imread(argv[1]).copyTo(img);
if (img.empty())
Mat *m=new Mat;
if (argc == 2)
*m = imread(argv[1]);
if (m->empty())
{
cout << "File not found or empty image\n";
}
imshow("Original",img);
img = Ptr<Mat>(m);
imshow("Original",*img.get());
namedWindow( window_name, WINDOW_AUTOSIZE );
/// Create a Trackbar for user to enter threshold
......
modules/ximgproc/src/paillou_filter.cpp
View file @ 7360642d
/*
* 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 "precomp.hpp"
#include "opencv2/highgui.hpp"
#include <math.h>
#include <vector>
#include <iostream>
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
Equation xx pyyy mean equation xx at page yyy in previous paper
*/
namespace cv {
namespace ximgproc {
class ParallelGradientPaillouYCols: public ParallelLoopBody
template<typename T> static void
VerticalIIRFilter(Mat &img, Mat &dst, const Range &range, double a,double w)
{
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<<" ("<<range.end-range.start<<" loops)" << std::endl;
float *f2;
int tailleSequence=(img.rows>img.cols)?img.rows:img.cols;
Mat matYp(1,tailleSequence,CV_64FC1), matYm(1,tailleSequence,CV_64FC1);
double *yp=matYp.ptr<double>(0), *ym=matYm.ptr<double>(0);
int rows=img.rows,cols=img.cols;
int tailleSequence = (img.rows>img.cols) ? img.rows : img.cols;
Mat matYp(1, tailleSequence, CV_64FC1), matYm(1, tailleSequence, CV_64FC1);
double *yp = matYp.ptr<double>(0), *ym = matYm.ptr<double>(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<range.end;j++)
{
// Equation 26 p194
uchar *c1 = img.ptr(0)+j;
f2 = dst.ptr<float>(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;i++,c1+=cols)
yp[i] = *c1-b1*yp[i-1]-b2*yp[i-2];
// Equation 27 p194
c1 = img.ptr(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<rows;i++,f2+=cols)
*f2 = (float)(a1*(ym[i]-yp[i]));
}
break;
case CV_16S :
for (int j = range.start; j<range.end; j++)
double b1 = -2 * exp(-a)*cosh(w);
double a1 = 2 * exp(-a)*cosh(w) - exp(-2 * a) - 1;
double b2 = exp(-2 * a);
for (int j = range.start; j < range.end; j++)
{
// Equation 26 p194
short *c1 = img.ptr<short>(0) + j;
T *c1 = (T *)img.ptr(0) + j;
f2 = dst.ptr<float>(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; i++, c1 += cols)
for (int i = 2; i < rows; i++, c1 += cols)
yp[i] = *c1 - b1*yp[i - 1] - b2*yp[i - 2];
// Equation 27 p194
c1 = img.ptr<short>(rows - 1) + j;
c1 = (T *)img.ptr(rows - 1) + j;
border = *c1;
ym[rows - 1] = *c1;
c1 -= cols;
......@@ -98,36 +83,88 @@ public:
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<rows; i++, f2 += cols)
for (int i = 0; i < rows; i++, f2 += cols)
*f2 = (float)(a1*(ym[i] - yp[i]));
}
break;
case CV_16U :
for (int j = range.start; j<range.end; j++)
}
template<typename T> static void
HorizontalIIRFilter(Mat &img, Mat &dst, const Range &range, double a, double w)
{
int tailleSequence = (img.rows>img.cols) ? img.rows : img.cols;
Mat matYp(1, tailleSequence, CV_64FC1), matYm(1, tailleSequence, CV_64FC1);
double *yp = matYp.ptr<double>(0), *ym = matYm.ptr<double>(0);
int cols = img.cols;
float *f2;
// 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);
for (int i = range.start; i<range.end; i++)
{
// Equation 26 p194
ushort *c1 = img.ptr<ushort>(0) + j;
f2 = dst.ptr<float>(0) + j;
T *c1 = img.ptr<T>(i);
double border = *c1;
yp[0] = *c1;
c1 += cols;
c1++;
yp[1] = *c1 - b1*yp[0] - b2*border;
c1 += cols;
for (int i = 2; i<rows; i++, c1 += cols)
yp[i] = *c1 - b1*yp[i - 1] - b2*yp[i - 2];
c1++;
for (int j = 2; j<cols; j++, c1++)
yp[j] = *c1 - b1*yp[j - 1] - b2*yp[j - 2];
// Equation 27 p194
c1 = img.ptr<ushort>(rows - 1) + j;
c1 = img.ptr<T>(i) + cols - 1;
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];
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 i = 0; i<rows; i++, f2 += cols)
*f2 = (float)(a1*(ym[i] - yp[i]));
f2 = dst.ptr<float>(i);
for (int j = 0; j<cols; j++, f2++)
*f2 = (float)(a1*(ym[j] - yp[j]));
}
}
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<<" ("<<range.end-range.start<<" loops)" << std::endl;
switch(img.depth()){
case CV_8U:
VerticalIIRFilter<uchar>(img, dst, range, a, w);
break;
case CV_8S:
VerticalIIRFilter<char>(img, dst, range, a, w);
break;
case CV_16S :
VerticalIIRFilter<short>(img, dst, range, a, w);
break;
case CV_16U :
VerticalIIRFilter<short>(img, dst, range, a, w);
break;
default :
return ;
......@@ -301,125 +338,20 @@ public:
{
if (verbose)
std::cout << getThreadNum()<<"# :Start from row " << range.start << " to " << range.end-1<<" ("<<range.end-range.start<<" loops)" << std::endl;
float *f2;
int tailleSequence = (img.rows>img.cols) ? img.rows : img.cols;
Mat matYp(1,tailleSequence,CV_64FC1), matYm(1,tailleSequence,CV_64FC1);
double *yp=matYp.ptr<double>(0), *ym=matYm.ptr<double>(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<range.end; i++)
{
// Equation 26 p194
uchar *c1 = img.ptr(i);
double border = *c1;
yp[0] = *c1;
c1++;
yp[1] = *c1 - b1*yp[0] - b2*border;
c1++;
for (int j = 2; j<cols; j++, c1++)
yp[j] = *c1 - b1*yp[j - 1] - b2*yp[j - 2];
// Equation 27 p194
c1 = img.ptr(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<float>(i);
for (int j = 0; j<cols; j++, f2 ++)
*f2 = (float)(a1*(ym[j] - yp[j]));
}
HorizontalIIRFilter<uchar>(img,im1,range,a,w);
break;
case CV_8S :
for (int i = range.start; i<range.end; i++)
{
// Equation 26 p194
char *c1 = img.ptr<char>(i);
double border = *c1;
yp[0] = *c1;
c1++;
yp[1] = *c1 - b1*yp[0] - b2*border;
c1++;
for (int j = 2; j<cols; j++, c1++)
yp[j] = *c1 - b1*yp[j - 1] - b2*yp[j - 2];
// Equation 27 p194
c1 = img.ptr<char>(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<float>(i);
for (int j = 0; j<cols; j++, f2 ++)
*f2 = (float)(a1*(ym[j] - yp[j]));
}
HorizontalIIRFilter<char>(img, im1, range, a, w);
break;
case CV_16S :
for (int i = range.start; i<range.end; i++)
{
// Equation 26 p194
short *c1 = img.ptr<short>(i);
f2 = im1.ptr<float>(i);
double border = *c1;
yp[0] = *c1;
c1++;
yp[1] = *c1 - b1*yp[0] - b2*border;
c1++;
for (int j = 2; j<cols; j++, c1++)
yp[j] = *c1 - b1*yp[j - 1] - b2*yp[j - 2];
// Equation 27 p194
c1 = img.ptr<short>(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<cols; j++, f2++)
*f2 = (float)(a1*(ym[i] - yp[i]));
}
HorizontalIIRFilter<short>(img, im1, range, a, w);
break;
case CV_16U :
for (int i = range.start; i<range.end; i++)
{
// Equation 26 p194
ushort *c1 = img.ptr<ushort>(i);
f2 = im1.ptr<float>(i);
double border = *c1;
yp[0] = *c1;
c1++;
yp[1] = *c1 - b1*yp[0] - b2*border;
c1++;
for (int j = 2; j<cols; j++, c1++)
yp[j] = *c1 - b1*yp[j - 1] - b2*yp[j - 2];
// Equation 27 p194
c1 = img.ptr<ushort>(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<cols; j++, f2++)
*f2 = (float)(a1*(ym[i] - yp[i]));
}
HorizontalIIRFilter<ushort>(img, im1, range, a, w);
break;
default :
return ;
......@@ -432,53 +364,38 @@ public:
void GradientPaillouY(InputArray _op, OutputArray _dst, double alpha, double omega)
{
Mat tmp(_op.size(),CV_32FC(_op.channels()));
_dst.create( _op.size(),CV_32FC(tmp.channels()) );
cv::Mat opSrc = _op.getMat();
std::vector<Mat> planSrc;
split(opSrc,planSrc);
split(_op,planSrc);
std::vector<Mat> planTmp;
split(tmp,planTmp);
std::vector<Mat> planDst;
split(_dst,planDst);
for (int i = 0; i < static_cast<int>(planSrc.size()); i++)
{
if (planSrc[i].isContinuous() && planTmp[i].isContinuous() && planDst[i].isContinuous())
for (int i = 0; i <static_cast<int>(planSrc.size()); i++)
{
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());
ParallelGradientPaillouYCols x(planSrc[i],planTmp[i],alpha,omega);
parallel_for_(Range(0,opSrc.cols), x,getNumThreads());
parallel_for_(Range(0, planSrc[i].cols), x,getNumThreads());
ParallelGradientPaillouYRows xr(planTmp[i],planDst[i],alpha,omega);
parallel_for_(Range(0,opSrc.rows), xr,getNumThreads());
}
else
std::cout << "PB";
parallel_for_(Range(0, planTmp[i].rows), xr,getNumThreads());
}
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<Mat> planSrc;
split(opSrc,planSrc);
split(_op,planSrc);
std::vector<Mat> planTmp;
split(tmp,planTmp);
std::vector<Mat> planDst;
split(_dst,planDst);
for (int i = 0; i < static_cast<int>(planSrc.size()); i++)
{
if (planSrc[i].isContinuous() && planTmp[i].isContinuous() && planDst[i].isContinuous())
for (int i = 0; i <static_cast<int>(planSrc.size()); i++)
{
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());
ParallelGradientPaillouXRows x(planSrc[i],planTmp[i],alpha,omega);
parallel_for_(Range(0,opSrc.rows), x,getNumThreads());
parallel_for_(Range(0, planSrc[i].rows), x,getNumThreads());
ParallelGradientPaillouXCols xr(planTmp[i],planDst[i],alpha,omega);
parallel_for_(Range(0,opSrc.cols), xr,getNumThreads());
}
else
std::cout << "PB";
parallel_for_(Range(0, planTmp[i].cols), xr,getNumThreads());
}
merge(planDst,_dst);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment