fams.cpp

mean-shift算法的实现算例

/////////////////////////////////////////////////////////////////////////////
// Name:        fams.cpp
// Purpose:     fast adaptive meanshift implementation
// Author:      Ilan Shimshoni
// Modified by: Bogdan Georgescu
// Created:     08/14/2003
// Version:     v0.1
/////////////////////////////////////////////////////////////////////////////

#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <math.h>
#include <memory.h>
#include <time.h>
#include "fams.h"


FAMS::FAMS(int no_lsh)
{
   hasPoints_ = 0;
   nsel_ = 0;
   npm_ = 0;
   hashCoeffs_ = NULL;
   noLSH_=no_lsh;

   time_t tt1;
   time(&tt1);
   srand((unsigned)tt1);
//   srand(100);

   nnres1_ = 0;
   nnres2_ = 0;
}

FAMS::~FAMS()
{
   CleanData();
   
}

void FAMS::CleanData()
{
   CleanPoints();
   CleanSelected();
   CleanPrunedModes();
   CleanHash();
}

void FAMS::CleanPoints()
{
   if (hasPoints_)
   {
      delete [] points_;
      delete [] data_;
      delete [] rr_;
      hasPoints_ = 0;
   }
}

void FAMS::CleanSelected()
{
   if (nsel_ > 0)
   {
      delete [] psel_;
      delete [] modes_;
      delete [] hmodes_;
      nsel_ = 0;
   }
}

void FAMS::CleanPrunedModes()
{
   if (npm_ > 0)
   {
      delete [] prunedmodes_;
      delete [] nprunedmodes_;
      npm_ = 0;
   }
}

void FAMS::CleanHash()
{
   if (hashCoeffs_ != NULL)
   {
      delete [] hashCoeffs_;
      hashCoeffs_ = NULL;
   }
}

int FAMS::LoadPoints(char* filename)
{
   bgLog("Load data points from %s...",filename);
   CleanPoints();
   FILE* fd;
   fd = fopen(filename, "r");

   if (!fd)
   {
      bgLog("Error opening %s\n", filename);
      return 1;
   }
   fscanf(fd, "%d %d", &n_, &d_);
   if ((n_<1) || (d_<1))
   {
      bgLog("Error reading %s\n", filename);
      fclose(fd);
      return 1;
   }

   // allocate data
   float* pttemp;
   pttemp = new float[n_*d_];
   int i;
   for (i=0; (i<(n_*d_)) && (fscanf(fd, "%g", &pttemp[i]) == 1); i++);
   fclose(fd);

   if (i!= (n_*d_))
   {
      bgLog("Error reading %s\n", filename);
      delete [] pttemp;
      return 1;
   }

   // allocate and convert to integer
   for (i=0, minVal_=pttemp[0], maxVal_=pttemp[0]; i<(n_*d_); i++)
   {
      if (minVal_>pttemp[i])
         minVal_ = pttemp[i];
      else if (maxVal_<pttemp[i])
         maxVal_ = pttemp[i];
   }
   data_ = new unsigned short[n_*d_];
   rr_ = new double[d_];
   hasPoints_ = 1;
   float deltaVal = maxVal_-minVal_;
   if (deltaVal == 0) deltaVal = 1;
   for (i=0; i<(n_*d_); i++)
      data_[i] = (unsigned short) (65535.0*(pttemp[i]-minVal_)/deltaVal);
   delete [] pttemp;
   dataSize_ = d_*sizeof(unsigned short);

   points_ = new fams_point[n_];
   unsigned short* dtempp;
   for (i=0, dtempp=data_; i<n_; i++, dtempp+=d_)
   {
      points_[i].data_ = dtempp;
      points_[i].usedFlag_ = 0;
   }
   bgLog("done\n");
   return 0;
}


// Choose a subset of points on which to perform the mean shift operation 
void FAMS::SelectMsPoints(double percent,int jump)
{
   if (hasPoints_ == 0)
      return;
   int i;
   int tsel;
   if(percent){
      tsel = (int)(n_*percent/100.0);
      if (tsel != nsel_)
      {
         CleanSelected();
         nsel_ = tsel;
         psel_ = new int[nsel_];
         modes_ = new unsigned short[nsel_*d_];
         hmodes_ = new unsigned int[nsel_];
      }
      for(i=0; i<nsel_;  i++)
         psel_[i] = MyMin(n_-1,(int)(drand48()*n_));
   }
   else
   {
      tsel = (int)ceil(n_/(jump+0.0));
      if (tsel != nsel_)
      {
         CleanSelected();
         nsel_ = tsel;
         psel_ = new int[nsel_];
         modes_ = new unsigned short[nsel_*d_];
         hmodes_ = new unsigned int[nsel_];
      }
      for(i=0; i<nsel_; i++)
         psel_[i]=i*jump;
   }
}

