completelinkage.cpp

这是一个GPS相关的程序

				class1 = pClassGrid->asInt(x,y); 
				edist = pMinDist->asDouble(x,y);
				it1 = classCntMap.find(class1);
				if( it1 == classCntMap.end() ){
					classCntMap.insert( map<long,long>::value_type( class1, 1 ) );
				}else{
					(*it1).second++;
				}
				it2 = classMinPixelMap.find(class1);
				if( it2 == classMinPixelMap.end() ){
					classMinPixelMap.insert( map<long,long>::value_type( class1, (y*Get_NY()+x) ) );
				}else{
					if( (edist > 0) && (edist < pMinDist->asDouble((*it2).second) )) 
						(*it2).second = (y*Get_NY()+x);
				}
			}
		}


		// construct request of class mappings
		for( it1 = classCntMap.begin(); it1 != classCntMap.end(); ++it1 ){
			class1 = (*it1).first;
			cnt =  (*it1).second;
			if( cnt < minClassCnt ){
				it2 = classMinPixelMap.find( class1 );
				pi = (*it2).second;
				class2 = getDestClass(pi);
				if( class2 == class1 ){
					// need to recalculate minPixel
					class2 = class2;
				}
				// check, if destination class is going to change
				it4 = classesToChange.find( class2 );
				if( it4 != classesToChange.end() ){
					class2 = (*it4).second;
				}
				// do nothing on cyclic reference (dest class changed to this class)
				it3 = classCntMap.find( class2 );
				if( class1 != class2 ){ 
					// add new lass change request
					classesToChange.insert( map<long,long>::value_type( class1, class2 ));
					destClasses.insert( class2 );
					// increment PixelCount of destination class
					(*it3).second += cnt;
				} // class1 != class2
				// check if target classes pixel count is still smaller than mininum 
				// occurs only on cyclic references: class1 -> class2 -> class1
				if( (*it3).second < minClassCnt ){
					remainingOrphants.insert( remainingOrphants.end(), class2 );
				}
			}// cnt < minClassCnt
		}//Iterator over all classes

		// perform class change
		changedPixels.clear();
		for(y=0; y<Get_NY() && Set_Progress(y); y++){
			for(x=0; x<Get_NX(); x++){
				class1 = pClassGrid->asInt(x,y); 
				it1 = classesToChange.find( class1 );
				if( it1 != classesToChange.end() ){
					class2 = (*it1).second;
					pClassGrid->Set_Value(x,y, class2 );
					// get minimal distance of class before aggregation and put to outputGrid
					if( (pOrphantsGrid != NULL) && pOrphantsGrid->asDouble(x,y) == pOrphantsGrid->Get_NoData_Value() ){
						it1 = classMinPixelMap.find( class1); 
						pi = (*it1).second;
						edist = pMinDist->asDouble(pi);
						pOrphantsGrid->Set_Value(x,y, edist );
					}
					// rember pixels for recalculating distances
					changedPixels.push_back( (y*Get_NY()+x) );
				}// found in request
				itS = destClasses.find( class1 );
				if( itS != destClasses.end() )
					changedPixels.push_back( (y*Get_NY()+x) );
			}//for x
		}//for y

		//recalculate distances and direction of updated pixels
		distMapT classDistMap;
		for( itV = changedPixels.begin(); itV!=changedPixels.end(); ++itV ){
			pi = (*itV);
			x = pi % Get_NX();
			y = pi / Get_NX();
			dist_dir = calcPixelsClassDistance(x,y,classDistMap);
			pMinDist->Set_Value(x,y, dist_dir[0] );
			pMinDistDir->Set_Value(x,y, dist_dir[1] );
			delete[] dist_dir;
		}//itV

	}while( remainingOrphants.size()>0 );
	//XXX loop if remaining oprhants
	*/

}

