completelinkage.cpp

这是一个GPS相关的程序

{
	int xmin = (x > 0) ? x-1 : 0;
	int xmax = (x < (Get_NX()-1) ) ? x+1 : Get_NX()-1;
	int ymin = (y > 0) ? y-1 : 0;
	int ymax = (y < (Get_NY()-1) ) ? y+1 : Get_NY()-1;
	double dist, minDist = DBL_MAX, minClass;
	int dir; 
	long class1, class2;
	class1 = pClassGrid->asInt(x,y);
	if( !pClassGrid->is_NoData_Value(class1) ){
		for( int lx = xmin; lx <= xmax; lx++ ){
			for( int ly = ymin; ly <= ymax; ly++ ){
				class2 = pClassGrid->asInt(lx,ly);
				if( (class2 != class1) && !pClassGrid->is_NoData_Value(class2) ){
					//dist = minDist-1;
					try{
						dist = getClassesDist( class1, class2, cache, remClass1, remClass2 );
						if( dist < minDist ){
							minDist = dist;
							dir = (lx+1-x) +3*(ly+1-y);
							minClass = class2;
						}
					}catch( CNoDataValueError ){ 
						minDist = minDist;
						/* just do not account direction */ 
					}
				}// class1 != class2
			}//x
		}//y
	}// no data value
	if( minDist == DBL_MAX ){
		minDist = 0;
		dir = 4;
		minClass = class1;
	}
	double *ret = new double[3];
	ret[0] = minDist;
	ret[1] = dir;
	ret[2] = minClass;
	return ret;
}


// returns the distance of both classes
// uses and fills a cache (classDistMap), if cache fails it calls calcClassesDist
// see clearClassDistMap, for clearing cahce
// remClass1, remClass2 : Classes to remove from Distance-Cache
double CCompleteLinkage::getClassesDist( long class1, long class2, distMapT& cache, long remClass1, long remClass2 ){
	// same class has no distance 
	if(  class1 == class2 ) return 0;
	if( pClassGrid->is_NoData_Value(class1) || pClassGrid->is_NoData_Value(class2) ) 
		throw CNoDataValueError("getClassesDist with one class as NoData");
	// distances are symetric, so evt switch classes: class2 > class1
	if( class1 > class2 ){
		long tmp = class1;
		class1 = class2;
		class2 = tmp;
	}
	if( class2 < Get_NCells() ){
		prof.start("getClassesDist_singlePixel");
		// both classes are just single pixels
		double x1,y1,x2,y2;
		x1 = class1 % Get_NX(); y1 = class1 / Get_NX();
		x2 = class2 % Get_NX(); y2 = class2 / Get_NX();
		// NoDataValueError must be caught at higher level, and insert NoDataValue in ClassesGrid
		return calculateEucledianDistance(x1,y1,x2,y2);
		prof.stop("getClassesDist_singlePixel");
	}
	// cash of calculated class differences
	// mapping class2(long) -> class1(long) -> distance(double)
	double dist;
	distPtrT pDistMap;
	distMapT::iterator it = cache.find( class2 );
	distT::iterator it2;
	if( it == cache.end() ){
		pDistMap = new distT();
		cache.insert( distMapT::value_type( class2, pDistMap) );
	}else{
		pDistMap = (*it).second;
		// remove former classes from distance cache
		if( remClass1 >= Get_NCells() ){
			it2 = pDistMap->find( remClass1 );
			if( it2 != pDistMap->end() ) 
				pDistMap->erase( it2 );
		}
		if( remClass2 >= Get_NCells() ){
			it2 = pDistMap->find( remClass2 );
			if( it2 != pDistMap->end() ) 
				pDistMap->erase( it2 );
		}
	}
	it2 = pDistMap->find( class1 );
	if( it2 == pDistMap->end() ){
		prof.start("getClassesDist_calc");
		try{
			dist = calcClassesDist(class1, class2);
			pDistMap->insert( distT::value_type( class1, dist ) );
		}catch( CNoDataValueError ){
			// should not occure
			throw runtime_error("cauth noDataError of calcClassesDist of valid classes");
		}
		prof.stop("getClassesDist_calc");
	}else{
		prof.start("getClassesDist_cache");
		dist = (*it2).second;
		prof.stop("getClassesDist_cache");
	}
	return dist;
}

