Sunday, July 14, 2013

Implicit surfaces and Level Set

Things are more complicated in discretized form. Say a line when discretized is not actually a line, Right(remembering the bresenham line algorithm ) ? Same thing goes to circle or any shapes. So for computing the associated properties we need more techniques. Implicit surface helps us to compute such properties.

Before I explain about implicit surfaces, you tell be what is the gradient of a scalar surface/curve.
Say you define a 3D surface by the equation x*x + y*y + z*z = 9, we can easily see that this surface is nothing but a sphere with radius 3. Right ?
So Q(x,y,z) = x*x + y*y + z*z - 9 = 0.
What is the gradient of Q then ? it represents the normal at any point on this surface , and it is (2x,2y,2z) (not normalized).

Why i said this was to show you how easy it is compute the gradient of a surface with an explicit equation. We can also easily compute other properties related with that surface like tangent,curvature,etc.

But what can you do if you don't have such explicit equations. In practical things will be like this.
In Implicit form we can define a shape implicitly. Our shape must be closed and non-self intersecting. With this agreement we can define the shape with following definitions.

Let P ( { Xi,Yi } ) be our point set which denotes the boundaries of our shape. We can define shape implicitly based on the following conditions.

1. For all points in shape boundaries  (Pi) = 0
2. For all other points outside the shape ∅(Pi) must be > 0 
3. For all other points inside the shape ∅(Pi) must be < 0 . (Conditions 2,3 can be interchanged though) 

Based on the earlier definitions consider the above picture. Pixel's with green boundary is our shape where ∅(Pi) will be 0. pixels having red color will have negative value, and rest of the pixels (blue) will have positive value. This is how we define implicit functions for complicated shapes. In the next post I will show you how we can numerically compute the properties of these shapes from this definitions also will introduce about level sets. It is not a big deal( Actually i had intention to write more about this, but I lost my mood so stopping now )

Saturday, June 15, 2013

Curvature Flow & Smoothing curves

I am getting addicted to curves. They are the perfect beautiful representation which we can numerically compute. I don't want to express more my feelings towards it which may be boring to you.

Coming to the topic, curvature flow is a kind of method to modify the curve.
Take any planar non intersecting curves , find the curvature at each point and multiply with it the normal there, then move the curve along that direction. This is the concept. It is analogs to the heat exchange. Heat will eventually spread uni-formally , no matter how you wrap it.

So take a curve {X(s),Y(s)} and it's curvature and normal , Say  {K(s)} and {N(s)}.
Then curvature flow vector is defined as {K(s)*N(s) }.

Using this technique we can smooth the curve from noise. Eventually this curve will become circular and will vanish. it is possible to know the exact point where it will vanish.

I just made a quick demo of this with matlab. See the images below

Original curve.
We can see that some edges are not smooth (think of it like created by noise).These spikes(non-smooth) in edges are not influencing our capability to detect the shape.We humans normally pick up a smooth shape from a contour. Now we want to remove some noise (or say smooth it) using curvature flow.
After 10 iteration
See that curve is more smooth now.

Curve after 50 Iterations
Curve is getting circular. Yeah, it will become circular at one point , because circle is the only one shape with a non-zero uniform curvature.

This is type of smoothing is also possible with Gaussian filtering along the curve,but with lesser accuracy. This type of technique can also be extended to 3D, Say you are having an object and you want to add a higher layer of layer/cover to it. Rather than just extending the surface along normal , use curvature flow.Then the surface will be more appealing and natural (I Guess).

Wednesday, May 22, 2013

Simple 2 Dimensional Curve Matching

In my previous post I have explained about curvature in depth. Now to the practical side, I created a simple application which will use these curvature to compare two simple planar curves. In the following video you can see two feature vectors indicates the similarity of curves;

See the video here.
The angle difference between the feature vectors indicates how similar those shapes are.In the video, you can also see that this matching is invariant to rotation and scale(when shape gets bigger, curvature will becomes lower). Right now the algorithm I used for computing the curvature 'Feature vector' is based on centroid. It needs to be refined further,But the underlying theory is very solid.Also the first impression giving me a very good hope on the concept.

However I am stopping my work on this concept, I don't have time to refine it. Next my target is 'level set methods' or solving the thin plate spline equation. The second one is duper super hard to fully understand , I already attempted it and lost my mind and motivation.Whenever I take it , suddenly everything becomes complicated, even my life (incidents!). So it is like the book of Amun-ra  But after looking it, I knew that i need to improve my 'calculus of variation' skills , and that topic is very nice.The same thing which helps to solve missile guidance problems!.

Curvature: Deriving the formula based on non-length parameter.

Matching curves is an important operation in computer vision. Several methods are there. But the important thing is how we are defining the curves.

