completelinkage.cpp

这是一个GPS相关的程序

				edist = calculateEucledianDistance( x,y, x+1, y );
				prevDist = edist;
				if( edist < minDist ){
					dir = 5;
					minDist = edist;
				}
			}catch( CNoDataValueError ){ prevDist = DBL_MAX; }
			setInitDist(x,y,minDist,dir, ndvCnt);
		}//x [1..n-1]
		x = Get_NX()-1; {
			dir = 0;
			minDist = nextRowDist[x][0]; 
			edist = nextRowDist[x][1];
			if( edist < minDist ){
				dir = 1;
				minDist = edist;
			}
			edist = prevDist;
			if( edist < minDist ){
				dir = 3;
				minDist = edist;
			}
			setInitDist(x,y,minDist,dir, ndvCnt);
		} // x = Get_NX)() -1
	}
	delete[] nextRowDist;

	prof.stop("initClassesGrid");
 	return ndvCnt;
}

void CCompleteLinkage::setInitDist(int x, int y, double minDist, int dir, long& ndvCnt){
	long pi = XY2PI(x,y);
	if( minDist == DBL_MAX ){
		// could not calculate Class-Distance to any direction, may be itself No_Data in one of the grids
		pClassGrid->Set_NoData(x,y);
		pMinDist->Set_NoData(x,y);
		pMinDistDir->Set_NoData(x,y);
		ndvCnt++;
		incrementRegionsNdvCnt(x,y);
	}else{
		pClassGrid->Set_Value(x,y,pi);
		pMinDist->Set_Value(x,y,minDist);
		pMinDistDir->Set_Value(x,y,dir);
	}
}

// returns the pixelSet for given class
// and the Iterator to class in classPixelMap or end() if constructed from single pixel 
pair<pixelSetPtrT,classPixelSetMapT::iterator> CCompleteLinkage::getPixelSetItForClass(long cl, byte tmpStorageNumber){
	classPixelSetMapT::iterator it;
	pixelSetPtrT pPixSet;
	if( cl < Get_NCells() ){
		it = classPixelMap.end();
		if( !(tmpStorageNumber <= MAX_SINGLE_PIXEL_ARR) || (tmpStorageNumber==0)) 
			throw runtime_error( ToString(tmpStorageNumber).append(" to large tmpStorageNumber") );
		pPixSet = &(singlePixelSetArr[tmpStorageNumber-1]);
		//pixelSetT::iterator itp = pPixSet->begin();
        //(*itp) = cl;
		(*pPixSet->begin()) = cl;
	}else{
		it = classPixelMap.find(cl);
		if( it != classPixelMap.end() ){
			pPixSet = (*it).second;
		}else{
			string s = "no pixel found for class ";
			s.append( ToString(cl) );
			throw runtime_error(s);
		}
	}
	return pair<pixelSetPtrT,classPixelSetMapT::iterator>( pPixSet, it );
}


// recalculates Minimum of all regions, found in set affectedRegionSet
// stores Minimum in 
void CCompleteLinkage::recalcRegionMins(RegionSetT& affectedRegionSet){
	for( RegionSetT::iterator it = affectedRegionSet.begin(); it != affectedRegionSet.end(); /* empty increment within loop */ ){
		long regpi = (*it);
		int regx = REGPI2REGX(regpi);
		int regy = REGPI2REGY(regpi);
		int minx, miny;
		double minDist;
		long minPixel;
		calcRegionsMinPixel(regx,regy,minx,miny,minDist);
		if( minDist == DBL_MAX ){ 
			minDist = pRegionMinDistGrid->Get_NoData_Value(); // whole region without a valid distance
			minPixel = pRegionMinPixelGrid->Get_NoData_Value();
			if( regy > pRegionMinDistGrid->Get_NY() )
				throw runtime_error( "to large index for grid pRegionMinDistGrid" );
		}else{
			minPixel = XY2PI(minx,miny);
		}
		pRegionMinDistGrid->Set_Value(regx,regy,minDist);
		pRegionMinPixelGrid->Set_Value(regx,regy,minPixel);
		// erase and increment iterator -> see http://groups.google.de/groups?hl=de&lr=&threadm=bxCa7.1518%24h94.176381%40paloalto-snr1.gtei.net&rnum=1&prev=/groups%3Fq%3Dstl%2Biterator%2Berase%26hl%3Dde%26lr%3D%26selm%3DbxCa7.1518%2524h94.176381%2540paloalto-snr1.gtei.net%26rnum%3D1
		affectedRegionSet.erase(it++);
	}//changedRegionSet_it
}

