|| || || || || || || ||

---

ColorSegmentation.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_ColorSegmentation_h
#define Impala_Core_Array_ColorSegmentation_h

#include <list>
#include "Core/Matrix/MatrixTem.h"
#include "Core/Array/ArrayListExport.h"
#include "Core/Array/ArrayListDelete.h"
#include "Core/Array/ColorSegmentationInvariant.h"
#include "Core/Array/ColorSegmentationAlgorithm.h"
#include "Core/Array/ColorGaussResponses.h"
#include "Core/Array/ColorGaborResponses.h"
#include "Basis/Timer.h"
#include "Core/Matrix/MatKMeans.h"
#include "Core/Array/WriteRaw.h"
#include "Core/Matrix/MatTranspose.h"
#include "Core/Matrix/MatKLM.h"
#include "Core/Matrix/MatFeatureNormalization.h"
#include "Core/Array/GetRgbPixels.h"

namespace Impala
{
namespace Core
{
namespace Array
{


typedef Impala::Core::Vector::VectorTem<double> VectorDouble;
typedef Impala::Core::Matrix::MatrixTem<double> MatrixDouble;
using namespace Impala::Core::Matrix;

// define a region as a group of pixels
class pixgroup{
public:

    int size; // set to number of pixels in image, but not used
    int nrPixInGroup;
    int label;
    int dim; // dim after PCA
    MatrixDouble cov; //covariance matrix
    VectorDouble m; //mean
    VectorDouble idx; // index for every pixel, index is actually an offset
};


void labelimg(unsigned char* dst, const VectorDouble &imap, int width, int height) ;

void updateMean(MatrixDouble &m, pixgroup &rg, long inew) ;

void extractRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap,
                    const MatrixDouble &m) ;

double regionDistance(pixgroup &u, pixgroup &v) ;

void rgmerge(pixgroup &toMerge, pixgroup &region, const MatrixDouble &m) ;

void mergeRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap, 
                  const MatrixDouble &m, int width, int height, double threshold) ;

void relabel(VectorDouble &clusterMap, int width, int height, pixgroup *rg,
             int &nrClusters, int numpixrm, long sizethreshold) ;

void refinery(VectorDouble &clusterMap, MatrixDouble &m, pixgroup *rg, int nrClusters,
              int width, int height) ;


inline void
ColorSegmentation(Array2dVec3UInt8*& dst, Array2dVec3Real64* im,
                    ColorSegmentationAlgorithm segAlg, 
                    ColorSegmentationInvariant invariantType, 
                    double minRegionFraction, double threshold,
                    bool useGauss = false,
                    std::vector<Array2dVec3UInt8*>* dispList = 0)
{
bool vv = false;
if (vv) std::cout << "start ColorSegmentation" << std::endl;
    if (dst == 0)
        dst = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
    Timer timer(1);
    timer.Start();
    std::vector<Array2dScalarReal64*> respList;
    if (useGauss)
        respList = ColorGaussResponses(im);
    else
        respList = ColorGaborResponses(im, segAlg, invariantType);
    //WriteRawList(respList, std::string("resp.raw"), true);

if (vv) std::cout << " making response matrix" << std::endl;

    typedef Array2dScalarReal64 Mat;
    typedef VecScalarReal64 Vec;
    int nrFeat = respList.size();
    int nrPixels = respList[0]->CW() * respList[0]->CH();
    int nrPcaDim = 4;

    double* data = ArrayListExport(respList);
    ArrayListDelete(&respList);

    Mat* mTmp = MatCreate<Mat>(nrFeat, nrPixels, data);
    //WriteRaw(mTmp, "abc_feat.raw", true);
    Mat* m = MatTranspose(mTmp);
    delete mTmp;

if (vv) std::cout << "time: " << timer.SplitTime() << " PCA " << std::endl;
    // reduce to nrPcaDim dimensions with PCA (if necessary)
    if (MatNrCol(m) > nrPcaDim)
    {
        mTmp = MatKLM(m, nrPcaDim);
        delete m;
        m = mTmp;
    }

if (vv) std::cout << "time: " << timer.SplitTime() << " feature_normalization " << std::endl;
        // normalize 
     MatFeatureNormalization(m); 

if (vv) std::cout << "time: " << timer.SplitTime() << " kmeans " << std::endl;

        // k-means
    int nrClusters = nrPixels/2500;
    Mat* clusters = MatCreate<Mat>(nrClusters, MatNrCol(m));
    Vec* clusterMap = VecCreate<Vec>(nrPixels);
    VecScalarInt32* vectorsInC = VecCreate<VecScalarInt32>(nrClusters);
    int clustersFound;
    MatKMeans(m, clusters, clusterMap, 25, false, clustersFound, vectorsInC);
    delete clusters;
if (vv)
for (int i=0 ; i<nrClusters ; i++)
{
    int nr = *VecE(vectorsInC, i);
    std::cout << "cluster " << i << " has " << nr//*VecE(vectorsInC, i) 
          << " vectors" << std::endl;
}
    delete vectorsInC;

// mMat is "wrapper" since Matrix dtor commented
MatrixDouble mMat(MatNrRow(m), MatNrCol(m), MatE(m));

    VectorDouble clusterMapMat(VecNrElem(clusterMap), VecE(clusterMap), false);
    delete clusterMap;

if (vv) std::cout << "time: " << timer.SplitTime() << " Doing region stuff " << std::endl;
if (dispList != 0)
{
    Array2dVec3UInt8* dIm = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
    labelimg(dIm->PB(), clusterMapMat, im->CW(), im->CH());
    dispList->push_back(dIm);
}

        // get the regions
    pixgroup *rg = new pixgroup[nrClusters];
    extractRegions(rg, nrClusters, clusterMapMat, mMat);

        // merge regions below threshold
    mergeRegions(rg, nrClusters, clusterMapMat, mMat, im->CW(), im->CH(),
                 threshold) ;
if (dispList != 0)
{
    Array2dVec3UInt8* dIm = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
    labelimg(dIm->PB(), clusterMapMat, im->CW(), im->CH());
    dispList->push_back(dIm);
}

        // UN-label region borders, and the image border
        // minRegionFraction indicates the minimal size of a region in a percentage
    relabel(clusterMapMat, im->CW(), im->CH(), rg, nrClusters, 3,
            nrPixels * minRegionFraction);
if (dispList != 0)
{
    Array2dVec3UInt8* dIm = ArrayCreate<Array2dVec3UInt8>(im->CW(), im->CH());
    labelimg(dIm->PB(), clusterMapMat, im->CW(), im->CH());
    dispList->push_back(dIm);
}

        // region growing
    refinery(clusterMapMat, mMat, rg, nrClusters, im->CW(), im->CH());

        // create colors for a region
    labelimg(dst->CPB(), clusterMapMat, im->CW(), im->CH());

    delete m;
    delete[] rg;
if (vv) std::cout << "time: " << timer.SplitTime() << " Done " << std::endl;
}



void labelimg(unsigned char* dst, const VectorDouble &imap, int width, int height)
{
    //color palette
    unsigned char pl[10][2] = { {255,0},{191,63},{127,0},{255,63},{191,0},
                                {255,127},{63,0}, {191,127},{255,191},{127,63} 
                              };
    unsigned char red[60],green[60],blue[60];

    //create palette
    int id = 0;
    long i;
    for (i=0; i<10; i++)
    {
        red[id] = pl[i][0]; green[id] = pl[i][1]; blue[id] = pl[i][1]; id++;
        red[id] = pl[i][0]; green[id] = pl[i][0]; blue[id] = pl[i][1]; id++;
        red[id] = pl[i][0]; green[id] = pl[i][1]; blue[id] = pl[i][0]; id++;
        red[id] = pl[i][1]; green[id] = pl[i][0]; blue[id] = pl[i][1]; id++;
        red[id] = pl[i][1]; green[id] = pl[i][0]; blue[id] = pl[i][0]; id++;
        red[id] = pl[i][1]; green[id] = pl[i][1]; blue[id] = pl[i][0]; id++;
    }

    long imsize = width*height;
    id = 0;
    
    for (i=0; i<imsize; i++)
    {
        if (imap[i] < 0)
        {
            dst[id++] = 0;
            dst[id++] = 0;
            dst[id++] = 0;
        }
        else
        {
            int cid = ((int)imap[i]) % 60;
            dst[id++] = red[cid];
            dst[id++] = green[cid];
            dst[id++] = blue[cid];
        }
    }
}


void 
extractRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap, const MatrixDouble &m) {

    int nrPixInImage = m.nRow() ;
    int i;

    // init
    for (i=0; i<nrClusters; i++) {

        rg[i].size = nrPixInImage ;

            // number of pixels in the group
        rg[i].nrPixInGroup = 0; 

            // for each pixel in the group, 
            // keep a reference to its place in the original image
        rg[i].idx = VectorDouble(nrPixInImage); 

        rg[i].dim = m.nCol();

    }

        // get regions for each pixel in the image
    for (i=0; i<nrPixInImage; i++) {
            // get the cluster number for pixel i
        int l = clusterMap[i];
            // add a reference to the pixel i
        int n = rg[l].nrPixInGroup ;
        rg[l].idx[n] = i; //add pixel's index
            // added pixel i to group l
        rg[l].nrPixInGroup++;
    }

        // set covariance and mean values for each region
    for (i=0; i<nrClusters; i++) {
        rg[i].label = i;

        MatrixDouble rgPixMatrix(rg[i].nrPixInGroup, m.nCol(), 0.0) ;
            // create a new matrix of the pixels in group i
        for(int j=0;j<rg[i].nrPixInGroup;j++) 
            rgPixMatrix.setRow(j, m.getRow(rg[i].idx[j]) ) ;

            // set covariance matrix, and mean vector
        VectorDouble means ;
        rg[i].cov = rgPixMatrix.covariance(means) ;

        rg[i].m = means ;
    }
}

