reference :
C. Liu. Beyond Pixels: Exploring New Representations and Applications for Motion Analysis. Doctoral Thesis. Massachusetts Institute of Technology. May 2009.

Algorithm
The core of the algorithm is based on [2] and [3]. The major difference is that I used conjugate gradient for solving large linear systems instead of Gauss-Seidel or SOR. Why conjugate gradient? The solver becomes very simple and easy to code. If in every inner step we need to solve a large linear equation Ax=b, where x is the flow field, then we may find that matrix A can be decomposed to concatenations of filtering and weighting. So in a conjugate gradient solver, we never need to formally write down matrix A, but only need to write a function that consists of filtering and weighting to apply A to x.

In the appendix of my thesis, I also showed how to derive optical flow solver from the angle of iterative reweighted least square (IRLS) instead of Euler-Lagrange. IRLS and Euler-Lagrange end up with the same equations, but IRLS is much easier to understand and simpler to derive.

Although there can be multiple extentions to the optical flow as shown in Appendix A of my thesis, in this package only two-frame optical flow code is attached.

http://people.csail.mit.edu/celiu/OpticalFlow/