Revision history [back]

1-Execute four times dilate using a cross operator of size 3 2-Execute seven times erode using same cross

1-Execute four times dilate using a cross operator of size 3

2-Execute seven times erode using same cross

0- binarize image 1-Execute four times dilate using a cross operator of size 3

~~2-Execute~~ 2-With previous result Execute seven times erode using same cross

3- With previous result execute 9 times erode using same cross

0- binarize image

1-Execute four times dilate using a cross operator of size 3

2-With previous result Execute seven times erode using same cross

3- With previous result execute 9 times erode using same cross

0- binarize image

1-Execute four times ~~dilate~~ erode using a cross operator of size 3 5

2-With previous result Execute ~~seven~~ 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

~~3- With previous~~ 4- inverse result ~~execute 9 times erode using same cross~~

and look for contour and filter results..

with this program :

int main (int argc,char **argv)
{
Mat     statComposante; /*< http://docs.opencv.org/trunk/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html */
Mat     centreGComposante; 

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 5, 5 ),cv::Point( 2, 2 ) );


Mat  m=imread("outputImage.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
nbIter=8;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
nbIter=4;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
    mu[i] = moments( contours[i], false ); 
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mu[i].m00>800)
    cout <<"cut in : ("<< mu[i].m10/mu[i].m00 <<","<<mu[i].m01/mu[i].m00<< ")\n";
}

imwrite("test.jpg",m4);
return 0;
};

0- binarize image

1-Execute four times erode using a cross operator of size 5

2-With previous result Execute 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

4- inverse result and look for contour and filter results..

One cut is missing but I think you can start with this program :

int main (int argc,char **argv)
{
Mat     statComposante; /*< http://docs.opencv.org/trunk/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html */
Mat     centreGComposante; 

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 5, 5 ),cv::Point( 2, 2 ) );


Mat  m=imread("outputImage.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
nbIter=8;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
nbIter=4;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
    mu[i] = moments( contours[i], false ); 
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mu[i].m00>800)
    cout <<"cut in : ("<< mu[i].m10/mu[i].m00 <<","<<mu[i].m01/mu[i].m00<< ")\n";
}

imwrite("test.jpg",m4);
return 0;
};

0- binarize image

1-Execute four times erode using a cross operator of size 5

2-With previous result Execute 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

4- inverse result and look for contour and filter results..

One cut is missing but I think you can start with this program :

int main (int argc,char **argv)
{
Mat     statComposante; /*< http://docs.opencv.org/trunk/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html */
Mat     centreGComposante; 

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 5, 5 7, 7 ),cv::Point( 2, 2 3, 3 ) );


Mat  m=imread("outputImage.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
nbIter=8;
nbIter=7;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
nbIter=4;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
    mu[i] = moments( contours[i], false ); 
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mu[i].m00>800)
(mu[i].m00>1000)
    cout <<"cut in : ("<< mu[i].m10/mu[i].m00 <<","<<mu[i].m01/mu[i].m00<< ")\n";
}

imwrite("test.jpg",m4);
return 0;
};

0- binarize image

1-Execute four times erode using a cross operator of size 5

2-With previous result Execute 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

4- inverse result and look for contour and filter results..

One cut is missing but I think you can start with this program :

int main (int argc,char **argv)
{
Mat     statComposante; /*< http://docs.opencv.org/trunk/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html */
Mat     centreGComposante; 

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 7, 7 ),cv::Point( 3, 3 ) );