// calulates distance betwenn two classes
// this is the maximum the distances between each pixel of class1 
// related to each pixel of class2
double CCompleteLinkage::calcClassesDist(long class1, long class2){
	int x1,y1,x2,y2, pi;
	double dist, maxDist = -1;

	pixelSetPtrT pC1pix = getPixelSetItForClass(class1,1).first;
	pixelSetPtrT pC2pix = getPixelSetItForClass(class2,2).first;

	for( pixelSetT::iterator it1 = pC1pix->begin(); it1 != pC1pix->end(); ++it1 ){
		pi = (*it1);
		x1 = (*it1) % Get_NX();
		y1 = (*it1) / Get_NX();
		for( pixelSetT::iterator it2 = pC2pix->begin(); it2 != pC2pix->end(); ++it2 ){
			pi = (*it2);
			x2 = (*it2) % Get_NX();
			y2 = (*it2) / Get_NX();
			try{			
				dist = calculateEucledianDistance(x1,y1,x2,y2);
				if( dist > maxDist ) maxDist = dist;
			}catch( CNoDataValueError ){} // do not account for distances to missing pixels
		}
	}

	/*
	// outer loop: all pixels of class1
	for(y1=0; y1<Get_NY(); y1++)	{
		for(x1=0; x1<Get_NX(); x1++)		{
			if( pClassGrid->asInt(x1,y1) == class1 ){

				// inner loop: all pixels of class2
				for(y2=0; y2<Get_NY(); y2++)			{
					for(x2=0; x2<Get_NX(); x2++)				{
						if( pClassGrid->asInt(x2,y2) == class2 ){
							dist = calculateEucledianDistance(x1,y1,x2,y2);
							if( dist > maxDist ) maxDist = dist;
						}// if class2 inner loop
					}// x inner loop
				}// y inner loop

			}// if class1 outer loop
		}// x outer loop
	}//y outer loop
	*/

	return maxDist;
}


// returns the Distance in parameter-space for given two pixel
double CCompleteLinkage::calculateEucledianDistance(int x1, int y1, int x2, int y2){
	double c1,c2,dist, sumDist = 0;
	double w;
	for( int i = 0; i < MAX_INPUT_GRIDS; i++ ){
		if( pInputGrid[i] != NULL){
			w = inputWeight[i];
			c1 = pInputGrid[i]->asDouble(x1,y1);
			c2 = pInputGrid[i]->asDouble(x2,y2);
			if( pInputGrid[i]->is_NoData_Value(c1) ) throw CNoDataValueError("calculateEucledianDistance");
			if( pInputGrid[i]->is_NoData_Value(c2) ) throw CNoDataValueError("calculateEucledianDistance");
			dist =  c1 - c2;
			dist = w * dist * dist; //weight and square
			sumDist += dist;
		}
	}
	//sqare root not needed, because its monotone
	return sumDist; 
}

// copies on grid to another grid, targetGrid is allocated
CSG_Grid *CCompleteLinkage::copyGrid(CSG_Grid *fromGrid, TSG_Grid_Type gridType){
	int x,y;
	CSG_Grid *targetGrid;
	targetGrid = SG_Create_Grid(gridType, Get_NX(), Get_NY(), Get_Cellsize());
	for(y=0; y<Get_NY();  y++){
		for(x=0; x<Get_NX(); x++){
			switch( gridType ) {
				case GRID_TYPE_Byte:
					targetGrid->Set_Value(x,y, fromGrid->asByte(x,y) );
					break;
				case GRID_TYPE_Char:
					targetGrid->Set_Value(x,y, fromGrid->asChar(x,y)  );
					break;
				case GRID_TYPE_Short:
					targetGrid->Set_Value(x,y, fromGrid->asShort(x,y)  );
					break;
				case GRID_TYPE_Int:
				case GRID_TYPE_Word:
				case GRID_TYPE_DWord:
					targetGrid->Set_Value(x,y, fromGrid->asInt(x,y)  );
					break;
				case GRID_TYPE_Float:
					targetGrid->Set_Value(x,y, fromGrid->asFloat(x,y)  );
				/*
					*/
				default:
					targetGrid->Set_Value(x,y, fromGrid->asDouble(x,y) );
					break;
			}// switch
 		}//x
	}//y
	return targetGrid;
}


