Log-Gabor濾波器構造，opencv實現
阿新 • • 發佈：2018-12-30
參照連結：https://github.com/carandraug/PeterKovesiImage/blob/master/PhaseCongruency/gaborconvolve.m點選開啟連結
畢設要用Log-Gabor濾波器來實現視網膜血管增強，這東西真是折騰了我好一段時間。。。。Matlab程式碼轉到C++程式碼沒動什麼腦子，基本就是按邏輯翻譯過來的
Log-Gabor濾波器的概念本文暫不表述了，待後期新增
首先是一些與濾波器無關的matlab轉C++用 的程式碼：
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); }
下文需要用到的是meshgrid函式、fftshift、ifftshift函式，這兩個函式是幹嘛用的自行百度
mathlab裡[x,y] = meshgrid(xrange, yrange); 等價於上述的 meshgrid（xrange.lower,xrange.upper,yrange.lower,yrange.upper）//大意虛擬碼
然後是Log-Gabor核構造需要的低通濾波器
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; }
返回的mat即為構造的低通濾波器核
然後是一些矩陣運算操作
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); }
注意，之所以是Mat[2]是因為 傳入的應是複數，mat[0] [1]分別為實部和虛部
上述三個函式分別是計算 × ÷ 和 強度的（平方和開根）
最後是返回結果的結構體的定義
struct GaborConvolveResult
{
vector<vector<Mat>>EO;
vector<Mat>BP;
};EO為 EO[nOrient,nScale]，分別代表了不同角度和不同尺度的處理結果
BP為BP[nScale] ，分別表示了不同尺度下的處理結果（方向已融合）
下面就是Log-Gabor核的具體構造了
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;
}\