CIMPL Code Examples and Demo

CIMPL Performance Library Home | Baris Sumengen's Home | Discussion Forum

Examples and demonstration of CIMPL

2-D DFT
Singular value decomposition (SVD)
2-D Gabor filters and texture analysis

2-D DFT Example:

The following code first creates a 5x3 complex matrix and fills it with random numbers between 0 and 1. After that, 2-D discrete fourier transform of this matrix is calculated and displayed.

// Create a 5x3 complex matrix initialized with random values between 0 and 1
Matrix<ComplexFloat> X(5,3);
Rand(X, 1.0);
cout << "X: " << X;

Matrix<ComplexFloat> Y = FFT2( X );
cout << "Y: " << Y;

After running this code, the following is displayed on screen:

X: class CIMPL::Matrix<class std::complex<float> > of size 5x3
=======================
ROW 1|    (0.0831019,0.789972)     (0.503098,0.612995)     (0.899411,0.820185)
ROW 2|     (0.431562,0.810633)    (0.0930815,0.341411)   (0.850856,0.00231941)
ROW 3|     (0.696799,0.784936)     (0.121708,0.160802)     (0.103336,0.795618)
ROW 4|     (0.141728,0.683523)     (0.024659,0.849666)     (0.286874,0.415754)
ROW 5|     (0.717795,0.255562)     (0.539201,0.350383)      (0.50264,0.202918)
=======================

Y: class CIMPL::Matrix<class std::complex<float> > of size 5x3
=======================
ROW 1|       (5.99585,7.87667)      (0.176505,2.22758)    (0.0406033,-0.13038)
ROW 2|    (1.54851,-0.0657822)     (0.349803,0.600561)   (-0.856106,-0.902587)
ROW 3|     (-0.225331,2.44667)   (-0.362703,-0.681963)     (-0.485415,1.46728)
ROW 4|       (-1.02634,1.1038)    (-1.92253,-0.143387)     (0.496551,-1.90531)
ROW 5|     (1.13537,-0.245603)      (-2.229,0.0802035)     (-1.38924,0.121821)
=======================

Singular Value Decomposition:

The following code first calculates the singular value decomposition of matrix A. Then matrix A is reconstructed from its singular values and the error between the original and the reconstruction is displayed.

// SVD Calculation

// Create a 3x5 matrix A that will be decomposed
Matrix<double> A(3,5);
// Initialize matrix with random values between 0 and 10
A.Rand(10.0);
cout << "A:\n\n" << A << endl << endl;

// Apply Singular Value Decomposition
Matrix<double> U;
Matrix<double> V;
Vector<double> s = Gesdd(A, U, V);

// Display Singular Value Decomposition
cout << "Singular Values:\n\n" << s << endl;
cout << "U:\n\n" << U << endl;
cout << "V:\n\n" << V << endl;

// Recreate A2 from U, s, and V

// Create diagonal matrix S from s
Matrix<double> S(3,5, 0.0);
S(0,0) = s(0);
S(1,1) = s(1);
S(2,2) = s(2);
cout << "S:\n\n" << S << endl;

// Calculate A2 = U*S*transpose(V) , where S is a 3x5 diagonal matrix
// & operator is overloaded for matrix multiplication
// ~ operator is overloaded as matrix transpose
Matrix<double> A2 = (U & S & ~V);
cout << "A2:\n\n" << A2 << endl;

// Display the calculation error A-A2
cout << "Calculation error A-A2:\n\n" << (A-A2) << endl;

After running this code, the following is displayed on screen:

A:

class CIMPL::Matrix<double> of size 3x5
=======================
ROW 1|   7.99768   6.98538   6.64876    3.1016   3.36711
ROW 2|   1.09317   3.91675   6.42872   6.08509   8.61843
ROW 3|    3.3964   7.44377   9.64568   4.06598   9.74303
=======================



Singular Values:

class CIMPL::Vector<double> of size 3
----------------------
ROW 1|   23.7417 |
ROW 2|   7.03635 |
ROW 3|   2.31459 |
----------------------


U:

class CIMPL::Matrix<double> of size 3x3
=======================
ROW 1|   -0.507154    0.809934    -0.29462
ROW 2|   -0.518296   -0.559738   -0.646578
ROW 3|   -0.688596   -0.175214    0.703659
=======================


V:

class CIMPL::Matrix<double> of size 5x5
=======================
ROW 1|    -0.293214     0.749054    -0.290845   -0.0680344      0.51354
ROW 2|    -0.450619     0.307133     0.279687     -0.49869     -0.61294
ROW 3|     -0.56213    0.0137285      0.29022     0.770742   -0.0745057
ROW 4|    -0.317024    -0.228298    -0.858557    0.0646527    -0.325693
ROW 5|    -0.542655    -0.540628     0.125842    -0.385301     0.498953
=======================


S:

class CIMPL::Matrix<double> of size 3x5
=======================
ROW 1|   23.7417         0         0         0         0
ROW 2|         0   7.03635         0         0         0
ROW 3|         0         0   2.31459         0         0
=======================


A2:

class CIMPL::Matrix<double> of size 3x5
=======================
ROW 1|   7.99768   6.98538   6.64876    3.1016   3.36711
ROW 2|   1.09317   3.91675   6.42872   6.08509   8.61843
ROW 3|    3.3964   7.44377   9.64568   4.06598   9.74303
=======================


Calculation error A-A2:

class CIMPL::Matrix<double> of size 3x5
=======================
----------------------
Columns 1 through 4
----------------------
ROW 1|   -2.66454e-015               0    1.77636e-015    8.88178e-016
ROW 2|   -8.88178e-016   -1.33227e-015    8.88178e-016   -8.88178e-016
ROW 3|   -8.88178e-015   -1.77636e-015               0               0

