Help, cannot identify where problem is!

[darelet]

Hello,
I'm researching in image processing and need to develop a Gabor filter in C++ and openCV. I found a good sample on which I'm basing my code. Problem is, while I've tried to emulate much as possible the other code, mine is not producing correct results. I changed a few things to suit my context but I don't see why it is not working! Please help, because I've labored on this for this last four days! I'll paste the three functions that are of concern, as well as the code I'm emulating below mine. Thanks in advance..!

(the following is the problematic code)

unsigned char *ideas::gabor(dataStructure *data)
{
	int height=data->nRows, width=data->nCols;
	float lambda=5.0; float bandwidth=01.0;
	//float PI = (float) System::Math::PI;
	#define PI 3.14159265358979323846
	float sigma = (float) (1.0/PI * sqrt(log(2.0)/2.0) *
		 (pow(2.0f,bandwidth) + 1.0)/(pow(2.0f,bandwidth) - 1.0) * lambda);
	
	//Gabor filter
	int filterHalfWidth = (int) (4.0*sigma + 0.5); // half width of gaussian enveloppe
	int filterWidth = 2*filterHalfWidth + 1;
	float norm = 0.0f;
	float *gx_r=new float[filterWidth], *gx_i=new float[filterWidth], 
		*gy=new float[filterWidth];// all of size [filterWidth]
	initializeArray(gx_r,filterWidth);	initializeArray(gx_i,filterWidth);	initializeArray(gy,filterWidth);
	
	for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i) { 
		float g = (float) exp( -(i*i) / (2.0*sigma*sigma) );
		gx_r[filterHalfWidth+i] = (float) (g * cos(2.0*PI/lambda * i + 0.0));
		gx_i[filterHalfWidth+i] = (float) (g * cos(2.0*PI/lambda * i + PI/2.0));
		gy[filterHalfWidth+i] = g;	norm += g;
	}

	float toto=0.0f;
	for (int y = 0; y < filterWidth; ++y){
		gx_r[y] /= norm; 
		gx_i[y] /= norm;	
		gy[y] /= norm;
		//printf("\n gx_r: %3.5f  gx_i: %3.5f  gy:%3.5f ",gx_r[y],gx_i[y],gy[y]);
		toto+=(gx_r[y]+gx_i[y]+gy[y]);
	}

	printf("\n filterHalfWidth:%3.3f, filterWidth:%3.5f",filterHalfWidth, filterWidth);
	printf("\n toto:%3.3f, sigma:%3.5f, norm:%3.5f, lambda:%3.5f \n",toto, sigma, norm, lambda);
	float *real=new float[width*height];		real=convolve2D(data,gx_r,gy,filterHalfWidth);
	float *imaginary=new float[width*height];	imaginary=convolve2D(data,gx_i,gy,filterHalfWidth); 
	float *finalImage=new float[width*height];	initializeArray(finalImage, width*height);
	
	for (int i=0; i<width*height;i++) printf("  %f",real[i]);

	////////////////////////////////////////
	// Combine real and imaginary component ( out = r^2 + i^2 )
	////////////////////////////////////////
	
	float max = -1.0;
	for (int y = filterHalfWidth; y < height-filterHalfWidth; ++y) {
		for (int x = filterHalfWidth; x < width-filterHalfWidth; ++x) {
			finalImage[y*width+x]  = real[y*width+x]*real[y*width+x];
			//finalImage[y*width+x] += imaginary[y*width+x]*imaginary[y*width+x];
			if (finalImage[y*width+x] > max) max = finalImage[y*width+x];
		}
		//printf(" %f",real[y*width+x]);
	}
	unsigned char* out=new unsigned char[width*height];
	initializeArray(out, width*height);

	for (int y = 0; y < height; ++y)for (int x = 0; x < width; ++x) {
			finalImage[y*width+x] = 255.0/max;
			out[y*width+x] = (uchar) finalImage[y*width+x];
	}
return out;
}
unsigned int var1=0;
float *ideas::convolve2D(dataStructure *data, const float* fX, const float* fY, int filterHalfWidth) {
	int height=data->nRows, width=data->nCols;
	double *tmp = new double[height*width]; initializeArray(tmp, height*width); double ttol=0;
	for (int y = 0; y < height; ++y){
		for (int x = filterHalfWidth; x < width-filterHalfWidth; ++x){
			for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i){
				tmp[y*width+x] += ((int) data->data[y*width+(x-i)]) * fX[i+filterHalfWidth];
				ttol+= (int) data->data[y*width+(x-i)] * fX[i+filterHalfWidth];
			}
			//printf(" %f",tmp[y*width+x]);
			//ttol=0;
		}
		//printf("\n");
	}


	float *output=new float[height*width]; initializeArray(tmp, height*width);
	for (int y = filterHalfWidth; y < height-filterHalfWidth; ++y)
		for (int x = filterHalfWidth; x < width-filterHalfWidth; ++x)
			for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i)
				output[y*width+x] += tmp[(y-i)*width+x] * fY[i+filterHalfWidth];
	return output;
}

