segmenterms.cpp

这是一个分水岭程序

// SegmenterMS.cpp: implementation of the SegmenterMS class.
// from《Robust Analysis of Feature Spaces》
// author：Dorin Comaniciu & Peter Meer
// Source is available from
// http://www.caip.rutgers.edu/~meer/RIUL/uploads.html
// which is an alpha version of the code and works with ppm images only.
// I edited it a little, so that it can now run in VC.
// and a single independent class, and support some other format, etc.
// by dzj, 04.07.02
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "SegmenterMS.h"
#include "include\imageobject.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
/*

  SegmenterMS::SegmenterMS()
  {
  
	}
	
	  SegmenterMS::~SegmenterMS()
	  {
	  
		}
*/

//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 <memory.h>
#include <time.h>

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 = (float) 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] = { { (float)3.2405, (float)-1.5371, (float)-0.4985 },
{(float)-0.9693,  (float)1.8760,  (float)0.0416 },
{ (float)0.0556, (float)-0.2040,  (float)1.0573 } };

// Constants for LUV transformation 
static const float     Xn = (float)0.9505;
static const float     Yn = (float)1.0;
static const float     Zn = (float)1.0888;
static const float     Un_prime = (float)0.1978;
static const float     Vn_prime = (float)0.4683;
static const float     Lt = (float)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;

/*
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;
	result_contour = NULL;//by dzj;
}

// 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;
	}
	if (result_contour!=NULL)
	{
		delete [] result_contour; result_contour = NULL;
	}
}

// LUV (final_T[]) to RGB (TI[]) conversion
void SegmenterMS::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 = (float) ( Yn * final_T[0] / 903.3 );
			else
				my_y = (float) ( Yn * pow((final_T[0] + 16.0) / 116.0, 3) );
			
			float u_prime = (float) ( final_T[1] / (13 * final_T[0]) + Un_prime ); 
			float v_prime = (float) ( 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]= (float) ( 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 = (float) ( y / 2550000.0 ); //Yn==1
		if ( tval >  Lt)  l_star = my_pow[(int)(tval*255+0.5)];
        else  l_star = (float) ( 903.3 * tval );
		
        my_temp = x + 15 * y + 3 * z;
		if(my_temp)
		{
			u_prime = (float) ( (float)(x << 2) / (float)(my_temp) );
			v_prime = (float) ( (float)(9 * y) / (float)(my_temp) );
		}
		else
		{
			u_prime = (float) 4.0;
			v_prime = (float) (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];