图像笔记（二）——特征提取及跟踪匹配

谢欣燕

于 2021-05-07 17:33:03 发布

阅读量2.3k

点赞数 4

分类专栏：笔记

本文链接：https://blog.csdn.net/weixin_43941538/article/details/116484851

版权

笔记专栏收录该内容

70 篇文章 6 订阅

订阅专栏

2.图像特征提取及匹配

2.1 基于灰度模板匹配算法

Blocking Matching 是根据移植模板图到另一幅图像中寻找与模板图像相似的自图像。基于灰度图像的匹配算法也称作相关匹配算法，用空间二维滑动模板进行匹配，不同的匹配算法相关准则不同。

2.1.1 平均绝对差算法（MAD）

思路简单
匹配精度较高
$D(i,j)=\frac{1}{M*N}\sum^m_{S=1}\sum^N_{t=1}|S(i+s-1,j+t-1)-T(s,t)| \\ m-M+1≥i≥1，n-N+1≥j≥1$

void MAD(Mat src,Mat in,int & col ,int & row,int step)
{
    int srow=src.rows;
    int scol=src.cols;
    int irow=in.rows;
    int icol=in.cols;
    assert(irow<=srow);
    assert(icol<=scol);
    Mat tmp_mat=Mat(Size(irow,icol),CV_8UC1);
    Rect rect;
    Point point1(0,0);
    Point point2(irow,icol);
    double minMean(1000);
    int finali,finalj;
    for(int i=0;i<srow-irow+1;i+=step)
    {
         for(int j=0;j<scol-icol+1;j+=step)
         {
             point1=Point(j,i);
             point2=Point(j+icol,i+irow);
             rect=Rect(point1,point2);
             rect=rect & Rect(Point(0,0),Point(scol,srow));
             tmp_mat=src(rect);
             if(tmp_mat.rows==irow && tmp_mat.cols==icol)
             {
                 Mat diff_image;
                 absdiff(tmp_mat,in,diff_image);
                 Mat mat_mean,mat_stddev;
                 meanStdDev(diff_image,mat_mean,mat_stddev);
                 double m,s;
                 m=mat_mean.at<double>(0,0);
                 s=mat_stddev.at<double>(0,0);
                 if(m<minMean)
                 {
                     minMean=m;
                     finali=i,finalj=j;
                 }
             }

         }
    }
    col=finalj;
    row=finali;
}

2.1.2 绝对误差和算法（SAD）

$D(i,j)=\sum^m_{S=1}\sum^N_{t=1}|S(i+s-1,j+t-1)-T(s,t)| \\ m-M+1≥i≥1，n-N+1≥j≥1$

2.1.3 误差平方和算法（SSD）

$D(i,j)=\frac{1}{M*N}\sum^m_{S=1}\sum^N_{t=1}|S(i+s-1,j+t-1)-T(s,t)|^2 \\ m-M+1≥i≥1，n-N+1≥j≥1$

2.1.4 平均误差平方和算法（MSD）

$D(i,j)=\frac{1}{M*N}\sum^m_{S=1}\sum^N_{t=1}|S(i+s-1,j+t-1)-T(s,t)|^2 \\ m-M+1≥i≥1，n-N+1≥j≥1$

2.1.5 归一化积算法（NCC，normalized cross correlation）

对原始图像内任意一个像素点（Px,Py）构建一个n*n的邻域作为匹配窗口。然后对于目标像素位置(px+d,py)同样构建一个n*n大小的匹配串口，对两个窗口进行相似度度量，注意这里的 d d d有一个取值范围。对于两幅图像来说，在进行 N C C NCC NCC计算之前要对图像处理，也就是将两帧图像校正到水平位置，即光心处于同一水平线上，此时极线是水平的，否则匹配过程只能在倾斜的极线方向上完成，这将消耗更多的计算资源。

$NCC(p,d)=\frac{\sum_{(x,y)\in W_p }(I_1(x,y)-\bar I_1(p_x,p_y))(I_2(x+d,y)-\bar I_2(p_x+d,p_y))}{\sqrt{\sum_{(x,y)\in W_p }(I_1(x,y)-\bar I_1(p_x,p_y))^2(I_2(x+d,y)-\bar I_2(p_x+d,p_y))^2}}$

$\\ Wp为之前提到的匹配窗口。 \\ I_1(x,y)为原始图像的像素值。 \\ \bar I_1(p_x,p_y)为原始窗口内像素的均值。 \\ I_2(x_d,y)为原始图像在目标图像上对应点位置在x方向上偏移d后的像素值。 \\ \bar I_2( p_x+d,p_y)为目标图像匹配窗口像素均值。$

$r(X,Y)=\frac{cov(X,Y)}{\sqrt {Var(X)Var(Y)}} \\ cov(X,Y)=E[(x-E(x))(y-E(y))]=E(xy)-E(x)E(y)$

若 NCC=-1，则表示两个匹配窗口完全不相关，相反，若NCC=1时，表示两个匹配窗口相关程度非常高。

# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
import cv2
from numpy import *
from numpy.ma import array
from scipy.ndimage import filters
def plane_sweep_ncc(im_l,im_r,start,steps,wid):
    """ 使用归一化的互相关计算视差图像 """
    m,n = im_l.shape
    # 保存不同求和值的数组
    mean_l = zeros((m,n))
    mean_r = zeros((m,n))
    s = zeros((m,n))
    s_l = zeros((m,n))
    s_r = zeros((m,n))
    # 保存深度平面的数组
    dmaps = zeros((m,n,steps))
    # 计算图像块的平均值
    filters.uniform_filter(im_l,wid,mean_l)
    filters.uniform_filter(im_r,wid,mean_r)
    # 归一化图像
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    # 尝试不同的视差
    for displ in range(steps):
        # 将左边图像移动到右边，计算加和
        filters.uniform_filter(np.roll(norm_l, -displ - start) * norm_r, wid, s) # 和归一化
        filters.uniform_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, s_l)
        filters.uniform_filter(norm_r*norm_r,wid,s_r) # 和反归一化
        # 保存 ncc 的分数
        dmaps[:,:,displ] = s / sqrt(s_l * s_r)
        # 为每个像素选取最佳深度
    return np.argmax(dmaps, axis=2)

def plane_sweep_gauss(im_l,im_r,start,steps,wid):
 """ 使用带有高斯加权周边的归一化互相关计算视差图像 """
 m,n = im_l.shape
 # 保存不同加和的数组
 mean_l = zeros((m,n))
 mean_r = zeros((m,n))
 s = zeros((m,n))
 s_l = zeros((m,n))
 s_r = zeros((m,n))
 # 保存深度平面的数组
 dmaps = zeros((m,n,steps))
 # 计算平均值
 filters.gaussian_filter(im_l,wid,0,mean_l)
 filters.gaussian_filter(im_r,wid,0,mean_r)
 # 归一化图像
 norm_l = im_l - mean_l
 norm_r = im_r - mean_r
 # 尝试不同的视差
 for displ in range(steps):
     # 将左边图像移动到右边，计算加和
     filters.gaussian_filter(np.roll(norm_l, -displ - start) * norm_r, wid, 0, s) # 和归一化
     filters.gaussian_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, 0, s_l)
     filters.gaussian_filter(norm_r*norm_r,wid,0,s_r) # 和反归一化
     # 保存 ncc 的分数
     dmaps[:,:,displ] = s / np.sqrt(s_l * s_r)
 # 为每个像素选取最佳深度
 return np.argmax(dmaps, axis=2)
im_l = array(Image.open(r'1.png').convert('L'), 'f')
im_r = array(Image.open(r'2.png').convert('L'),'f')

# 开始偏移，并设置步长
steps = 12
start = 4
# ncc 的宽度
wid = 10
res = plane_sweep_gauss(im_l,im_r,start,steps,wid)
import scipy.misc
scipy.misc.imsave('depth1.png',res)
show()

2.1.6 序贯相似性（SSDA，normalized cross correlation

定义绝对误差

$\epsilon(i,j,s,t)=|S_{i,j}(s,t)-\bar{S_{i,j}}-T(s,t)+\bar T| \\ \bar{S_{i,j}}=E(S_{i,j})=\frac{1}{M*N}\sum^M_{s=1}\sum^N_{t=1}S_{i,j}(s,t) \\ \bar{T}=E(T)=\frac{1}{M*N}\sum^M_{s=1}\sum^N_{t=1}T(s,t)$

设定阈值Th
在模板图中随机选取不重复的像素点，计算与子图的绝对误差，将误差累加，当误差累加超过了Th时，记下累加次数H，所有子图的累加次数H用一个表R(i,j)表示。
$R(i,j)=\{H|min_{1≤H≤M*N}[\sum^H_{h=1}\epsilon(i,j,s,t)≥Th]|\}$
找出累加次数最多的R，表示为最相近的区域。
为了提高速度，可以线进行粗匹配，隔行选取子图，定位到子图后在将子图分割，提高计算速度。

void SSDA(Mat src,Mat in,int & col, int& row ,float th)
{
    int srow=src.rows;
    int scol=src.cols;
    int irow=in.rows;
    int icol=in.cols;
    assert(irow<=srow);
    assert(icol<=scol);
    Mat tmp_mat=Mat(Size(irow,icol),CV_8UC1);
    Rect rect;
    Point point1(0,0);
    Point point2(irow,icol);
    double minMean(1000);
    int finali,finalj;
    Mat mat_mean,mat_stddev;
    meanStdDev(in,mat_mean,mat_stddev);
    double m_in;
    m_in=mat_mean.at<double>(0,0);
    int R[scol-icol][srow-irow];
    for(int i=0;i<scol-icol;i++)
        for(int j=0;j<srow-irow;j++)
            R[i][j]=0;
    for(int i=0;i<srow-irow;i+=10)
    {
         for(int j=0;j<scol-icol;j+=10)
         {
             point1=Point(j,i);
             point2=Point(j+icol,i+irow);
             rect=Rect(point1,point2);
             rect=rect & Rect(Point(0,0),Point(scol,srow));
             tmp_mat=src(rect);
//             imshow("tmp_mat",tmp_mat);
//             waitKey(1);
             meanStdDev(tmp_mat,mat_mean,mat_stddev);
             double m_temp;
             m_temp=mat_mean.at<double>(0,0);
             double r_ij(0.0);
             //std::cout<<"i is "<<i<<"j is "<<j<<std::endl;
             int h=0;
             for(int c=0;c<icol;)
                 for(int r=0;r<irow;)
                 {
                     h++;
                     r_ij+=fabs((float)tmp_mat.at<uchar>(c,r)-m_temp-(float)in.at<uchar>(c,r)+m_in);
                     c++;
                     r++;
                     if(r_ij>th)
                     {
                         R[i][j]=h;
                         c=icol;
                         r=irow;
                     }
                 }
         }
    }
    int h_max=0;
    for(int i=0;i<scol-icol;i+=10)
    {
        for(int j=0;j<srow-irow;j+=10)
        {
            std::cout<<R[i][j]<<"  ";
            if(R[i][j]>h_max)
            {
                h_max=R[i][j];

                col=j;
                row=i;
            }
        }
        std::cout<<std::endl;
    } std::cout<<col<<"  "<<row<<"  "<<h_max<<std::endl;
}

