# LogGabor function problem

• (4 Pages)
• « First
• 2
• 3
• 4

## 52 Replies - 1596 Views - Last Post: 29 July 2013 - 02:19 AMRate Topic: //<![CDATA[ rating = new ipb.rating( 'topic_rate_', { url: 'http://www.dreamincode.net/forums/index.php?app=forums&module=ajax&section=topics&do=rateTopic&t=324141&amp;s=f24713e21cd8bf1dc9a3a1c80d27ebb6&md5check=' + ipb.vars['secure_hash'], cur_rating: 0, rated: 0, allow_rate: 0, multi_rate: 1, show_rate_text: true } ); //]]>

### #46 jimblumberg

Reputation: 4631
• Posts: 14,541
• Joined: 25-December 09

## Re: LogGabor function problem

Posted 17 July 2013 - 06:38 PM

After a bit of detective work I see that in the following snippet you're accessing your Y array out of bounds.

```void pre_filter_Computations(vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows){
...

for(int z=0;z<cols;z++){
for(int j=0;j<rows;j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2)));
}
Y.push_back(y);
y.clear();
}

for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){

X[a][b] = pow(X[a][b],2);
Y[a][b] = pow(Y[a][b],2); // PROBLEM HERE.
X[a][b] = X[a][b] + Y[a][b];
}
}
```

For the picture you posted:
```Rows: 534 Cols: 696
X.size() 534 X[0].size() 696
Y.size() 696 Y[0].size() 534
```

So every time you access Y with a value of greater than 533 you are accessing the vector out of range. Why did you reverse the cols and rows for your Y vector in the first set of loops?

By the way I found this issue by using the .at() vector member function instead of the array notation [].

```          Y.at(a).at(B)/> = pow(Y.at(a).at(B)/>, 2);
```

The at() member function throws an exception if you try to access the vector out of range.

Jim

### #47 residentelvio

Reputation: 0
• Posts: 30
• Joined: 26-June 13

## Re: LogGabor function problem

Posted 18 July 2013 - 04:05 AM

Just thinking over that for my picture :

```//Rows: 534 Cols: 696
//x is a vector of cols --> x[696]
//y is a vector of rows --> y[534]
//X is a vector of vectors that for meshgrid will be --> X[696][534] ( x elements goes in X's rows elements)
//Y is a vector of vectors that for meshgrid will be --> Y[696][534] ( y elements goes in Y's cols elements)
//radius[696][534] so it' s possible to sum XandY to create it.

```

Now just changed into fuction:

``` for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){
X[a][b] = pow(X[a][b],2);
Y[b][a] = pow(Y[b][a],2); // CHANGED HERE.
}
}

```

That' s the only way it works. In general it s wrong. Beacause X and Y must be egual about size. In this way it comes back to the main and enter to the next function logGabor. And there it' s stopped. Another change I need to do it 's radius and theta's declaration that must be the opposite for the size for what I was thinking.
Howewer using what you said I suppose that's the way.

Do have I to use push_back for radius bacause it' s empty till this for innested loop?

```for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){
X.at(a).at(B)/>/> = pow(x.at(a).X(B)/>/>, 2);
Y.at(a).at(B)/>/> = pow(Y.at(a).at(B)/>/>, 2);
Y[b][a] = pow(Y[b][a],2); // CHANGED HERE.
X[a][b] = X[a][b] + Y[a][b];
}
}

```

### #48 jimblumberg

Reputation: 4631
• Posts: 14,541
• Joined: 25-December 09

## Re: LogGabor function problem

Posted 18 July 2013 - 05:45 AM

Quote

Do have I to use push_back for radius bacause it' s empty till this for innested loop?

Unless you've changed things radius is not "empty".
```       vector< vector<double> >radius(rows,vector <double> (cols));
```

This created a vector with blank elements with a size of [row][cols].

But this will still access your Y out of range:
```            radius[a][b] = sqrt(X[a][b]+Y[a][b]);
```