This time i am writing about some concept which uses the curvature to do this matching operation. In some of my previous posts i had written about parameterizing the curve based on length(parameter 's').It is important to understand such concepts which will help us to apply some imagination.

Curvature 'k' is the magnitude of second derivative of a curve which is parameterized using length parameter('s')
In real scenarios(live)we will not have the luxury of getting curves which follows some strict equations.Most of the cases what we get is some points(of-course again corrupted by some noise).
We need to use a numeric approximation here to find the derivatives. That is simple,the hardest part is making a list with ordered point pairs in the correct order which will closely matches the property of the original curve
Once we have such a list, it is possible to use curvature for matching. Now you may think like curvature is the magnitude of second derivative based on 's' parameter. So do we need to again parametrize the curve to 's' ? (i had this thought ). No need , we can find a curvature based on non-length parameter too.
Here it is.
Say our curve is r(s) = { x(s),y(s) }, s is from 0 to length
Now we need to find dr/dt. ( 't' is non-s parameter )
dr(s)/dt = (dr/ds) * (ds/dt) [ using chain rule) ]
= t * v (lets call (dr/ds) as tangent (t vector, not same as parameter 't' which is scalar) w.r to s, and ds/dt as 'v' which will tell how wast length is changing with respect to 't']
Now we need to find d^2r(s)/dt^2
= d(t)/dt * v + dv/dt * t
= d(dr/ds)/ds * (ds/dt) * v + (dv/dt) * t [again using chain rule of diff., first differentiate with s, because r is a function of 's' , then differentiate 's' w.r to 't']
= d(dr/ds)/ds * v * v + (dv/dt) * t
= dt/ds * v^2 + (dv/dt) * t

Here dt/ds is actually the second derivative of 'r' with respect to parameter 's'. This quantity can be written as k*n,where k is a constant and n is a unit normal vector. So the equation will become
= k*n * v^2 + (dv/dt) * t

Now take cross product between dr(s)/dt , d^2r(s)/dt^2

dr(s)/dt X d^2r(s)/dt^2 = (t*v) X (k*n*v^2 + (dv/dt) *t )
 = (t*v) X (k*n*v^2) [because tXt is zero ]
= t X (k*n*v^3) [ in this only t, and n are vectors.It is very much clear.]
= (k*t*v^3) X (n)

So dr(s)/dt X d^2r(s)/dt^2 = kv^3t x n ---------------- (1)

vector t in the previous equation (1) is dr/ds, tangent vector of curve which is parameterized over 's' and is always unit length. so that means we can find 't' just by normalizing the vector 't'.

t = <dx(t)/dt,dy(t)/dt,0> / Norm[ <dx(t)/dt,dy(t)/dt,0> ]

the vectors 't' and 'n' lies in the x,y plane and also perpendicular. so take a vector 'z' which is also perpendicular to the x-y plane and its value is (0,0,1). so these three vectors (t,n,e) form a 3 dimensional coordinate system. 

By using normal cross product logic we can interpret vector 'n' as cross product between e and t.

n = e x t 

lets call dr(s)/dt simply as dr
lets call d^r(s)/dt^2 simply as drr  

dr x drr = kv^3t X (e x t )
=(dr x drr)/ v^3 = k* t x (e x t )
=(dr x drr)/ v^3 = k*( e(t.t) - t(t.e))
=(dr x drr)/ v^3 = k*( e(t.t) - t(t.e))
=(dr x drr)/ v^3 = k*( e(t.t) - t(0))
=(dr x drr)/ v^3 = ke
k = (e. (dr x drr)) / v^3
k = (<0,0,1> . (dr x drr)) / v^3

Now lets find dr x drr

if we know the curve 'r' , we can easily find dr and drr. first parameterize curve to 't' , mostly all curve's natural parameterization will be t which is non-arc length parameter. Say you are drawing some a curve with mouse, then points will be placed on the order you place it right ? it may not be possible for you to place points with equal interval considering the total length of curve.
    dr = < dx(t)/dt ,dy(t)/dt ,0 >
    drr = < d^2x(t)/dt^2 , d^2y(t)/dt^2 ,0 >

dr X drr = <0,0,(dx(t)/dt) * d^2y(t)/dt^2 - dy(t)/dt * d^2x(t)/dt^2> = <0,0,dx*dyy - dy*dxx>

so k = (<0,0,1>. <0,0,dx*dyy - dy*dxx>) / v^3
    k = (dx*dyy - dy*dxx>) / v^3

we need to get ride of v also. 'v' is ds/dt . what is the value of 's' , it is a function of length
s = Integrate [ Sqrt[ |dr/dt| ] ]  over the entire curve. Differentiating it again with respect to t , we get
ds/dt = Sqrt[ |dr/dt| ] = v

k = (dx*dyy - dy*dxx>) / Sqrt[ |dr/dt| ] ^3
k = (dx*dyy - dy*dxx>) / (dx^2 + dy^2)  ^(3/2)

This is the final equation for finding curvature of plane curve which is not parameterized on arc length parameter.This k is invariant to rotation, only problem lies in how efficiently we can parameterize the curve which follows the properties which we used to derive this 'k'.

see the images which shows the original curve, its curvature('k') plot, rotated image and its curvature plot. We can see that the curvature remains same.

original curve
curvature plot of original curve

Rotated curve

rotated curve's curvature plot

From the figures we can see that curvature remains same. we can use this feature for matching planar curves. I will post some video of that in next post.

Thursday, April 4, 2013

Wiener Filter + Inverse Filter, Contd.

In the last post i derived the formula for wiener filter. I am not getting enough time to write something here. This night i decided to write something.  If you carefully examine the wiener filter formula it can be seen that when the K is zero ( that is no noise),it act just an inverse filter.
     This means we can just divide the Fourier transform of the input signal(degraded) with the Fourier transform of the degrade function. This will produce the original data. You may be wondering how this works. It is based on convolution theorem.

In my pc mathematica stopped working. So i am switching to matlab from now on

Following mat lab snippet proves that convolution in spatial domain and frequency domain is same

%% in this example convolution signal is not centered ( so the source)
%% by this way we can avoid the one shifting(ifftshift) at the end.

data1 = [1,2,3,2,1,0,0];
kernel = [1,1,1];
kernel_padded = [1,1,1,0,0,0,0];
'convolution spatial'
d = fft(data1) .* fft(kernel_padded);

When you run it this will produce the result  1  3    6   7  6  3  1. (conv function will out two more zeros, but it is not a part of input data)

Using this idea inverse filter works. Since inverse filter uses a division problem can happen when the denominator is zero or say too high. because of this it is not stable as itself also with noise.

Here is the results and code

a = imread('f:\\ip\\imtests\\exotic.jpg');
d = rgb2gray(a);
degrade = [.1 .1 .1;.1 .5 .1;.1 .1 .1];
dim = conv2(d,degrade,'valid');
figure('Name','degraded image with some filter')

img_size = size(d) ;
img_size = size(dim); 
conv_mat = zeros(img_size(1),img_size(2));

% create conv matrix for inverse.
pw = floor(img_size(1))-3
ph = floor(img_size(2))-3
conv_mat = padarray(degrade,[pw,ph],'post');

fo = fft2(dim);
fc = fft2(conv_mat);
% the final division.
myimg = fo ./ (fc )

corrected = ifft2(myimg);
figure('Name','corrected with inverse filter')
imshow(real(corrected),[0 255])

Oringal Image

Corrected with Inverse

Don't expect this kind of correction all cases , because of that division things can go totally wrong . So one idea is to add some constant with denominator before division and play with it until you find some satisfactory results Or you need to check it against zero or non-numeric value.
Ok its time for me to stop. I just made this because of me :) . I was worried about not writing anything for sometime. That's why this quick post.  More interesting this are coming(to my mind)

Friday, March 8, 2013

Wiener Deconvolution: Deriving the final formula

Recently i watched a video which talks about noise and filtering.That why i am posting this.After searching the same in wiki, i didn't see a detailed derivation of the formula. so I thought to write this and post.

Wiener deconvolution is a frequency domain technique which helps to remove the noise if we know the degradation function(H(t)) in advance. The degradation function degrades the image quality. It can be blurring like Gaussian blur or it can be motion blur.

X(x,y) -> Passes through degradation function ->  Some noise gets added here(N) -> our final signal Y(x,y)

If we multiply the Y with Wiener filter then it will provide an approximation of X. Lets call this approximation as X$ .Now am going to write down the full derivation. a part of this is available in wiki. But those who are not remembering complex number maths may find it difficult.

Mean square error between approximation and original data is = E | X(x,y) - X$(x,y) | ^2

but we know X$ is  can be obtained by G(x,y) * Y(x,y)  [ * is convolution operator ]
 = E( | X(x,y) -  G(x,y) * Y(x,y)  |^2)
But if convert all to frequency domain convolution becomes multiplication, ie a(t)*b(t) = a(f)b(f)
Also we know Y is X * H + N
= E(| X -  G ( XH + N)) |^2)
= E(| X -  GXH - GN) |^2)
= E(| (1 -  G H)X - GN) |^2)

if A and B are two complex numbers | A B |^2 = (A B)(AB)* , here * is complex conjugate.

= E(  ((1-GH)X - GN)((1-GH)X - GN)* )
= E(  ((1-GH)X - GN)((1-G*H*)X* - G*N*)) , here all * means complex conjugate

Expanding this is a painful thing. but i am not giving up. So are you ?