// returns the class of the pixel with minimal distance
// pMinDistDir must have been calculated
long CCompleteLinkage::getDestClass(int x, int y){
	int dir = (int) pMinDistDir->asDouble(x,y);
	if( dir == pMinDistDir->Get_NoData_Value() ) 
		return pClassGrid->Get_NoData_Value();
	int dx = DIRX(dir);
	int dy = DIRY(dir);
	return pClassGrid->asInt( x+DIRX(dir), y+DIRY(dir) );
}

/*
void CCompleteLinkage::removeOrphantClasses(long minClassCnt, CSG_Grid *pOrphantsGrid){
	// calulate number of pixels per class and minimum distance
	classPixelSetMapT classPixelMap;	//class -> List of Pixel
	classPixelSetMapT::iterator it1;

	classPixelSetMapT::value_type pixelVec;
	classPixelSetMapT::value_type::iterator it2;

	long class1, class2, cnt, pi;
	double edist;
	int x,y;
	for(y=0; y<Get_NY() && Set_Progress(y); y++){
		for(x=0; x<Get_NX(); x++){
			class1 = pClassGrid->asDouble(x,y); 
			it1 = classPixelMap.find(class1);
			if( it1 == classPixelMap.end() ){
				pixelVec.clear();
				pixelVec.insert( pixelVec.end(), (y*Get_NX()+x) );
				//inserts copy, so we do not need to create new pointer
				classPixelMap.insert( classPixelSetMapT::value_type( class1, pixelVec ) );
			}else{
				//(*it1).second  is pixelVec
				(*it1).second.insert( (*it1).second.end(), (y*Get_NX()+x) );
			}
		}
	}

	// reclass orphant pixels and remember in output
	double avgDist;
	double *dist_array;
	for( it1 = classPixelMap.begin(); it1 != classPixelMap.end(); ++it1 ){
		class1 = (*it1).first;
		cnt = (*it1).second.size();
		if( cnt < minClassCnt ){
			// calculate class distances to each neighboring class
			clearClassDistMap();
			for( it2 = (*it1).second.begin(); it2 != (*it1).second.end(); ++it2 ){
				pi = (*it2);
				//x = pi % Get_NX();
				//y = pi / Get_NX();
				dist_array = calcPixelsClassDistance(pi);
				// class distances anyway
			}

		}
	}
}
*/


void CCompleteLinkage::searchSetsMinPixel(pixelSetPtrT pPix, int &minx, int &miny){
	double edist, minDist = DBL_MAX;
	long pi;
	for( pixelSetT::iterator it = pPix->begin(); it != pPix->end(); it++ ){
		pi = (*it);
		edist = pMinDist->asDouble(pi);
		if( (edist > 0) && (edist < minDist)){
			minDist = edist;
			minx = PI2X(pi);
			miny = PI2Y(pi);
		}
	}
}


// remove orphant classes
// all classes, consisting of fewer than minClassCnt pixels are attached to 
// classes with minimal distance
// the distance of these pixels is put to pOrphantsGrid
// pClassGrid, pMinDist, pMinDistDir must have been calculated
void CCompleteLinkage::removeOrphantClasses(
		const long minClassCnt, CSG_Grid *pOrphantsGrid,
		long& i, 
		/* inout */ distMapT& classDistMap,	// mapping class -> ptr(map class -> distance)
		/* out */ RegionSetT& affectedRegionSet // set of region pixels, thats minimum value has been changed
		){
	long class1;
	double edist;
	pixelSetPtrT pPix, pPixN;
	vector<long> littleSets;
	//vector<pixelSetPtrT>::iterator itp;
	int minx, miny;
	classPixelSetMapT::iterator itclass;

	for( int y = 0; y < Get_NY(); y++ ){
		for( int x = 0; x < Get_NX(); x++ ){
			class1 = pClassGrid->asInt(x,y);
			if( !pClassGrid->is_NoData_Value(class1) && class1 < Get_NCells() ){
				edist = pMinDist->asDouble(x,y);
				aggregatePixel(i++, x, y, classDistMap, affectedRegionSet );
				if( pOrphantsGrid != NULL ) 
					pOrphantsGrid->Set_Value(x,y,edist);
			}else if(pOrphantsGrid != NULL ) {
				pOrphantsGrid->Set_NoData(x,y);
			}	
		}
	}

	if( minClassCnt > 1 ){
		//search all Sets and remember litte sets
		for( classPixelSetMapT::iterator it = classPixelMap.begin(); it != classPixelMap.end(); ++it ){
			pPix = (*it).second;
			if( pPix->size() < minClassCnt )
				littleSets.push_back( (*it).first );
		}
		// continue to merge pixelSets until merged Sets are large enough
		while( littleSets.size() > 0 ){
			class1 = littleSets[ littleSets.size()-1 ];
			littleSets.pop_back();
			itclass = classPixelMap.find(class1);
			// do search for class, because it could have been aggregated away in meantime
			if( itclass != classPixelMap.end() ){
				pPix = (*itclass).second;
				searchSetsMinPixel( pPix, minx, miny );
				if( pOrphantsGrid != NULL ){
					edist = pMinDist->asDouble(minx,miny);
					for( pixelSetT::iterator it = pPix->begin(); it != pPix->end(); it++ )
						pOrphantsGrid->Set_Value( (*it), edist );
				}
				pPixN = aggregatePixel(i++, minx,miny, classDistMap, affectedRegionSet ); 
				// if new Class has still to few pixels, remember its class
				if( pPixN->size() < minClassCnt )
					littleSets.push_back( Get_NCells()+i ); 
			}// class yet exists
		}// loop littleSets
	}// minClassCnt > 0

	/* XXXX unite smaller classes */
	/*
	// calulate number of pixels per class and minimum distance
	map<long,long> classCntMap;	//class -> cnt
	map<long,long>::iterator it1;
	map<long,long>::iterator it3;
	
	map<long,long> classMinPixelMap; // class -> pixel, that has minimal distance
	map<long,long>::iterator it2;

	map<long,long> classesToChange; // mapping of class->class to change
	map<long,long>::iterator it4;

	set<long> remainingOrphants; // classes that remain orphants after uniting with neighbor
	set<long> destClasses; // classes that are destination of aggregation
							// this list is needed because, pixels of target Classes must be recalculated too
	set<long>::iterator itS;
							
	
	vector<long> changedPixels; // pixels that have been united
	vector<long>::iterator itV;

	double *dist_dir;

	long class1, class2, cnt, pi;
	double edist;
	int x,y;
	
	// clear oprhants Grid
	//pOrphantsGrid = SG_Create_Grid(GRID_TYPE_Double, Get_NX(), Get_NY(), Get_Cellsize()); input parameter
	if( pOrphantsGrid != NULL ){
		double ndv = pOrphantsGrid->Get_NoData_Value();
		for(y=0; y<Get_NY(); y++){
			for(x=0; x<Get_NX(); x++){
				pOrphantsGrid->Set_Value(x,y,ndv);
			}
		}
	}

	//unite classes, until size of united classes reaches minimum pixel size
	do{

		remainingOrphants.clear();
		classesToChange.clear();
		destClasses.clear();
		classCntMap.clear();
		classMinPixelMap.clear();

		// construct histogram of classes and minimum points
		for(y=0; y<Get_NY() && Set_Progress(y); y++){
			for(x=0; x<Get_NX(); x++){