From ed04d106579ea65189dfca9c9fe112aaa62fe91f Mon Sep 17 00:00:00 2001
From: Florian Wolf <florian.wolf@idiv.de>
Date: Thu, 6 Oct 2016 15:11:54 +0200
Subject: [PATCH] changes from Michael

- includes elevation in distance calculation
---
 makefile         |   2 +-
 makefile.mpi     |   6 +-
 src/DrawArea.cpp | 170 ++++++++++++++++++-------
 src/DrawArea.h   |  10 +-
 src/Error.h      |  35 ++++--
 src/Graph.cpp    |  26 ++--
 src/ParseArgs.h  |  48 ++++---
 src/pnm.h        | 319 ++++++++++++-----------------------------------
 8 files changed, 285 insertions(+), 331 deletions(-)

diff --git a/makefile b/makefile
index d1f3936..af176c5 100644
--- a/makefile
+++ b/makefile
@@ -1,5 +1,5 @@
 CC = g++
-CFLAGS = -Wall -O3 -c
+CFLAGS = -Wall -O3 -fpermissive -c
 LFLAGS = -s -lm
 
 PROG = DrawArea
diff --git a/makefile.mpi b/makefile.mpi
index e4967fd..e08e8ca 100644
--- a/makefile.mpi
+++ b/makefile.mpi
@@ -1,5 +1,5 @@
 CC = mpiCC -DDRAWAREA_MPI
-CFLAGS = -Wall -O3 -c
+CFLAGS = -Wall -O3 -fpermissive -c
 LFLAGS = -s -lm
 
 PROG = DrawArea-mpi
@@ -13,10 +13,10 @@ HDR = $(HDR_FILES:%=$(SRC_DIR)/%)
 OBJ = $(OBJ_FILES:%=$(OBJ_DIR)/%)
 
 $(PROG): $(OBJ) $(HDR)
-	$(CC) $(LDFLAGS) $(LFLAGS) -o $(PROG) $(OBJ)
+	$(CC) $(LFLAGS) -o $(PROG) $(OBJ)
 
 $(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp object
-	$(CC) $(CPPFLAGS) $(CFLAGS) $< -o $@
+	$(CC) $(CFLAGS) $< -o $@
 
 object:
 	mkdir -p $(OBJ_DIR)
diff --git a/src/DrawArea.cpp b/src/DrawArea.cpp
index f830068..0dc403a 100644
--- a/src/DrawArea.cpp
+++ b/src/DrawArea.cpp
@@ -22,10 +22,7 @@
 /*************************************************
 
     TODO-List:
-
 	- implement multi-threading
-	- change temporal distance calculation to
-	  topographic maps (incl. heuristic splitting)
 
 **************************************************/
 
@@ -33,15 +30,12 @@
 	#include "mpi.h"
 	#define MPI_MSG(msg) "Proc %d " msg, ID
 	#define MPI_ROOT_PROC 0
-	#define MPI_STDOUT stderr
 #else
 	#define MPI_MSG(msg) msg
-	#define MPI_STDOUT stdout
 #endif
 
 #include <stdio.h>
 #include <stdlib.h>
-#include <stdint.h>
 #include <time.h>
 #include <math.h>
 #include <string.h>
@@ -56,27 +50,28 @@ using namespace std;
 int ID;
 
 Node* ReadScopeFromMatrix(const char* fileName, Scope& scope);
-int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, size_t nColumns, size_t nRows, float limitTime);
+int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, size_t nColumns, size_t nRows, float limitTime, float pixelSize);
 int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope, float limitTime, float minTime, float areaScaling, int nBinning);
-Scope& FitBinning(Scope& scope, int nBinning, int xOffset, int yOffset);
-size_t* SplitToEqualParts(size_t nValues, double* values, size_t nParts, size_t *parts=NULL);
+Scope& Fit2Binning(Scope& scope, int nBinning, int xOffset, int yOffset);
+size_t* Split2EqualParts(size_t nValues, double* values, size_t nParts, size_t *parts=NULL);
+float MeanTravelTime(size_t *origin, size_t dataWidth, float pixelSize, int& validNeighbors);
 
 double PixelWeight(double time, double limit)
 {
 	return pow(max(1.0-time/limit,0.0),4);
 }
 
-char helpText[]="Usage: %s [OPTIONS] -l TIME -m TIME input\n"
+char helpText[]="Usage: %s [OPTIONS] -l TIME -p PIXELSIZE input\n"
 		"\nMandatories:\n"
 		"  -l, --limit TIME\n"
-			"\ttravel time limit\n"
-		"  -m, --min-time TIME\n"
-			"\tminimum time to pass single pixel\n"
+			"\ttravel time limit in seconds\n"
+		"  -p, --pixel-size PIXELSIZE\n"
+			"\tinput image pixel size in meters\n"
 		"\nOptions:\n"
 		"  -h, --help\n"
 			"\tprint this message\n"
 		"  -a, --area AREA\n"
-			"\tarea scaling for output data in pixels / AREA\n"
+			"\toutput data will be scaled to AREA square kilometers\n"
 		"  -b, --binning N\n"
 			"\tcalculate only every Nth point in each row and column\n"
 		"  -o, --output FILE\n"
@@ -97,7 +92,7 @@ char helpText[]="Usage: %s [OPTIONS] -l TIME -m TIME input\n"
 			"\tsplit map calculation into N threads\n"
 			"\tIf the number of columns or rows and threads is given, the requested number of threads is ignored.\n"
 #endif