2.1.7 hadamard变换算法（SATD）

hadamarddiff_imagehadamard

2.1.8 CENSUS变换

census变换保留了窗口中像素的位置特征，并且对亮度偏差较为鲁棒，可以减少光照差异引起的误匹配。
步骤：

census变换是使用像素邻域内的局部灰度差异将像素灰度转化为比特串，通过将邻域窗口（窗口大小为n*m，n和m都是奇数）内的像素灰度值与窗口中心像素的灰度值进行比较，将比较得到的布尔值映射到一个比特串中，最后用比特串的值作为中心像素的census变换值Cs。
对于欲求取视差的左右图，要比较两个视图㕜两点的相似度，可将此两点的census值逐位进行异或运算，然后计算结果为1的个数，记为两点之间的汉明值，汉明值是两点间相似度的一种体现，汉明越小，两点相似度越大。

Mat CensusTransform(Mat src,int window_sizex,int window_sizey)
{
    int image_height=src.rows;
    int image_width=src.cols;
    Mat dst=Mat::zeros(image_height,image_width,CV_64F);
    //census transform
    int offsetx=(window_sizex-1)/2;
    int offsety=(window_sizey-1)/2;
    for(int j=0;j<image_width-window_sizex;j++)
    {
        for(int i=0;i<image_height-window_sizey;i++)
        {
            unsigned long census=0;
            uchar current_pixel=src.at<uchar>(i+offsety,j+offsetx);
            Rect roi(j,i,window_sizex,window_sizey);
            Mat window(src,roi);
            for(int a = 0; a <window_sizey; a++)
            {
                for(int b = 0; b < window_sizex; b++)
                {
                    if(!(a==offsety && b==offsetx))//中心像素不做判断
                    {
                        census = census << 1;//左移1位
                    }
                    uchar temp_value = window.at<uchar>(a, b);
                    if(temp_value <= current_pixel ) //当前像素小于中心像素 01
                    {
                        census += 1;
                    }
                }
            }
            dst.at<double>(i+offsety, j+offsetx) = census;
        }
    }
    return dst;
}

2.2 基于特征的匹配算法（点，线，面| 区域，全局）

区域特征主要用点特征和边缘特征表示：
点特征：harris，Mooravec，KLT,SITF,SURF,BRIEF,SUSAN,FAST,CENSUS,FREAK,BRISK,ORB,光流法，A-KAZE
边缘特征：LoG，Robert，Sobel，Prewitt，Canny

2.2.1 harris角点&&& [[KLT]]

轮廓之间的交点
视角发生变化时，具有稳定的特性
该点附近的像素值在梯度方向和幅值上都有较大的变化
$E(u,v)=\sum_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)]^2$
w(x,y)是窗口函数，最简单情形就是窗口内的所有像素所对应的w权重系数均为1。但有时候，我们会将w(x,y)函数设定为以窗口中心为原点的二元正态分布。如果窗口中心点是角点时，移动前与移动后，该点的灰度变化应该最为剧烈，所以该点权重系数可以设定大些，表示窗口移动时，该点在灰度变化贡献较大；而离窗口中心(角点)较远的点，这些点的灰度变化几近平缓，这些点的权重系数，可以设定小点，以示该点对灰度变化贡献较小，那么我们自然想到使用二元高斯函数来表示窗口函数。
$E(u,v)=\sum_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)]^2 \\ Taylor 展开:f(x+u,y+v)≈f(x,y)+uf_x(x,y)+vf_y(x,y)\\ \sum_{x,y}[I(x+u,y+v)-I(x,y)]^2\\ ≈ \sum_{x,y}[I(x,y)+uI_x+vI_y-I(x,y)]^2 \\ = \sum u^2I_x^2+2uvI_xI_y+v^2I_y^2 \\ =\sum [u \quad v][\begin{matrix}I_x^2&I_xI_y\\I_xI_y& I_y^2\end{matrix}][\begin{matrix}u \\ v\end{matrix}] \\ M=\sum_{x,y}w(x,y)[\begin{matrix}I_x^2&I_xI_y\\I_xI_y& I_y^2\end{matrix}] \\ E(u,v)=\sum [u \quad v]M[\begin{matrix}u \\ v\end{matrix}] \\$
通过引入协方差矩阵，计算图片的主成分分析图片中是否存在角点。
特征值都比较大时，即窗口中含有角点
特征值一个较大，一个较小，窗口中含有边缘
特征值都比较小，窗口处在平坦区域
最后通过度量角点相应

$R=det(M)-k(trace(M))^2 \\ det(M)=\lambda_1\lambda_2 \\ trace(M)=\lambda_1+\lambda_2$

度量公式

void cornerHarris_demo( int, void* )
{
    Mat dst, dst_norm, dst_norm_scaled;
    dst = Mat::zeros( src.size(), CV_32FC1 );

    /// Detector parameters
    int blockSize = 2;
    int apertureSize = 3;
    double k = 0.04;

    /// Detecting corners
    cornerHarris( src_gray, dst, blockSize, apertureSize, k, BORDER_DEFAULT );

    /// Normalizing
    normalize( dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat() );
    convertScaleAbs( dst_norm, dst_norm_scaled );

    /// Drawing a circle around corners
    for( int j = 0; j < dst_norm.rows ; j++ )
    { for( int i = 0; i < dst_norm.cols; i++ )
        {
            if( (int) dst_norm.at<float>(j,i) > thresh )
            {
                circle( dst_norm_scaled, Point( i, j ), 5,  Scalar(0), 2, 8, 0 );
            }
        }
    }
    /// Showing the result
    namedWindow( corners_window, CV_WINDOW_AUTOSIZE );
    imshow( corners_window, dst_norm_scaled );
}

2.2.2 SITF

具有尺度不变性，旋转角度，亮度，拍摄角度对sift特征影响较小。

尺度空间极值检测
使用高斯差分函数来计算并搜索所有尺度上的图像位置，用于识别对尺度和方向不变的潜在兴趣点。
关键点定位
通过一个拟合精细的模型在每个候选位置上确定位置和尺度，关键点的选择依赖于它们的稳定程度。
方向匹配
基于局部图像的梯度方向，为每个关键点位置分配一个或多个方向，后续所有对图像数据的操作都是相对于关键点的方向、尺度和位置进行变换，从而而这些变换提供了不变形。
关键点描述
这个和HOG算法有点类似之处，在每个关键点周围的区域内以选定的比例计算局部图像梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。

import cv2
import numpy as np
#import pdb
#pdb.set_trace()#turn on the pdb prompt

#read image
img = cv2.imread('1.png',cv2.IMREAD_COLOR)
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
cv2.imshow('origin',img);

#SIFT
detector = cv2.xfeatures2d.SIFT_create()
keypoints = detector.detect(gray,None)
img = cv2.drawKeypoints(gray,keypoints,cv2.CV_8UC1)
#img = cv2.drawKeypoints(gray,keypoints,flags = cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow('test',img);
cv2.waitKey(0)
cv2.destroyAllWindows()

2.2.3 SURF

surf 和 sift 区别
 surf 和 sift 分析

