Matrix Multiplication: 3x3 and 3x1

Formula:

| a11 a12 a13 |    | b1 |    | a11*b1 + a12*b2 + a13*b3 |
| a21 a22 a23 | x | b2 | = | a21*b1 + a22*b2 + a23*b3 |
| a31 a32 a33 |    | b3 |    | a31*b1 + a32*b2 + a33*b3 |

For arrays:

| a11 a12 a13 |      | a[0][0] a[0][1] a[0][2] |
| a21 a22 a23 |  =  | a[1][0] a[1][1] a[1][2] |
| a31 a32 a33 |      | a[2][0] a[2][1] a[2][2] |

| b1 |    | b[0] |
| b2 | = | b[1] |
| b3 |    | b[2] |

| c1 |    | a11 a12 a13 |    | b1 |    | a[0][0]*b[0] + a[0][1]*b[1] + a[0][2]*b[2] |
| c2 | = | a21 a22 a23 | x | b2 | = | a[1][0]*b[0] + a[1][1]*b[1] + a[1][2]*b[2] |
| c3 |    | a31 a32 a33 |    | b3 |    | a[2][0]*b[0] + a[2][1]*b[1] + a[2][2]*b[2] |

Code:   

c[0] = a[0][0]*b[0] + a[0][1]*b[1] + a[0][2]*b[2]
c[1] = a[1][0]*b[0] + a[1][1]*b[1] + a[1][2]*b[2]
c[2] = a[2][0]*b[0] + a[2][1]*b[1] + a[2][2]*b[2]

Matrix Inverse - 3x3

Operations of matrices of 3x3 are very common in day to day applications. Because the dimensions are already known, its better to avoid using a for-loop because of the extra computation involved in finding the subscript values.

Following is a quick formula for the inverse of a 3x3 matrix (implemented as a 2D array).

Formula

The inverse of a 3x3 matrix:
| a11 a12 a13 |-1                   |   a33a22-a32a23  -(a33a12-a32a13)   a23a12-a22a13   | 
| a21 a22 a23 |    =  1/DET *  | -(a33a21-a31a23)   a33a11-a31a13  -(a23a11-a21a13) |
| a31 a32 a33 |                      |   a32a21-a31a22  -(a32a11-a31a12)   a22a11-a21a12   |

where DET is the determinant of the matrix, i.e.

DET  =  a11(a33a22 - a32a23) - a21(a33a12 - a32a13) + a31(a23a12 - a22a13)

For 2D array

Moving on from matrices to 2D arrays

| a11 a12 a13 |        | mat[0][0] mat[0][1] mat[0][2] |

| a21 a22 a23 |   =   | mat[1][0] mat[1][1] mat[1][2] |

| a31 a32 a33 |        | mat[2][0] mat[2][1] mat[2][2] |

and,

| a11 a12 a13 |-1     | inv[0][0] inv[0][1] inv[0][2] |

| a21 a22 a23 |   =   | inv[1][0] inv[1][1] inv[1][2] |

| a31 a32 a33 |        | inv[2][0] inv[2][1] inv[2][2] |

Hence for code, the assignments are:

inv[0][0] = mat[2][2] * mat[1][1] - mat[2][1] * mat[1][2]

inv[0][1] = mat[2][1] * mat[0][2] - mat[2][2] * mat[0][1]

inv[0][2] = mat[1][2] * mat[0][1] - mat[1][1] * mat[0][2]

inv[1][0] = mat[2][0] * mat[1][2] - mat[2][2] * mat[1][0]

inv[1][1] = mat[2][2] * mat[0][0] - mat[2][0] * mat[0][2]

inv[1][2] = mat[1][0] * mat[0][2] - mat[1][2] * mat[0][0]

inv[2][0] = mat[2][1] * mat[1][0] - mat[2][0] * mat[1][1]

inv[2][1] = mat[2][0] * mat[0][1] - mat[2][1] * mat[0][0]

inv[2][2] = mat[1][1] * mat[0][0] - mat[1][0] * mat[0][1] 

DET =   mat[0][0]*inv[0][0] + mat[1][0]*inv[1][0] + mat[2][0]*inv[0][2]

Forgotten rules of C/C++: Part 1

Following is a list of important points people find very convenient to forget. This is the first article of a series

·   For bitwise operations, operand is promoted to "int" before evaluation.