double regionDistance(pixgroup &u, pixgroup &v) {

    //the statistic distance between region u and v (multivariate)
    //see Hieu's thesis, chapter 3

    MatrixDouble covSums = (u.cov + v.cov) ;
    MatrixDouble invCovSums = covSums.i() ;

    int nDim = u.m.Size() ;

    VectorDouble vec(nDim) ;

    for (int i=0; i<nDim; i++) {
        vec[i] = v.m[i] - u.m[i]; //mean1-mean2
    }

    double res = vec * invCovSums * vec ;

    return res ;
}

// merge regions toMerge and region to the region 'toMerge'
void rgmerge(pixgroup &toMerge, pixgroup &region, const MatrixDouble &m) {
    
    int i ;
    int nn = toMerge.nrPixInGroup+region.nrPixInGroup;

        // add references 
    int j=0 ;
    for (i=toMerge.nrPixInGroup;i<nn; i++) 
        toMerge.idx[i] = region.idx[j++];

    toMerge.nrPixInGroup = nn;

    region.nrPixInGroup = 0;

        //recompute the covariance matrix
    MatrixDouble rgPixMatrix(nn, m.nCol(), 0.0) ;
    
        // create a matrix of the pixels 
    for(i=0;i<nn;i++) 
        rgPixMatrix.setRow(i, m.getRow(toMerge.idx[i]) ) ;

        // set covariance and mean vector
    VectorDouble mean ;
    toMerge.cov = rgPixMatrix.covariance(mean) ;

    toMerge.m = mean ;
}