-		"\nTIME and AREA must be positive numbers; N must be a positive integer; WIDTH, HEIGHT, LEFT, and TOP must be integers.\n";
+		"\nTIME, PIXELSIZE and AREA must be positive numbers; N must be a positive integer; WIDTH, HEIGHT, LEFT, and TOP must be integers.\n";
 
 
 
@@ -182,10 +177,10 @@ int main(int argc, char* argv[])
 
 		// split map into params.nRows*params.nColumns scopes
 		// heuristic will split in scopes with approx. equal calculation time
-		if ((error=SplitHeuristic(params.inputName, params.scope, scopeList, params.nColumns, params.nRows, params.limitTime))!=ERROR_SUCCESS)
+		if ((error=SplitHeuristic(params.inputName, params.scope, scopeList, params.nColumns, params.nRows, params.limitTime, params.pixelSize))!=ERROR_SUCCESS)
 			return error;
 		for (int i=0; i<params.nThreads; i++)
-			FitBinning(scopeList[i],params.nBinning,xOffsetBinning,yOffsetBinning);
+			Fit2Binning(scopeList[i],params.nBinning,xOffsetBinning,yOffsetBinning);
 
 #ifdef DRAWAREA_MPI
 
@@ -209,9 +204,9 @@ int main(int argc, char* argv[])
 	// send evaluated cmd args to all processes
 	MPI_Bcast(&params.inputName, 256, MPI_UNSIGNED_CHAR, MPI_ROOT_PROC, MPI_COMM_WORLD);
 	MPI_Bcast(&params.outputName, 256, MPI_UNSIGNED_CHAR, MPI_ROOT_PROC, MPI_COMM_WORLD);
