Harris Corners--00 : principle

    Harris Corner is an famous algorithm to locate the useful features(corners) in the image.Now there are two questions for us to understand.

1 : What are features?   

• Exactly, there are no formal definition of feature yet(atleast I can't find it), in the world of computer vision, you could treat it as the information you use for tracking, discern the objects(just my intuition, if you find any errors, please give me some comments).

2 : Why choose corners as features?

• Because in many cases, corners are good features.They provide enough messages for the program to discern the location of the objects.

3 : How do we find the corners?
     This is the main part of this post, now let me try to explain it.

In a nutshell, we use the lambdas of Hessian Matrix(1) to measure the change rate of the vectors.There are three conditions

1. lambda one and lambda two are small, that means this pixel belong to a flat region.
2. lambda one or lambda two are big, that means this pixel belong to edges.
3. lambda one and lambda two are big, that means this pixel belong to corners.

So far,so good, but why do this work? From graph_00, we discovered that the region around the corners are changing rapidly, that means if the pixel belong the the corner, the pixels surrounded it have a high chance to change rapidly for every direction.

graph_00

Computation of the lambdas are very expensive, so Harris and Stephens proposed a function Mc(2) to approximate the it.

graph_01(I take it from the original paper) explain what are the Mc doing.We could decide which Mc qualify the condition as corner.

graph_01

4 : Where are the Hessian matrix come from?

The Hessian Matrix are deduced from the sum of square difference(3), (4) is approximated by the taylor series(only take the zero and first derivative terms).Substitude (4) to (3), we get (5), and we can express (5) to the new equation--(6). The A of (6) is the Harris matrix we care about.

Segmenting leaf by openCV2--02 : Segment leafs to single leaf part 1

  Last post I  use watershed to segment the leafs, this post I would like to segment those leafs into single leaf.Let us begin from simple case(graph_00).graph_00 is easier to segment than graph_01 and graph_02 because the leafs are not stick together as close as graph_01 and graph_02. If they stick(merge) together like graph_01 and graph_02 do, we will need to pay much more efforts to separate those leafs.

graph_00

graph_01

graph_02

  Separate the leafs of graph_00 could be done by simple algorithms as following.

 

/**
 * @brief cut the leafs to single leaf
 * @param input  input image
 * @param name   name of the image want to be saved
 */
void cut_to_single_leaf_simple(cv::Mat const &input, 
                               std::string const &name)
{
    cv::Mat mask;
    //step 1 : convert the input to gray image
    cv::cvtColor(input, mask, CV_BGR2GRAY);
    //step 2 : binarize it
    cv::threshold(mask, mask, 128, 255, CV_THRESH_BINARY);
    //step 3 : do serious erosion
    cv::erode(mask, mask, 
             cv::getStructuringElement(cv::MORPH_RECT, 
             cv::Size(5, 5)), 
             cv::Point(-1, -1), 5);
    //step 4 : find and remove the outliers contours
    ContoursType contours;
    cv::findContours(mask, contours, CV_RETR_EXTERNAL, 
                     CV_CHAIN_APPROX_NONE);
    OCV::remove_contours(contours, 500, 120000);
    cv::Mat result;
    //step 5 : draw the contours one by one
    for(size_t i = 0; i != contours.size(); ++i){
        mask.setTo(0);
        cv::drawContours(mask, contours, i, cv::Scalar(255), 
                         CV_FILLED);
        //step 6 : dilate and close the contours
        cv::dilate(mask, mask, 
                   cv::getStructuringElement(cv::MORPH_RECT, 
                   cv::Size(5, 5)), 
                   cv::Point(-1, -1), 6);
        cv::morphologyEx(mask, mask, CV_MOP_CLOSE, 
                         cv::getStructuringElement(cv::MORPH_RECT, 
                         cv::Size(15, 15)));
        //step 6 : save the results
        input.copyTo(result, mask);
        cv::imwrite(name + std::to_string(i) + ".png", result);
        result.setTo(0);
    }
}

I think the codes and the comment already show us the whole story.

As usual, the codes can download from github.

Next post of Recognize leaf by openCV2 

  Try my best to segment the graph_01 and graph_02, I think I would use machine learning to solve this problem.

Segmenting leaf by openCV2--01 : By watershed and contours

In the last post(recognize leaf--00) I use threshold and contours to cut the leafs from the images, this time I will give watershed a try.

graph_00

graph_01

graph_02

Step 0

  Read the images(graph_00~graph_02) and convert it to gray image.

cv::Mat color_image = cv::imread(Folder + "leaf" + num + ".png");
if(color_image.empty()){
  std::cerr<<"input is empty"<<std::endl;
  return -1;
}

cv::Mat image;
cv::cvtColor(color_image, image, CV_BGR2GRAY);

Step 1

Binarize the image.

cv::Mat binary;
cv::threshold(image, binary, 140, 255, cv::THRESH_BINARY); 

Step 2

Remove small objects and noise.

auto const structure = cv::getStructuringElement(cv::MORPH_ELLIPSE, 
                       cv::Size(5, 5));
cv::Mat fore_ground;
cv::erode(binary, fore_ground, structure, cv::Point(-1, -1), 4);

We severely erode the images in order to retain only pixels belonging to the foreground.

Step 3

Some of the objects can't removed by erosion, so we find out the contours of the image and remove unimportant objects from the images in order to get the foreground(255), they will be considered to correspond to objects of interest.

std::vector<std::vector<cv::Point>> contours;
cv::findContours(fore_ground, contours, CV_RETR_EXTERNAL, 
      CV_CHAIN_APPROX_NONE);
OCV::remove_contours(contours, 8000, 50000);
fore_ground.setTo(0);
cv::drawContours(fore_ground, contours, -1, cv::Scalar(255), CV_FILLED);

remove_contours was explained in the previous post.

graph_03

graph_04

graph_05

Step 4

  Mark the background as 128 and mark unknown pixels(neither foreground nor background) as 0.

cv::Mat back_ground;
cv::dilate(binary, back_ground, structure, cv::Point(-1, -1), 4);
cv::threshold(back_ground, back_ground, 1, 128, cv::THRESH_BINARY_INV);

graph_06

graph_07

graph_08

Step 5

  Combined foreground and background to get marker and use the marker to get the mask we need.

cv::Mat markers = back_ground + fore_ground;       
cv::Mat mask;
markers.convertTo(mask, CV_32S);
cv::watershed(color_image, mask);
mask.convertTo(mask, CV_8U);
cv::threshold(mask, mask, 150, 255, CV_THRESH_BINARY);

graph_09

graph_10

graph_11

Step 6

Generate final result.

cv::Mat result;
color_image.copyTo(result, mask);

graph_12

graph_13

graph_14

As usual, the codes can download from github.

Next post of Recognize leaf by openCV2

 Next post I will try to use K-means to cluster the images. I decided to put k-means as a new topic, next post of this series I will try to segment the leafs into independent leaf.

Segmenting leaf by openCV2--00 : By threshold and findContours

  These series of post will try to cut the leaf from the images(graph_00~graph_02) with different solutions.The first solution is based on threshold and contours.

graph_00

graph_01

graph_02

Step 0

read the image and covert it to gray image

cv::Mat color_image = cv::imread(Folder + "leaf" + num + ".png");
if(color_image.empty()){
   std::cerr<<"input is empty"<<std::endl;
   return -1;
}

cv::Mat image;
cv::cvtColor(color_image, image, CV_BGR2GRAY);

Step 1

 
  Apply threshold.

 

cv::Mat binary;
cv::threshold(image, binary, 140, 255, cv::THRESH_BINARY);

graph_03

graph_04

graph_05

Step 2 

Use morphology operations(erosion and close) to remove discontinuous region(erosion) and close the holes(close).

auto const structure = cv::getStructuringElement(cv::MORPH_ELLIPSE, 
                                                 cv::Size(5, 5));
cv::erode(binary, binary, structure);
cv::morphologyEx(binary, binary, cv::MORPH_CLOSE, structure);

graph_06

graph_07

graph_08

 Step 3

find contours of the images.

std::vector<std::vector<cv::Point>> contours;
cv::findContours(binary, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);

graph_09

graph_10

graph_11

Step 4

From step 3, we could see that graph_09~graph_11 still contains a lot of undesired blobs, we could remove them by the area of contours(the pixels in the contours).

void remove_contours(std::vector<std::vector<cv::Point>>&contours, 
                     double cmin, double cmax)
{
    auto it = std::partition(std::begin(contours), std::end(contours), 
    [=](std::vector const &data)
    {
       auto const size = cv::contourArea(data);
       return !(size < cmin || size > cmax);
    });
    contours.erase(it, std::end(contours));
} 

After that, we invert the mask(change 0 to 255, 255 to 0).

graph_12

graph_13

graph_14

Final results

 Use the masks(graph_12~graph_14) to locate the region of interset, we get the results.

graph_15

graph_16

graph_17

 As usual, the codes can download from github.

Next post of Recognize leaf by openCV2

 Next post I will try to use watershed to segment the images, maybe it could give us better results.

Hough transform and openCV2--00 : implementation of standard hough-transform

    OpenCV2 provide different implementation of Hough-transform, in this post I want to analyze the implementation of the standard Hough-transform--HoughLinesStandard.You could find the codes of HoughLinesStandard at here(houghLinesStandard).

   Before we step into the codes of HoughLinesStandard, better take a look at the tutorial of openCV2(tutorial).

Codes 

struct LinePolar
{
    float rho;
    float angle;
};


struct hough_cmp_gt
{
    hough_cmp_gt(const int* _aux) : aux(_aux) {}
    bool operator()(int l1, int l2) const
    {
        return aux[l1] > aux[l2] || (aux[l1] == aux[l2] && l1 < l2);
    }
    const int* aux;
};

static void
HoughLinesStandard( const Mat& img, float rho, float theta,
                    int threshold, std::vector& lines, int linesMax )
{
    int i, j; //#ra0
    float irho = 1 / rho;

    CV_Assert( img.type() == CV_8UC1 ); //#ra1

    const uchar* image = img.data; //#rb1
    int step = (int)img.step; //#rb2
    int width = img.cols; //#ra2
    int height = img.rows; //#ra3

    int numangle = cvRound(CV_PI / theta);
    int numrho = cvRound(((width + height) * 2 + 1) / rho); //#ra4, #1

    AutoBuffer _accum((numangle+2) * (numrho+2)); //#ra5
    std::vector _sort_buf; //#ra6
    AutoBuffer _tabSin(numangle);
    AutoBuffer _tabCos(numangle);
    int *accum = _accum; //#ra7
    float *tabSin = _tabSin, *tabCos = _tabCos;

    memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );

    float ang = 0;
    for(int n = 0; n < numangle; ang += theta, n++ ) //#rd1
    {
        tabSin[n] = (float)(sin((double)ang) * irho); //#rc1
        tabCos[n] = (float)(cos((double)ang) * irho); //#rc2
    }

    // stage 1. fill accumulator
    for( i = 0; i < height; i++ ) //#rd2
        for( j = 0; j < width; j++ ) //#rd3
        {
            if( image[i * step + j] != 0 ) //#rb3
                for(int n = 0; n < numangle; n++ ) //#rd4
                {
                    int r = cvRound( j * tabCos[n] + 
                                   i * tabSin[n] ); //#2a
                    r += (numrho - 1) / 2; //#2b
                    accum[(n+1) * (numrho+2) + r+1]++; //#rd5, #3
                }
        }

    // stage 2. find local maximums
    for(int r = 0; r < numrho; r++ ) //#rd6
        for(int n = 0; n < numangle; n++ ) //#rd7
        {
            int base = (n+1) * (numrho+2) + r+1;
            if( accum[base] > threshold &&
                accum[base] > accum[base - 1] && 
                accum[base] >= accum[base + 1] &&
                accum[base] > accum[base - numrho - 2] && 
                accum[base] >= accum[base + numrho + 2] )
                _sort_buf.push_back(base);
        }

    // stage 3. sort the detected lines by accumulator value
    std::sort(_sort_buf.begin(), _sort_buf.end(), hough_cmp_gt(accum));

    // stage 4. store the first min(total,linesMax) lines to 
    // the output buffer
    linesMax = std::min(linesMax, (int)_sort_buf.size()); //#re1
    double scale = 1./(numrho+2);
    for( i = 0; i < linesMax; i++ ) //#rd8
    {
        LinePolar line; //#ra8
        int idx = _sort_buf[i];
        int n = cvFloor(idx*scale) - 1; //#4
        int r = idx - (n+1)*(numrho+2) - 1; //#5
        line.rho = (r - (numrho - 1)*0.5f) * rho; //#6
        line.angle = n * theta;
        lines.push_back(Vec2f(line.rho, line.angle)); //#rf1
    }
}

 

