00001 #ifndef Impala_Core_Array_ColorSegmentation_h
00002 #define Impala_Core_Array_ColorSegmentation_h
00003
00004 #include <list>
00005 #include "Core/Matrix/MatrixTem.h"
00006 #include "Core/Array/ArrayListExport.h"
00007 #include "Core/Array/ArrayListDelete.h"
00008 #include "Core/Array/ColorSegmentationInvariant.h"
00009 #include "Core/Array/ColorSegmentationAlgorithm.h"
00010 #include "Core/Array/ColorGaussResponses.h"
00011 #include "Core/Array/ColorGaborResponses.h"
00012 #include "Basis/Timer.h"
00013 #include "Core/Matrix/MatKMeans.h"
00014 #include "Core/Array/WriteRaw.h"
00015 #include "Core/Matrix/MatTranspose.h"
00016 #include "Core/Matrix/MatKLM.h"
00017 #include "Core/Matrix/MatFeatureNormalization.h"
00018 #include "Core/Array/GetRgbPixels.h"
00019
00020 namespace Impala
00021 {
00022 namespace Core
00023 {
00024 namespace Array
00025 {
00026
00027
00028 typedef Impala::Core::Vector::VectorTem<double> VectorDouble;
00029 typedef Impala::Core::Matrix::MatrixTem<double> MatrixDouble;
00030 using namespace Impala::Core::Matrix;
00031
00032
00033 class pixgroup{
00034 public:
00035
00036 int size;
00037 int nrPixInGroup;
00038 int label;
00039 int dim;
00040 MatrixDouble cov;
00041 VectorDouble m;
00042 VectorDouble idx;
00043 };
00044
00046
00047 void labelimg(unsigned char* dst, const VectorDouble &imap, int width, int height) ;
00048
00049 void updateMean(MatrixDouble &m, pixgroup &rg, long inew) ;
00050
00051 void extractRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap,
00052 const MatrixDouble &m) ;
00053
00054 double regionDistance(pixgroup &u, pixgroup &v) ;
00055
00056 void rgmerge(pixgroup &toMerge, pixgroup ®ion, const MatrixDouble &m) ;
00057
00058 void mergeRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap,
00059 const MatrixDouble &m, int width, int height, double threshold) ;
00060
00061 void relabel(VectorDouble &clusterMap, int width, int height, pixgroup *rg,
00062 int &nrClusters, int numpixrm, long sizethreshold) ;
00063
00064 void refinery(VectorDouble &clusterMap, MatrixDouble &m, pixgroup *rg, int nrClusters,
00065 int width, int height) ;
00066
00068
00069 inline void
00070 ColorSegmentation(Array2dVec3UInt8*& dst, Array2dVec3Real64* im,
00071 ColorSegmentationAlgorithm segAlg,
00072 ColorSegmentationInvariant invariantType,
00073 double minRegionFraction, double threshold,
00074 bool useGauss = false,
00075 std::vector<Array2dVec3UInt8*>* dispList = 0)
00076 {
00077 bool vv = false;
00078 if (vv) std::cout << "start ColorSegmentation" << std::endl;
00079 if (dst == 0)
00080 dst = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
00081 Timer timer(1);
00082 timer.Start();
00083 std::vector<Array2dScalarReal64*> respList;
00084 if (useGauss)
00085 respList = ColorGaussResponses(im);
00086 else
00087 respList = ColorGaborResponses(im, segAlg, invariantType);
00088
00089
00090 if (vv) std::cout << " making response matrix" << std::endl;
00091
00092 typedef Array2dScalarReal64 Mat;
00093 typedef VecScalarReal64 Vec;
00094 int nrFeat = respList.size();
00095 int nrPixels = respList[0]->CW() * respList[0]->CH();
00096 int nrPcaDim = 4;
00097
00098 double* data = ArrayListExport(respList);
00099 ArrayListDelete(&respList);
00100
00101 Mat* mTmp = MatCreate<Mat>(nrFeat, nrPixels, data);
00102
00103 Mat* m = MatTranspose(mTmp);
00104 delete mTmp;
00105
00106 if (vv) std::cout << "time: " << timer.SplitTime() << " PCA " << std::endl;
00107
00108 if (MatNrCol(m) > nrPcaDim)
00109 {
00110 mTmp = MatKLM(m, nrPcaDim);
00111 delete m;
00112 m = mTmp;
00113 }
00114
00115 if (vv) std::cout << "time: " << timer.SplitTime() << " feature_normalization " << std::endl;
00116
00117 MatFeatureNormalization(m);
00118
00119 if (vv) std::cout << "time: " << timer.SplitTime() << " kmeans " << std::endl;
00120
00121
00122 int nrClusters = nrPixels/2500;
00123 Mat* clusters = MatCreate<Mat>(nrClusters, MatNrCol(m));
00124 Vec* clusterMap = VecCreate<Vec>(nrPixels);
00125 VecScalarInt32* vectorsInC = VecCreate<VecScalarInt32>(nrClusters);
00126 int clustersFound;
00127 MatKMeans(m, clusters, clusterMap, 25, false, clustersFound, vectorsInC);
00128 delete clusters;
00129 if (vv)
00130 for (int i=0 ; i<nrClusters ; i++)
00131 {
00132 int nr = *VecE(vectorsInC, i);
00133 std::cout << "cluster " << i << " has " << nr
00134 << " vectors" << std::endl;
00135 }
00136 delete vectorsInC;
00137
00138
00139 MatrixDouble mMat(MatNrRow(m), MatNrCol(m), MatE(m));
00140
00141 VectorDouble clusterMapMat(VecNrElem(clusterMap), VecE(clusterMap), false);
00142 delete clusterMap;
00143
00144 if (vv) std::cout << "time: " << timer.SplitTime() << " Doing region stuff " << std::endl;
00145 if (dispList != 0)
00146 {
00147 Array2dVec3UInt8* dIm = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
00148 labelimg(dIm->PB(), clusterMapMat, im->CW(), im->CH());
00149 dispList->push_back(dIm);
00150 }
00151
00152
00153 pixgroup *rg = new pixgroup[nrClusters];
00154 extractRegions(rg, nrClusters, clusterMapMat, mMat);
00155
00156
00157 mergeRegions(rg, nrClusters, clusterMapMat, mMat, im->CW(), im->CH(),
00158 threshold) ;
00159 if (dispList != 0)
00160 {
00161 Array2dVec3UInt8* dIm = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
00162 labelimg(dIm->PB(), clusterMapMat, im->CW(), im->CH());
00163 dispList->push_back(dIm);
00164 }
00165
00166
00167
00168 relabel(clusterMapMat, im->CW(), im->CH(), rg, nrClusters, 3,
00169 nrPixels * minRegionFraction);
00170 if (dispList != 0)
00171 {
00172 Array2dVec3UInt8* dIm = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
00173 labelimg(dIm->PB(), clusterMapMat, im->CW(), im->CH());
00174 dispList->push_back(dIm);
00175 }
00176
00177
00178 refinery(clusterMapMat, mMat, rg, nrClusters, im->CW(), im->CH());
00179
00180
00181 labelimg(dst->CPB(), clusterMapMat, im->CW(), im->CH());
00182
00183 delete m;
00184 delete[] rg;
00185 if (vv) std::cout << "time: " << timer.SplitTime() << " Done " << std::endl;
00186 }
00187
00188
00190
00191 void labelimg(unsigned char* dst, const VectorDouble &imap, int width, int height)
00192 {
00193
00194 unsigned char pl[10][2] = { {255,0},{191,63},{127,0},{255,63},{191,0},
00195 {255,127},{63,0}, {191,127},{255,191},{127,63}
00196 };
00197 unsigned char red[60],green[60],blue[60];
00198
00199
00200 int id = 0;
00201 long i;
00202 for (i=0; i<10; i++)
00203 {
00204 red[id] = pl[i][0]; green[id] = pl[i][1]; blue[id] = pl[i][1]; id++;
00205 red[id] = pl[i][0]; green[id] = pl[i][0]; blue[id] = pl[i][1]; id++;
00206 red[id] = pl[i][0]; green[id] = pl[i][1]; blue[id] = pl[i][0]; id++;
00207 red[id] = pl[i][1]; green[id] = pl[i][0]; blue[id] = pl[i][1]; id++;
00208 red[id] = pl[i][1]; green[id] = pl[i][0]; blue[id] = pl[i][0]; id++;
00209 red[id] = pl[i][1]; green[id] = pl[i][1]; blue[id] = pl[i][0]; id++;
00210 }
00211
00212 long imsize = width*height;
00213 id = 0;
00214
00215 for (i=0; i<imsize; i++)
00216 {
00217 if (imap[i] < 0)
00218 {
00219 dst[id++] = 0;
00220 dst[id++] = 0;
00221 dst[id++] = 0;
00222 }
00223 else
00224 {
00225 int cid = ((int)imap[i]) % 60;
00226 dst[id++] = red[cid];
00227 dst[id++] = green[cid];
00228 dst[id++] = blue[cid];
00229 }
00230 }
00231 }
00232
00233
00234 void
00235 extractRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap, const MatrixDouble &m) {
00236
00237 int nrPixInImage = m.nRow() ;
00238 int i;
00239
00240
00241 for (i=0; i<nrClusters; i++) {
00242
00243 rg[i].size = nrPixInImage ;
00244
00245
00246 rg[i].nrPixInGroup = 0;
00247
00248
00249
00250 rg[i].idx = VectorDouble(nrPixInImage);
00251
00252 rg[i].dim = m.nCol();
00253
00254 }
00255
00256
00257 for (i=0; i<nrPixInImage; i++) {
00258
00259 int l = clusterMap[i];
00260
00261 int n = rg[l].nrPixInGroup ;
00262 rg[l].idx[n] = i;
00263
00264 rg[l].nrPixInGroup++;
00265 }
00266
00267
00268 for (i=0; i<nrClusters; i++) {
00269 rg[i].label = i;
00270
00271 MatrixDouble rgPixMatrix(rg[i].nrPixInGroup, m.nCol(), 0.0) ;
00272
00273 for(int j=0;j<rg[i].nrPixInGroup;j++)
00274 rgPixMatrix.setRow(j, m.getRow(rg[i].idx[j]) ) ;
00275
00276
00277 VectorDouble means ;
00278 rg[i].cov = rgPixMatrix.covariance(means) ;
00279
00280 rg[i].m = means ;
00281 }
00282 }
00283
00284 double regionDistance(pixgroup &u, pixgroup &v) {
00285
00286
00287
00288
00289 MatrixDouble covSums = (u.cov + v.cov) ;
00290 MatrixDouble invCovSums = covSums.i() ;
00291
00292 int nDim = u.m.Size() ;
00293
00294 VectorDouble vec(nDim) ;
00295
00296 for (int i=0; i<nDim; i++) {
00297 vec[i] = v.m[i] - u.m[i];
00298 }
00299
00300 double res = vec * invCovSums * vec ;
00301
00302 return res ;
00303 }
00304
00305
00306 void rgmerge(pixgroup &toMerge, pixgroup ®ion, const MatrixDouble &m) {
00307
00308 int i ;
00309 int nn = toMerge.nrPixInGroup+region.nrPixInGroup;
00310
00311
00312 int j=0 ;
00313 for (i=toMerge.nrPixInGroup;i<nn; i++)
00314 toMerge.idx[i] = region.idx[j++];
00315
00316 toMerge.nrPixInGroup = nn;
00317
00318 region.nrPixInGroup = 0;
00319
00320
00321 MatrixDouble rgPixMatrix(nn, m.nCol(), 0.0) ;
00322
00323
00324 for(i=0;i<nn;i++)
00325 rgPixMatrix.setRow(i, m.getRow(toMerge.idx[i]) ) ;
00326
00327
00328 VectorDouble mean ;
00329 toMerge.cov = rgPixMatrix.covariance(mean) ;
00330
00331 toMerge.m = mean ;
00332 }
00333
00334
00335 void
00336 mergeRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap,
00337 const MatrixDouble &m, int width, int height, double threshold)
00338 {
00339
00340 MatrixTem<int> adjacencyTable(nrClusters, nrClusters, 0) ;
00341
00342
00343
00344 for (int y=1; y<height-1; y++) {
00345 for (int x=1; x<width-1; x++) {
00346
00347 int i = y*width + x;
00348 for (int u=-1; u<=1; u++)
00349 for (int v=-1; v<=1; v++) {
00350 int nb = i+u*width+v;
00351 if (clusterMap[i] != clusterMap[nb])
00352 adjacencyTable[(int)clusterMap[i]][(int)clusterMap[nb]] =
00353 adjacencyTable[(int)clusterMap[nb]][(int)clusterMap[i]] = 1;
00354 }
00355 }
00356 }
00357
00358 double minval = 1E10;
00359 int umin, vmin;
00360
00361
00362 MatrixDouble distanceTable(nrClusters, nrClusters,0.0) ;
00363 for (int u=0; u<nrClusters; u++) {
00364 for (int v=u+1; v<nrClusters; v++) {
00365 if (adjacencyTable[u][v]>0) {
00366 double d = regionDistance(rg[u],rg[v]) ;
00367 distanceTable[u][v] = distanceTable[v][u] = d;
00368
00369 if (d < minval) {
00370 minval = d; umin = u; vmin = v;
00371 }
00372 }
00373 }
00374 }
00375
00376 bool ismerged = (minval <= threshold);
00377
00378
00379 while (ismerged) {
00380
00381
00382 rgmerge(rg[umin],rg[vmin], m);
00383
00384
00385
00386 adjacencyTable[umin][vmin] = adjacencyTable[vmin][umin] = 0;
00387 for (int i=0; i<nrClusters; i++) {
00388 if(adjacencyTable[vmin][i]>0) {
00389
00390 adjacencyTable[i][umin] = adjacencyTable[umin][i] = 1;
00391 adjacencyTable[i][vmin] = adjacencyTable[vmin][i] = 0;
00392 }
00393
00394
00395 }
00396
00397
00398
00399 bool foundmin = false;
00400 for (int i=0; i<nrClusters; i++) {
00401 if (i!=umin && adjacencyTable[umin][i]>0) {
00402 double d = regionDistance(rg[umin],rg[i]) ;
00403 distanceTable[umin][i] = distanceTable[i][umin] = d;
00404 if (d <= minval) {
00405 minval = d; vmin = i;
00406
00407
00408 }
00409 }
00410 }
00411
00412 if (!foundmin)
00413 {
00414 minval = 1E10;
00415 for (int u=0; u<nrClusters; u++)
00416 for (int v=u+1; v<nrClusters; v++)
00417 if (adjacencyTable[u][v] && (distanceTable[u][v]<minval)) {
00418 minval = distanceTable[u][v]; umin = u; vmin = v;
00419 }
00420 }
00421 if (minval > threshold) ismerged = false;
00422
00423 }
00424
00425
00426 int count = 0;
00427
00428 for (int c=0; c<nrClusters; c++) {
00429 if (rg[c].nrPixInGroup!=0) {
00430
00431 for (int i=0; i<rg[c].nrPixInGroup; i++)
00432 clusterMap[rg[c].idx[i]] = count;
00433
00434 rg[c].label = count++;
00435 }
00436 }
00437 }
00438
00439
00440
00441
00442 void updateMean(MatrixDouble &m, pixgroup &rg, long inew)
00443 {
00444
00445 int dim = m.nCol();
00446 long n = rg.nrPixInGroup - 1;
00447 int i ;
00448
00449
00450
00451
00452
00453
00454 for (i=0; i<dim; i++) rg.m[i] *= n;
00455
00456 for (i=0; i<dim; i++)
00457 {
00458 rg.m[i] += m[inew][i];
00459
00460
00461
00462
00463 }
00464
00465 for (i=0; i<dim; i++) rg.m[i] /= rg.nrPixInGroup;
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475 }
00476
00477 #define BORDER -4
00478 #define PROCESSED -3
00479 #define INQUEUE -2
00480 #define UNLABELED -1
00481
00482
00483 void relabel(VectorDouble &clusterMap, int width, int height, pixgroup *rg,
00484 int &nrClusters, int numpixrm, long sizethreshold)
00485 {
00486 long imsize = width*height;
00487 int x,y;
00488 long i;
00489 int d[4] = {1, -width, -1, width};
00490
00491
00492 for (int l=numpixrm; l>0; l--)
00493 {
00494
00495 VectorDouble pixstate(clusterMap) ;
00496
00497
00498 i = width ;
00499 for (y=1; y<height-1; y++) {
00500
00501 i++;
00502 for (x=1; x<width-1; x++) {
00503 if (clusterMap[i]!=UNLABELED) {
00504
00505 bool boundary = false;
00506 for (int u=0; u<4; u++) {
00507
00508
00509 if (clusterMap[i+d[u]]!=clusterMap[i]) {
00510 boundary = true; break;
00511 }
00512 }
00513 if (boundary) pixstate[i] = UNLABELED;
00514 }
00515 i++;
00516 }
00517
00518
00519 i++;
00520 }
00521
00522
00523 clusterMap = pixstate;
00524 }
00525
00526
00527 for (i=0; i<nrClusters; i++) rg[i].nrPixInGroup = 0;
00528 for (i=0; i<imsize; i++)
00529 if (clusterMap[i] != UNLABELED) {
00530 int lb = clusterMap[i];
00531 rg[lb].idx[rg[lb].nrPixInGroup] = i;
00532 rg[lb].nrPixInGroup++;
00533 }
00534
00535
00536 int count = 0;
00537 for (int c=0; c<nrClusters; c++) {
00538 if (rg[c].nrPixInGroup >= sizethreshold) {
00539 for (i=0; i<rg[c].nrPixInGroup; i++) clusterMap[rg[c].idx[i]] = count;
00540 count++;
00541 }
00542 else for (i=0; i<rg[c].nrPixInGroup; i++) clusterMap[rg[c].idx[i]] = UNLABELED;
00543 }
00544
00545
00546 nrClusters = count;
00547 }
00548
00549
00550
00551 void refinery(VectorDouble &clusterMap, MatrixDouble &m, pixgroup *rg, int nrClusters,
00552 int width, int height)
00553 {
00554
00555
00556 int imsize = width*height;
00557 int i;
00558 int x,y,k;
00559 int d[4] = {1, -width, -1, width};
00560 VectorDouble pixstate(imsize);
00561
00562 VectorDouble vec(m.nCol());
00563 double val ;
00564 std::list<long> boundarylist;
00565 std::list<long> *mylist = new std::list<long>[imsize];
00566
00567
00568 for (i=0; i<nrClusters; i++) rg[i].nrPixInGroup = 0;
00569
00570 for (i=0; i<width; i++) pixstate[i] = BORDER;
00571
00572 for (y=1; y<height-1; y++)
00573 {
00574 pixstate[i++] = BORDER;
00575 for (x=1; x<width-1; x++)
00576 {
00577 if (clusterMap[i]!=UNLABELED)
00578 {
00579 int l = clusterMap[i];
00580 long n = rg[l].nrPixInGroup;
00581 rg[l].idx[n] = i;
00582 rg[l].nrPixInGroup++;
00583 pixstate[i] = PROCESSED;
00584 }
00585 else
00586 {
00587 int adjflag = 0;
00588 for (int u=0; u<4; u++) if (clusterMap[i+d[u]]!=UNLABELED)
00589 {
00590 boundarylist.push_back(i);
00591 pixstate[i] = clusterMap[i+d[u]];
00592 adjflag = 1;
00593 break;
00594 }
00595 if (!adjflag) pixstate[i] = UNLABELED;
00596 }
00597 i++;
00598 }
00599 pixstate[i++] = BORDER;
00600 }
00601 for (;i<imsize; i++) pixstate[i] = BORDER;
00602
00603 for (i=0; i<nrClusters; i++)
00604 {
00605 rg[i].label = i;
00606
00607 MatrixDouble rgPixMatrix(rg[i].nrPixInGroup,m.nCol());
00608 for(int j=0;j<rg[i].nrPixInGroup;j++)
00609 rgPixMatrix.setRow(j, m.getRow(rg[i].idx[j]) ) ;
00610
00611 VectorDouble means ;
00612 rg[i].cov = rgPixMatrix.covariance(means) ;
00613 rg[i].m = means ;
00614
00615 }
00616
00617
00618 int counter=0 ;
00619 while(!boundarylist.empty())
00620 {
00621 i = boundarylist.front();
00622
00623 boundarylist.pop_front();
00624 clusterMap[i] = pixstate[i];
00625
00626 int l = clusterMap[i];
00627 for (k=0; k<m.nCol(); k++) {
00628 vec[k] = rg[l].m[k] - m[i][k];
00629
00630 }
00631 val = vec * vec ;
00632 int b = (int)(val*200);
00633 if (b>=imsize) b=imsize-1;
00634 mylist[b].push_back(i);
00635 counter++;
00636 }
00637
00638 int nmin = 0;
00639 while (1)
00640 {
00641 for (; nmin<imsize && mylist[nmin].empty();nmin++);
00642 if (nmin==imsize) break;
00643
00644 while(!mylist[nmin].empty())
00645 {
00646 i = mylist[nmin].front();
00647 mylist[nmin].pop_front();
00648 clusterMap[i] = pixstate[i];
00649 pixstate[i] = PROCESSED;
00650
00651 int l = clusterMap[i];
00652 rg[l].idx[rg[l].nrPixInGroup] = i;
00653 rg[l].nrPixInGroup++;
00654 updateMean(m, rg[l],i);
00655 for (int u=0; u<4; u++) if (pixstate[i+d[u]]==UNLABELED)
00656 {
00657 int newid = i + d[u];
00658 for (k=0; k<m.nCol(); k++) {
00659 vec[k] = rg[l].m[k] - m[newid][k];
00660 }
00661
00662 val = vec * vec ;
00663 int b = (int)(val*200);
00664 if (b>=imsize) b=imsize-1;
00665 mylist[b].push_back(newid);
00666 pixstate[newid] = clusterMap[i];
00667 if (b<nmin) nmin = b;
00668 }
00669 }
00670 }
00671
00672 delete[] mylist;
00673
00674 }
00675
00676
00677 #undef UNLABELED
00678 #undef INQUEUE
00679 #undef PROCESSED
00680 #undef BORDER
00681
00682 }
00683 }
00684 }
00685
00686 #endif