Mean shift implementation in C++

Question

Can anyone recommend a lightweight mean shift clustering implementation in C++? I am already using OpenCV, however their mean shift implementation is for tracking, not clustering. I have seen EDISON, however, this is for image segmentation and not clustering.

I could implement it myself, however would rather not invest the time, and not take the risk of bugs.

Thanks

Answer 1

Here is my C++ version of mean shift object tracking method. To run the code successfully, you need to add the tf.h header file to the main code directory.

#include "tf.h" // include the header file
using namespace cv;
using namespace std;
#include <stdio.h>

makerect mkr;  // rectangle for encompassing target 
// rcs for row coordination of target window center
//ccs for column coordination of target window center
double rcs=0,ccs=0; 
// w for width of target window 
// l for length of target window
double W=70,L=60;

const int mySizes[3]={16,16,16}; // vector for number of histogram bins
cv::Mat q4hsv = Mat::zeros(3,mySizes,CV_64F); // initializing histogram variable
uchar nbins=16; // var for num of histogram bins
int main(void){     
    printf("enter 14 or 36: \t"); // enter number of input video name
    uint flag; // var for video flag or number
    cin>>flag;
    Mat ref4frame; // reference frame which is used for initializing mean shift parameters
char filename [50]; 
sprintf(filename,"im%d.avi",flag);
VideoCapture capture(filename);
if( !capture.isOpened() )
    throw "Error when reading steam_avi";

unsigned long int f4counter=0; // frame counter var
histhsv hsv; // instantiating histhsv class
for (f4counter=1;f4counter<=40000000;f4counter++){ // a big number to support many frames
    capture >> ref4frame; //reading input image from specified directory
    if( !(ref4frame.data ))   // checking the read status
    {       
       printf("Cannot load file image %s\n", filename);    
        break;    } 
    uchar ndiv = uchar(256/nbins); // division value to which intesity values are divided

if (f4counter==1)   { // special for 1st frame
    char modelname[20];
    if(flag==36){
    sprintf(modelname,"data%d.png",flag);       
    }
    else if (flag==14){
        sprintf(modelname,"data%d.jpg",flag);   
    }
    // imread is defined in tf.h
    Mat img = imread(modelname,1);//reads 1st image
    if( !(img.data ))//check if file loading is ok
    {       
       printf("Cannot load file image %s\n", modelname);       
        break;    } 
    hsv.img=img;//assign new image to hsv object (calculates hsv or rgb hist)
    hsv.run();//run the histogram calculator
// assign temporary hsv object to reference hist q4hsv object
for (int i=0;i<16;i++){
        for(int j=0;j<16;j++){
            for(int k=0;k<16;k++){
                q4hsv.at<double>(i,j,k)=hsv.q[i][j][k];
            }
        }
}
}
if(f4counter<5){averageglobalhsv(ref4frame,q4hsv,rcs,ccs);}
averagelocalhsv(ref4frame,q4hsv,rcs,ccs);//normalizing histogram values (0-1)

    Point pt1; pt1.x=ccs; pt1.y=rcs;
    int thickness=4;//thickness of marker - here is a circle
    int lineType=8;
    int shift=0;
    RNG rng(0xFFFFFFFF);
    cv::circle(ref4frame, pt1, 5, randomColor(rng), thickness, lineType,
shift); //marking object center with a circle
    myimshow("reference frame",ref4frame);//show current image
    waitKey(1); 
}

int c=waitKey(0);
//release memory
ref4frame.release();
destroyAllWindows();
    return 0;
}

here is the tf.h header file contents:

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc_c.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2\opencv.hpp"
#include "core\core.hpp"
#include <cstdio>
#include <iostream>
#include <fstream>
#include <math.h>
using namespace cv;
using namespace std;