Anatomy of codes

I do not sure that my anatomy are 100% correct, if you find out anything wrong, please give me some advices, thanks a lot.

 #1 

numrho is the number of distance.According to Digital Image Processing 3rd, the range of the distance is between equation(1).

According to equation(1) and Pythagoras' Theorem, we know that the number of distance should have 2D + 1 units(at least).

It is 2D + 1 but not 2D because we need to include 0.

 

#2b 

r need to do some offset to make sure the value is positive.

#3 

 accum[(n + 1) * (numrho + 2) + r + 1]++;

(numrho + 2) is the number of columns of the accumalate table

(n + 1) make sure the first row would not be used

 (r + 1) make sure the first column would not be used

that is why the dimensions of the accumulate table(_accum) are (numangle+2) * (numrho+2)

because it need two more rows and columns to avoid the pointer access to invalid data.

#4 

 n is the row number of the table, which corresponding to the theta

#5 

r is same as the r at #3, from #3 we know that   
idx == accum[(n+1) * (numrho+2) + r+1]  => r = idx -(n+1) * (numrho+2) -1  

#6 

From #5, we regain the r same as #2b, but the r we need is not the #2b provided but the same as #2a, so we need to subtract it.

Refine codes 

   Thanks to openCV2, we know the details of how to translate the algorithm to codes, this save us a lot of troubles.But I believe we could make the codes become better.

  #ra1~#ra8