long CCompleteLinkage::getGridsMinPixel(){
	double gridMinDist = DBL_MAX, regMinDist;
	int regMinX = -1, regMinY = -1;
	for( int regx = 0; regx < pRegionMinDistGrid->Get_NX(); regx++ ){
		for( int regy = 0; regy < pRegionMinDistGrid->Get_NY(); regy++ ){
			regMinDist = pRegionMinDistGrid->asDouble(regx, regy);
			if( (regMinDist > 0) && (regMinDist < gridMinDist) ){
				gridMinDist = regMinDist;
				regMinX = regx;
				regMinY = regy;
			}
		}
	}
	if( gridMinDist == DBL_MAX ) 
		throw runtime_error("no min-pixel found");
	return getRegionsMinPixel( regMinX, regMinY );
}

long CCompleteLinkage::getRegionsMinPixel(int regx, int regy){
	return pRegionMinPixelGrid->asInt(regx, regy);
}


void CCompleteLinkage::searchNeighbourPixels( const long nClass, RegionSetT& affectedRegionSet, pixelSetT& neighbourPixels ){
		//neighbourPixels.reserve( neighbourPixels.size() + 2*cNpix.size() ); // allocate memory in range of pixels of referenced class
		neighbourPixels.reserve( neighbourPixels.size() + 64 );
		{for( RegionSetT::iterator it = affectedRegionSet.begin(); it != affectedRegionSet.end(); ++it ){
			long regpi = (*it);
			int regx = REGPI2REGX(regpi);
			int regy = REGPI2REGY(regpi);
			int mx = min( Get_NX(), (regx+1)*regionExtent );
			int my = min( Get_NY(), (regy+1)*regionExtent );
			long fclass;
			for( int y = (regy*regionExtent); y < my; y++ ){
				for( int x = (regx*regionExtent); x < mx; x++ ){
					// get pixels, that point to one of these classes
					fclass = getDestClass(x,y);
					if(fclass == nClass){ 
						long pi = XY2PI(x,y);
						neighbourPixels.push_back( pi );
					}
				}//x
			}//y
		}}
}


/* returns the pixel-Sets for given classes
pCNpix may be the same as pC1Pix of pC2pix
the Iterators point to the entry of pC1pix in given classPixelMap
or they point to end(), if the pixelSets are temporary single pixels
*/
void CCompleteLinkage::getPixelSets(
			classPixelSetMapT& classPixelMap, 
			const long class1, const long class2, 
			pixelSetPtrT& pC1pix, classPixelSetMapT::iterator& itc1pix, 
			pixelSetPtrT& pC2pix, classPixelSetMapT::iterator& itc2pix, 
			pixelSetPtrT& pCNpix
			){
	pair<pixelSetPtrT,classPixelSetMapT::iterator> res;
	res = getPixelSetItForClass( class1, 1 );
	pC1pix = res.first; // pixel set
	itc1pix = res.second; // iterator to entry in pixelMap
	/*
	if( class1 >= Get_NCells()) classPixelMap.erase( itc1pix ); // iterator
	if( class2 >= Get_NCells()) classPixelMap.erase( itc2pix ); // iterator
	*/
	res = getPixelSetItForClass( class2, 2 );
	pC2pix = res.first;
	itc2pix = res.second;

	// classes are symetric, so unite smaller to bigger class, 
	// assign new class the number of the larger class
	// create new class if former class consists of just single pixel
	// let pCNpix point to larger vector or create new one
	if( pC1pix->size() >= pC2pix->size() ){
		if( class1 < Get_NCells() ){
			pCNpix = new pixelSetT;
		}else{
			pCNpix = pC1pix;
		}
	}else{
		if( class2 < Get_NCells() ){
			pCNpix = new pixelSetT;
		}else{
			pCNpix = pC2pix;
		}
	}
}

// update map class->pixelSet, append smaller set 2 to larger set
// copy pixels from smaller vector to new Vector, delete other ones
// but other vectors are not erased from classPixelMap
void CCompleteLinkage::updatePixelSet(const long class1, const long class2, pixelSetPtrT pC1pix, pixelSetPtrT pC2pix, pixelSetPtrT pCNpix){
	long pi;
	if( pCNpix != pC1pix ){
		pCNpix->reserve( pCNpix->size() + pC1pix->size() ); // allocate enough memory
		for( pixelSetT::iterator itn1 = pC1pix->begin(); itn1 != pC1pix->end(); ++itn1){
			pi = (*itn1);
			pCNpix->push_back( pi ); // used in conjunction with reserve

		}
		if( class1 >= Get_NCells()){
			delete pC1pix; // original Classes Pixel Vector are not constructed
		}
	}
	if( pCNpix != pC2pix ){
		pCNpix->reserve( pCNpix->size() + pC2pix->size() );
		for( pixelSetT::iterator itn = pC2pix->begin(); itn != pC2pix->end(); ++itn){
			pi = (*itn);
			pCNpix->push_back( pi );
		}
		if( class2 >= Get_NCells()){
			delete pC2pix;
		}
	}
}