// makerect class: to create the desired size encompassing window
class makerect
{
public:
    double rcs,ccs,w,l; // ctl row, ctl column, width, length
    uint rmin,rmax,cmin,cmax,height,length;//min max coordination vars
    void run(void);
    makerect();
};
makerect::makerect(){
    rmin=0;rmax=0;cmin=0;cmax=0;rcs=0;ccs=0;w=0;l=0;
}
void makerect::run(){
//not that all points must be either integer or floating point
    rmin=uint(rcs-floor(w/2));//min row coordination
    rmax=uint(rmin+w);//max row coordination
    cmin=uint (ccs-floor(l/2));//min column coordination
    cmax=uint(cmin+l);//max column coordination
    if(cmax>length){cmax=length;cmin=cmax-l;// checking column to be inside image}
    if(rmax>height){rmax=height;rmin=rmax-w;//checking row to be inside image}
}

static Scalar randomColor(RNG& rng)
{
    int icolor = (unsigned)rng;
    return Scalar(icolor&255, (icolor>>8)&255, (icolor>>16)&255);
}
//myimshow: is a function to to show image
void myimshow(char* name4window,Mat &tmp4image){
    namedWindow(name4window,1); imshow(name4window,tmp4image);
}


void averageglobalhsv(Mat ref4frame,Mat &f,double &rcs,double &ccs)
{
    Mat img;
    cvtColor(ref4frame,img,CV_BGR2HSV);//rgb2hsv conversion
    uint n4rowsref=ref4frame.rows;// num of rows
    uint n4colsref=ref4frame.cols;// num of cols
    double *im4bp = new double [n4rowsref*n4colsref];//1-D dynamic array equal to image pixels
    uint nbins=16;// num of histogram bins
    uint ndiv=256/nbins; //division value to which intensities are divided
//vars for image moments (zero to second moments)
double m00=0,m01=0,m10=0,m20=0,m02=0,m11=0,w=0; 

uint R=0,G=0,B=0; //red bin num, green bin num, blue bin num
    for(uint i=0;i<n4rowsref;i++){ //loop through all image rows
        for(uint j=0;j<n4colsref;j++){//loop through all image columns
            Vec3b inten=img.at<Vec3b>(i,j);//a vector of three element
            R=inten.val[2]; G=inten.val[1];  B=inten.val[0];
            R/=ndiv; G/=ndiv;  B/=ndiv;//calculating the bin to which current pixel intensity belong

                im4bp[i*n4colsref+j]=1.3*sqrt(f.at<double>(R,G,B));//calculating spatial weighted kernel histogram formula

        }
    }

    for(uint i=0;i<n4rowsref;i++){//loop through all image rows
        for(uint j=0;j<n4colsref;j++){//loop through all image columns
            w=im4bp[j+n4colsref*i];// weight for pixel at (i,j)
            m01=m01+double(i)*w;//first moment on y/row coordination
            m10=m10+double(j)*w;//first moment on x/column coordination
            m00=m00+w;//zeroth moment which is sum of all moments
        }
    }
if(m00>0){  
rcs=ceil(m01/m00);//central point for row
ccs=ceil(m10/m00);}//central point for column
delete im4bp;//cleaning memory
}


void averagelocalhsv(Mat ref4frame,Mat &f,double &rcs,double &ccs)
{
    Mat img;
cvtColor(ref4frame,img,CV_BGR2HSV); 
makerect mkr;
int sz[]={2,2};
    uint n4rowsref=ref4frame.rows;
    uint n4colsref=ref4frame.cols;
    double *im4bp = new double [n4rowsref*n4colsref];
    uint nbins=16;
    uint ndiv=256/nbins;
double m00=0,m01=0,m10=0,m20=0,m02=0,m11=0,w=0,Dxx,Dyy,Dxy;

uint R=0,G=0,B=0;
    for(uint i=0;i<n4rowsref;i++){
        for(uint j=0;j<n4colsref;j++){
            Vec3b inten=img.at<Vec3b>(i,j);
            R=inten.val[2]; G=inten.val[1];  B=inten.val[0];
            R/=ndiv; G/=ndiv;  B/=ndiv;
            im4bp[i*n4colsref+j]=1.3*sqrt(f.at<double>(R,G,B));
        }
    }


for(int iter=1;iter<5;iter++){
    mkr.rcs=rcs;mkr.ccs=ccs;mkr.w=100;mkr.l=100;mkr.height=ref4frame.rows;
    mkr.length=ref4frame.cols;  mkr.run();

    m00=0;m01=0;m10=0;m20=0;m02=0;m11=0;
    for(uint i=mkr.rmin;i<mkr.rmax;i=i+1){
        for(uint j=mkr.cmin;j<mkr.cmax;j=j+1){
            w=im4bp[j+n4colsref*i];
            m01=m01+double(i)*w;
            m10=m10+double(j)*w;
            m00=m00+w;
        }
    }   
}
if(m00>0){
rcs=ceil(m01/m00);
ccs=ceil(m10/m00);
}
delete im4bp;
}
class histhsv{
public:
    histhsv();
    void run(void);
    Mat img;    
    double q[16][16][16];

};
histhsv::histhsv(){};
void histhsv::run(void){
//Mat hsv4image;
double sum4hist=0;
uchar nbins=16;
uchar ndiv=256/nbins;
double wmax =0;;    
double r_center=0;
double c_center =0;
r_center = img.rows/2;
c_center = img.cols/2;
for (int i=0;i<nbins;i++){for(int j=0;j<nbins;j++){for(int k=0;k<nbins;k++){
q[i][j][k]=0;   }   }   };
    cvtColor(img,img,CV_BGR2HSV);       
    int H=0,S=0,V=0;    
    for(int r=0;r<img.rows;r++){
        for(int c=0;c<img.cols;c++){
            Vec3b intensity = img.at<Vec3b>(r,c);
            H=intensity.val[0]/ndiv;
            S=intensity.val[1]/ndiv;
            V=intensity.val[2]/ndiv;
            q[H][S][V]+=wmax-(pow(r-r_center,2)+pow(c-c_center,2));
            sum4hist+=q[H][S][V];
        }
    }
    for (int i=0;i<nbins;i++){
        for(int j=0;j<nbins;j++){
            for(int k=0;k<nbins;k++){
                q[i][j][k]/=sum4hist;
            }
        }
    }       

}

Answer 2

This is old, but I am working with mean shift right now so I thought it best to answer.

I think I understand the distinction you are making here, but when you say you are looking for mode detection this is vague in the technical sense as from the point of view of the algorithm as the algorithm inherently is for searching for "modes", which are the local minima or maxima depending on how you frame the optimization problem (Gradient descent or ascent).

This source, which was found on the EDISON site, claims to be a c++ implementation of the mean shift clustering algorithm, but as discussed above, clustering is the main implementation of the mode seeking behavior that all other uses of mean shift is based on, especially segmentation, so you can certainly use the EDISON source to find a clustering implementation, even if you have to search through it a bit.

I also found this Github project, for what it is worth, but I haven't worked with it before.

LAST NOTE: I also noticed you said "lightweight" implementation. Note that mean shift is not a very efficient algorithm (i think it is something like O(N^3), but I will check that). That said, it can still be efficiently implemented, though how that should be gauged is more ambiguous. Needless to say, Quick Shift, an attempt by UCLA researchers to solve the issues of the more efficient medoid shift, a similar non-parametric mode seeking algorithm, might be more like what you are looking for in a "lightweight" algorithm.