H-S光流法

一、实验目的

计算出图像序列的光流场，并用有向箭头在图中标出。

二、算法实现

1、构建图像金字塔

金字塔的每一层包括：图像的灰度图、x和y方向的梯度图、相邻两帧灰度图的像素灰度差( $I_t$ )。

a、计算图像x、y、t方向的梯度。 $I_x、I_y$ 计算方法：每个像素沿x、y方向的两个相邻像素的灰度差的一半。 $I_t$ 计算方法：当前帧与上一帧对应位置的像素灰度差。代码如下：

//x方向梯度
void GradImg::calculateDerivativeX()
{
  gradx = cv::Mat::zeros(im.size(),CV_64F);  
  for(int y = 0; y < im.rows; ++y)
  {
    for(int x = 0; x < im.cols; ++x)
    {
      int prev = std::max(x - 1, 0);
      int next = std::min(x + 1, im.cols - 1);
      gradx.at<double>(y, x) = (double) (im.at<double>(y, next) - im.at<double>(y, prev)) * 0.5f;
    }
  }
}

//y方向梯度
void GradImg::calculateDerivativeY()
{
  grady = cv::Mat::zeros(im.size(),CV_64F);  
  for(int y = 0; y < im.rows; ++y)
  {
    for(int x = 0; x < im.cols; ++x)
    {
      int prev = std::max(y - 1, 0);
      int next = std::min(y + 1, im.rows - 1);
      grady.at<double>(y, x) = (double) (im.at<double>(next, x) - im.at<double>(prev, x)) * 0.5f;
    }
  }
}
//t方向梯度
void GradImg::calculateDerivativeT(GradImg& ref)
{
  gradt = cv::Mat::zeros(im.size(),CV_64F);  
  for(int y = 0; y < gradt.rows; ++y)
  {
    for(int x = 0; x < gradt.cols; ++x)
    {
      gradt.at<double>(y, x) = (double) (im.at<double>(y, x) - ref.getImg().at<double>(y, x));
    }
  }
}

b、计算图像金字塔，下采样方法为均匀采样，代码如下：

//均匀采样
cv::Mat HSOpticalFlow::pyrDownMeanSmooth(const cv::Mat& in)
{
  cv::Mat out;
  out = cv::Mat::zeros(cv::Size(in.size().width / 2, in.size().height / 2),in.type()); 
  for(int y = 0; y < out.rows; ++y)
  {
    for(int x = 0; x < out.cols; ++x)
    {
      int x0 = x * 2;
      int x1 = x0 + 1;
      int y0 = y * 2;
      int y1 = y0 + 1;
      out.at<double>(y, x) = (double) ( (in.at<double>(y0, x0) + in.at<double>(y0, x1) + in.at<double>(y1, x0) + in.at<double>(y1, x1)) / 4.0f );      
    }
  }
  return out.clone();
}

//建图像金字塔
void HSOpticalFlow::buildpyr()
{
  //清空存放金字塔的容器
  curlevels.clear();
  reflevels.clear();

  reflevels.push_back(ref);
  curlevels.push_back(cur);
  cv::Mat refcopy;
  cv::Mat curcopy;
  refcopy = ref.getImg();
  curcopy = cur.getImg();
  GradImg gradref;
  GradImg gradcur;
  for (int i = 0 ; i < py_layers ; i++)
  {
    cv::Mat refdown = pyrDownMeanSmooth(refcopy);
    gradref.setImg(refdown);
    reflevels.push_back(gradref);
    refcopy = refdown.clone();
    cv::Mat curdown = pyrDownMeanSmooth(curcopy);

    gradcur.setImg(curdown);
    gradcur.calculateDerivativeX();
    gradcur.calculateDerivativeY();
    gradcur.calculateDerivativeT(gradref);
    curlevels.push_back(gradcur);
    curcopy = curdown.clone();
  }
}

2、计算光流

沿着金字塔，从上往下，依次计算光流。上一层的光流通过插值法传到下一层，作为下一层光流的初值。
这样做的好处：由于光流法基于灰度一致性假设，只能处理帧间运动较小的情况，使用金字塔，可以由粗到精计算大运动情况下的光流。
每个像素的光流计算方法如下：

u = u ⎯ ⎯ - I x I x u ⎯ ⎯ + I y v ⎯ ⎯ - I t a + I 2 x + I 2 y v = v ⎯ ⎯ - I y I x u ⎯ ⎯ + I y v ⎯ ⎯ - I t a + I 2 x + I 2 y

$u=\overline{u}-I_x\frac{I_x\overline{u}+I_y\overline{v}-I_t}{a+I_x^2+I_y^2}\\ v=\overline{v}-I_y\frac{I_x\overline{u}+I_y\overline{v}-I_t}{a+I_x^2+I_y^2}$

其中 $a$ 为常数，反映了对图像数据及平滑约束的可信度，当图像数据本身含有较大噪声时，此时需要加大 $a$ 的值，相反，当输入图像含有较少的噪声时，此时可减小a<script type="math/tex" id="MathJax-Element-177">a</script>的值。
通过在每一层使用上式迭代计算每个像素的光流，再将该层光流传到下一层，作为下一层光流初值，计算下一层的光流。代码如下：

while(iter < maxiter) 
{
   iter++; 
   for(int y = 0; y < im.rows; ++y)
   {
      for(int x = 0; x < im.cols; ++x)
      {
      double mm = optow_ub.at<double>(y,x)*dx.at<double>(y,x) +optow_vb.at<double>         (y,x)*dy.at<double>(y,x) + dt.at<double>(y,x);
      double nn = thea + dx.at<double>(y,x)*dx.at<double>(y,x) + dy.at<double>(y,x)*dy.at<double>(y,x);
      optow_u.at<double>(y,x) = optow_ub.at<double>(y,x) - dx.at<double>(y,x)*(mm/nn);
      optow_v.at<double>(y,x) = optow_vb.at<double>(y,x) - dy.at<double>(y,x)*(mm/nn);    
    }}
    if ( iter < maxiter )  
    {
      optow_ub = smooth<double>(optow_u);
      optow_vb = smooth<double>(optow_v);
    }}

3、光流可视化

用箭头表示图像像素的光流，箭头方向表示光流方向，长度表示光流大小。设定一个阈值区间，只显示在该区间内的光流。代码如下：

//在图片的相应像素上画出光流箭头
for (int y = 0 ; y < u.rows ; y++)
for (int x = 0 ; x < u.cols ; x++)
  {
   if ((abs(u.at<double>(y,x))>th1 && abs(u.at<double>(y,x))<th2) && 
       (abs(v.at<double>(y,x))>th1 && abs(v.at<double>(y,x))<th2))
    {
      cv::Point pStart = Point(x, y);
      cv::Point pEnd = Point(x+u.at<double>(y,x) *1.8f, y+v.at<double>(y,x) *1.8f);   
      drawArrow(rgbcur, pStart, pEnd);
    }
  }
cv::imshow ( rgb_file, rgbcur.clone() );