void CCompleteLinkage::aggregateClasses(const long nClass, pixelSetPtrT pCNpix, RegionSetT &affectedRegionSet){
	//determine Regions, in which to search f, or pixels to change
	//look at environment
	long pi, regpi;
	int x,y,regx,regy,modx,mody;
	for( pixelSetT::iterator it = pCNpix->begin(); it != pCNpix->end(); ++it ){
		pi = (*it);
		x = PI2X(pi);
		y = PI2Y(pi);
		regx = X2REGX(x);
		regy = Y2REGY(y);
		// do change class of given pixel
		pClassGrid->Set_Value(x,y, nClass );
		regpi = REGXY2REGPI( regx, regy );
		affectedRegionSet.insert( regpi );
		modx = x % regionExtent;
		mody = y % regionExtent;
		if((modx == 0) && (x>0)){ // left edge of region
			regpi = REGXY2REGPI( regx-1, regy );
			affectedRegionSet.insert( regpi );
			if((mody == 0) && (y>0)){ // top edge of region left corner
				regpi = REGXY2REGPI( regx-1, regy-1 );
				affectedRegionSet.insert( regpi );
			}else if((mody == regionExtent-1) && (y < Get_NY()-1)){ // bottom edge of region left corner
				regpi = REGXY2REGPI( regx-1, regy+1 );
				affectedRegionSet.insert( regpi );
			}
		}else if((modx == regionExtent-1) && (x < Get_NX()-1)){ // right edge of region
			regpi = REGXY2REGPI( regx+1, regy );
			affectedRegionSet.insert( regpi );
			if((mody == 0) && (y>0)){
				regpi = REGXY2REGPI( regx+1, regy-1 ); // top edge of region right corner
				affectedRegionSet.insert( regpi );
			}else if((mody == regionExtent-1) && (y < Get_NY()-1)){ // bottom edge of region right corner
				regpi = REGXY2REGPI( regx+1, regy+1 );
				affectedRegionSet.insert( regpi );
			}
		}
		if((mody == 0) && (y>0)){ // top edge of region
			regpi = REGXY2REGPI( regx, regy-1 );
			affectedRegionSet.insert( regpi );
		}else if((mody == regionExtent-1) && (y < Get_NY()-1)){ // bottom edge of region
			regpi = REGXY2REGPI( regx, regy+1 );
			affectedRegionSet.insert( regpi );
		}
	}
}

void CCompleteLinkage::recalcClassDistances(pixelSetPtrT pCNpix, pixelSetT &neighbourPixels, const long remClass1, const long remClass2, distMapT &classDistMap){
	// delete old classes from first mapping
	distMapT::iterator it;
	it = classDistMap.find( remClass1 );
	if( it != classDistMap.end() ){
		classDistMap.erase( it );
	}
	it = classDistMap.find( remClass2 );
	if( it != classDistMap.end() ){
		classDistMap.erase( it );
	}

	// recalculate minimal distances of cells,
	// that originate or point to classes of aggregated cells (stored in neighbourPixels)
	// clearClassDistMap();
	// pixel of new Class
	{for( pixelSetT::iterator it = pCNpix->begin(); it != pCNpix->end(); ++it ){
		long pi = (*it);
		int x = PI2X(pi);
		int y = PI2Y(pi);
		double *dist_dir = calcPixelsClassDistance(x,y, classDistMap, remClass1, remClass2);
		pMinDist->Set_Value(x,y, dist_dir[0] );
		pMinDistDir->Set_Value(x,y, dist_dir[1] );
		delete[] dist_dir;
	}}
	// all pixel that point to pixels of new class
	{for( vector<long>::iterator it = neighbourPixels.begin(); it != neighbourPixels.end(); ++it ){
		long pi = (*it);
		int x = PI2X(pi);
		int y = PI2Y(pi);
		double *dist_dir = calcPixelsClassDistance(x,y, classDistMap, remClass1, remClass2);
		pMinDist->Set_Value(x,y, dist_dir[0] );
		pMinDistDir->Set_Value(x,y, dist_dir[1] );
		delete[] dist_dir;
	}}
}

