2.92 Bigger, faster LIDAR

CHANGELOG

@@ -1,5 +1,9 @@
 SIGNAL SERVER CHANGELOG 
 
+2.92 - 29 September 2016
+Overhauled LIDAR functions to read in (more) native format tiles.
+Replaced 'super tile' code with more scalable multi-tile code: currently limited to 16.
+
 2.9 - 01 September 2016
 MODIS Landcover support in ASCII Grid format. Works for 1200/3600 SRTM only.
 

inputs.cc

@@ -6,6 +6,7 @@
 #include "common.h"
 #include "main.hh"
 
+
 int loadClutter(char *filename, double radius, struct site tx)
 {
 	/* This function reads a MODIS 17-class clutter file in ASCII Grid format.
@@ -124,37 +125,29 @@ int loadClutter(char *filename, double radius, struct site tx)
 	return 0;
 }
 
-void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx,
+void readLIDAR(FILE *fd, int h, int w, int indx,double n, double e, double s, double west)
-	       double n, double e, double s, double west)
 {
-	int x = 0, y = 0, reads = 0;
+	int x = 0, y = 0, reads = 0, a=0, b=0, avg=0, tWidth = 0, tHeight = 0;
 	char line[150000];
 	char *pch;
 
-	// use offsets to determine middle lat/lon for 5 x 10k tiles
-	// TALL
-	if (voffset == 0 && h == 10000) {
-		s = (s+n)/2;
-		h = 5000;
-	}
-	if (voffset == 5000 && h == 10000) {
-		n = (s+n)/2;
-	}
-	// WIDE
-	if (hoffset == 0 && w == 10000) {
-		e = (e+west)/2;
-		w = 5000;
-	}
-	if (hoffset == 5000 && w == 10000) {
-		west = (e+west)/2;
-		w = 5000;
-	}
-
 	dem[indx].max_north=n;
 	dem[indx].min_west=e;
 	dem[indx].min_north=s;
 	dem[indx].max_west=west;
 
+	if (max_north == -90)
+		max_north = dem[indx].max_north;
+
+	else if (dem[indx].max_north > max_north)
+		max_north = dem[indx].max_north;
+
+	if (min_north == 90)
+		min_north = dem[indx].min_north;
+
+	else if (dem[indx].min_north < min_north)
+		min_north = dem[indx].min_north;
+
 	if (dem[indx].max_west > max_west)
 		max_west = dem[indx].max_west;
 	if (dem[indx].min_west < min_west)
@@ -186,18 +179,12 @@ void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx,
 
 	for (y = h-1; y > -1; y--) {
 		x = w-1;
+	
 		if (fgets(line, 150000, fd) != NULL) {
-			// do nothing until n rows have passed
+				pch = strtok(line, " "); // split line into values
-			if (y < voffset || voffset == 0) {
-				pch = strtok(line, " ");
-
-				//dummy reads until we reach offset
-				// for 5000 offset, width must be 10e3
-				for (n = 0; n < hoffset; n++)
-					pch = strtok(NULL, " ");
 
 				while (pch != NULL && x > -1) {
-					if (atoi(pch) < -999)
+					if (atoi(pch) <= -9999)
 						pch = "0";
 					dem[indx].data[y][x] = atoi(pch);
 					dem[indx].signal[x][y] = 0;
@@ -210,46 +197,32 @@ void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx,
 						dem[indx].min_el = atoi(pch);
 						min_elevation = dem[indx].min_el;
 					}
-					x--;
+
+   					x--;
 					pch = strtok(NULL, " ");
 				}//while
-			}//voffset
+
+
 		} else {
 			fprintf(stdout, "LIDAR error @ x %d y %d indx %d\n",
 					x, y, indx);
 		}//if
 	}//for
+
 }
 
 int loadLIDAR(char *filenames)
 {
-	/*
-	 * This function reads either 9 LIDAR tiles of n rows and n columns
-	 * in ASCII grid format OR a single super tile composed of 2 or more
-	 * tiles. The tile must have WGS84 bounds in the header in the order:
-	 * WEST,SOUTH,EAST,NORTH
-	 *	ncols        5000
-	 *	nrows        5000
-	 *	xllcorner    -2.291359
-	 *	yllcorner    51.788295
-	 *	xurcorner    -2.146674
-	 *	yurcorner    51.878474
-	 *	cellsize     2
-	 *	NODATA_value  -9999
-	 *
-	 *	Tiles must be entered in the format
-	 *
-	 *	    -lid tile1.asc,tile2.asc,tile3.asc
-	 */
 	char *filename;
-	char *files[4]; // 4 tiles
+	char *files[25]; // 25 tiles
 	int x, y, indx = 0, fc = 0, hoffset = 0, voffset = 0, pos,
 	    dem_alloced = 0;
-	double xll, yll, xur, yur, cellsize;
+	double xll, yll, xur, yur, cellsize, avgCellsize;
 	char found, free_page = 0, jline[20], lid_file[255],	
 	     path_plus_name[255], *junk = NULL;
 	char line[50000];
 	char * pch;
+    	double TO_DEG = (180 / PI);
 	FILE *fd;
 
 	// test for multiple files
@@ -265,10 +238,10 @@ int loadLIDAR(char *filenames)
 
 		if (fd != NULL) {
 
-			if (fgets(line, 19, fd) != NULL) {
+			if (fgets(line, 255, fd) != NULL) {
 				pch = strtok (line," ");
 				pch = strtok (NULL, " ");
-				width = atoi(pch);
+				width = atoi(pch); // ncols
 
 				if (debug) {
 					fprintf(stdout, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, width);
@@ -276,63 +249,47 @@ int loadLIDAR(char *filenames)
 				}
 
 				if (!dem_alloced) {
-					IPPD = width;
+					IPPD = width * 1.1; // +10%
 					ARRAYSIZE = (MAXPAGES * IPPD) + 10;
 					do_allocs();
 					dem_alloced = 1;
 				}
 
 			}
-			if (fgets(line, 19, fd) != NULL)
+			if (fgets(line, 255, fd) != NULL)
-				height = atoi(pch);
+				height = atoi(pch); // nrows
-			fgets(line, 24, fd);
 
-			if (fgets(line, 24, fd) != NULL) {
+			if (fgets(line, 255, fd) != NULL) {
-				//xll=atof(pch);
+				sscanf(pch, "%lf", &xll); // xll
-				sscanf(pch, "%lf", &xll);
-			}
-			fgets(line, 24, fd);
-			if (fgets(line, 24, fd) != NULL) {
-				//yll=atof(pch);
-				sscanf(pch, "%lf", &yll);
 			}
 
-			fgets(line, 24, fd);
+			if (fgets(line, 255, fd) != NULL) {
-
+				sscanf(pch, "%lf", &yll); // yll
-			if (fgets(line, 24, fd) != NULL) {
-				//xur=atof(pch);
-				sscanf(pch, "%lf", &xur);
 			}
 
-			fgets(line, 24, fd);
+			if (fgets(line, 255, fd) != NULL)
+				sscanf(pch, "%lf", &cellsize); 
+			avgCellsize=avgCellsize+cellsize;
 
-			if (fgets(line, 24, fd) != NULL) {
+			if(cellsize>=0.5){ // 50cm LIDAR?
-				//yur=atof(pch);
+				// compute xur and yur with inverse haversine if cellsize in *metres*
-				sscanf(pch, "%lf", &yur);
+				double roundDistance = (width*cellsize)/6371000;
+				yur = asin(sin(yll*DEG2RAD) * cos(roundDistance) + cos(yll * DEG2RAD) * sin(roundDistance) * cos(0)) * TO_DEG;
+				xur = ((xll*DEG2RAD) + atan2(sin(90*DEG2RAD) * sin(roundDistance) * cos(yll*DEG2RAD), cos(roundDistance) - sin(yll * DEG2RAD) * sin(yur*DEG2RAD))) * TO_DEG;
+			}else{
+				// Degrees with GDAL option: -co "FORCE_CELLSIZE=YES"
+				xur = xll+(cellsize*width);
+				yur = yll+(cellsize*height);
 			}
-
+	
-			fgets(line, 15, fd);
+			if (xur > eastoffset)
-
-			if (fgets(line, 15, fd) != NULL)
-				cellsize = strtod(pch, NULL);
-
-			// LIDAR 10m @ 10800 PPD
-			if (cellsize == 10.0)
-				MAXRAD = 30;
-			// LIDAR 2m @ 54000 PPD
-			if (cellsize == 2.0)
-				MAXRAD = 15;
-			// LIDAR 1m @ 108000 PPD!
-			if (cellsize == 1.0)
-				MAXRAD = 10;
-
-			if (xur < eastoffset)
 				eastoffset = xur;
-			if (xll > westoffset)
+			if (xll < westoffset)
 				westoffset = xll;
 
 			if (debug)
-				fprintf(stdout, "PRE yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta);
+				fprintf(stdout,"%d, %d, %.7f, %.7f, %.0f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur);
+
 
 			// Greenwich straddling hack
 			if (xll < 0 && xur > 0) {
@@ -353,58 +310,16 @@ int loadLIDAR(char *filenames)
 			if (debug)
 				fprintf(stdout, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta);
 
-			if (yll < min_north)
-				min_north = yll;
-			if (yur > max_north)		
-				max_north = yur;
 
-			fgets(line, 30, fd); // NODATA
+			fgets(line, 255, fd); // NODATA
 			pos = ftell(fd);
 
 			// tile 0 [x| ]
 			if (debug)
-				fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", 0, 0, height, width, indx, yur, xur, yll, xll);
+				fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", height, width, indx, yur, xur, yll, xll);
 
-			readLIDAR(fd, 0, 0, height, width, indx, yur, xur, yll,	xll);
+			readLIDAR(fd, height, width, indx, yur, xur, yll, xll);
-			//rewind
-			fseek(fd, pos, SEEK_SET);
 
-			// tile 1 [ |x]
-			if (width == 2000 ||
-			    width == 10000 ||
-			    width == 10082) {
-				indx++;
-				if (debug)
-					fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", width / 2, 0, height, width, indx, yur, xur, yll, xll);
-
-				readLIDAR(fd, width / 2, 0, height, width, indx, yur, xur, yll, xll);
-			}
-			//rewind
-			fseek(fd, pos, SEEK_SET);
-
-			// tile 2 [x | ]
-			if (height == 2000 ||
-			    height == 10000 ||
-			    height == 10082) {
-				indx++;
-				if (debug)
-					fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", 0, height / 2, height, width, indx, yur, xur, yll, xll);
-
-				readLIDAR(fd, 0, height / 2, height, width, indx, yur, xur, yll, xll);
-			}
-			//rewind
-			fseek(fd, pos, SEEK_SET);
-
-			// tile 3 [ |x]
-			if ((width == 2000 && height == 2000) ||
-			    (width == 10000 && height == 10000) ||
-			    (width == 10082 && height == 10082)) {
-				indx++;
-				if (debug)
-					fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", width / 2, height / 2, height, width, indx, yur, xur, yll, xll);
-
-				readLIDAR(fd, width / 2, height / 2, height, width, indx, yur, xur, yll, xll);
-			}
 			fclose(fd);
 			if (debug)
 				fprintf(stdout, "LIDAR LOADED %d x %d\n", width, height);
@@ -414,6 +329,15 @@ int loadLIDAR(char *filenames)
 		indx++;
 	} // filename(s)
 
+	// TESTING!
+	IPPD=width;
+	ippd=width;
+	avgCellsize=avgCellsize/fc;
+	height = (unsigned)((max_north-min_north) / avgCellsize);
+	width = (unsigned)((max_west-min_west) / avgCellsize);
+
+	if (debug)
+		fprintf(stdout, "fc %d WIDTH %d HEIGHT %d ippd %d minN %.5f maxN %.5f minW %.5f maxW %.5f\n", fc, width, height, ippd,min_north,max_north,min_west,max_west);
 	return 0;
 }
 

main.cc

@@ -1,4 +1,4 @@
-double version = 2.91;
+double version = 2.92;
 /****************************************************************************\
 *  Signal Server: Radio propagation simulator by Alex Farrant QCVS, 2E0TDW   *
 ******************************************************************************
@@ -45,7 +45,7 @@ char string[255], sdf_path[255], udt_file[255], opened = 0, gpsav =
 double earthradius, max_range = 0.0, forced_erp, dpp, ppd, yppd,
     fzone_clearance = 0.6, forced_freq, clutter, lat, lon, txh, tercon, terdic,
     north, east, south, west, dBm, loss, field_strength,
-	min_north = 90, max_north = -90, min_west = 360, max_west = -1, westoffset=-180, eastoffset=180, delta=0, rxGain=0;
+	min_north = 90, max_north = -90, min_west = 360, max_west = -1, westoffset=180, eastoffset=-180, delta=0, rxGain=0;
 
 int ippd, mpi,
     max_elevation = -32768, min_elevation = 32768, bzerror, contour_threshold,
@@ -1024,7 +1024,7 @@ int main(int argc, char *argv[])
 	unsigned char LRmap = 0, txsites = 0, topomap = 0, geo = 0, kml =
 	    0, area_mode = 0, max_txsites, ngs = 0;
 
-	char mapfile[255], udt_file[255], ano_filename[255], lidar_tiles[512], clutter_file[255];
+	char mapfile[255], udt_file[255], ano_filename[255], lidar_tiles[1024], clutter_file[255];
 
 	double altitude = 0.0, altitudeLR = 0.0, tx_range = 0.0,
 	    rx_range = 0.0, deg_range = 0.0, deg_limit = 0.0, deg_range_lon;
@@ -1036,9 +1036,9 @@ int main(int argc, char *argv[])
 	}
 
 	if (strstr(argv[0], "signalserverLIDAR")) {
-		MAXPAGES = 4; 
+		MAXPAGES = 25; // 5x5
 		lidar = 1;
-		IPPD = 5000;
+		IPPD = 500; // will be overridden based upon file header...
 	}
 
 	strncpy(ss_name, "Signal Server\0", 14);
@@ -1227,7 +1227,7 @@ int main(int argc, char *argv[])
 			z = x + 1;
 			lidar=1;
 			if (z <= y && argv[z][0] && argv[z][0] != '-')
-				strncpy(lidar_tiles, argv[z], 253);
+				strncpy(lidar_tiles, argv[z], 1022);
 		}
 
 		if (strcmp(argv[x], "-res") == 0) {
@@ -1616,14 +1616,18 @@ int main(int argc, char *argv[])
 				"correct and try again.\n");
 			exit(EXIT_FAILURE);
 		}
-		ippd = IPPD;
+
 
 		if(debug){
-			fprintf(stdout,"%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west);
+			fprintf(stdout,"%.4f,%.4f,%.4f,%.4f,%d x %d\n",max_north,min_west,min_north,max_west,width,height);
 		}
 		ppd=rint(height / (max_north-min_north));
 		yppd=rint(width / (max_west-min_west));
 
+		//ppd=(double)ippd;
+		//yppd=ppd; 
+
+
 		if(delta>0){
 			tx_site[0].lon+=delta;
 		}

outputs.cc

@@ -567,18 +567,12 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml,
 
 				x0 = (int)rint((ppd *
 					      (lat -
-						(double)dem[indx].min_north))); // +4549 fix
+						(double)dem[indx].min_north))); 
 				y0 = mpi -
-				    (int)rint(ppd * // ONLY exception. All others are yppd
+				    (int)rint(ppd * 
 					      (LonDiff
 					       ((double)dem[indx].max_west,lon)));
 
-				// fix for multi-tile lidar
-				if(width==10000 && (indx==1 || indx==3)){
-					if(y0 >= 3432){ //3510,3535
-						y0=y0-3432;
-					}
-				}
 
 				if (x0 >= 0 && x0 <= mpi && y0 >= 0
 				    && y0 <= mpi)