rainfalling.cpp

刚上传内容的相关CODEC不能单独上传。于是

// Rainfalling.cpp: implementation of the CRainfalling class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Rainfalling.h"

#include "IEIntImage.h"
#include "IEFloatImage.h"

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

//////////////////////////////////////////////////////////////////////
// Global Varibles
//////////////////////////////////////////////////////////////////////

static float * south;
static float * north;
static float * south_east;
static float * south_west;
static float * north_east;
static float * north_west;

static float * south_org;
static float * north_org;
static float * south_east_org;
static float * south_west_org;
static float * north_east_org;
static float * north_west_org;

//////////////////////////////////////////////////////////////////////
// Steepest 8 Macros
//////////////////////////////////////////////////////////////////////

#define SINGLE_PIXEL_SEGMENT -1

#define SWITCH_PTRS_8() \
	temp_ptr = north_org; \
	north = north_org = south_org; south = south_org = temp_ptr; \
	temp_ptr = north_east_org; \
	north_east = north_east_org = south_west_org; \
	south_west = south_west_org = temp_ptr; \
	temp_ptr = north_west_org; \
	north_west = north_west_org = south_east_org; \
	south_east = south_east_org = temp_ptr;

// Merge only if the central pixel is NOT the local minimum
#define MERGE(BR, BD, A, B) \
	if( minimum < 0.0f ) { \
		label_cp = out[BR ][BD ]; \
		label_sd = out[BR+A][BD+B]; \
		MergeSegments(label_cp, label_sd); \
	}

// 
#define DIFF_MIN_CMP(VAL, K, L) \
	if(VAL < minimum) { \
		minimum = VAL; r = K; c = L; \
	}

// Also inits 'minimum' and 'r' and 'c', so it
// should always be the first one used !!!!
#define SOUTH() *south++ = minimum = sp - cp; r = 1; c = 0;

// Also inits 'minimum' and 'r' and 'c', so it
// should always be the first one used !!!!
#define F_SOUTH() *south++ = sp - cp;

// Also inits 'minimum' and 'r' and 'c', so it
// should always be the first one used !!!!
#define SOUTH_RB() *south++ = minimum = sep - ep; r = 1; c = 0;

#define SOUTH_EAST() \
	*south_east++ = difference = (sep - cp)/1.4142136f; \
	DIFF_MIN_CMP(difference, 1, 1);
	
#define F_SOUTH_EAST() *south_east++ = (sep - cp)/1.4142136f;

#define SOUTH_WEST() \
	*south_west++ = difference = ( swp - cp )/1.4142136f; \
	DIFF_MIN_CMP(difference, 1, -1);

#define F_SOUTH_WEST() *south_west++ = (swp - cp)/1.4142136f;

#define SOUTH_WEST_RB() \
	*south_west++ = difference = ( sp - ep )/1.4142136f; \
	DIFF_MIN_CMP(difference, 1, -1);

#define EAST() \
	east = difference = ep - cp ; \
	DIFF_MIN_CMP(difference, 0, 1);

#define F_EAST() east = ep - cp; 	

// should only be used before new east calculation
#define WEST()          DIFF_MIN_CMP( -east, 0, -1) ;  

#define NORTH()         DIFF_MIN_CMP( -(*north), -1, 0); north++;  

#define F_NORTH()       north++;  

#define NORTH_WEST()    DIFF_MIN_CMP( -(*north_west), -1, -1); north_west++;

#define F_NORTH_WEST()  north_west++;

#define NORTH_EAST()    DIFF_MIN_CMP( -(*north_east), -1, +1); north_east++;

#define F_NORTH_EAST()  north_east++;


#define STEEPEST8_UPPER_LEFT_CORNER() \
	SOUTH(); \
	EAST(); \
	SOUTH_EAST(); \
	MERGE(0, 0, r, c); 

#define STEEPEST8_1st_LINE_INTERNAL() \
	SOUTH(); \
	SOUTH_EAST(); \
	SOUTH_WEST(); \
	WEST(); \
	EAST(); \
	MERGE(0, l, r, c); 

#define STEEPEST8_UPPER_RIGHT_CORNER() \
	SOUTH_RB(); \
	WEST(); \
	SOUTH_WEST_RB(); \
	MERGE(0, w - 1, r, c); 

#define STEEPEST8_LEFT_BORDER_PIXEL() \
	SOUTH(); \
	SOUTH_EAST(); \
	NORTH(); \
	NORTH_EAST(); \
	EAST(); \
	MERGE(k, 0, r, c); 

