Photometric stereo 2009-06-07

To do: Surface shape reconstruction
: Least squares surface fitting

ref.
http://www.cs.washington.edu/education/courses/csep576/05wi/projects/project3/project3.htm

Forsyth and Ponce (84-86p) : 5.4.2 Shape from Normals : Shape by Integration

http://en.wikipedia.org/wiki/Lucas%E2%80%93Kanade_Optical_Flow_Method

- path integration
- least squares optimization
- Lucas-Kanade algorithm

http://mathworld.wolfram.com/LeastSquaresFitting.html

 /* We want to solve Mz = v in a least squares sense.  The
    solution is M^T M z = M^T v.  We denote M^T M as A and
    M^T v as b, so A z = b. */

 CMatrixSparse<double> A(M.mTm());
    assert(A.isSymmetric());
    CVector<double> r = A*z;  /* r is the "residual error" */
    CVector<double> b(v*M);

 // solve the equation A z = b
    solveQuadratic<double>(A,b,z,300,CGEPSILON);

 // copy the depths back from the vector z into the image depths
    copyDepths(z,zind,depths);

 

template <class T>
double solveQuadratic(const CMatrixSparse<T> & A, const CVector<T> & b,
       CVector<T> & x,int i_max, double epsilon)
{
  //my conjugate gradient solver for .5*x'*A*x -b'*x, based on the
  // tutorial by Jonathan Shewchuk  (or is it +b'*x?)
 
  printf("Performing conjugate gradient optimization\n");

  int numvars = x.Length();
  assert(b.Length() == numvars && A.rows() == numvars &&
  A.columns() == numvars);

  int i =0;

  CVector<T> r = b-A*x;
  CVector<T> d= r;
  double delta_new = r.dot(r);
  double delta_0 = delta_new;

  int numrecompute = (int)floor(sqrt(float(numvars)));
  //int numrecompute = 1;
  printf("numrecompute = %d\n",numrecompute);

  printf("delta_new = %f\n", delta_new);
  while (i < i_max && delta_new > epsilon)//epsilon*epsilon*delta_0)
    {
      printf("Step %d, delta_new = %f      \r",i,delta_new);
     
      CVector<T> q = A*d;
      double alpha = delta_new/d.dot(q);
      x.daxpy(d,alpha); //      x += d*alpha;
      if (i % numrecompute == 0)
 {
   //   printf(" ** recompute\n");
   r = b-A*x;
 }
      else
 r.daxpy(q,-alpha); //  r = r-q*alpha;
      double delta_old = delta_new;
      delta_new = r.dot(r);
      double beta = delta_new/delta_old;
      d = r+d*beta;
      i++;
    }

  return delta_new;
  //  return delta_new <= epsilon;
  //  return !(delta_new > epsilon*epsilon*delta_0);
}

surfNorm.cpp

lights.txt