= E( ((1-GH)X(1-GH)* XX* - (1-GH)XG*N* - GN(1-GH)*X* + GNG*N* )
= ( (1-GH)(1-GH)* E{|X|^2} - (1-GH)G* E{XN*} - G(1-GH)* E{NX*} + GG* E{|N|^2}
Lets assume noise and input signal has no relation,That why it is noise. So we can just ignore the terms E{XN*} ,E{NX*}.I think this is because at a particular frequency we can't find any relation between input data and noise. So the probability is null here. That why taking E( expectation) which sums the probability of values. here since we can't find any probability we can assume it 0.

We can denote S(f) = E{|X|^2} , it is the power spectrum of the complex values of the input data
We can denote N(f) = E{|N|^2} , it is the power spectrum of the complex values of the noise  data

so final equation becomes

minError(f) =  (1-GH)(1-GH)* S + GG*N
Now we need to find the minimum error which we can obtain for G. So following the basic calculus (like how we will find the min value of a function ) . careful about G* ,it is not G. So we can treat that as constant.

d(minError(f)/dG = -H(1-GH)*S + G*N = 0 ,
Now we need to find G from this. we will find it in a minute.

That is  G*N - H[1-GH]*S = 0

G*N = HS[1-GH]*
G*N = HS[1-G*H*]
G*N = S[H-G*HH*]
G*N = S[H-G*|H|^2 ]
G*N = SH-SG*|H|^2 ]
G*N = (G*) [ (1/G*)SH-S|H|^2 ]
N = [ (1/G*)SH - S|H|^2 ]
N + S|H|^2 = SH/G*
G*/SH  = 1/ (N+ S|H|^2 )
G* = SH / ( N+  S HH* ).         Take out s
G* =  SH / S( N/S+ HH* )
N/S is noise to signal ratio, We can replace it with a constant K

G* =  H / (K+ HH* )
G   = H*/ (K+H*H), H*H is |H|^2|
G   = H*/ (K+|H|^2)

Yes we reached the final answer , G = H*/(K+|H|^2).
if you know H, then you can find G, Then just multiply this with Fourier transformed input data. you will get a better picture. Thus you can remove some noise from the input data. In next post i will try to post some samples after applying G or wiener filter. Be ready for it.

Sunday, February 24, 2013

Histogram Equalization : A simple technique for improving low contrast images

I have been busy. Not getting much time to give enough attention to this. I am really worried about this and this makes me unhappy. But life must goes on. I can do what i want but i am not able to do.. It is an interesting thought. I am sure many of us still have this.(just had a wine)

What is histogram , it is a frequency showing some information in image. Due to some factors based on the equipment or lighting or whatever it is possible sometimes taken image will contain colors that uses near by gray values. That is is the image's gray levels are not really stretched to required level.

See some images here
These are images taken on low light or artificially generated. if you examine this these uses intensities that are around near by spectrum.That is why our eye is not able to perceive the amount of detail we want to interpret.  Histogram equalization is a very simple technique to stretch out this contrast by using a cumulative distributive function.

This CDF function is very simple, it means for the pixel with highest value it assigns a probability of 1, because it is cumulative  That is for the brightest pixel value( say 255) it sums probabilities of all pixels from 0 to 255 and it must be 1 it will be (for other pixel values it sum ups probabilities up to that pixel's value). So our CDF function is a monotonically increasing function.

I am attaching the mathematica code and results. so i don't have to explain.

yimage = Import["f:\\imtests\\lowContrast1.jpg"];
grayImg = ColorConvert[myimage, "GrayScale"];
mydata = ImageData[grayImg, "byte"];
temp = ImageDimensions[grayImg];
width = temp[[1]];
height = temp[[2]];
numberOfPixels = width*height;
histo = Array[0 &, 256];

(* compute histogram *)
m = mydata[[2]][[138]];
For[ i = 1 , i <= height, i++,
  For[ j = 1 , j <= width, j++,
    m = mydata[[i]][[j]];
    m = m + 1;
    histo[[m]] ++;

(* create cummlative distributive function to map pixels *)
cdf = Array[0 &, 256];
For[ i = 1, i <= 256, i++,
  sum = 0.0;
  For[ j = 1, j <= i, j++,
   sum = sum + histo[[j]] / numberOfPixels;
  sum = sum * 255;
  cdf[[i]] = sum;

(* plot cdf *)

(* equalize histogram *)

orignalData = mydata;
(* change pixels in original image based on cdf function *)
For[ i = 1 , i <= height, i++,
  For[ j = 1 , j <= width, j++,
    m = mydata[[i]][[j]] + 1;
    mydata[[i]][[j]] = cdf[[m]]
Image[orignalData, "byte"]
Image[mydata, "byte"]


CDF  function

Original image

Output image

See how good is the equalized image. but do not expect Histogram Equalization to perform this in all conditions. In general It is good for x-ray images. One other problem is it will also amplify the noise. So it better to remove the noise before passing to this.