void FAMS::MakeCuts(fams_partition* cuts)
{
   int i;
   for(i=0; i<L_; i++)
      MakeCutL(cuts[i]);
}

// data-driven uniform partition

void FAMS::MakeCutL(fams_partition& cut)
{
   int n1 = (int)floor(K_/(1.0*d_));
   int i,j;
   int ncu=0;
   int w;
   for(i=0; i<d_; i++)
   {
      for(j=0; j<n1; j++)
      {
         cut[ncu].which_ = i;
         w = MyMin((int)(drand48()*n_),n_-1);
         cut[ncu].where_ = points_[w].data_[i];
         ncu++;
      }
   }
   /*
   int *which = new int[d_];
   for (i=0; i<d_; i++)
      which[i]=i;
   int wh, ndleft, itmp;
   ndleft = d_;
   for (i=0; i<(K_-ncu); i++)
   {
      wh = MyMin((int)(drand48()*ndleft),ndleft-1);
      itmp = which[wh+i];
      which[wh+i] = which[i];
      which[i] = itmp;
      ndleft--;
   }

   for(i=0 ; ncu < K_; ncu++, i++)
   {
      w = MyMin((int)(drand48()*n_),n_-1);
      cut[ncu].which_ = which[i];
      cut[ncu].where_ = points_[w].data_[which[i]];
   }
   delete [] which;
   */
   int *which = new int[d_];
   memset(which,0,sizeof(int)*d_);
   for (; ncu<K_; )
   {
      int wh = MyMin((int)(drand48()*d_),d_-1);
      if(which[wh]!=0)
         continue;
      which[wh]=1;
      int w = MyMin((int)(drand48()*n_), n_-1);
      cut[ncu].which_ = wh;
      cut[ncu].where_ = points_[w].data_[wh];
      ncu++;
   }
   delete [] which;
}

void FAMS::InitHash(int nk)
{
   CleanHash();
   hashCoeffs_ = new int[nk];
   for(int i=0; i<nk; i++)
      hashCoeffs_[i] = rand();
}

/* compute the hash key and and the double hash key if needed 
It is possible to give M the the size of the hash table and 
hjump for the double hash key 
*/

int FAMS::HashFunction(int *cutVals, int whichPartition, int kk,int M,int *hjump)
{
   int i;
   int res = whichPartition;
   for(i=0; i<kk; i++)
   {
      res += cutVals[i]*hashCoeffs_[i];
   }
   if(M)
   {
      res = abs(res);
      if(hjump)
         *hjump = (res%(M-1))+1;
      res = res%M;
   }
   return res;
}


// Add a point to the LSH hash table using double hashing 
void FAMS::AddDataToHash(block HT[],int hs[],fams_point& pt,int where,int Bs,int M,int which,int which2,
                      int hjump)
{
   int nw=0;
   for(;;where = (where+hjump)%M)
   {
      nw++;
      if(nw > M)
      {
         bgLog("LSH hash table overflow exiting\n");
         exit(-1);
      }
      if(hs[where] == Bs)
         continue;
      HT[where][hs[where]].pt_ = &pt;
      HT[where][hs[where]].whichCut_ = which;
      HT[where][hs[where]].which2_ = which2;
      hs[where]++;
      break;
   }
}

// compute the pilot h_i's for the data points
void FAMS::ComputePilot(block *HT,int *hs, fams_partition *cuts, char* pilot_fn)
{
   const int win_j = 10,max_win=7000;
   int i,j;
   unsigned int nn;
   unsigned int wjd = (unsigned int)(win_j*d_);
   int num_l[1000];
   fams_res_cont res(n_);
   if(LoadBandwidths(pilot_fn)==0)
   {
      bgLog("compute bandwidths...");
      for(j=0; j<n_; j++)
      { 
         int numn=0;
         int numns[max_win/win_j];
         memset(numns,0,sizeof(numns));
         int nel;
         if(noLSH_)
         {
            nel = n_;
            for(i=0; i<nel; i++)
            {
               fams_pointp pt = &points_[i];
               nn = DistL1(points_[j],*pt) / wjd;
               if(nn <max_win/win_j)
                  numns[nn]++;	
            }
         }
         else
         {
            GetNearestNeighbours(points_[j],HT,hs,cuts,res,0,num_l);
            nel = res.nel_;      
            for(i=0; i<nel; i++)
            {
               fams_pointp pt = res.vec_[i];
               nn = DistL1(points_[j],*pt) / wjd;
               if(nn <max_win/win_j)
                  numns[nn]++;	
            }
         }
         for(nn=0; nn<max_win/win_j; nn++)
         {
            numn+=numns[nn];
            if(numn>k_)
            {
               break;
            }
         }	
         points_[j].window_=(nn+1)*wjd;
      }
      SaveBandwidths(pilot_fn);
   }
   else
      bgLog("load bandwidths...");
   for(j=0; j<n_; j++){
      points_[j].weightdp2_ = (float) pow(FAMS_FLOAT_SHIFT/points_[j].window_, (d_+2)*FAMS_ALPHA);
//      points_[j].weightdp2_ = (float) (1.0/pow(points_[j].window_, (d_+2)*FAMS_ALPHA));
   }
}

