Thread Subject:

Mutual information(Image Registration)

Hello everybody
I wrote a code for image registration . In order to optimization I used several norms like sum of squared difference(SSD) or cross-corelation(CC). but the optimization function dose not work with mutual information code anr the reason is (the exit flag),"the grdient of mutual information is 0(zero). it is so strange beacuse it use finite diff to calcuate the grdient and the value of MI function changes. sb can help me please?
thanks

here is my MI function:
function h=registration_error_Mutual_Information(image_1,image_2)
% Takes a pair of images and returns the mutual information Ixy using joint entropy function JOINT_H.m

a=joint_h(image_1,image_2); % calculating joint histogram for two images
[r,c] = size(a);
b= a./sum(a(:)); % normalized joint histogram
y_marg=sum(b,1); %sum of the rows of normalized joint histogram
x_marg=sum(b,2);%sum of columns of normalized joint histogran

Hy=0;
   for i=1:c; % col
      if( y_marg(i)==0 )
         %do nothing
      else
         Hy = Hy + (y_marg(i)*(log2(y_marg(i)))); %marginal entropy for image 1
      end
   end
Hy=-Hy;
Hx=0;
   for i=1:r; %rows
      if( x_marg(i)==0 )
         %do nothing
      else
         Hx = Hx +(x_marg(i)*(log2(x_marg(i)))); %marginal entropy for image 2
      end
   end
   Hx=-Hx;
h_xy = -sum(sum(b.*(log2(b+(b==0))))); % joint entropy

h = (Hx + Hy)/h_xy;% Mutual information



function h=joint_h(image_1,image_2)
% takes a pair of images of equal size and returns the 2d joint histogram.
% used for MI calculation
image_1=round(image_1);image_2=round(image_2);
rows=size(image_1,1);
cols=size(image_1,2);
N1=max(image_1(:))+1;
N2=max(image_2(:))+1;

h=zeros(N1,N2);
%h=zeros(256,256);
for i=1:rows; % col
  for j=1:cols; % rows
    h(image_1(i,j)+1,image_2(i,j)+1)= h(image_1(i,j)+1,image_2(i,j)+1)+1;
  end
end

"eli " <elaheh_kohan@yahoo.com> wrote in message <jrc8g5$m20$1@newscl01ah.mathworks.com>...
> Hello everybody
> I wrote a code for image registration . In order to optimization I used several norms like sum of squared difference(SSD) or cross-corelation(CC). but the optimization function dose not work with mutual information code anr the reason is (the exit flag),"the grdient of mutual information is 0(zero).
============

What optimizer are you using? FMINUNC?
If so, it doesn't sound like it failed. The gradient is supposed to be zero at the solution. Did you try calculating the gradient at the solution returned? Was it not zero?

Or do you mean that it stops at Iteration #0? If so, that usually happens because you have quantization operations like round, floor, ceil, etc... in your objective function. See this recent thread for a very similar discussion:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/320855

Indeed, you appear to have round operations in your joint_h function, which is a likely rendering your function piecewise constant (and hence has zero gradient almost everywhere). . Most people use a differentiable Parzen window to avoid this.

Matt Thank you for your suggestion
My optimization function is "fminlbfgs" from mathwork by 'Dirk-Jan Kroon '.
It works when I did chang the Maximum and Minimum stepsize used for finite difference in order to compute the gradient. Now my problem is it works good sometimes when I change it and sometimes when I use the defualt(with other measures too like SSD and CC).
So do you think I have to make my iamges smooth(by gaussian or averaging)?
It is so shame but my function is so complicate and I do not know how I can compute the grdient of this function.
Thank you very very much for your helps
 



"Matt J" wrote in message <jrcrc0$bic$1@newscl01ah.mathworks.com>...
> "eli " <elaheh_kohan@yahoo.com> wrote in message <jrc8g5$m20$1@newscl01ah.mathworks.com>...
> > Hello everybody
> > I wrote a code for image registration . In order to optimization I used several norms like sum of squared difference(SSD) or cross-corelation(CC). but the optimization function dose not work with mutual information code anr the reason is (the exit flag),"the grdient of mutual information is 0(zero).
> ============
>
> What optimizer are you using? FMINUNC?
> If so, it doesn't sound like it failed. The gradient is supposed to be zero at the solution. Did you try calculating the gradient at the solution returned? Was it not zero?
>
> Or do you mean that it stops at Iteration #0? If so, that usually happens because you have quantization operations like round, floor, ceil, etc... in your objective function. See this recent thread for a very similar discussion:
>
> http://www.mathworks.com/matlabcentral/newsreader/view_thread/320855
>
> Indeed, you appear to have round operations in your joint_h function, which is a likely rendering your function piecewise constant (and hence has zero gradient almost everywhere). . Most people use a differentiable Parzen window to avoid this.