void SITF(cv::Mat src1,cv::Mat src2)
{
    //sift feature detect
    int numfEATURES=400;
    //cv::Ptr<cv::xfeatures2d::SIFT> detector = cv::xfeatures2d::SIFT::create(numfEATURES);
    cv::Ptr<cv::xfeatures2d::SURF> detector = cv::xfeatures2d::SURF::create(numfEATURES);

    std::vector<cv::KeyPoint> kp1, kp2;
    detector->detect( src1, kp1 );
    detector->detect( src2, kp2 );

    ///【4】 使用SIFT算子提取特征（计算特征向量）
    Ptr<xfeatures2d::SurfDescriptorExtractor> extractor = xfeatures2d::SurfDescriptorExtractor::create();
    Mat descriptors1, descriptors2;
    extractor->compute( src1, kp1, descriptors1 );
    extractor->compute( src2, kp2, descriptors2 );
    std::cout<<kp1.size()<<"  "<<kp2.size()<<std::endl;
   // drawKeypoints(src1,kp1,src1,Scalar::all(-1),cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS);//在内存中画出特征点
    //drawKeypoints(src2,kp2,src2,Scalar::all(-1),cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS);

    //【5】使用BruteForce进行匹配
    // 实例化一个匹配器
    Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce");
    std::vector< DMatch > matches;
    //匹配两幅图中的描述子（descriptors）
    matcher->match( descriptors1, descriptors2, matches );
    double max_dist=0; double min_dist=100;
    for(int i=0;i<descriptors1.rows;i++)
    {
        double dist=matches[i].distance;
        min_dist=dist>min_dist?min_dist:dist;
        max_dist=dist<max_dist?max_dist:dist;
    }
    std::cout<<"mindist is "<<min_dist<<"max dist is "<<max_dist<<std::endl;
    std::vector< DMatch > good_matches;
    for(int i=0;i<descriptors1.rows;i++)
    {
        if(matches[i].distance<2*min_dist)
        {
            good_matches.push_back(matches[i]);
        }
    }
//    FlannBasedMatcher matcher;


    Mat img_match;
    drawMatches(src1,kp1,src2,kp2,good_matches,img_match, Scalar::all(-1), Scalar::all(-1),
                vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
    cout<<"number of matched points: "<<matches.size()<<endl;
    for( int i = 0; i < good_matches.size(); i++ )
    {
        printf( "matches [%d] quertidx: %d  -- trainidx: %d  \n",
                i, good_matches[i].queryIdx, good_matches[i].trainIdx );
    }

    imshow("matches",img_match);

    cvWaitKey();
    cvDestroyAllWindows();


}

2.2.4 SUSAN

CENSURE
当一个圆形模板在图像上滑动时，可以利用模板所覆盖区域内像素与中心点处像素灰度值的差异寻找角点。将模板内与中心点像素差异小于某个阈值的像素点组成的区域定义为USAN(univalue Segment assimilating nucleus)区域。可以发现，当模板中心处于平坦区域时，USAN面积最大，当模板中心处于边界上时，则USAN面积约为最大值的1/2，当模板中心位于角点处上时，USAN面积约为最大值的1/4。也就是说USAN面积越接近1/4，其为角点的概率就越大，因此将这种算法称为SUSAN（Smallest Univalue Segment Assimilating Nucleus）。
SUSAN角点检测步骤：

使用一个圆形模板在图像上滑动，目标像素与模板中心重叠，比较模板内各像素灰度值与目标像素灰度值的差异，构成USAN区域。

$c(x,y)=\{\begin{matrix}1, \quad |I(x,y)-I(x_0,y_0)|<t \\ 0 \quad |I(x,y)-I(x_0,y_0)|≥t\end{matrix}$

计算SUSAN区域面积

$n(x_0,y_0)=\sum_{x,y}c(x,y)$

3. 计算角点响应值

$R(x_0,y_0)=\{\begin{matrix}g-n(x_0,y_0), \quad n(x_0,y_0)<g \\ 0, \quad n(x_0,y_0)≥g\end{matrix}$
阈值可取g=nmax/2，USAN面积达到最小时，角点响应值达到最大
4. 对角点响应值做非极大值抑制
计算处角点响应值中非0值对应的位置即为候选角点，在以每个候选角点为中心的一定大小的窗口中，保留局部极大值点作为最终角点。

void SUSAN(cv::Mat src)
{
    uchar *data0, *data1, *data2;
    int same, max, min, thresh;
    Mat img=Mat(src.size(),CV_8SC1); //中间图像
    Mat dst=Mat(src.size(),CV_8SC1); //结果图像
    int height = src.rows;
    int width = src.cols;
    int step = 3;
    int channels = src.channels();
    data0 = (uchar*)src.data;
    data1 = (uchar*)img.data;
    data2 = (uchar*)dst.data;
    memset(data1,0,width*height*sizeof (char));
    memset(data2,0,width*height*sizeof (char));
    int  g = 20; //核值相似区域中点个数的阈值

    //模版 x 和 y的坐标的偏移量
    int OffSetX[37] =
    {
        -1, 0, 1,
        -2, -1, 0, 1, 2,
        -3, -2, -1, 0, 1, 2, 3,
        -3, -2, -1, 0, 1, 2, 3,
        -3, -2, -1, 0, 1, 2, 3,
        -2, -1, 0, 1, 2,
        -1, 0, 1
    };
    int OffSetY[37] =
    {
        -3, -3, -3,
        -2, -2, -2, -2, -2,
        -1, -1, -1, -1, -1, -1, -1,
        0, 0, 0, 0, 0, 0, 0,
        1, 1, 1, 1, 1, 1, 1,
        2, 2, 2, 2, 2,
        3, 3, 3
    };

    //求阈值
    max = min = data0[0];
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j<width; j++)
        {
            if (data0[i*step + j]>max) max = data0[i*step + j];
            if (data0[i*step + j]<min) min = data0[i*step + j];
        }
    }
//    thresh = (max - min) / 10;//可选取其他方法
    thresh=10;
    for (int i = 3; i < height - 3; i++)
    {
        for (int j = 3; j < width - 3; j++)
        {
            same = 0;
            for (int k = 0; k < 37; k++)
            {
                if (abs((int)data0[(i + OffSetY[k])*step + (j + OffSetX[k])] - (int)data0[i*step + j]) < thresh)
                    same++;
            }
            if (same < g)//g值可改
                data1[i*step + j] = g-same;
            else
                data1[i*step + j] = 0;
        }
    }


    //非极大值抑制
    int i_s[8] = { -1, -1, -1, 0, 0, 1, 1, 1 };
    int j_s[8] = { -1, 0, 1, -1, 1, -1, 0, 1 };
    int flag;
    for (int i = 4; i < height - 4; i++)
    {
        for (int j = 4; j < width -4; j++)
        {
            flag = 0;
            for (int k = 0; k < 8; k++)
            {
                if (data1[i*step + j] <= data1[ (i + i_s[k]) * step + (j + j_s[k]) ] )
                {
                    flag = 1;
                    break;
                }
            }
            if (flag == 0)
            {
                data2[i*step + j] = 255;
            }
            else
            {
                data2[i*step + j] = 0;
            }
        }
    }
    cv::imshow("des",dst);
    cv::waitKey();
}

2.2.5 FAST

fast角点检测1
fast角点检测2
若某个像素点与周围足够多的像素点处于不同的区域，则该像素点可能为角点。

void FAST(int thred, cv::Mat src)
{
    std::vector<KeyPoint> keypoints;
    Mat dst = src.clone();
    Ptr<FastFeatureDetector> detector = FastFeatureDetector::create(thred);
    detector->detect(src,keypoints);
    drawKeypoints(dst, keypoints, dst, Scalar::all(-1), DrawMatchesFlags::DRAW_OVER_OUTIMG);
    imshow("output", dst);

}

2.2.6 MORAVEC

Moravec 角点检测算法是最早的角点检测算法之一。该算法将角点定义为具有低“自相关性”的点。算法会检测图像的每一个像素，将像素周边的一个邻域作为一个patch，并检测这个patch和周围其他patch的相关性。这种相关性通过两个patch间的平方差之和(SSD)来衡量，SSD 值越小则相似性越高。
如果像素位于平滑图像区域内，周围的patch都会非常相似。如果像素在边缘上，则周围的patch在与边缘正交的方向上会有很大差异，在与边缘平行的方向上则较为相似。而如果像素是各个方向上都有变化的特征点，则周围所有的patch都不会很相似。
Moravec 会计算每个像素patch和周围patch的SSD最小值作为强度值，取局部强度最大的点作为特征点。

Moravec 角点检测

Mat MoravecCorners(Mat srcImage, int kSize, int threshold)
{
    //角点检测的结果图
    Mat resMorMat = srcImage.clone();
    int r = kSize / 2;
    //获取图像的高和宽
    const int nRows = srcImage.rows;
    const int nCols = srcImage.cols;
    int nCount = 0;
    //保存角点的坐标
    CvPoint *pPoint = new CvPoint[nRows*nCols];
    //遍历图像
    for (int i = r; i < srcImage.rows - r; i++)
    {
        for (int j = r; j < srcImage.cols - r; j++)
        {
            int wV1, wV2, wV3, wV4;
            wV1 = wV2 = wV3 = wV4 = 0;
            //计算水平方向窗内的兴趣值
            for (int k = -r; k <= r; k++)
            {
                for (int m = -r; m <= r; m++)
                {
                    //判断移动的过程中是否越界，越界的话就跳过当前的循环，以免出错
                    int a = i + k ;
                    int b = j + m + 1;
                    if (b>=srcImage.cols)
                    {
                        continue;
                    }
                    wV1 += (srcImage.at<uchar>(i + k, j + m + 1) - srcImage.at<uchar>(i + k, j + m))
                            *(srcImage.at<uchar>(i + k, j + m + 1) - srcImage.at<uchar>(i + k, j + m));
                }
            }

            //计算垂直方向窗内的兴趣值
            for (int k = -r; k <= r; k++)
            {
                for (int m = -r; m <= r; m++)
                {
                    int a = i + k + 1;
                    int b = j + m ;
                    if (a >=srcImage.rows)
                    {
                        continue;
                    }
                    wV2 += (srcImage.at<uchar>(i + k + 1, j + m) - srcImage.at<uchar>(i + k, j + m))
                            *(srcImage.at<uchar>(i + k + 1, j + m) - srcImage.at<uchar>(i + k, j + m));
                }
            }

            //计算45°方向窗内的兴趣值
            for (int k = -r; k <= r; k++)
            {
                for (int m = -r; m <= r; m++)
                {
                    int a = i + k + 1;
                    int b = j + m + 1;
                    if (a >=srcImage.rows || b >=srcImage.cols)
                    {
                        continue;
                    }
                    wV3 += (srcImage.at<uchar>(i + k + 1, j + m + 1) - srcImage.at<uchar>(i + k, j + m))
                            *(srcImage.at<uchar>(i + k + 1, j + m + 1) - srcImage.at<uchar>(i + k, j + m));
                }
            }

            //计算135°方向窗内的兴趣值
            for (int k = -r; k <= r; k++)
            {
                for (int m = -r; m <= r; m++)
                {
                    int a = i + k + 1;
                    int b = j + m - 1;
                    if (a >=srcImage.rows || b < 0)
                    {
                        continue;
                    }

                    wV4 += (srcImage.at<uchar>(a, b) - srcImage.at<uchar>(i + k, j + m))
                            *(srcImage.at<uchar>(a, b) - srcImage.at<uchar>(i + k, j + m));
                }
            }
            int  value = min(min(wV1, wV2), min(wV3, wV4));
            //如果兴趣值大于阈值，那么将坐标存入数组中
            if (value > threshold)
            {
                pPoint[nCount] = cvPoint(j, i);
                nCount++;
            }
        }
    }
    //绘制兴趣点
    for (int i = 0; i < nCount; i++)
    {
        circle(resMorMat, pPoint[i], 5, Scalar(255, 0, 0));
    }
    return resMorMat;
}

2.2.6 BRISK

brisk算法描述
BRISK算法主要利用FAST9-16进行特征点检测。BRISK算法中构造了图像金字塔进行多尺度表达。
构建尺度空间
特征点检测
非极大值抑制
亚像素插值
特征点描述

        Ptr<Feature2D> detector = BRISK::create();
        vector<KeyPoint> keypoints_obj;
        vector<KeyPoint> keypoints_scene;
        Mat descriptor_obj, descriptor_scene;
        detector->detectAndCompute(img1, Mat(), keypoints_obj, descriptor_obj);
        detector->detectAndCompute(img2, Mat(), keypoints_scene, descriptor_scene);
        // matching
        BFMatcher matcher(NORM_L2);
        vector<DMatch> matches;
        matcher.match(descriptor_obj, descriptor_scene, matches);
        double t2 = (double)getTickCount();
        double t = (t2 - t1) / getTickFrequency();
        cout << "spend time : " << t << "s" << endl;

2.2.7 ORB(Oriented FAST and Rotated BRIEF)

ORB
ORB

特征点检测
ORB采用FAST算法来检测特征点。FAST核心思想是找出与周围点不同的点。
特征点描述
ORB采用优化的BRIEF算法计算一个特征点的描述子。BRIEF算法的核心思想是在关键点P的周围以一定模式选取N个点对，把这N个点对的比较结果组合起来作为描述子。ORB中以关键点为圆心，计算一定范围内质心为特征点，这样得到的特征点具有旋转不变性。
特征点匹配

void TORB(cv::Mat scene, cv::Mat box,float ratio)
{
    std::vector <cv::KeyPoint> keypoints_obj,keypoints_sence;
    Mat descriptors_box, descriptors_sence;
    Ptr<ORB> detector = ORB::create();
    detector->detectAndCompute(scene, Mat(), keypoints_sence, descriptors_sence);
    detector->detectAndCompute(box, Mat(), keypoints_obj, descriptors_box);
    vector<DMatch> matches;
    // 初始化flann匹配
    // Ptr<FlannBasedMatcher> matcher = FlannBasedMatcher::create(); // default is bad, using local sensitive hash(LSH)
    Ptr<DescriptorMatcher> matcher = makePtr<FlannBasedMatcher>(makePtr<flann::LshIndexParams>(12, 20, 2));
    matcher->match(descriptors_box, descriptors_sence, matches);
    // 发现匹配
    vector<DMatch> goodMatches;
    printf("total match points : %d\n", matches.size());
    float maxdist = 0;
    for (unsigned int i = 0; i < matches.size(); ++i)
    {
        printf("dist : %.2f \n", matches[i].distance);
        maxdist = max(maxdist, matches[i].distance);
    }
    for (unsigned int i = 0; i < matches.size(); ++i)
    {

        if(matches[i].distance < maxdist*ratio)
            goodMatches.push_back(matches[i]);
    }

    Mat dst;
    drawMatches(box, keypoints_obj, scene, keypoints_sence, goodMatches, dst);
    imshow("output", dst);
    waitKey(10);
}

2.2.8 A-KAZE

A-KAZE
由于SURF，SIFT算法采用高斯尺度空间，高斯模糊不保留对象边界信息并且在所有尺度上平滑到相同程度的细节与噪声，影响定位的准确性和独特性。
因此，KAZE采用非线性滤波构建尺度空间：双边滤波、非线性扩散滤波。
KAZE

非线性尺度空间构建
KAZE算法中，作者通过非线性扩散滤波和加性算子分裂（AOS）算法来构造非线性尺度空间。
1.1 非线性扩散滤波

非线性扩散滤波方法是将图像亮度（L）在不同尺度上的变化视为某种形式的流动函数的散度，可以通过非线性偏微分方程来描述：
$\frac{\partial L}{\partial t}=div(c(x,y,t)\cdot \Delta L) \\ c(x,y,t)=g(|\Delta L_\sigma(x,y,t)|) \quad 传导函数 \\ g_1=exp(-\frac{|\Delta L_\sigma|^2}{k^2}) \quad 能够保留高对比度的边缘\\ g_2=\frac{1}{1+\frac{|\Delta L_\sigma|^2}{k^2}} \quad 能够保留宽度较大的区域\\ g_3=\{\begin{matrix} 1 & ,|\Delta L_\sigma|^2=0 \\ 1-exp(-\frac{3.315}{(|\Delta L_\sigma|/k)^8}) &,|\Delta L_\sigma|^2>0\end{matrix} \quad 能够有效平滑区域内部而保留边界信息$

1.2 AOS算法

非线性扩散滤波中的偏微分方程没有解析界。因此，需要使用数值方法来逼近微分方程。将上述偏微分方程离散化：
$\frac {L^{i+1}-L^i}{\tau}=\sum^m_{l=1}A_l(L^i)L^{i+1}) \\ A_L表示图像在各个维度（l）上传导性的矩阵，\tau 位时间步长 \\ L^{i+1}=(I-\tau\sum^m_{l=1}A_l(L^i)^{-1}L^i$

特征点的检测与精确定位
通过寻找不同尺度归一化后的Hessian行列式的局部极大值（或者极小值）点来实现
特征点主方向的确定
这一过程和SURF一样。若特征点的尺度参数为σi，则搜索半径设为6σi。在这个圆形领域内做一个60度的扇形区域，统计这个扇形区域内的haar小波特征总和，然后转动扇形区域，再统计小波特征总和。小波特征总和最大的方向为主方向
特征描述子的生成
对于尺度参数为σi的特征点，在梯度图像上以特征点为中心取一个24σi×24σi的窗口，并将窗口划分为4×4个子区域，每个子区域大小为9σi×9σi，相邻的子区域有宽度为2σi的交叠带（此处我认为应该是相邻的子区域有宽度为4σi的交叠带，不然24σi不够划分为4×4个子区域）。每个子区域都用一个高斯核（σ1 =2.5σi）进行加权，然后计算出长度为4的子区域描述向量：
再通过另一个大小为4×4的高斯窗口（σ2 =1.5σi）对每个子区域的向量dv进行加权，最后进行归一化处理。这样就得到了4×4×4=64维的描述向量。

void AKZE(cv::Mat img1,cv::Mat img2)
{
    Ptr<KAZE> detector = KAZE::create();
    vector<KeyPoint> keypoints1,keypoints2;
    Mat descriptors_box, descriptors_sence;
//    double t1 = getTickCount();//获取当前的内部滴答数（计算机从一开始就会一直计数）
    detector->detect(img1, keypoints1, descriptors_box);
//    double t2 = getTickCount();//获取当前的内部滴答数
//    //算法执行的时间
//    double tkaze = 1000 * (t2 - t1) / getTickFrequency();//getTickFrequency()每一秒的滴答数
//    printf("KAZE Time consume(ms):%f", tkaze);

    detector->detect(img2,keypoints2,descriptors_sence);
    vector<DMatch> matches;
    // 初始化flann匹配
    // Ptr<FlannBasedMatcher> matcher = FlannBasedMatcher::create(); // default is bad, using local sensitive hash(LSH)
    Ptr<DescriptorMatcher> matcher = makePtr<FlannBasedMatcher>(makePtr<flann::LshIndexParams>(12, 20, 2));
    matcher->match(descriptors_box, descriptors_sence, matches);
    // 发现匹配
    vector<DMatch> goodMatches;
    printf("total match points : %d\n", matches.size());
    float maxdist = 0;
    for (unsigned int i = 0; i < matches.size(); ++i)
    {
        printf("dist : %.2f \n", matches[i].distance);
        maxdist = max(maxdist, matches[i].distance);
    }
    for (unsigned int i = 0; i < matches.size(); ++i)
    {

        if(matches[i].distance < maxdist*0.4)
            goodMatches.push_back(matches[i]);
    }

    Mat dst;
    drawMatches(img1, keypoints1, img2, keypoints2, goodMatches, dst);
    imshow("output", dst);
    waitKey(10);
}

2.2.9 LoG（Laplacian of Gaussian）

Dog：一阶边缘提取
LoG：二阶边缘提取，采用laplacian-Gaussian边缘提取算法，先高斯滤波后拉普拉斯边缘提取。
$\Delta G(x,y)^2=\frac{\partial G(x,y)^2}{\partial x^2}+\frac{\partial G(x,y)^2}{\partial y^2}=[\frac{x^2}{\delta ^2}-\frac{1}{\delta ^2}]e^{-\frac{x^2+y^2}{2\delta ^2}} \\ =[\frac{x^2+y^2-2\delta ^2}{\delta^4}]e^{-\frac{x^2+y^2}{2\delta ^2}}$

GaussianBlur(img, imgGaussian, Size(3, 3), 1);//高斯模糊
Laplacian(imgGaussian, img16S, 3);
convertScaleAbs(img16S, imgLoG, 1);

2.2.10 Robert

Mat kernel_x = (Mat_<int>(2, 2) << 1, 0, 0, -1);
    filter2D(img, img, -1, kernel_x, Point(-1, -1), 0.0);  //-1是图像的深度，-1程序自动判断

Mat kernel_Y = (Mat_<int>(2, 2) << 0,1,-1,0);
    filter2D(img, img, -1, kernel_Y, Point(-1, -1), 0.0);

2.2.11 Sobel

//Sobel X方向
Mat Sobel_X = (Mat_<int>(3, 3) << -1,0,1,-2,0,2,-1,0,1);
filter2D(src, sobel_x, -1, Sobel_X, Point(-1, -1), 0.0);
imshow("Sobel_x", sobel_x);

//Sobel y方向
Mat Sobel_Y = (Mat_<int>(3, 3) << -1,-2,-1,0,0,0,1,2,1);
filter2D(src, dst_y, -1, Sobel_Y, Point(-1, -1), 0.0);
imshow("Sobel_y", dst_y);

Sobel(img, img16S, 3, 1, 0);
    convertScaleAbs(img16S, imgSobelx, 1);
    Sobel(img, img16S, 3, 0, 1);
    convertScaleAbs(img16S, imgSobely, 1);
    add(imgSobelx, imgSobely, imgSobel);
    namedWindow("Sobel");
    imshow("Sobel", imgSobel);
    waitKey(0);

2.2.12 Prewitt

边缘检测算子

2.2.13 Canny

 Canny(img, imgCanny, 100, 200);

2.3 BRIEF

由于SIFT和SURF具有庞大的特征向量结构，而且在特征匹配中通常不需要这么多，所以可以通过
（1）PCA,LDA等特征降维的方法来压缩特征描述子的维度（2）LSH,SIFT的特征描述子转化为一个二值码串，然后计算码串的汉明距离进行特征点之间的匹配。
在现代计算机中计算汉明距离可以用异或操作，然后计算二进制位数来实现，BRIEF就提供了一种计算二值串的捷径，具体计算步骤如下：
（1）平滑图像
（2）在特征点周围选择一个Patch，在这个Patch内通过一种选定的方法来挑选出来nd个点对。然后对于每一个点对（p,q）比较两个点的亮度值，如果I§>I(q),则这个点对生成的二值串值为1，如果I§<I(q),则这个点对生成的二值串值为-1，否则为0.所以nd个点对，会生成nd长度的二进制串。
（3）nd=256[default] | 128 | 512，在这种情况下，非匹配的汉明距离呈现均值为nd特征的高斯分布。一旦维数选定了，我们就可以用汉明距离来匹配这些描述子了。
（4）BRIEF是一种特征描述符，它不提供特提取特征点的方法。所以必须使用一种特征点定位的方法，如FAST,SURF,SITF，CENSURE等。
点对的选择方法

在图像快内平均采样
p,q都符合（0，1/25S*S）的高斯分布
p符合（0，1/25S*S），q符合（0，1/100S*S）的高斯分布
在空间量化极坐标下的离散位置随机采样
p固定在(0,0),q在周围随机采样

void BRIEF(cv::Mat img_1,cv::Mat img_2)
{
    // -- Step 1: Detect the keypoints using STAR Detector
       std::vector<KeyPoint> keypoints_1,keypoints_2;
       cv::Ptr<cv::xfeatures2d::StarDetector> detector = cv::xfeatures2d::StarDetector::create();

       detector->detect(img_1, keypoints_1);
       detector->detect(img_2, keypoints_2);
       std::cout<<keypoints_1.size()<<"  "<<keypoints_2.size()<<std::endl;
       // -- Stpe 2: Calculate descriptors (feature vectors)
       Ptr<xfeatures2d::BriefDescriptorExtractor> brief = xfeatures2d::BriefDescriptorExtractor::create();
       Mat descriptors_1, descriptors_2;
       brief->compute(img_1, keypoints_1, descriptors_1);
       brief->compute(img_2, keypoints_2, descriptors_2);
       //-- Step 3: Matching descriptor vectors with a brute force matcher
       BFMatcher matcher(NORM_HAMMING);
       std::vector<DMatch> mathces;
       matcher.match(descriptors_1, descriptors_2, mathces);
       std::nth_element(mathces.begin(), // initial position
                        mathces.begin()+29, // position of the sorted element
                        mathces.end()); // end position
       // remove all elements after the 30th
       mathces.erase(mathces.begin()+30, mathces.end());

       // -- dwaw matches
       Mat img_mathes;
       drawMatches(img_1, keypoints_1, img_2, keypoints_2, mathces, img_mathes);
       // -- show
       imshow("Mathces", img_mathes);
       waitKey(0);
}

2.4 立体匹配

左右相机的匹配，用灰度模板匹配，对所有像素点匹配
1.区域立体匹配算法
给定一幅图像上的一点，选取该点邻域内的一个子窗口，在另一幅图像中的一个区域内，根据某种相似性，寻找与子窗口图像最相似的窗口，得到的匹配窗口中对应的像素点就为该像素的匹配点。可获得稠密视差图。
2.基于特征的立方匹配算法
基于几何特征信息（边缘、线、轮廓、兴趣点、角点和几何基元等）提取图像的几何特征点，针对几何特征点进行视差估计，利用得到的视差信息重建三维空间场景。可获得稀疏视差图，通过插值可获得稠密视差图。
3.基于相位立体匹配算法
假设在图像对应点中，其频率范围内局部相位相等，在频率范围内进行视差估计。

算法分类	算法分类	算法描述	优点	缺点
局部法	区域匹配	在两幅图像给定区域内搜索中心点最大相关值或最小偏差值	直观，简单	对噪声敏感，弱文理区域易造成匹配
	特征匹配	寻找特征（点，线，边缘等）之间的相似度，而不是亮度本身	速度快，精度高，克服区域法对弱文理区的敏感性	不适用于缺少明显特征的环境
	相位匹配	对带通滤波信号的相位信息进行处理，寻找局部相位相等的对应点	可直接得到亚像素级竞速，抑制高频噪声	奇异点不稳定
全局法	动态规划	在左右图像对应的扫描线上寻找最小匹配代价路径	提高局部无文理区域匹配精度	视差图上易出现条纹瑕疵
	图割法	根据能量函数构造合适的图，求其最小割（最大流）	能有效去除条纹瑕疵	较动态规划算法的计算复杂度高
	置信度传播	用马尔科夫随机场（MRF）进行建模，求MRF最优解来获得最佳试图差	精度高	计算复杂度高

立体匹配（Stereo matching）的步骤如下：
1：预处理：亮度归一化，去噪，图像增强，滤波等等
2：匹配Cost计算Cost aggregation
每个像素点的matching cost可用下图所示的两种方式表示
将所有像素的Cost加起来，选择和最小的方案。
3：全局最优化
4：后处理：Refinement 精化
一般基于像素点的匹配含有较大的噪音，所以人们更愿意使用基于局部或全局的匹配算法。常用的立体匹配算法可分为基于局部的匹配方法和基于全局的匹配方法。
（1）基于局部的匹配算法：块匹配，模板匹配，SIFT[1]，SURF[2]，AKAZE[3]
（2）基于全局的匹配算法：首先定义Cost Function，然后根据约束条件优化Cost Function.

Graph Cut：将每一个像素看作Graph中的一个Node，然后假设该图有L个不同的深度值，则添加L个新的Node作为Label，然后求视差图的问题就变成了多重Labelling的问题了，具体Graph Cut参照文献。

Yuri Boykov and Vladimir Kolmogorov (2003), “Computing Geodesics and
Minimal Surfaces via Graph Cuts”

Belief Propagation: 将图像看作马尔科夫随机场，则求视差问题变成了最大化联合后验概率的问题，每一个像素对其深度值的猜测来自于其邻居给的信息，然后再把信息传给其他的邻居。每一个信息有一个可靠度（概率）。最后经过多次循环深度值会收敛。

2.4.1 BM算法（Block Matching）

四种block match算法

 cv::Mat imgL=cv::imread("1.png");
    cv::Mat imgR=cv::imread("2.png");
    cvtColor(imgL,imgL,CV_BGR2GRAY);
    cvtColor(imgR,imgR,CV_BGR2GRAY);

    int numberOfDisparities = 48;
    cv::Ptr<cv::StereoBM> bm = cv::StereoBM::create(16, 9);
    cv::Rect roi1, roi2;
    bm->setROI1(roi1);
    bm->setROI2(roi2);
    bm->setPreFilterCap(31);
    bm->setBlockSize(9);
    bm->setMinDisparity(0);
    bm->setNumDisparities(numberOfDisparities);
    bm->setTextureThreshold(10);
    bm->setUniquenessRatio(15);
    bm->setSpeckleWindowSize(100);
    bm->setSpeckleRange(32);
    bm->setDisp12MaxDiff(1);
    Mat disp;
    bm->compute(imgL, imgR, disp);
    disp.convertTo(disp, CV_8U, 255 / (numberOfDisparities*16.));   //将16位符号整形的视差矩阵转换为8位无符号整形矩阵
    cv::imshow("disparity", disp);
    cv::waitKey();

    return 0;

2.4.2 SGBM（Semi-Global Block matching）

SGBM

 cv::Mat imgL=cv::imread("1.png");
    cv::Mat imgR=cv::imread("2.png");
    cvtColor(imgL,imgL,CV_BGR2GRAY);
    cvtColor(imgR,imgR,CV_BGR2GRAY);


    int numberOfDisparities = 48;
    int SADWindowSize=11;
    int uniquenessRatio = 15;
    int speckleWindowSize = 50;
    int speckleRange = 32;
    cv::Ptr<cv::StereoSGBM> sgbm = cv::StereoSGBM::create(0, numberOfDisparities, SADWindowSize);
    int cn = imgL.channels();
    sgbm->setP1(8 * cn * SADWindowSize * SADWindowSize);
    sgbm->setP2(32 * cn * SADWindowSize * SADWindowSize);
    sgbm->setPreFilterCap(63);
    sgbm->setUniquenessRatio(uniquenessRatio);
    sgbm->setSpeckleWindowSize(speckleWindowSize);
    sgbm->setSpeckleRange(speckleRange);
    sgbm->setDisp12MaxDiff(1);

    //计算并显示视差
    Mat disp;
    sgbm->compute(imgL, imgR, disp);
    disp.convertTo(disp, CV_8U, 255 / (numberOfDisparities*16.));   //将16位符号整形的视差矩阵转换为8位无符号整形矩阵
    cv::imshow("disparity", disp);
    cv::waitKey();

2.4.3 GC(Graph Cut)

int main()
{
 
	//IplImage * img1 = cvLoadImage("left.png",0);
	//IplImage * img2 = cvLoadImage("right.png",0);
	//IplImage * img1 = cvLoadImage("tsukuba_l.png",0);
	//IplImage * img2 = cvLoadImage("tsukuba_r.png",0);
	IplImage * img1 = cvLoadImage("left.png",0);
	IplImage * img2 = cvLoadImage("right.png",0);
	CvStereoGCState* GCState=cvCreateStereoGCState(64,3);
	assert(GCState);
	cout<<"start matching using GC"<<endl;
	CvMat* gcdispleft=cvCreateMat(img1->height,img1->width,CV_16S);
	CvMat* gcdispright=cvCreateMat(img2->height,img2->width,CV_16S);
	CvMat* gcvdisp=cvCreateMat(img1->height,img1->width,CV_8U);
	int64 t=getTickCount();
	cvFindStereoCorrespondenceGC(img1,img2,gcdispleft,gcdispright,GCState);
	t=getTickCount()-t;
	cout<<"Time elapsed:"<<t*1000/getTickFrequency()<<endl;
	//cvNormalize(gcdispleft,gcvdisp,0,255,CV_MINMAX);
	//cvSaveImage("GC_left_disparity.png",gcvdisp);
	cvNormalize(gcdispright,gcvdisp,0,255,CV_MINMAX);
	cvSaveImage("GC_right_disparity.png",gcvdisp);
 
 
	cvNamedWindow("GC_disparity",0);
	cvShowImage("GC_disparity",gcvdisp);
	cvWaitKey(0);
	cvReleaseMat(&gcdispleft);
	cvReleaseMat(&gcdispright);
	cvReleaseMat(&gcvdisp);
	return 0;
}

立体视觉匹配
GC-SGBM-DP-BM

2.4.4 DP(Dynamic Programming)

2.4.5 preFilterCap（）匹配图像预处理

static void prefilterXSobel(const cv::Mat& src, cv::Mat& dst, int ftzero)
{
    int x, y;
    const int OFS = 256 * 4, TABSZ = OFS * 2 + 256;
    uchar tab[TABSZ];
    cv::Size size = src.size();

    for (x = 0; x < TABSZ; x++)
        tab[x] = (uchar)(x - OFS < -ftzero ? 0 : x - OFS > ftzero ? ftzero * 2 : x - OFS + ftzero);
    uchar val0 = tab[0 + OFS];

    for (y = 0; y < size.height - 1; y += 2)
    {
        const uchar* srow1 = src.ptr<uchar>(y);
        const uchar* srow0 = y > 0 ? srow1 - src.step : size.height > 1 ? srow1 + src.step : srow1;
        const uchar* srow2 = y < size.height - 1 ? srow1 + src.step : size.height > 1 ? srow1 - src.step : srow1;
        const uchar* srow3 = y < size.height - 2 ? srow1 + src.step * 2 : srow1;
        uchar* dptr0 = dst.ptr<uchar>(y);
        uchar* dptr1 = dptr0 + dst.step;

        dptr0[0] = dptr0[size.width - 1] = dptr1[0] = dptr1[size.width - 1] = val0;
        x = 1;
        for (; x < size.width - 1; x++)
        {
            int d0 = srow0[x + 1] - srow0[x - 1], d1 = srow1[x + 1] - srow1[x - 1],
                d2 = srow2[x + 1] - srow2[x - 1], d3 = srow3[x + 1] - srow3[x - 1];
            int v0 = tab[d0 + d1 * 2 + d2 + OFS];
            int v1 = tab[d1 + d2 * 2 + d3 + OFS];
            dptr0[x] = (uchar)v0;
            dptr1[x] = (uchar)v1;
        }
    }

    for (; y < size.height; y++)
    {
        uchar* dptr = dst.ptr<uchar>(y);
        x = 0;
        for (; x < size.width; x++)
            dptr[x] = val0;
    }
}

void mySobelX(cv::Mat srcImg, cv::Mat dstImg, int preFilterCap)
{
    assert(srcImg.channels() == 1);
    int radius = 1;
    int width = srcImg.cols;
    int height = srcImg.rows;
    uchar *pSrcData = srcImg.data;
    uchar *pDstData = dstImg.data;
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            int idx = i*width + j;
            if (i >= radius && i < height - radius && j >= radius && j < width - radius)
            {
                int diff0 = pSrcData[(i - 1)*width + j + 1] - pSrcData[(i - 1)*width + j - 1];
                int diff1 = pSrcData[i*width + j + 1] - pSrcData[i*width + j - 1];
                int diff2 = pSrcData[(i + 1)*width + j + 1] - pSrcData[(i + 1)*width + j - 1];

                int value = diff0 + 2 * diff1 + diff2;
                if (value < -preFilterCap)
                {
                    pDstData[idx] = 0;
                }
                else if (value >= -preFilterCap && value <= preFilterCap)
                {
                    pDstData[idx] = uchar(value + preFilterCap);
                }
                else
                {
                    pDstData[idx] = uchar(2 * preFilterCap);
                }

            }
            else
            {
                pDstData[idx] = 0;
            }
        }
    }
}

