Почему код отслеживания исчезающих точек с помощью фильтров Gabor дает неверные результаты?

Я пытаюсь определить точку схождения изображения с помощью фильтров Габора. Я получаю 4 фильтра Габора с ориентациями {0, 45, 90, 135} и применяю к изображению, получая для отфильтрованных изображений. Затем я получаю их энергетические отклики и подавляю их, вычитая полученную матрицу с ее ортогональным откликом (отклик энергии 0 ° минус отклик энергии 90 ° и т. Д.), Получая 4 новые матрицы, 2 из подавленных откликов и 2 ориентации, соответствующие таковым у Габора. Фильтрует в зависимости от того, был ли результат предыдущего вычитания положительным или отрицательным. Я использую эти матрицы для вычисления доминирующей ориентации каждого пикселя в исходном изображении путем умножения подавленной реакции каждого пикселя на вектор, представляющий унарное комплексное число на основе матриц ориентации.
Доминирующие ориентации используются в схеме голосования, взвешенной функцией расстояния, чтобы получить точку схода, которая будет областью с наибольшим количеством голосов. Процесс показан в этом алгоритме:

Отслеживание точки схода с помощью фильтров Gabor

И это реализация, которую я имею в C ++:

#include <opencv2/opencv.hpp>
#include <string>
#include <vector>
#include <cmath>

cv::Mat image, imreduc, imvect, imeq, imblur, imedge, immap, imvert;
int theta = 4;
int ksize = 5;
int vlambda = 5.7;
double vsigma = 0.5;
double vgamma = 0.5;
double vpsi = 0;
double imscale = 0.8;

std::vector<cv::Mat> gabor_createR(int gsize, double stdev, double *orien, double wave, double asrat, double offset);
std::vector<cv::Mat> gabor_createI(int gsize, double stdev, double *orien, double wave, double asrat, double offset);
std::vector<cv::Mat> imfilteringR(std::vector<cv::Mat> fgaborR, cv::Mat imorigin);
std::vector<cv::Mat> imfilteringI(std::vector<cv::Mat> fgaborI, cv::Mat imorigin);
std::vector<cv::Mat> eresponse(std::vector<cv::Mat> fimagesR, std::vector<cv::Mat> fimagesI);
std::vector<cv::Mat> suppresor(std::vector<cv::Mat> energy);
std::vector<cv::Mat> suppdir(std::vector<cv::Mat> energy);
cv::Mat oldom(std::vector<cv::Mat> suppenergy, std::vector<cv::Mat> offseter, double *angles);
cv::Mat maxdist(cv::Mat fmatx, cv::Mat fmaty, cv::Mat dominor, cv::Mat dangles);
cv::Mat voter(cv::Mat dmax, cv::Mat dominor, cv::Mat dangles);