I really think the major problem is in this snippet.
```  for(int i=0;i<rows;i++){
for(int k=0;k<cols;k++){
x.push_back(((double)k-(cols/2))/(double)(cols/2));
}
X.push_back(x);
x.clear();
}

for(int z=0;z<cols;z++){
for(int j=0;j<rows;j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2)));
}
Y.push_back(y);
y.clear();
}
```

Why did you reverse the cols and rows in the second loop? The size of your picture didn't change it's still [rows][cols] not the [cols][rows] that you're creating in the second loop. It should probably be this instead.

```  for(int z = 0; z < rows; z++){
for(int j = 0; j < cols; j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2)));
}
Y.push_back(y);
y.clear();
}
```

Jim

### #49 residentelvio

Reputation: 0
• Posts: 30
• Joined: 26-June 13

## Re: LogGabor function problem

Posted 18 July 2013 - 07:32 AM

If I did as you say that snippet block the compilation inside this function. In my version it s blocked in the next function.
All code if you want to try as you say:

```#include <stdio.h>
#include <cxcore.h>
#include <cv.h>
#include <highgui.h>
#include <math.h>
#include <vector>

using namespace cv;
using namespace std;

#define PI 3.14159265

struct ts{
Mat data;
}tiese[4];

Mat colorDecorrelation(Mat img,int numColComp);
Mat FourierTrasform2D(Mat img, double Lambda_max,Mat imfft,int *padSize,int *sz2filt_r,int *sz2filt_c);
int roundIt(double num);
void pre_filter_Computations(vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows);
Mat logGabor(vector< vector<double> > &filter,Mat imfft,double r_o, double theta_o,double sigma_theta,vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows,double sigma_r,int *padSize);
Mat whitening_sqr_meanDistance(Mat X);
void ftshifting(vector< vector<double> > &filter,int rows,int cols);

int main(){

double achrScales,achrOrient,achrLambda_min,achrScFactor,achrDeltaAng,achrBwF,achrLambda_max,chrScales,chrOrient,chrLambda_min,chrScFactor,chrDeltaAng,chrBwF,chrLambda_max,sigma_r,wavelengths,sigma_theta,sigma,r_o,theta_o;
Mat rgbDecorrImg,imfft,imfft2,imgDecorr,conspic1,conspic,conspic2,Tsquared,inCompImg,Salmap,SumSalmap,maxT,GaussKernel;
vector<Mat>stack;

if(img.empty()){
cout<<"The image is empty \n";
return -1;
}

namedWindow("Image Display", CV_WINDOW_AUTOSIZE);
imshow("Image Display", img);
waitKey(0);
destroyWindow("Image Display");

sz2filt_r = new int;
sz2filt_c = new int;

rows=img.rows;
cols= img.cols;
nChannels= img.channels();
numColComp = nChannels ;

inCompImg = img.reshape(numColComp,rows*cols);

imgDecorr = colorDecorrelation(inCompImg,numColComp);

rows=imgDecorr.rows;
cols=imgDecorr.cols;
nChannels= imgDecorr.channels();

imgDecorr.convertTo(rgbDecorrImg,CV_8UC3);
rgbDecorrImg = rgbDecorrImg.reshape(numColComp,img.rows);

namedWindow("Image Display", CV_WINDOW_AUTOSIZE);
imshow("Image Display", rgbDecorrImg);
waitKey(0);
destroyWindow("Image Display");

if(numColComp==1){
achrScales=7.0;
achrOrient=4.0;
achrLambda_min=3.0;
achrScFactor=2.0;
achrDeltaAng = 1.0;
achrBwF=1.1;
achrLambda_max= (achrLambda_min) * (pow(achrScFactor,achrScales-1));
}
else {
chrScales=5.0;
chrOrient=4.0;
chrLambda_min=6.0;
chrScFactor=2.0;
chrDeltaAng = 1.7;
chrBwF=1.1;
chrLambda_max= (chrLambda_min) * (pow(chrScFactor,chrScales-1));
}

if(numColComp==1){

rows = *sz2filt_r;
cols = *sz2filt_c;

vector< vector<double> >theta(rows,vector<double>(cols));
vector< vector<double> >filter(rows,vector<double>(cols));

sigma_r = exp(-(pow(2,achrBwF -1))/(pow(2,achrBwF + 1)));
sigma_theta = (PI/achrOrient)/achrDeltaAng;

for(int i_orient=1;i_orient<=achrOrient;i_orient++){
for(int i_scale=1;i_scale<=achrScales;i_scale++){

wavelengths = achrLambda_min * (pow(achrScFactor,i_scale-1));
r_o = 2/wavelengths;
theta_o = (i_orient-1)*(PI/achrOrient);

}
}
conspic1 = Mat::zeros(img.rows,img.cols,CV_8UC3);

for(int i_orient=1;i_orient<=achrOrient;i_orient++){

inCompImg = Mat(stack).reshape(1).t();
inCompImg = inCompImg.reshape(0,rows*cols);
Tsquared=whitening_sqr_meanDistance(inCompImg);
conspic = Tsquared.reshape(0,rows);

divide(conspic,max(0,conspic),conspic);
maxhw = max(0,(int)floor(min(rows/2,cols/2)-1));
sigma = max(rows,cols)*(0.025);
GaussianBlur(conspic,conspic,Size(maxhw,maxhw),sigma,sigma,BORDER_DEFAULT);

conspic1 = conspic1 + conspic;
}

divide(conspic1,max(0,conspic1),conspic1);
GaussianBlur(conspic,conspic,Size(maxhw,maxhw),sigma,sigma,BORDER_DEFAULT);

}

if(numColComp!=1){

conspic2 = Mat::zeros(img.rows,img.cols,CV_8UC3);

for(int i_c_comp=2;i_c_comp<numColComp;i_c_comp++){

rows=*sz2filt_r;
cols=*sz2filt_c;

imfft.convertTo(imfft2,CV_8UC3);
imfft2= imfft2.reshape(numColComp,rows);

namedWindow("Image Display", CV_WINDOW_AUTOSIZE);
imshow("Image Display", imfft2);
waitKey(0);
destroyWindow("Image Display");

vector< vector<double> >theta(cols,vector <double> (rows));
vector< vector<double> >filter(cols,vector <double> (rows));

sigma_r = exp(-(pow(2,chrBwF) -1)/(pow(2,chrBwF) + 1));
sigma_theta = (PI/chrOrient)/chrDeltaAng;

for(int i_orient=1;i_orient<=chrOrient;i_orient++){
for(int i_scale=1;i_scale<=chrScales;i_scale++){
Mat prova;
wavelengths = chrLambda_min * (pow(chrScFactor,i_scale -1));
r_o = 2/wavelengths;
theta_o = (i_orient-1)*(PI/chrOrient);
cout<< "Before Loggabor\n";

}
}

conspic = Mat::zeros(img.rows,img.cols,CV_8UC3);

for(int i_orient=1;i_orient<=achrOrient;i_orient++){

inCompImg = Mat(stack).reshape(numColComp).t();
inCompImg = inCompImg.reshape(0,rows*cols);
Tsquared=whitening_sqr_meanDistance(inCompImg);
Tsquared.convertTo(img,CV_8UC3);
img = img.reshape(numColComp,rows);

tiese[i_orient-1].data= Tsquared.reshape(0,rows);

for(int a=1;a<=chrOrient;a++){
maxT=max(0,tiese[a-1].data);
}

tiese[i_orient-1].data= tiese[i_orient-1].data/ maxT;

maxhw = max(0,(int)floor(min(rows/2,cols/2)-1));
sigma = max(rows,cols)*(0.025);
GaussianBlur(tiese[i_orient-1].data,tiese[i_orient-1].data,Size(maxhw,maxhw),sigma,sigma,BORDER_DEFAULT);

conspic=conspic + tiese[i_orient-1].data;

}
divide(conspic,max(0,conspic),conspic);

GaussianBlur(conspic,conspic,Size(maxhw,maxhw),sigma,sigma,BORDER_DEFAULT);

conspic2 = conspic2 + conspic;
}

multiply(conspic1,2,Salmap);
Salmap= Salmap + conspic2;
SumSalmap= sum(Salmap);
Salmap = Salmap / SumSalmap;

namedWindow("Salmap Display", CV_WINDOW_AUTOSIZE);
imshow("Salmap Display", Salmap);
waitKey(0);
destroyWindow("Salmap Display");

return 0;
}

}

Mat colorDecorrelation(Mat img,int numColComp){

Mat V,D,A,Xdecorr,img_tr,img_f;
Scalar mu;
SVD svd;

if(numColComp>1){

img.convertTo(img_f,CV_32FC3);
img_f=img_f.reshape(1,img.rows*img.cols);

mu=mean(img_f);
subtract(img_f,mu,img_f);
transpose(img_f,img_tr);
A = img_tr * img_f;
V = svd(A).u;
Xdecorr = img_f*V;

}
else
Xdecorr=img;

return Xdecorr;
}

Mat FourierTrasform2D(Mat img, double Lambda_max,Mat imfft,int *padSize,int *sz2filt_r,int *sz2filt_c){

namedWindow("Image Display", CV_WINDOW_AUTOSIZE);
waitKey(0);
destroyWindow("Image Display");

return imfft;

}

int roundIt(double num){
int numero;
double diff = fabs(num - (int)num);

if((diff >= 0.5 && num < 0) || (diff < 0.5 && num >= 0))
numero=(int)num;

else if(diff >=0.5 && num >= 0)
numero=(int)num + 1;

else if(diff >0.5 && num < 0)
numero=(int)num - 1;

return numero;
}

void pre_filter_Computations(vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows){

vector<double>x;
vector<double>y;
vector< vector<double> > X;
vector< vector<double> > Y;
double epsilon = 0.0001;

for(int i=0;i<rows;i++){
for(int k=0;k<cols;k++){
x.push_back(((double)k-(cols/2))/(double)(cols/2));
}
X.push_back(x);
x.clear();
}

for(int z=0;z<cols;z++){
for(int j=0;j<rows;j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2)));
}
Y.push_back(y);
y.clear();
}

for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){

X[a][b] = pow(X[a][b],2);
Y[b][a] = pow(Y[b][a],2);
}
}

for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){
theta[a][b] = atan2(Y[a][b],X[a][b])*180/PI;
}
}
X.clear();
Y.clear();
}

Mat whitening_sqr_meanDistance(Mat X){

Mat X_f,X_tr,A,V,D,S,Z,Xwh,whMat;
Scalar mu,mu_wh;
SVD svd;
double epsilon = 0.0001;
int r,c;

X.convertTo(X_f,CV_32FC3);
X_f=X_f.reshape(1,X.rows*X.cols);

mu = mean(X_f);
subtract(X_f,mu,X_f);
transpose(X_f,X_tr);
A=X_tr * X_f;

V = svd(A).u;
D = svd(A).w;

r= X.rows;
c= 1;

S= Mat::eye((int)r*epsilon,(int)c*epsilon,CV_32FC3);
S=S*epsilon;
Z=(D+S);
Z=Z.inv(DECOMP_LU);
sqrt(Z,S);
whMat=V*S;

Xwh=X*whMat;
mu_wh=mean(Xwh);
subtract(Xwh,mu_wh,Xwh);

Z= Mat::zeros(r,c,CV_32FC3);

sqrt(Xwh,Z);

return Z;
}

Mat logGabor(vector< vector<double> > &filter,Mat imfft,double r_o, double theta_o,double sigma_theta,vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows,double sigma_r,int *padSize){

vector< vector<double> >sintheta(rows,vector <double> (cols));
vector< vector<double> >costheta(rows,vector <double> (cols));
vector< vector<double> >dalpha(rows,vector <double> (cols));
vector< vector<double> >ds(rows,vector <double> (cols));
vector< vector<double> >dc(rows,vector <double> (cols));
vector< vector<double> >divis2(rows,vector <double> (cols));
vector< vector<double> >divis(rows,vector <double> (cols));

cout<< "Inside Loggabor\n";

for(int a=0; a<rows; a++){
for(int b=0; b<cols; b++){

sintheta[a][b]=sin(theta[a][b]);
costheta[a][b]=cos(theta[a][b]);
ds[a][b] = (sintheta[a][b] * cos(theta_o) - costheta[a][b] * sin(theta_o));
dc[a][b] = (costheta[a][b] * cos(theta_o) + sintheta[a][b] * sin(theta_o));
dalpha[a][b] = abs(atan2(ds[a][b],dc[a][b]));
divis[a][b] = -(dalpha[a][b]*dalpha[a][b])/(2*(sigma_theta*sigma_theta));
filter[a][b]= exp(divis2[a][b]);

}
}

ftshifting(filter,rows,cols);

for(int m=0;m<rows;m++){
for(int n=0;n<cols;n++) {

imfft =  imfft.mul(filter[m][n]);

}
}

}

void ftshifting(vector< vector<double> > &filter,int rows,int cols){

double tmp1_3,tmp2_4;
int r2,c2;

r2 = (rows+1)/2;
c2 = (cols+1)/2;

for (int i = 0; i<rows/2; i++){
for (int k = 0; k<cols/2; k++){

tmp1_3 = filter[i][k];
filter[i][k] = filter[i+r2][k+c2];
filter[i+r2][k+c2] = tmp1_3;

tmp2_4 = filter[i+r2][k];
filter[i+r2][k] = filter[i][k+c2];
filter[i][k+c2] = tmp2_4;
}
}
}

```

