2D Vector Class Wrapper SSE Optimized for Math Operations

• Download source - 7.02 KB

Introduction

As you all probably know, arrays are evil. Therefore, I've written a nice 1D/2D vector wrapper with matrix operations support, optimized with SSE for fast processing. It is useful in applications with matrix math where you have convolution, addition, subtraction, multiplication etc... of matrices with floating point precision. For example, I'm using it in neural networks, to support vector machine classifiers for computer vision programs. The class has been written with proper coding style, supporting read-only access, implemented with const-correctness.

Background

You should be familiar with matrices and the math behind them. Have a look at SOS math for a quick start. Have a look at the various SSE programming articles present here in the CodeProject as well.

Using the code

First, make sure you've included the following headers in stdafx.h:

#include <mm3dnow.h>
#include <float.h>
#include <time.h>
#include <math.h>
#include <windows.h>
#include <stdio.h>

The vector wrapper is presented with the vec2D class. You may have an N x M (N by M) 2D vector with N rows and M columns, or a simple 1D vector - either a row vector 1xM or a column vector Nx1.

• vec2D(const wchar_t* file);
• vec2D(unsigned int ysize = 1, unsigned int xsize = 1, int yoffset = 0, int xoffset = 0, const float* data = 0);

With the first constructor, you read the array from the file on the disk in text format:

matrix.txt

2 3
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0

The matrix.txt file provides a 2 by 3 matrix with all zero elements.

Or, you can initialize the vector with the second constructor. If you do not provide any arguments to it, you have a 1 by 1 matrix, which is a scalar. The ysize denotes the number of rows, and xsize, the number of columns. With yoffset and xoffset, you may choose the negative or positive index with which you access the first element in the matrix. Thus, you can build convolution filters like Gaussian for smoothing, or Sobel, Prewitt for edge detection, where you have a 3 by 3 matrix and the zero index element is situated exactly in the center of the matrix. You can also supply an external buffer, data, to initialize the matrix; otherwise, it is filled with zeros.

vec2D v1(1, 10);    //is the 1x10 row vector
vec2D v2(10, 1);    //is the 10x1 column vector

//this one initializes gaus_filter matrix with read-only access
const float gblur[] = { 0.0625f,  0.125f,  0.0625f,
                        0.125f,   0.25f,   0.125f,
                        0.0625f,  0.125f,  0.0625f };
const vec2D gaus_filter(3, 3, -1, -1, gblur);
float val = gaus_filter(0, 0);    //val = 0.25f

With the following functions, you can access the first and the last available indices of the row and the column for the matrix:

• inline int xfirst() const; //return first column index into array
• inline int yfirst() const; //return first row index into array
• inline int xlast() const; //return last col index into array
• inline int ylast() const;//return last row index into array

You can access all the elements of the matrix this way:

//Let's print out our gaus_filter;
for (int y = gaus_filter.yfirst(); y <= gaus_filter.ylast(); y++) {
        for (int x = gaus_filter.xfirst(); x <= gaus_filter.xlast(); x++) 
                wprintf(L" %g", gaus_filter(y, x));        
        wprintf(L"\n");
}

Or you can do the same with:

• void print(const wchar_t* file = 0) const;

If the file is not NULL, you save the matrix to a file on the disk; otherwise, print it to stdout.

To get the sizes of the matrix, use the following functions:

• inline unsigned int width() const; //width of array, number of columns
• inline unsigned int height() const; //height size of array, number of rows
• inline unsigned int length() const; //total size of array width*height

To access the elements of the matrix, there are four overloaded operators present, with both read-write and read-only access:

• inline float& operator()(int y, int x);
• inline float operator()(int y, int x) const;
• inline float& operator[](unsigned int i);
• inline float operator[](unsigned int i) const;

With the [] operators, you can access all the elements in a Matlab like fashion. But, the first row and column indices in your matrix should start with 0.

//In Matlab we define the array like:
// v = zeros(3, 5);
//and we can print its contents with [] operator
// v[:]

//The same applies to [] overloaded operator
vec2D v(3, 5);   //we have xfirst() and yfirst() equals to 0 as ordinary C array
for (unsigned int i = 0; i < v.length(); i++)
        wprintf(L"%g\n", v[i]);

You can access the elements of the matrix also with:

• inline float get(int y, int x) const;

It allows the indices to exceed the matrix dimensions. In that case, you will have a periodic boundary extension. So, if you have your first row and column indices equal to 0, and try to access vec(-1, -1), you get the vec(1, 1) element.