// compute real bandwiths for selected points
void FAMS::ComputeRealBandwidths(unsigned int h)
{
   const int win_j = 10,max_win=7000;
   int i,j;
   unsigned int nn;
   unsigned int wjd;
   wjd = (unsigned int) (win_j*d_);
   int who;
   if (h==0)
   {
      for(j=0; j<nsel_; j++)
      {
         who = psel_[j];
         int numn=0;
         int numns[max_win/win_j];
         memset(numns,0,sizeof(numns));
         for(i=0; i<n_; i++)
         {
            fams_pointp pt = &points_[i];
            nn = DistL1(points_[who],*pt) / wjd;
            if(nn <max_win/win_j)
               numns[nn]++;	
         }
         for(nn=0; nn<max_win/win_j; nn++)
         {
            numn+=numns[nn];
            if(numn>k_)
            {
               break;
            }
         }	
         points_[who].window_=(nn+1)*win_j;
      }
   } else
   {
      for(j=0; j<nsel_; j++)
      {
         who = psel_[j];
         points_[who].window_=h;
      }
   }
}

// compute the pilot h_i's for the data points
void FAMS::ComputeScores(block *HT,int *hs, fams_partition *cuts, float* scores)
{
   const int win_j = 10,max_win=7000;
   int i,j,who;
   unsigned int nn;
   unsigned int wjd = (unsigned int)(win_j*d_);
   int num_l[1000];
   memset(scores, 0, L_*sizeof(float));
   fams_res_cont res(n_);
   for(j=0; j<nsel_; j++)
   { 
      who = psel_[j];
      int numn;
      int nl=0;
      int numns[max_win/win_j];
      memset(numns,0,sizeof(numns));
      int nel;
      if(noLSH_)
      {
         nel = n_;
         num_l[L_]=n_+1;
         for(i=0; i<nel; i++)
         {
            fams_pointp pt = &points_[i];
            nn = DistL1(points_[who],*pt) / wjd;
            if(nn <max_win/win_j)
               numns[nn]++;	
            if(i == (num_l[nl]-1))
            {
               numn=0;
               for(nn=0; nn<max_win/win_j; nn++)
               {
                  if(numn>k_)
                     break;
               }
               for(; (num_l[nl]-1) == i; nl++)
                  scores[nl] += (float) (((nn+1.0)*win_j)/points_[who].window_);
            }
         }
      }
      else
      {
         GetNearestNeighbours(points_[who],HT,hs,cuts,res,0,num_l);
         nel = res.nel_;
         num_l[L_] = n_+1;
         for(i=0; i<nel; i++)
         {
            fams_pointp pt = res.vec_[i];
            nn = DistL1(points_[who],*pt) / wjd;
            if(nn <max_win/win_j)
               numns[nn]++;
            if(i == (num_l[nl]-1))
            {
               numn=0;
               for(nn=0; nn<max_win/win_j; nn++)
               {
                  numn+=numns[nn];
                  if(numn>k_)
                     break;
               }
               for(; (num_l[nl]-1) == i; nl++)
                  scores[nl] += (float) (((nn+1.0)*win_j)/points_[who].window_);
            }
         }
      }
      numn=0;
      for(nn=0; nn<max_win/win_j; nn++)
      {
         numn+=numns[nn];
         if(numn>k_)
            break;
      }	
   }
   for(j=0; j<L_; j++)
      scores[j]/=nsel_;

}

/*
Insert an mean-shift result into a second hash table so when another mean shift computationis 
performed about the same C_intersection region, the result can be retreived without further 
computation
*/

void FAMS::InsertIntoHash(block2* HT2,int *hs2,int where,int which,
                    unsigned short **solution,int M,int hjump)
{
   int nw=0;
   for(;;where = (where+hjump)%M)
   {
      nw++;
      if(nw==M)
      {
         bgLog("Second Hash Table Full\n");
         exit(-1);
      }