int main(int argc, char** argv){

std::string imageName("/home/jorge/testimage/testo0.jpg");
if(argc > 1){
imageName = argv[1];
}

image = cv::imread(imageName, CV_8UC1);
if(image.empty()){
std::cout << "Could not open or find the image" << std::endl;
return -1;
}
//cv::equalizeHist(image,imeq);
//cv::GaussianBlur(imeq, imblur, cv::Size(ksize,ksize), vsigma);
cv::Laplacian(image,imedge,-1);
//cv::subtract(255, imedge, imvert);
cv::resize(imedge, imreduc, cv::Size(), imscale, imscale, cv::INTER_AREA); //Image scaling for speed purposes
//Image dimensions for reference
int height = imreduc.rows;
int width = imreduc.cols;

//Array of orientations
double artheta[theta];
artheta[0] = 0;
for(int i = 1;i<theta;i++){
artheta[i] = i*CV_PI/theta;
}

//Function calls to create real and imagibary Gabor Filter banks
std::vector<cv::Mat> filtersR = gabor_createR(ksize,vsigma,artheta,vlambda,vgamma,vpsi);
std::vector<cv::Mat> filtersI = gabor_createI(ksize,vsigma,artheta,vlambda,vgamma,vpsi);

//Functions calls to convolve image with filter banks
std::vector<cv::Mat> imagebankR = imfilteringR(filtersR, imreduc);
std::vector<cv::Mat> imagebankI = imfilteringI(filtersI, imreduc);

//Function call to compute Gabor energy response
std::vector<cv::Mat> gaborenergy = eresponse(imagebankR, imagebankI);

//Function call to compute the suppressed energy response
std::vector<cv::Mat> suppgabor = suppresor(gaborenergy);

//Function call to compute angle of suppressed responses P.S: May be scrapped
std::vector<cv::Mat> predir = suppdir(gaborenergy);

//Function call for computing dominant orientations
cv::Mat dominor = oldom(suppgabor, predir, artheta);
cv::Mat compor = dominor*(180/CV_PI); //Conversion of orientations to degrees

cv::Mat matx = cv::Mat::zeros(compor.rows, compor.cols, CV_32S);
cv::Mat maty = cv::Mat::zeros(compor.rows, compor.cols, CV_32S);
for(int i = 0;i < compor.rows;i++){
for(int j = 0;j < compor.cols;j++){
matx.at<int>(i,j) = j;
maty.at<int>(i,j) = i;
}
}

int px,py, downsampler = 0;
imreduc.copyTo(imvect);
for(int i = 0;i < compor.rows;i++){
for(int j = 0;j < compor.cols;j++){
if(downsampler == 100){
px = (int) matx.at<int>(i,j) + 10 * cos(dominor.at<float>(i,j));
py = (int) maty.at<int>(i,j) + 10 * sin(dominor.at<float>(i,j));
cv::line(imvect,cv::Point(matx.at<int>(i,j),maty.at<int>(i,j)),cv::Point(px,py),cv::Scalar(100),1,8,0);
downsampler = 0;
}else{
downsampler++;
}
}
}
cv::namedWindow("direction test", CV_WINDOW_AUTOSIZE );
cv::imshow("direction test", imvect );
cv::waitKey(0);

cv::Mat fmatx, fmaty;
matx.convertTo(fmatx,CV_32F);
maty.convertTo(fmaty,CV_32F);

cv::Mat dmax = maxdist(fmatx, fmaty, dominor, compor);

cv::Mat urna = voter(dmax, dominor, compor);
cv::applyColorMap(urna,immap,cv::COLORMAP_JET);
cv::namedWindow("Display window", cv::WINDOW_AUTOSIZE); // Create a window for display.
cv::imshow("Display window", immap);                // Show our image inside it.
//cv::imshow("Display window", urna);                // Show our image inside it.
cv::waitKey(0); // Wait for a keystroke in the window
return 0;
}

std::vector<cv::Mat> gabor_createR(int gsize, double stdev, double *orien, double wave, double asrat, double offset){
std::vector<cv::Mat> filters;
cv::Mat kern;
cv::Mat aux;
float scaler;
double min, max;
for(int i = 0;i < theta;i++){
kern = cv::getGaborKernel(cv::Size(gsize, gsize), stdev, orien[i], wave, asrat, offset, CV_32F);
aux = cv::Mat::ones(kern.rows, kern.cols, CV_32F);
cv::minMaxLoc(kern, &min, &max);
cv::subtract(kern,min,kern);
cv::divide(kern,(max-min),kern);
//scaler = (float) max;
//scaler = 1.5*cv::sum(kern)[0];
//scaler = (float) 1;
//cv::multiply(kern, aux, kern, 1.0/(scaler));
filters.push_back(kern);
}

return filters;
}

std::vector<cv::Mat> gabor_createI(int gsize, double stdev, double *orien, double wave, double asrat, double offset){
std::vector<cv::Mat> filters;
cv::Mat kern;
cv::Mat aux;
float scaler;
double min, max;
for(int i = 0;i < theta;i++){
kern = cv::getGaborKernel(cv::Size(gsize, gsize), stdev, orien[i], wave, asrat, (offset + (CV_PI/2)), CV_32F);
aux = cv::Mat::ones(kern.rows, kern.cols, CV_32F);
cv::minMaxLoc(kern, &min, &max);
cv::subtract(kern,min,kern);
cv::divide(kern,(max-min),kern);
//scaler = (float) max;
//scaler = 1.5*cv::sum(kern)[0];
//scaler = (float) 1;
//cv::multiply(kern, aux, kern, 1.0/(scaler));
filters.push_back(kern);
}

return filters;
}

std::vector<cv::Mat> imfilteringR(std::vector<cv::Mat> fgaborR, cv::Mat imorigin){
std::vector<cv::Mat> imagebank;
cv::Mat aux, fimage;
cv::Mat bridge = cv::Mat::zeros(imorigin.rows, imorigin.cols, CV_8U);
for(int i = 0;i < theta;i++){
cv::filter2D(imorigin, aux, CV_32F, fgaborR[i]);
double min, max;
cv::minMaxLoc(aux, &min, &max);
cv::subtract(aux,min,aux);
cv::divide(aux,(max-min),aux);
aux *= 255;
aux.convertTo(bridge,CV_8U);
imagebank.push_back(bridge*1.0);
cv::namedWindow("Display window", cv::WINDOW_AUTOSIZE); // Create a window for display.
cv::imshow("Display window", bridge);                // Show our image inside it.
cv::waitKey(0); // Wait for a keystroke in the window
}
return imagebank;
}

std::vector<cv::Mat> imfilteringI(std::vector<cv::Mat> fgaborI, cv::Mat imorigin){
std::vector<cv::Mat> imagebank;
cv::Mat aux, fimage;
cv::Mat bridge = cv::Mat::zeros(imorigin.rows, imorigin.cols, CV_8U);
for(int i = 0;i < theta;i++){
cv::filter2D(imorigin, aux, CV_64F, fgaborI[i]);
double min, max;
cv::minMaxLoc(aux, &min, &max);
cv::subtract(aux,min,aux);
cv::divide(aux,(max-min),aux);
aux *= 255;
aux.convertTo(bridge,CV_8U);
imagebank.push_back(bridge*1.0);
cv::namedWindow("Display window", cv::WINDOW_AUTOSIZE); // Create a window for display.
cv::imshow("Display window", bridge);                // Show our image inside it.
cv::waitKey(0); // Wait for a keystroke in the window
}
return imagebank;
}

std::vector<cv::Mat> eresponse(std::vector<cv::Mat> fimagesR, std::vector<cv::Mat> fimagesI){
std::vector<cv::Mat> gaborenergy;
cv::Mat aux = cv::Mat::zeros(fimagesR[0].rows, fimagesR[0].cols, CV_32F);
cv::Mat real = cv::Mat::zeros(fimagesR[0].rows, fimagesR[0].cols, CV_8U);
cv::Mat imag = cv::Mat::zeros(fimagesR[0].rows, fimagesR[0].cols, CV_8U);
cv::Mat bridge = cv::Mat::zeros(fimagesR[0].rows, fimagesR[0].cols, CV_8U);
for(int i = 0;i < theta;i++){
fimagesR[i].convertTo(aux,CV_32F);
cv::pow(aux, 2.0, real);
fimagesI[i].convertTo(aux,CV_32F);
cv::pow(aux, 2.0, imag);
cv::sqrt((real + imag),aux);
aux.convertTo(bridge,CV_8U);
gaborenergy.push_back(aux*1.0);
cv::namedWindow("Display window", cv::WINDOW_AUTOSIZE); // Create a window for display.
cv::imshow("Display window", bridge);                // Show our image inside it.
cv::waitKey(0); // Wait for a keystroke in the window
}

return gaborenergy;
}

std::vector<cv::Mat> suppresor(std::vector<cv::Mat> energy){
std::vector<cv::Mat> suppgabor;
cv::Mat aux = cv::Mat::zeros(energy[0].rows, energy[0].cols, CV_32F);
cv::Mat bridge = cv::Mat::zeros(energy[0].rows, energy[0].cols, CV_8U);
for(int i = 0;i < theta/2;i++){
aux = cv::abs((cv::abs(energy[i])-cv::abs(energy[i+theta/2])));
//aux = cv::abs((cv::abs(energy[i])-cv::abs(energy[theta-(i+1)])));
//aux = cv::abs(energy[i])-cv::abs(energy[i+theta/2]);
aux.convertTo(bridge,CV_8U);
suppgabor.push_back(aux*1.0);
cv::namedWindow("Display window", cv::WINDOW_AUTOSIZE); // Create a window for display.
cv::imshow("Display window", bridge);                // Show our image inside it.
cv::waitKey(0); // Wait for a keystroke in the window
}
return suppgabor;
}

std::vector<cv::Mat> suppdir(std::vector<cv::Mat> energy){
std::vector<cv::Mat> endmat;
cv::Mat aux = cv::Mat::zeros(energy[0].rows, energy[0].cols, CV_32F);
cv::Mat mask, angmat;
for(int i = 0;i < theta/2;i++){
angmat = (cv::Mat::ones(energy[0].rows, energy[0].cols, CV_32F))*(i*CV_PI/theta);
mask = cv::Mat::zeros(energy[0].rows, energy[0].cols, CV_32F);
aux = cv::abs(energy[i])-cv::abs(energy[i+theta/2]);
//aux = cv::abs(energy[i])-cv::abs(energy[theta-(i+1)]);
cv::compare(aux,0,mask,cv::CMP_LT);
cv::add(angmat,(CV_PI/2),angmat,mask);
endmat.push_back(angmat);
}
return endmat;
}