#define STEEPEST8_INTERNAL_PIXEL() \
	SOUTH(); \
	SOUTH_WEST(); \
	SOUTH_EAST(); \
	WEST(); \
	EAST(); \
	NORTH(); \
	NORTH_EAST(); \
	NORTH_WEST(); \
	MERGE(k, l, r, c); 

#define F_STEEPEST8_INTERNAL_PIXEL() \
	F_SOUTH(); \
	F_SOUTH_WEST(); \
	F_SOUTH_EAST(); \
	F_EAST(); \
	F_NORTH(); \
	F_NORTH_EAST(); \
	F_NORTH_WEST();

#define STEEPEST8_RIGHT_BORDER_PIXEL() \
	SOUTH_RB(); \
	WEST(); \
	SOUTH_WEST_RB(); \
	NORTH(); \
	NORTH_WEST(); \
	MERGE(k, w - 1, r, c); 

#define STEEPEST8_LAST_LINE() \
	difference = minimum = ep - cp ; r = 0; c = 1; \
	WEST(); \
	east = difference; \
	NORTH(); \
	NORTH_EAST(); \
	NORTH_WEST(); \
	MERGE(h - 1, l, r, c); 

#define NORTH_LL() \
	minimum = -(*north); r = -1; c = 0; north++;

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CRainfalling::CRainfalling()
{
	m_pResult = m_pBuffer = NULL;
	m_nSegments = 0;
}

CRainfalling::~CRainfalling()
{
	if( m_pResult )
		delete m_pResult;

	if( m_pBuffer )
		delete m_pBuffer;
}

//////////////////////////////////////////////////////////////////////
// Member Functions
//////////////////////////////////////////////////////////////////////

// Watershed of a (gradient type) image in 8 connectivity   
// see the relevant papers/reports for a full explanation on how
// this routine works; some comments:
// major parts (3 in total):
// 111111111111111111111111  : "upper row" part
// 222222222222222222222222  : "internal image" part starts here
// 222222222222222222222222  
// 222222222222222222222222  
// 222222222222222222222222  
// 222222222222222222222222  
// 222222222222222222222222  
// 222222222222222222222222  
// 222222222222222222222222  : "internal image" part stops here
// 333333333333333333333333  : "lowest row" part

// within each part there is an additional distinction between
// the first (leftmost), the "internal" and the last (rightmost) 
// pixel in each line of the image being processed