Now my question is:

```for(int i=0;i<rows;i++){
for(int k=0;k<cols;k++){
x.push_back(((double)k-(cols/2))/(double)(cols/2)); //creating vector x(cols)
}
X.push_back(x); //Creating X(cols,rows)
x.clear();
}

for(int z=0;z<cols;z++){
//maybe here I have to built y vector
for(int j=0;j<rows;j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2))); //it creates vector y(rows);
}
Y.push_back(y);   // does it creates vector of vectors Y(cols,rows)? I m just thinking to have Y[cols][rows] maybe y creation must be just inside first loop because it s a copy of second part of vector of vectors
y.clear();
}
}

```

This post has been edited by residentelvio: 18 July 2013 - 07:33 AM

### #50 jimblumberg

Reputation: 4631
• Posts: 14,541
• Joined: 25-December 09

## Re: LogGabor function problem

Posted 18 July 2013 - 09:22 AM

Quote

If I did as you say that snippet block the compilation inside this function. In my version it s blocked in the next function.

It doesn't matter where it is crashing the problem is the same, your trying to access one of your vectors out of bounds.

I'll ask again why are you reversing the cols and rows when you compute your y vector values?

Did your picture magically transform it's self into different dimensions?

```
for(int i=0;i<rows;i++){
for(int k=0;k<cols;k++){
x.push_back(((double)k-(cols/2))/(double)(cols/2));
}
X.push_back(x);
x.clear();
}
// Did you picture change here?
for(int z=0;z<cols;z++){
for(int j=0;j<rows;j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2)));
}
Y.push_back(y);
y.clear();
}

```

