📄 completelinkage.cpp
字号:
///////////////////////////////////////////////////////////
// //
// SAGA //
// //
// System for Automated Geoscientific Analyses //
// //
// Module Library //
// //
// contrib_t_wutzler //
// //
//-------------------------------------------------------//
// //
// CompleteLinkage.cpp //
// //
// Copyright (C) 2003 Your Name //
// //
//-------------------------------------------------------//
// //
// This file is part of 'SAGA - System for Automated //
// Geoscientific Analyses'. SAGA is free software; you //
// can redistribute it and/or modify it under the terms //
// of the GNU General Public License as published by the //
// Free Software Foundation; version 2 of the License. //
// //
// SAGA is distributed in the hope that it will be //
// useful, but WITHOUT ANY WARRANTY; without even the //
// implied warranty of MERCHANTABILITY or FITNESS FOR A //
// PARTICULAR PURPOSE. See the GNU General Public //
// License for more details. //
// //
// You should have received a copy of the GNU General //
// Public License along with this program; if not, //
// write to the Free Software Foundation, Inc., //
// 59 Temple Place - Suite 330, Boston, MA 02111-1307, //
// USA. //
// //
//-------------------------------------------------------//
// //
// e-mail: your@e-mail.abc //
// //
// contact: Your Name //
// And Address //
// //
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#include "CompleteLinkage.h"
#include <math.h>
#include<LIMITS>
#include<list>
#include<set>
#include<iostream>
#include"DebugStream.h"
#include <stdexcept>
#include <sstream>
#include "NoDataValueError.h"
using std::runtime_error;
using std::deque;
#define MIN(X,Y) (((X)<(Y))?(X):(Y))
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
CCompleteLinkage::CCompleteLinkage(void)
{
//-----------------------------------------------------
// Place information about your module here...
Set_Name (_TL("CompleteLinkage"));
Set_Author (_TL("Copyrights (c) 2004 by Thomas Wutzler"));
Set_Description (_TW(
"CompleteLinkage\n"
));
//-----------------------------------------------------
//Parameters.Add_List(NULL, SG_T("InputGrids") , SG_T("Input Grids") , SG_T("Grids to input"), PARAMETER_INPUT, PARAMETER_TYPE_Grid);
//Parameters.Add_List(NULL, SG_T("InputWeights") , SG_T("Input Weights") , SG_T("Weights for each input Grid"), PARAMETER_INPUT, PARAMETER_TYPE_Int); does not work
//Parameters.Add_Grid( NULL, SG_T("Curvature") , SG_T("Curvature") , SG_T("Input Curvate grid.") , PARAMETER_INPUT);
CSG_Parameter *parentPar;
Parameters.Add_Grid( NULL, SG_T("InputGrid1") , SG_T("Input Grid1") , SG_T("Input Grid1 (e.g. profile curvatue)") , PARAMETER_INPUT);
parentPar = Parameters.Add_Grid( NULL, SG_T("InputGrid2") , SG_T("Input Grid2") , SG_T("Input Grid2 (e.g. curvatue)") , PARAMETER_INPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("InputWeight2"), SG_T("weighting factor"), SG_T("weighting factor of this grid compared to InputGrid1=1"), PARAMETER_TYPE_Double, 1.0 );
parentPar = Parameters.Add_Grid( NULL, SG_T("InputGrid3") , SG_T("Input Grid3") , SG_T("Input Grid3 (e.g. slope)") , PARAMETER_INPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("InputWeight3"), SG_T("weighting factor"), SG_T("weighting factor of this grid compared to InputGrid1=1"), PARAMETER_TYPE_Double, 1.0 );
parentPar = Parameters.Add_Grid( NULL, SG_T("InputGrid4") , SG_T("Input Grid4") , SG_T("Input Grid4") , PARAMETER_INPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("InputWeight4"), SG_T("weighting factor"), SG_T("weighting factor of this grid compared to InputGrid1=1"), PARAMETER_TYPE_Double, 1.0 );
parentPar = Parameters.Add_Grid( NULL, SG_T("InputGrid5") , SG_T("Input Grid5") , SG_T("Input Grid5") , PARAMETER_INPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("InputWeight5"), SG_T("weighting factor"), SG_T("weighting factor of this grid compared to InputGrid1=1"), PARAMETER_TYPE_Double, 1.0 );
Parameters.Add_Grid( NULL, SG_T("ClassesGrid") , SG_T("Classes") , SG_T("Output grid of aggregated classes") , PARAMETER_OUTPUT);
Parameters.Add_Grid( NULL, SG_T("OrphantsGrid"), SG_T("Orphants") , SG_T("Output grid of removed orphants") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Grid( NULL, SG_T("MinDistGrid"), SG_T("MinDist") , SG_T("Output grid of minimum distances to neighbor pixels class in parameter space") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Grid( NULL, SG_T("MinDistDirGrid"), SG_T("MinDistDir") , SG_T("Output grid of direction of minimum distances to neighbor pixels class in parameter space") , PARAMETER_OUTPUT_OPTIONAL);
//Parameters.Add_Value( NULL, SG_T("BOOLEAN") , SG_T("Boolean") , SG_T("A value of type boolean.") , PARAMETER_TYPE_Bool , true);
//Parameters.Add_Value( NULL, SG_T("INTEGER") , SG_T("Integer") , SG_T("A value of type integer.") , PARAMETER_TYPE_Int , 200);
//Parameters.Add_Value( NULL, SG_T("DOUBLE") , SG_T("Double") , SG_T("A floating point value.") , PARAMETER_TYPE_Double , 3.145);
//Parameters.Add_Value( NULL, SG_T("NoDataValue") , SG_T("NoDataValue") , SG_T("Value for 'no data'") , PARAMETER_TYPE_Double , -9999);
Parameters.Add_Value( NULL, SG_T("MaxOrphantsSize") , SG_T("MaxOrphantsSize") , SG_T("class of size up to this will be joined to next neighbouring class") , PARAMETER_TYPE_Int, 5);
Parameters.Add_Value( NULL, SG_T("GenFactor") , SG_T("GenFactor") , SG_T("Factor of Generalisation (range 0-1)") , PARAMETER_TYPE_Double , 0.99);
Parameters.Add_Value( NULL, SG_T("RegionExtent"), SG_T("RegionExtent"), SG_T("Extend of one dimension of algorithms subgrid") , PARAMETER_TYPE_Int , 30);
parentPar = Parameters.Add_Grid( NULL, SG_T("ClassesGrid1") , SG_T("Classes add1") , SG_T("additional output of classes at specified generalisation") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("GenFactor1") , SG_T("GenFactor add1") , SG_T("Factor of Generalisation (range 0-genFactor)") , PARAMETER_TYPE_Double , 0.75);
parentPar = Parameters.Add_Grid( NULL, SG_T("ClassesGrid2") , SG_T("Classes add2") , SG_T("additional output of classes at specified generalisation") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("GenFactor2") , SG_T("GenFactor add2") , SG_T("Factor of Generalisation (range 0-genFactor)") , PARAMETER_TYPE_Double , 0.9);
parentPar = Parameters.Add_Grid( NULL, SG_T("ClassesGrid3") , SG_T("Classes add3") , SG_T("additional output of classes at specified generalisation") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("GenFactor3") , SG_T("GenFactor add3") , SG_T("Factor of Generalisation (range 0-genFactor)") , PARAMETER_TYPE_Double , 0.95);
parentPar = Parameters.Add_Grid( NULL, SG_T("ClassesGrid4") , SG_T("Classes add4") , SG_T("additional output of classes at specified generalisation") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("GenFactor4") , SG_T("GenFactor add4") , SG_T("Factor of Generalisation (range 0-genFactor)") , PARAMETER_TYPE_Double , 0.975);
parentPar = Parameters.Add_Grid( NULL, SG_T("ClassesGrid5") , SG_T("Classes add5") , SG_T("additional output of classes at specified generalisation") , PARAMETER_OUTPUT_OPTIONAL);
Parameters.Add_Value( parentPar, SG_T("GenFactor5") , SG_T("GenFactor add5") , SG_T("Factor of Generalisation (range 0-genFactor)") , PARAMETER_TYPE_Double , 0.987);
/*
Parameters.Add_Choice( NULL, SG_T("METHOD") , SG_T("Method") , SG_T("Choose a method from this select option.",
"First Method|"
"Second Method|",
0
);
*/
//----------------------------------------
// Other init
for( int i = 0; i < MAX_SINGLE_PIXEL_ARR; i++ ){
//singlePixelSetArr[i].insert( -1 );
singlePixelSetArr[i].push_back( -1 );
}
}
//---------------------------------------------------------
CCompleteLinkage::~CCompleteLinkage(void)
{}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
// The only thing left to do is to realize your ideas whithin
// the On_Execute() function (which will be called automatically
// by the framework). But that's really your job :-)
template <class T>
std::string ToString ( const T& t ){
std::ostringstream oss;
oss << t;
return oss.str();
}
bool CCompleteLinkage::On_Execute(void){
// ----------------------------- declarations ---------------------
CSG_Grid *pOrphantsGrid;
long minPixelCnt, ndvCnt, nPi;
double genFactor;
double genFactorRegion = 0.8; // ratio of overall steps performed only in subgrid without calculating grids minimum
long i,nAll, ireg, nRegion; //aggStep, maximum Step for Region and overall
long nRegionStd, nRegionNX, nRegionNY, nRegionNXY; // number of steps for edge regions
int minx, miny; // pixel to be aggretate in one step
double minDist; // minimum pixels distance
distMapT classDistMap; // mapping class -> ptr(map class -> distance)
RegionSetT affectedRegionSet; // set of region pixels, thats minimum value has been changed
pixelSetT neighbourPixels; // set of pixels pointing to changed class
DebugStream dout;
//------------------------------- user input ----------------------
pInputGrid[0] = Parameters("InputGrid1")->asGrid();
pInputGrid[1] = Parameters("InputGrid2")->asGrid();
pInputGrid[2] = Parameters("InputGrid3")->asGrid();
pInputGrid[3] = Parameters("InputGrid4")->asGrid();
pInputGrid[4] = Parameters("InputGrid5")->asGrid();
inputWeight[0] = 1;
inputWeight[1] = Parameters("InputWeight2")->asDouble();
inputWeight[2] = Parameters("InputWeight3")->asDouble();
inputWeight[3] = Parameters("InputWeight4")->asDouble();
inputWeight[4] = Parameters("InputWeight5")->asDouble();
pOutputGridAdd[0] = Parameters("ClassesGrid1")->asGrid();
pOutputGridAdd[1] = Parameters("ClassesGrid2")->asGrid();
pOutputGridAdd[2] = Parameters("ClassesGrid3")->asGrid();
pOutputGridAdd[3] = Parameters("ClassesGrid4")->asGrid();
pOutputGridAdd[4] = Parameters("ClassesGrid5")->asGrid();
genFactor = Parameters("GenFactor")->asDouble();
regionExtent= Parameters("RegionExtent")->asInt();
minPixelCnt = Parameters("MaxOrphantsSize")->asInt() +1;
pClassGrid = Parameters("ClassesGrid")->asGrid();
pOrphantsGrid = Parameters("OrphantsGrid")->asGrid();
pMinDist = Parameters("MinDistGrid")->asGrid();
pMinDistDir = Parameters("MinDistDirGrid")->asGrid();
if( pMinDist == NULL ){
pMinDist = SG_Create_Grid(GRID_TYPE_Double, Get_NX(), Get_NY(), Get_Cellsize());
pMinDist->Set_Name(_TL("MinDist"));
}
if( pMinDistDir == NULL ){
pMinDistDir = SG_Create_Grid( GRID_TYPE_Double, Get_NX(), Get_NY(), Get_Cellsize());
pMinDistDir->Set_Name(_TL("MinDistDir"));
}
//------------- Region Grids ---------------
{int rnx = ceil( Get_NX() / (double) regionExtent ); // without (double) its an integer division
int rny = ceil( Get_NY() / (double) regionExtent );
pRegionMinDistGrid = SG_Create_Grid(GRID_TYPE_Double, rnx, rny, Get_Cellsize()*regionExtent);
pRegionMinDistGrid->Set_Name(_TL("MinDist Region"));
pRegionMinPixelGrid = SG_Create_Grid(GRID_TYPE_Int, rnx, rny, Get_Cellsize()*regionExtent);
pRegionMinPixelGrid->Set_Name(_TL("Regions MinPixel"));
pRegionNdvCntGrid = SG_Create_Grid(GRID_TYPE_Int, rnx, rny, Get_Cellsize()*regionExtent);
pRegionNdvCntGrid->Set_Name(_TL("Regions Count of No-DataValue-Pixels"));
}
try{
prof.resetAll();
// ---------------- initalize
// and calculate pixels and grids minimal distance
Message_Add(_TL("calculating initial distances\n"));
prof.start("init");
ndvCnt = initClassesGrid();
prof.stop("init");
// number of steps until output
nPi = Get_NCells()-ndvCnt;
nAll = (long) ( nPi * genFactor); // number of originial classes to remove
nAll = min(nAll,30000);
outputStep[0] = Parameters("GenFactor1")->asDouble() * nPi;
//pOutputGridAdd[0]->Set_Name( string("Classes ").assign( ToString( Parameters("GenFactor1")->asDouble())).c_str() );
outputStep[1] = Parameters("GenFactor2")->asDouble() * nPi;
//pOutputGridAdd[1]->Set_Name( string("Classes ").assign( ToString( Parameters("GenFactor2")->asDouble())).c_str() );
outputStep[2] = Parameters("GenFactor3")->asDouble() * nPi;
//pOutputGridAdd[2]->Set_Name( string("Classes ").assign( ToString( Parameters("GenFactor3")->asDouble())).c_str() );
outputStep[3] = Parameters("GenFactor4")->asDouble() * nPi;
//pOutputGridAdd[3]->Set_Name( string("Classes ").assign( ToString( Parameters("GenFactor4")->asDouble())).c_str() );
outputStep[4] = Parameters("GenFactor5")->asDouble() * nPi;
//pOutputGridAdd[4]->Set_Name( string("Classes ").assign( ToString( Parameters("GenFactor5")->asDouble())).c_str() );
nRegionStd = (long) regionExtent*regionExtent * genFactor * genFactorRegion;
//edge-Regions only can afford fewer runs
nRegionNX = (long) (Get_NX() % regionExtent) * regionExtent * genFactor * genFactorRegion;
nRegionNY = (long) regionExtent * (Get_NY() % regionExtent) * genFactor * genFactorRegion;
nRegionNXY = (long) (Get_NX() % regionExtent) * (Get_NY() % regionExtent) * genFactor * genFactorRegion;
// ------------------- region loops -------------------------
Message_Add(_TL("aggregating inside regions\n"));
i = 0;
for( int regy = 0; regy < pRegionMinDistGrid->Get_NY() && Set_Progress(regy); regy++ ){
for( int regx = 0; regx < pRegionMinDistGrid->Get_NX() ; regx++ ){
nRegion = nRegionStd;
if( regy == pRegionMinDistGrid->Get_NY()-1 )
nRegion = nRegionNY;
if( regx == pRegionMinDistGrid->Get_NX()-1 ){
if( regy == pRegionMinDistGrid->Get_NY()-1 )
nRegion = nRegionNXY;
else
nRegion = nRegionNX;
}
// decrease number of steps within region if there are no-data-values
ndvCnt = pRegionNdvCntGrid->asDouble(regx,regy);
nRegion -= (ndvCnt / genFactorRegion);
for( ireg = 0; ireg < nRegion; ireg++ ){
prof.start("regionLoop");
calcRegionsMinPixel(regx,regy, /*out*/ minx, miny, minDist );
if( minDist == DBL_MAX ){
prof.stop("regionLoop");
break;
}
if( minx == -1 )
throw new runtime_error("could not determind minPixel");
aggregatePixel(i++, minx,miny, classDistMap, affectedRegionSet);
prof.stop("regionLoop");
}
int regpi = REGXY2REGPI(regx,regy);
affectedRegionSet.insert( regpi );
}
}
Message_Add(_TL("Recalculating all Regions Minimum Distances\n"));
prof.start("recalcRegionMins_1");
recalcRegionMins(affectedRegionSet);
prof.stop("recalcRegionMins_1");
// ----------------- grids loop --------------------------------
// needs recalculation of several Regions
// and minimum loop over all regions
Message_Add(_TL("aggregating over all the grid\n"));
try{
for( i = i; (i<nAll) && Set_Progress_NCells(i-Get_NCells()); i++ ){
prof.start("mainLoop");
prof.start("checkOutput");
checkOutput(i);
prof.stop("checkOutput");
prof.start("calcMinDist");
recalcRegionMins(affectedRegionSet);
prof.stop("calcMinDist");
prof.start("getGridsMinPixel");
long minpi = getGridsMinPixel();
minx = PI2X(minpi);
miny = PI2Y(minpi);
prof.stop("getGridsMinPixel");
prof.start("aggregatePixel");
aggregatePixel(i, minx,miny, classDistMap, affectedRegionSet);
prof.stop("aggregatePixel");
prof.stop("mainLoop");
}
}catch( runtime_error e){
// no minpixel found
dout << e.what();
}
// -------------- remove Orphants ---------------------
//(classes with only few cells)
if( minPixelCnt > 1){
Message_Add(_TL("removing orphants\n"));
prof.start("removeOrphantClasses");
removeOrphantClasses(minPixelCnt, pOrphantsGrid, i, classDistMap, affectedRegionSet);
prof.stop("removeOrphantClasses");
}
// ----------------- rename classes --------------------
//prof.start("renameClasses");
Message_Add(_TL("renaming classes\n"));
renameClasses();
//prof.stop("renameClasses");
}catch( std::runtime_error e){
string str = e.what();
}
// -------------------- cleanup --------------
Message_Add(_TL("clean up temporary data structures\n"));
prof.start("cleanup");
{for( distMapT::iterator it = classDistMap.begin(); it != classDistMap.end(); ++it){
delete (*it).second;
}}
classDistMap.clear();
{for( classPixelSetMapT::iterator it = classPixelMap.begin(); it != classPixelMap.end(); ++it){
delete (*it).second;
}}
classPixelMap.clear();
//grids
delete pRegionMinDistGrid;
delete pRegionMinPixelGrid;
delete pRegionNdvCntGrid;
if( Parameters("MinDistGrid")->asGrid() == NULL)
delete pMinDist;
if( Parameters("MinDistDirGrid")->asGrid() == NULL)
delete pMinDistDir;
//DataObject_Add(pRegionMinDistGrid, true);
prof.stop("cleanup");
//-----------------------------------------------------
#ifdef _DEBUG
prof.output( dout );
#endif
return( true );
}
// for given pixel it returns the
// - the minimum distance to next distant cell (0),
// - the direction to the next distant cell (1)
// - the class of the next distant cell (2)
// returns 0, if there are no other neighboring classes, direction then is undefined
// do not forget to clear up the allocated double array by delete[]
double *CCompleteLinkage::calcPixelsClassDistance(int x, int y, distMapT& cache, long remClass1, long remClass2)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -