Thread Subject:

Implementing Mean Curvature Flow (MCF) for an image in Matlab

Hi I need to write a Matlab program in which I need to be able to calulate the gaussian gradient of the image, followed by computation of the Mean Curvature Flow (MCF) using this formula :

k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2)(Iy^2)]^ 3/2}

I am able to calculate the gaussian gradient but I need some help to understand this formula and suggestions on how I could implement this using Matlab. Any guidance/ help that I can get on this is much appreciated!

Thankyou!

"Padmavathy Subramanian" wrote in message <itq9tr$4a3$1@newscl01ah.mathworks.com>...
> Hi I need to write a Matlab program in which I need to be able to calulate the gaussian gradient of the image, followed by computation of the Mean Curvature Flow (MCF) using this formula :
>
> k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2)(Iy^2)]^ 3/2}
>
> I am able to calculate the gaussian gradient but I need some help to understand this formula and suggestions on how I could implement this using Matlab. Any guidance/ help that I can get on this is much appreciated!
>
> Thankyou!
- - - - - - - - - - - -
  That formula doesn't look right to me. Curvature should be inversely proportional to a distance measure, whereas your expression looks proportional to distance squared if I interpret its symbols correctly. My guess is that at least an addition sign has been lost in the denominator expression.

  Could you please give us a reference to the formula.

Roger Stafford

"Padmavathy Subramanian" wrote in message <itq9tr$4a3$1@newscl01ah.mathworks.com>...
> Hi I need to write a Matlab program in which I need to be able to calulate the gaussian gradient of the image, followed by computation of the Mean Curvature Flow (MCF) using this formula :
>
> k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2)(Iy^2)]^ 3/2}
>
> I am able to calculate the gaussian gradient but I need some help to understand this formula and suggestions on how I could implement this using Matlab. Any guidance/ help that I can get on this is much appreciated!
>
> Thankyou!
- - - - - - - - -
  Yes, that formula you gave is undoubtedly lacking an addition sign, Padmavathy. It should be (in my opinion):

 k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2) + (Iy^2)]^3/2}

or if Ix, Iy, Ixx, Ixy, and Iyy are matrices, the matlab format would be:

 k = (Ix.^2.*Iyy-2*Ix.*Iy.*Ixy+Iy.^2.*Ixx)./(Ix.^2+Iy.^2).^(3/2);

where Ix, Iy are partial first derivatives (gradient components) and Ixx, Ixy, and Iyy are partial second derivatives.

  This can be interpreted as the curvature of contours of constant intensity:

 I(x,y) = c.

The sign of k will be positive if the contour bends toward the area of lower intensity and negative otherwise.

  You can derive the formula by the appropriate total and partial differentiation of I(x,y) = c and applying the well-known formula in calculus for the curvature of a curve.

  I presume the above is what is meant by "mean curvature flow" as applied to images, though I am unfamiliar with it by that name.

Roger Stafford

Hi Roger,

Sorry about the typo error in the formula. The formula you gave is correct.

Below are the steps that I calculated in my matlab program
1. Compute Ix and Iy
Ix=imfilter(I,hx,'replicate','conv');
Iy=imfilter(I,hy,'replicate','conv');
where hx and hy are 2D gaussian kernels along x direction and y direction

2. Compute Ixx, Ixy and Iyy
I get stuck in calulating these three matrices and I need advise on how to do it.