2.4.6 filterSpeckles（）视差图后处理

OpenCV中有对应的API函数，
void filterSpeckles(InputOutputArray img, double newVal, int maxSpeckleSize, double maxDiff, InputOutputArray buf=noArray() )

2.4.7 双目立体匹配获得深度图步骤

双目立体匹配步骤

相机标定
(1)内参标定，获得fx, fy, cx, cy, k1, k2, p1, p2, k3。
(2)外参标定，获得摄像机坐标系和世界坐标系之间的旋转R和平移T关系。
如果两个相机的内参均已知，并且知道各自与世界坐标系之间的R1、T1和R2，T2，就可以算出这两个相机之间的Rotation和Translation，也就找到了从一个相机坐标系到另一个相机坐标系之间的位置转换关系。
或者：
（1）OpenCV中的光流法提取匹配特征点对，pts1和pts2。
（2）利用特征点对pts1和pts2，以及内参矩阵camK，解算出本质矩阵E：

cv::Mat E = cv::findEssentialMat(tmpPts1, tmpPts2, camK, CV_RANSAC);

（3）利用本质矩阵E解算出两个摄像机之间的Rotation和Translation，也就是两个摄像机之间的外参。

    cv::Mat R1, R2;
    cv::decomposeEssentialMat(E, R1, R2, t);

    R = R1.clone();
    t = -t.clone();