With overloaded assignment operators, you can make a copy of the matrix, or initialize it with an external data buffer:

• inline const vec2D& operator=(const vec2D& v);
• inline const vec2D& operator=(const float* pf);

const float gblur[] = { 0.0625f,  0.125f,  0.0625f,
                        0.125f,   0.25f,   0.125f,
                        0.0625f,  0.125f,  0.0625f };
vec2D v1(3, 3);
v1 = gblur;
vec2D v2(3, 3);
vec2D v3(10, 5);
v2 = v1;
v3 = v2;    //v3 will be resized to 3 by 3 matrix

With set element functions, the matrix contents can be initialized to specific values:

• void set(float scalar); //set all matrix elements to scalar
• void set(float scalar, RECT& r); //set matrix elements in the RECT to scalar
• void setrand(); //set matrix elements to random values in the range -1.0 ... 1.0

The RECT is just a Windows structure.

SSE optimized operations

The following functions deal with SSE optimized math operations:

• void add(float scalar);
• void sub(float scalar);
• void mul(float scalar);
• void div(float scalar);

You can add, subtract, multiply, or divide all the matrix elements with a scalar.

vec2D v(3, 5);
v.print();
v.add(3.4f);
v.print();

The functions for arithmetic operations on matrices are:

• void add(const vec2D& a, const vec2D& b); //c = a.+b sse optimized
• void sub(const vec2D& a, const vec2D& b); //c = a.-b sse optimized
• void mul(const vec2D& a, const vec2D& b); //c = a*b
• void div(const vec2D& a, const vec2D& b); //c = a./b sse optimized

Where c is the matrix whose function you are calling, and its elements will be the result of the operation. add, sub, and div are the element wise operations, so the a and b matrices along with c should be of the same size. The mul (multiplication) operation is not SSE optimized, as during multiplication, every row in a is multiplied by every column in b.

vec2D a(10, 10);
vec2D b = a;
vec2D c = a;
a.setrand();
b.setrand();
c.add(a, b);      //c = a + b
a.print();
b.print();
c.print();

vec2D a(3, 5);
vec2D b(5, 2);
vec2D c(3, 2);
a.setrand();
b.setrand();
c.mul(a, b);     //c = a * b

The SSE optimized matrices multiplication goes here:

• void mult(const vec2D& a, const vec2D& b); //c = a*b' sse optimized
• void mule(const vec2D& a, const vec2D& b); //c = a.*b sse optimized

The difference between mult and mul is that the b matrix is transposed, so during multiplication, we multiply every row from a by every row from b, which is easily SSE optimized.

//arrange transposed version of b for mult
vec2D a(3, 5);
vec2D b(2, 5);
vec2D c(3, 2);
a.setrand();
b.setrand();
c.mult(a, b);     //c = a * b';   SSE optimized
a.print();
b.print();
c.print();

And, mule is just element wise multiplication of equally sized a and b matrices.

2D convolution, normalization, and histogram equalization

The other function useful in signal processing is convolution with a filter:

• void conv2D(const vec2D& a, const vec2D& filter);
• void conv2D(const vec2D& a, const vec2D& re, const vec2D& im);

The first one is just convolution with a filter, and the second is convolution with a complex filter consisting of real and imaginary parts.

//gaussian smoothing
vec2D v_original(100, 100);
vec2D v_smoothed = v_original;
v_smoothed.conv2D(v_original, gaus_filter);


//edge detection
const float Re[] = { 0.5f, 0.0f, -0.5f,
                     0.5f, 0.0f, -0.5f,
                     0.5f, 0.0f, -0.5f };
const float Im[] = { 0.5f,  0.5f,  0.5f,
                     0.0f,  0.0f,  0.0f,
                    -0.5f, -0.5f, -0.5f};
vec2D re(3, 3, -1, -1, Re);
vec2D im(3, 3, -1, -1, Im);
vec2D v_edges(100, 100);

//set the rectangle in the matrix to 0.0 
//so there will be edges on the brim of it
RECT r(10, 10, 50, 50);
v_smoothed.set(0.0f, r);

v_edges.conv2D(v_smoothed, re, im);

The normalization and histogram equalization are used in image processing, for example, during a face detection process:

• void normalize(float a, float b); //normalize to a...b range
• void histeq(vec2D& hist); //histogram normalization

For normalization, the matrix is just normalized so its minimum value becomes a and maximum value is b.