-	MPI_Bcast(&params.areaScaling, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
 	MPI_Bcast(&params.limitTime, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
-	MPI_Bcast(&params.minTime, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
+	MPI_Bcast(&params.pixelSize, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
+	MPI_Bcast(&params.areaScaling, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
 	MPI_Bcast(&params.nBinning, 1, MPI_INT, MPI_ROOT_PROC, MPI_COMM_WORLD);
 	MPI_Bcast(&xOffsetBinning, 1, MPI_INT, MPI_ROOT_PROC, MPI_COMM_WORLD);
 	MPI_Bcast(&yOffsetBinning, 1, MPI_INT, MPI_ROOT_PROC, MPI_COMM_WORLD);
@@ -241,8 +236,8 @@ int main(int argc, char* argv[])
 
 	// perform local calculation
 	clock_t start=clock();
-	CalculateAreaMap(params.inputName, params.outputName, params.scope, params.limitTime, params.minTime, params.areaScaling, params.nBinning);
-	fprintf(MPI_STDOUT, MPI_MSG("map calculation took %.2f s\n"),float(clock()-start)/CLOCKS_PER_SEC);
+	CalculateAreaMap(params.inputName, params.outputName, params.scope, params.limitTime, params.pixelSize, params.areaScaling, params.nBinning);
+	Info(MPI_MSG("map calculation took %.2f s\n"),float(clock()-start)/CLOCKS_PER_SEC);
 
 #ifdef DRAWAREA_MPI
 	MPI_Finalize();
@@ -252,7 +247,7 @@ int main(int argc, char* argv[])
 }
 
 // shrink scope to matches binning offsets
-Scope& FitBinning(Scope& scope, int nBinning, int xOffset, int yOffset)
+Scope& Fit2Binning(Scope& scope, int nBinning, int xOffset, int yOffset)
 {
 	size_t last;
 	last=scope.Last(COLUMN)-(nBinning+scope.Last(COLUMN)%nBinning-xOffset)%nBinning;
@@ -267,7 +262,7 @@ Scope& FitBinning(Scope& scope, int nBinning, int xOffset, int yOffset)
 }
 
 // split a list of nValues values into nParts with approx. equal sum(values)
-size_t* SplitToEqualParts(size_t nValues, double* values, size_t nParts, size_t *parts)
+size_t* Split2EqualParts(size_t nValues, double* values, size_t nParts, size_t *parts)
 {
 	// get total weight
 	double totalWeight=0.0;
@@ -308,11 +303,51 @@ size_t* SplitToEqualParts(size_t nValues, double* values, size_t nParts, size_t
 	return parts;
 }
 
+// calculate the mean travel time for one pixel as mean of travel times to all valid neighbors
+float MeanTravelTime(size_t *origin, size_t dataWidth, float pixelSize, int& validNeighbors)
+{
+	validNeighbors=0;
+
+	// skip if pixel is invalid
+	if (*origin==INVALID_ALTITUDE)
+		return 0.0;
+	
+	float distance=0.0;
+	float diagPixelSize=SQRT_2*pixelSize;
+	size_t *neighbor;
+
+	// sum up travel times to all valid neighbors
+	for (int i=-1; i<2; i++)
+	{
+		neighbor=origin+i*dataWidth-1;
+		for (int j=-1; j<2; j++)
+		{
+			if ( (i!=0 || j!=0) && *neighbor!=INVALID_ALTITUDE )
+			{
+				distance+=Edge::Distance(*origin,*neighbor,(i*j!=0?diagPixelSize:pixelSize));
+				validNeighbors++;
+			}
+			neighbor++;
+		}
+	}
+
+	// calculate mean value
+	if (validNeighbors>0)
+		distance/=validNeighbors;
+
+	return distance;
+}
+
+
 // read input file and split into nColumns x nRows equally weighted scopes
-// 	- the weight of each scope is the sum of the evaluated pixel times to the power of 4
-// 	- the evaluated pixel time is MAX(1.0 - PIXEL_TIME / TIME_LIMIT, 0.0)
-int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, size_t nColumns, size_t nRows, float limitTime)
+// 	- the weight of each scope is the sum of the normalized approximated pixel times to the power of 4
+// 	- the approximated pixel time is the mean time needed to get to a neighboring pixel
+//	- normalized pixel time is MAX(1.0 - PIXEL_TIME / TIME_LIMIT, 0.0) times the number of valid neighbors
+int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, size_t nColumns, size_t nRows, float limitTime, float pixelSize)
 {
+	if (nColumns>1 || nRows>1)
+		Info("heuristic map splitting ...\n");
+
 	// check if scopeList is initialized
 	if (!scopeList)
 		return Error(-1, MPI_MSG("No memory allocated for scope list.\n"));
@@ -325,60 +360,109 @@ int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, si
 	inputScope.Fit(inputMap.GetSize());
 
 	// allocate memory for map evaluation
-	size_t *pixs=new size_t[inputScope.Width()];
+	size_t *invalid=new size_t[inputScope.Width()+2];
+	size_t *pixs=new size_t[3*inputScope.Width()+6];
 	double *weights=new double[std::max(inputScope.Height(),inputScope.Width())];
 	size_t *parts=new size_t[std::max(nRows,nColumns)+1];
-	if (!pixs || !weights || !parts)
+	if (!invalid || !pixs || !weights || !parts)
 		return Error(-1, MPI_MSG("Unable to allocate memory for heuristics.\n"));
+	
+	// initialize memory for map evaluation
+	for (size_t i=0; i<inputScope.Width()+2; i++) 
+		invalid[i]=INVALID_ALTITUDE;
+	size_t *prevPixs=pixs+1;
+	size_t *currPixs=pixs+inputScope.Width()+3;
+	size_t *nextPixs=pixs+2*inputScope.Width()+5;
+	memcpy(prevPixs-1,invalid,(inputScope.Width()+2)*sizeof(*pixs));
+	memcpy(currPixs-1,invalid,(inputScope.Width()+2)*sizeof(*pixs));
+	memcpy(nextPixs-1,invalid,(inputScope.Width()+2)*sizeof(*pixs));
+	int validNeighbors;
 
 	// row splitting first:
 	if (nRows>1)
 	{
+		Info("\tanalyzing rows ");
+
 		// get weight of each line in map and total map weight
 		Scope readScope=inputScope;
 		readScope.Height()=1;
+		inputMap.ReadScope(readScope,nextPixs);
 		for (size_t i=0; i<inputScope.Height(); i++)
 		{
 			weights[i]=0.0;
-			size_t num=inputMap.ReadScope(readScope,pixs);
-			for (size_t j=0; j<num; j++)
-				weights[i]+=PixelWeight(pixs[j],limitTime);
 			readScope.First(ROW)++;
+			memcpy(prevPixs,currPixs,inputScope.Width()*sizeof(*pixs));
+			memcpy(currPixs,nextPixs,inputScope.Width()*sizeof(*pixs));
+
+			size_t num=inputMap.ReadScope(readScope,nextPixs);
+			if (num<inputScope.Width())
+				memcpy(nextPixs+num,invalid,(inputScope.Width()-num)*sizeof(*pixs));
+
+			// get weight of each pixel by calculating travel times to all neighbors
+			for (size_t j=0; j<inputScope.Width(); j++)
+			{
+				float time=MeanTravelTime(currPixs+j,inputScope.Width(),pixelSize,validNeighbors);
+				weights[i]+=PixelWeight(time,limitTime)*validNeighbors;
+			}
+			if (i%1000==0) Info(".");
 		}
+		Info(" done\n");
 	}
 
-	SplitToEqualParts(inputScope.Height(),weights,nRows,parts);
+	// split map into equally weighted rows
+	Split2EqualParts(inputScope.Height(),weights,nRows,parts);
 	for (size_t row=0; row<nRows; row++)
 		for (size_t col=0; col<nColumns; col++)
 			scopeList[row*nColumns+col].GetRange(ROW)=Range(parts[row+1]-parts[row],inputScope.First(ROW)+parts[row]);
 
+	// reinitialize memory for map evaluation
+	memcpy(prevPixs-1,invalid,(inputScope.Width()+2)*sizeof(*pixs));
+	memcpy(currPixs-1,invalid,(inputScope.Width()+2)*sizeof(*pixs));
+	memcpy(nextPixs-1,invalid,(inputScope.Width()+2)*sizeof(*pixs));
+
 	// loop through all row bundles to split into columns
 	for (size_t row=0; row<nRows; row++)
 	{
 		if (nColumns>1)
 		{
+			Info("\tanalyzing columns of row %d ",row+1);
+
 			// get weight of each column in row bundle and total bundle weight
 			Scope readScope=inputScope;
 			readScope.Height()=1;
 			readScope.First(ROW)=scopeList[row*nColumns].First(ROW);
+			inputMap.ReadScope(readScope,nextPixs);
 			for (size_t i=0; i<inputScope.Width(); i++)
 				weights[i]=0.0;
 			for (size_t i=0; i<scopeList[row*nColumns].Height(); i++)
 			{
-				size_t num=inputMap.ReadScope(readScope,pixs);
-				for (size_t j=0; j<num; j++)
-					weights[j]+=PixelWeight(pixs[j],limitTime);
 				readScope.First(ROW)++;
+				memcpy(prevPixs,currPixs,inputScope.Width()*sizeof(*pixs));
+				memcpy(currPixs,nextPixs,inputScope.Width()*sizeof(*pixs));
+
+				size_t num=inputMap.ReadScope(readScope,nextPixs);
+				if (num<inputScope.Width())
+					memcpy(nextPixs+num,invalid,(inputScope.Width()-num)*sizeof(*pixs));
+
+				// get weight of each pixel by calculating travel times to all neighbors
+				for (size_t j=0; j<num; j++)
+				{
+					float time=MeanTravelTime(currPixs+j,inputScope.Width(),pixelSize,validNeighbors);
+					weights[j]+=PixelWeight(time,limitTime)*validNeighbors;
+				}
+				if (i%1000==0) Info(".");
 			}
+			Info(" done\n");
 		}
 
 		// get splitting and save in scopes
-		SplitToEqualParts(inputScope.Width(),weights,nColumns,parts);
+		Split2EqualParts(inputScope.Width(),weights,nColumns,parts);
 		for (size_t col=0; col<nColumns; col++)
 			scopeList[row*nColumns+col].GetRange(COLUMN)=Range(parts[col+1]-parts[col],parts[col]+inputScope.First(COLUMN));
 	}
 
 	inputMap.Close();
+	delete [] invalid;
 	delete [] pixs;
 	delete [] weights;
 	delete [] parts;
@@ -386,7 +470,6 @@ int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, si
 	return ERROR_SUCCESS;
 }
 
-
 // read Node array scope from map file
 // if scope is out of total map boundaries, adjust scope
 Node* ReadScopeBinaryMatrix(const char* fileName, Scope& scope)
@@ -410,7 +493,7 @@ Node* ReadScopeBinaryMatrix(const char* fileName, Scope& scope)
 		return NULL;
 	}
 
-	fprintf(MPI_STDOUT, MPI_MSG("reading nodes from scope %u x %u + %u + %u\n"),scope.Width(),scope.Height(),scope.First(COLUMN),scope.First(ROW));
+	Info(MPI_MSG("reading nodes from scope %u x %u + %u + %u\n"),scope.Width(),scope.Height(),scope.First(COLUMN),scope.First(ROW));
 
 	// read nodes scope from input map
 	input.ReadScope(scope,nodes),scope.Size();
@@ -426,7 +509,7 @@ Node* ReadScopeBinaryMatrix(const char* fileName, Scope& scope)
 //	  input map file
 //	- calculate area for each node in scope
 //	  and write it to output map file
-int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope, float limitTime, float minTime, float areaScaling, int nBinning)
+int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope, float limitTime, float pixelSize, float areaScaling, int nBinning)
 {
 	// open map file to update pixel data
 	PortableAnyMap output(outputName,PNM_MODE_UPDATE);
@@ -438,7 +521,7 @@ int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope
 
 	// enlarge scope by the maximal travel distance in TIME_LIMIT
 	// to get read scope (all nodes that may be needed)
-	Scope readScope=scope+int(limitTime/minTime);
+	Scope readScope=scope+int(limitTime/(pixelSize*MIN_TIME));
 
 	// read nodes from binary file
 	Node* nodes=ReadScopeBinaryMatrix(inputName,readScope);
@@ -446,8 +529,7 @@ int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope
 		return -1;
 
 	// connect nodes
-	Edge::minTime=minTime;
-	CreateGraphEdges(nodes,readScope.Width(),readScope.Height());
+	CreateGraphEdges(nodes,readScope.Width(),readScope.Height(),pixelSize);
 
 	// adjust scope if necessary
 	scope.Fit(readScope);
@@ -462,7 +544,7 @@ int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope
 	if (!area)
 		return Error(-1, MPI_MSG("No memory allocated for row data.\n"));
 
-	fprintf(MPI_STDOUT, MPI_MSG("calculating area in scope %u x %u + %u + %u\n"),scope.Width(),scope.Height(),scope.First(COLUMN),scope.First(ROW));
+	Info(MPI_MSG("calculating area in scope %u x %u + %u + %u ...\n"),scope.Width(),scope.Height(),scope.First(COLUMN),scope.First(ROW));
 
 	clock_t start;
 
@@ -480,7 +562,7 @@ int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope
 		output.UpdateScope(writeScope,area);
 		writeScope.First(ROW)++;
 
-		fprintf(MPI_STDOUT, "\n\trow %u (%u/%u) done in %.2f s\n",i+scope.First(ROW),(i/nBinning)+1,(scope.Height()-1)/nBinning+1,float(clock()-start)/CLOCKS_PER_SEC);
+		Info("\trow %u (%u/%u) done in %.2f s\n",i+scope.First(ROW),(i/nBinning)+1,(scope.Height()-1)/nBinning+1,float(clock()-start)/CLOCKS_PER_SEC);
 	}
 
 	output.Close();
