segm.cc

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

  for ( i = 0; i < p; i++ )
    TT[i] = 0;
  for ( i = 0; i < N; i++ )
    {
      k = J[i];
      for ( j = 0; j < p; j++ )
        TT[j] += data[j][k];
     }
   for ( i = 0; i < p; i++ )
      T[i] = (float)TT[i] / (float)N;
}

// Build a subsample set of 9 points
int SegmenterMS::subsample(float *Xmean )
{
  int J[9];
  register int my_contor=0, uj, i0;
  if(auto_segm)
    i0=J[my_contor]=
       gen_remain[int(float(n_remain)*float(rand())/float(SHRT_MAX))]; 
  else
    i0=J[my_contor]=int(float(_n_points)*float(rand())/float(SHRT_MAX));
  my_contor++;
  for(register i=0;i<8;i++){
    uj=i0 + my_neigh_r[i];
    if(uj>=0 && uj<_n_points)
      {
        if((auto_segm && gen_class[uj]!=255)) break;
	else
          {        
            J[my_contor] = uj;
            my_contor++;
	  }      
      }
    }
  mean_s(my_contor, _p, J, _data, Xmean);
  return 1;
} 

// Sampling routine with all needed tests
float SegmenterMS::my_sampling( int rect, float T[Max_rects][p_max])
{ 
  register int  k, c;
  register float L,U,V,Res;
  register float my_dist=max_dist, my_sqrt_dist=fix_RADIUS[0];
  float	TJ[Max_J][p_max];
  int	l =  0;       //contor of number of subsample sets
  int	ll = 0;       //contor of trials
  float	Xmean[p_max];
  float	Obj_fct[Max_J];
 
  //Max_trials = max number of failed trials
  //_NJ = max number of subsample sets
   
  while ( (ll < Max_trials) && (l < _NJ ) )
    {
      if ( subsample(Xmean) )    // the subsample procedure succeeded
	{ 
	      ll = 0; c=0;

	      // Save the mean
	      for ( k = 0; k < _p; k++ ) TJ[l][k] = Xmean[k];

	      // Compute the square residuals (Euclid dist.)
	      if(auto_segm)
                  for ( register int p = 0; p < _n_col_remain; p++ )
		    {
		     k=_col_remain[p];
                     L=_col0[k]-Xmean[0]; if(my_abs(L)>=my_sqrt_dist) continue;
		     U=_col1[k]-Xmean[1]; if(my_abs(U)>=my_sqrt_dist) continue;
		     V=_col2[k]-Xmean[2]; if(my_abs(V)>=my_sqrt_dist) continue;
		     if(L*L+U*U+V*V<my_dist) c+=_m_col_remain[k];
		    }
	      else
	          for ( k = 0; k < _n_points; k++ )
		    {
                     L=_data[0][k]-Xmean[0]; 
                     if(my_abs(L)>=my_sqrt_dist) continue; 
		     U=_data[1][k]-Xmean[1]; 
                     if(my_abs(U)>=my_sqrt_dist) continue; 
		     V=_data[2][k]-Xmean[2]; 
                     if(my_abs(V)>=my_sqrt_dist) continue; 
                     if(L*L+U*U+V*V<my_dist) c++;
		    }
	      
              //  Objective functions
	      Obj_fct[l]=c;
              l++;
	    }
      else ++ll;	
    }
  if ( ll == Max_trials && l < 1) return( BIG_NUM ); // Cannot find a kernel

  //  Choose the highest density  
  L = -BIG_NUM; c=0;
  for ( k = 0; k < _NJ; k++ )
      if ( Obj_fct[k] > L)
	{
	  L = Obj_fct[k];
	  c = k;
	}
  if(Obj_fct[c]>0)
    for(k=0;k<_p;k++)
       T[rect][k]=TJ[c][k];
  else return -BIG_NUM; // Not enough points
  return ( 0 );
}