void CRainfalling::Watershed(IEFloatImage* pInput, float drowning_level)
{
	// auxiliary variables
	int k, l; 
	int r, c;
	int label_cp, label_sd;

	int w = pInput->GetWidth();
	int h = pInput->GetHeight();
	
	float * temp_ptr;
	
	register float difference, minimum;
	register float cp, ep, sp, sep, swp;
	register float east;
	
	// Output image creation
	AllocBuffer(w, h);
	float** in = pInput->GetPixels2D();
	int** out  = m_pResult->GetPixels2D();

	int alloc_space = w + 1; // +1 to be safe for '*ptr++'
	south = south_org = new float[sizeof(float) * alloc_space];
	north = north_org = new float[sizeof(float) * alloc_space];
	south_east = south_east_org = new float[sizeof(float) * alloc_space];
	north_east = north_east_org = new float[sizeof(float) * alloc_space];
	south_west = south_west_org = new float[sizeof(float) * alloc_space];
	north_west = north_west_org = new float[sizeof(float) * alloc_space];

	// ############## process the upper row ##############
	// ---------upper left corner---------
	cp = in[ 0 ][ 0 ] ; ep  = in[ 0 ][ 1 ]; 
	sp = in[ 1 ][ 0 ] ; sep = in[ 1 ][ 1 ];
	if( cp <= drowning_level )
	{
		if( sp <= drowning_level ) 
			MergeSegments(out[ 0 ][ 0 ], out[ 1 ][ 0 ]);
		if( ep <= drowning_level ) 
			MergeSegments(out[ 0 ][ 0 ], out[ 0 ][ 1 ]);
		if( sep <= drowning_level ) 
			MergeSegments(out[ 0 ][ 0 ], out[ 1 ][ 1 ]);
	}
	STEEPEST8_UPPER_LEFT_CORNER() ;

	// ----------internal line ----------- 
	for( l = 1 ; l < w - 1 ; l++ )
	{
		cp = ep; swp = sp; sp = sep; 
		ep = in[ 0 ][l+1]; sep = in[ 1 ][l+1];
		if( cp <= drowning_level )  
		{
			if( sp <= drowning_level ) 
				MergeSegments(out[ 0 ][ l ], out[ 1 ][ l ]);
			if( ep <= drowning_level ) 
				MergeSegments(out[ 0 ][ l ], out[ 0 ][l+1]);
			if( sep <= drowning_level ) 
				MergeSegments(out[ 0 ][ l ], out[ 1 ][l+1]);
			if( swp <= drowning_level ) 
				MergeSegments(out[ 0 ][ l ], out[ 1 ][l-1]);
		}
		STEEPEST8_1st_LINE_INTERNAL() ;
	}

	// ----------upper right corner---------
	// l = w - 1
	if( ep <= drowning_level )
	{
		if( sp <= drowning_level )
			MergeSegments(out[ 0 ][ l ], out[ 1 ][l-1]);
		if( sep <= drowning_level )
			MergeSegments(out[ 0 ][ l ], out[ 1 ][ l ]);
	}
	STEEPEST8_UPPER_RIGHT_CORNER();
	SWITCH_PTRS_8();
	
	// ################### general case; internal image ################
	for( k = 1 ; k < h - 1 ; k++ )
	{
		// -------left border pixel-------
		cp  = in[ k ][ 0 ]; ep = in[ k ][ 1 ]; 
		sep = in[k+1][ 1 ]; sp = in[k+1][ 0 ] ;
		if( cp <= drowning_level )
		{
			if( sp <= drowning_level )
				MergeSegments(out[ k ][ 0 ],out[k+1][ 0 ]);
			if( ep <= drowning_level )
				MergeSegments(out[ k ][ 0 ],out[ k ][ 1 ]);
			if( sep <= drowning_level )
				MergeSegments(out[ k ][ 0 ],out[k+1][ 1 ]);
		}
		STEEPEST8_LEFT_BORDER_PIXEL();
		
		// --------internal line---------- 
		for( l = 1 ; l < w - 1 ; l++ )
		{
			cp = ep; swp = sp; sp = sep ; // switch 
			ep = in[ k ][l+1]; sep = in[k+1][l+1];
			if( cp <= drowning_level )  
			{
				if( sp <= drowning_level )
					MergeSegments(out[ k ][ l ], out[k+1][ l ]);
				if( swp <= drowning_level )
					MergeSegments(out[ k ][ l ], out[k+1][l-1]);

label_ep_red :  // jump label

				if( sep <= drowning_level )
					FastMergeSegments(out[ k ][ l ], out[k+1][l+1]);

				if( ep <= drowning_level ) 
				{
					MergeSegments(out[ k ][ l ], out[ k ][l+1]);
					STEEPEST8_INTERNAL_PIXEL();
					if( l >= w - 2 )	continue;
					l++;
					cp = ep; swp = sp; sp = sep ; // switch
					ep = in[ k ][l+1]; sep= in[k+1][l+1];
					goto label_ep_red;
				}
				else
				{
					if( l < w - 2 )
					{
						STEEPEST8_INTERNAL_PIXEL();
						l++;
						cp = ep; swp = sp; sp = sep; // switch
						ep = in[ k ][l+1]; sep= in[k+1][l+1];
					}
				}
			}
			STEEPEST8_INTERNAL_PIXEL() ;
		} 

		// --------right border pixel-------- 
		l = w - 1;
		if( ep <= drowning_level ) 
		{
			if( sp <= drowning_level )
				MergeSegments(out[ k ][ l ], out[k+1][l-1]);
			if( sep <= drowning_level )
				MergeSegments(out[ k ][ l ], out[k+1][ l ]);
		}
		STEEPEST8_RIGHT_BORDER_PIXEL();
		SWITCH_PTRS_8();
	} 

	// ############### process the lowest row ############## 
	// ---------lower left corner--------- 
	k = h - 1;
	cp = in[ k ][ 0 ];  ep = in[ k ][ 1 ]; 
	if( ( cp <= drowning_level ) && ( ep <= drowning_level ) )
		MergeSegments(out[ k ][ 0 ], out[ k ][ 1 ]) ;

	NORTH_LL();
	EAST();
	NORTH_EAST();

	// ---------internal  lower line---------
	for( l = 1 ; l < w - 2 ; l++ ) 
	{
		cp = ep; ep = in[ k ][l+1]; 
		if( ( cp <= drowning_level ) && ( ep <= drowning_level ) )
			MergeSegments(out[ k ][ l ], out[ k ][l+1]);

		STEEPEST8_LAST_LINE();
	}

	// ---------lower right corner---------
	NORTH_LL();
	WEST();
	NORTH_WEST();
	
	delete [] south_org;
	delete [] north_org;
	delete [] south_east_org;
	delete [] north_east_org;
	delete [] south_west_org;
	delete [] north_west_org;

	//
	/*
	int count = 0;
	int size = w * h;
	int* segnr = m_pBuffer->GetPixels();
	for( l = 0 ; l < size ; l++ )
	{
		if( segnr[l] != SINGLE_PIXEL_SEGMENT )
			count++;
	}
	m_nSegments = size - count;
	*/
}