diff --git a/src/DrawArea.h b/src/DrawArea.h
index 24a3565..409a35a 100644
--- a/src/DrawArea.h
+++ b/src/DrawArea.h
@@ -25,8 +25,12 @@
 #include <vector>
 #include <set>
 #include <stdio.h>
+#include <stdint.h>
 
-#define INVALID_ALTITUDE 0
+#define INVALID_ALTITUDE 55537
+#define MIN_TIME               0.72	// s/m
+#define POSITIVE_SLOPE_FACTOR  5	// s/m/100%
+#define NEGATIVE_SLOPE_FACTOR -2 	// s/m/-100%
 #define SQRT_2 1.41421356237
 
 // Node is the basic structure of the graph
@@ -35,18 +39,18 @@ struct Node {
 	std::vector<struct Edge> neighbors;		// vector of edges to all nearest and next-nearest neighbors
 	Node(size_t alt=0);
 	void AddNeighbour(int shift, float distance);
+	float MeanDistance();
 };
 
 // Edge is the basic structure to connect nodes of the graph
 struct Edge {
-	static float minTime;			// minimum pixel pass time
 	int shift;				// pointer shift to get from node to neighbor
 	float distance;				// temporal distance between the neighboring nodes
 	Edge(Node* start, int shift, float spacialDistance);
 	static float Distance(float start, float finish, float spacialDistance);	// calculating temporal edge length
 };
 