template <class t>
void ideas::initializeArray(t *arr, int size){
	//arr = new t [size]; 
	int i=0;
	while (i<size) {
		arr[i]= (t) 0; ++i;
	}
}

(The following is what I'm trying to emulate)

#define M_PI 3.14159265358979323846

CImg<float> convolve2D(const CImg<float> &input, const CImg<float> &fX, const CImg<float> &fY);


int main(int argc, char **argv) {
	////////////////////////////////////////
	// read command line arguments
	////////////////////////////////////////
	string infile = "c:\\tit\\lena.bmp";
	string outfile = "c:\\tit\\savedimage.jpg";
	
	
	float lambda = 5.0;
	float bandwidth = 1.0;
	
	/*argstream as(argc, argv);
	as >> parameter("i", infile,    "infile ", true);
	as >> parameter("o", outfile,   "outfile ", false);
	as >> parameter("l", lambda,    "wavelength lambda [pixel]", false);
	as >> parameter("b", bandwidth, "bandwidth [octaves]", false);
	as >> help();
	as.defaultErrorHandling();*/
	
	
	////////////////////////////////////////
	// load input image
	////////////////////////////////////////
	CImg<float> imgInput(infile.c_str());
	const int width = imgInput.dimx();
	const int height = imgInput.dimy();
	if ( width == 0 || height == 0 ) {
		cerr << "Error when loading input image." << endl;
		return -1;
	}
	imgInput = imgInput.get_norm_pointwise(1)/(float)imgInput.dimv(); // to grey
	
	cout << "Resolution: "<< width << " x " << height << endl;
	
	
	////////////////////////////////////////
	// calculate spread of spatial gaussian
	// based on bandwidth
	////////////////////////////////////////
	float sigma = 1.0/M_PI * sqrt(log(2.0)/2.0) *
		(pow(2.0f,bandwidth) + 1.0)/(pow(2.0f,bandwidth) - 1.0) * lambda;
	
	
	
	////////////////////////////////////////
	// Create separated Gabor filter for real and imaginary component
	////////////////////////////////////////
	int filterHalfWidth = (int) (4.0*sigma + 0.5); // half width of gaussian enveloppe
	int filterWidth = 2*filterHalfWidth + 1;
	float norm = 0.0f;
	CImg<float> gx_r(filterWidth), gx_i(filterWidth), gy(filterWidth);
	for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i) { 
		float g = exp( -(i*i) / (2.0*sigma*sigma) );
		gx_r[filterHalfWidth+i] = g * cos(2.0*M_PI/lambda * i + 0.0);
		gx_i[filterHalfWidth+i] = g * cos(2.0*M_PI/lambda * i + M_PI/2.0);
		gy[filterHalfWidth+i] = g;
		norm += g;
	}
	// normalize the gaussian
		for (int y = 0; y < filterWidth; ++y){
		//printf("\n gx_r: %f  gx_i: %f  gy:%f ",gx_r[y],gx_i[y],gy[y]);
	}
	gx_r /= norm; gx_i /= norm; gy /= norm;
	
	//printf("\n filterHalfWidth:%3.3f, filterWidth:%3.5f",filterHalfWidth, filterWidth);
	printf("\n norm: %3.5f lambda: %3.5f sigma: %3.5f",norm, lambda,sigma);
	
	////////////////////////////////////////
	// 1. Convolve input image with gx_r and gy --> r_component
	// 2. Convolve input image with gx_i and gy --> i_component
	////////////////////////////////////////
	CImg<float> r_component = convolve2D(imgInput, gx_r, gy);
	CImg<float> i_component = convolve2D(imgInput, gx_i, gy);
	
	for (int y = 0; y < height; ++y){for (int x = 0; x < width; ++x) {
		printf(" r: %f, \t c: %f", r_component,i_component );
	}printf("\n");}
	////////////////////////////////////////
	// Combine real and imaginary component ( out = r^2 + i^2 )
	////////////////////////////////////////
	CImg<float> imgOut(width, height, 1, 1); imgOut.fill(0.0f);
	float max = -1.0;
	for (int y = filterHalfWidth; y < height-filterHalfWidth; ++y) {
		for (int x = filterHalfWidth; x < width-filterHalfWidth; ++x) {
			imgOut(x,y)  = r_component(x,y)*r_component(x,y);
			imgOut(x,y) += i_component(x,y)*i_component(x,y);
			if (imgOut(x,y) > max) max = imgOut(x,y);
		}
	}
	
	
	////////////////////////////////////////
	// display the filtered image
	////////////////////////////////////////
	imgOut *= 255.0/max; // scale output range to [0;255]
	CImgl<float> (imgInput, imgOut).display ("Input Image and Normalized Gabor response",'y');
	
	if (outfile.size() > 0) imgOut.save(outfile.c_str());
	return 0;
}


