segm.cc

A general technique for the recovery of signi cant image features is presented. The technique is ba

//Color Image Segmentation
//This is the implementation of the algorithm described in 
//D. Comaniciu, P. Meer, 
//Robust Analysis of Feature Spaces: Color Image Segmentation,
//http://www.caip.rutgers.edu/~meer/RIUL/PAPERS/feature.ps.gz
//appeared in Proceedings of CVPR'97, San Juan, Puerto Rico.
// ===========================================================================
// =====      Module: segm.cc
// ===== -------------------------------------------------------------- ======
// =====      Version 01   Date: 04/22/97
// ===== -------------------------------------------------------------- ======
// ===========================================================================
// =====      Written by Dorin Comaniciu
// =====      e-mail:  [e-mail]
// ===========================================================================
// Permission to use, copy, or modify this software and its documentation
// for educational and research purposes only is hereby granted without
// fee, provided that this copyright notice appear on all copies and
// related documentation.  For any other uses of this software, in original
// or modified form, including but not limited to distribution in whole
// or in part, specific prior permission must be obtained from
// the author(s).
//
// THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND,
// EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
// WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
//
// IN NO EVENT SHALL RUTGERS UNIVERSITY BE LIABLE FOR ANY SPECIAL,
// INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, OR ANY
// DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
// WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY
// THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE
// OR PERFORMANCE OF THIS SOFTWARE.
// ===========================================================================
//

#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <fstream.h>
#include <strstream.h>
#include <memory.h>
#include <time.h>
#include "segm.hh"


int option = 2;
static const int rect_gen_contor=3;

// Radius of the searching window
static float gen_RADIUS[3]={2, 3, 4};
static float rect_RADIUS[rect_gen_contor]={8, 6, 4};
static float fix_RADIUS[rect_gen_contor];
static float final_RADIUS;
static float RADIUS2;
static float RADIUS;

static float max_dist;

static int my_threshold[3]={50, 100, 400};
static int my_rap[3]={4, 2, 1};
static int act_threshold;

#define my_abs(a) ((a) > 0 ? (a): (-a)) 
#define SQRT2 1.4142
#define SQRT3 1.7321
static const float	BIG_NUM      =	1.0e+20;

// Coefficient matrix for xyz and rgb spaces
static const int    XYZ[3][3] = { { 4125, 3576, 1804 },
                                  { 2125, 7154,  721 },
                                  {  193, 1192, 9502 } };
static const float  RGB[3][3] = { { 3.2405, -1.5371, -0.4985 },
                                  {-0.9693,  1.8760,  0.0416 },
                                  { 0.0556, -0.2040,  1.0573 } };
 
// Constants for LUV transformation 
static const float     Xn = 0.9505;
static const float     Yn = 1.0;
static const float     Zn = 1.0888;
static const float     Un_prime = 0.1978;
static const float     Vn_prime = 0.4683;
static const float     Lt = 0.008856;

// # of samples
static const int	Max_J	     =	25;

// Limit of the number of failed trials
static const int       Max_trials    =  50;

// Defaults values for the parameters.
static const int	sam_max = 60;

// Few more trials at the end
static const int        MAX_TRIAL=10;

// Used in 3-D histogram computation
static const int FIRST_SIZE=262144; // 2^18
static const int SEC_SIZE=64;  // 2^6

// Make coordinate computation faster
// my_neigh is for auto_segm, my_neigh_r works with region
static int my_neigh[8];
static int my_neigh_r[8];


// Results
static int ORIG_COLORS;
static int SEGM_COLORS;

Boolean my_write_PGM_file( char*, Octet*,int,int);
void    covariance_w(const int N, int M, const int p, int ** data, 
                   int *w, float T[],  float C[p_max][p_max]);
void    mean_s(const int N, const int p, int J[], int **data, float T[]);
void    my_convert(int, float *, int *);
void    reverse_map(Octet *, Octet *, int *, Octet *, float T[][p_max]);

// Class constructor
SegmenterMS::SegmenterMS( )
{
  _p = 0;
  _p_ptr = 0;
  _rrows  = 0;
  _rcolms = 0;
  _data_all = nil;
  _data = nil;
  _NJ = Max_J;
}