In my opinion this should be done in one loop.

```void pre_filter_Computations(vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows)
{

vector<double>x;
vector<double>y;
vector< vector<double> > X;
vector< vector<double> > Y;
double epsilon = 0.0001;

for(int i=0; i<rows; i++)
{
for(int k=0; k<cols; k++)
{
x.push_back(((double)k-(cols/2))/(double)(cols/2));
y.push_back(-(((double)i-(rows/2))/(double)(rows/2)));
}
X.push_back(x);
Y.push_back(y);
x.clear();
y.clear();
}

for(int a=0; a<rows; a++)
{
for(int b=0; b<cols; b++)
{
X[a][b] = pow(X[a][b],2);
Y[a][b] = pow(Y[a][b],2);
}
}

for(int a=0; a<rows; a++)
{
for(int b=0; b<cols; b++)
{
theta[a][b] = atan2(Y[a][b],X[a][b])*180/PI;
}
}
}

```

With the above code you program now moves on to the next part of the program. However I don't understand what you're trying to do in your logGabor() function. But this seems to be a problem:

```	for(int m=0; m<rows; m++)
{
for(int n=0; n<cols; n++)
{
imfft =  imfft.mul(filter[m][n]);

}
}

```

This computation seems to be taking a very long time, at least on my machine (about 20 min).

