Add resampling code

inputs.cc

@@ -9,6 +9,78 @@
 #include "main.hh"
 #include "tiles.hh"
 
+/*
+ * resample_data
+ * This is used to resample tile data. It is particularly designed for
+ * use with LIDAR tiles where the resolution can be anything up to 2m.
+ * This function is capable of merging neighbouring pixel values
+ * The scaling factor is the distance to merge pixels.
+ * NOTE: This means that new resolutions can only increment in multiples of the original
+ * (ie 2m LIDAR can be 4/6/8/... and 20m can be 40/60)
+ */
+int resample_data(int scaling_factor){
+	short **new_data;
+	/* Tile width and height is guaranteed to be IPPD.
+	 * We now need to calculate the new width and height */
+	size_t new_ippd = ((IPPD+1) / scaling_factor);
+	max_elevation = -32768;
+	min_elevation = 32768;
+
+	for (int indx = 0; indx < MAXPAGES; indx++) {
+		/* Check if this is a valid tile. Uninitialized tiles have this default value */
+		if (dem[indx].max_el == 32768)
+			continue;
+		
+		/* Now we allocate the replacement data arrays */
+		if ((new_data = new short *[new_ippd]) == NULL) {
+			return ENOMEM;
+		}
+		for (size_t i = 0; i < new_ippd; i++) {
+			if ((new_data[i] = new short[new_ippd]) == NULL) {
+				return ENOMEM;	// If this happens we will leak but thats ok because we must bail anyway
+			}
+		}
+
+		/* Nearest neighbour normalization. For each subsample of the original, simply
+		 * assign the value in the top left to the new pixel */
+		for (size_t x = 0, i = 0; i < new_ippd; x += scaling_factor, i++) {
+			for (size_t y = 0, j = 0; j < new_ippd; y += scaling_factor, j++){
+				// fprintf(stderr,"[%d,%d,%d,%d] => %d\n", dem[indx].data[x][y], dem[indx].data[x][y+1], dem[indx].data[x+1][y], dem[indx].data[x+1][y+1], dem[indx].data[x][y]);
+				new_data[i][j] = dem[indx].data[x][y];
+			}
+		}
+
+		/* We need to fixup the min/max elevations */
+		for (size_t i = 0; i < new_ippd; i++) {
+			for (size_t j = 0; j < new_ippd; j++) {
+				if (new_data[i][j] > dem[indx].max_el)
+					dem[indx].max_el = new_data[i][j];
+				if (new_data[i][j] > max_elevation)
+					max_elevation = new_data[i][j];
+
+				if (new_data[i][j] < dem[indx].min_el)
+					dem[indx].min_el = new_data[i][j];
+				if (new_data[i][j] < min_elevation)
+					min_elevation = new_data[i][j];
+			}
+		}
+
+		/* Resampling complete, delete the original and assign the new values */
+		for (size_t i = 0; i < IPPD; i++)
+			delete [] dem[indx].data[i];
+		delete [] dem[indx].data;
+		dem[indx].data = new_data;
+	}
+
+	/* Report to the user */
+	fprintf(stderr, "Resampling IPPD %d->%d min/max el %d/%d\n", IPPD, new_ippd, min_elevation, max_elevation);
+
+	/* Finally, set the new IPPD value */
+	IPPD = new_ippd;
+	ippd = IPPD;
+	ARRAYSIZE = (MAXPAGES * IPPD) + 50;
+	return 0;
+}
 
 int loadClutter(char *filename, double radius, struct site tx)
 {
@@ -291,7 +363,7 @@ int loadLIDAR(char *filenames)
 			return errno;
 		fseek(fd, tiles[indx].datastart, SEEK_SET);
 		if (debug)
-			fprintf(stderr, "readLIDAR(fd,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", height, width, indx, yur, xur, yll, xll);
+			fprintf(stderr, "readLIDAR(fd,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", tiles[indx].height, tiles[indx].width, indx, tiles[indx].yur, tiles[indx].xur, tiles[indx].yll, tiles[indx].xll);
 		readLIDAR(fd, tiles[indx].height, tiles[indx].width, indx, tiles[indx].yur, tiles[indx].xur, tiles[indx].yll, tiles[indx].xll);
 	}
 

inputs.hh

@@ -3,6 +3,9 @@
 
 #include "common.h"
 
+/* Resample input tiles to new resolution */
+int resample_data(int scaling_factor);
+
 int LoadSDF_SDF(char *name, int winfiles);
 int LoadSDF(char *name, int winfiles);
 int LoadPAT(char *az_filename, char *el_filename);

main.cc

@@ -53,7 +53,7 @@ double earthradius, max_range = 0.0, forced_erp, dpp, ppd, yppd,
 int ippd, mpi,
     max_elevation = -32768, min_elevation = 32768, bzerror, contour_threshold,
     pred, pblue, pgreen, ter, multiplier = 256, debug = 0, loops = 100, jgets =
-    0, MAXRAD, hottest = 10, height, width;
+    0, MAXRAD, hottest = 10, height, width, resample;
 
 unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0;
 
@@ -1145,6 +1145,7 @@ int main(int argc, char *argv[])
 	max_txsites = 30;
 	fzone_clearance = 0.6;
 	contour_threshold = 0;
+	resample = -1;
 
 	ano_filename[0] = 0;
 	earthradius = EARTHRADIUS;
@@ -1322,6 +1323,14 @@ int main(int argc, char *argv[])
 			}
 		}
 
+		if (strcmp(argv[x], "-resample") == 0) {
+			z = x + 1;
+
+			if(!lidar)
+				fprintf(stderr, "[!] Warning, this should only be used with LIDAR tiles. Trying anyway\n");
+			sscanf(argv[z], "%d", &resample);
+		}
+
 		if (strcmp(argv[x], "-lat") == 0) {
 			z = x + 1;
 
@@ -1694,6 +1703,14 @@ int main(int argc, char *argv[])
 			exit(result);
 		}
 
+		/* If we have been asked to resample the input data; do it now. */
+		if (resample != -1 ){
+			if ((result = resample_data(resample)) != 0) {
+				fprintf(stderr, "Error resampling data\n");
+				return result;
+			}
+		}
+
 		if(debug){
 			fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%d x %d\n",max_north,min_west,min_north,max_west,width,height);
 		}
@@ -1710,9 +1727,9 @@ int main(int argc, char *argv[])
 
 	}else{
 		// DEM first
-                if(debug){
-                        fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west,max_lon,min_lon);
-                }
+		if(debug){
+			fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west,max_lon,min_lon);
+		}
 
 		max_lon-=3;
 

tiles.hh

@@ -23,6 +23,6 @@ typedef struct _tile_t{
 	int		*data;
 } tile_t, *ptile_t;
 
-int tile_load_lidar(tile_t*, char *, int);
+int tile_load_lidar(tile_t*, char *);
 
 #endif