// Class destructor
SegmenterMS::~SegmenterMS( )
{
  if ( _data_all ) {
    for ( register int i = 0; i < _p; i++ )
      if ( _data_all[i] )   delete [] _data_all[i];
    delete [] _data_all;
  }
  if ( _data ) {
    for ( register int i = 0; i < _p; i++ )
      if ( _data[i] )   delete [] _data[i];
    delete [] _data;
  }
}

// LUV (final_T[]) to RGB (TI[]) conversion
void my_convert(int selects, float *final_T, int *TI)
{
  // this condition is always true
  if ( selects & Lightness && selects & Ustar && selects & Vstar )
    {
      if(final_T[0]<0.1)
        {
	  TI[0]=0;TI[1]=0;TI[2]=0;
	}
      else
	{
          float my_x, my_y, my_z;
	  if(final_T[0]< 8.0)
            my_y = Yn * final_T[0] / 903.3;
          else
            my_y = Yn * pow((final_T[0] + 16.0) / 116.0, 3);

          float u_prime = final_T[1] / (13 * final_T[0]) + Un_prime;
          float v_prime = final_T[2] / (13 * final_T[0]) + Vn_prime;

          my_x = 9 * u_prime * my_y / (4 * v_prime);
          my_z = (12 - 3 * u_prime - 20 * v_prime) * my_y / (4 * v_prime);

          TI[0] =int((RGB[0][0]*my_x + RGB[0][1]*my_y + RGB[0][2]*my_z)*255.0);
          TI[1] =int((RGB[1][0]*my_x + RGB[1][1]*my_y + RGB[1][2]*my_z)*255.0);
          TI[2] =int((RGB[2][0]*my_x + RGB[2][1]*my_y + RGB[2][2]*my_z)*255.0);
 
          if(TI[0]>255) TI[0]=255;
          else if(TI[0]<0) TI[0]=0;
	  
	  if(TI[1]>255) TI[1]=255;
          else if(TI[1]<0) TI[1]=0;
	  
	  if(TI[2]>255) TI[2]=255;
	  else if(TI[2]<0) TI[2]=0;
	}
    }
  else
    {
       TI[0]=(int)final_T[0];
       TI[1]=(int)final_T[1];
       TI[2]=(int)final_T[2];
    }
}