双目图像校正
（1）畸变校正
（2）立体校正

得到两个摄像机之间的 Rotation和Translation之后，要用下面的API对两幅图像进行立体对极线校正，这就需要算出两个相机做对极线校正需要的R和T，用R1,T1, R2, T2表示，以及透视投影矩阵P1，P2：

cv::stereoRectify(camK, D, camK, D, imgL.size(), R, -R*t,  R1, R2, P1, P2, Q);

得到上述参数后，就可以使用下面的API进行对极线校正操作了，并将校正结果保存到本地：

　　cv::initUndistortRectifyMap(P1(cv::Rect(0, 0, 3, 3)), D, R1, P1(cv::Rect(0, 0, 3, 3)), imgL.size(), CV_32FC1, mapx, mapy);
cv::remap(imgL, recImgL, mapx, mapy, CV_INTER_LINEAR);
cv::imwrite("data/recConyL.png", recImgL);

cv::initUndistortRectifyMap(P2(cv::Rect(0, 0, 3, 3)), D, R2, P2(cv::Rect(0, 0, 3, 3)), imgL.size(), CV_32FC1, mapx, mapy);
cv::remap(imgR, recImgR, mapx, mapy, CV_INTER_LINEAR);
cv::imwrite("data/recConyR.png", recImgR);