"eli " <elaheh_kohan@yahoo.com> wrote in message <js92fa$g6d$1@newscl01ah.mathworks.com>...
> Matt Thank you for your suggestion
> My optimization function is "fminlbfgs" from mathwork by 'Dirk-Jan Kroon '.
> It works when I did chang the Maximum and Minimum stepsize used for finite difference in order to compute the gradient. Now my problem is it works good sometimes when I change it and sometimes when I use the defualt(with other measures too like SSD and CC).
================

That's an artificial solution. Obviously if you've accidentally made your function piece-wise constant, taking larger steps can help you avoid getting trapped on the constant piece that your initial point is currently sitting on. However, that can be pretty unreliable and also give you a pretty crude approximation to the gradient. Gradient computations are meant to be done with small steps, not large ones.


> So do you think I have to make my iamges smooth(by gaussian or averaging)?
===============

No, as I said before, I think you have to make your histogram calculations differentiable, by using a Parzen window.

Matt
Thanks for your suggestions.
I think I know why I have piece-wise constant function. I do the registration in 3 resolution(Pyramid gussian), so at the coarse resolution (1/4 of the real resolution) the image is going to be more like a piece-wise constant function.
Do you think if in each resolution I use different step-size to calculate the gradient I will have "artificial solution" yet?
Tanx
Ella

"Matt J" wrote in message <js9qr9$pje$1@newscl01ah.mathworks.com>...
> "eli " <elaheh_kohan@yahoo.com> wrote in message <js92fa$g6d$1@newscl01ah.mathworks.com>...
> > Matt Thank you for your suggestion
> > My optimization function is "fminlbfgs" from mathwork by 'Dirk-Jan Kroon '.
> > It works when I did chang the Maximum and Minimum stepsize used for finite difference in order to compute the gradient. Now my problem is it works good sometimes when I change it and sometimes when I use the defualt(with other measures too like SSD and CC).
> ================
>
> That's an artificial solution. Obviously if you've accidentally made your function piece-wise constant, taking larger steps can help you avoid getting trapped on the constant piece that your initial point is currently sitting on. However, that can be pretty unreliable and also give you a pretty crude approximation to the gradient. Gradient computations are meant to be done with small steps, not large ones.
>
>
> > So do you think I have to make my iamges smooth(by gaussian or averaging)?
> ===============
>
> No, as I said before, I think you have to make your histogram calculations differentiable, by using a Parzen window.

"eli " <elaheh_kohan@yahoo.com> wrote in message <jsbsde$3eo$1@newscl01ah.mathworks.com>...
> Matt
> Thanks for your suggestions.
> I think I know why I have piece-wise constant function. I do the registration in 3 resolution(Pyramid gussian), so at the coarse resolution (1/4 of the real resolution) the image is going to be more like a piece-wise constant function.
==============

No, I don't think that would cause piece-wise constant-ness.
I still think it is piece-wise constant because of the rounding operations you have in joint_h. You also haven't shown how you're deforming image1 to image2 and what kind of interpolation you're using. If you're not using a smooth interpolator, you will have non-differentiability problems there as well.



> Do you think if in each resolution I use different step-size to calculate the gradient I will have "artificial solution" yet?
======

I think the need to manually enlarge step-size is a sign that something is wrong/bad in the formulation of the problem.

"Matt J" wrote in message <jsch3b$16f$1@newscl01ah.mathworks.com>...
> "eli " <elaheh_kohan@yahoo.com> wrote in message <jsbsde$3eo$1@newscl01ah.mathworks.com>...
> > Matt
> > Thanks for your suggestions.
> > I think I know why I have piece-wise constant function. I do the registration in 3 resolution(Pyramid gussian), so at the coarse resolution (1/4 of the real resolution) the image is going to be more like a piece-wise constant function.
> ==============
>
> No, I don't think that would cause piece-wise constant-ness.
> I still think it is piece-wise constant because of the rounding operations you have in joint_h. You also haven't shown how you're deforming image1 to image2 and what kind of interpolation you're using. If you're not using a smooth interpolator, you will have non-differentiability problems there as well.
>
>
>
> > Do you think if in each resolution I use different step-size to calculate the gradient I will have "artificial solution" yet?
> ======
>
> I think the need to manually enlarge step-size is a sign that something is wrong/bad in the formulation of the problem.

