光流法 HS与LK

1、HS算法

#include "cv.h"
#include "highgui.h"
#include <math.h>
#include <stdio.h>
#include "opencv2/legacy/legacy.hpp"
#define CVX_GRAY50 cvScalar(100)
#define CVX_WHITE  cvScalar(255)

int main(int argc, char** argv)
{
    // Initialize, load two images from the file system, and
    // allocate the images and other structures we will need for
    // results.

    // exit if no input images
	 IplImage *imgA = 0, *imgB = 0;
    //imgA = cvLoadImage("OpticalFlow0.jpg",0);
    //imgB = cvLoadImage("OpticalFlow1.jpg",0);
    //if(!(imgA)||!(imgB)){ printf("One of OpticalFlow0.jpg and/or OpticalFlow1.jpg didn't load\n"); return -1;}
	imgA = cvLoadImage("OpticalFlow0.jpg",0);
    imgB = cvLoadImage("OpticalFlow1.jpg",0);
    if(!(imgA)||!(imgB)){ printf("One of OpticalFlow0.jpg and/or OpticalFlow1.jpg didn't load\n"); return -1;}

    IplImage* velx = cvCreateImage(cvGetSize(imgA),IPL_DEPTH_32F,1);
    IplImage* vely = cvCreateImage(cvGetSize(imgA),IPL_DEPTH_32F,1);

    IplImage* imgC = cvCreateImage(cvGetSize(imgA),IPL_DEPTH_8U,3);

    cvNamedWindow( "OpticalFlow0" );
    cvNamedWindow( "OpticalFlow1" );
    cvNamedWindow( "Flow Results" );

    cvShowImage( "OpticalFlow0",imgA );
    cvShowImage( "OpticalFlow1",imgB );

    // Call the actual Horn and Schunck algorithm
    //
    cvCalcOpticalFlowHS( 
        imgA, 
        imgB, 
        0,
        velx,
        vely,
        .10,
        cvTermCriteria( 
            CV_TERMCRIT_ITER | CV_TERMCRIT_EPS,
            imgA->width,
            1e-6
        )
    );

    // Now make some image of what we are looking at:
    //
    cvZero( imgC );
    int step = 4;
    for( int y=0; y<imgC->height; y += step ) {
        float* px = (float*) ( velx->imageData + y * velx->widthStep );
        float* py = (float*) ( vely->imageData + y * vely->widthStep );
        for( int x=0; x<imgC->width; x += step ) {
            if( px[x]>1 && py[x]>1 ) {
                cvCircle(
                    imgC,
                    cvPoint( x, y ),
                    2,
                    CVX_GRAY50,
                    -1
                );
                cvLine(
                    imgC,
                    cvPoint( x, y ),
                    cvPoint( x+px[x]/2, y+py[x]/2 ),
                    CV_RGB(255,0,0),
                    1,
                    0
                );
            }
        }
    }
    // show tracking
    cvShowImage( "Flow Results",imgC );
    
    cvWaitKey(0);

    // destroy windows
    cvDestroyWindow( "OpticalFlow0" );
    cvDestroyWindow( "OpticalFlow1" );
    cvDestroyWindow( "Flow Results" );
    // release memory
    cvReleaseImage( &imgA );
    cvReleaseImage( &imgB );
    cvReleaseImage( &imgC );

    return 0;
}

2、LK算法

#include "cv.h"
#include "cxcore.h"
#include "highgui.h"
#include "stdio.h"

const int MAX_CORNERS = 500;
int main(int argc, char** argv) {
   // Initialize, load two images from the file system, and
   // allocate the images and other structures we will need for
   // results.
	//IplImage* imgA = cvLoadImage("OpticalFlow0.jpg",CV_LOAD_IMAGE_GRAYSCALE);
	//IplImage* imgB = cvLoadImage("OpticalFlow1.jpg",CV_LOAD_IMAGE_GRAYSCALE);
	IplImage* imgA = cvLoadImage("3.jpg",CV_LOAD_IMAGE_GRAYSCALE);
	IplImage* imgB = cvLoadImage("4.jpg",CV_LOAD_IMAGE_GRAYSCALE);
	//IplImage* imgA = cvLoadImage("10.png",CV_LOAD_IMAGE_GRAYSCALE);
	//IplImage* imgB = cvLoadImage("11.png",CV_LOAD_IMAGE_GRAYSCALE);
	CvSize      img_sz    = cvGetSize( imgA );
	int         win_size = 10;
	//IplImage* imgC = cvLoadImage("OpticalFlow1.jpg",CV_LOAD_IMAGE_UNCHANGED);
	IplImage* imgC = cvLoadImage("4.jpg",CV_LOAD_IMAGE_UNCHANGED);
	//IplImage* imgC = cvLoadImage("11.png",CV_LOAD_IMAGE_UNCHANGED);
	
	// The first thing we need to do is get the features
	// we want to track.
	//
	IplImage* eig_image = cvCreateImage( img_sz, IPL_DEPTH_32F, 3 );
	IplImage* tmp_image = cvCreateImage( img_sz, IPL_DEPTH_32F, 3 );
	int              corner_count = MAX_CORNERS;
	CvPoint2D32f* cornersA        = new CvPoint2D32f[ MAX_CORNERS ];
	cvGoodFeaturesToTrack(
		imgA,
		eig_image,
		tmp_image,
		cornersA,
		&corner_count,
		0.01,
		5.0,
		0,
		3,
		0,
		0.04
	);
	cvFindCornerSubPix(
		imgA,
		cornersA,
		corner_count,
		cvSize(win_size,win_size),
		cvSize(-1,-1),
		cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,20,0.03)
	);
	// Call the Lucas Kanade algorithm
	//
	char features_found[ MAX_CORNERS ];
	float feature_errors[ MAX_CORNERS ];
	CvSize pyr_sz = cvSize( imgA->width+8, imgB->height/3 );
	IplImage* pyrA = cvCreateImage( pyr_sz, IPL_DEPTH_32F,3 );
    IplImage* pyrB = cvCreateImage( pyr_sz, IPL_DEPTH_32F,3 );
  CvPoint2D32f* cornersB        = new CvPoint2D32f[ MAX_CORNERS ];
  cvCalcOpticalFlowPyrLK(
     imgA,
     imgB,
     pyrA,
     pyrB,
     cornersA,
     cornersB,
     corner_count,
     cvSize( win_size,win_size ),
     5,
     features_found,
     feature_errors,
     cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3 ),
     0
  );
  // Now make some image of what we are looking at:
  //
  for( int i=0; i<corner_count; i++ ) {
     if( features_found[i]==0|| feature_errors[i]>550 ) {
 //       printf("Error is %f/n",feature_errors[i]);
        continue;
     }
 //    printf("Got it/n");
     CvPoint p0 = cvPoint(
        cvRound( cornersA[i].x ),
        cvRound( cornersA[i].y )
     );
     CvPoint p1 = cvPoint(
        cvRound( cornersB[i].x ),
        cvRound( cornersB[i].y )
     );
     cvLine( imgC, p0, p1, CV_RGB(255,0,0),2 );
  }
  cvNamedWindow("ImageA",0);
  cvNamedWindow("ImageB",0);
  cvNamedWindow("LKpyr_OpticalFlow",0);
  cvShowImage("ImageA",imgA);
  cvShowImage("ImageB",imgB);
  cvShowImage("LKpyr_OpticalFlow",imgC);
  cvWaitKey(0);
  return 0;
}