-void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows);
+void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows, float distance);
 
 // MapPoint is basically just a wrapper
 // to list, sort and compare WayPoint pointers
diff --git a/src/Error.h b/src/Error.h
index db35a0e..efae218 100644
--- a/src/Error.h
+++ b/src/Error.h
@@ -37,6 +37,7 @@
 void InitErrorHandling(const char* appName, const char* heltText, bool consoleErrorOut);
 int Error(const char* str, ...);
 int Error(int errCode, const char* str, ...);
+int Info(const char* str, ...);
 int Help(const char* msg=NULL);
 int Help(int code, const char* msg=NULL);
 
@@ -52,21 +53,17 @@ inline void InitErrorHandling(const char* appName, const char* helpText, bool co
 }
 
 // print error message
-inline int Error(int code, const char* str, ...)
+inline int ErrorVA(int code, const char* str, va_list pArg)
 {
-	// set prefix to error if error code != 0
-	va_list pArg;
-
+	va_list pArgCopy;
+	va_copy(pArgCopy, pArg);
 	// get length of error msg to be printed
-	va_start(pArg,str);
-	int len=vsnprintf(NULL,0,str,pArg)+1;
-	va_end(pArg);
+	int len=vsnprintf(NULL,0,str,pArgCopy)+1;
+	va_end(pArgCopy);
 
 	// allocate memory and write error message to string
-	va_start(pArg,str);
 	char* msg=(char*)malloc(len);
 	vsprintf(msg,str,pArg);
-	va_end(pArg);
 
 	// output error message
 	if (CONSOLE_ERROR_OUT)
@@ -79,12 +76,28 @@ inline int Error(int code, const char* str, ...)
 	return code;
 }
 
-// print error message with default error code -1
+// print error message wrapper
+inline int Error(int code, const char* str, ...)
+{
+	va_list pArg;
+	va_start(pArg,str);
+	return ErrorVA(code,str,pArg);
+}
+
+// print error message wrapper with default error code 0
+inline int Info(const char* str, ...)
+{
+	va_list pArg;
+	va_start(pArg,str);
+	return ErrorVA(ERROR_SUCCESS,str,pArg);
+}
+
+// print error message wrapper with default error code -1
 inline int Error(const char* str, ...)
 {
 	va_list pArg;
 	va_start(pArg,str);
-	return Error(-1,str,pArg);
+	return ErrorVA(-1,str,pArg);
 }
 
 // create help message
diff --git a/src/Graph.cpp b/src/Graph.cpp
index 9de40b5..dcf1b7b 100644
--- a/src/Graph.cpp
+++ b/src/Graph.cpp
@@ -21,7 +21,6 @@
 
 #include "DrawArea.h"
 
-
 // construct Node
 Node::Node(size_t alt)
 {
@@ -38,7 +37,6 @@ void Node::AddNeighbour(int shift, float distance)
 }
 
 // construct Edge as pointer-shift and temporal distance
-float Edge::minTime;
 Edge::Edge(Node* start, int shift, float spacialDistance)
 {
 	this->shift=shift;
@@ -48,14 +46,17 @@ Edge::Edge(Node* start, int shift, float spacialDistance)
 // calculate travel time along the edge based on altitude difference and spacial distance
 float Edge::Distance(float startAltitude, float endAltitude, float spacialDistance)
 {
-	return 0.5*spacialDistance*(startAltitude+endAltitude);
+	float dHeight=int16_t(endAltitude)-int16_t(startAltitude);
+	return MIN_TIME*spacialDistance+(dHeight>0.0?POSITIVE_SLOPE_FACTOR:NEGATIVE_SLOPE_FACTOR)*dHeight;
 }
 
 // connect each node with its 8 nearest neighbors (if existing)
-void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows)
+void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows, float distance)
 {
+	float diagDistance=SQRT_2*distance;
 	Node* node=nodes;
 	for (size_t i=0; i<numRows; i++)
+	{
         	for (size_t j=0; j<numCols; j++)
         	{
     			node->neighbors.clear();
@@ -67,27 +68,28 @@ void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows)
 
 			if (i>0)
 			{
-				node->AddNeighbour(-numCols,1.0);		// northern neighbor
+				node->AddNeighbour(-numCols,distance);		// northern neighbor
 				if (j>0)
-					node->AddNeighbour(-numCols-1,SQRT_2);	// north-western
+					node->AddNeighbour(-numCols-1,diagDistance);	// north-western
 				if (j+1<numCols)
-					node->AddNeighbour(-numCols+1,SQRT_2);	// north-eastern
+					node->AddNeighbour(-numCols+1,diagDistance);	// north-eastern
 			}
 
 			if (i+1<numRows)
 			{
-				node->AddNeighbour(numCols,1.0);		// southern neighbor
+				node->AddNeighbour(numCols,distance);		// southern neighbor
 				if (j>0)
-					node->AddNeighbour(numCols-1,SQRT_2);	// south-western
+					node->AddNeighbour(numCols-1,diagDistance);	// south-western
 				if (j+1<numCols)
-					node->AddNeighbour(numCols+1,SQRT_2);	// south-eastern
+					node->AddNeighbour(numCols+1,diagDistance);	// south-eastern
 			}
 
 			if (j>0)
-				node->AddNeighbour(-1,1.0);			// western neighbor
+				node->AddNeighbour(-1,distance);			// western neighbor
 			if (j+1<numCols)
-				node->AddNeighbour(1,1.0);			// eastern neighbor
+				node->AddNeighbour(1,distance);			// eastern neighbor
 
 			node++;
 		}
+	}
 }
diff --git a/src/ParseArgs.h b/src/ParseArgs.h
index 3f9ceb0..49ece6c 100644
--- a/src/ParseArgs.h
+++ b/src/ParseArgs.h
@@ -81,9 +81,9 @@ inline bool getCmdOption(char** begin, char** end, const char* option, float& va
 	return *error==0;
 }
 
-inline char* ExtractAppName(char* cmd)
+inline char* ExtractAppName(const char* cmd)
 {
-	char *fileName=cmd;
+	const char *fileName=cmd;
 	char *path=strrchr(fileName,PATH_SEPARATOR);		// get start pointer of filename w/o path
 	if (path) fileName=path+1;				// set filename start
 	int cCount=strlen(fileName);				// get length of file name
@@ -114,7 +114,7 @@ struct CmdArgStruct
 	int nBinning;
 	int nColumns;
 	float limitTime;
-	float minTime;
+	float pixelSize;
 	bool noCalc;
 	int nRows;
 	int nThreads;
@@ -127,16 +127,38 @@ struct CmdArgStruct
 // evaluated command line arguments
 inline int CmdArgStruct::Evaluate(int argc, char** argv)
 {
+
+/* mandatory arguments*/
+
 	// set input file
 	if (argc<2)
 		return Help(-1,"Missing input file.");
 	strcpy(inputName,argv[argc-1]);
 
+	// set time limit
+	if (getCmdOption(argv+1,argv+argc,"-l",limitTime) || getCmdOption(argv+1,argv+argc,"--limit",limitTime))
+	{
+		if (!limitTime>0)
+			return Help(-1,"Travel time limit must be a positive.");
+	}
+	else return Help(-1,"Missing travel time limit.");
+
+	// set pixel size
+	if (getCmdOption(argv+1,argv+argc,"-p",pixelSize) || getCmdOption(argv+1,argv+argc,"--pixel-size",pixelSize))
+	{
+		if (!pixelSize>0)
+			return Help(-1,"Pixel size must be a positive.");
+	}
+	else return Help(-1,"Missing image pixel size.");
+
+/* optional arguments*/
+
 	// set area scaling
 	areaScaling=1.0;
 	if (getCmdOption(argv+1,argv+argc,"-a",areaScaling) || getCmdOption(argv+1,argv+argc,"--area",areaScaling))
 		if (!areaScaling>0)
 			return Help(-1,"Area scaling factor must be a positive.");
+	areaScaling*=1E6/pixelSize/pixelSize;
 
 	// set binning
 	nBinning=1;
@@ -155,22 +177,6 @@ inline int CmdArgStruct::Evaluate(int argc, char** argv)
 				return Help(-1,"Number of columns must be a positive.");
 	}
 
-	// set time limit
-	if (getCmdOption(argv+1,argv+argc,"-l",limitTime) || getCmdOption(argv+1,argv+argc,"--limit",limitTime))
-	{
-		if (!limitTime>0)
-			return Help(-1,"Travel time limit must be a positive.");
-	}
-	else return Help(-1,"Missing travel time limit.");
-
-	// set minimum pixel pass time
-	if (getCmdOption(argv+1,argv+argc,"-m",minTime) || getCmdOption(argv+1,argv+argc,"--min-time",minTime))
-	{
-		if (!minTime>0)
-			return Help(-1,"Minimum time to pass single pixel must be a positive.");
-	}
-	else return Help(-1,"Missing minimum pixel time.");
-
 	// user requested only scope splitting but no map calculation
 	noCalc=false;
 	if (getCmdOption(argv+1,argv+argc,"-n") || getCmdOption(argv+1,argv+argc,"--no-calc"))
@@ -188,7 +194,7 @@ inline int CmdArgStruct::Evaluate(int argc, char** argv)
 	if (getCmdOption(argv+1,argv+argc,"-r") || getCmdOption(argv+1,argv+argc,"--rows"))
 	{
 		nRows=-1;
-		if (getCmdOption(argv+1,argv+argc,"-r",nRows) || getCmdOption(argv+1,argv+argc,"--cols",nRows))
+		if (getCmdOption(argv+1,argv+argc,"-r",nRows) || getCmdOption(argv+1,argv+argc,"--rows",nRows))
 			if (nRows<1)
 				return Help(-1,"Number of rows must be a positive.");
 	}
@@ -223,6 +229,8 @@ inline int CmdArgStruct::Evaluate(int argc, char** argv)
 		}
 	}
 
+	printf("cols: %d\nrows: %d\n",nColumns,nRows);
+
 	return ERROR_SUCCESS;
 }
 
diff --git a/src/pnm.h b/src/pnm.h
index 0b0e7f3..a471d6e 100644
--- a/src/pnm.h
+++ b/src/pnm.h
@@ -29,6 +29,7 @@
 // reading, writing, and (if binary P4-P6) updating
 // portable anymap images
 