----------------------
Columns 5 through 5
----------------------
ROW 1|    2.22045e-015
ROW 2|               0
ROW 3|               0
=======================

2-D Gabor filters and texture analysis:

Note that you need to download and install CIMPL performance library first to be able to compile the following code.

The following is a complete program (adapted from the C code written by Wei Ma at Vision Research Lab at UCSB) to
1) Generate complex Gabor filters at certain number of orientations and scales (given the target frequency range),
2) Write them to *.gif files for visualization (after scaling to 0 to 255--might be misleading),
3) Filter an image with these Gabor filters,
4) Write magnitude of the filter outputs to *.gif files (after scaling to 0 to 255--might be misleading).

Image filename is a commandline argument.

Most of the job is handled by the Conv2() function that comes with CIMPL (in Analysis namespace). This function convolves the input image with a Gabor filter. By default symmetric (mirror) boundary conditions are used and the middle part of the filter output (matching the input image) is taken.

#define PI 3.141592653589793115997963468544185161590576171875

// Creates a complex Gabor filter at a certain orientation and scale
Matrix<ComplexFloat> Gabor(int side, int s, int n, float Ul, float Uh, int scale, int orientation, bool flag);

// Filters the input image with Gabor filters at various scales and orientations. 
// Outputs is a list of Complex Matrices, which correspond to the filter outputs.
MatrixList<ComplexFloat> GaborFiltering(char *imagename, int side, float Ul, float Uh, int scale, int orientation, bool flag);

#include "cimpl.h"
using namespace CIMPL;

#include "cimpltoolboxes.h"
using namespace IO;
using namespace MathCore;
using namespace Analysis;

int main(int argc, char* argv[])
{
	int scale = 3;
	int orientation = 4;
	float Ul = 0.1f;
	float Uh = 0.35f;
	int side = 40;
	MatrixList<ComplexFloat> fo = GaborFiltering(argv[1], side, Ul, Uh, scale, orientation, true);

	for (int s=0;s<scale;s++)
	{
		for (int n=0;n<orientation;n++)
		{
			int index = s*orientation+n;
			Matrix<float> filt(fo[index].Rows(), fo[index].Columns());
			for (int k=0;k<filt.Length();k++)
			{
				filt.ElemNC(k) = real(fo[index].ElemNC(k))*real(fo[index].ElemNC(k))+imag(fo[index].ElemNC(k))*imag(fo[index].ElemNC(k));
			}

			char filename[100]; sprintf(filename, "fo.%02d.%02d.gif", s+1, n+1);
			ImWrite(filt, filename);
		}
	}

	return 0;
}


MatrixList<ComplexFloat> GaborFiltering(char *imagename, int side, float Ul, float Uh, int scale, int orientation, bool flag)
{
	MatrixList<unsigned char> image = ImRead(imagename);
	Matrix<float> im;
	ConvertType(image[0],im);

	MatrixList<ComplexFloat> fo;

	for (int s=0;s<scale;s++)
	{
		for (int n=0;n<orientation;n++)
		{
			cout << s+1 << ":" << n+1 << endl;
			Matrix<ComplexFloat> Gb = Gabor(side, s, n, Ul, Uh, scale, orientation, true);
			char filename[100]; 
			sprintf(filename, "gb_real.%02d.%02d.gif", s+1, n+1);
			ImWrite(Real(Gb), filename);
			sprintf(filename, "gb_imag.%02d.%02d.gif", s+1, n+1);
			ImWrite(Imag(Gb), filename);

			Matrix<ComplexFloat> filtered = Conv2( ToComplexFloat(im), Gb);
			
			if(s==0 && n==0)
			{
				fo = MatrixList<ComplexFloat>(filtered);
			}
			else
			{
				fo.PushBack(filtered);
			}
		}
	}

	return fo;
}


Matrix<ComplexFloat> Gabor(int side, int s, int n, float Ul, float Uh, int scale, int orientation, bool flag)
{
	Matrix<float> Gr(2*side+1, 2*side+1);
	Matrix<float> Gi(2*side+1, 2*side+1);
	double base = Uh/Ul;
	double a = pow(base, 1.0/(double)(scale-1));
	double u0 = Uh/pow(a, (double) scale-s);
	double var = pow(0.6/Uh*pow(a, (double) scale-s), 2.0);
	double t1 = cos((double) PI/orientation*(n-1.0));
	double t2 = sin((double) PI/orientation*(n-1.0));

	for (int x=0;x<2*side+1;x++) 
	{
		for (int y=0;y<2*side+1;y++) 
		{
			double X = (double) (x-side)*t1+ (double) (y-side)*t2;
			double Y = (double) -(x-side)*t2+ (double) (y-side)*t1;
			double G = 1.0/(2.0*PI*var)*pow(a, (double) scale-s)*exp(-0.5*(X*X+Y*Y)/var);
			Gr.ElemNC(x,y) = float( G*cos(2.0*PI*u0*X) );
			Gi.ElemNC(x,y) = float( G*sin(2.0*PI*u0*X) );
		}
	}

	/* if flag == true, then remove the DC from the real part of Gabor */
	if (flag == true)
	{
		double m = 0.0;
		for (int x=0;x<2*side+1;x++)
		{
			for (int y=0;y<2*side+1;y++)
			{
				m += Gr.ElemNC(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.ElemNC(x,y) -= (float)m;
			}
		}
	}

	return Complex(Gr, Gi);
}

After running this program, the following images are generated:

Gabor Filters: