Mathematics for Computer Graphics

GFD(Generalized Fourier Descriptor), Part 1

2012-01-23T15:10:00.000Z

it is interesting.. GFD is a rotation invariant. In GFD the orginal image is transformed to polar cordinate system. This provides the rotation invariant ability. I will explain how this happens. When we map normal image in cartesian coodinate system to polar coordinates, a rotation in Cartesian coordinate will cause translation/shifting in polar system. Remeber in polar system r, and theta are the axis. So when you rotate original image, R remain same so in effect the coresponding image is shifted(to left/right according to the direction of rotation) in polar system.

Still you may be wondering even if that is the case, why Fourier transform output remains same ? because we are operating on different data.For each rotation , the image is shifted in polar coordinate system.
If you thought like that, you are thinking.. Good.

I will explain the answer to the above puzzle.

If you remember one dimensional Fourier transform you can see that we are finding the dot product of our Image data with a number of Cos,Sin vectors (N diamensional, same as image size). The sin and cos are orthonormal , so the magnitude remains same. Please read my previous post http://cgmath.blogspot.com/2011/12/image-recognition-using-phase-only.html to get an intuitive grasp on DFT.

that is

Magnitude Of DFT[ {1,2,3,4,5,6}] is equal to Magnitude of DFT[ {4,5,6,1,2,3}]

Magnitude is Sqrt( i*i + j*j ) of the complex number.

Thus we get rotation in-variance. Now the next thing is Fourier descriptor. I will post about that after some days. It is simple( Really ? )

Just started with generalized Fourier Descriptor

2012-01-09T19:17:00.001Z

I started working on a new algorithm for image search. This is based on Fourier descriptor. Image is first converted in to a corresponding polar representation.This representation allows to get complete rotation in-variance. This is Because a rotation in Cartesian plane corresponds to a angle shifts in polar domain. So the DFT remains same. This is a simple and great idea.

I also started working on my Engine to add MD2 animation support. The image processing things takes too much thinking time also sometimes makes me dull. Thats why i started this to feed my interests.MD2 animation is simple to implement also to understand. Bone animation is difficult to understand if you are a beginner/intermediate in computer graphics. You may be wondering about which format should be used like that.. I had a looked that before some years. But at that time i didn't got time to implement it. Now sadly i am not remembering much.I should have implemented it that time.

Image recognition using phase-only correlation

2011-12-13T17:56:00.005Z

I have been working on image matching program which based on frequency analysis.Phase correlation principle is used here.

The concept is like this

Get different frequency information from image using DFT
Do phase correlation of source and templates frequencies (got from source and template images, both are n-dimensional)
After phase only correlation do inverse Fourier and find the peaks in the real part. that peaks shows the matching positions..

To really understand how this works , you need to know how DFT works. Not just that magic equation.
How it could transform the signal(here image) to freqency data ?infact we can make an intutive represenation of DFT equation in mind. Here it is in without any formulas.

You need to know the vector projetion operation first. If you don't know that my humble opinion is you better learn some basic vector maths right now. (vectors are everywhere.. Beware!!).
In DFT we have an image, which can be represented as an N-Dimensional vector.Then we take the dot product of this N-dimensional vector with an another N-Dimensional vector , lets call it 'Sin' vector, we also take dot product with an another vector, lets call it 'Cos' vector

that is (ImageN . Sin) + i (ImageN . Cos ) : ". is Dot product between vectors "

(Dot product is actually gives the distance in terms of vector, you can conceive it as projecting the image data to n-dimensional sine and cosine vectors which gives the real and imaginary parts). We do this operation for a number of frequencies, so the Sin,Cos Vectors changes giving different n-dimensional vectors.That's how you get the output N complex numbers (Now look refer the original DFT equation)
So what we just said is, we project(DOT) the image data to a number of vectors which are obtained from different sin,cos signals. Another point worth to remember is Sin and Cos vectors are always ortho normal ). I hopes i just wrote enough theory for DFT so that you can understand it intutively.Visualization is very important.

So back to matching thing., we do DFT of source and template image( both are same size, template can be padded with zero ). After this apply phase correlation using equation

(a+ib)(p+iq)* / (| (a+ib)(p+iq)* | )

(a+ib): complex number got from source DFT (so for NxN dimensional image , there will be N*N complex numbers,represented as 2D array NxN.)
(p+iq): complex number got from template DFT, (p+iq)* is conjugate operation.

What this equation does is it gives high values for signals where where peaks and bottoms correctly
aligns each other.

Now we do inverse DFT on the phase only correlated datas, and apply some threshold on the real part inorder to get the peaks which indicates the matched positions

Images from my tests.

Original image (after edge detection)

Template image (template image size must be same as source, so pad remaining area with zeros)

Result after phase only Correlation, see the peaks, You can see some invalid peaks also there (that is another story). I did this with mathematics. With 'mathematica' or 'matlab' implementation is easy. But i had to spend months to really understand the theory.

This method has advantage of being small invariant to smaller rotation and scale. But drawback is it needs higher time for processing due to DFT.This algorithm doesn't consider any shape information for matching, so this can give false results when too much edges present in source image.

That's all for now,this was a quick post. Good Luck with your projects.

A Board game within 3 hours

2011-06-26T17:04:00.001Z

Yes, I made a board game within 3 hours. I don't know what made me to code it. May be after playing some board games in mobile, I wished to create my own, especially those crystals graphics. I learnt to create some nice looking glass buttons with Gimp, it is just silly task but it is fun.

See game image shown below(the red text indicates the connected cells count to it), the game idea is like whoever first able to make 5 coins in a line(vertical,horizontal or diagonal) will win the game.

I coded the AI for this small game, it was fun and so easy, just count the number of opponent(player1) connected coins and calculate probability for each empty cell. Next is to calculate the probability of player2's connected cell, based on these determined the cell where the computer must put the coin. It works!!, Can be made better by adding some other ideas too(like considering the distance , so that if same probability comes it will choose the cell which is more connected).

[if you need source code mail me.]

Shape Context Matching

2011-05-11T18:04:00.001Z