+namespace {
 // Range describes a numbered interval
 class Range {
 
@@ -37,26 +38,57 @@ class Range {
 	size_t length;
 
     public:
-	Range(size_t length=0, size_t start=0);
-	size_t& Size();
-	size_t& First();
-	size_t Last();
-	Range& Fit(const Range limits);
-	bool IsInside(size_t value);
-
-	Range operator<<(int n);
-	Range operator>>(int n);
-	Range operator+(int n);
-	Range operator-(int n);
-	Range operator*(double x);
-	Range operator/(double x);
-
-	Range& operator<<=(int n);
-	Range& operator>>=(int n);
-	Range& operator+=(int n);
-	Range& operator-=(int n);
-	Range& operator*=(double x);
-	Range& operator/=(double x);
+	Range(size_t length=0, size_t start=0) { this->length=length; this->start=start; }
+	size_t& Size() { return length; }
+	size_t& First() { return start; }
+	size_t Last() { return start+length-1; }
+
+	Range& Fit(const Range limits)
+		{
+			size_t last=std::min(start+length-1,limits.start+limits.length-1);
+			start=std::max(start,limits.start);
+			length=(last<start?(size_t)0:last-start+1);
+			return *this;
+		}
+	bool IsInside(size_t value) { return (value>=start && value<start+length); }
+
+	Range& operator<<=(int n) { if (n<0) return *this>>=n; start=(start>(size_t)n?start-n:0); return *this; }
+	Range& operator>>=(int n) { if (n<0) return *this<<=n; start+=n; return *this; }
+	Range& operator+=(int n)
+		{
+			if (n<0) return *this-=n;
+			size_t oldLast=this->Last();
+			*this<<=n;
+			length=oldLast-start+n+1;
+			return *this;
+		}
+	Range& operator-=(int n)
+		{
+			if (n<0) return *this+=n;
+			if (length-2*n<0)
+			{
+				start=start+length/2;
+				length=0;
+			}
+			else
+			{
+				*this>>=n;
+				length-=2*n;
+			}
+			return *this;
+		}
+	Range& operator*=(double x) { if (x<0) x=-x; length*=x; return *this; }
+	Range& operator/=(double x) { if (x<0) x=-x; length/=x; return *this; }
+
+	Range operator<<(int n) { Range r=*this; return r<<=n; }
+	Range operator>>(int n) { Range r=*this; return r>>=n; }
+	Range operator+(int n) { Range r=*this; return r+=n; }
+	Range operator-(int n) { Range r=*this; return r-=n; }
+	Range operator*(double x) { Range r=*this; return r*=x; }
+	Range operator/(double x) { Range r=*this; return r/=x; }
+
+	bool operator==(Range r) { return start==r.start && length==r.length; }
+	bool operator!=(Range r) { return start!=r.start || length!=r.length; }
 };
 
 // Scope is a 2D Range
@@ -67,28 +99,33 @@ class Scope {
 	Range x,y;
 
     public:
-	Scope(Range x, Range y);
-	Scope(size_t xLength=0, size_t yLength=0, size_t xStart=0, size_t yStart=0);
-	Range& GetRange(DIMENSION d);
-	size_t Size();
-	size_t& First(DIMENSION d);
-	size_t& Height();
-	size_t& Width();
-	size_t Last(DIMENSION d);
-	Scope& Fit(const Scope limits);
-	bool IsInside(size_t x, size_t y);
-
-	Scope operator+(int n);
-	Scope operator-(int n);
-	Scope operator*(double x);
-	Scope operator/(double x);
-
-	Scope& operator+=(int n);
-	Scope& operator-=(int n);
-	Scope& operator*=(double x);
-	Scope& operator/=(double x);
+	Scope(Range x, Range y) { this->x=x; this->y=y; }
+	Scope(size_t xLength=0, size_t yLength=0, size_t xStart=0, size_t yStart=0) { x=Range(xLength,xStart); y=Range(yLength,yStart); }
+
+	Range& GetRange(DIMENSION d) { if (d==ROW) return y; return x; }
+	size_t Size() { return x.Size()*y.Size(); }
+	size_t& Height() { return y.Size(); }
+	size_t& Width() { return x.Size(); }
+	size_t& First(DIMENSION d) { if (d==ROW) return y.First(); return x.First(); }
+	size_t Last(DIMENSION d) { if (d==ROW) return y.Last(); return x.Last(); }
+
+	Scope& Fit(const Scope limits) { x.Fit(limits.x); y.Fit(limits.y); return *this; }
+	bool IsInside(size_t x, size_t y) { return this->x.IsInside(x) && this->y.IsInside(y); }
+
+	Scope& operator+=(int n) { this->x+=n; this->y+=n; return *this; }
+	Scope& operator-=(int n) { this->x-=n; this->y-=n; return *this; }
+	Scope& operator*=(double x) { this->x*=x; this->y*=x; return *this; }
+	Scope& operator/=(double x) { this->x/=x; this->y/=x; return *this; }
+
+	Scope operator+(int n) { Scope s=*this; return s+=n; }
+	Scope operator-(int n) { Scope s=*this; return s-=n; }
+	Scope operator*(double x) { Scope s=*this; return s*=x; }
+	Scope operator/(double x) { Scope s=*this; return s/=x; }
+
+	bool operator==(Scope s) { return x==s.x && y==s.y; }
+	bool operator!=(Scope s) { return x!=s.x || y!=s.y; }
 };
-
+}
 
 namespace {
 
@@ -610,6 +647,7 @@ size_t PortableAnyMap::ReadScope(Scope& readScope, T* data, T (*explicitCast)(si
 			if (feof(file))
 			{
 				SetError(8);
+				delete [] row;
 				return 0;
 			}
 
@@ -637,6 +675,7 @@ size_t PortableAnyMap::ReadScope(Scope& readScope, T* data, T (*explicitCast)(si
 				}
 			}
 		}
+		delete [] row;
 	}
 
 	return dataCount;
@@ -748,10 +787,12 @@ size_t PortableAnyMap::UpdateScope(Scope& updateScope, T* data, size_t (*explici
 		if (feof(file))
 		{
 			SetError(8);
+			delete [] row;
 			return 0;
 		}
 	}
 
+	delete [] row;
 	updateScope=writeScope;
 	return dataCount;
 }
@@ -867,200 +908,4 @@ inline bool PortableAnyMap::FillUp(size_t data_r, size_t data_g, size_t data_b)
 	return true;
 }
 }