Mat  m=imread("outputImage.jpg",CV_LOAD_IMAGE_GRAYSCALE);
m=imread("C:/Users/Laurent.PC-LAURENT-VISI/Downloads/simplifiedInput.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
imwrite("m1.png",m1);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
imwrite("m2.png",m2);
nbIter=7;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
nbIter=4;
imwrite("m3.png",m3);
nbIter=5;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
imwrite("m4.png",m4);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mu(contours.size() mo(contours.size() );
vector<Mat> hu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
    mu[i] {
    mo[i] = moments( contours[i], false ); 
    HuMoments((const Moments)mo[i], hu[i]);
}
for( int i = 0; i < contours.size(); i++ )
{
    for (int j=0;j<7;j++)
            cout << hu[i].at<double>(j, 0)<<"\t";
    cout << "\n";
}
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mu[i].m00>1000)
    cout <<"cut (mo[i].m00>1000 && mo[i].m00<5000)
        cout<<i <<" cut in : ("<< mu[i].m10/mu[i].m00 <<","<<mu[i].m01/mu[i].m00<< ")\n";
mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")"<<1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02))<<"\n";
}
int nb=0;
double pi = acos(-1.0);
for( int i = 0; i < contours.size(); i++ )
{ 
    if (fabs(hu[i].at<double>(0, 0) - 0.2) < 0.02)
    {
        cout<<nb<<"   "<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")" <<hu[i]<<"\n";
        double theta=1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02));
        double xg=mo[i].m10/mo[i].m00,yg=mo[i].m01/mo[i].m00;
        double xf,yf;
        double minorAxis2= sqrt((mo[i].mu20+mo[i].mu02-sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
        if (mo[i].mu20>mo[i].mu02)
            theta  += pi/2;
        if (theta<0)
            theta += pi;
        double majorAxis2= sqrt((mo[i].mu20+mo[i].mu02+sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
        xf = xg + sqrt(minorAxis2)*cos(theta+pi/2);
        yf = yg + sqrt(minorAxis2)*sin(theta+pi/2);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(0),5);
        xf = xg - sqrt(minorAxis2)*cos(theta+pi/2);
        yf = yg - sqrt(minorAxis2)*sin(theta+pi/2);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(0),5);
        xf = xg + sqrt(majorAxis2)*cos(theta);
        yf = yg + sqrt(majorAxis2)*sin(theta);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(128),5);
        xf = xg - sqrt(majorAxis2)*cos(theta);
        yf = yg - sqrt(majorAxis2)*sin(theta);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(128),5);
        nb++;
    }
}
imwrite("Fin.jpg",m);

imwrite("test.jpg",m4);
return 0;
};

with this result image description

0- binarize image

1-Execute four times erode using a cross operator of size 5

2-With previous result Execute 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

4- inverse result and look for contour and filter results..

One cut is missing but I think you can start with this program :

int main (int argc,char **argv)
{

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 7, 7 ),cv::Point( 3, 3 ) );


Mat  m=imread("C:/Users/Laurent.PC-LAURENT-VISI/Downloads/simplifiedInput.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
imwrite("m1.png",m1);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
imwrite("m2.png",m2);
nbIter=7;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
imwrite("m3.png",m3);
nbIter=5;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
imwrite("m4.png",m4);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mo(contours.size() );
vector<Mat> hu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
{
    mo[i] = moments( contours[i], false ); 
    HuMoments((const Moments)mo[i], hu[i]);
}
for( int i = 0; i < contours.size(); i++ )
{
    for (int j=0;j<7;j++)
            cout << hu[i].at<double>(j, 0)<<"\t";
    cout << "\n";
}
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mo[i].m00>1000 && mo[i].m00<5000)
        cout<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")"<<1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02))<<"\n";
}
int nb=0;
double pi = acos(-1.0);
for( int i = 0; i < contours.size(); i++ )
{ 
    if (fabs(hu[i].at<double>(0, 0) - 0.2) < 0.02)
    {
        cout<<nb<<"   "<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")" <<hu[i]<<"\n";
        double theta=1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02));
        double xg=mo[i].m10/mo[i].m00,yg=mo[i].m01/mo[i].m00;
        double xf,yf;
        double minorAxis2= sqrt((mo[i].mu20+mo[i].mu02-sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
        if (mo[i].mu20>mo[i].mu02)
            theta  += pi/2;
        if (theta<0)
            theta += pi;
        double majorAxis2= sqrt((mo[i].mu20+mo[i].mu02+sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
        xf = xg + sqrt(minorAxis2)*cos(theta+pi/2);
sqrt(minorAxis2)*cos(theta);
        yf = yg + sqrt(minorAxis2)*sin(theta+pi/2);
sqrt(minorAxis2)*sin(theta);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(0),5);
        xf = xg - sqrt(minorAxis2)*cos(theta+pi/2);
sqrt(minorAxis2)*cos(theta);
        yf = yg - sqrt(minorAxis2)*sin(theta+pi/2);
sqrt(minorAxis2)*sin(theta);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(0),5);
        xf = xg + sqrt(majorAxis2)*cos(theta);
sqrt(majorAxis2)*cos(theta+pi/2);
        yf = yg + sqrt(majorAxis2)*sin(theta);
sqrt(majorAxis2)*sin(theta+pi/2);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(128),5);
        xf = xg - sqrt(majorAxis2)*cos(theta);
sqrt(majorAxis2)*cos(theta+pi/2);
        yf = yg - sqrt(majorAxis2)*sin(theta);
sqrt(majorAxis2)*sin(theta+pi/2);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(128),5);
        nb++;
    }
}
imwrite("Fin.jpg",m);

imwrite("test.jpg",m4);
return 0;
};

