title: Image registration and focus stacking
body: <p><strong>Background:</strong>
  I am looking to align images for a focus stacking application using a smartphone.  Links to images:</p>
  
  <p>First in stack: <a href="https://i.stack.imgur.com/QEFjJ.jpg" rel="nofollow noreferrer">1</a>, Last in stack: <a href="https://i.stack.imgur.com/nZJbo.jpg" rel="nofollow noreferrer">2</a>, Final stacked images: <a href="https://i.stack.imgur.com/qm1bv.png" rel="nofollow noreferrer">3</a></p>
  
  <p>I.e. images are nominally the same, BUT contain:</p>
  
  <ol>
  <li>Systematic change in FOCUS as the focal plane shifts between images</li>
  <li>Magnification changes slightly (smartphone feature as focus changes!)</li>
  <li>Camera moves slightly due to random vibrations.</li>
  </ol>
  
  <p>Images need to be aligned for the focus-stacking APP to work.    </p>
  
  <p><strong>Progress to date:</strong>
  I use OpenCV's findTransformECC() to get alignment.  It works well after some experimentation i.e. see <a href="https://stackoverflow.com/questions/45997891/cv2-motion-euclidean-for-the-warp-mode-in-ecc-image-alignment-method">cv2.MOTION_EUCLIDEAN for the warp_mode in ECC image alignment method</a> which was useful to improve the initialization of the Warp matrix: </p>
  
  <ol>
  <li>Images aligned at pixel level</li>
  <li>60secs to process 8Mpix image (1sec for 0.5Mpix image) (on 3 year old portable PC with OpenCV release libraries)</li>
  </ol>
  
  <p>See stacked image link above.</p>
  
  <p>I briefly investigated a feature detector (SIFT).  It did not align the images well, presumably due to the change in focus between images.</p>
  
  <p><strong>Code:</strong></p>
  
  <pre><code> int scale = 1;
   int scaleSmall = 4;
   float scaleDiff = scaleSmall / scale;
  
       for (i = 0; i&lt; numImages; i++) {
          file = dir + image + to_string(i) + ".jpg";
          col[i] = imread(file);
  
          resize(col[i], z[i], Size(col[i].cols/scale, col[i].rows/scale));       
          cvtColor(z[i], zg[i], CV_BGR2GRAY);
          resize(zg[i], zgSmall[i],  Size(col[i].cols / scaleSmall, col[i].rows / scaleSmall));
       }
  
      // Set a 2x3 or 3x3 warp matrix depending on the motion model.
      // See https://www.learnopencv.com/image-alignment-ecc-in-opencv-c-python/
      // Define the motion model
      const int warp_mode = MOTION_HOMOGRAPHY;
  
      // Initialize the matrix to identity
      if (warp_mode == MOTION_HOMOGRAPHY) {
          warp_init = Mat::eye(3, 3, CV_32F);
          warp_matrix = Mat::eye(3, 3, CV_32F);
          warp_matrix_prev = Mat::eye(3, 3, CV_32F);
          scaleTX = (Mat_&lt;float&gt;(3, 3) &lt;&lt; 1, 1, scaleDiff, 1, 1, scaleDiff, 1 / scaleDiff, 1 / scaleDiff, 1);
      }
      else {
          warp_init = Mat::eye(2, 3, CV_32F);
          scaleTX = Mat::eye(2, 3, CV_32F);
          warp_matrix = Mat::eye(2, 3, CV_32F);
          warp_matrix_prev = Mat::eye(2, 3, CV_32F);
          scaleTX = (Mat_&lt;float&gt;(2, 3) &lt;&lt; 1, 1, scaleDiff, 1, 1, scaleDiff);
      }
  
      // Specify the number of iterations.
      int number_of_iterations = 5000;
  
      // Specify the threshold of the increment
      // in the correlation coefficient between two iterations
      double termination_eps = 1e-8;
  
      // Define termination criteria
      TermCriteria criteria(TermCriteria::COUNT + TermCriteria::EPS, number_of_iterations, termination_eps);
  
      for (i = 1; i &lt; numImages; i++) {
          // Check images right size
          if (zg[0].rows &lt; 10 || zg[1].rows &lt; 10)
              return;
  
          // Run the ECC algorithm at start to get an initial guess. The results are stored in warp_matrix.
          if (i == 1) {
              findTransformECC(zgSmall[0], zgSmall[i], warp_init, warp_mode, criteria     );
  
              // See https://stackoverflow.com/questions/45997891/cv2-motion-euclidean-for-the-warp-mode-in-ecc-image-alignment-method
              warp_matrix = warp_init * scaleTX;
          }
  
          // Warp Matrix from previous iteration is used as initialisation  
          findTransformECC(zg[0], zg[i], warp_matrix, warp_mode,  criteria);
  
          if (warp_mode != MOTION_HOMOGRAPHY) {
              warpAffine(zg[i], ag[i], warp_matrix, zg[i].size(), INTER_LINEAR + WARP_INVERSE_MAP);
              warpAffine(z[i], acol[i], warp_matrix, zg[i].size(), INTER_LINEAR + WARP_INVERSE_MAP);
          }
          else {
              // Use warpPerspective for Homography
              warpPerspective(z[i], acol[i], warp_matrix, z[i].size(), INTER_LINEAR + WARP_INVERSE_MAP);
              warpPerspective(zg[i], ag[i], warp_matrix, zg[i].size(), INTER_LINEAR + WARP_INVERSE_MAP);
             }
          }
      }
  </code></pre>
  
  <p><strong>Question:</strong>
  Can the image registration speed be significantly improved (using the same hardware)?</p>