查看对极线校正结果是否准确，可以通过观察若干对应点是否在同一行上粗略估计得出

立体匹配算法获得视差图
（1）SGBM,BM,GC,等获得视差图
（2）视差图中视差值不可靠的视差大多数是由于遮挡引起，或者光照不均匀引起。既然牛逼如SGBM也觉得不可靠，那与其留着做个空洞，倒不如用附近可靠的视差值填充一下。
（3）视差图转换为深度图：视差的单位是像素（pixel），深度的单位往往是毫米（mm）表示。而根据平行双目视觉的几何关系（此处不再画图推导，很简单），可以得到下面的视差与深度的转换公式：
$d e p t h = (f * b a s e l i n e) / d i s p$
上式中，depth表示深度图；f表示归一化的焦距，也就是内参中的fx； baseline是两个相机光心之间的距离，称作基线距离；disp是视差值。等式后面的均已知，深度值即可算出。

 /*
函数作用：视差图转深度图
输入：
　　dispMap ----视差图，8位单通道，CV_8UC1
　　K       ----内参矩阵，float类型
输出：
　　depthMap ----深度图，16位无符号单通道，CV_16UC1
*/

void disp2Depth(cv::Mat dispMap, cv::Mat &depthMap, cv::Mat K)
{
    int type = dispMap.type();

    float fx = K.at<float>(0, 0);
    float fy = K.at<float>(1, 1);
    float cx = K.at<float>(0, 2);
    float cy = K.at<float>(1, 2);
    float baseline = 65; //基线距离65mm

    if (type == CV_8U)
    {
        const float PI = 3.14159265358;
        int height = dispMap.rows;
        int width = dispMap.cols;

        uchar* dispData = (uchar*)dispMap.data;
        ushort* depthData = (ushort*)depthMap.data;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                int id = i*width + j;
                if (!dispData[id])  continue;  //防止0除
                depthData[id] = ushort( (float)fx *baseline / ((float)dispData[id]) );
            }
        }
    }
    else
    {
        cout << "please confirm dispImg's type!" << endl;
        cv::waitKey(0);
    }
}

视差图和深度图的空洞填充步骤如下：

① 以视差图dispImg为例。计算图像的积分图integral，并保存对应积分图中每个积分值处所有累加的像素点个数n（空洞处的像素点不计入n中，因为空洞处像素值为0，对积分值没有任何作用，反而会平滑图像）。

② 采用多层次均值滤波。首先以一个较大的初始窗口去做均值滤波（积分图实现均值滤波就不多做介绍了，可以参考我之前的一篇博客），将大区域的空洞赋值。然后下次滤波时，将窗口尺寸缩小为原来的一半，利用原来的积分图再次滤波，给较小的空洞赋值（覆盖原来的值）；依次类推，直至窗口大小变为3x3，此时停止滤波，得到最终结果。

③ 多层次滤波考虑的是对于初始较大的空洞区域，需要参考更多的邻域值，如果采用较小的滤波窗口，不能够完全填充，而如果全部采用较大的窗口，则图像会被严重平滑。因此根据空洞的大小，不断调整滤波窗口。先用大窗口给所有空洞赋值，然后利用逐渐变成小窗口滤波覆盖原来的值，这样既能保证空洞能被填充上，也能保证图像不会被过度平滑。