with this result

image description

0- binarize image

1-Execute four times erode using a cross operator of size 5

2-With previous result Execute 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

4- inverse result and look for contour and filter results..

One cut is missing but I think you can start with this program :

int main (int argc,char **argv)
{

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 7, 7 ),cv::Point( 3, 3 ) );


Mat  m=imread("C:/Users/Laurent.PC-LAURENT-VISI/Downloads/simplifiedInput.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
imwrite("m1.png",m1);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
imwrite("m2.png",m2);
nbIter=7;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
imwrite("m3.png",m3);
nbIter=5;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
imwrite("m4.png",m4);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mo(contours.size() );
vector<Mat> hu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
{
    mo[i] = moments( contours[i], false ); 
    HuMoments((const Moments)mo[i], hu[i]);
}
for( int i = 0; i < contours.size(); i++ )
{
    for (int j=0;j<7;j++)
            cout << hu[i].at<double>(j, 0)<<"\t";
    cout << "\n";
}
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mo[i].m00>1000 && mo[i].m00<5000)
        cout<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")"<<1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02))<<"\n";
}
int nb=0;
double pi = acos(-1.0);
for( int i = 0; i < contours.size(); i++ )
{ 
 if (fabs(hu[i].at<double>(0, 0) - 0.2) < 0.02)
 {
     cout<<nb<<"   "<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")" <<hu[i]<<"\n";
     double theta=1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02));
     double xg=mo[i].m10/mo[i].m00,yg=mo[i].m01/mo[i].m00;
        double xf,yf;
    double xf1,yf1;
    double xf2,yf2;
    double minorAxis2= sqrt((mo[i].mu20+mo[i].mu02-sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
     if (mo[i].mu20>mo[i].mu02)
         theta  += pi/2;
     if (theta<0)
         theta += pi;
     double majorAxis2= sqrt((mo[i].mu20+mo[i].mu02+sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
        xf = xg + sqrt(minorAxis2)*cos(theta);
        yf = yg + sqrt(minorAxis2)*sin(theta);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(0),5);
        xf = xg - sqrt(minorAxis2)*cos(theta);
        yf = yg - sqrt(minorAxis2)*sin(theta);
        line(m, Point(xg, yg), Point(xf, yf), Scalar(0),5);
        xf xf1 = xg + sqrt(majorAxis2)*cos(theta+pi/2);
        yf yf1 = yg + sqrt(majorAxis2)*sin(theta+pi/2);
    xf2 = xg - sqrt(majorAxis2)*cos(theta+pi/2);
    yf2 = yg - sqrt(majorAxis2)*sin(theta+pi/2);
    line(m, Point(xg, yg), Point(xf, yf), Scalar(128),5);
        xf Point(xf1, yf1), Point(xf2, yf2), Scalar(255),2*sqrt(majorAxis2));
    xf1 = xg - + sqrt(minorAxis2)*cos(theta)-2* sqrt(majorAxis2)*cos(theta+pi/2);
        yf yf1 = yg - + sqrt(minorAxis2)*sin(theta)-2* sqrt(majorAxis2)*sin(theta+pi/2);
    xf2 = xg - sqrt(minorAxis2)*cos(theta)-2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf2 = yg - sqrt(minorAxis2)*sin(theta)-2* sqrt(majorAxis2)*sin(theta+pi/2);
    line(m, Point(xg, yg), Point(xf, yf), Scalar(128),5);
    Point(xf1, yf1), Point(xf2, yf2), Scalar(0),5);
    xf1 = xg + sqrt(minorAxis2)*cos(theta)+2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf1 = yg + sqrt(minorAxis2)*sin(theta)+2* sqrt(majorAxis2)*sin(theta+pi/2);
    xf2 = xg - sqrt(minorAxis2)*cos(theta)+2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf2 = yg - sqrt(minorAxis2)*sin(theta)+2* sqrt(majorAxis2)*sin(theta+pi/2);
    line(m, Point(xf1, yf1), Point(xf2, yf2), Scalar(0),5);
    nb++;
    }
}
imwrite("Fin.jpg",m);

imwrite("test.jpg",m4);
return 0;
};

with this result

0- binarize image

1-Execute four times erode using a cross operator of size 5

2-With previous result Execute 8 times dilate using same cross

3- With previous result execute 4 times erode using same cross

4- inverse result and look for contour and filter results..

One cut is missing but I think you can start with this program :