// RGB to LUV conversion
// To gain speed the conversion works on a table of colors (_col_RGB[])
// rather than on the whole image
void SegmenterMS::convert_RGB_LUV( RasterIpChannels* signal, long selects )
{
  int x, y, z, my_temp;
  
  float l_star, u_star, v_star;
  float u_prime, v_prime;
  register int temp_col, temp_index, temp_ind;
  register int j,k;
  
  int a00=XYZ[0][0], a01=XYZ[0][1], a02=XYZ[0][2];
  int a10=XYZ[1][0], a11=XYZ[1][1], a12=XYZ[1][2];
  int a20=XYZ[2][0], a21=XYZ[2][1], a22=XYZ[2][2];
  
  int *A00 = new int[MAXV]; int *A01 = new int[MAXV]; int *A02 = new int[MAXV];
  int *A10 = new int[MAXV]; int *A11 = new int[MAXV]; int *A12 = new int[MAXV];
  int *A20 = new int[MAXV]; int *A21 = new int[MAXV]; int *A22 = new int[MAXV];
  
  for(j=0; j<MAXV;j++)
    {
      A00[j]=a00*j; A01[j]=a01*j; A02[j]=a02*j;
      A10[j]=a10*j; A11[j]=a11*j; A12[j]=a12*j;
      A20[j]=a20*j; A21[j]=a21*j; A22[j]=a22*j;
    }
 
  float *my_pow = new float[MAXV];
  for(j=0; j<MAXV;j++)
    my_pow[j]= 116.0 * pow(j/255.0, 0.3333333) - 16;
 
  Octet* temp_ch0 = signal->chdata(0);
  Octet* temp_ch1 = signal->chdata(1);
  Octet* temp_ch2 = signal->chdata(2);

  int pp;    
  int *temp0, *temp1, *temp2;
  pp = _p_ptr;
  if ( selects & Lightness ) temp0 = _data_all[pp++];
  if ( selects & Ustar ) temp1 = _data_all[pp++];
  if ( selects & Vstar ) temp2 = _data_all[pp++]; 
  _p_ptr=pp;

  for ( j = 0; j < _n_colors; j++)
    {
        temp_col=_col_RGB[j];
	int R=temp_col>>16; int G=(temp_col>>8) & 255; int B=temp_col & 255;

        x = A00[R] + A01[G] + A02[B];
        y = A10[R] + A11[G] + A12[B];
        z = A20[R] + A21[G] + A22[B];
	
        float  tval = y / 2550000.0; //Yn==1
	if ( tval >  Lt)  l_star = my_pow[(int)(tval*255+0.5)];
        else  l_star = 903.3 * tval;

        my_temp = x + 15 * y + 3 * z;
	if(my_temp)
	  {
	    u_prime = (float)(x << 2) / (float)(my_temp);
	    v_prime = (float)(9 * y) / (float)(my_temp);
	  }
	else
	  {
	    u_prime=4.0;
	    v_prime=9.0/15.0;
	  }
	
	tval=13*l_star;
        u_star = tval* (u_prime - Un_prime); // Un_prime = 0.1978
        v_star = tval* (v_prime - Vn_prime); // Vn_prime = 0.4683

	_col0[j] = (int)(l_star+0.5);
	if(u_star>0) _col1[j] = (int)(u_star+0.5);
        else _col1[j] = (int)(u_star-0.5);
  
	if(v_star>0) _col2[j] = (int)(v_star+0.5);
        else _col2[j] = (int)(v_star-0.5);
    }
  for(j=0;j<_ro_col;j++)
    {
      temp_col=(((((int)temp_ch0[j])<<8)+(int)temp_ch1[j])<<8)+(int)temp_ch2[j];
      temp_ind=_col_misc[temp_col>>6];
      for(k=temp_ind;k<temp_ind+SEC_SIZE;k++)
	if(_col_RGB[k]==temp_col)
	  {
           temp_index=_col_index[j]=k;
           break;
	 }
      temp0[j]=_col0[temp_index];
      temp1[j]=_col1[temp_index];
      temp2[j]=_col2[temp_index];
    }   
  delete [] my_pow;
  delete [] A22; delete [] A21; delete [] A20;
  delete [] A12; delete [] A11; delete [] A10;
  delete [] A02; delete [] A01; delete [] A00;
  delete [] _col_misc;
  delete [] _col_RGB;
  cerr<<":";
}

// 3-D Histogram computation
// Implement a trade-off between speed and required memory
void SegmenterMS::my_histogram(RasterIpChannels* signal, long selects)
{
  int *first_tab= new int[FIRST_SIZE];
  _col_misc= new int[FIRST_SIZE]; 
  int **third_tab;
  int *fourth_tab;
  int *fifth_tab=new int[SEC_SIZE];
  _n_colors=0;
  
  register int k,j,p,r;
  int temp_ind, sec_ind, third_ind;
  
  int first_contor=0, third_contor=0;
  
  memset(first_tab,0,sizeof(int)*FIRST_SIZE);
  memset(_col_misc,0,sizeof(int)*FIRST_SIZE);
  
  register Octet* ch0 = signal->chdata(0);
  register Octet* ch1 = signal->chdata(1);
  register Octet* ch2 = signal->chdata(2);
  
  //first_tab -> how many
  for(k=0;k<_ro_col;k++)
    {
      temp_ind=(((ch0[k]<<8)+ch1[k])<<2)+(ch2[k]>>6);
      first_tab[temp_ind]++;
    }
  //_col_misc -> memo position
  for(k=0;k<FIRST_SIZE;k++)
    if(first_tab[k])
      {
        _col_misc[k]=first_contor;
	first_contor++;
      }
  //contors
  fourth_tab=new int[first_contor];
  memset(fourth_tab,0,sizeof(int)*first_contor);
  //tab of pointers to reduced colors
  third_tab=new int *[first_contor]; 
  first_contor=0;
  for(k=0;k<FIRST_SIZE;k++)
    if(first_tab[k])
      { 
        third_tab[first_contor]=new int[first_tab[k]];       
	first_contor++;
      }
        
  for(k=0;k<_ro_col;k++)
    {
      temp_ind=(((ch0[k]<<8)+ch1[k])<<2)+(ch2[k]>>6);
      sec_ind=ch2[k] & 63;
      third_ind=_col_misc[temp_ind];      
      third_tab[third_ind][fourth_tab[third_ind]]=sec_ind;
      fourth_tab[third_ind]++;      
    }
  for(k=0;k<first_contor;k++)
    {
      memset(fifth_tab,0,sizeof(int)*SEC_SIZE);
      for(j=0;j<fourth_tab[k];j++)
        fifth_tab[third_tab[k][j]]++;
      for(j=0;j<SEC_SIZE;j++)
        if(fifth_tab[j])
	   _n_colors++;
      }
  _col_RGB=new int[_n_colors];
  _m_colors=new int[_n_colors];
  
  k=0;p=0;
  for(r=0;r<FIRST_SIZE;r++)
    if(first_tab[r])
     {
       memset(fifth_tab,0,sizeof(int)*SEC_SIZE);
       for(j=0;j<fourth_tab[k];j++)
         fifth_tab[third_tab[k][j]]++;
       _col_misc[r]=p;
       for(j=0;j<SEC_SIZE;j++)
         if(fifth_tab[j])
	   {	     
	     _col_RGB[p]=(r<<6)+j;
	     _m_colors[p]=fifth_tab[j];
             p++;
	   }
       delete [] third_tab[k];
       k++;
     }
  delete [] third_tab;
  delete [] fourth_tab;
  delete [] fifth_tab;
  delete [] first_tab;
  
  _col_all = new int*[3];
  _col0=_col_all[0] = new int[_n_colors];
  _col1=_col_all[1] = new int[_n_colors];
  _col2=_col_all[2] = new int[_n_colors];
  _col_index = new int[_ro_col];
  cerr<<":";
}    