"Roger Stafford" wrote in message <itqvji$et2$1@newscl01ah.mathworks.com>...
> "Padmavathy Subramanian" wrote in message <itq9tr$4a3$1@newscl01ah.mathworks.com>...
> > Hi I need to write a Matlab program in which I need to be able to calulate the gaussian gradient of the image, followed by computation of the Mean Curvature Flow (MCF) using this formula :
> >
> > k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2)(Iy^2)]^ 3/2}
> >
> > I am able to calculate the gaussian gradient but I need some help to understand this formula and suggestions on how I could implement this using Matlab. Any guidance/ help that I can get on this is much appreciated!
> >
> > Thankyou!
> - - - - - - - - -
> Yes, that formula you gave is undoubtedly lacking an addition sign, Padmavathy. It should be (in my opinion):
>
> k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2) + (Iy^2)]^3/2}
>
> or if Ix, Iy, Ixx, Ixy, and Iyy are matrices, the matlab format would be:
>
> k = (Ix.^2.*Iyy-2*Ix.*Iy.*Ixy+Iy.^2.*Ixx)./(Ix.^2+Iy.^2).^(3/2);
>
> where Ix, Iy are partial first derivatives (gradient components) and Ixx, Ixy, and Iyy are partial second derivatives.
>
> This can be interpreted as the curvature of contours of constant intensity:
>
> I(x,y) = c.
>
> The sign of k will be positive if the contour bends toward the area of lower intensity and negative otherwise.
>
> You can derive the formula by the appropriate total and partial differentiation of I(x,y) = c and applying the well-known formula in calculus for the curvature of a curve.
>
> I presume the above is what is meant by "mean curvature flow" as applied to images, though I am unfamiliar with it by that name.
>
> Roger Stafford

"Padmavathy Subramanian" wrote in message <itvecg$3pn$1@newscl01ah.mathworks.com>...
> Below are the steps that I calculated in my matlab program
> 1. Compute Ix and Iy
> Ix=imfilter(I,hx,'replicate','conv');
> Iy=imfilter(I,hy,'replicate','conv');
> where hx and hy are 2D gaussian kernels along x direction and y direction
>
> 2. Compute Ixx, Ixy and Iyy
> I get stuck in calulating these three matrices and I need advise on how to do it.
- - - - - - - - - -
  For Ixx, Ixy, and Iyy I would try using the three second partial derivatives of the 2D gaussian distribution function to make gaussian kernels hxx, hxy, and hyy in a manner analogous to the way you created hx and hy. Guanglei Xiong has a FEX file written in 2005 which finds the gaussian gradient with hx and hy, probably much as you have done. If you follow that example but use second partial derivatives instead of first partial derivatives, I would think you could create Ixx, Ixy, and Iyy without any great difficulty. I would guess you could use the same sigma as in Ix and Iy. Of course you will have to arrange that everything is properly scaled.

  You probably already have this FEX link, but in case you don't, it is located at:

http://www.mathworks.com/matlabcentral/fileexchange/8060-gradient-using-first-order-derivative-of-gaussian/content/gaussgradient/gaussgradient.m

Roger Stafford

Hi Roger,
Thankyou for the suggestion!!
I understand how the gaussian function works by itself but i'm not able to relate it to the matlab codes and how do they then derive at the gaussian kernel ?
Any explanation or help is much appreciated.. !

"Roger Stafford" wrote in message <iu1a35$k5o$1@newscl01ah.mathworks.com>...
> "Padmavathy Subramanian" wrote in message <itvecg$3pn$1@newscl01ah.mathworks.com>...
> > Below are the steps that I calculated in my matlab program
> > 1. Compute Ix and Iy
> > Ix=imfilter(I,hx,'replicate','conv');
> > Iy=imfilter(I,hy,'replicate','conv');
> > where hx and hy are 2D gaussian kernels along x direction and y direction
> >
> > 2. Compute Ixx, Ixy and Iyy
> > I get stuck in calulating these three matrices and I need advise on how to do it.
> - - - - - - - - - -
> For Ixx, Ixy, and Iyy I would try using the three second partial derivatives of the 2D gaussian distribution function to make gaussian kernels hxx, hxy, and hyy in a manner analogous to the way you created hx and hy. Guanglei Xiong has a FEX file written in 2005 which finds the gaussian gradient with hx and hy, probably much as you have done. If you follow that example but use second partial derivatives instead of first partial derivatives, I would think you could create Ixx, Ixy, and Iyy without any great difficulty. I would guess you could use the same sigma as in Ix and Iy. Of course you will have to arrange that everything is properly scaled.
>
> You probably already have this FEX link, but in case you don't, it is located at:
>
> http://www.mathworks.com/matlabcentral/fileexchange/8060-gradient-using-first-order-derivative-of-gaussian/content/gaussgradient/gaussgradient.m
>
> Roger Stafford

