segmenterms.cpp

这是一个分水岭程序

	register REGION *first_region, *prev_region, *current_region;
	taken = new Octet[_ro_col];
	memset(taken,0,_ro_col);
	conn_selected = new int[_ro_col];
	int *ptr_tab=new int[_ro_col];
	
	for(k=0;k<_ro_col;k++)
		if(!taken[k])
		{	
			current_region=new REGION;
			current_region->my_contor=0;
			current_region->my_class=gen_class[k];
			current_region->my_label=local_label;
			if(k!=0) prev_region->next_region_str=current_region;
			if(k==0){ first_region=current_region;}
			
			local_label++;
			conn_selected[current_region->my_contor++]=k;
			taken[k]=1;
			find_other_neigh(k,ptr_tab,current_region);
			if(change_type==0)
				if(my_max_region[current_region->my_class]<current_region->my_contor)
					my_max_region[current_region->my_class]=current_region->my_contor; 
				current_region->my_region=new int[current_region->my_contor];
				
				memcpy(current_region->my_region,conn_selected,sizeof(int)*current_region->my_contor);    
				prev_region=current_region;
		}  
		current_region->next_region_str=0;
		
		delete [] ptr_tab; delete [] taken; delete [] conn_selected;
		return first_region;
}  

// Find connected components and remove small regions of classes 
// with small regions 
void SegmenterMS::conn_comp(Octet *my_class, int *n_rects, Octet *inv_map, float T[][p_max],int my_lim, int change_type)
{
	REGION *first_region; 
	int *my_max_region;
	if(change_type==0)
    {
		my_max_region = new int[(*n_rects)+1];
		memset(my_max_region,0,sizeof(int)*((*n_rects)+1));
    }
	first_region=create_region_list(my_max_region, change_type);
	if(change_type==0)  //elliminate classes with small regions
		eliminate_class(my_class,my_max_region,n_rects,my_lim,inv_map,T,first_region); 
	else if(change_type==1)  //elliminate small regions
		eliminate_region(n_rects,my_lim,T,first_region);
	destroy_region_list(first_region);
	if(change_type==0) delete [] my_max_region;
	cerr<<":";
}


// Cut a rectangle from the entire input data
// Deletes the previous rectangle, if any
void SegmenterMS::cut_rectangle( sRectangle* rect )
{
	if ( _data ) {
		for ( register int i = 0; i < _p; i++ )
			if ( _data[i] )   delete [] _data[i];
			delete [] _data;
	}
	
	// Set the dimensions of the currently processed region.
	_rrows  = rect->height;
	_rcolms = rect->width;
	_data = new int*[_p];
	
	register int my_x = rect->x;
	register int my_y = rect->y;
	register int i, j, d;
	for ( i = 0; i < _p; i++ )
		_data[i] = new int[_rcolms*_rrows];
	
	if(auto_segm)
		for ( d = 0; d < _p; d++ )
			memcpy(_data[d], _data_all[d],sizeof(int)*_ro_col);
		else
		{ 
			int	idx1 = my_y * _colms + my_x;
			int	idx2 = 0;
			for ( j = my_y, d;
            j < my_y + _rrows; j++, idx1 += _colms - _rcolms )
				for ( i = my_x; i < my_x + _rcolms; i++, idx1++, idx2++ )
				{
					for ( d = 0; d < _p; d++ )
						_data[d][idx2] = _data_all[d][idx1];
				}
		}
		cerr<<":";
}

// Compute the mean of N points given by J[]
void SegmenterMS::mean_s(const int N, const int p, int J[], int **data, float T[])
{
	int TT[p_max];
	register int k, i, j;
	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;by dzj;
	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 //auto_segm;
			{
				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++;
				}
			}//end of if(auto_segm)
			//  Objective functions
			Obj_fct[l]=(float) c;
			l++;
		}else //subsample(Xmean);
		{
			++ll;
		} //end of if ( subsample(Xmean) )
    }//end while;
	
	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 SegmenterMS::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]=(float) ( 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]=(float) ( rect_RADIUS[i]*sqrt(my_th/100) );
    }
	final_RADIUS=(float) ( 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];//_p为维数=3;
	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;//break for, next rect;
			}
		}
		
		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 //RAD2<0.1;
					{
						T[rect][0] = final_T[0];
						T[rect][1] = final_T[1];
						T[rect][2] = final_T[2];
					} //end of if RAD2<0.1;
				}else //auto_segm;
				{
					my_contor = 0;
					memset(my_class,0,_ro_col);