diff --git a/makefile b/makefile index d1f3936815888d8824d83806bed2899d76dd160d..af176c5f489df00a983de522ca100e4ec16a39bd 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 e4967fdf5fad8b67af41e0066cd9de6fa41dc679..e08e8cab908a3e2ace42a63a34bce7eeb42e0ee6 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 f8300687f1c42df77801bb48a342138732fb2975..0dc403ab279ca01c5e098228e579ad466001ee0e 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(¶ms.inputName, 256, MPI_UNSIGNED_CHAR, MPI_ROOT_PROC, MPI_COMM_WORLD); MPI_Bcast(¶ms.outputName, 256, MPI_UNSIGNED_CHAR, MPI_ROOT_PROC, MPI_COMM_WORLD); - MPI_Bcast(¶ms.areaScaling, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD); MPI_Bcast(¶ms.limitTime, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD); - MPI_Bcast(¶ms.minTime, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD); + MPI_Bcast(¶ms.pixelSize, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD); + MPI_Bcast(¶ms.areaScaling, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD); MPI_Bcast(¶ms.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 24a35658bb1e6df53728ddcf741bb12926170a0b..409a35ae9c2623d595bf198610ab256c5dcec510 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 db35a0e21e1c179e2d47d18edd5fb953cb53b648..efae218406ca207567cf51b60cb7cfd22de9b305 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 9de40b54e76489d75ee3a9f52089b49658bb3ea2..dcf1b7b00c1f0b873a45ce466e20dedae66143d7 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 3f9ceb033a3cb017c8b754501d798d940e9d6100..49ece6c122737cbddbb9390ec1daadf3179f4e2e 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 0b0e7f351ecb8196146af5dafe787d749eb733e4..a471d6ed80d0f1633cc81a8e84af210b44a046d6 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