CImg<float> convolve2D(const CImg<float> &input, const CImg<float> &fX, const CImg<float> &fY) {
	const int width  = input.dimx();
	const int height = input.dimy();
	int filterHalfWidth = (fX.dimx()-1) / 2;

	CImg<float> tmp(width, height, 1, 1); tmp.fill(0.0f); float ttol=0;
	for (int y = 0; y < height; ++y){
		for (int x = filterHalfWidth; x < width-filterHalfWidth; ++x){
			for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i){
				tmp(x,y) += input(x-i,y) * fX[i+filterHalfWidth];
				ttol+= input(x-i,y)  * fX[i+filterHalfWidth];
			}
			//printf(" %f",ttol);
			//ttol=0;
		} 
		//printf("\n");
	}

	
	CImg<float> output(width, height, 1, 1); output.fill(0.0f);
	for (int y = filterHalfWidth; y < height-filterHalfWidth; ++y)
		for (int x = filterHalfWidth; x < width-filterHalfWidth; ++x)
			for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i)
				output(x,y) += tmp(x,y-i) * fY[i+filterHalfWidth];
	
	return output;
}

[daviddoria]

I think you'd have a lot better luck on the OpenCV forum.

Dave

[darelet]

I think you'd have a lot better luck on the OpenCV forum.

Dave

Thanks for the reply. I think the Problem is more of C++ than it is an OpenCV one. Furthermore, I only use OpenCV for reading images (into my matrices) and saving back into the files, so the dependency is only at those two levels - beginning and end. Anyway, I found a work-around on it. Someone please advise me where I can get help on the Gabor filter, from an implementation point of view!

[Nick Evan]

From wikipedia perhaps?

Next time you ask a question about a code, you should probably mention what is wrong with it.
- Does it compile? If not: post error messages
- Does it produce expected output? If no: what is the output you were expecting, and what is your current output.
That way we can be of more help.

[Aranarth]

A glaring mistake I noticed when skimming the code is

float *real=new float[width*height]; real=convolve2D(data,gx_r,gy,filterHalfWidth);
float *imaginary=new float[width*height]; imaginary=convolve2D(data,gx_i,gy,filterHalfWidth);

You're requesting memory for real and imaginary, then immediately overwrite the pointer, creating a huge memory leak.
And while I see an abundance of new[]s in other parts of your code, I don't see a single delete[]. You shouldn't do manual memory management unless you know exactly what you're doing. Instead, use classes like std::vector and boost::scoped_array.

[mitrmkar]

I think the Problem is more of C++ than it is an OpenCV one.

Looking at the code you've posted, you have huge problems with memory overruns, and forgetting to delete the allocated memory (which may be hiding these overruns).

One option would be to use vector<float> for all of your arrays. That way, the memory is managed automatically, plus you'll be alerted to out-of-bounds read/writes at the moment they occur.

Try the following program and see what happens

#include <iostream>
#include <vector>

int main()
{
  using namespace std;

  int filterHalfWidth = 3;
  vector<float> gx_v(filterHalfWidth, 0.0f);

  try 
  {
    // Note:
    // This is exactly the same loop construct that you are using
    //
    for (int i = -filterHalfWidth; i <= filterHalfWidth; ++i)
    {
      cout << "accessing index: " << filterHalfWidth+i;

      // Use .at() rather than [], because it gives you
      // an out_of_range C++ exception on illegal access ..
      gx_v.at(filterHalfWidth+i) = i;

      // Doing good
      cout << "-> OK" << endl;
    }
  }
  catch(out_of_range & e)
  {
    // .. which you can then catch
    cout << "-> FAIL" << endl << "exception: " << e.what() << endl;
  }

  // Done, the memory reserved by the vector will be 
  // automatically released in a moment or two ..

  return 0;
}

[darelet]

Mit, Nick, Aran and David <-- Thanx for your help. While I had "patched" my way thro' as far as the problem is concerned, I think I need to follow your advice if I'm to be a disciplined programmer (and forum member, Nick ;-)). I was very aware of some of the bad practices, especially not deallocating memory, but I normally work by getting my code right then "cleaning" up later!
@Mitrmkar: I've been avoiding vectors (and most libraries) all along since I'll need to parallelize my code later, once the serial version is running; so apart from knowing the nitty gritties of all functions I use for optimal solution design, I also need to own them ...