# Log-Gabor滤波器构造，opencv实现

5 篇文章 0 订阅
3 篇文章 0 订阅

Log-Gabor滤波器的概念本文暂不表述了，待后期添加

void circshift(Mat& out, const Point& delta)
{
Size sz = out.size();

// 错误检查
assert(sz.height > 0 && sz.width > 0);
// 此种情况不需要移动
if ((sz.height == 1 && sz.width == 1) || (delta.x == 0 && delta.y == 0))
return;

// 需要移动参数的变换，这样就能输入各种整数了
int x = delta.x;
int y = delta.y;
if (x > 0) x = x % sz.width;
if (y > 0) y = y % sz.height;
if (x < 0) x = x % sz.width + sz.width;
if (y < 0) y = y % sz.height + sz.height;

// 多维的Mat也能移动
vector<Mat> planes;
split(out, planes);

for (size_t i = 0; i < planes.size(); i++)
{
// 竖直方向移动
Mat tmp0, tmp1, tmp2, tmp3;
Mat q0(planes[i], Rect(0, 0, sz.width, sz.height - y));
Mat q1(planes[i], Rect(0, sz.height - y, sz.width, y));
q0.copyTo(tmp0);
q1.copyTo(tmp1);
tmp0.copyTo(planes[i](Rect(0, y, sz.width, sz.height - y)));
tmp1.copyTo(planes[i](Rect(0, 0, sz.width, y)));

// 水平方向移动
Mat q2(planes[i], Rect(0, 0, sz.width - x, sz.height));
Mat q3(planes[i], Rect(sz.width - x, 0, x, sz.height));
q2.copyTo(tmp2);
q3.copyTo(tmp3);
tmp2.copyTo(planes[i](Rect(x, 0, sz.width - x, sz.height)));
tmp3.copyTo(planes[i](Rect(0, 0, x, sz.height)));
}

merge(planes, out);
}

void fftshift(Mat& out)
{
Size sz = out.size();
Point pt(0, 0);
pt.x = (int)floor(sz.width / 2.0);
pt.y = (int)floor(sz.height / 2.0);
circshift(out, pt);
}

void ifftshift(Mat& out)
{
Size sz = out.size();
Point pt(0, 0);
pt.x = (int)ceil(sz.width / 2.0);
pt.y = (int)ceil(sz.height / 2.0);
circshift(out, pt);
}

void meshgrid(double xStart, double xEnd, double yStart, double yEnd,
Mat& X, Mat& Y)
{
std::vector<double> t_x, t_y;
while (xStart <= xEnd)
{
t_x.push_back(xStart);

++xStart;
}
while (yStart <= yEnd)
{
t_y.push_back(yStart);

++yStart;
}

repeat(Mat(t_x).t(), t_y.size(), 1, X);
repeat(Mat(t_y), 1, t_x.size(), Y);
}

mathlab里[x,y] = meshgrid(xrange, yrange); 等价于上述的 meshgrid（xrange.lower,xrange.upper,yrange.lower,yrange.upper）//大意伪代码

Mat lowpassFilter(Size size, double cutOff, int n)
{
int rows, cols;

double xStart, xEnd, yStart, yEnd;

if (cutOff < 0 || cutOff > 0.5)
throw std::exception("cutoff frequency must be between 0 and 0.5");

if (size.width == 1)
{
rows = size.width;
cols = size.width;
}
else
{
rows = size.height;
cols = size.width;
}

Mat x(cols, rows, CV_64FC1);
Mat y(cols, rows, CV_64FC1);

bool isRowsOdd = (rows % 2 == 1);
bool isColsOdd = (cols % 2 == 1);
if (isRowsOdd)
{
yStart = -(rows - 1) / 2.0;
yEnd = (rows - 1) / 2.0;
}
else
{
yStart = -(rows / 2.0);
yEnd = (rows / 2.0 - 1);
}

if (isColsOdd)
{
xStart = -(cols - 1) / 2.0;
xEnd = (cols - 1) / 2.0;
}
else
{
xStart = -(cols / 2.0);
xEnd = (cols / 2.0 - 1);
}

meshgrid(xStart, xEnd, yStart, yEnd, x, y);
if (isColsOdd)
x /= (cols - 1);
else
x /= cols;

if (isRowsOdd)
y /= (rows - 1);
else
y /= rows;

Mat radius;
Mat x2;
Mat y2;

pow(x, 2, x2);
pow(y, 2, y2);
sqrt(x2 + y2, radius);
Mat f;
Mat temp;
pow((radius / cutOff), 2 * n, temp);
divide(Mat::ones(radius.rows, radius.cols, CV_64F), 1.0 + temp, f);
ifftshift(f);
return f;
}