void CRainfalling::AllocBuffer(int w, int h)
{
	// Memory allocation
	if( m_pResult )
		delete m_pResult;
	m_pResult = new IEIntImage(w, h);

	if( m_pBuffer )
		delete m_pBuffer;
	m_pBuffer = new IEIntImage(w, h);
	
	// Initialization
	register int i;
	int size = m_pBuffer->GetSize();

	int* pixels = m_pResult->GetPixels();
	int* segnr = m_pBuffer->GetPixels();
	for( i = 0 ; i < size ; i++ )
	{
		pixels[i] = i;
		segnr[i] = SINGLE_PIXEL_SEGMENT;
	}

	m_nSegments = size;
}

// merge the two segments with label first and next
void CRainfalling::MergeSegments(int first, int next)
{
	int offset;

	// these names make it easier to understand the algorithm
	// but read them as smallest SEGMENT, largest SEGMENT
	int smallest, largest;
	
	// if they already belong to/are the same segment
	if( first == next )
		return; // ALREADY_MERGED; //do nothing
	
	// normally the bigger identifier is also a smaller segment
	// because of the way they were initialized (*) and the way
	// the segments are progressively grown/merged (also (*) in
	// video scanning order)
	if( first <= next )
	{
		smallest = next ; largest = first;
		//ReturnValue = MERGED_IN_SWITCHED_ORDER;
	}
	else 
	{
		smallest = first; largest = next;
		//ReturnValue = MERGED_IN_ORDER;
	}

	//
	int* pixels = m_pResult->GetPixels();
	int* segnr = m_pBuffer->GetPixels();

	// search end of smallest SLL (singly linked list) and 
	// relabel all pixels belonging to it while doing so.
	offset = smallest;
	while( segnr[offset] != SINGLE_PIXEL_SEGMENT )  
	{
		pixels[offset] = largest;    //relabel
		offset = segnr[offset];      //next in SLL
	}
	pixels[offset] = largest;        //last relabel
	
	// prefix smallest SLL to largest SLL
	segnr[offset ] = segnr[largest];
	segnr[largest] = smallest;
	
	/* 
	What happens is:
	a--->b--->c--->X     D--->E--->Y      
	(X and Y end pointers of lists a,b,c and D,E)
	(so X and Y are both == SINGLE_PIXEL_SEGMENT)
	changes to:
	   __________________
	  /                  \
	 /                    \/
	a    b--->c--->X       D--->E--->Y      
	     /\                         /
	       \_______________________/
	*/
	
	// update the number of remaining segments
	m_nSegments--;
}

void CRainfalling::FastMergeSegments(int first, int next)
{
	// these names make it easier to understand the algorithm
	// but read them as smallest SEGMENT, largest SEGMENT
	int smallest, largest;

	// if they already belong to/are the same segment
	if( first == next )
		return; // ALREADY_MERGED; //do nothing
	
	// normally the bigger identifier is also a smaller segment
	// because of the way they were initialized (*) and the way
	// the segments are progressively grown/merged (also (*) in
	// video scanning order)
	if( first <= next )
	{
		smallest = next ; largest = first;
		//ReturnValue = MERGED_IN_SWITCHED_ORDER;
	}
	else 
	{
		smallest = first; largest = next;
		//ReturnValue = MERGED_IN_ORDER;
	}

	int* pixels = m_pResult->GetPixels();
	int* segnr = m_pBuffer->GetPixels();

	pixels[smallest] = largest;     // last relabel

	segnr[smallest] = segnr[largest];
	segnr[largest ] = smallest;

	/* 
	What happens is:
	a--->b--->c--->X#     Y#      
	(X and Y end pointers of lists a,b,c and Y)
	(X and Y are both == SINGLE_PIXEL_SEGMENT )
	changes to:
	   __________________
	  /                  \
	 /                    \/
	a    b--->c--->X#      Y#
	     /\              /
	       \____________/
	*/

	//update the number of remaining segments
	m_nSegments--;
}