Hello Mat
I have (sometimes) this problem with other measures like 'ssd' and 'cc' too.(so the problem is not just related to joint_h).
I am using rigid transformation(6 DOF) for transformation, and interp3 with linear method for interpolation, I can not use the interp3 with spline method beacause my data is too big and it returns me "out of memory" error.
do you know another function for spline interpolation?
Tanx alot
Ella

"eli " <elaheh_kohan@yahoo.com> wrote in message <jsefec$lh2$1@newscl01ah.mathworks.com>...
>
> Hello Mat
> I have (sometimes) this problem with other measures like 'ssd' and 'cc' too.(so the problem is not just related to joint_h).
============

I couldn't comment on that without seeing your code, but I would still imagine that somewhere you are quantizing using ROUND, CEIL, FLOOR, etc....

To investigate what's going on, you should make 1D plots your function through your initial point, as a function of one of your 6DOF variables for example. If the plot is locally flat, you know you have quantization issues. If it isn't locally flat, the plot would still give you a graphical insight into why the optimizer thinks it is at a minimum.
 


> I am using rigid transformation(6 DOF) for transformation, and interp3 with linear method for interpolation, I can not use the interp3 with spline method beacause my data is too big and it returns me "out of memory" error.
> do you know another function for spline interpolation?
===============

I can't imagine why splines would give you an out of memory error. The disadvantage I would expect is that it is slower.

In any case, if you have large data, you probably shouldn't be using MATLAB to do registrations anyway. You should be using some high performance registration package like ITK.

"Matt J" wrote in message <jsch3b$16f$1@newscl01ah.mathworks.com>...
> "eli " <elaheh_kohan@yahoo.com> wrote in message <jsbsde$3eo$1@newscl01ah.mathworks.com>...
> > Matt
> > Thanks for your suggestions.
> > I think I know why I have piece-wise constant function. I do the registration in 3 resolution(Pyramid gussian), so at the coarse resolution (1/4 of the real resolution) the image is going to be more like a piece-wise constant function.
> ==============
>
> No, I don't think that would cause piece-wise constant-ness.
> I still think it is piece-wise constant because of the rounding operations you have in joint_h. You also haven't shown how you're deforming image1 to image2 and what kind of interpolation you're using. If you're not using a smooth interpolator, you will have non-differentiability problems there as well.
>
>
>
> > Do you think if in each resolution I use different step-size to calculate the gradient I will have "artificial solution" yet?
> ======
>
> I think the need to manually enlarge step-size is a sign that something is wrong/bad in the formulation of the problem.

"Matt J" wrote in message <jsch3b$16f$1@newscl01ah.mathworks.com>...
> "eli " <elaheh_kohan@yahoo.com> wrote in message <jsbsde$3eo$1@newscl01ah.mathworks.com>...
> > Matt
> > Thanks for your suggestions.
> > I think I know why I have piece-wise constant function. I do the registration in 3 resolution(Pyramid gussian), so at the coarse resolution (1/4 of the real resolution) the image is going to be more like a piece-wise constant function.
> ==============
>
> No, I don't think that would cause piece-wise constant-ness.
> I still think it is piece-wise constant because of the rounding operations you have in joint_h. You also haven't shown how you're deforming image1 to image2 and what kind of interpolation you're using. If you're not using a smooth interpolator, you will have non-differentiability problems there as well.
>
>
>
> > Do you think if in each resolution I use different step-size to calculate the gradient I will have "artificial solution" yet?
> ======
>
> I think the need to manually enlarge step-size is a sign that something is wrong/bad in the formulation of the problem.


Matt
 I did plot my function when my measure is mutual information it is extremely not smooth. I tried to use " Parzen window " but I have problem with it. when we say we have to sample 5% of total pixels in the overlap region it means that I will have 600 gray scales in my hostogram?
so what we mean by 64 bins?
do you know any good article or code(I do not know it good)
Ella