void insertDepth32f(cv::Mat& depth)
{
    const int width = depth.cols;
    const int height = depth.rows;
    float* data = (float*)depth.data;
    cv::Mat integralMap = cv::Mat::zeros(height, width, CV_64F);
    cv::Mat ptsMap = cv::Mat::zeros(height, width, CV_32S);
    double* integral = (double*)integralMap.data;
    int* ptsIntegral = (int*)ptsMap.data;
    memset(integral, 0, sizeof(double) * width * height);
    memset(ptsIntegral, 0, sizeof(int) * width * height);
    for (int i = 0; i < height; ++i)
    {
        int id1 = i * width;
        for (int j = 0; j < width; ++j)
        {
            int id2 = id1 + j;
            if (data[id2] > 1e-3)
            {
                integral[id2] = data[id2];
                ptsIntegral[id2] = 1;
            }
        }
    }
    // 积分区间
    for (int i = 0; i < height; ++i)
    {
        int id1 = i * width;
        for (int j = 1; j < width; ++j)
        {
            int id2 = id1 + j;
            integral[id2] += integral[id2 - 1];
            ptsIntegral[id2] += ptsIntegral[id2 - 1];
        }
    }
    for (int i = 1; i < height; ++i)
    {
        int id1 = i * width;
        for (int j = 0; j < width; ++j)
        {
            int id2 = id1 + j;
            integral[id2] += integral[id2 - width];
            ptsIntegral[id2] += ptsIntegral[id2 - width];
        }
    }
    int wnd;
    double dWnd = 2;
    while (dWnd > 1)
    {
        wnd = int(dWnd);
        dWnd /= 2;
        for (int i = 0; i < height; ++i)
        {
            int id1 = i * width;
            for (int j = 0; j < width; ++j)
            {
                int id2 = id1 + j;
                int left = j - wnd - 1;
                int right = j + wnd;
                int top = i - wnd - 1;
                int bot = i + wnd;
                left = max(0, left);
                right = min(right, width - 1);
                top = max(0, top);
                bot = min(bot, height - 1);
                int dx = right - left;
                int dy = (bot - top) * width;
                int idLeftTop = top * width + left;
                int idRightTop = idLeftTop + dx;
                int idLeftBot = idLeftTop + dy;
                int idRightBot = idLeftBot + dx;
                int ptsCnt = ptsIntegral[idRightBot] + ptsIntegral[idLeftTop] - (ptsIntegral[idLeftBot] + ptsIntegral[idRightTop]);
                double sumGray = integral[idRightBot] + integral[idLeftTop] - (integral[idLeftBot] + integral[idRightTop]);
                if (ptsCnt <= 0)
                {
                    continue;
                }
                data[id2] = float(sumGray / ptsCnt);
            }
        }
        int s = wnd / 2 * 2 + 1;
        if (s > 201)
        {
            s = 201;
        }
        cv::GaussianBlur(depth, depth, cv::Size(s, s), s, s);
    }
}

利用视差图，进行虚拟视点拟合

正向映射：简单的利用左视点原图和视差图进行视点合成，取每一个像素点处的视差值，然后计算新图像中像素点位置，然后赋值。前向映射，单点赋值代码如下。

#include <iostream>
#include <string>
#include <opencv.hpp>

using namespace std;
using namespace cv;

void main()
{
    string imgPath="data/source_images/teddy/";
    Mat srcImgL=imread(imgPath+"imgL.png");
    Mat dispL=imread(imgPath+"dispL.png",0);
    dispL=dispL/4;
    
    int imgHeight=srcImgL.rows;
    int imgWidth=srcImgL.cols;
    int channels=srcImgL.channels();
    
    Mat dstImgL=Mat::zeros(imgHeight,imgWidth, CV_8UC3);
    
    uchar* pImgDataL=(uchar*)srcImgL.data;
    uchar* pDispDataL=(uchar*)dispL.data;
    uchar* pDstDataL=(uchar*)dstImgL.data;
    
    VideoWriter writer("video.avi", CV_FOURCC('D','I','V','X'), 30, Size(imgWidth, imgHeight), 1);
    int cnt=0;
    int viewCnt=50;
    while (cnt !=4)
    {
        for (int k=0; k<viewCnt; k++)
        {
            dstImgL.setTo(0);
            float interp;
            if (cnt%2==0)    interp=(float)k/viewCnt;
            else  interp=float(viewCnt-k)/viewCnt;
            for (int j=0; j<imgHeight; j++)
            {
                for (int i=0; i<imgWidth; i++)
                {
                    uchar dispL=pDispDataL[j*imgWidth+i];
                    float offsetL=dispL* interp;
                    int idL=(int)(offsetL+0.5);  //计算视差值
                    
                    if (idL+i>=imgWidth) continue;
                    //插值结果
                    int idxResult=(j*imgWidth+i)*channels;
                    int idx=(j*imgWidth+i+idL)*channels;
                    for (int chan=0; chan<channels; chan++)
                    {
                        pDstDataL[idxResult+chan]=pImgDataL[idx+chan];
                    }
                }
            }
            namedWindow("show");
            imshow("show", dstImgL);
            waitKey(10);
            writer<<dstImgL;
        }
        cnt++;
    }
    writer.release();
}

反向映射：先根据左视点视差图生成虚拟视点的视差图，然后反向映射得到每一个像素点在原图像中的浮点位置，利用线性插值获取最终颜色值。（虚拟视点位置视差图没有填充空洞版本），可见有很多裂纹。

void main(void)
{
    string imgPath="data/source_images/teddy/";
    Mat srcImgL=imread(imgPath+"imgL.png");
    Mat dispL=imread(imgPath+"dispL.png",0);
    dispL=dispL/4;
    
    int imgHeight=srcImgL.rows;
    int imgWidth=srcImgL.cols;
    
    Mat dstImgL=Mat::zeros(imgHeight,imgWidth, CV_8UC3);
    Mat dstImg=Mat::zeros(imgHeight,imgWidth, CV_8UC3);
    Mat dstNewDispImg=Mat::zeros(imgHeight,imgWidth, CV_32FC1);
    
    uchar* pImgDataL=(uchar*)srcImgL.data;
    uchar* pDispDataL=(uchar*)dispL.data;
    uchar* pDstDataL=(uchar*)dstImgL.data;
    
    VideoWriter writer("video.avi", CV_FOURCC('D','I','V','X'), 30, Size(imgWidth, imgHeight), 1);
    int cnt=0;
    int viewCnt=50;
    while (cnt!=4)
    {
        float interp;
        for (int k=0; k<viewCnt; k++)
        {
            dstNewDispImg.setTo(255);
            dstImgL.setTo(0);
            if (cnt%2==0)   interp=(float)k/viewCnt;
            else interp=(float)(viewCnt-k)/viewCnt;
            
            obtainNewDispMap(dispL, dstNewDispImg, interp);
            
            float* pNewDispData=(float*)dstNewDispImg.data;
            for (int j=0; j<imgHeight; j++)
            {
                for (int i=0; i<imgWidth; i++)
                {
                    float disp=pNewDispData[j*imgWidth+i];
                    float id=i+disp;
                    
                    int id0=floor(id);
                    int id1=floor(id+1);
                    
                    float weight1=1-(id-id0);
                    float weight2=id-id0;
                    
                    id0=index(id0, imgWidth);
                    id1=index(id1, imgWidth);
                    
                    //插值结果
                    pDstDataL[j*imgWidth*3+i*3+0]=weight1*pImgDataL[j*imgWidth*3+id0*3+0]+weight2*pImgDataL[j*imgWidth*3+id1*3+0];
                    pDstDataL[j*imgWidth*3+i*3+1]=weight1*pImgDataL[j*imgWidth*3+id0*3+1]+weight2*pImgDataL[j*imgWidth*3+id1*3+1];
                    pDstDataL[j*imgWidth*3+i*3+2]=weight1*pImgDataL[j*imgWidth*3+id0*3+2]+weight2*pImgDataL[j*imgWidth*3+id1*3+2];
                }
            }
            namedWindow("virImg");
            imshow("virImg", dstImgL);
            waitKey(10);
            writer<<dstImgL;
        }
        cnt++;
    }
    writer.release();
}

反向映射+空洞填充+双线性插值:上面生成虚拟视点位置的视差图时没有填充空洞，生成的虚拟视点会有很多裂纹存在。加上空洞填充能够有效消除裂纹。

2.5 跟踪匹配

使用特征匹配，对前后帧进行比较
openCV中匹配函数：

Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce");
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("FlannBased");
BFMatcher matcher(NORM_HAMMING);
BFMatcher matcher(NORM_L2);
matcher->knnMatch(first_desc, desc, matches, 2);

opencv跟踪样例

2.6 optical_flow

光流是由对象或相机的移动引起的两个连续帧之间的图像对象的明显运动的模式，它是2d矢量场，其中每个矢量是位移矢量，表示从前一帧到当前帧点的移动。
原理基础：

目标像素强度在连续帧之间不变。
相邻像素具有相似的运动
第一帧的像素I(x,y,t),在dt时间之后的下一帧中移动距离(dx，dy),因为这些像素是相同的,而且亮度不变,所以：

$\\ Taylor级数展开，删除常用项病除以dt得到： \\ I(x,y,t)=I(x+dx,y+dy,t+dt) = I(x,y,t)+\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t+H.O.T \\ H.O.T--高阶无穷小 \\ 由于假设1： \\ \frac{\partial I}{\partial x}\frac{\Delta x}{\Delta t}+\frac{\partial I}{\partial y}\frac{\Delta y}{\Delta t}+\frac{\partial I}{\partial t}\frac{\Delta t}{\Delta t}= f_xu+f_yv+f_t=0 \\ 其中：f_x和f_y分别为图像的梯度，f_t为图像沿着时间的梯度。但是u，v是未知的，lucas-kanade算法就是解决这个问题的。$

nextPts,status,err = cv2.calcOpticalFlowPyrLK(prevImg,   #上一帧图片
                                              nextImg,   #当前帧图片
                                              prevPts,   #上一帧找到的特征点向量 
                                              nextPts    #与返回值中的nextPtrs相同
                                              [, status[, err[, winSize
                                              [, maxLevel[, criteria
                                              [, flags[, minEigThreshold]]]]]]])

(1)LK-lucas kanade稀疏光流跟踪
.Lucas-Kanade 方法需要设置窗口大小,所有9个点都有相同的动作,我们可以找到这9个点的（f x，f y，f t）. 所以现在我们的问题变成解决了9个方程式，其中两个未知变量是过度确定的.解的个数大于未知数的个数，这是个超定方程，使用最小二乘的方法来求解最优值.
$[\begin{matrix} u \\ v\end{matrix}]=[\begin{matrix} \sum_if_{x_i}^2 & \sum_if_{x_i}f_{y_i} \\ \sum_if_{x_i}f_{y_i} & \sum_if_{y_i}^2 \end{matrix}]^{-1}[\begin{matrix}-\sum_if_{x_i}f_{t_i} \\ -\sum_if_{y_i}f_{t_i}\end{matrix}]$