// Compute the weighted covariance of N points 
void covariance_w(const int N, int M, const int p,  int **data, 
               int *w, float T[], float C[p_max][p_max])
{
  register int i, j, k, l;
  int TT[p_max];
  for ( i = 0; i < p; i++ )
    TT[i] = 0;
  for ( i = 0; i < M; i++ )
    for ( j = 0; j < p; j++ )
      TT[j] += w[i]*data[j][i]; 
  for ( i = 0; i < p; i++ )
      T[i] = (float) TT[i] / (float)N;

  for ( i = 0; i < p; i++ )
    for ( j = i; j < p; j++ )
      C[i][j] = 0.0;
  for ( i = 0; i < M; i++ )
    {
      for ( k = 0; k < p; k++ )
	for ( l = k; l < p; l++ )
	   C[k][l]+=w[i]*(data[k][i]-T[k])*(data[l][i]-T[l]);
    }
  for ( k = 0; k < p; k++ )
    {
    for ( l = k; l < p; l++ )
      C[k][l] /= (float)(N-1);
    for ( l = 0; l < k; l++ )
      C[k][l] = C[l][k];
    }
}

// initialization
void SegmenterMS::init_neigh(void)
{
  my_neigh[0]= -_colms-1;  my_neigh[1]= -_colms;
  my_neigh[2]= -_colms+1;  my_neigh[3]= +1;
  my_neigh[4]= +_colms+1;  my_neigh[5]= +_colms;
  my_neigh[6]= +_colms-1;  my_neigh[7]= -1;

  my_neigh_r[0]= -_rcolms-1;  my_neigh_r[1]= -_rcolms;
  my_neigh_r[2]= -_rcolms+1;  my_neigh_r[3]= +1;
  my_neigh_r[4]= +_rcolms+1;  my_neigh_r[5]= +_rcolms;
  my_neigh_r[6]= +_rcolms-1;  my_neigh_r[7]= -1;
}

// Init matrices parameters
void SegmenterMS::init_matr(void)
{
  // General statistic parameters for X.
  float	Mg[p_max];         //sample mean of X
  float	C[p_max][p_max];   //sample covariance matrix of X
  covariance_w(_ro_col, _n_colors, _p, _col_all, _m_colors, Mg, C);

  // Adaptation
  float my_th=C[0][0]+C[1][1]+C[2][2];
  int active_gen_contor=1;
  
  if(auto_segm)
     fix_RADIUS[0]=gen_RADIUS[option]*sqrt(my_th/100);
  else
    {  
      active_gen_contor=rect_gen_contor;
      for(int i=0;i<active_gen_contor;i++)
        fix_RADIUS[i]=rect_RADIUS[i]*sqrt(my_th/100);
    }
  final_RADIUS=fix_RADIUS[active_gen_contor-1]*1.26;
  max_dist=fix_RADIUS[0]*fix_RADIUS[0];
#ifdef TRACE
  printf("\n %.2f %.2f ", fix_RADIUS[0], final_RADIUS);
#endif
  act_threshold=(int)((my_threshold[option]>sqrt(_ro_col)/my_rap[option])? 
                 my_threshold[option]:sqrt(_ro_col)/my_rap[option]);
  cerr<<":";
} 

// Init
void SegmenterMS::initializations(RasterIpChannels* pic, sRectangle rects[], int *n_rects, long selects, int *active_gen_contor)
{
  register int i; 
  XfRaster::Info        info;
  pic->raster_info(info);
  _colms = info.columns; _rows  = info.rows; _ro_col = _rows * _colms;

  _data_all = new int*[_p];
  for ( i = 0; i < _p; i++ )
      _data_all[i] = new int[_ro_col];
  _data0=_data_all[0]; _data1=_data_all[1]; _data2=_data_all[2];
  init_neigh();
  my_histogram(pic, selects);
  convert_RGB_LUV( pic, selects );
  gen_class = new Octet[_ro_col];
  memset(gen_class,255,_ro_col);
  if(!(*n_rects))
    {
      auto_segm=1;
      *n_rects=Max_rects;
      n_remain=_ro_col;
      _n_col_remain=_n_colors;
      gen_remain = new int[_ro_col];
      _col_remain = new int[_n_colors];
      _m_col_remain = new int[_n_colors];
      for ( i = 0; i< _ro_col ; i++ )
        gen_remain[i] = i;
      for ( i = 0; i< _n_colors ; i++ )
        _col_remain[i] = i;
      memcpy(_m_col_remain,_m_colors,sizeof(int)*_n_colors);

      for ( i = 0; i < Max_rects ; i++ )
        {
          rects[i].width=_colms; rects[i].height=_rows;
          rects[i].x = 0;        rects[i].y = 0;
        }
      *active_gen_contor=1;
    }
  else
    {
      auto_segm=0;
      n_remain=0;
      *active_gen_contor=rect_gen_contor;
      option=2;
    }
  init_matr();
  delete [] _m_colors;
}

