title: iris to screen calculation for eye tracking
body: <p>I'm currently experimenting with Eye tracking I've successfully built an iris tracking algorithm using OpenCV with contours and Hough transform. But the next step is unclear for me. I want to know if the calculations i'm doing are correct for translating the center of an eye to the screen. The head of the user has an fixed position.</p>
  
  <p><a href="https://i.stack.imgur.com/eMxLm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/eMxLm.png" alt="result"></a></p>
  
  <p>What I want is an algorithm that works on all eyes off course. Is there like an angle calculation? So when the user is looking more to the right, linear?   </p>
  
  <p>What I do right now is:
   First I let the user look at specific points and use RANSAC to detect the iris position that's closest to the position on the screen. I do that with four 2D points on the screen and iris. I'm using Homography for this to get the correct calculation. </p>
  
  <pre><code>void gaussian_elimination(float *input, int n){
  // ported to c from pseudocode in
  // http://en.wikipedia.org/wiki/Gaussian_elimination
  
  float * A = input;
  int i = 0;
  int j = 0;
  int m = n-1;
  while (i &lt; m &amp;&amp; j &lt; n){
      // Find pivot in column j, starting in row i:
      int maxi = i;
      for(int k = i+1; k&lt;m; k++){
          if(fabs(A[k*n+j]) &gt; fabs(A[maxi*n+j])){
              maxi = k;
          }
      }
      if (A[maxi*n+j] != 0){
          //swap rows i and maxi, but do not change the value of i
          if(i!=maxi)
              for(int k=0;k&lt;n;k++){
                  float aux = A[i*n+k];
                  A[i*n+k]=A[maxi*n+k];
                  A[maxi*n+k]=aux;
              }
          //Now A[i,j] will contain the old value of A[maxi,j].
          //divide each entry in row i by A[i,j]
          float A_ij=A[i*n+j];
          for(int k=0;k&lt;n;k++){
              A[i*n+k]/=A_ij;
          }
          //Now A[i,j] will have the value 1.
          for(int u = i+1; u&lt; m; u++){
              //subtract A[u,j] * row i from row u
              float A_uj = A[u*n+j];
              for(int k=0;k&lt;n;k++){
                  A[u*n+k]-=A_uj*A[i*n+k];
              }
              //Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
          }
  
          i++;
      }
      j++;
  }
  
  //back substitution
  for(int i=m-2;i&gt;=0;i--){
      for(int j=i+1;j&lt;n-1;j++){
          A[i*n+m]-=A[i*n+j]*A[j*n+m];
          //A[i*n+j]=0;
      }
  }
  }
  
  
  
  ofMatrix4x4 findHomography(ofPoint src[4], ofPoint dst[4]){
  ofMatrix4x4 matrix;
  
  // create the equation system to be solved
  //
  // from: Multiple View Geometry in Computer Vision 2ed
  //       Hartley R. and Zisserman A.
  //
  // x' = xH
  // where H is the homography: a 3 by 3 matrix
  // that transformed to inhomogeneous coordinates for each point
  // gives the following equations for each point:
  //
  // x' * (h31*x + h32*y + h33) = h11*x + h12*y + h13
  // y' * (h31*x + h32*y + h33) = h21*x + h22*y + h23
  //
  // as the homography is scale independent we can let h33 be 1 (indeed any of the terms)
  // so for 4 points we have 8 equations for 8 terms to solve: h11 - h32
  // after ordering the terms it gives the following matrix
  // that can be solved with gaussian elimination:
  
  float P[8][9]={
      {-src[0].x, -src[0].y, -1,   0,   0,  0, src[0].x*dst[0].x, src[0].y*dst[0].x, -dst[0].x }, // h11
      {  0,   0,  0, -src[0].x, -src[0].y, -1, src[0].x*dst[0].y, src[0].y*dst[0].y, -dst[0].y }, // h12
  
      {-src[1].x, -src[1].y, -1,   0,   0,  0, src[1].x*dst[1].x, src[1].y*dst[1].x, -dst[1].x }, // h13
      {  0,   0,  0, -src[1].x, -src[1].y, -1, src[1].x*dst[1].y, src[1].y*dst[1].y, -dst[1].y }, // h21
  
      {-src[2].x, -src[2].y, -1,   0,   0,  0, src[2].x*dst[2].x, src[2].y*dst[2].x, -dst[2].x }, // h22
      {  0,   0,  0, -src[2].x, -src[2].y, -1, src[2].x*dst[2].y, src[2].y*dst[2].y, -dst[2].y }, // h23
  
      {-src[3].x, -src[3].y, -1,   0,   0,  0, src[3].x*dst[3].x, src[3].y*dst[3].x, -dst[3].x }, // h31
      {  0,   0,  0, -src[3].x, -src[3].y, -1, src[3].x*dst[3].y, src[3].y*dst[3].y, -dst[3].y }, // h32
  };
  
  gaussian_elimination(&amp;P[0][0],9);
  
  matrix(0,0)=P[0][8];
  matrix(0,1)=P[1][8];
  matrix(0,2)=0;
  matrix(0,3)=P[2][8];
  
  matrix(1,0)=P[3][8];
  matrix(1,1)=P[4][8];
  matrix(1,2)=0;
  matrix(1,3)=P[5][8];
  
  matrix(2,0)=0;
  matrix(2,1)=0;
  matrix(2,2)=0;
  matrix(2,3)=0;
  
  matrix(3,0)=P[6][8];
  matrix(3,1)=P[7][8];
  matrix(3,2)=0;
  matrix(3,3)=1;
  
  return matrix;
  </code></pre>
  
  <p>}</p>