void CCompleteLinkage::calcRegionsMinPixel(int regx, int regy, int &minx, int &miny, double &minDist){
	int mx = min( Get_NX(), (regx+1)*regionExtent );
	int my = min( Get_NY(), (regy+1)*regionExtent );
	minDist = DBL_MAX;
	minx = -1; miny = -1;
	for( int y = (regy*regionExtent); y < my; y++ ){
		for( int x = (regx*regionExtent); x < mx; x++ ){
			double edist = pMinDist->asDouble(x,y);
			if(  (edist > 0) && (edist < minDist) ){
				minDist = edist;
				minx = x;
				miny = y;
			}//minDist
		}//x
	}//y
}

// basic step of the algorithm, aggregates specified pixel to its 
// closest neighbor and recalculates all class distances and caches
pixelSetPtrT CCompleteLinkage::aggregatePixel(int i, int minx, int miny,
		/* inout */ distMapT& classDistMap,	// mapping class -> ptr(map class -> distance)
		/* out */ RegionSetT& affectedRegionSet // set of region pixels, thats minimum value has been changed
		){
	long nClass, class1, class2; // new and former class-numbers
	pixelSetT neighbourPixels;	// set of pixels pointing to changed class

	class1 = pClassGrid->asInt(minx,miny);
	class2 = getDestClass(minx,miny);
	nClass = Get_NCells() + i;

	// ------ get PixelSets --------
	//prof.start("searchPixelSets");
	pixelSetPtrT pC1pix, pC2pix; // vector of pixels of old classes 
	pixelSetPtrT pCNpix; // vector of pixels of new class
	classPixelSetMapT::iterator itc1pix, itc2pix;
	getPixelSets( classPixelMap, class1, class2, pC1pix, itc1pix, pC2pix, itc2pix, pCNpix );
	//prof.stop("searchPixelSets");

	// ----- copy pixels -----
	//into larger vector, delete other one and update classPixelMap
	//prof.start("updatePixelSets");
	updatePixelSet( class1, class2, pC1pix, pC2pix, pCNpix );
	// erase old pixelSets from classPixelMap
	if( itc1pix != classPixelMap.end() ) classPixelMap.erase( itc1pix );
	if( itc2pix != classPixelMap.end() ) classPixelMap.erase( itc2pix );
	// insert new pixel-Set
	pair<classPixelSetMapT::iterator, bool> ret =
		classPixelMap.insert( classPixelSetMapT::value_type( nClass, pCNpix ) );
	if( !ret.second )
		throw runtime_error( ToString(nClass).append(" Class: could not insert in classPixelMap") );
	//prof.stop("updatePixelSets");

	// ------ aggregate -----------
	// all pixels of two classes to new class
	// remember regions, that are affected
	prof.start("aggregatePixels");
	affectedRegionSet.clear(); // will be output from next method
	aggregateClasses( nClass, pCNpix, affectedRegionSet ); 
	prof.stop("aggregatePixels");

	//------- search neighbours: --------
	//iterate over all pixels in found regions, remember pixels, that point to new class for recalculating distances
	prof.start("searchNeighbours");
	searchNeighbourPixels( nClass, affectedRegionSet, neighbourPixels );
	prof.stop("searchNeighbours");

	// -------- recalculate class distances ----
	//, update class distance cache		
	prof.start("recalcDistanceAll");
	recalcClassDistances( pCNpix, neighbourPixels, class1, class2, classDistMap );
	prof.stop("recalcDistanceAll");

	return pCNpix;
}


void CCompleteLinkage::checkOutput(long step){
	for( int i = 0; i < MAX_OUTPUT_GRIDS; i++ ){
		if( (pOutputGridAdd[i] != NULL) && (step == outputStep[i])){
			pOutputGridAdd[i]->Assign( pClassGrid );
			CSG_Colors colors;
			DataObject_Get_Colors( pOutputGridAdd[i], colors );
			colors.Random();
			DataObject_Set_Colors( pOutputGridAdd[i], colors );
		}
	}
}


long CCompleteLinkage::incrementRegionsNdvCnt(int x, int y){
	int regx, regy;
	regx = X2REGX(x);
	regy = Y2REGY(y);
	long cnt = pRegionNdvCntGrid->asInt(regx,regy);
	pRegionNdvCntGrid->Set_Value(regx,regy, ++cnt );
	return cnt;
}