mean shift and openCV2--01 : implementation of meanShift

  In the last post(mean shift and openCV2--00) I try to track the region of interest by mean-shift and Camshift, but I haven't tried to understand how to implement mean-shift. In this post, I would like to analyze the implementation of the mean-shift from openCV2(openCV2 mean-shift).

  Mean-shift start from an initial location and iteratively around,  find the centroid location and repeats the procedure until the window center converges to a stable point. The new centroid could be found by the equation(1) come from wikipedia(mean-shift equation).(2) is the kernal function, which determine the weight of nearby points.OpenCV2 do not choose gaussian distribution as the kernel but using spatial moments as kernal.

    Here is the core of the mean-shift.

for( i = 0; i < niters; i++ )
{
   //step 1 : make sure cur_rect has reasonable position and size
   cur_rect = cur_rect & Rect(0, 0, mat.cols, mat.rows); #1
   if( cur_rect == Rect() )
   {
      cur_rect.x = mat.cols/2;
      cur_rect.y = mat.rows/2;
   }
   cur_rect.width = std::max(cur_rect.width, 1);
   cur_rect.height = std::max(cur_rect.height, 1);
        
   //step 2 : calculate the moments
   Moments m = moments(mat(cur_rect));
        
   //step 3 : Calculating center of mass
   if( fabs(m.m00) < DBL_EPSILON ) //#2
      break;
        
   int dx = cvRound( m.m10/m.m00 - window.width*0.5 ); //#3
   int dy = cvRound( m.m01/m.m00 - window.height*0.5 ); 
        
   int nx = std::min(std::max(cur_rect.x + dx, 0), 
                     mat.cols - cur_rect.width); //#4
   int ny = std::min(std::max(cur_rect.y + dy, 0), 
                     mat.rows - cur_rect.height); 
        
   dx = nx - cur_rect.x; //#5
   dy = ny - cur_rect.y;
   cur_rect.x = nx; #6
   cur_rect.y = ny;
        
   // step 4 : Check for coverage centers mass & window
   if( dx*dx + dy*dy < eps ) //#7
      break;
}

Anatomy of codes 

  The purpose of #1 is make sure the cur_rect(initial position when i == 0) is in the range of the input cv::Mat(mat). If the cur_rect do not in the range of the mat, then cur_rect will equal to Rect().

   Purpose of #2 is make sure the m00 big enough.m00, m01, m10 are found by the equation(3).

  If we express it by codes, they would be the same as

template<typename Func>
float moment(cv::Mat const &input, Func func)
{
    float sum = 0;
    for(int row = 0; row != input.rows; ++row){
        auto ptr = input.ptr<float>(row);
        for(int col = 0; col != input.cols; ++col){
            sum += func(ptr[col], row, col);
        }
    }

    return sum;
}

int main()
{
    cv::Mat input = (cv::Mat_(3,3) <<
                     1, 2, 3,
                     4, 5, 6,
                     7, 8, 9);

   float const m00 = moment(input, [](float data, int, int)
   {
      return data;
   });

   float const m10 = moment(input, [](float data, int, int col)
   {
      return data * col;
   });

   float const m01 = moment(input, [](float data, int row, int)
   {
      return data * row;
   });
}

#3 estimate how far the horizontal position of the new centroid could be.

#4 measure the new centroid position of horizontal.

std::max(cur_rect.x + dx, 0)

make sure the new centroid horizontal position would not smaller than 0.

std::min(std::max(cur_rect.x + dx, 0), 
                     mat.cols - cur_rect.width);

outer std::min make sure the new centroid horizontal position greater or equal to mat.cols - cur_rect.width.

#5 is the distance between new centroid horizontal position and old centroid horizontal position.

#6 is the new horizontal position of centroid(the center of the window).

#7 measure the distance between old centroid and new centroid are small enough or not, if they
are small enough, the program will end the iteration, else continue until one of the criteria meet(
iteration counts or the distance between old centroid and new centroid).The distance is measure by
euclidean distance(4).The program don't need to take the square root of (dx * dx) + (dy * dy), because the eps already multiply with itself(eps = cvRound(eps*eps)).