Given
Ix=imfilter(I,hx,'replicate','conv');
Iy=imfilter(I,hy,'replicate','conv');

For Ixx Iyy, i use below formulas
Ixx = imfilter(Ix,hx,'replicate','conv');
Iyy = imfilter(Iy,hy,'replicate','conv');

For Ixy, i use two formulas
Ixy1 = imfilter(Ix,hy,'replicate','conv');
Ixy2 = imfilter(Iy,hx,'replicate','conv');

But Ixy1 and Ixy2 are different. Please advise if i use the formula correctly.

Thanks


"Padmavathy Subramanian" wrote in message <iui8kf$k1h$1@newscl01ah.mathworks.com>...
> Hi Roger,
> Thankyou for the suggestion!!
> I understand how the gaussian function works by itself but i'm not able to relate it to the matlab codes and how do they then derive at the gaussian kernel ?
> Any explanation or help is much appreciated.. !
>
> "Roger Stafford" wrote in message <iu1a35$k5o$1@newscl01ah.mathworks.com>...
> > "Padmavathy Subramanian" wrote in message <itvecg$3pn$1@newscl01ah.mathworks.com>...
> > > Below are the steps that I calculated in my matlab program
> > > 1. Compute Ix and Iy
> > > Ix=imfilter(I,hx,'replicate','conv');
> > > Iy=imfilter(I,hy,'replicate','conv');
> > > where hx and hy are 2D gaussian kernels along x direction and y direction
> > >
> > > 2. Compute Ixx, Ixy and Iyy
> > > I get stuck in calulating these three matrices and I need advise on how to do it.
> > - - - - - - - - - -
> > For Ixx, Ixy, and Iyy I would try using the three second partial derivatives of the 2D gaussian distribution function to make gaussian kernels hxx, hxy, and hyy in a manner analogous to the way you created hx and hy. Guanglei Xiong has a FEX file written in 2005 which finds the gaussian gradient with hx and hy, probably much as you have done. If you follow that example but use second partial derivatives instead of first partial derivatives, I would think you could create Ixx, Ixy, and Iyy without any great difficulty. I would guess you could use the same sigma as in Ix and Iy. Of course you will have to arrange that everything is properly scaled.
> >
> > You probably already have this FEX link, but in case you don't, it is located at:
> >
> > http://www.mathworks.com/matlabcentral/fileexchange/8060-gradient-using-first-order-derivative-of-gaussian/content/gaussgradient/gaussgradient.m
> >
> > Roger Stafford

"Padmavathy Subramanian" wrote in message <iui8kf$k1h$1@newscl01ah.mathworks.com>...
> Hi Roger,
> Thankyou for the suggestion!!
> I understand how the gaussian function works by itself but i'm not able to relate it to the matlab codes and how do they then derive at the gaussian kernel ?
> Any explanation or help is much appreciated.. !
- - - - - - - - - - -
  The following is my version of what I believe you need for the five gaussian kernels in your curvature problem.

  The notion of a derivative implies the presence of a distance measure. In this code I have assumed that the distance between pixels is to be regarded as a unit distance. You should be able to continue with that assumption up to the evaluation of curvature. At this point if you wish a different unit of distance, you can easily rescale that initial curvature value, which is inversely proportional to distance.

  You will notice that initial evaluations below of the five kernels, hxx, hxy, hyy, hx, and hy are subjected to a rescaling and adjustment process. This is because the theoretical forms will possess a certain amount of error due to its use with discrete points over a finite span. After these corrections the kernels should have the property that any image intensity function of the quadratic form

 I(x,y) = A*x^2/2 + B*x*y + C*y^2/2 + D*x + E*y + F