unsigned char I = 0x80;
printf("%d", i<<1 nbsp="nbsp">  256

·   printf(5 + "intelligent");  => “ligent”

      Number specifies displacement to string pointer. This is the same as doing…

char *s = "intelligent";
s += 5;
printf("%s",s);

·    In a switch statement, initializations are allowed in the beginning, but they are NOT EXECUTED. The control is passed directly to an executable statement, i.e. matching case statement.

switch(1)
{
     printf("hello");                         => NOT PRINTED
     case 1:printf("case 1");break;
     case 2:printf("case 2");break;
}

·   The ++ operator means ``add one to a variable'' and DOES NOT work with constants.

int i = ++ 3 ;
printf("%d", i);                              =>  GARBAGE VALUE

·   Static arrays are constant! The base address cannot be modified. Any pointer arithmetic that attempts to do so causes a compilation error.

int arr[ 10 ] ;                     => “arr” is a constant. It is defined to be the address, &arr[ 0 ].

·   Array pointer logic for an array, “int a[10][10];”

      a+1
=> a[1]
=> &a[1][0]

·    In pointer arithmetic, addition and subtraction are valid operations, (as long as they don’t try to change the base address of the pointer) but pointer division and pointer multiplication are INVALID.

Code Snippet - Summed Area Table

SUMMED AREA TABLE
Using the method described here, I have written the following code to calculate a summed area table. A sample image table of 5x5 dimensions has been assumed for simplicity. The code snippet and its results are given below.

CODE
// assumed height and width of the input table
#define height 6
#define width 6

// input matrix - 5x5
long matrix[ height - 1][width -1] = {{5,2,3,4,1},{1,5,4,2,3},{2,2,1,3,4},{3,5,6,4,5},{4,1,3,2,6}};
// output matrix - 6x6
long sat[height][width];

void sat_matrix(){
 // formula variables
      int a=0, b=0, c=0, m=0;
      // matrix traversal loop for calculating the SAT
for(int i = 0; i < height; i++){
for(int j = 0; j < width; j++){
                     // following code picks up array elements within bounds and picks "zero"
                     // for values outside bounds.
a = (i-1>=0)?sat[i-1][j]:0;
b = (j-1>=0)?sat[i][j-1]:0;
c = ((i-1>=0)&&(j-1>=0))?sat[i-1][j-1]:0;
m = ((i-1>=0)&&(j-1>=0))?matrix[i-1][j-1]:0;
                      // ACTUAL FORMULA FOR SUMMED AREA TABLE
sat[i][j] = m + a + b - c;
}
}
}

PROCEDURE:
Use the function as it is and write supporting code. The code written was compiled using g++ compiler in cygwin environment

OUTPUT:

Summed Area Table

Definition

Summed Area Table is both an algorithm and data structure used in reference with the concept of Integral images. SAT is a name used to refer to both the method and the result of conversion of an image to an integral image.

SAT – The Algorithm

Summed Area Table (or Integral Image) is an algorithm applied on a 2-dimensional array of elements. It’s a simple, single pass algorithm to obtain the integral image values from the given pixel values of image.

SAT – The Data Structure

SAT also refers to the table of values generated after applying the conversion to Integral Image. This table of values is then used as input for improving the speed of more complicated operations.


Procedure
The algorithm takes as input, a table of order nxn and returns a table of order (n+1)x(n+1). The fundamental operation done here is to apply the following formula;


I(x,y) = i(x,y) + I(x-1,y) + I(x,y-1) - I(x-1,y-1)
Where, 
  i(x,y) = Element of image array i[x][y]
  I(x,y) =  Element of integral image array I[x][y]


The procedure can be better explained by the following code snippet

Integral Image

Definition

An Integral image is one that is conducive to frequent summation operations. Any rectangular subset of such an image can be evaluated in constant time. Such an image is achieved by converting its pixel values to a SAT (summed area table) by simple interpolation of known pixel values.

An Integral Image is defined as,

I(x,y) = i(x,y) + I(x-1,y) + I(x,y-1) - I(x-1,y-1)

Where,

i(x,y) = Pixel value of base image at (x,y)

I(x,y) = Pixel value of integral image at (x,y)

How this helps

Once the integral image has been calculated, the sum of any rectangular subset can be calculated by just 4 array points, i.e.,

SUM = (Bottom right + top left – top right – bottom left)

Application

Despite being a very simple transform (mathematically), calculating integral images helps perform many complex calculations with ease. Processing an integral image is much simpler than processing a normal pixel table. Hence the prerequisite step of several image processing algorithms is the conversion to an integral image.

Integral Images are useful for calculating HAAR wavelets, gradients, means and other measures. These are also used in image blurring, face recognition and other similar algorithms.

The concept of integral images can be easily extended to continuous domain (using limits) and multidimensional images. Also, the nature of summation operations can be modified to suit the algorithm, such as summation over non rectangular areas.

Method

Integral images are calculated using the following method.

How I Installed CUDA on my PC

Platform: Windows 7 (32 bit)

CUDA Hardware: None yet, will use simulator till then.

1 Prologue

      Before you start, get a clear idea about CUDA and its features from here.

2 Installation

2.1 Software setups

Basic step is to download the latest SDK and to choose an appropriate Toolkit according to your hardware’s compute capability. Download links are given here.

If you don’t have the hardware yet, you will need to use the emulator to compile and run programs (covered here). The emulator was however deprecated in the 3.x updates. Hence you need to download CUDA Toolkit 2.3 from here.

2.2  Installation Steps

Thankfully for windows, no post-installation configuration is required. Just run the setups and install-away.

3      Choose your language: C\C++

3.1 C++

Visual Studio will be used for C++ development for convenience. Express versions can be downloaded and registered for no cost. Or check with your countries IEEE MSDNAA alliance website if you are a member.

(If other, better IDE’s now have compatibility with CUDA; please notify me in the comments)

Downloads:

 Visual Studio2008

CUDA VS WIZARD (Win 32) 2.00 (or latest version)

  

3.2 C

C is best used on Linux or a native Linux environment. For windows users, the limited capability offered by the command line is satisfactory in this context. No separate setup is needed as the environment variables are already in place and CUDA’s compiler; nvcc can be invoked from the command line directly.

3.3 Java

I am, as of now, not committed to the idea of using Java for CUDA programs. But for reference, please go to this site.

After these steps have been completed, you are now ready to compile and execute CUDA programs.