Single flow direction (SFD) patch from Andrew Danner

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@32067 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.terraflow/IOStream/include/ami_stream.h

@@ -178,7 +178,7 @@ unsigned int AMI_STREAM<T>::get_block_length() {
 /**********************************************************************/
 /* creates a random file name, opens the file for reading and writing
    and and returns a file descriptor */
-int ami_single_temp_name(char *base, char* tmp_path);
+int ami_single_temp_name(const std::string& base, char* tmp_path);
 
 
 /**********************************************************************/

raster/r.terraflow/IOStream/include/mem_stream.h

@@ -84,7 +84,7 @@ off_t MEM_STREAM<T>::stream_len(void) {
 template<class T>
 AMI_err MEM_STREAM<T>::name(char **stream_name)  {
 
-  char *path = "dummy";
+  char const* path = "dummy";
 
   *stream_name = new char [strlen(path) + 1];
   strcpy(*stream_name, path);

raster/r.terraflow/IOStream/include/mm_utils.h

@@ -22,12 +22,12 @@
 
 
 #include "mm.h"
-
+#include <string>
 
 void  LOG_avail_memo();
 
 size_t getAvailableMemory();
 
-void  MEMORY_LOG(char* str);
+void  MEMORY_LOG(std::string str);
 
 #endif

raster/r.terraflow/IOStream/lib/src/ami_stream.cc

@@ -32,7 +32,7 @@
 /* creates a random file name, opens the file for reading and writing
    and and returns a file descriptor */
 int
-ami_single_temp_name(char *base, char* tmp_path) {
+ami_single_temp_name(const std::string& base, char* tmp_path) {
  
   char *base_dir;
   int fd;
@@ -41,7 +41,7 @@ ami_single_temp_name(char *base, char* tmp_path) {
   base_dir = getenv(STREAM_TMPDIR);
   assert(base_dir);
 
-  sprintf(tmp_path, "%s/%s_XXXXXX", base_dir, base);
+  sprintf(tmp_path, "%s/%s_XXXXXX", base_dir, base.c_str());
 #ifdef __MINGW32__
   fd = mktemp(tmp_path) ? open(tmp_path, O_CREAT|O_EXCL|O_RDWR, 0600) : -1;
 #else

raster/r.terraflow/IOStream/lib/src/mm_utils.cc

@@ -46,8 +46,7 @@ getAvailableMemory() {
   return fmem;
 }
 
-void 
-MEMORY_LOG(char* str) {
-  printf("%s", str);
+void MEMORY_LOG(std::string str) {
+  printf("%s", str.c_str());
   fflush(stdout);
 }

raster/r.terraflow/direction.cc

@@ -19,9 +19,10 @@
 
 #include "direction.h"
 #include "nodata.h"
+#include "common.h" /*for global opt->d8 flag*/
+
+#define TF_ROOTTWO 1.4142135623
 
-/***************************************************************/
-/* returns the direction corresponding to the window */
 /* directions:
    32 64 128
    16 *   1
@@ -33,6 +34,21 @@ encodeDirection(const genericWindow<elevation_type>& elevwin,
 				dimension_type row, 
 				dimension_type col) {
   
+  /*check global opt flag and call appropriate model*/
+	if(opt->d8){
+    return encodeDirectionSFD(elevwin, nrows, ncols, row, col);
+	}
+  return encodeDirectionMFD(elevwin, nrows, ncols, row, col);
+}
+
+/***************************************************************/
+/* returns the direction corresponding to the window using MFD */
+direction_type 
+encodeDirectionMFD(const genericWindow<elevation_type>& elevwin,
+				const dimension_type nrows, const dimension_type ncols,
+				dimension_type row, 
+				dimension_type col) {
+  
   direction_type dir = DIRECTION_UNDEF;
   
   if(!is_nodata(elevwin.get())) {
@@ -67,7 +83,71 @@ encodeDirection(const genericWindow<elevation_type>& elevwin,
   return dir;
 }
 
+/***************************************************************/
+/* returns the direction corresponding to the window using SFD */
+direction_type 
+encodeDirectionSFD(const genericWindow<elevation_type>& elevwin,
+				const dimension_type nrows, const dimension_type ncols,
+				dimension_type row, 
+				dimension_type col) {
 
+  direction_type dir = DIRECTION_UNDEF;
+	float tdrop, max_drop;
+	int max_indx;
+  
+  if(!is_nodata(elevwin.get())) {
+    dir = 0;
+		max_drop = 0;
+		max_indx = -1;
+		for(int i=0; i<9; i++){
+			if(i%2==1){ /*cardinal direction*/
+				tdrop=elevwin.get()-elevwin.get(i);
+				if(tdrop>max_drop){ max_drop=tdrop; max_indx=i; }
+			}
+			else if(i!=4){ /*diagonal*/
+				tdrop=(elevwin.get()-elevwin.get(i))/TF_ROOTTWO;
+				if(tdrop>max_drop){ max_drop=tdrop; max_indx=i; }
+			}
+		}
+		switch(max_indx){
+			case 0:
+			case 1:
+			case 2:
+				dir=32<<max_indx; break;
+		  case 3:
+				dir=16; break;
+			case 5:
+				dir=1; break;
+			case 6:
+			case 7:
+			case 8:
+				dir=8>>(max_indx-6); break;
+			default:
+				dir=0; break;
+		}
+  }
+  
+  /* if no direction, check for boundary */
+  if(dir==0 || dir==DIRECTION_UNDEF) {
+    if(row==0) {
+      dir = 64;
+    }
+    if(row==nrows-1) {
+      dir = 48;
+    }
+    if(col==0) {
+      if(row==0) dir = 32;
+      else if(row==nrows-1) dir = 8;
+      else dir = 16;
+    }
+    if(col==ncols-1) {
+      if(row==0) dir = 128;
+      else if(row==nrows-1) dir = 2;
+      else dir = 1;
+    }
+  }
+  return dir;
+}
 
 direction_type
 findDominant(direction_type dir) {
@@ -143,8 +223,25 @@ findDominant(direction_type dir) {
 	return 128;
   }
 
-
-  return dir;
+  /* Otherwise, there is no dominant direction.
+	 * SFD must output a single direction 1,2,4,8,16,32,64, or 128
+	 * pick the first matching one, with preference to cardinal direction
+	 */
+	if(dir & 85 ){ 
+		/* at least one cardinal direction (1+4+16+64=E+S+W+N) matches */
+		if(dir & 1){ return 1;}
+		if(dir & 4){ return 4;}
+		if(dir & 16){ return 16;}
+		if(dir & 64){ return 64;}
+	}
+	else{ 
+		/* 2 8 32 128 = SE SW NW NE*/
+		if(dir & 2){ return 2;}
+		if(dir & 8){ return 8;}
+		if(dir & 32){ return 32;}
+		if(dir & 128){ return 128;}
+  }
+  return dir; /* shouldn't get here unless dir <= 0 */
 }
 
 
@@ -152,7 +249,7 @@ char
 directionSymbol(direction_type dir) {
   char c='?';
   int cnt=0;
-  char *symbols = ">\\v/<\\^/";
+  char symbols[] = ">\\v/<\\^/";
 
   if(dir == 0) return '.';
   
@@ -201,3 +298,4 @@ directionSymbol(direction_type dir) {
 } 
 
 
+#undef TF_ROOTTWO

raster/r.terraflow/direction.h

@@ -141,6 +141,16 @@ direction_type encodeDirection(const genericWindow<elevation_type>& elevwin,
 			       const dimension_type ncols,
 			       dimension_type row, dimension_type col);
 
+direction_type encodeDirectionMFD(const genericWindow<elevation_type>& elevwin,
+			       const dimension_type nrows, 
+			       const dimension_type ncols,
+			       dimension_type row, dimension_type col);
+
+direction_type encodeDirectionSFD(const genericWindow<elevation_type>& elevwin,
+			       const dimension_type nrows, 
+			       const dimension_type ncols,
+			       dimension_type row, dimension_type col);
+
 direction_type findDominant(direction_type dir);
 char directionSymbol(direction_type dir);
 

raster/r.terraflow/fill.cc

@@ -20,8 +20,8 @@
 #include <ctype.h>
 #include <time.h>
 
+#include <string> 
 #include "fill.h"
-#include "option.h"
 #include "common.h"
 #include "water.h"
 #include "sortutils.h"
@@ -138,9 +138,9 @@ public:
 
 
 char *
-verbosedir(char *s) {
+verbosedir(std::string s) {
   static char buf[BUFSIZ];
-  sprintf(buf, "dump/%s", s);
+  sprintf(buf, "dump/%s", s.c_str());
   return buf;
 }
 

raster/r.terraflow/main.cc

@@ -111,7 +111,7 @@ parse_args(int argc, char *argv[]) {
   d8cut->key  = "d8cut";
   d8cut->type = TYPE_DOUBLE;
   d8cut->required = NO;
-  d8cut->answer = "infinity"; /* default value */
+  d8cut->answer = G_store("infinity"); /* default value */
   d8cut->description =
     _("If flow accumulation is larger than this value it is routed using "
       "SFD (D8) direction \n \t\t (meaningfull only  for MFD flow)");
@@ -122,7 +122,7 @@ parse_args(int argc, char *argv[]) {
   mem->key         = "memory";
   mem->type        = TYPE_INTEGER;
   mem->required    = NO;
-  mem->answer      = "300"; /* 300MB default value */
+  mem->answer      = G_store("300"); /* 300MB default value */
   mem->description = _("Maximum runtime memory size (in MB)");
 
   /* temporary STREAM path */
@@ -131,7 +131,7 @@ parse_args(int argc, char *argv[]) {
   streamdir->key        = "STREAM_DIR";
   streamdir->type       = TYPE_STRING;
   streamdir->required   = NO;
-  streamdir->answer     = "/var/tmp"; 
+  streamdir->answer     = G_store("/var/tmp"); 
   streamdir->description=
      _("Directory to hold temporary files (they can be large)");
 
@@ -150,7 +150,7 @@ parse_args(int argc, char *argv[]) {
   stats_opt->type       = TYPE_STRING;
   stats_opt->required   = NO;
   stats_opt->description= _("Name of file containing runtime statistics");
-  stats_opt->answer     = "stats.out";
+  stats_opt->answer     = G_store("stats.out");
 
 
   if (G_parser(argc, argv)) {