-
-
-/* ---------------------------------------------------------------
-**********     Here comes the Range and Scope part!     **********
---------------------------------------------------------------- */
-
-inline Range::Range(size_t length, size_t start)
-{
-	this->length=length;
-	this->start=start;
-}
-
-inline Scope::Scope(Range x, Range y)
-{
-	this->x=x;
-	this->y=y;
-}
-
-inline Scope::Scope(size_t xLength, size_t yLength, size_t xStart, size_t yStart)
-{
-	x=Range(xLength,xStart);
-	y=Range(yLength,yStart);
-}
-
-inline size_t& Range::Size() { return length; }
-inline size_t Scope::Size() { return x.Size()*y.Size(); }
-inline size_t& Scope::Height() { return y.Size(); }
-inline size_t& Scope::Width() {	return x.Size(); }
-
-inline size_t& Range::First() { return start; }
-inline size_t Range::Last() { return start+length-1; }
-
-inline bool Range::IsInside(size_t value) { return (value>=start && value<start+length); }
-inline bool Scope::IsInside(size_t x, size_t y) { return this->x.IsInside(x) && this->y.IsInside(y); }
-
-inline Range& Range::Fit(const Range limits)
-{
-	size_t last=std::min(start+length-1,limits.start+limits.length-1);
-	start=std::max(start,limits.start);
-	length=(last<start?(size_t)0:last-start+1);
-	return *this;
-}
-
-inline Scope& Scope::Fit(const Scope limits)
-{
-	x.Fit(limits.x);
-	y.Fit(limits.y);
-	return *this;
-}
-
-inline Range& Scope::GetRange(DIMENSION d)
-{
-	if (d==ROW)
-		return y;
-	return x;
-}
-
-inline size_t& Scope::First(DIMENSION d)
-{
-	if (d==ROW)
-		return y.First();
-	return x.First();
-}
-
-inline size_t Scope::Last(DIMENSION d)
-{
-	if (d==ROW)
-		return y.Last();
-	return x.Last();
-}
-
-Range Range::operator<<(int n)
-{
-	Range r=*this;
-	return r<<=n;
-}
-Range Range::operator>>(int n)
-{
-	Range r=*this;
-	return r>>=n;
-}
-Range Range::operator+(int n)
-{
-	Range r=*this;
-	return r+=n;
-}
-Range Range::operator-(int n)
-{
-	Range r=*this;
-	return r-=n;
-}
-Range Range::operator*(double x)
-{
-	Range r=*this;
-	return r*=x;
-}
-Range Range::operator/(double x)
-{
-	Range r=*this;
-	return r/=x;
-}
-Range& Range::operator<<=(int n)
-{
-	if (n<0) return *this>>=n;
-	start=(start>(size_t)n?start-n:0);
-	return *this;
-}
-Range& Range::operator>>=(int n)
-{
-	if (n<0) return *this<<=n;
-	start+=n;
-	return *this;
-}
-Range& Range::operator+=(int n)
-{
-	if (n<0) return *this-=n;
-	size_t oldLast=this->Last();
-	*this<<=n;
-	length=oldLast-start+n+1;
-	return *this;
-}
-Range& Range::operator-=(int n)
-{
-	if (n<0) return *this+=n;
-	if (length-2*n<0)
-	{
-		start=start+length/2;
-		length=0;
-	}
-	else
-	{
-		*this>>=n;
-		length-=2*n;
-	}
-	return *this;
-}
-Range& Range::operator*=(double x)
-{
-	if (x<0)
-		x=-x;
-	length*=x;
-	return *this;
-}
-Range& Range::operator/=(double x)
-{
-	if (x<0)
-		x=-x;
-	length/=x;
-	return *this;
-}
-
-Scope Scope::operator+(int n)
-{
-	Scope s=*this;
-	return s+=n;
-}
-Scope Scope::operator-(int n)
-{
-	Scope s=*this;
-	return s-=n;
-}
-Scope Scope::operator*(double x)
-{
-	Scope s=*this;
-	return s*=x;
-}
-Scope Scope::operator/(double x)
-{
-	Scope s=*this;
-	return s/=x;
-}
-Scope& Scope::operator+=(int n)
-{
-	this->x+=n;
-	this->y+=n;
-	return *this;
-}
-Scope& Scope::operator-=(int n)
-{
-	this->x-=n;
-	this->y-=n;
-	return *this;
-}
-Scope& Scope::operator*=(double x)
-{
-	this->x*=x;
-	this->y*=x;
-	return *this;
-}
-Scope& Scope::operator/=(double x)
-{
-	this->x/=x;
-	this->y/=x;
-	return *this;
-}
-
 #endif // PGM_H_INCLUDED
-- 
GitLab