For histogram equalization, provide an external row vector hist of size 1 by 256. The matrix before equalization should contain values in the range of 0 ... 255 as the normal bitmap image. After equalization, its values are mapped to the range of 0 ... 1.0.

vec2D hist(1, 256);
vec2D v_heq(100, 100);
v_heq.setrand();
v_heq.norm(0.0f, 255.0f);
v_heq.histeq(hist);

You must Sign In to use this message board.

FAQ  Noise Tolerance Very HighHighMediumLowVery Low  Search Messages    Layout Expand All MessagesThread ViewNo Javascript (slow)No Javascript Preview    Per page 102550

Msgs 1 to 11 of 11 (Total in Forum: 11) (Refresh) FirstPrevNext

i dun understand why u need an offset in init function f2 3:04 16 Oct '08   init(unsigned int ysize, unsigned int xsize, int yoffset, int xoffset) pls advice from, -= aLbert =- Sign In·View Thread·PermaLink Re: i dun understand why u need an offset in init function Chesnokov Yuriy 21:53 16 Oct '08   To have different indexed arrays. In Matlab you access arrays from 1 to N indexes as in C++ from 0 to N - 1. There are some functions that needs to access them say from -5, -4, ... 3 so at 0 index you have a center of that function chesnokov Sign In·View Thread·PermaLink Re: i dun understand why u need an offset in init function f2 8:28 17 Oct '08   thx for ur real fast response. so u are currently working in Confidential Japanese firm. that's interesting. should be challenging. do u explore into SIFT algo b4, if yes, how do u rate his method compare to ur face detection method. from, -= aLbert =- Sign In·View Thread·PermaLink Re: i dun understand why u need an offset in init function Chesnokov Yuriy 21:47 17 Oct '08   Yes confidential, so other competitive companies will not divert me or something. So I can only comment on my current articles at codeproject. Hence I do not know if I explore into SIFT 'May be yes, may be no' - Little Britan, Scotish hotel sketch. chesnokov Sign In·View Thread·PermaLink Re: i dun understand why u need an offset in init function f2 22:58 17 Oct '08   u can try, it gives some new approach of getting features point. SIFT by David Lowe. Scale Invariant Feature Transform. but i think ur Haar wavelet implementation, it was claim better. from, -= aLbert =- modified on Saturday, October 18, 2008 5:22 AM Sign In·View Thread·PermaLink Bad Kuryn 8:38 30 May '08   I found your code badly designed. I maybe a year or so late. But most of your code. Like mul, div. could have been put into operator overloading and made into a friend class. And you could have also put it into templates also. Iono if your trying to keep C compatiability or what. But rewrite it to something better. Sign In·View Thread·PermaLink Re: Not Bad but Great! Chesnokov Yuriy 4:05 1 Jun '08   I found you blind looking at my code. There is a mul in overloaded operator. Taking just once a mul or div worth operator invokation, but consider taking 1000000000 of muls (convolution with a filter), preallocation is a much faster rather than allocating every time you do the mul. The code is specialized one for floats only with SSE support. I'm not going it to use in C also. Providing templates will be a bunch of work to support every operation in MMX, SSE, SSE2. You can rewrite it yourself for templates supporting every MMX, SSE, SSE2 so I will be able to put your modes here chesnokov Sign In·View Thread·PermaLink A question Ali Khanlarkhani 5:58 31 Oct '07   Can I simply change "float" to "double"? without loosing performance Sign In·View Thread·PermaLink 4.00/5 (2 votes) Re: A question Chesnokov Yuriy 8:12 31 Oct '07   No you can't. The operations are SSE optimized for floats. It will fail if just put doubles, you have to change it to SSE2. Have a look at my 'convolution experiment' article for speed comparison for floats and doubles chesnokov Sign In·View Thread·PermaLink No SSE command was found tyjiang 22:18 16 Oct '07   No SSE command was found in your program. Only a few ms intrinsic command was used. Sign In·View Thread·PermaLink Re: No SSE command was found Chesnokov Yuriy 1:54 17 Oct '07   Yes, they are provided by microsoft so you do not have to delve into SSE assembly, the code becomes more concise. It does not need a lot of intrinsics, as the optimization goes with addition, division, multiplication, convolution and these are the main operations you are involved in neural networks, support vector machines etc... chesnokov Sign In·View Thread·PermaLink 3.50/5 (2 votes)

Last Visit: Monday, December 1, 2008 8:57 PM     Last Update:Monday, December 1, 2008 8:57 PM