//matrix faces: a column is a "face"-image represented in greyscale-pixels (rows=numFaces)

  pcabuildbasis(faces, numPixels, numFaces, info, coeffs, eigenFaces);

  // reproduce faces by scalar-multiplying coeffs with eigenFaces and adding mean:
  for (int face = 0; face<numFaces; ++face) {
    for (int pixel = 0; pixel<numPixels; ++pixel) {
      double reconstruct = 0.0;
      // reconstruct pixel using all basisfunctions
      for (int eigenface = 0; eigenface<numPixels; ++eigenface) {         
        reconstruct += eigenFaces[pixel][eigenface] * coeffs[eigenface];
      }
      reconstruct += means[face];
      matrix[pixel][face] = reconstruct;     
    }
  }
  // Why is matrix!=faces ?

// here I assume that all elements were correctly initialized
pcabuildbasis(faces, numPixels, numFaces, info, coeffs, eigenFaces);
for (int face_to_reconstruct = 0; face_to_reconstruct<numFaces; ++face_to_reconstruct)
{
    for(pixel = 0; pixel<numPixels; ++pixel)
        reconstructed[pixel] = 0.0;
    for (int eigenface = 0; eigenface<numPixels; ++eigenface)
    {
        // coefficient before eigenface
        double coefficient = 0.0;

        // calculate coefficient
        for (int pixel = 0; pixel<numPixels; ++pixel)
            coefficient += eigenFaces[pixel][eigenface]*faces[face_to_reconstruct];

        // update reconstructed face
        for (int pixel = 0; pixel<numPixels; ++pixel)
            reconstructed[pixel] += eigenFaces[pixel][eigenface]*coefficient;
    }
}