// Update _col_remain[], _m_col_remain, and _n_col_remain
void SegmenterMS::my_actual(Octet *my_class)
{ 
  register int i;
  int temp_contor=n_remain;
  register int *temp_rem= new int[_ro_col];
  memcpy(temp_rem,gen_remain,sizeof(int)*temp_contor);
  n_remain=0;	
  for(i=0;i<temp_contor;i++)
    if(!my_class[temp_rem[i]])
      gen_remain[n_remain++]=temp_rem[i];
  delete [] temp_rem;
  memset(_col_remain,0,sizeof(int)*_n_col_remain); 
  memset(_m_col_remain,0,sizeof(int)*_n_colors);  
  _n_col_remain=0; 
  for(i=0;i<n_remain;i++)
    _m_col_remain[_col_index[gen_remain[i]]]++;
  for(i=0;i<_n_colors;i++)
     if(_m_col_remain[i])
       {
	 _col_remain[_n_col_remain]=i;
         _n_col_remain++;
       }
}

// if more than "how_many"  neighbors, consider the point
void SegmenterMS::test_neigh(Octet* my_class, int *selected, int* my_contor, int how_many)
{ 
  register int i,j,p,k;
  register Octet* local_class=my_class;
  register int temp_contor=*my_contor;
  register int my_index;
  if(auto_segm) my_index=n_remain;
  else          my_index=_ro_col;
  for ( p = 0, i; p < my_index; p++ )
    {
      if(auto_segm) i=gen_remain[p];
      else          i=p;
      if(!local_class[i])
	{
          int neigh_contor=0, no_neigh=1;
          for(j=0;j<8;j++)
            {
              k=i+my_neigh[j];
              if(k>=0 && k<_ro_col && local_class[k])
	        {
                  if(auto_segm  && gen_class[k]!=255) continue;
	          neigh_contor++;
                  if(neigh_contor>how_many)
		    {
		      no_neigh=0;
                      break;
		    }
	        }
	    }
        if(!no_neigh)
	  {
            if(auto_segm) selected[*my_contor]=i;
	    *my_contor=*my_contor+1;
	  }
	}
    }
  for(i=temp_contor;i<*my_contor;i++)
    local_class[selected[i]]=1;
}


// Find the feature vectors inside the given window
// Use Improved Absolute Error Inequality Criterion 
// when computing Euclidean distance
// See J.S.Pan el al, Fast Clustering Alg. for VQ, Pattern Recognition,
// Vol. 29, No. 3, pp. 511-518, 1996