int main (int argc,char **argv)
{

Mat mTest,mThresh,mConnected;
int nbIter;
int borderType=cv::BORDER_REPLICATE;
Scalar borderValue=cv::morphologyDefaultBorderValue();
Point ancrage=cv::Point(-1,-1);
cv::Mat element = cv::getStructuringElement( MORPH_CROSS,cv::Size( 7, 7 ),cv::Point( 3, 3 ) );


Mat  m=imread("C:/Users/Laurent.PC-LAURENT-VISI/Downloads/simplifiedInput.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat m1,m2,m3,m4,m5;
nbIter=4;
threshold(m,m1,50,255,THRESH_BINARY);
imwrite("m1.png",m1);
erode(m1, m2, element,ancrage,nbIter,borderType,borderValue);
imwrite("m2.png",m2);
nbIter=7;
dilate(m2, m3, element,ancrage,nbIter,borderType,borderValue);
imwrite("m3.png",m3);
nbIter=5;
erode(m3, m4, element,ancrage,nbIter,borderType,borderValue);
imwrite("m4.png",m4);
bitwise_not(m4,m5);
Mat mClone = m5.clone();
vector<vector<cv::Point> > contours;

findContours( mClone,  contours,cv::RETR_EXTERNAL,cv::CHAIN_APPROX_NONE);   
vector<Moments> mo(contours.size() );
vector<Mat> hu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
{
    mo[i] = moments( contours[i], false ); 
    HuMoments((const Moments)mo[i], hu[i]);
}
for( int i = 0; i < contours.size(); i++ )
{
    for (int j=0;j<7;j++)
            cout << hu[i].at<double>(j, 0)<<"\t";
    cout << "\n";
}
for( int i = 0; i < contours.size(); i++ )
{ 
    if (mo[i].m00>1000 && mo[i].m00<5000)
        cout<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")"<<1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02))<<"\n";
}
int nb=0;
double pi = acos(-1.0);
for( int i = 0; i < contours.size(); i++ )
{ 
//    if (fabs(hu[i].at<double>(0, 0) - 0.2) < 0.02)
if (mo[i].m00>1000 && mo[i].m00<5000)
    {
    cout<<nb<<"   "<<i <<" cut in : ("<< mo[i].m10/mo[i].m00 <<","<<mo[i].m01/mo[i].m00<< ")" <<hu[i]<<"\n";
    double theta=1/2.0*atan(2*mo[i].mu11/(mo[i].mu20-mo[i].mu02));
    double xg=mo[i].m10/mo[i].m00,yg=mo[i].m01/mo[i].m00;
    double xf1,yf1;
    double xf2,yf2;
    double minorAxis2= sqrt((mo[i].mu20+mo[i].mu02-sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
    if (mo[i].mu20>mo[i].mu02)
        theta  += pi/2;
    if (theta<0)
        theta += pi;
    double majorAxis2= sqrt((mo[i].mu20+mo[i].mu02+sqrt(pow(mo[i].mu20-mo[i].mu02,2.0)+4*mo[i].mu11*mo[i].mu11))/2);
    xf1 = xg + sqrt(majorAxis2)*cos(theta+pi/2);
    yf1 = yg + sqrt(majorAxis2)*sin(theta+pi/2);
    xf2 = xg - sqrt(majorAxis2)*cos(theta+pi/2);
    yf2 = yg - sqrt(majorAxis2)*sin(theta+pi/2);
    line(m, Point(xf1, yf1), Point(xf2, yf2), Scalar(255),2*sqrt(majorAxis2));
    xf1 = xg + sqrt(minorAxis2)*cos(theta)-2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf1 = yg + sqrt(minorAxis2)*sin(theta)-2* sqrt(majorAxis2)*sin(theta+pi/2);
    xf2 = xg - sqrt(minorAxis2)*cos(theta)-2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf2 = yg - sqrt(minorAxis2)*sin(theta)-2* sqrt(majorAxis2)*sin(theta+pi/2);
    line(m, Point(xf1, yf1), Point(xf2, yf2), Scalar(0),5);
    xf1 = xg + sqrt(minorAxis2)*cos(theta)+2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf1 = yg + sqrt(minorAxis2)*sin(theta)+2* sqrt(majorAxis2)*sin(theta+pi/2);
    xf2 = xg - sqrt(minorAxis2)*cos(theta)+2* sqrt(majorAxis2)*cos(theta+pi/2);
    yf2 = yg - sqrt(minorAxis2)*sin(theta)+2* sqrt(majorAxis2)*sin(theta+pi/2);
    line(m, Point(xf1, yf1), Point(xf2, yf2), Scalar(0),5);
    nb++;
}
imwrite("Fin.jpg",m);

imwrite("test.jpg",m4);
return 0;
};

with this result