after convolution with these respective five kernels will result in the five partial derivatives, Ixx = A, Ixy = B, Iyy = C, Ix = D, and Iy = E at the point x = 0, y = 0, as you would expect, provided the kernel does not extend beyond the image boundary. The quantities x and y here are considered to be in units of pixel distance.

  I don't know the number of kernel pixels you have in mind, but I recommend that the standard deviation, sigma, you use be one-third or less of the pixel distance, m, from the center to the edge of the square kernel. The kernels will have a 2*m+1 by 2*m+1 size.

  Here is the kernel generation code:

 sigma = 4; % Standard deviation - chosen by the user
 m = 12; % Form kernels of size 2*m+1 by 2*m+1 - also chosen by user
 [x,y] = meshgrid(-m:m);
 h = exp(-(x.^2+y.^2)/(2*sigma^2)); % Start with 2D gaussian distrib.
 hx = x.*h;
 hx = hx/(-sum(sum(x.*hx))); % Rescale hx appropriately
 hy = hx.';
 hxx = x.^2.*h;
 t = [sum(sum(x.^2.*hxx)),sum(sum(x.^2.*h)); ...
      sum(sum(y.^2.*hxx)),sum(sum(y.^2.*h))]\[2;0];
 hxx = t(1)*hxx+t(2)*h; % Use approp. lin. comb. of x.^2.*h and h
 hyy = hxx.';
 hxy = x.*y.*h;
 hxy = hxy/sum(sum(x.*y.*hxy)); % Rescale hxy

Notice that the normalizing factor, 2*pi*sigma^2, was not used in h above because subsequent rescaling made that superfluous.

  I would recommend you use surf(x,y,hx), surf(x,y,hxx), etc. to visualize these five kernels. You will notice that if sigma is chosen too large, these surface kernel plots will begin to show an increasing departure from an ideal shape.

Roger Stafford

"Roger Stafford" wrote in message <ium59t$iaj$1@newscl01ah.mathworks.com>...
> sigma = 4; % Standard deviation - chosen by the user
> m = 12; % Form kernels of size 2*m+1 by 2*m+1 - also chosen by user
> [x,y] = meshgrid(-m:m);
> h = exp(-(x.^2+y.^2)/(2*sigma^2)); % Start with 2D gaussian distrib.
> hx = x.*h;
> hx = hx/(-sum(sum(x.*hx))); % Rescale hx appropriately
> hy = hx.';
> hxx = x.^2.*h;
> t = [sum(sum(x.^2.*hxx)),sum(sum(x.^2.*h)); ...
> sum(sum(y.^2.*hxx)),sum(sum(y.^2.*h))]\[2;0];
> hxx = t(1)*hxx+t(2)*h; % Use approp. lin. comb. of x.^2.*h and h
> hyy = hxx.';
> hxy = x.*y.*h;
> hxy = hxy/sum(sum(x.*y.*hxy)); % Rescale hxy
- - - - - - - - - - -
  The following code produces the same kernels, hx, hy, hxx, hyy, and hxy, (within round off error) as in the previous code but does so with fewer computations, in case you are interested in efficiency in generating kernels.

 sigma = 4; % Standard deviation - chosen by the user
 m = 12; % Form kernels of size 2*m+1 by 2*m+1 - also chosen by user
 x = -m:m;
 h = exp(-x.^2/(2*sigma^2));
 a = sum(h);
 b = sum(x.^2.*h);
 c = sum(x.^4.*h);
 d = 2/(a*c-b^2);
 e = b/a*d;
 t = x.*h/b;
 hxy = t.'*t;
 hx = h.'*(t/(-a));
 hy = hx.';
 hxx = h.'*(d*x.^2.*h-e*h);
 hyy = hxx.';

Roger Stafford

Thank you so much for your code.

I have verified with my code and both can create a similar Gaussian kernels hx,hy,hxx,hxy,hyy.

Now, i can use these kernels to convolute with the input image I to generate Ix,Iy,Ixx,Ixy,Iyy.