// Mean shift segmentation applied on the selected region(s) of an image or
// on the whole image
Boolean SegmenterMS::ms_segment( RasterIpChannels* pic, sRectangle rects[],
				     int n_rects, long selects , 
                                     unsigned int seed_default, Boolean block)
{
  int t=time(0);
  if (n_rects > Max_rects) return false;
  if ( selects & Lightness || selects & Ustar || selects & Vstar ) 
    _p=3;
  else return false;

  int contor_trials=0, active_gen_contor;
  float	T[Max_rects][p_max];
  int TI[Max_rects][p_max];
  float L,U,V,RAD2,q;
  register int i,k;
  
  srand( seed_default ); cerr<<":";
  initializations(pic, rects, &n_rects, selects, &active_gen_contor);

  // Mean shift algorithm and removal of the detected feature
  int rect;
  for ( rect = 0; rect < n_rects; rect++ )
    {
      _n_points = rects[rect].width * rects[rect].height;
      cut_rectangle( &rects[rect] );
     
    MS:
      if(auto_segm && contor_trials) _NJ=1;
      
      q =  my_sampling( rect, T);
      if(q == -BIG_NUM || q == BIG_NUM)
	  if(contor_trials++<MAX_TRIAL)  goto MS;
	  else                           break;
      
      float final_T[p_max];
      for ( i = 0; i < _p; i++ ) final_T[i]=T[rect][i];

      int *selected = new int[_ro_col];
      Octet *sel_col = new Octet[_n_colors];
      Octet *my_class = new Octet[_ro_col];
       
      int my_contor=0, gen_gen_contor=-1, how_many=10; 
      while(gen_gen_contor++<active_gen_contor-1)
        {
          set_RADIUS(gen_gen_contor, 0);
	  int gen_contor=0; 
          Boolean last_loop=0;
	 
          while(gen_contor++<how_many)
            { 
              if(auto_segm)
		{
                  memset(sel_col,0,_n_colors);
                  new_auto_loop(final_T,sel_col); 
                  L=T[rect][0]-final_T[0]; U=T[rect][1]-final_T[1];
                  V=T[rect][2]-final_T[2]; RAD2=L*L+U*U+V*V; 
                  if(RAD2<0.1){
                    my_contor=0;
                    memset(my_class,0,_ro_col);     
                    for(k=0;k<n_remain;k++)
		      {
                        register int p=gen_remain[k];
			if(sel_col[_col_index[p]])
			  {
			    selected[my_contor++]=p;
                            my_class[p]=1;
			  }
		      }
                    break;
		    }
		  else
                  {
	            T[rect][0]=final_T[0];T[rect][1]=final_T[1];
                    T[rect][2]=final_T[2];
                   }
	         }
	       else
		 {
		    my_contor=0;
                    memset(my_class,0,_ro_col);     
	            nauto_loop(final_T, selected, my_class,&my_contor);
                    mean_s(my_contor, _p, selected, _data, final_T);
                    L=T[rect][0]-final_T[0]; U=T[rect][1]-final_T[1];
                    V=T[rect][2]-final_T[2]; RAD2=L*L+U*U+V*V; 
    
                    if(RAD2<0.1)
                      break;
                    else
                     {
	               T[rect][0]=final_T[0];T[rect][1]=final_T[1];
                       T[rect][2]=final_T[2];
                     } 
		  }
	    }
	}
      if(auto_segm)
	{
          if(option==0) test_neigh(my_class, selected, &my_contor,2);
          else if(option==1) test_neigh(my_class, selected, &my_contor,1);
          else if(option==2) test_neigh(my_class, selected, &my_contor,0);
	}      
#ifdef TRACE   
      printf("my_contor = %d contor_trials = %d",my_contor,contor_trials);  
#endif	
      if(!auto_segm)
        if(test_same_cluster(rect,T))
	  {
             delete [] my_class; delete [] sel_col; delete [] selected;
             continue;
	  }
      if(auto_segm && my_contor<act_threshold)
        {
          delete [] my_class; delete [] sel_col; delete [] selected;
          if(contor_trials++<MAX_TRIAL) goto MS;
          else                          break;
        }
      
      if(auto_segm) my_actual(my_class);

      my_convert(selects, final_T, TI[rect]);  
      for ( k = 0; k < _ro_col; k++ )
        if(my_class[k])
          gen_class[k]=rect;
     
      delete [] my_class; delete [] sel_col; delete [] selected;
    }     			
  n_rects=rect;
  Octet** isegment = new Octet*[3];	
  register int j;
  for ( j = 0; j < 3; j++ )
    {
      isegment[j] = new Octet[_ro_col];
      memset( isegment[j], 0, _ro_col );
    }
  Octet *isegment0=isegment[0]; Octet *isegment1=isegment[1];
  Octet *isegment2=isegment[2];
 
  if(auto_segm)
    {
      Octet *my_class = new Octet[_ro_col];
      for ( k = 0; k < _ro_col; k++ )
         if(gen_class[k]==255) gen_class[k]=n_rects;

      Octet inv_map[Max_rects];
      for(k=0;k<n_rects;k++) inv_map[k]=k;

      if(option==0)  conn_comp(my_class,&n_rects,inv_map,T,10,0);
      else if(option==1) conn_comp(my_class,&n_rects,inv_map,T,20,0); 
      else conn_comp(my_class,&n_rects,inv_map,T,act_threshold,0); 
      optimize(T,n_rects);
      for(k=0;k<n_rects;k++) my_convert(selects, T[k], TI[k]);
      
      // Postprocessing
      if(option==1 || option ==2)
          for(k=2;k<10;k+=2) conn_comp(my_class,&n_rects,inv_map,T,k,1);
#ifdef TRADE
      printf("\nTime = %d seconds\n",time(0)-t);
#endif 
      register Octet my_max_ind; 
      for ( k = 0; k < _ro_col; k++ )
        {
          if((my_max_ind=gen_class[k])!=n_rects)
	    {
              int *ti=TI[my_max_ind];
	      isegment0[k]=ti[0]; isegment1[k]=ti[1]; isegment2[k]=ti[2];
	    } 
	}

      // Save the borders
      memset(my_class,255,_ro_col); 
      for ( k = 0; k < _ro_col; k++ )
	{
          int pres_class=gen_class[k];
          int pos_colms=k%_colms;
	      
          if(pos_colms==0 || pos_colms==(_colms-1) ||
             k<(_colms-1) || k>(_ro_col-_colms))
             my_class[k] = 0;
          else
             for(j=1;j<8;j+=2)
               {
                 int u=k+my_neigh[j];
                 if(u>=0 && u<_ro_col &&  (pres_class<gen_class[u]))
		   {
		     my_class[k] = 0; break;
		   }
	        }
	}
      my_write_PGM_file( "result.pgm", my_class,_rows,_colms);
      delete [] my_class;
    }
  else   // if not auto_segmentation
    {
      optimize(T,n_rects);
      for(k=0;k<n_rects;k++)
        my_convert(selects, T[k], TI[k]);
      register int temp_class;
      for (k=0;k<_ro_col;k++)
	  {
            if((temp_class=gen_class[k])<n_rects)
	      {
                 isegment0[k] = /*TI[temp_class][0];*/pic->chdata(0)[k];
	         isegment1[k] = /*TI[temp_class][1];*/pic->chdata(1)[k];
	         isegment2[k] = /*TI[temp_class][2];*/pic->chdata(2)[k];
	       }
	    } 
     } 
  if(auto_segm){
    delete [] _m_col_remain; 
    delete [] _col_remain;
    delete [] gen_remain;
  }
  delete [] gen_class;
  delete [] _col_index;
  delete [] _col0; delete [] _col1; delete [] _col2; 
  delete [] _col_all;
  
  XfRaster::Info	info;
  pic->raster_info(info);
  result_ras_ = new RasterIpChannels( info, 3, eDATA_OCTET,
				      isegment, true );
  cerr << endl << "Original Colors: " << _n_colors << endl;
  cerr << "Segment. Colors: " << n_rects << endl;
  
  return true;  
}