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); vec2D v2(10, 1);
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);
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:
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.
vec2D v(3, 5); 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;
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); 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);
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.
vec2D a(3, 5);
vec2D b(2, 5);
vec2D c(3, 2);
a.setrand();
b.setrand();
c.mult(a, b); 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.
vec2D v_original(100, 100);
vec2D v_smoothed = v_original;
v_smoothed.conv2D(v_original, gaus_filter);
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);
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);