What i notice is Ix And Iy contains a lot of 0s. For example at i=5,j=5, Ix(5,5)=Iy(5,5)=0. Therefore, what puzzles me is how the curvature can be computed, given the below formula:
 k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2) + (Iy^2)]^3/2}

because now the denominator is 0 (Ix^2+Iy^2).

I really need an explanation for this. Please help!


 

"Roger Stafford" wrote in message <iuos54$ddp$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <ium59t$iaj$1@newscl01ah.mathworks.com>...
> > sigma = 4; % Standard deviation - chosen by the user
> > m = 12; % Form kernels of size 2*m+1 by 2*m+1 - also chosen by user
> > [x,y] = meshgrid(-m:m);
> > h = exp(-(x.^2+y.^2)/(2*sigma^2)); % Start with 2D gaussian distrib.
> > hx = x.*h;
> > hx = hx/(-sum(sum(x.*hx))); % Rescale hx appropriately
> > hy = hx.';
> > hxx = x.^2.*h;
> > t = [sum(sum(x.^2.*hxx)),sum(sum(x.^2.*h)); ...
> > sum(sum(y.^2.*hxx)),sum(sum(y.^2.*h))]\[2;0];
> > hxx = t(1)*hxx+t(2)*h; % Use approp. lin. comb. of x.^2.*h and h
> > hyy = hxx.';
> > hxy = x.*y.*h;
> > hxy = hxy/sum(sum(x.*y.*hxy)); % Rescale hxy
> - - - - - - - - - - -
> The following code produces the same kernels, hx, hy, hxx, hyy, and hxy, (within round off error) as in the previous code but does so with fewer computations, in case you are interested in efficiency in generating kernels.
>
> sigma = 4; % Standard deviation - chosen by the user
> m = 12; % Form kernels of size 2*m+1 by 2*m+1 - also chosen by user
> x = -m:m;
> h = exp(-x.^2/(2*sigma^2));
> a = sum(h);
> b = sum(x.^2.*h);
> c = sum(x.^4.*h);
> d = 2/(a*c-b^2);
> e = b/a*d;
> t = x.*h/b;
> hxy = t.'*t;
> hx = h.'*(t/(-a));
> hy = hx.';
> hxx = h.'*(d*x.^2.*h-e*h);
> hyy = hxx.';
>
> Roger Stafford

"Padmavathy Subramanian" wrote in message <iv1pjj$62o$1@newscl01ah.mathworks.com>...
> Thank you so much for your code.
> I have verified with my code and both can create a similar Gaussian kernels hx,hy,hxx,hxy,hyy.
> Now, i can use these kernels to convolute with the input image I to generate Ix,Iy,Ixx,Ixy,Iyy.
> What i notice is Ix And Iy contains a lot of 0s. For example at i=5,j=5, Ix(5,5)=Iy(5,5)=0. Therefore, what puzzles me is how the curvature can be computed, given the below formula:
> k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2) + (Iy^2)]^3/2}
> because now the denominator is 0 (Ix^2+Iy^2).
> I really need an explanation for this. Please help!
- - - - - - - - - -
  I have no way of advising you on that because I have not seen the data you are processing. However I find it rather strange that Ix and Iy should both be exactly zero at the same point in the image. Are you sure about that? Perhaps they simply look zero in a display because they are small in comparison to other values being displayed along with them. Try using "format long" to display individual values of Ix and Iy at those same points. I venture to guess they are not exactly zero.

  Also bear in mind that because you found it necessary to do filtering with gaussian kernels, the results are bound to only be approximate. You may have to experiment with choosing optimum sigma values and best values for the range of the kernels (what I called 'm') to get the best filtering.

  In evaluating the accuracy of these gaussian kernels, I urge you to make that test I described earlier. Create an artificial image of a quadratic form:

 I(x,y) = A(x-x0)^2/2+B*(x-x0)*(y-y0)+C*(y-y0)^2/2+D*(x-x0)+E*(y-y0)+F