cv::Mat oldom(std::vector<cv::Mat> suppenergy, std::vector<cv::Mat> predir, double *angles){
cv::Mat votimr = cv::Mat::zeros(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
cv::Mat votimi = cv::Mat::zeros(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
cv::Mat testo = cv::Mat::ones(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
cv::Mat aux, esubnr, esubni;
cv::Mat dominor = cv::Mat::zeros(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
for(int i = 0;i < theta/2;i++){
esubnr = cv::Mat::zeros(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
esubni = cv::Mat::zeros(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
aux = cv::Mat::zeros(suppenergy[0].rows, suppenergy[0].cols, CV_32F);
for(int j = 0;j < suppenergy[0].rows;j++){
for(int k = 0;k < suppenergy[0].cols;k++){
esubnr.at<float>(j,k) = cos(predir[i].at<float>(j,k));
esubni.at<float>(j,k) = sin(predir[i].at<float>(j,k));
}
}
cv::multiply(suppenergy[i],esubnr,aux);
//aux *= cos(angles[i]);
votimr += aux;
cv::multiply(suppenergy[i],esubni,aux);
//aux *= sin(angles[i]);
votimi += aux;
}
for(int i = 0;i < suppenergy[0].rows;i++){
for(int j = 0;j < suppenergy[0].cols;j++){
dominor.at<float>(i,j) = atan2(votimi.at<float>(i,j), votimr.at<float>(i,j));
dominor.at<float>(i,j) *= sin(dominor.at<float>(i,j));
}
}
return dominor;
}

cv::Mat maxdist(cv::Mat fmatx, cv::Mat fmaty, cv::Mat dominor, cv::Mat dangles){
cv::Mat cosdominor = cv::Mat::zeros(dangles.rows, dangles.cols, CV_32F);
cv::Mat sindominor = cv::Mat::zeros(dangles.rows, dangles.cols, CV_32F);
for(int i = 0;i < dangles.rows;i++){
for(int j = 0;j < dangles.cols;j++){
cosdominor.at<float>(i,j) = cos(dominor.at<float>(i,j));
sindominor.at<float>(i,j) = sin(dominor.at<float>(i,j));
}
}

cv::Mat ywidth, xheight, yxzero, xyzero, drborder, dlborder, duborder, ddborder, sub1, sub2, sum1;

cv::subtract(dangles.cols - 1,fmatx,sub1);
cv::divide(sub1,cosdominor,sub2);
drborder = cv::abs(sub2);
cv::multiply(drborder,sindominor,sum1);
cv::add(fmaty,sum1,ywidth);

cv::divide(fmatx,cosdominor,sub2);
dlborder = cv::abs(sub2);
cv::multiply(dlborder,sindominor,sum1);
cv::add(fmaty,sum1,yxzero);

cv::subtract(dangles.rows - 1,fmaty,sub1);
cv::divide(sub1,sindominor,sub2);
ddborder = cv::abs(sub2);
cv::multiply(ddborder,cosdominor,sum1);
cv::add(fmatx,sum1,xheight);

cv::divide(fmaty,sindominor,sub2);
duborder = cv::abs(sub2);
cv::multiply(duborder,cosdominor,sum1);
cv::add(fmatx,sum1,xyzero);

for(int i = 0;i < drborder.rows;i++){
for(int j = 0;j < drborder.cols;j++){
if(ywidth.at<float>(i,j) >= drborder.rows || ywidth.at<float>(i,j) < 0 || (dangles.at<float>(i,j) >= 90 && dangles.at<float>(i,j) <= 270)){
drborder.at<float>(i,j) = -1;
}
if(yxzero.at<float>(i,j) >= drborder.rows || yxzero.at<float>(i,j) < 0 || ((dangles.at<float>(i,j) >= 0 && dangles.at<float>(i,j) <= 90) || (dangles.at<float>(i,j) >= 270 && dangles.at<float>(i,j) <= 360))){
dlborder.at<float>(i,j) = -1;
}
if(xheight.at<float>(i,j) >= drborder.cols || xheight.at<float>(i,j) < 0 || ((dangles.at<float>(i,j) >= 180 && dangles.at<float>(i,j) <= 360) || dangles.at<float>(i,j) == 0)){
ddborder.at<float>(i,j) = -1;
}
if(xyzero.at<float>(i,j) >= drborder.cols || xyzero.at<float>(i,j) < 0 || (dangles.at<float>(i,j) >= 0 && dangles.at<float>(i,j) <= 180)){
duborder.at<float>(i,j) = -1;
}
}
}

cv::Mat dmax;
cv::max(drborder, cv::max(dlborder,cv::max(ddborder,duborder)),dmax);

return dmax;
}

cv::Mat voter(cv::Mat dmax, cv::Mat dominor, cv::Mat dangles){
cv::Mat urna = cv::Mat::zeros(dmax.rows, dmax.cols, CV_8U);
cv::Mat aux = cv::Mat::zeros(dmax.rows, dmax.cols, CV_32F);
int tx, ty;
int oldty;
double dist, shindist, exparg, fundist, suppor, vote;
for(int i = 0;i < dmax.rows;i++){
for(int j = 0;j < dmax.cols;j++){
suppor = fabs(sin(dominor.at<double>(i,j)));
//suppor = 1.0;
oldty = 10000;
if(dangles.at<float>(i,j) == 0 || dangles.at<float>(i,j) == 360){
for(tx = j;tx < dmax.cols;tx++){
dist = sqrt(pow(j - tx, 2.0));
shindist = dist/dmax.at<float>(i,j);
exparg = (-1*pow(shindist,2.0))/(2*pow(vsigma,2.0));
fundist = exp(exparg);
aux.at<float>(i,tx) += suppor*fundist;
}
}else if(dangles.at<float>(i,j) == 180){
for(tx = 0;tx < j;tx++){
dist = sqrt(pow(j - tx, 2.0));
shindist = dist/dmax.at<float>(i,j);
exparg = (-1*pow(shindist,2.0))/(2*pow(vsigma,2.0));
fundist = exp(exparg);
aux.at<float>(i,tx) += suppor*fundist;
}
}else if(dangles.at<float>(i,j) == 90){
for(ty = i;ty < dmax.rows;ty++){
dist = sqrt(pow(i - ty, 2.0));
shindist = dist/dmax.at<float>(i,j);
exparg = (-1*pow(shindist,2.0))/(2*pow(vsigma,2.0));
fundist = exp(exparg);
aux.at<float>(ty,j) += suppor*fundist;
}
}else if(dangles.at<float>(i,j) == 270){
for(ty = 0;ty <= i;ty++){
dist = sqrt(pow(i - ty, 2.0));
shindist = dist/dmax.at<float>(i,j);
exparg = (-1*pow(shindist,2.0))/(2*pow(vsigma,2.0));
fundist = exp(exparg);
aux.at<float>(ty,j) += suppor*fundist;
}
}else if((dangles.at<float>(i,j) > 0 && dangles.at<float>(i,j) < 90) || (dangles.at<float>(i,j) > 270 && dangles.at<float>(i,j) < 360)){
for(tx = j;tx < dmax.cols;tx++){
ty = (int)(tan(dominor.at<float>(i,j))*(tx - j) + i);
if(ty >= dmax.rows || ty < 0 || oldty == ty){
break;
}
dist = sqrt(pow(j - tx, 2.0) + pow(i - ty, 2.0));
shindist = dist/dmax.at<float>(i,j);
exparg = (-1*pow(shindist,2.0))/(2*pow(vsigma,2.0));
fundist = exp(exparg);
aux.at<float>(ty,tx) += suppor*fundist;
oldty = ty;
}
}else if(dangles.at<float>(i,j) > 90 && dangles.at<float>(i,j) < 270){
for(tx = 0;tx <= j;tx++){
ty = (int)(i - tan(dominor.at<float>(i,j))*(j - tx));
if(ty >= dmax.rows || ty < 0 || oldty == ty){
}else{
dist = sqrt(pow(j - tx, 2.0) + pow(i - ty, 2.0));
shindist = dist/dmax.at<float>(i,j);
exparg = (-1*pow(shindist,2.0))/(2*pow(vsigma,2.0));
fundist = exp(exparg);
aux.at<float>(ty,tx) += suppor*fundist;
oldty = ty;
}
}
}
}
}
double min, max;
cv::minMaxLoc(aux, &min, &max);
cv::subtract(aux,min,aux);
cv::divide(aux,(max-min),aux);
aux *= 255;
cv::Mat aux2;
cv::resize(aux,aux2,cv::Size(),(1.0/imscale),(1.0/imscale), cv::INTER_CUBIC);
aux2.convertTo(urna, CV_8U);

return urna;
}

* Закомментированные строки кода — это несколько вариантов, которые я пытался комбинировать в разных конфигурациях, чтобы увидеть, какая комбинация работала лучше всего

С этим алгоритмом я должен получить изображения, подобные этим трем:

Правильные результаты алгоритма

На этих изображениях есть одна относительно четко определенная небольшая область с большим количеством голосов, чем остальная часть изображения.

С моей реализацией я получаю что-то вроде этого:

Мои полученные изображения

На этом изображении я получаю несколько полос, которые предположительно представляют точку схода и идут по исходному изображению (Исходное изображение), они не

Основная проблема, которую я получаю, заключается в том, что доминирующие ориентации, которые я получаю, слишком однородны (слишком много пикселей с одинаковой ориентацией), чтобы не упомянуть, что они кажутся одной из ориентаций фильтров Габора. Это была часть кода, с которой я больше всего возился, но в результате доминирующие ориентации всегда слишком однородны.

Какая часть (или части) в моем коде нуждается в изменениях для получения правильных результатов?