void 
mergeRegions(pixgroup* rg, int nrClusters, const VectorDouble &clusterMap, 
             const MatrixDouble &m, int width, int height, double threshold)
{

    MatrixTem<int> adjacencyTable(nrClusters, nrClusters, 0) ;

    // initialization of adjacency table
    // if cluster i and j are adjacent to each other then adjacencyTable[i][j] = 1 
    for (int y=1; y<height-1; y++) {
        for (int x=1; x<width-1; x++)   {   
                // convert y and x to pixelnumber i
            int i = y*width + x;
            for (int u=-1; u<=1; u++)
                for (int v=-1; v<=1; v++)   {
                    int nb = i+u*width+v;
                    if (clusterMap[i] != clusterMap[nb])
                        adjacencyTable[(int)clusterMap[i]][(int)clusterMap[nb]] =
 adjacencyTable[(int)clusterMap[nb]][(int)clusterMap[i]] = 1;
                }
        }
    }
    
    double minval = 1E10; //min distance between regions
    int umin, vmin;

        //initialize the joint distance table
    MatrixDouble distanceTable(nrClusters, nrClusters,0.0) ;
    for (int u=0; u<nrClusters; u++) {
        for (int v=u+1; v<nrClusters; v++)  {
            if (adjacencyTable[u][v]>0) { //adjacent
                double d = regionDistance(rg[u],rg[v]) ;
                distanceTable[u][v] = distanceTable[v][u] = d;
                    // find the most similar adjacent regions u and v
                if (d < minval) {
                    minval = d; umin = u; vmin = v;
                }
            }
        }
    }

    bool ismerged = (minval <= threshold);

    //main loop for merging
    while (ismerged) {

            // merge vmin to umin
        rgmerge(rg[umin],rg[vmin], m);

        //int nrRegionsLeft = 0;
            //update the adjacencyTable
        adjacencyTable[umin][vmin] = adjacencyTable[vmin][umin] = 0;
        for (int i=0; i<nrClusters; i++) {
            if(adjacencyTable[vmin][i]>0) {
                    //vmin has been merged to umin
                adjacencyTable[i][umin] = adjacencyTable[umin][i] = 1; 
                adjacencyTable[i][vmin] = adjacencyTable[vmin][i] = 0;
            }
            //if (rg[i].nrPixInGroup > 0)
            //  nrRegionsLeft++;
        }
        //std::cout << "nrRegionsLeft : " << nrRegionsLeft << std::endl;

        //now recompute the statistic distance between new region u and others
        bool foundmin = false;
        for (int i=0; i<nrClusters; i++) {
            if (i!=umin && adjacencyTable[umin][i]>0) {
                double d = regionDistance(rg[umin],rg[i]) ;
                distanceTable[umin][i] = distanceTable[i][umin] = d;
                if (d <= minval) {
                    minval = d; vmin = i; //found another min value
                    // DK : ignore this one so we always find the real min (?)
                    //foundmin = true;
                }
            }
        }

        if (!foundmin) //min value has not been found, find it
        {
            minval = 1E10;
            for (int u=0; u<nrClusters; u++) 
                for (int v=u+1; v<nrClusters; v++)
                    if (adjacencyTable[u][v] && (distanceTable[u][v]<minval)) {
                        minval = distanceTable[u][v]; umin = u; vmin = v;
                    }
        }
        if (minval > threshold) ismerged = false; //no merge anymore

    } //while

    //regenerate new labels
    int count = 0;

    for (int c=0; c<nrClusters; c++) {
        if (rg[c].nrPixInGroup!=0)  {

            for (int i=0; i<rg[c].nrPixInGroup; i++) 
                    clusterMap[rg[c].idx[i]] = count;

            rg[c].label = count++;
        }
    }
}