By the way what exactly is a Mat?

Jim

### #51 residentelvio

Reputation: 0
• Posts: 30
• Joined: 26-June 13

## Re: LogGabor function problem

Posted 18 July 2013 - 12:29 PM

Hi, ok. I ll try to explain, but the problem is my English is not so good in technical expressions. I m reversing Y because of meshgrid.

```[X,Y] = meshgrid(x,y) transforms the domain specified by vectors x and y into matrices X and Y that can be used for the evaluation of functions of two variables and three-dimensional mesh/surface plots.

The rows of the output matrix X are copies of the vector x      // x contains cols so the rows of X contains cols
the columns of the output matrix Y are copies of the vector y.  // y contains rows so the cols of Y contains rows

```

If you see after in LogGabor, radius and theta will build filter that will use fftshifting and the inverse of Fourier and after I take off the pads

### #52 jimblumberg

Reputation: 4631
• Posts: 14,541
• Joined: 25-December 09

## Re: LogGabor function problem

Posted 18 July 2013 - 12:37 PM

Well that doesn't really explain anything to me, so I'll let someone else try to help.

But I really suggest you look up the documentation for Mat, I suspect you should be using this class instead of all the vectors.

Jim

### #53 residentelvio

Reputation: 0
• Posts: 30
• Joined: 26-June 13