I recently completed my shape context image matching project without much success :(.if you don't know about shape contexts , check out this link SC

The shape context can be created easily, For matching operation a bipartite based graph matching algorithm can be used. It can be done using Hungarian algorithm with complexity is less than the brute force method( O(n^2) ).

But after implementing Hungarian algorithm i am not able to find enough match point pairs between shapes. The problem with Hungarian method is that it can't give you the answer in fixed amount of time. It tries to optimize, sometimes never ending optimizations.

Still I believe it can be corrected(using a different approch to SC creation) , I need to work more. Not getting enough time.

A nice lecture about Hungarian Algorithm can be found at here

Face Detection using PCA.

2011-02-13T17:03:00.001Z

I have been working on my image processing library to support face detection. I started with basic method PCA ( Principal component analysis). Basically you need to have a set of images(20-50) with different lighting conditions etc. The next step is to create covariance matrix out of it and find the eigen vector. Then simply project the the image you need to check in to Eigen vectors and find the distance between them.Do some thresholding to classify it.One important thing is you don't have to take all eigen vectors,may be its better to sort (descending ) based on Eigen value and take only first 'N' vectors.

See the video

The difficult part in PCA may be to find the eigne vectors , QA algorithm seems a good choice. The current problem with running time.It takes almost 1 second to process 200X201 image. Roughly O(n^3) complexity. I am plan to implement it in CUDA,so that it can be used for real time detection.

Number Plate region extraction

2010-12-14T15:36:00.001Z

After a long break I started writing about my image processing studies.This time KD came with a project to extract number plate region from an image. The advantage of my method over other methods are the following.
1. fast processing
2. It can give you multiple regions in image if more than one number plate present
3. It can handle image rotation up to a certain degree( +- 35 ) .I used Eigen vectors.
4. No third party libraries like openCV or aforge ( yes some times i like to reinvent the wheels again )

See the video to see the project in action.

Although the number plate extraction parts works pretty good , I don't have a good OCR module. So i am having troubles to extract numbers from image. I tried using a simple back propagation neural network, its quality of recognition is not that great.Now i am trying to develop a rotation,scale invariant recognizer. It may take another 9 or eight months to do that. But if it works i think that would be a great achievement. I will try to post more updates here..

If you know any good optical character recognition library, please let me know.

Calculating the reflected ray/vector

2010-06-12T09:49:00.006Z

In computer graphics applications its often needed to calculate the reflection ray for example if you are writing a ray tracer, a shader for some advanced lighting , or environment mapping etc.

If you are writing shaders , there are standard library function to do that . In cg shading language there is a function reflect(also its more efficient than writing our own).

In this post rather than just giving the vector formula for reflection ray , i am trying to explain the simple mathematics behind that.

See the following image, I is the original ray, and R is the reflected ray which we need to found.N is the normal of the incident plane. P is the line perpendicular from normal to both rays. it is obvious that at both ends the length of P will be same.

Dot product between two unit vectors gives the cosine of the angle between them. So using this idea we can find

R = DotProduct[ I, N ] * N + P . ---> Eq(1)

We don't know P now. But I + P = N * DotProduct[ I,N].

So by rearranging P = N * DotProduct[ I,N] - I. Substituting the value of P now in equation(1) gives the final equation.

Here it is the final equation R = 2 * N * ( DotProduct[ I,N] ) - I

Color Image Segmentation using Meanshift Algorithm

2010-03-16T06:19:00.009Z

The meanshift method can be used to segment color image. In this method image pixels is treated as points in color space.In each iteration the meanshift vector is calculated for points which are inside the kernel radius.After that the the old kernel location is changed to meanshift vector's position. Color is also updated. This process continues untill both converge.

Original Image

Meanshift Filterd Image

As the iteration count increases , the same color segment which has same type of colors will get merged together. You can see the effect of meanshift filter on sachin's photo.Its like water painting (not exactly,there are other filters for that. )

How to check a Point inside a Triangle

2010-03-11T17:06:00.007Z

When i started to learn vector mathematics , i had found many methods to determine whether a point is inside a triangle or not.Easy method is not the fastest. This type of calculations are very important for making a fast graphics library.
Following are the methods which i remember now.
1. Using Cross Products

Take cross product between each triangle edge and the point to check. If the direction of cross product result is for the three edges that point is inside the triangle , otherwise not

2. Using the angle between the vectors made by point and triangle coordinates.If the sum of angle is 360 the point is inside , otherwise not.

3. Using Point and Plane test.See the figure. You need to know plane equation to understand it.

4. Using Bary centric Coordinates

This is the fastest of above. It involves checking the point in Bary centric coordinates system.

Real time motion segmentation

2010-02-19T06:52:00.007Z

Motion segmentation is one of the greatest challenges in computer vision. Accurate real time segmentation is difficult to achieve with ordinary camera. Using 3D camera we will get much better approximation .but 3D camera is costly.

i had done edge based motion detection , now i just did chain method of clustering on the edges to group the edges in to different groups. Its not accurate yet. with ordinary 2d camera it is difficult to cluster objects due to the perspective effect. But it can give better results if we could place the camera in somewhat higher position compared the the objects.

see the videos here.. ( edges in a group have same colors )

Updated motion segmentation , see the results in video below. I need to remove the noisy segments. Noise segments are created due to the illumination changes. Hopes it can be done with some pixel modeling schemes.

Boundary Approximation

2010-02-03T07:07:00.000Z

I have been playing with snakes( aka active contours ) for a quite long time.
Basic snake is not able to detect boundary correctly. But it is very fast.

I tried to extend the gradient vector field of the image using Laplace equation. it works. But it takes a few seconds to generate the vector field.. :( , now i am thinking about using distance transform to generate the vector field.

See the video below to see snake in action.

Tracking using Active shape models

2009-12-04T05:15:00.001Z

A few months back i had done a program which uses active shape models for tracking objects. Following video shows the result of it. Still I need to work on it to make it suitable for commercial use( but i don't have time :( , i need to study some other things related to gender classification ) , but the first looks very much promising to me.

Ray - Circle intersection test.

2009-10-09T08:06:00.001Z

If you are writing a game or any computer graphics applications you may want to check for intersections like triangle-triangle, triangle-ray, sphere-ray-sphere-sphere,Box-Triangle etc...There exists different algorithms for finding this.. Separating Axis theorem (SAT) seems to be a good option.

Following code written in c++ find intersection between a circle and line ( or ray-circle ). it will return true if there is intersection.Other function parameters are described in next lines.

vOrginCircle : center of the circle , r is the radius , vtLineStart: Line start point, vtLineEnd : Line end point, outIp : intersecting point..

As you can see it is easy to extend this method for finding ray-sphere intersection. just add a 'z' coordinate :)

bool IntersectLineCircle(vector vOrginCircle, double r,vector vtLineStart, vector vtLineEnd, vector& outIp)

{

vector vOrginNew = vOrginCircle - vtLineStart;

vector vtRay = (vtLineEnd - vtLineStart).normalize();

double b = vOrginNew.x* vtRay.x + vtRay.y * vOrginNew.y;

double c = r*r - ( vOrginNew.Norm2() );

double a = vtRay.Norm2(); // Norm2() is x*x + y*y

double diva = 1.0F/a;

// Negative doesn't have square root.

if( (b*b + a * c ) <= 0 ) return false;

double b4ac = sqrt( (b*b + a * c ) );

double l1 = (b - b4ac) * diva;

double l2 = (b + b4ac) * diva;

// we need the closest intersection point.. so find smaller l.

// To get the other end , just change the comparison operator.

double l = l1 <>

// find the intersecting point

outIp = (vtRay * l) + vtLineStart;

return true;

}

2009-05-08T04:40:00.000Z

Figure below shows the city renderd with my graphics engine. It has around 4000000 triangles.

Still selection is very fast. You can see a red trianlge on some building which is the selected one using mouse.I used Octree + Ray checking for selection.

2X2 Matrix inversion.

2009-04-24T04:28:00.001Z

Sometimes it is needed to find the inverse of a simple 2x2 matrix. For this there is a simple method exitst , we don't have to apply any gauss jordan method or crammers etc.

The Matrix A = [ a00 a01 ]

[ a10 a11 ]

then the inverse is 1 [ a11 -a01 ]

--- x [ ]

| A | [ -a10 a00 ]

| A | is determinent of matrix A which is a00 * a11 - a10 * a01

State estimation using Kalman filter

2009-04-20T06:18:00.001Z

Kalman filter helps to find the correct state of a linear system , provided the approximate state of the system at time 't'. It has wide application in the following fields Aerospace,computer vision, signal processing etc.

Let 'St' be the state variable at time t

St1 = A * St0 + V ,

V is the random noise with mean = 0

A is the gain Matrix of the system.

Let Mt be the measurment variable. Measurment is the approxiamte state at time t. Kalman filter can predict the correct state from this measurment( also depending on other params etc ).

Mt1 = H * St1 + W

H is the gain Matrix.

W is the random noise with mean = 0.

The interseting thing is usually Mt1 is known at time t1 and we want to find the st1 which is the correcte state at time t1 from this Mt1.

Kalman filter uses kalman gain to compute the state St.

St = St + Kt*( Mt - H*St ) , Kt is the gain matrix. and the new St is calculatd using the kalman gain Kt.

The noise is represented in this system using a covariace matrix. Let P represent the noise of function V ( system noise ) and R represent the noise during measurment.

I am omitting the detailed derivation of equations here. The final equation and steps are shown below.

At time 't' do the follwing to estimate the correct state St.

Update the state vector and measurment covariance matrix.

St = A * St

P = A * P * Transpose[P]

find the kalman gain for this.

K = A*P0*Transpose[H]* Inverse[ ( H*P0*Transpose[H] + R) ]

Mt = Do measurment here and assign to Mt.

Calcualte the estimated state St.

St = St + K* ( Mt - H*St );

Finally update the state updation covariance matrix P

P = P - K*(H*P )

Repeat these steps for each measurment.Following figure shows the result.

The measurment data 'Mt' , which is affected by noise.and rectified data St which is the output of Kalman filter and it gives us better prediction of the state.

Tracking results using meanshift and bhattacharya distance minimization technique.

2009-03-29T14:09:00.000Z

In my last post I had said about mean shift algorithm . I implemented tacking in video frames usiung meanshift algorithm and bhattacharya distance minimaization techniques( ask me if you need more details) . See the video below. ( Video resolution is poor.Please bear with me ).

In my opinion meanshift tracking is not so good,But still it can be useful on some occasions , where you don't want to know the shape of tracking object. It doesn't know anything about the shape of the object.There are other methods exists which are working in shape space and uses the edge detection for shape trasnformation. I have implemented one part of it, need to study some probabilistic modeling for completeing it.

Gradient vector calculation for non parametric data

2009-03-19T13:44:00.001Z

Gradient vector points to the nearst position where a high rate of change can b e seen.

Normally for a function it is obtained by applying the partial detivative operator. For example if we could represent the temperature inside a room by a function ( with parameter hight or time any conceiveing things ) , we could find the the gradient vector by applying partial derivative operator on that function. Using this technique we can go to the position where the temperature is highest than the current position. We just need to simply go along the direction of gradient vector.

Okay .. Things are simple when we can paramterize things. Practically parametrization of physical things is not easy. In such situation we can find the gradient density vector using kernal density gradiant estimator. In the follwing figure illustrates this mathematic technique. You can see that the circle is moved from it orginal position to the densest region (where more points present) . I calculated gradient vector and moved the circle through that direction , thus finally the circle reached at the densest region .

Used Epanechnikov Kernel for density estimation.

Arc ball rotation using Quaternion

2009-03-04T09:53:00.001Z

Some times we need to rotate the scene for a better view. Rather than providing separate slider ball control or similar mechanisam to rotate around X,Y,Z axis it is better to do rotation with mouse and it will give a intutive feel to user.

This is not a complex task. With some simple mathematical calculation we can achieve this. A unit Quaternion can represent any arbitary a 3D orientation. Quaternions are nothing but 4Diamensional complex number. Visualizing quaternion is difficult when compared to affine matrices.

Rotation Using Quaternion

When user clicks on the scene we will get the points in window coordinates. we need to convert it into the correct world coordinates. If you are familar with opengl you can use gluUnProject API for this. Otherwise you can simply mul

tiply with corresponding matrices to get this.

So When user dragging the mouse we will get two points in world space ( since dragging is not a continus operation on computer unlike the physical draging process). We can generate a axis from this two points ( vectors ) by taking the cross product between them. The angle of rotation is nothing but arc cosine of dot prodcut of first and second vector.

Now we can create a quaternion which represents the rotation around that axis.

Following code ilustrates this process. I used my Quaternion & vector classes to simplify the process.

Math3D::vector vtFrom;

Math3D::vector vtTo;

gluUnProject(m_PrePoint.x,m_PrePoint.y,0,fModeView,fProjection,iViewPort,&vtFrom.x,&vtFrom.y,&vtFrom.z);

gluUnProject(point.x,point.y,0,fModeView,fProjection,iViewPort,&vtTo.x,&vtTo.y,&vtTo.z);

vtFrom.y = -vtFrom.y;

vtTo.y = -vtTo.y;

vtFrom.normalize();

vtTo.normalize();

Math3D::vector vAxis = vtFrom.Cross( vtTo );

vAxis.normalize();

float fTheta = acos(vtFrom.Dot(vtTo));

Math3D::Quaternion qTmp;

qTmp.FromAxisAngle( vAxis, fTheta *2 );

qTmp.Normalize();

// multiply the current quaternion with the new one.

m_Quat *= qTmp;

m_PrePoint = point;

// normalise it.

m_Quat.Normalize();

Now the scene can be rotated with quaternion m_Quat and the scene will be rotated according the mouse movement made by the user.

Pseudo Inverse calculation

2009-03-03T05:19:00.001Z

Pseudo inverse matrix helps to get a satisfactory solution to linear equation which don't have any exact solution. This technique is used for curve fitting .

The pseudo inverse can be calculated like this . If A is the matrix

The pseudo inverse is this : Inverse[Transpose[a] * A ] * Transpose[A]

Normal inverse exists only for square matrix. So we can call this pseudo inverse as a general matrix inverse method.

Consider the following linear equations

5X1 + 4 X2 = 9

3X1 + 2X2 = 5

X1 + x2 = 3

Actually this set of equations doesn't have any exact solution.

A X R

5 4 x1 9

3 2 * x2 = 5

1 1 3

Clearly we can't calculate the normal inverse of A because it is not square matrix.

The pseudoinverse of A is

-.5 1.5 -1.0

.833 -1.833 1.333

So [Pseudoinverse]* [R] give the values for X1 and X2, which is the nearst solution for the equations.

Eigen Vector

2009-01-23T11:43:00.001Z

Today i just created one application to demonstrate the use of eigen vectors.

It is only for 2 diamensional things. Beacuse in 3D calculating eigen vector is a complex process ( It depends), involving linear equation solving , and cubic equations' root finding.

You canse see the eigen axis can be used for finding the OBB( oriented bounding boxes) for a set of points. This has significance inmportant in Computer graphics.

One more interesting thing is eigen vector of a rotaion matrix is actually its axis of rotation (similar to curl of a vector field , but curl calculation is easier).

Accelerated 3D rendering.

2008-12-01T08:31:00.000Z

Following screen shots are from my rendering engine. It uses octree and occluion culling to accelerate the rendering process. It took almost 1 month to complete. Still some optimaization works are pending like use SIMD instructions etc.

In first image you can see some buildings(of course the building are not realistic, but it doesn't matter ) and the camera is placed near a wall

When we move camera a little downwards , the camera got occluded totaly by that wall . This cause to cull a lot of polygons in the scene. The second image shows that scene (scene is a bit zoomed out to get a clear picture). Also you can see the frame rate in title bar(for the second image it is 360). Frame rate is increased by a great amount.

Now I have plan to write a tool for manually selecting occluder polygons for a scene.

Projecting vectors

2008-11-07T12:06:00.000Z

Sometimes we need to project one vector over another. See the picture below , we can see a vector u. Now we need to project this u on to x axis, .

This can be done with this equation : p = x * u . (x / Norm[x] ) . p is the result

"." means Dot product. What that equation means is we need to normalize x first ( in general to which we are projecting u vector ) and find the dot product between u and normalized x , and mulitply our x(destingation ) vector by this scalar value. Now we will get the result.

2008-11-04T04:30:00.000Z

Finding rotation axis of a rotation matrix.

A rotation matrix can be represented as

x1 x2 x3
y1 y2 y3
z1 z2 z3

Where (x1, y1, z1) are the first vector (x2, y2, z2), (x3, y3, z3) are the second and third vectors.
Note this matrix is column major oriented.

When you multiply this matrix by an arbitrary vector in space you will get the rotated vector coordinates.

In some situations it is important to know the axis of rotation. Here I am going to present one method to get the rotation axis from rotation matrix.

Consider this matrix which represents a rotation around axis ( 0.7017,0.7017,0 ) , this actually a vector in x-y plane rotated by 45 degree, Like this

The rotation matrix for rotating around this vector by 45 degree is shown below (which I calculated using quaternion)

0.8536 0.146 0.500
0.1460 0.850 -0.500
-0.500 0.500 0.707

When we multiply a arbitrary vector (x,y,z) by this matrix we gets the vector function F as below

(X 0.8536 + y 0.146 + z 0.5) i + ( x 0.146 + y 0.85 – z 0.5 ) j + ( -x 0.5 + y 0.5 + y 0.707 ) k

I,j,k are the base vectors.

Taking curl of this vector function F, we get the axis of rotation.

I calculated the curl of that function and it is (1,-1, 0). After normalizing this vector we can see the result as (.707,-.707, 0),. Earlier we had created the rotation matrix to rotate around the axis ( .707,.707,0). Our answer curl is in the opposite direction of (.707,.707,0) , but it is still correct.

Only problem is with the direction of rotation that is clockwise or anti clock wise..