用户给出一些点用来追踪，从而获得点的光流向量。但是有另外一个问题需要解决，目前讨论的运动都是小步长的运动，如果有幅度大的运动出现，本算法就会失效。
解决办法是利用图像金字塔。在金字塔顶端的小尺寸图片当中，大幅度的运动就变成了小幅度的运动。所以使用LK算法，可以得到尺度空间上的光流。

import numpy as np
import cv2
import matplotlib.pyplot as plt

cap = cv2.VideoCapture('1.avi')
# params for ShiTomasi corner detection
feature_params = dict( maxCorners = 100,
                       qualityLevel = 0.3,
                       minDistance = 7,
                       blockSize = 7 )

# Parameters for lucas kanade optical flow
lk_params = dict( winSize  = (15,15),
                  maxLevel = 2,
                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

# Create some random colors
color = np.random.randint(0,255,(100,3))

# Take first frame and find corners in it
ret, old_frame = cap.read()
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
p0 = cv2.goodFeaturesToTrack(old_gray, mask = None, **feature_params)

# Create a mask image for drawing purposes
mask = np.zeros_like(old_frame)

while(1):
    ret,frame = cap.read()
    frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    # calculate optical flow
    p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)

    # Select good points
    good_new = p1[st==1]
    good_old = p0[st==1]

    # draw the tracks
    for i,(new,old) in enumerate(zip(good_new,good_old)):
        a,b = new.ravel()
        c,d = old.ravel()
        mask = cv2.line(mask, (a,b),(c,d), color[i].tolist(), 2)
        frame = cv2.circle(frame,(a,b),5,color[i].tolist(),-1)
    img = cv2.add(frame,mask)

    cv2.imshow('frame',img)
    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break

    # Now update the previous frame and previous points
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1,1,2)

cv2.destroyAllWindows()
cap.release()

（2）稠密光流跟踪

import numpy as np
import cv2
import matplotlib.pyplot as plt

cap = cv2.VideoCapture('1.avi')
ret, frame1 = cap.read()
prvs = cv2.cvtColor(frame1,cv2.COLOR_BGR2GRAY)
hsv = np.zeros_like(frame1)
hsv[...,1] = 255

while(1):
    ret, frame2 = cap.read()
    next = cv2.cvtColor(frame2,cv2.COLOR_BGR2GRAY)

    flow = cv2.calcOpticalFlowFarneback(prvs,next, None, 0.5, 3, 15, 3, 5, 1.2, 0)

    mag, ang = cv2.cartToPolar(flow[...,0], flow[...,1])
    hsv[...,0] = ang*180/np.pi/2
    hsv[...,2] = cv2.normalize(mag,None,0,255,cv2.NORM_MINMAX)
    bgr = cv2.cvtColor(hsv,cv2.COLOR_HSV2BGR)

    cv2.imshow('frame2',bgr)
    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break
    elif k == ord('s'):
        cv2.imwrite('opticalfb.png',frame2)
        cv2.imwrite('opticalhsv.png',bgr)
    prvs = next

cap.release()
cv2.destroyAllWindows()

(3)反向检测光流跟踪法

import numpy as np
import cv2 as cv
import math
cap = cv.VideoCapture('car_flow.mp4')

#角点检测参数
feature_params = dict(maxCorners=100, qualityLevel=0.1, minDistance=7, blockSize=7)

#KLT光流参数
lk_params = dict(winSize=(15, 15), maxLevel=2, criteria=(cv.TERM_CRITERIA_EPS | cv.TERM_CRITERIA_COUNT, 10, 0.02))

height = cap.get(cv.CAP_PROP_FRAME_HEIGHT)
width = cap.get(cv.CAP_PROP_FRAME_WIDTH)
fps = cap.get(cv.CAP_PROP_FPS)
#out = cv.VideoWriter("reslut.avi", cv.VideoWriter_fourcc('D', 'I', 'V', 'X'), fps,
                     #(np.int(width), np.int(height)), True)

tracks = []
track_len = 15
frame_idx = 0
detect_interval = 5
while True:

    ret, frame = cap.read()
    if ret:
        frame_gray = cv.cvtColor(frame, cv.COLOR_BGR2GRAY)
        vis = frame.copy()

        if len(tracks)>0:
            img0 ,img1 = prev_gray, frame_gray
            p0 = np.float32([tr[-1] for tr in tracks]).reshape(-1,1,2)
            # 上一帧的角点和当前帧的图像作为输入来得到角点在当前帧的位置  
            p1, st, err = cv.calcOpticalFlowPyrLK(img0, img1, p0, None, **lk_params)

            # 反向检查,当前帧跟踪到的角点及图像和前一帧的图像作为输入来找到前一帧的角点位置  
            p0r, _, _ = cv.calcOpticalFlowPyrLK(img1, img0, p1, None, **lk_params)

            # 得到角点回溯与前一帧实际角点的位置变化关系 
            d = abs(p0-p0r).reshape(-1,2).max(-1)

            #判断d内的值是否小于1，大于1跟踪被认为是错误的跟踪点
            good = d < 1

            new_tracks = []

            for i, (tr, (x, y), flag) in enumerate(zip(tracks, p1.reshape(-1, 2), good)):

                # 判断是否为正确的跟踪点
                if not flag:
                    continue

                # 存储动态的角点
                tr.append((x, y))

                # 只保留track_len长度的数据，消除掉前面的超出的轨迹
                if len(tr) > track_len:
                    del tr[0]
                # 保存在新的list中
                new_tracks.append(tr)

                cv.circle(vis, (x, y), 2, (0, 255, 0), -1)

            # 更新特征点    
            tracks = new_tracks

            # #以上一振角点为初始点，当前帧跟踪到的点为终点,画出运动轨迹
            cv.polylines(vis, [np.int32(tr) for tr in tracks], False, (0, 255, 0), 1)


        # 每隔 detect_interval 时间检测一次特征点
        if frame_idx % detect_interval==0:
            mask = np.zeros_like(frame_gray)
            mask[:] = 255

            if frame_idx !=0:
                for x,y in [np.int32(tr[-1]) for tr in tracks]:
                    cv.circle(mask, (x, y), 5, 0, -1)

            p = cv.goodFeaturesToTrack(frame_gray, mask=mask, **feature_params)
            if p is not None:
                for x, y in np.float32(p).reshape(-1,2):
                    tracks.append([(x, y)])

        frame_idx += 1
        prev_gray = frame_gray

        cv.imshow('track', vis)
        #out.write(vis)
        ch = cv.waitKey(1)
        if ch ==27:
            cv.imwrite('track.jpg', vis)
            break
    else:
        break

cv.destroyAllWindows()
cap.release()

（4）删除静止点的光流分析

import numpy as np
import cv2 as cv
import math
cap = cv.VideoCapture('car_flow.mp4')

#角点检测参数
feature_params = dict(maxCorners=100, qualityLevel=0.1, minDistance=7, blockSize=7)

#KLT光流参数
lk_params = dict(winSize=(15, 15), maxLevel=2, criteria=(cv.TERM_CRITERIA_EPS | cv.TERM_CRITERIA_COUNT, 10, 0.02))


# 随机颜色
color = np.random.randint(0,255,(100,3))

# 读取第一帧
ret, old_frame = cap.read()
old_gray = cv.cvtColor(old_frame, cv.COLOR_BGR2GRAY)
p0 = cv.goodFeaturesToTrack(old_gray, mask=None, **feature_params,useHarrisDetector=False,k=0.04)
good_ini=p0.copy()

def caldist(a,b,c,d):
    return abs(a-c)+abs(b-d)

mask = np.zeros_like(old_frame)
# 光流跟踪
while True:
    ret, frame = cap.read()
    if ret is False:
        break
    frame_gray = cv.cvtColor(frame, cv.COLOR_BGR2GRAY)
    # 计算光流
    p1, st, err = cv.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
    # 根据状态选择

    good_new = p1[st == 1]
    good_old = p0[st == 1]
    #删除静止点
    k=0
    for i, (new0, old0) in enumerate(zip(good_new,good_old)):
        a0,b0 = new0.ravel()
        c0,d0 = old0.ravel()
        dist=caldist(a0,b0,c0,d0)
        if dist>2:
            good_new[k]=good_new[i]
            good_old[k]=good_old[i]
            good_ini[k]=good_ini[i]
            k=k+1

    # 提取动态点
    good_ini=good_ini[:k]
    good_new=good_new[:k]
    good_old=good_old[:k]  

    # 绘制跟踪线
    for i, (new, old) in enumerate(zip(good_new,good_old)):
        a,b = new.ravel()
        c,d = old.ravel()
        mask = cv.line(mask, (a,b),(c,d), color[i].tolist(), 2)
        frame = cv.circle(frame,(a,b),5,color[i].tolist(),-1)
    cv.imshow('frame',cv.add(frame,mask))
    k = cv.waitKey(30) & 0xff
    if k == 27:
        cv.imwrite("flow.jpg", cv.add(frame,mask))
        break

    # 更新
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1, 1, 2)

    if good_ini.shape[0]<40:
        p0 = cv.goodFeaturesToTrack(old_gray, mask=None, **feature_params)
        good_ini=p0.copy()

cv.destroyAllWindows()
cap.release()

optical_flow

谢欣燕

关注

4
点赞
踩
20

收藏

觉得还不错? 一键收藏
0
评论
图像笔记（二）——特征提取及跟踪匹配

1.图像目标检测2.图像特征提取及匹配2.1 基于灰度模板匹配算法Blocking Matching 是根据移植模板图到另一幅图像中寻找与模板图像相似的自图像。基于灰度图像的匹配算法也称作相关匹配算法，用空间二维滑动模板进行匹配，不同的匹配算法相关准则不同。2.1.1 平均绝对差算法（MAD）思路简单匹配精度较高D(i,j)=1M∗N∑S=1m∑t=1N∣S(i+s−1,j+t−1)−T(s,t)∣m−M+1≥i≥1，n−N+1≥j≥1D(i,j)=\frac{1}{M*N}\sum^m_
复制链接

扫一扫