## Re: LogGabor function problem

Posted 29 July 2013 - 02:19 AM

A mat is a Matrix. The problem with Mat is I can not do something like :

```Mat a  ;
for(int x=0;x<10;x++)
a[x]= x;

```

Using a pointer to Mat is it possible?
Howewew now I have an error in the famous function I didn t have before compiling, maybe it s the problem:

```void pre_filter_Computations(vector< vector<double> > &radius,vector< vector<double> > &theta,int cols,int rows){

vector<double>x;
vector<double>y;
vector< vector<double> > X;
vector< vector<double> > Y;
double epsilon = 0.0001;

for(int i=0;i<rows;i++){
for(int k=0;k<cols;k++){
x.push_back(((double)k-(cols/2))/(double)(cols/2));
}
X.push_back(x);
x.clear();
}

for(int z=0;z<cols;z++){
for(int j=0;j<rows;j++){
y.push_back(-(((double)j-(rows/2))/(double)(rows/2)));
}
Y.push_back(y);
y.clear();
}

for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){

X[a][b] = pow(X[a][b],2);
Y[b][a] = pow(Y[b][a],2);
} // no matching function for call to ‘std::vector<std::vector<double> >::push_back(double&)’
}     // no matching function for call to ‘std::vector<std::vector<double> >::push_back(double&)’

for(int a=0;a<rows;a++){
for(int b=0;b<cols;b++){
theta[a][b] = atan2(Y[a][b],X[a][b])*180/PI;
}
}
X.clear();
Y.clear();
}

```