//***************************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "image.h"
#include "matrix.h"
#include "gabor.h"
//***************************************************************************
// The Gabor routine is modified from the following code:
// http://vivaldi.ece.ucsb.edu/users/wei/code_gabor/
//***************************************************************************
// The Gabor function generates a Gabor filter with the selected index 's'
// and 'n' (scale and orientation, respectively) from a Gabor filter bank.
// This filter bank is designed by giving the range of spatial frequency
// (Uh and Ul) and the total number of scales and orientations used to
// partition the spectrum.
//
// The returned filter is stored in 'Gr' (real part) and 'Gi' (imaginary part).
//
// Symbol translation from paper (pg 19, 20 of extended paper):
// http://vivaldi.ece.ucsb.edu/users/wei/mypapers/edgeflow.ps.gz
// S = scale
// k = orientation
//***************************************************************************
// The GaborFilteredImg provides the outputs of the Gabor filter bank
void GaborFilteredImg(Matrix **OutR, Matrix **OutI,
Matrix &img, int side, double Ul, double Uh, int scale, int orientation, int flag)
{
int border = side;
int height = img.nx;
int width = img.ny;
/* FFT2 */
int xs = (int) pow(2.0, ceil(log2((double)(height+2.0*border))));
int ys = (int) pow(2.0, ceil(log2((double)(width+2.0*border))));
Matrix IMG(xs,ys);
int r1, r2, r3, r4;
r1 = width+border;
r2 = width+2*border;
for (int w = 0; w < border; w++)
{
for (int h = 0; h < border; h++)
IMG.p[w][h] = img.p[border-1-w][border-1-h];
for (int h = border; h < r1; h++)
IMG.p[w][h] = img.p[border-1-w][h-border];
for (int h = r1; h < r2; h++)
IMG.p[w][h] = img.p[border-1-w][2*width-h+border-1];
}
r1 = height+border;
r2 = width+border;
r3 = width+2*border;
for (int w = border; w < r1; w++)
{
for (int h = 0; h < border; h++)
IMG.p[w][h] = img.p[w-border][border-1-h];
for (int h = border; h < r2; h++)
IMG.p[w][h] = img.p[w-border][h-border];
for (int h = r2; h < r3; h++)
IMG.p[w][h] = img.p[w-border][2*width-h+border-1];
}
r1 = height+border;
r2 = height+2*border;
r3 = width+border;
r4 = width+2*border;
for (int w = r1; w < r2; w++)
{
for (int h = 0; h < border; h++)
IMG.p[w][h] = img.p[2*height-w+border-1][border-1-h];
for (int h = border; h < r3; h++)
IMG.p[w][h] = img.p[2*height-w+border-1][h-border];
for (int h = r3; h < r4; h++)
IMG.p[w][h] = img.p[2*height-w+border-1][2*width-h+border-1];
}
Matrix F_real(xs,ys), F_imag(xs,ys), IMG_imag(xs,ys);
MatrixFFT2D(&F_real, &F_imag, &IMG, &IMG_imag);
// ----------- compute the Gabor filtered output -------------
Matrix Gr(2*side+1, 2*side+1), Gi(2*side+1, 2*side+1);
Matrix Tmp_1(xs,ys), Tmp_2(xs,ys);
Matrix F_1(xs,ys), F_2(xs,ys);
Matrix G_real(xs,ys), G_imag(xs,ys);
cout << " GaborFilteredImg (S=" << scale << ",N=" << orientation
<< "): Processing Image 0/" << scale*orientation;
for (int s = 0; s < scale; s++)
{
for (int n = 0; n < orientation; n++)
{
cout << "\r GaborFilteredImg (S=" << scale << ",N=" << orientation
<< "): Processing Image " << s*orientation+n+1 << "/" << scale*orientation;
cout.flush();
Gabor(&Gr, &Gi, s+1, n+1, Ul, Uh, scale, orientation, flag);
MatCopy(&F_1, &Gr, 0, 0, 0, 0, 2*side, 2*side);
MatCopy(&F_2, &Gi, 0, 0, 0, 0, 2*side, 2*side);
MatrixFFT2D(&G_real, &G_imag, &F_1, &F_2);
// ^ stands for pixel by pixel multiplication
Tmp_1 = G_real ^ F_real;
Tmp_2 = G_imag ^ F_imag;
IMG = Tmp_1 - Tmp_2;
// ^ stands for pixel by pixel multiplication
Tmp_1 = G_real ^ F_imag;
Tmp_2 = G_imag ^ F_real;
IMG_imag = Tmp_1 + Tmp_2;
MatrixIFFT2D(&Tmp_1, &Tmp_2, &IMG, &IMG_imag);
MatCopy(&OutR[s][n], &Tmp_1, 0, 0, 2*border, 2*border, height+2*border-1, width+2*border-1);
MatCopy(&OutI[s][n], &Tmp_2, 0, 0, 2*border, 2*border, height+2*border-1, width+2*border-1);
}
}
cout << "...done!" << endl;
cout.flush();
}
//***************************************************************************
void Gabor(Matrix *Gr, Matrix *Gi, int s, int n, double Ul, double Uh, int scale, int orientation, int flag)
{
double base = Uh/Ul;
double a = pow(base, 1.0/(double)(scale-1));
double u0 = Uh/pow(a, (double) (scale-s));
double Uvar = (a-1.0)*u0/((a+1.0)*sqrt(2.0*log(2.0)));
double z = -2.0*log(2.0)*(Uvar*Uvar)/u0;
double Vvar = tan(Pi/(2*orientation))*(u0+z)/sqrt(2.0*log(2.0)-z*z/(Uvar*Uvar));
double Xvar = 1.0/(2.0*Pi*Uvar);
double Yvar = 1.0/(2.0*Pi*Vvar);
double t1 = cos(Pi/orientation*(n-1.0));
double t2 = sin(Pi/orientation*(n-1.0));
int side = (Gr->nx-1)/2;
for (int x = 0; x < 2*side+1; x++)
{
for (int y = 0; y < 2*side+1; y++)
{
double X = (double) ( (x-side)*t1 + (y-side)*t2);
double Y = (double) (-(x-side)*t2 + (y-side)*t1);
double G = 1.0/(2.0*Pi*Xvar*Yvar)*pow(a, (double) (scale-s))
*exp(-0.5*((X*X)/(Xvar*Xvar)+(Y*Y)/(Yvar*Yvar)));
Gr->p[x][y] = G*cos(2.0*Pi*u0*X);
Gi->p[x][y] = G*sin(2.0*Pi*u0*X);
}
}
/* if flag = 1, then remove the DC from the filter */
if (flag == 1)
{
double m = 0.0;
for (int x = 0; x < 2*side+1; x++)
for (int y = 0; y < 2*side+1; y++)
m += Gr->p[x][y];
m /= pow((double) (2.0*side+1), 2.0);
for (int x = 0; x < 2*side+1; x++)
for (int y = 0; y < 2*side+1; y++)
Gr->p[x][y] -= m;
}
}
//***************************************************************************