void complexMul(Mat left[2], Mat right[2], Mat* result)
{
Mat real1;
Mat real2;
Mat imag1;
Mat imag2;

real1 = left[0].mul(right[0]);
real2 = (left[1].mul(right[1]));
imag1 = left[1].mul(right[0]);
imag2 = left[0].mul(right[1]);

result[0] = real1 - real2;
result[1] = imag1 + imag2;
}

void complexDiv(Mat left[2], Mat right[2], Mat* result)
{
Mat negRight[2];
negRight[0] = right[0];
negRight[1] = -right[1];

Mat upper[2];
Mat lower[2];
complexMul(left, negRight, upper);
complexMul(right, negRight, lower);

divide(upper[0], lower[0], result[0]);
divide(upper[1], lower[0], result[1]);
}

void complexAbs(Mat mat[2], Mat& result)
{
Mat left, right;
cv::pow(mat[0], 2, left);
pow(mat[1], 2, right);
sqrt(left + right, result);
}

struct GaborConvolveResult
{

vector<vector<Mat>>EO;
vector<Mat>BP;
};

EO为 EO[nOrient,nScale]，分别代表了不同角度和不同尺度的处理结果

BP为BP[nScale] ，分别表示了不同尺度下的处理结果（方向已融合）

GaborConvolveResult garborConvolve(const Mat& mat, int nScale, int nOrient, double minWaveLength, double mult, double sigmaOnf,
double dThetaSigma, int Lnorm = 0, double feedback = 0)
{
//mat应已确认是灰度图,并且是double
int rows = mat.rows;
int cols = mat.cols;

Mat matDft;
toFre(mat, matDft);

vector<vector<Mat>>EO(nOrient, vector<Mat>(nScale, Mat(cols, rows, CV_64F,Scalar(0))));
vector<Mat>BP(nScale, Mat(cols, rows, CV_64F, Scalar(0)));

Mat x;
Mat y;
double xStart, xEnd, yStart, yEnd;
bool isRowsOdd = (rows % 2 == 1);
bool isColsOdd = (cols % 2 == 1);
if (isRowsOdd)
{
yStart = -(rows - 1) / 2.0;
yEnd = (rows - 1) / 2.0;
}
else
{
yStart = -(rows / 2.0);
yEnd = (rows / 2.0 - 1);
}

if (isColsOdd)
{
xStart = -(cols - 1) / 2.0;
xEnd = (cols - 1) / 2.0;
}
else
{
xStart = -(cols / 2.0);
xEnd = (cols / 2.0 - 1);
}

meshgrid(xStart, xEnd, yStart, yEnd, x, y);

if (isColsOdd)
x /= (cols - 1);
else
x /= cols;

if (isRowsOdd)
y /= (rows - 1);
else
y /= rows;

Mat radius;
Mat x2;
Mat y2;

pow(x, 2, x2);
pow(y, 2, y2);
sqrt(x2 + y2, radius);

Mat theta(rows, cols, CV_64FC1);

for (int i = 0; i < rows; ++i)
for (int j = 0; j < cols; ++j)
theta.at<double>(i, j) = atan2(y.at<double>(i, j), x.at<double>(i, j));

ifftshift(radius);
ifftshift(theta);

radius.at<double>(0, 0) = 1;

Mat sinTheta(rows, cols, CV_64FC1);
Mat cosTheta(rows, cols, CV_64FC1);

for (int i = 0; i < rows; ++i)
for (int j = 0; j < cols; ++j)
sinTheta.at<double>(i, j) = sin(theta.at<double>(i, j));

for (int i = 0; i < rows; ++i)
for (int j = 0; j < cols; ++j)
cosTheta.at<double>(i, j) = cos(theta.at<double>(i, j));

Mat lp = lowpassFilter(Size(cols,rows), 0.45, 15);

//vector<Mat>logGabor(nScale, Mat(rows, cols, CV_64F,Scalar(0,0)));
vector<Mat>logGabor;
for(int  s = 0;s<nScale;++s)
{
logGabor.push_back(Mat(rows, cols, CV_64F));
double waveLength = minWaveLength* pow(mult, s );
double fo = 1.0 / waveLength;

Mat tempUpper;

log(radius / fo, tempUpper);
pow(tempUpper, 2, tempUpper);

double tempLower = pow(log(sigmaOnf), 2);

double factory = -1 / 2.0;
tempUpper = tempUpper / tempLower * factory;

exp(tempUpper, logGabor[s]);

logGabor[s] = logGabor[s].mul(lp);
logGabor[s].at<double>(0, 0) = 0;

double L = 0;
switch (Lnorm)
{
case 0:
L = 1;
break;
case 1:
{
Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
Mat complex;
idft(logGabor[s], complex, DFT_COMPLEX_OUTPUT + DFT_SCALE);
split(complex, planes);
Mat realPart = planes[0];

L = sum(abs(realPart))[0];
break;
}

case 2:
{
Mat temp;
pow(logGabor[s], 2, temp);

L = sqrt(sum(temp)[0]);

}

break;
default:
break;
}

logGabor[s] = logGabor[s] / L;
//cout << logGabor[s] << endl;
//cout << curLogGabor;

Mat complex;
Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
split(matDft, planes);

planes[0] = planes[0].mul(logGabor[s]);
planes[1] = planes[1].mul(logGabor[s]);

//idft(complex, EO, DFT_COMPLEX_OUTPUT + DFT_SCALE);

//cout << temp.depth() << "  " << temp.channels();
//split(temp, planes);

Mat complexd;
merge(planes, 2, complexd);
idft(complexd,BP[s], DFT_COMPLEX_OUTPUT + DFT_SCALE);

}
//cout << logGabor[0] << endl;

for(int o = 0 ; o < nOrient;++o)
{
double angl = o*CV_PI / nOrient;
double waveLength = minWaveLength;

Mat ds = sinTheta * cos(angl) - cosTheta * sin(angl);
Mat dc = cosTheta * cos(angl) + sinTheta * sin(angl);

Mat dTheta(rows, cols, CV_64F);
for (int i = 0; i < rows; ++i)
for (int j = 0; j < cols; ++j)
dTheta.at<double>(i, j) = abs(atan2(ds.at<double>(i, j), dc.at<double>(i, j)));

Mat temp;
pow(dTheta, 2, temp);
temp = -temp;
Mat spread;
double thetaSigma = CV_PI / nOrient / dThetaSigma;
exp(temp / (2 * pow(thetaSigma, 2)), spread);

for(int s =0 ;s<nScale;++s)
{
Mat filter = spread.mul(logGabor[s]);
double L = 0;
switch (Lnorm)
{
case 0: L = 1;
break;
case 1:
{
Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
Mat complex;
idft(filter, complex, DFT_COMPLEX_OUTPUT + DFT_SCALE);
split(complex, planes);
Mat realPart = planes[0];
L = sum(abs(realPart))[0];
}
break;
case 2:
{

Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };

split(temp, planes);
Mat realPart = planes[0];

Mat imagPart = planes[1];
pow(realPart, 2, realPart);
pow(imagPart, 2, imagPart);

L = sqrt(sum(realPart)[0] + sum(imagPart)[0]);
}
break;
default:
break;
}
filter = filter / L;

Mat complex;
Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
cv::split(matDft, planes);

planes[0] = planes[0].mul(filter);
planes[1] = planes[1].mul(filter);

merge(planes,2, complex);

//here
//Mat  multed = matDft.mul(filter);
//cout << filter << endl;
idft(complex, EO[o][s], DFT_COMPLEX_OUTPUT + DFT_SCALE);
//cout << EO[s][o].cols << " " << EO[s][o].rows << EO[s][o].channels() << " " << EO[s][o].depth() << endl;

Mat EOPlanes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
split(EO[o][s], EOPlanes);

waveLength = waveLength *mult;
}

}

GaborConvolveResult result;
result.BP = BP;
result.EO = EO;
return result;
}

\

• 0
点赞
• 5
评论
• 4
收藏
• 打赏
• 扫一扫，分享海报

04-02
11-13 5657
10-26 1079
09-21
11-04 1253
03-23 2541
06-26 2044
06-24
©️2022 CSDN 皮肤主题：编程工作室 设计师：CSDN官方博客

MaloryVer9

¥2 ¥4 ¥6 ¥10 ¥20

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、C币套餐、付费专栏及课程。