// renames classes, so that classes are numbered from 0 to max-1
// classPixelMap is not longer valid afterwards
void CCompleteLinkage::renameClasses(){
	long nr = 0, pi;
	pixelSetPtrT pPixSet;
	for( classPixelSetMapT::iterator it = classPixelMap.begin(); it!=classPixelMap.end(); ++it ){
		pPixSet = (*it).second;
		// constant (*it).first = nr;
		for( pixelSetT::iterator itp = pPixSet->begin(); itp != pPixSet->end(); ++itp ){
			pi = (*itp);
			pClassGrid->Set_Value(pi,nr);
		}
		nr++;
	}
	/*
	long maxC = (long) pClassGrid->Get_ZMax(); // maximum class number
	long class1, class2, cnum;
	map<long,long> cmap;
	map<long,long>::iterator it;
	int x,y;

	// calculate new numbers above maximum
	cnum = 0;
	for(y=0; y<Get_NY() && Set_Progress(y); y++){
		for(x=0; x<Get_NX(); x++){
			class1 = pClassGrid->asInt(x,y); 
			it = cmap.find( class1 );
			if( it == cmap.end() ){
				cnum++;
				class2 = cnum + maxC;
				cmap.insert( map<long,long>::value_type( class1, class2 ) );
			}else{
				class2 = (*it).second;
			}
			pClassGrid->Set_Value(x,y, class2 );
		} // x
	} // y

	// substract maximum
	for(y=0; y<Get_NY() && Set_Progress(y); y++){
		for(x=0; x<Get_NX(); x++){
			pClassGrid->Set_Value(x,y, pClassGrid->asInt(x,y) - maxC );
		} // x
	} // y
	*/
}


