00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "precomp.hpp"
00043 #include "gcgraph.hpp"
00044 #include "grabcut_3d/grabcut_3d.h"
00045 #include <limits>
00046 #include <opencv2/highgui/highgui.hpp>
00047 #include <opencv2/core/core.hpp>
00048
00049 using namespace cv;
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 class GMM
00061 {
00062 public:
00063 static const int componentsCount = 5;
00064
00065 GMM( Mat& _model );
00066 float operator()( Vec3f color ) const;
00067 float operator()( int ci, Vec3f color ) const;
00068 int whichComponent( Vec3f color ) const;
00069
00070 void initLearning();
00071 void addSample( int ci, Vec3f color );
00072 void endLearning();
00073 private:
00074 void calcInverseCovAndDeterm( int ci );
00075 Mat model;
00076 float* coefs;
00077 float* mean;
00078 float* cov;
00079
00080 float inverseCovs[componentsCount][3][3];
00081 float covDeterms[componentsCount];
00082
00083 float sums[componentsCount][3];
00084 float prods[componentsCount][3][3];
00085 int sampleCounts[componentsCount];
00086 int totalSampleCount;
00087 };
00088
00089 GMM::GMM( Mat& _model )
00090 {
00091 const int modelSize = 3 + 9 + 1;
00092 if( _model.empty() )
00093 {
00094 _model.create( 1, modelSize*componentsCount, CV_32FC1 );
00095 _model.setTo(Scalar(0));
00096 }
00097 else if( (_model.type() != CV_32FC1) || (_model.rows != 1) || (_model.cols != modelSize*componentsCount) )
00098 CV_Error( CV_StsBadArg, "_model must have CV_32FC1 type, rows == 1 and cols == 13*componentsCount" );
00099
00100 model = _model;
00101
00102 coefs = model.ptr<float>(0);
00103 mean = coefs + componentsCount;
00104 cov = mean + 3*componentsCount;
00105
00106 for( int ci = 0; ci < componentsCount; ci++ )
00107 if( coefs[ci] > 0 )
00108 calcInverseCovAndDeterm( ci );
00109 }
00110
00111 float GMM::operator()( Vec3f color ) const
00112 {
00113 float res = 0;
00114 for( int ci = 0; ci < componentsCount; ci++ )
00115 res += coefs[ci] * (*this)(ci, color );
00116 return res;
00117 }
00118
00119 float GMM::operator()( int ci, Vec3f color ) const
00120 {
00121 float res = 0;
00122 if( coefs[ci] > 0 )
00123 {
00124 if( covDeterms[ci] > std::numeric_limits<float>::epsilon() )
00125 {
00126 Vec3f diff = color;
00127 float* m = mean + 3*ci;
00128 diff[0] -= m[0]; diff[1] -= m[1]; diff[2] -= m[2];
00129 float mult = diff[0]*(diff[0]*inverseCovs[ci][0][0] + diff[1]*inverseCovs[ci][1][0] + diff[2]*inverseCovs[ci][2][0])
00130 + diff[1]*(diff[0]*inverseCovs[ci][0][1] + diff[1]*inverseCovs[ci][1][1] + diff[2]*inverseCovs[ci][2][1])
00131 + diff[2]*(diff[0]*inverseCovs[ci][0][2] + diff[1]*inverseCovs[ci][1][2] + diff[2]*inverseCovs[ci][2][2]);
00132 res = 1.0f/sqrt(covDeterms[ci]) * exp(-0.5f*mult);
00133 }
00134 }
00135 return res;
00136 }
00137
00138 int GMM::whichComponent( Vec3f color ) const
00139 {
00140 int k = 0;
00141 float max = 0;
00142
00143 for( int ci = 0; ci < componentsCount; ci++ )
00144 {
00145 float p = (*this)( ci, color );
00146 if( p > max )
00147 {
00148 k = ci;
00149 max = p;
00150 }
00151 }
00152 return k;
00153 }
00154
00155 void GMM::initLearning()
00156 {
00157 for( int ci = 0; ci < componentsCount; ci++)
00158 {
00159 sums[ci][0] = sums[ci][1] = sums[ci][2] = 0;
00160 prods[ci][0][0] = prods[ci][0][1] = prods[ci][0][2] = 0;
00161 prods[ci][1][0] = prods[ci][1][1] = prods[ci][1][2] = 0;
00162 prods[ci][2][0] = prods[ci][2][1] = prods[ci][2][2] = 0;
00163 sampleCounts[ci] = 0;
00164 }
00165 totalSampleCount = 0;
00166 }
00167
00168 void GMM::addSample( int ci, Vec3f color )
00169 {
00170 sums[ci][0] += color[0]; sums[ci][1] += color[1]; sums[ci][2] += color[2];
00171 prods[ci][0][0] += color[0]*color[0]; prods[ci][0][1] += color[0]*color[1]; prods[ci][0][2] += color[0]*color[2];
00172 prods[ci][1][0] += color[1]*color[0]; prods[ci][1][1] += color[1]*color[1]; prods[ci][1][2] += color[1]*color[2];
00173 prods[ci][2][0] += color[2]*color[0]; prods[ci][2][1] += color[2]*color[1]; prods[ci][2][2] += color[2]*color[2];
00174 sampleCounts[ci]++;
00175 totalSampleCount++;
00176 }
00177
00178 void GMM::endLearning()
00179 {
00180 for( int ci = 0; ci < componentsCount; ci++ )
00181 {
00182 int n = sampleCounts[ci];
00183 if( n == 0 )
00184 coefs[ci] = 0;
00185 else
00186 {
00187 coefs[ci] = (float)n/totalSampleCount;
00188
00189 float* m = mean + 3*ci;
00190 m[0] = sums[ci][0]/n; m[1] = sums[ci][1]/n; m[2] = sums[ci][2]/n;
00191
00192 float* c = cov + 9*ci;
00193 c[0] = prods[ci][0][0]/n - m[0]*m[0]; c[1] = prods[ci][0][1]/n - m[0]*m[1]; c[2] = prods[ci][0][2]/n - m[0]*m[2];
00194 c[3] = prods[ci][1][0]/n - m[1]*m[0]; c[4] = prods[ci][1][1]/n - m[1]*m[1]; c[5] = prods[ci][1][2]/n - m[1]*m[2];
00195 c[6] = prods[ci][2][0]/n - m[2]*m[0]; c[7] = prods[ci][2][1]/n - m[2]*m[1]; c[8] = prods[ci][2][2]/n - m[2]*m[2];
00196
00197 calcInverseCovAndDeterm(ci);
00198 }
00199 }
00200 }
00201
00202 void GMM::calcInverseCovAndDeterm( int ci )
00203 {
00204 if( coefs[ci] > 0 )
00205 {
00206 float *c = cov + 9*ci;
00207 float dtrm =
00208 covDeterms[ci] = c[0]*(c[4]*c[8]-c[5]*c[7]) - c[1]*(c[3]*c[8]-c[5]*c[6]) + c[2]*(c[3]*c[7]-c[4]*c[6]);
00209
00210 if( dtrm > std::numeric_limits<float>::epsilon() )
00211 {
00212 inverseCovs[ci][0][0] = (c[4]*c[8] - c[5]*c[7]) / dtrm;
00213 inverseCovs[ci][1][0] = -(c[3]*c[8] - c[5]*c[6]) / dtrm;
00214 inverseCovs[ci][2][0] = (c[3]*c[7] - c[4]*c[6]) / dtrm;
00215 inverseCovs[ci][0][1] = -(c[1]*c[8] - c[2]*c[7]) / dtrm;
00216 inverseCovs[ci][1][1] = (c[0]*c[8] - c[2]*c[6]) / dtrm;
00217 inverseCovs[ci][2][1] = -(c[0]*c[7] - c[1]*c[6]) / dtrm;
00218 inverseCovs[ci][0][2] = (c[1]*c[5] - c[2]*c[4]) / dtrm;
00219 inverseCovs[ci][1][2] = -(c[0]*c[5] - c[2]*c[3]) / dtrm;
00220 inverseCovs[ci][2][2] = (c[0]*c[4] - c[1]*c[3]) / dtrm;
00221 }
00222 }
00223 }
00224
00225 #define ALPHA 1.0
00226 #define BETA 1.0
00227
00228 Vec3f getDataTerm( const Mat& img, const Mat& depth_img, int row, int col )
00229 {
00230 const Vec3b& color = img.at<Vec3b>(row,col);
00231
00232
00233 return Vec3f(color[0],color[1],color[2]);
00234 }
00235
00236 Vec3f getDataTerm( const Mat& img, const Mat& depth_img, Point& p)
00237 {
00238 const Vec3b& color = img.at<Vec3b>(p);
00239
00240
00241 return Vec3f(color[0],color[1],color[2]);
00242 }
00243
00244
00245
00246
00247
00248 float calcBeta( const Mat& img, const Mat& depth_img )
00249 {
00250 float beta = 0;
00251 for( int y = 0; y < img.rows; y++ )
00252 {
00253 for( int x = 0; x < img.cols; x++ )
00254 {
00255 Vec3f data = getDataTerm(img,depth_img,y,x);
00256
00257 if( x>0 )
00258 {
00259 Vec3f diff = data - getDataTerm(img,depth_img,y,x-1);
00260 diff[2]*=ALPHA;
00261 beta += diff.dot(diff);
00262 }
00263 if( y>0 && x>0 )
00264 {
00265 Vec3f diff = data - getDataTerm(img,depth_img,y-1,x-1);
00266 diff[2]*=ALPHA;
00267 beta += diff.dot(diff);
00268 }
00269 if( y>0 )
00270 {
00271 Vec3f diff = data - getDataTerm(img,depth_img,y-1,x);
00272 diff[2]*=ALPHA;
00273 beta += diff.dot(diff);
00274 }
00275 if( y>0 && x<img.cols-1)
00276 {
00277 Vec3f diff = data - getDataTerm(img,depth_img,y-1,x+1);
00278 diff[2]*=ALPHA;
00279 beta += diff.dot(diff);
00280 }
00281 }
00282 }
00283 beta = 1.f / (2 * beta/(4*img.cols*img.rows - 3*img.cols - 3*img.rows + 2) );
00284 return beta;
00285 }
00286
00287
00288
00289
00290
00291 void calcNWeights( const Mat& img, const Mat& depth_img, Mat& leftW, Mat& upleftW, Mat& upW, Mat& uprightW, float beta, float gamma )
00292 {
00293 const float gammaDivSqrt2 = gamma / std::sqrt(2.0f);
00294 leftW.create( img.rows, img.cols, CV_32FC1 );
00295 upleftW.create( img.rows, img.cols, CV_32FC1 );
00296 upW.create( img.rows, img.cols, CV_32FC1 );
00297 uprightW.create( img.rows, img.cols, CV_32FC1 );
00298 for( int y = 0; y < img.rows; y++ )
00299 {
00300 for( int x = 0; x < img.cols; x++ )
00301 {
00302 Vec3f data = getDataTerm(img,depth_img,y,x);
00303
00304 if( x-1>=0 )
00305 {
00306 Vec3f diff = data - getDataTerm(img,depth_img,y,x-1);
00307 diff[2]*=ALPHA;
00308 leftW.at<float>(y,x) = gamma * exp(-beta*diff.dot(diff));
00309 }
00310 else
00311 leftW.at<float>(y,x) = 0;
00312 if( x-1>=0 && y-1>=0 )
00313 {
00314 Vec3f diff = data - getDataTerm(img,depth_img,y-1,x-1);
00315 diff[2]*=ALPHA;
00316 upleftW.at<float>(y,x) = gammaDivSqrt2 * exp(-beta*diff.dot(diff));
00317 }
00318 else
00319 upleftW.at<float>(y,x) = 0;
00320 if( y-1>=0 )
00321 {
00322 Vec3f diff = data - getDataTerm(img,depth_img,y-1,x);
00323 diff[2]*=ALPHA;
00324 upW.at<float>(y,x) = gamma * exp(-beta*diff.dot(diff));
00325 }
00326 else
00327 upW.at<float>(y,x) = 0;
00328 if( x+1<img.cols-1 && y-1>=0 )
00329 {
00330 Vec3f diff = data - getDataTerm(img,depth_img,y-1,x+1);
00331 diff[2]*=ALPHA;
00332 uprightW.at<float>(y,x) = gammaDivSqrt2 * exp(-beta*diff.dot(diff));
00333 }
00334 else
00335 uprightW.at<float>(y,x) = 0;
00336 }
00337 }
00338 }
00339
00340
00341
00342
00343 void checkMask( const Mat& img, const Mat& mask )
00344 {
00345 if( mask.empty() )
00346 CV_Error( CV_StsBadArg, "mask is empty" );
00347 if( mask.type() != CV_8UC1 )
00348 CV_Error( CV_StsBadArg, "mask must have CV_8UC1 type" );
00349 if( mask.cols != img.cols || mask.rows != img.rows )
00350 CV_Error( CV_StsBadArg, "mask must have as many rows and cols as img" );
00351 for( int y = 0; y < mask.rows; y++ )
00352 {
00353 for( int x = 0; x < mask.cols; x++ )
00354 {
00355 uchar val = mask.at<uchar>(y,x);
00356 if( val!=GC_BGD && val!=GC_FGD && val!=GC_PR_BGD && val!=GC_PR_FGD )
00357 CV_Error( CV_StsBadArg, "mask element value must be equel"
00358 "GC_BGD or GC_FGD or GC_PR_BGD or GC_PR_FGD" );
00359 }
00360 }
00361 }
00362
00363
00364
00365
00366 void initMaskWithRect( Mat& mask, Size imgSize, Rect rect )
00367 {
00368 mask.create( imgSize, CV_8UC1 );
00369 mask.setTo( GC_BGD );
00370
00371 rect.x = max(0, rect.x);
00372 rect.y = max(0, rect.y);
00373 rect.width = min(rect.width, imgSize.width-rect.x);
00374 rect.height = min(rect.height, imgSize.height-rect.y);
00375
00376 (mask(rect)).setTo( Scalar(GC_PR_FGD) );
00377 }
00378
00379
00380
00381
00382 void initGMMs( const Mat& img, const Mat& depth_img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM )
00383 {
00384 const int kMeansItCount = 10;
00385 const int kMeansType = KMEANS_PP_CENTERS;
00386
00387 Mat bgdLabels, fgdLabels;
00388 vector<Vec3f> bgdSamples, fgdSamples;
00389 Point p;
00390 for( p.y = 0; p.y < img.rows; p.y++ )
00391 {
00392 for( p.x = 0; p.x < img.cols; p.x++ )
00393 {
00394 if( mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD ) {
00395 bgdSamples.push_back(getDataTerm(img,depth_img,p));
00396 }
00397 else
00398 {
00399 fgdSamples.push_back(getDataTerm(img,depth_img,p));
00400 }
00401 }
00402 }
00403 CV_Assert( !bgdSamples.empty() && !fgdSamples.empty() );
00404 Mat _bgdSamples( (int)bgdSamples.size(), 3, CV_32FC1, &bgdSamples[0][0] );
00405 kmeans( _bgdSamples, GMM::componentsCount, bgdLabels,
00406 TermCriteria( CV_TERMCRIT_ITER, kMeansItCount, 0.0), 0, kMeansType);
00407 Mat _fgdSamples( (int)fgdSamples.size(), 3, CV_32FC1, &fgdSamples[0][0] );
00408 kmeans( _fgdSamples, GMM::componentsCount, fgdLabels,
00409 TermCriteria( CV_TERMCRIT_ITER, kMeansItCount, 0.0), 0, kMeansType);
00410
00411 bgdGMM.initLearning();
00412 for( int i = 0; i < (int)bgdSamples.size(); i++ )
00413 bgdGMM.addSample( bgdLabels.at<int>(i,0), bgdSamples[i] );
00414 bgdGMM.endLearning();
00415
00416 fgdGMM.initLearning();
00417 for( int i = 0; i < (int)fgdSamples.size(); i++ )
00418 fgdGMM.addSample( fgdLabels.at<int>(i,0), fgdSamples[i] );
00419 fgdGMM.endLearning();
00420 }
00421
00422
00423
00424
00425 void assignGMMsComponents( const Mat& img, const Mat& depth_img, const Mat& mask, const GMM& bgdGMM, const GMM& fgdGMM, Mat& compIdxs )
00426 {
00427 Point p;
00428 for( p.y = 0; p.y < img.rows; p.y++ )
00429 {
00430 for( p.x = 0; p.x < img.cols; p.x++ )
00431 {
00432 Vec3f data = getDataTerm(img,depth_img,p);
00433
00434 compIdxs.at<int>(p) = mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD ?
00435 bgdGMM.whichComponent(data) : fgdGMM.whichComponent(data);
00436 }
00437 }
00438 }
00439
00440
00441
00442
00443 void learnGMMs( const Mat& img, const Mat& depth_img, const Mat& mask, const Mat& compIdxs, GMM& bgdGMM, GMM& fgdGMM )
00444 {
00445 bgdGMM.initLearning();
00446 fgdGMM.initLearning();
00447 Point p;
00448 for( int ci = 0; ci < GMM::componentsCount; ci++ )
00449 {
00450 for( p.y = 0; p.y < img.rows; p.y++ )
00451 {
00452 for( p.x = 0; p.x < img.cols; p.x++ )
00453 {
00454 if( compIdxs.at<int>(p) == ci )
00455 {
00456 if( mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD ) {
00457 Vec3f data = getDataTerm(img,depth_img,p);
00458 bgdGMM.addSample( ci, data );
00459 }
00460 else {
00461 Vec3f data = getDataTerm(img,depth_img,p);
00462 fgdGMM.addSample( ci, data );
00463 }
00464
00465 }
00466 }
00467 }
00468 }
00469 bgdGMM.endLearning();
00470 fgdGMM.endLearning();
00471 }
00472
00473
00474
00475
00476 void constructGCGraph( const Mat& img, const Mat& depth_img, const Mat& mask, const GMM& bgdGMM, const GMM& fgdGMM, float lambda,
00477 const Mat& leftW, const Mat& upleftW, const Mat& upW, const Mat& uprightW,
00478 GCGraph<float>& graph )
00479 {
00480 int vtxCount = img.cols*img.rows,
00481 edgeCount = 2*(4*img.cols*img.rows - 3*(img.cols + img.rows) + 2);
00482 graph.create(vtxCount, edgeCount);
00483 Point p;
00484 for( p.y = 0; p.y < img.rows; p.y++ )
00485 {
00486 for( p.x = 0; p.x < img.cols; p.x++)
00487 {
00488
00489 int vtxIdx = graph.addVtx();
00490 Vec3f data = getDataTerm(img,depth_img,p);
00491
00492
00493 float fromSource, toSink;
00494 if( mask.at<uchar>(p) == GC_PR_BGD || mask.at<uchar>(p) == GC_PR_FGD )
00495 {
00496 fromSource = -log( bgdGMM(data) );
00497 toSink = -log( fgdGMM(data) );
00498 }
00499 else if( mask.at<uchar>(p) == GC_BGD )
00500 {
00501 fromSource = 0;
00502 toSink = lambda;
00503 }
00504 else
00505 {
00506 fromSource = lambda;
00507 toSink = 0;
00508 }
00509 graph.addTermWeights( vtxIdx, fromSource, toSink );
00510
00511
00512 if( p.x>0 )
00513 {
00514 float w = leftW.at<float>(p);
00515 graph.addEdges( vtxIdx, vtxIdx-1, w, w );
00516 }
00517 if( p.x>0 && p.y>0 )
00518 {
00519 float w = upleftW.at<float>(p);
00520 graph.addEdges( vtxIdx, vtxIdx-img.cols-1, w, w );
00521 }
00522 if( p.y>0 )
00523 {
00524 float w = upW.at<float>(p);
00525 graph.addEdges( vtxIdx, vtxIdx-img.cols, w, w );
00526 }
00527 if( p.x<img.cols-1 && p.y>0 )
00528 {
00529 float w = uprightW.at<float>(p);
00530 graph.addEdges( vtxIdx, vtxIdx-img.cols+1, w, w );
00531 }
00532 }
00533 }
00534 }
00535
00536
00537
00538
00539 void estimateSegmentation( GCGraph<float>& graph, Mat& mask )
00540 {
00541 graph.maxFlow();
00542 Point p;
00543 for( p.y = 0; p.y < mask.rows; p.y++ )
00544 {
00545 for( p.x = 0; p.x < mask.cols; p.x++ )
00546 {
00547 if( mask.at<uchar>(p) == GC_PR_BGD || mask.at<uchar>(p) == GC_PR_FGD )
00548 {
00549 if( graph.inSourceSegment( p.y*mask.cols+p.x ) )
00550 mask.at<uchar>(p) = GC_PR_FGD;
00551 else
00552 mask.at<uchar>(p) = GC_PR_BGD;
00553 }
00554 }
00555 }
00556 }
00557
00558 void cv::grabCut3D( const Mat& img, const Mat& depth_img, Mat& mask, Rect rect,
00559 Mat& bgdModel, Mat& fgdModel,
00560 int iterCount, int mode )
00561 {
00562 if( img.empty() )
00563 CV_Error( CV_StsBadArg, "image is empty" );
00564 if( img.type() != CV_8UC3 )
00565 CV_Error( CV_StsBadArg, "image mush have CV_8UC3 type" );
00566
00567 if( depth_img.empty() )
00568 CV_Error( CV_StsBadArg, "depth_img is empty" );
00569 if( depth_img.type() != CV_8UC3 )
00570 CV_Error( CV_StsBadArg, "depth_img mush have CV_8UC3 type" );
00571
00572 GMM bgdGMM( bgdModel ), fgdGMM( fgdModel );
00573 Mat compIdxs( img.size(), CV_32SC1 );
00574
00575 if( mode == GC_INIT_WITH_RECT || mode == GC_INIT_WITH_MASK )
00576 {
00577 if( mode == GC_INIT_WITH_RECT )
00578 initMaskWithRect( mask, img.size(), rect );
00579 else
00580 checkMask( img, mask );
00581 initGMMs( img, depth_img, mask, bgdGMM, fgdGMM );
00582 }
00583
00584 if( iterCount <= 0)
00585 return;
00586
00587 if( mode == GC_EVAL )
00588 checkMask( img, mask );
00589
00590 const float gamma = 50;
00591 const float lambda = 9*gamma;
00592 const float beta = calcBeta( img, depth_img );
00593
00594 Mat leftW, upleftW, upW, uprightW;
00595 calcNWeights( img, depth_img, leftW, upleftW, upW, uprightW, beta, gamma );
00596
00597 for( int i = 0; i < iterCount; i++ )
00598 {
00599 GCGraph<float> graph;
00600 assignGMMsComponents( img, depth_img, mask, bgdGMM, fgdGMM, compIdxs );
00601 learnGMMs( img, depth_img, mask, compIdxs, bgdGMM, fgdGMM );
00602 constructGCGraph(img, depth_img, mask, bgdGMM, fgdGMM, lambda, leftW, upleftW, upW, uprightW, graph );
00603 estimateSegmentation( graph, mask );
00604 }
00605 }
