title: Leveberg-Marquard Optimization for Intensity-based Image Alignment Matlab Failed
body: <p>I am trying to implement the intensity-based image alignment algorithm in the paper <a href="http://ieeexplore.ieee.org/abstract/document/826785/" rel="nofollow noreferrer">Efficient, Robust, and Fast Global Motion Estimation for Video Coding</a> on MATLAB. The problem is that the Levenberg-Marquardt (LM) optimization doesn't work properly, so the energy function cannot approach even a local minima due to very tiny update term computed from the LM approximation equation. I've searching for a week but still could not figure out the problem. Here is my implementation.</p>
  
  <p>I have two gray images <strong>I</strong> and <strong>T</strong> equal in size (W x H) and would like to estimate a rigid transformation <strong>F</strong> (scale + rotation + translation) that aligns <strong>I</strong> onto <strong>T</strong>. This could be obtained by minimizing the following energy function</p>
  
  <p><a href="https://i.stack.imgur.com/dKjNd.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/dKjNd.gif" alt="enter image description here"></a></p>
  
  <p>Where</p>
  
  <p><a href="https://i.stack.imgur.com/avRtH.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/avRtH.gif" alt="enter image description here"></a></p>
  
  <p>According to LM, I can derive the gradient term of E with respect to each parameter <strong>a_i</strong> (i=1,4) as follow (from equation (5) in the paper).</p>
  
  <p><a href="https://i.stack.imgur.com/G7hxx.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/G7hxx.gif" alt="enter image description here"></a></p>
  
  <p>And the terms of the Hessian matrix of <strong>E</strong> are (from equation (4) in the paper)</p>
  
  <p><a href="https://i.stack.imgur.com/4qux8.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4qux8.gif" alt="enter image description here"></a></p>
  
  <p>From the above equations, I derive the bellow MATLAB procedure to compute beta and alpha terms.</p>
  
  <pre><code>%========== 1. Compute gradient of T using Sobel filter
  gx = [-1 -2 -1; 0 0 0; 1 2 1] / 4;
  gy = [-1 0 1; -2 0 2; -1 0 1] / 4;
  
  Tx = conv2( T, gx, 'same' );
  Ty = conv2( T, gy, 'same' );
  
  %========== 2. Warp I using F to compute diff_image = I(x, y) - T(F(x,y,a)). F was previously initialized by an identity matrix of size 3x3
  tform = affine2d( F );
  I_warp = imwarp( I, tform, 'linear', 'OutputView', imref2d( size( T ) ) );
  diff_img = I_wapr - T;
  
  % create a mask for overlaping region between two aligned images
  mask = ones( size( T ) );
  mask_warp = imwarp( mask, tform, 'nearest', 'OutputView', imref2d( size( T ) ) );
  overlap_area = numel( mask_warp ); 
  diff_img = diff_img .* mask;
  error = sum( sum( diff_img ) ) / 2 / overlap_area;
  
  % create x, y grids 
  [y, x] = ndgrid( 0 : h-1, 0 : w-1 );
  x_warp = imwarp( x, t_form, 'nearest', 'OutputView', imref2d( size(T) ) );
  y_warp = imwarp( y, t_form, 'nearest', 'OutputView', imref2d( size(T) ) );
  
  %======== compute beta_i = - dE/da_i (i=1,4)
  sx = Tx .* diff_img;
  sy = Ty .* diff_img;
  
  beta_1 = sum( sum( x_warp .* sx + y_warp .* sy ) );
  beta_2 = sum( sum( y_warp .* sx - x_warp .* sy ) );
  beta_3 = sum( sum( sx ) );
  beta_4 = sum( sum( sy ) );
  beta = -[beta_1; beta_2; beta_3; beta_4] / overlap_area;
  
  %======= compute alpha_ij = (dE/da_i) * (dE/da_j)  i,j = 1,4
  Sxx = (Tx .^ 2) .* mask_warp;
  Syy = (Ty .^ 2) .* mask_warp;
  Sxy = (Tx.* Ty) .* mask_warp;
  
  xx = x_warp .^2;
  yy = y_warp .^2;
  xy = x_warp .* y_warp;
  
  alpha_11 = sum( sum( xx .* Sxx + 2 * xy .* Sxy + yy .* Syy ) );
  alpha_12 = sum( sum( (yy - xx) .* Sxy + xy .* (Sxx - Syy) ) );
  alpha_13 = sum( sum( x .* Sxx + y .* Sxy ) );
  alpha_14 = sum( sum( x .* Sxy + y .* Syy ) );
  alpha_22 = sum( sum( yy .* Sxx - 2 * xy .* Sxy + xx .* Syy ) );
  alpha_23 = sum( sum( y .* Sxx - x .* Sxy ) );
  alpha_24 = sum( sum( y .* Sxy - x .* Syy ) );
  alpha_33 = sum( sum( Sxx ) );
  alpha_34 = sum( sum( Sxy ) );
  alpha_44 = sum( sum( Syy ) );
  
  alpha = [alpha_11 alpha_12 alpha_13 alpha_14;
           alpha_12 alpha_22 alpha_23 alpha_24;
           alpha_13 alpha_23 alpha_33 alpha_34;
           alpha_14 alpha_24 alpha_34 alpha_44] / overlap_area;
  % With lamda was previously initialized by 0.0001
  for i = 1 : 4
      alpha(i, i) = alpha(i, i) * (lamda + 1);
  end
  
  %======== Find the update term: delta_a = alpha^(-1) * beta
  delta_a = pinv( alpha ) * beta
  
  % Or we can solve for delta_a using SVD
  %[U, S, V] = svd( alpha );
  %inv_S = S;
  %for ii = 1 : size(S, 1)
  %    if S(ii, ii)
  %        inv_S(ii, ii) = 1 / S(ii, ii);
  %    end
  %end
  %delta_a = V * inv_S * U' * beta;
  
  %======== Update a_i and new error value
  a = [   F(1, 1)-1;
          F(2, 1);
          F(3, 1);
          F(3, 2)];
  new_a = a + delta_a;
  new_F = [new_a(1)+1   -new_a(2)    0;
            new_a(2)    new_a(1)+1   0;
            new_a(3)     new_a(4)    1];
  
  tform = affine2d( new_F );
  mask_warp = imwarp( mask, tform, 'nearest', 'OutputView', imref2d(size(T)));
  new_I_warp = imwarp( I, tform, 'linear', 'OutputView', imref2d(size(T)));
  new_diff_img = (new_I_warp - T) .* mask_warp;
  new_error = sum( sum( new_diff_img .^2 ) ) / 2/ sum( sum( mask_warp ) );
  
  if( new_error &gt; error )
      lamda = lamda * 10;
  elseif new_error &lt; error
      lamda = lamda /10;
  </code></pre>
  
  <p>The above process is repeated until the iteration reach the limitation or the new_error value less than a predetermined threshold. However, after each iteration, the update terms are too small that parameters make trivial movement. For example,</p>
  
  <pre><code>delta a = 1.0e-04 * [0.0011   -0.0002    0.2186    0.2079]
  </code></pre>
  
  <p>Is there anybody know how to fixed it? Any help would be highly appreciated.</p>