// initialized ClassesGrid and calculates first Distances and Directions
// returns pixel, that has minimum distance
long CCompleteLinkage::initClassesGrid(){
	prof.start("initClassesGrid");
	int x,y;
	double minDist, dir, edist;
	long ndvCnt = 0;

	typedef double double3[3];
	double3 *nextRowDist; // [ Get_NX() ][3] distances of next row to directions 0,1,2 of next rows pixel
	double prevDist; // previous pixels distance to direction 5
	double nextRowDist0; // temporary store for value nextRowDist[x+1][0]

	// no new for dynamic 2dimensional arrays, Get_NX() not constant
	nextRowDist = new double3[Get_NX()];
	// first row
	y = 0;{
		x = 0;{
			dir = 5;
			try{
				minDist = calculateEucledianDistance( x,y, x+1, y );
				prevDist = minDist;
			}catch( CNoDataValueError ){prevDist = DBL_MAX;}
			try{
				edist = calculateEucledianDistance( x,y, x, y+1 );
				nextRowDist[x][1] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 7;
				}
			}catch( CNoDataValueError ){nextRowDist[x][1] =  DBL_MAX;}
			try{
				edist = calculateEucledianDistance( x,y, x+1, y+1 );
				nextRowDist[x+1][0] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 8;
				}
			}catch( CNoDataValueError ){ nextRowDist[x+1][0] = DBL_MAX; }
			setInitDist(x,y,minDist,dir, ndvCnt);
		} // x = 0
		for(x=1; x<Get_NX()-1; x++){
			dir = 3;
			minDist = prevDist; //calculateEucledianDistance( x,y, x-1, y );
			try{
				edist = calculateEucledianDistance( x,y, x+1, y );
				prevDist = edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 5;
				}
			}catch( CNoDataValueError ){ prevDist = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x-1, y+1 );
				nextRowDist[x-1][2] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 6;
				}
			}catch( CNoDataValueError ){ nextRowDist[x-1][2] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x, y+1 );
				nextRowDist[x][1] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 7;
				}
			}catch( CNoDataValueError ){ nextRowDist[x][1] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x+1, y+1 );
				nextRowDist[x+1][0] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 8;
				}
			}catch( CNoDataValueError ){ nextRowDist[x+1][0] = DBL_MAX; }
			setInitDist(x,y,minDist,dir, ndvCnt);
		} // x [1..n-1]
		x = Get_NX()-1; {
			dir = 3;
			minDist = prevDist; //calculateEucledianDistance( x,y, x-1, y );
			try{
				edist = calculateEucledianDistance( x,y, x-1, y+1 );
				nextRowDist[x-1][2] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 6;
				}
			}catch( CNoDataValueError ){ nextRowDist[x-1][2] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x, y+1 );
				nextRowDist[x][1] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 7;
				}
			}catch( CNoDataValueError ){ nextRowDist[x][1] = DBL_MAX; }
			setInitDist(x,y,minDist,dir, ndvCnt);
		} // x = Get_NX)() -1
	}// y = 0

	//------------------ center rows ------------------
	for(y=1; y<Get_NY()-1 && Set_Progress(y); y++)	{
		x = 0;{
			dir = 1;
			minDist = nextRowDist[x][1];
			edist = nextRowDist[x][2];
			if( edist < minDist ){
				dir = 2;
				minDist = edist;
			}
			try{
				edist = calculateEucledianDistance( x,y, x+1, y );
				prevDist = edist;
				if( edist < minDist ){
					dir = 5;
					minDist = edist;
				}
			}catch( CNoDataValueError ){ prevDist = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x, y+1 );
				nextRowDist[x][1] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 7;
				}
			}catch( CNoDataValueError ){nextRowDist[x][1] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x+1, y+1 );
				nextRowDist0 =  edist; //nextRowDist[x+1][0];
				if( edist < minDist ){
					minDist = edist;
					dir = 8;
				}
			}catch( CNoDataValueError ){ nextRowDist0 = DBL_MAX; }
			setInitDist(x,y,minDist,dir, ndvCnt);
		}//x = 0
		for(x=1; x<Get_NX()-1; x++)		{
			dir = 0;
			minDist = nextRowDist[x][0]; nextRowDist[x][0] = nextRowDist0;
			edist = nextRowDist[x][1];
			if( edist < minDist ){
				dir = 1;
				minDist = edist;
			}
			edist = nextRowDist[x][2];
			if( edist < minDist ){
				dir = 2;
				minDist = edist;
			}
			edist = prevDist;
			if( edist < minDist ){
				dir = 3;
				minDist = edist;
			}
			try{
				edist = calculateEucledianDistance( x,y, x+1, y );
				prevDist = edist;
				if( edist < minDist ){
					dir = 5;
					minDist = edist;
				}
			}catch( CNoDataValueError ){ prevDist = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x-1, y+1 );
				nextRowDist[x-1][2] = edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 6;
				}
			}catch( CNoDataValueError ){ nextRowDist[x-1][2] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x, y+1 );
				nextRowDist[x][1] = edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 7;
				}
			}catch( CNoDataValueError ){ nextRowDist[x][1] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x+1, y+1 );
				//nextRowDist[x+1][0]; do not overwrite before use
				nextRowDist0 = edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 8;
				}
			}catch( CNoDataValueError ){ nextRowDist0 = 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;
			}
			try{
				edist = calculateEucledianDistance( x,y, x-1, y+1 );
				nextRowDist[x-1][2] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 6;
				}
			}catch( CNoDataValueError ){ nextRowDist[x-1][2] = DBL_MAX; }
			try{
				edist = calculateEucledianDistance( x,y, x, y+1 );
				nextRowDist[x][1] =  edist;
				if( edist < minDist ){
					minDist = edist;
					dir = 7;
				}
			}catch( CNoDataValueError ){ nextRowDist[x][1] = DBL_MAX; }
			setInitDist(x,y,minDist,dir, ndvCnt);
		} // x = Get_NX)() -1
	}// y [1..n-1]

	//---------------- last row ----------------------------
	y = Get_NY() -1;{
		x = 0;{
			dir = 1;
			minDist = nextRowDist[x][1];
			edist = nextRowDist[x][2];
			if( edist < minDist ){
				dir = 2;
				minDist = edist;
			}
			try{
				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 = 0
		for(x=1; x<Get_NX()-1; x++)		{
			dir = 0;
			minDist = nextRowDist[x][0]; 
			edist = nextRowDist[x][1];
			if( edist < minDist ){
				dir = 1;
				minDist = edist;
			}
			edist = nextRowDist[x][2];
			if( edist < minDist ){
				dir = 2;
				minDist = edist;
			}
			edist = prevDist;
			if( edist < minDist ){
				dir = 3;
				minDist = edist;
			}
			try{