// update mean vector a region rg 
// when add new pixel to rg (don't have to recompute all
void updateMean(MatrixDouble &m, pixgroup &rg, long inew)
{
    //COVxy = SUM(x.y)/n - Mx.My
    int dim = m.nCol();
    long n = rg.nrPixInGroup - 1; //when not added yet
    int i ;
//  int j;
    //restore :-)
//  for (i=0; i<dim; i++) for (j=i; j<dim; j++)
//  {
//      rg.cov[i][j] = (rg.cov[i][j] + rg.m[i]*rg.m[j])*n;
//  }
    for (i=0; i<dim; i++) rg.m[i] *= n;

    for (i=0; i<dim; i++)
    {
        rg.m[i] += m[inew][i];
//      for (j=i; j<dim; j++)
//      {
//          rg.cov[i][j] += (m[inew][i]*m[inew][j]);
//      }
    }

    for (i=0; i<dim; i++) rg.m[i] /= rg.nrPixInGroup;

//  for (i=0; i<dim; i++)
//  {
//      for (j=i; j<dim; j++)
//      {
//          rg.cov[i][j] = rg.cov[i][j]/rg.nrPixInGroup - rg.m[i]*rg.m[j];
//          rg.cov[j][i] = rg.cov[i][j];
//      }
//  }
}