for some point (x0,y0) and for various arbitrary values of A, B, C, D, E, and F, and see how close Ixx(x0,y0), Ixy(x0,y0), Iyy(x0,y0), Ix(x0,y0), and Iy(x0,y0) are respectively to A, B, C, D, and E. I designed the code I wrote for you to do well on such a test.

Roger Stafford

On Jul 6, 9:57 am, "Padmavathy Subramanian" <sso...@gmail.com> wrote:
> Thank you so much for your code.
> I have verified with my code and both can create a similar Gaussian kernels hx,hy,hxx,hxy,hyy.
> Now, i can use these kernels to convolute with the input image I to generate Ix,Iy,Ixx,Ixy,Iyy.
> What i notice is Ix And Iy contains a lot of 0s. For example at i=5,j=5, Ix(5,5)=Iy(5,5)=0. Therefore, what puzzles me is how the curvature can be computed, given the below formula:
>  k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2) + (Iy^2)]^3/2}
--------------------------------------------------------------------------------------------
I've only briefly dived into understanding this thread so perhaps I
shouldn't say anything but...it appears that your lx and lx are 2D
arrays so make sure you use the "dot" form of the operations if you
want to perform the operation element-by-element, for example lx.^2
and not lx^2, and use ./ and not /. I'm also not sure why you're
stuffing the numerator and denominator into cells with {} rather than
using() but maybe that's just an equation like you'd see in a math
book and not real MATLAB code (is so, ignore what I said).
k = (Ix .^ 2 .* Iyy - 2 .* Ix .* Iy .* Ixy + Iy .^ 2 .* Ixx) / (Ix .^
2 + Iy .^ 2) .^ (3/2)

Hi Roger,

There is a region in the image with only 1 colour (white in this case) and that's why Ix = Iy = 0 exactly....


"Roger Stafford" wrote in message <iv22gm$4be$1@newscl01ah.mathworks.com>...
> "Padmavathy Subramanian" wrote in message <iv1pjj$62o$1@newscl01ah.mathworks.com>...
> > Thank you so much for your code.
> > I have verified with my code and both can create a similar Gaussian kernels hx,hy,hxx,hxy,hyy.
> > Now, i can use these kernels to convolute with the input image I to generate Ix,Iy,Ixx,Ixy,Iyy.
> > What i notice is Ix And Iy contains a lot of 0s. For example at i=5,j=5, Ix(5,5)=Iy(5,5)=0. Therefore, what puzzles me is how the curvature can be computed, given the below formula:
> > k = {[(Ix^2)(Iyy)] - [2IxIyIxy] + [(Iy^2)(Ixx)]} / {[(Ix^2) + (Iy^2)]^3/2}
> > because now the denominator is 0 (Ix^2+Iy^2).
> > I really need an explanation for this. Please help!
> - - - - - - - - - -
> I have no way of advising you on that because I have not seen the data you are processing. However I find it rather strange that Ix and Iy should both be exactly zero at the same point in the image. Are you sure about that? Perhaps they simply look zero in a display because they are small in comparison to other values being displayed along with them. Try using "format long" to display individual values of Ix and Iy at those same points. I venture to guess they are not exactly zero.
>
> Also bear in mind that because you found it necessary to do filtering with gaussian kernels, the results are bound to only be approximate. You may have to experiment with choosing optimum sigma values and best values for the range of the kernels (what I called 'm') to get the best filtering.
>
> In evaluating the accuracy of these gaussian kernels, I urge you to make that test I described earlier. Create an artificial image of a quadratic form:
>
> I(x,y) = A(x-x0)^2/2+B*(x-x0)*(y-y0)+C*(y-y0)^2/2+D*(x-x0)+E*(y-y0)+F
>
> for some point (x0,y0) and for various arbitrary values of A, B, C, D, E, and F, and see how close Ixx(x0,y0), Ixy(x0,y0), Iyy(x0,y0), Ix(x0,y0), and Iy(x0,y0) are respectively to A, B, C, D, and E. I designed the code I wrote for you to do well on such a test.
>
> Roger Stafford