title: My inverse compositional homograhy image align cannot converge?
body: <p>I implement the algorithm based on <a href="https://github.com/albertoCrive/homographyTrackingDemo" rel="nofollow noreferrer">homographyTrackingDemo</a> and <a href="http://www.ncorr.com/download/publications/bakerunify.pdf" rel="nofollow noreferrer">LK-20-years</a>. It's adaption of inverse compositional LK algorithm for estimating affine transformation. There are several steps I have tried: </p>
  
  <ol>
  <li><p>Derive the derivative of template image;</p></li>
  <li><p>Derive the derivative of homography transformation on template control points at identity transformation</p></li>
  <li>Compute the Steepest Gradient Matrix and Hessian matrix</li>
  </ol>
  
  <p>These steps are precomputed before stepping into iteration and we iterate:</p>
  
  <ol>
  <li>Warp target image by homography</li>
  <li>Compute the difference between warped target image and template image</li>
  <li>Solve the linear equations and then update homography parameters</li>
  </ol>
  
  <p>I hope some guy can give me a hint about whether my understanding about this algorithm is wrong. The c++ code is as follows:</p>
  
  <pre><code>AlignmentResults ICA::GaussNewtonMinimization(const StructOfArray2di &amp; pixelsOnTemplate, const vector&lt;Mat&gt; &amp; images, const vector&lt;Mat&gt; &amp; templates, const OptimizationParameters optParam, vector&lt;float&gt; &amp; parameters)
  {
      AlignmentResults alignmentResults;
      alignmentResults.nIter = 0;
      alignmentResults.exitFlag = 1e6;
  
      //parameters must contain the initial guess. It is updated during optimization
      uint nChannels(images.size());
      vector&lt;vector&lt;float&gt; &gt; templatePixelIntensities(nChannels,vector&lt;float&gt;(pixelsOnTemplate.size()));//, templatePixelDx(nChannels,vector&lt;float&gt;(pixelsOnTemplate.size())), templatePixelDy(nChannels,vector&lt;float&gt;(pixelsOnTemplate.size()));
      vector&lt;Mat&gt; imageDx(nChannels), imageDy(nChannels);
      uint nParam = parameters.size();
  
      //vector&lt;Eigen::Matrix&lt;float, 2, N_PARAM&gt;, Eigen::aligned_allocator&lt;Eigen::Matrix&lt;float, 2, N_PARAM&gt; &gt; &gt; warpJacobians(pixelsOnTemplate.size());
  
      Eigen::MatrixXf sdImages(pixelsOnTemplate.size(), nParam);
      vector&lt;float&gt;  errorImage(pixelsOnTemplate.size(), 0.0);
  
      StructOfArray2di warpedPixels;
      Eigen::Matrix&lt;float, N_PARAM, N_PARAM&gt; hessian;
      Eigen::Matrix&lt;float, N_PARAM, N_PARAM&gt; hessianInv;
      hessian.setZero();
  
  
      std::vector&lt;Eigen::Matrix&lt;float, N_PARAM, N_PARAM&gt;, Eigen::aligned_allocator&lt;Eigen::Matrix&lt;float, N_PARAM, N_PARAM&gt;&gt;&gt; hessians(nChannels);
      std::vector&lt;Eigen::MatrixXf&gt; JacosT(nChannels);
      Eigen::Matrix&lt;float, N_PARAM, 1&gt; rhs;
      Eigen::Matrix&lt;float, N_PARAM, 1&gt; deltaParam;
  
      std::vector&lt;float&gt; zeroPara(8, 0.);
  #pragma omp parallel for
      for(int iChannel = 0; iChannel&lt;nChannels; ++iChannel)
      {
          ComputeImageDerivatives(templates[iChannel], imageDx[iChannel], imageDy[iChannel]);
          hessians[iChannel].setZero();
          JacosT[iChannel] = Eigen::MatrixXf(N_PARAM, pixelsOnTemplate.size());
  
          for(int iPoint = 0; iPoint &lt; pixelsOnTemplate.size(); ++iPoint)
          {
              int pos = templates[iChannel].cols * pixelsOnTemplate.y[iPoint] + pixelsOnTemplate.x[iPoint];
              templatePixelIntensities[iChannel][iPoint] = ((float*)templates[iChannel].data)[pos];
  
              float Dx = ((float*)imageDx[iChannel].data)[pos];
              float Dy = ((float*)imageDy[iChannel].data)[pos];
              //printf("%f %f\n", Dx, Dy);
              Eigen::Matrix&lt;float, 2, N_PARAM&gt; warpJ;
              Homography::ComputeWarpJacobian(pixelsOnTemplate.x[iPoint], pixelsOnTemplate.y[iPoint], zeroPara, warpJ);
              //std::cout &lt;&lt; warpJ &lt;&lt; std::endl;
              Eigen::Matrix&lt;float, 1, N_PARAM&gt; J2H;
              Eigen::Matrix&lt;float, 1, 2&gt; ID;
              ID[0] = Dx;
              ID[1] = Dy;
              J2H = ID*warpJ;
  
              JacosT[iChannel].col(iPoint) = J2H.transpose();
              hessians[iChannel] += JacosT[iChannel].col(iPoint) * J2H;
          }
          hessian += hessians[iChannel];
      }
  
      while (alignmentResults.exitFlag == 1e6)
      {
          rhs.setZero();
  
          ComputeWarpedPixels(pixelsOnTemplate, parameters, warpedPixels);
  
  #pragma omp parallel for
          for(int iChannel = 0; iChannel&lt;images.size(); ++iChannel)
          {
              ComputeResiduals(images[iChannel], templatePixelIntensities[iChannel], warpedPixels, errorImage);
              for (int i = 0; i&lt;nParam; ++i)
              {
                  for(uint iPoint(0); iPoint&lt;pixelsOnTemplate.size(); ++iPoint)
                  {
                      float val = (errorImage[iPoint] == std::numeric_limits&lt;float&gt;::infinity() ? 0: errorImage[iPoint]); 
                      rhs(i,0) -= (JacosT[iChannel])(i, iPoint) * val;
                  }
              }
          }
  
  
          deltaParam = hessian.fullPivLu().solve(rhs);
  
          std::vector&lt;float&gt; vDelta(8, 0.), invDelta(8, 0.);
          for(size_t i = 0; i&lt;N_PARAM; ++i)
          {
              vDelta[i] = deltaParamd(i,0);
          }
  
          Homography::InverseWarpParameters(vDelta, invDelta);
          parameters = Homography::ParametersUpdateCompositional(parameters, invDelta);
  
  
          alignmentResults.poseIntermediateGuess.push_back(parameters);
          alignmentResults.residualNorm.push_back(ComputeResidualNorm(errorImage));
          if(alignmentResults.nIter &gt; 0)
              alignmentResults.exitFlag = CheckConvergenceOptimization(deltaParam.norm(), alignmentResults.nIter, abs(alignmentResults.residualNorm[alignmentResults.nIter] - alignmentResults.residualNorm[alignmentResults.nIter-1]), optParam);
          alignmentResults.nIter++;
  
      return alignmentResults;
  }
  </code></pre>