#define BORDER -4
#define PROCESSED -3
#define INQUEUE -2
#define UNLABELED -1

// unlabel the border of the image, and the pixels between regions 
void relabel(VectorDouble &clusterMap, int width, int height, pixgroup *rg,
             int &nrClusters, int numpixrm, long sizethreshold)
{
    long imsize = width*height;
    int x,y;
    long i;
    int d[4] = {1, -width, -1, width}; 

    //remove pixels at the boundary of a region (thickness = 2*numpixrm)
    for (int l=numpixrm; l>0; l--)
    {
        // make a copy of the clusterMap, to keep the state a pixel is in.
        VectorDouble pixstate(clusterMap) ;

        // start after the borderline (first line of pixels)
        i = width ;
        for (y=1; y<height-1; y++)  {
            // skip left-border pixel
            i++;
            for (x=1; x<width-1; x++) {
                if (clusterMap[i]!=UNLABELED) {

                    bool boundary = false;
                    for (int u=0; u<4; u++) {
                            // if right,top,left or bottom pixel is from another
                            // region, or unlabelled
                        if (clusterMap[i+d[u]]!=clusterMap[i])  {
                        boundary = true; break;
                        }
                    }
                    if (boundary) pixstate[i] = UNLABELED; //remove from region
                }
                i++;
            } // end for(x)

            // skip right-border pixel
            i++;
        } // end for(y)

            // set the result to clusterMap 
        clusterMap = pixstate; 
    }

    // renumber the leftover regions
    for (i=0; i<nrClusters; i++) rg[i].nrPixInGroup = 0;
    for (i=0; i<imsize; i++) 
        if (clusterMap[i] != UNLABELED) {
            int lb = clusterMap[i];
            rg[lb].idx[rg[lb].nrPixInGroup] = i;
            rg[lb].nrPixInGroup++;
        }

    // UNLABEL the regions that are too small 
    int count = 0;
    for (int c=0; c<nrClusters; c++) {
        if (rg[c].nrPixInGroup >= sizethreshold) {
            for (i=0; i<rg[c].nrPixInGroup; i++) clusterMap[rg[c].idx[i]] = count;
            count++;
        }
        else for (i=0; i<rg[c].nrPixInGroup; i++) clusterMap[rg[c].idx[i]] = UNLABELED;
    }

    // store the new number of regions.
    nrClusters = count;
}

