1. Preface
As metioned earliear, in this post I present briefly how to combine Open Computer Vision Library and Blitz++ in a easy way. The purpose of this combination is to take advantage of elegant matrix operations on Blitz++ (almost like Matlab language) and powerful implementation of computer vision and image filtering algorithms.
2. Motivation
Our object of interest is the data arrays stored in memory and referred as C++ pointers or references. Say, we have a 2D matrix with 
In Matlab
subwnd= A(10: (10+3), 12: (10+3));
In Blitz++
Array<int,2> subwnd = A(Range(10:13), Range(12:15)); // or TinyVector<int,2> lbound(10,12); TinyVector<int,2> ubound(13,15); RectDomain<2> rect(lbound,ubound); Array<int,2> subwnd = A(rect);
Another example is to extract the 3rd row in a matrix:
In Matlab
row = A(3,:);
In Blitz++
Array<int,1> row = A(3, Range::all());
In summary, if you work much with image and video, those operations are neccessary to improve coding speed, error-prone, and readability. Another linear algebra library also attract much my attention is Armadillo. Basically, Armadillo introduces more compact code syntax than Blitz++. We will talk about this library in the future.
3. The core idea
Let’s think OpenCV coding syntax and Blitz++ coding syntax are just programming interfaces. We, in fact, manipulate the underlying memory. The memory layout in C or C++ is identical no matter what library you use. Hence, when you would like to use Blitz++ operands, let’s point Array<T,N> to the data. And if you like OpenCV to takes part in, let’s point Mat header to the data.
4. A case study
/**
*@brief Gaussian smoothing for 2D array
*
*@param A[in] 2D array want to smooth
*@param C[out] 2D array out, must be allocated before
*@param sigmas[in] standard deviations along each dimension
*@param radius[in] mask size = 2*radius + 1
*
*@see Convolve3D, FilterGauss1D
*/
int GaussianSmooth2D (
Array<float,2>& src,
Array<float,2>& dst,
TinyVector<float,2>& sigmas,
TinyVector<float,2>& radius
)
{
using namespace cv;
// Compute radii for the Gaussian kernel
TinyVector<float,2> r = ceil (sigmas * radius);
Mat blurred;
// get the kernels of Gaussian
Mat kernelX = getGaussianKernel (r(0), sigmas(0));
Mat kernelY = getGaussianKernel (r(1), sigmas(1));
// create a OpenCV Mat header and point to the data owned by 'src' (Blitz++'s type)
Mat fsrc(src.rows(), src.cols(), CV_32FC1, src.data());
// take 2D filtering on that data
sepFilter2D (fsrc, blurred, CV_32FC1, kernelX, kernelY);
// allocate new data space for output
dst.resize(src.shape());
// create a Blitz++ Array<T,N> header and point to the data owned by 'blurred' (OpenCV's type)
Array<float,2> tmp( (float*)(blurred.data), // this is the pointer
shape(blurred.rows, // the number of rows
blurred.cols), // the number of cols
neverDeleteData); // just refer to it, do not release it
// copy data from 'tmp' to 'dst'
dst = tmp;
return 1;
// At the end of this function, data owned by 'src' still remain and unchanged
// because 'fscr' is just a Mat header, hence it does not release the data.
// Besides, both 'tmp' and 'blurred' point to the same data space. We might fall
// into "double" delete. But we defend by indicating 'neverDeleteData' in 'tmp'.
// Consequently, just 'blurred' release its data.
}
5. Notes
5.1 Memory allocation
Array memory allocation in C/C++ is continuous in its physical address. In other word, the next element in an array is also be the next element in memory. Therefore, one can use athrimethic operators such as 
5.2 Memory layouts
When using useful operators of Blitz++, one must take attention in the continuity property of the underlying data.
// following allocations are continuity in memory layout Array<int,1> arr1(10); Array<int,2> arr2(20,40); Array<int,3> arr3(10,20,30);
Array<float,2> matrix(10,10); matrix = 1; // the 'subm' is not continuous in memory layout Array<float,2> subm = matrix(Range(2,4), Range(3,5)); // the 'row' is continuous in memory layout Array<float,1> row = matrix(4,Range::all()); // the 'col' is not continuous in memory layout Array<float,1> col = matrix(Range::all(), 6); // the 'stride' is not continuous in memory layout Array<float,2> stride = matrix(Range(1,7,3), Range(1,5,2));
Assume that you might apply operations on such sub matrices that requires continuous property (i.e. block processing), how to overcome this limitation? Just copy them to another newly allocated array with the same size. After then, the new array has continuous memory layout and can be referred by OpenCV Mat structure.
The following program gets input a 3D volume and outputs a slice of interest. A slice can be in the xy plane, xt plane, or yt plane. The output must be contiguous memory layout to be processed by OpenCV block processing procedures.
/**
* @brief Extract a 2D slice along arbitrary dimension (x, y, or t) of a volume
*
* @param vol[in] a 3D array of Blitz++, says spatio-temporal space
* @param slice[out] a 2D slice indicated by idx in the dim_th dimension
* @param idx[in] index value indicate which slice needed to be extracted
* @param dim[in] the dimension at which we want to extract
* @param show[in] show the extracted slice onto HighGUI or not
*/
template<typename _SrcTp, typename _DstTp>
int slice2mat (Array<_SrcTp,3>& vol, Mat& slice, int idx, int dim, int show)
{
if (dim == firstDim) // get slice along y-t plane
{
// the 'vol' variable has dis-contiguous memory layout
Array<_SrcTp,2> tmp = vol( idx, blitz::Range::all(), blitz::Range::all());
// copy tmp's content to dump so that dump has contiguous memory layout
Array<_SrcTp,2> dump = tmp.copy ();
// create a Mat header and point to the dump's data
Mat mat_tmp(dump.rows(), dump.cols(), DataType<_SrcTp>::type, dump.data());
// apply arbitrary OpenCV operation
mat_tmp.convertTo (slice, DataType<_DstTp>::type);
}
else if (dim == secondDim) // get slice along x-t plane
{
Array<_SrcTp,2> tmp = vol( blitz::Range::all(), idx, blitz::Range::all());
Array<_SrcTp,2> dump = tmp.copy ();
Mat mat_tmp(dump.rows(), dump.cols(), DataType<_SrcTp>::type, dump.data());
mat_tmp.convertTo (slice, DataType<_DstTp>::type);
}
else if (dim == thirdDim) // get slice along x-y plane
{
Array<_SrcTp,2> tmp = vol( blitz::Range::all(), blitz::Range::all(), idx);
Array<_SrcTp,2> dump = tmp.copy ();
Mat mat_tmp(dump.rows(), dump.cols(), DataType<_SrcTp>::type, dump.data());
mat_tmp.convertTo (slice, DataType<_DstTp>::type);
}
else
return 0;
if (show)
{
namedWindow ("slice2mat",1);
imshow ("slice2mat", slice);
waitKey ();
}
return 1;
}
5.3 Memory layout
In fact, Blitz++ has its own complex memory layout types and customizable.
For the case of two dimensional matrix:
Array<float,2> A(4/*firstDim*/, 3/*secondDim*/); // firstDim: number of rows // secondDim: number of columns // accessing to the (row, col) element float elem = A(row,col);
In order to access this element by pointer operations, we have to think in the reverse order:
float* mem_offset = (float*)(A.data()) + row*(A.extent(secondDim)) + col;
In other word, C++ stores array elements using row major order: go through 
for (int i = 0 ; i < A.extent(secondDim) ; ++i)
for (int j = 0; j < A.extent(firstDim) ; ++j )
{
float* ptr = (float*)(A.data()) + i*(A.extent(secondDim)) + j;
}
For the case of three dimensional matrix:
Array<float,3> A(4/*firstDim*/, 3/*secondDim*/, 2/*thirdDim*/); // firstDim: number of rows // secondDim: number of columns // thirdDim: number of slices // accessing to the (row, col,slice) element float elem = A(row,col,slice);
In order to access this element by pointer operations, we have to think in the reverse order:
float* mem_offset = (float*)(A.data()) + slice*(A.extent(thirdDim)) + row*(A.extent(secondDim)) + col;
In other word, C++ stores array elements using row major order: go through 
ffor (int i = 0 ; i < A.extent(thirdDim) ; ++i )
for (int j = 0 ; j < A.extent(secondDim) ; ++j)
for (int k = 0; k < A.extent(firstDim) ; ++k )
{
float* ptr = (float*)(A.data()) + i*(A.extent(thirdDim)) + j*(A.extent(secondDim)) + j;
}
The following code verifies the memory layout order in Blitz++
Array<int,2> A(3,3);
A = 5, 5, 5,
2, 2, 2,
3, 3, 3;
// accessing the (0,0)th element
cout << A(0,0) << endl;
Array<int,3> cube(3,3,3);
cube = 1, 1, 1,
1, 1, 1,
1, 1, 1,
2, 2, 2,
2, 2, 2,
2, 2, 2,
3, 3, 3,
3, 3, 3,
3, 3, 3;
Array<int,2> slicexy = cube(blitz::Range::all(), blitz::Range::all(),1);
cout << "Is slicexy storage contiguous:" << slicexy.isStorageContiguous() << endl;
cout << slicexy << endl;
slicexy = A;
cout << slicexy << endl;
Array<int,2> slicext = cube(blitz::Range::all(), 1, blitz::Range::all());
cout << "Is sliceyt storage contiguous:" << slicext.isStorageContiguous() << endl;
cout << slicext << endl;
Array<int,2> sliceyt = cube(1, blitz::Range::all(), blitz::Range::all());
cout << "Is sliceyt storage contiguous:" << sliceyt.isStorageContiguous() << endl;
cout << sliceyt;
Generally speaking, if you have a 3D matrix, let’s think it as a stack of 2D matrices along the bottom-up dimension. For each 2D matrix, the first dimension (the number of rows) is the second dimension of the 3D matrix and the second dimension (the number of columns) is the third dimension of the 3D matrix.
4.4 More notes on using Blitz++
The copy constructor of Blitz++ does not duplicate the data, it just refer to the same data.
Array<float,2> arr2(arr1);
The copy constructor has the same effect as assignment at the declaration.
Array<float,2> arr2 = arr1;
The assignment operator between two variables having the same size and each of them has its own data is in fact a copy data operator.
Array<float,2> arr2(arr1.rows(), arr1.cols()); arr2 = arr1; Array<float,2> arr3; arr3 = arr1.copy ();
Posted by vodinhphong 