There are one technique which called lazy evaluation. In short, we could declare the variables we need just in time;In other words, we better avoid to declare the variables before we need it.#ra1~#ra7 declare the variables too early.

example : 

we could check the type before we declare any variable

original

int i, j;

float irho = 1.0/rho;

CV_Assert(img.type() == CV_8U);

after refine

CV_Assert(img.type() == CV_8U);

float const irho = 1.0/rho;

#rb1~#rb3  

We could access the pixels with more efficient solution

 original

const uchar* image = img.data; //#rb1

int step = (int)img.step; //#rb2

//.....a lot of codes 

 for( i = 0; i < height; i++ )

        for( j = 0; j < width; j++ ) 

        {             

          if( image[i * step + j] != 0 ) 

        } 

after refine

 for(int i = 0; i < height; ++i){

        uchar const *img_ptr = img.ptr<uchar>(i);

        for(int j = 0; j < width; ++j){

            //access the pixel with cheaper operation

            //original one need (i * img.step + j)

            if( img_ptr[j] != 0 )

 #rc1~#rc2

c-style sin only provide us the api which accept "double" because c do not support overload(this suck, really).That means we need to convert the float to double and convert it back float, this is

a waste.If we use std::sin this problem would not exist anymore, because c++ provide different overload for std::sin.

  original

tabSin[n] = (float)(sin(double(ang)) * irho);

after refine

tabSin[n] = std::sin(ang) * irho;

#rd1~#rd7 

We could declare the variable of for loop just in time, this way the codes would be easier to maintain and read.

  original

int i, j;

for(i = 0; i != 10; ++i)

after refine

for(int i = 0; i != 10; ++i)

 #re1

Use static_cast to replace old type casting.

 original

linesMax = std::min(linesMax, (int)_sort_buf.size());

after refine

linesMax = std::min(linesMax, static_cast<int>(_sort_buf.size()));

 #rf1

Call resize and use random access to replace push_back, this could save the cost of memory reallocation and the cost of checking the memory is big enough or not when you call push_back.

 

 original

lines.push_back(cv::Vec2f(lines.rho, lines.angle));

after refine

lines.resize(linesMax);

//.........

lines[i] = cv::Vec2f(lines.rho, lines.angle); 

Whatever, this is just my view about the codes, if you find out anything wrong, please give me some comments, thanks. 

As always, the codes(houghTransform/standardHough) can download from github

otsu threshold and openCV2--00 : implementation of otsu threshold

    First of all, otsu threshold was proposed by Nobuyuki otsu in their paper.This post will talk about the implementation of the otsu threshold of openCV2.

     Following post will introduce the basic idea of otsu threshold based on the text book Digital image processing 3rd, this could help us understand what are the codes trying to do, I would not try to explain the details of otsu threshold becaus there are many good resouces to study, you should read the paper or find a good text book like Digital image processing 3rd if you want to know the whole story.

     Objective of otsu threshold is maximize the between-class variance(9).The basic idea is simple, if the threshold is good enough, the intensity of different class should be pretty different;In other words, if we could make the intensities between classes have biggest different, the threshold should be the best one.

    Class is the component of pixels within specific range.In otsu threshold, we classify the pixels into two classes.If C1 range is within [0, k], then C2 is within [k + 1, L - 1](assume the pixels intensity is within [0,  L - 1]).
   

   I want to give some explanation of the equations which will be used by the implementation of openCV2, without some investigation to these equations, it is almost impossible to understand what are the codes doing.

(1) is the total probability of C1

(2) is the total probability of C2.

(3) and (4) is the mean intensity of C1 and C2.

(5) is the global mean of the whole image.

(7) is deduced by (6) and will be used to compute (4).

(8) is the total variance combine by two classes

(9) is the  between-class variance.

  

    Now we should be able to study the codes.The implementation of the otsu threshold of openCV2 could be found on otsu threshold.Below is the code from openCV2(I do some small modification to it ).

static double
getThreshVal_Otsu_8u( const Mat& _src )
{
    Size size = _src.size();
    if( _src.isContinuous() )
    {
        size.width *= size.height;
        size.height = 1;
    }
    int const N = 256;
    int h[N] = {0}; //histogram
    for(int i = 0; i < size.height; ++i ) //#1
    {
        uchar const *src = _src.data + _src.step*i;
        
        for(int j = 0; j < size.width; ++j )
            ++h[src[j]];
    }

    double mu = 0, scale = 1./(size.width*size.height);
    for(int i = 0; i < N; ++i )
        mu += i * static_cast<double>(h[i]);

    mu *= scale;  //#2
    double mu1 = 0, q1 = 0;
    double max_sigma = 0, max_val = 0;

    for(int i = 0; i < N; ++i )
    {       
        double const p_i = h[i] * scale; //#3   
        q1 += p_i; //#4
        double const q2 = 1. - q1; //#5

        
        if( std::min(q1,q2) < FLT_EPSILON || 
            std::max(q1,q2) > 1. - FLT_EPSILON )
            continue;  //#6

        mu1 = (mu1 * q1 + i * p_i) / q1; //#7
        double const mu2 = (mu - q1*mu1)/q2; //#8
        double const sigma = q1 * q2 *
             (mu1 - mu2)*(mu1 - mu2); //#9
        if( sigma > max_sigma )
        {
            max_sigma = sigma;
            max_val = i;
        }
    }

    return max_val;
}

 

The algorithm of openCV2 adopt is pretty straigtforward.  

Step 1 : build the histogram of the 8bit input image  

Step 2 : find out global mean of equation(4)  

Step 3 : iterate i from 0~255, classify the classes from [0, i], [i, L-1]. Find out the i which make the equation(9) become biggest.  
 

Anatomy of codes
 

  The loop of #1 build the histogram.  

  mu(#2) equal to equation(5).

  #3 is the probability of pixel i.

  #4 is equal to equation(1).

  #5 use equation(8) to find out the result of equation(2).

  #6 pass the loop if one of the class probability is too small, this make sense since m1-m2
   would be very small in this case. If you wonder why do the codes check the min and max but
   not just checking min or max, please give this site(otsu problem) a visit, I have the same
   question too.

   #7 find out the result of equation(3)

   #8 equal to equation(7)

   #9 equal to equation(9)

Next post of otsu
 

  There are an extension of otsu--multi otsu method, I would like to implement it by openCV2 next time(openCV2 haven't implemented this algorithm yet when I write this post).