// this is near-original code by Minh.
//relaxation: refine the result (VERY DIRTY for now and from here:-) 
void refinery(VectorDouble &clusterMap, MatrixDouble &m, pixgroup *rg, int nrClusters,
              int width, int height)
{

//  printf("refinery") ;
    int imsize = width*height;
    int i;
    int x,y,k;
    int d[4] = {1, -width, -1, width}; 
    VectorDouble pixstate(imsize);

    VectorDouble vec(m.nCol());
    double val ;
    std::list<long> boundarylist;
    std::list<long> *mylist = new std::list<long>[imsize];

    //now initialize for region growing
    for (i=0; i<nrClusters; i++) rg[i].nrPixInGroup = 0;

    for (i=0; i<width; i++) pixstate[i] = BORDER; //over the border

    for (y=1; y<height-1; y++) 
    {
        pixstate[i++] = BORDER;
        for (x=1; x<width-1; x++)
        {
            if (clusterMap[i]!=UNLABELED)
            {
                int l = clusterMap[i];
                long n = rg[l].nrPixInGroup;
                rg[l].idx[n] = i; //add pixel's index
                rg[l].nrPixInGroup++;
                pixstate[i] = PROCESSED;
            }
            else //check adjacency (just on boundary or not)
            {
                int adjflag = 0;
                for (int u=0; u<4; u++) if (clusterMap[i+d[u]]!=UNLABELED)
                {
                    boundarylist.push_back(i); //store the boundary pixel
                    pixstate[i] = clusterMap[i+d[u]];
                    adjflag = 1;
                    break;
                }
                if (!adjflag) pixstate[i] = UNLABELED;
            }
            i++;
        }
        pixstate[i++] = BORDER;
    }
    for (;i<imsize; i++) pixstate[i] = BORDER;

    for (i=0; i<nrClusters; i++)
    {
        rg[i].label = i;

        MatrixDouble rgPixMatrix(rg[i].nrPixInGroup,m.nCol());
        for(int j=0;j<rg[i].nrPixInGroup;j++) 
            rgPixMatrix.setRow(j, m.getRow(rg[i].idx[j]) ) ;

        VectorDouble means ;
        rg[i].cov = rgPixMatrix.covariance(means) ;
        rg[i].m = means ;

    }

    //now initialize for region growing
    int counter=0 ;
    while(!boundarylist.empty())
    {
        i = boundarylist.front();

        boundarylist.pop_front();
        clusterMap[i] = pixstate[i];

        int l = clusterMap[i];
        for (k=0; k<m.nCol(); k++) {
            vec[k] = rg[l].m[k] - m[i][k]; //mean[l]-this pix i

        }
        val = vec * vec ;
        int b = (int)(val*200); //scale
        if (b>=imsize) b=imsize-1;
        mylist[b].push_back(i);
        counter++;
    }
    //now grow the regions
    int nmin = 0;
    while (1)
    {
        for (; nmin<imsize && mylist[nmin].empty();nmin++);
        if (nmin==imsize) break; //done

        while(!mylist[nmin].empty())
        {
            i = mylist[nmin].front();
            mylist[nmin].pop_front();
            clusterMap[i] = pixstate[i];
            pixstate[i] = PROCESSED;
            //update region
            int l = clusterMap[i]; //i is added to region clusterMap[i]
            rg[l].idx[rg[l].nrPixInGroup] = i; //add pixel;
            rg[l].nrPixInGroup++;               
            updateMean(m, rg[l],i);
            for (int u=0; u<4; u++) if (pixstate[i+d[u]]==UNLABELED) //extended
            {
                int newid = i + d[u];
                for (k=0; k<m.nCol(); k++) {
                    vec[k] = rg[l].m[k] - m[newid][k]; //mean[l]-this pix i
                }

                val = vec * vec ;
                int b = (int)(val*200);
                if (b>=imsize) b=imsize-1;
                mylist[b].push_back(newid);
                pixstate[newid] = clusterMap[i];
                if (b<nmin) nmin = b;
            }
        } //while
    }

    delete[] mylist;

}


#undef UNLABELED
#undef INQUEUE
#undef PROCESSED
#undef BORDER

} // namespace Array
} // namespace Core
} // namespace Impala

#endif

---