From b01b30bba0a3cebdfe86ea6d308f029ddadb8580 Mon Sep 17 00:00:00 2001
From: alex <alex@minty>
Date: Tue, 6 Feb 2018 21:47:52 +0000
Subject: [PATCH] LIDAR edge cases and SUI enhancement

---
 CHANGELOG     |   7 +++
 inputs.cc     | 157 ++++++++++++++++++++++++++++++++++++--------------
 inputs.hh     |   2 +-
 main.cc       |  52 ++++++++---------
 models/sui.cc |   3 +-
 outputs.cc    |  21 +++----
 tiles.cc      |   7 +--
 7 files changed, 160 insertions(+), 89 deletions(-)

diff --git a/CHANGELOG b/CHANGELOG
index 1ae3d1f..046ea10 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -1,5 +1,12 @@
 SIGNAL SERVER CHANGELOG 
 
+3.08 - 17 Dec 2017
+Proper fix for nearfield void
+More LIDAR edge cases: 1x2, 2x1...
+Removed exit() if no colour file
+Polyfilla for gaps at tile edges  
+SUI model clutter factor
+
 3.07 - 2 Sep 2017
 Bugfixes for LIDAR edge cases with mismatched tiles, missing tiles etc.
 
diff --git a/inputs.cc b/inputs.cc
index 0f19cd4..57f9af3 100644
--- a/inputs.cc
+++ b/inputs.cc
@@ -125,6 +125,34 @@ int loadClutter(char *filename, double radius, struct site tx)
 	return 0;
 }
 
+int averageHeight(int height, int width, int x, int y){
+	
+	int total = 0;
+	int c=0;
+	if(dem[0].data[y-1][x-1]>0){
+		total+=dem[0].data[y-1][x-1];
+		c++;
+	}
+	if(dem[0].data[y+1][x+1]>0){
+		total+=dem[0].data[y+1][x+1];
+		c++;
+	}
+	if(dem[0].data[y-1][x+1]>0){
+		total+=dem[0].data[y-1][x+1];
+		c++;
+	}
+	if(dem[0].data[y+1][x-1]>0){
+		total+=dem[0].data[y+1][x-1];
+		c++;
+	}
+
+	if(c>0){
+		return (int)(total/c);
+	}else{
+		return 0;
+	}
+}
+
 int loadLIDAR(char *filenames, int resample)
 {
 	char *filename;
@@ -229,50 +257,77 @@ int loadLIDAR(char *filenames, int resample)
 	/* Now we need to rescale all tiles the the lowest resolution or the requested resolution. ie if we have
 	 * one 1m lidar and one 2m lidar, resize the 2m to fake 1m */
 	float desired_resolution = resample != 0 && smallest_res < resample ? resample : smallest_res;
-	if (desired_resolution > resample && resample != 0 && debug )
-		fprintf(stderr, "Warning: Unable to rescale to requested resolution\n");
-	for (size_t i = 0; i< fc; i++) {
-		float rescale = tiles[i].resolution / (float)desired_resolution;
-		if(debug)
-			fprintf(stderr,"res %.2f desired_res %.2f\n",tiles[i].resolution,(float)desired_resolution);
-		if (rescale != 1){
-			if( (success = tile_rescale(&tiles[i], rescale) != 0 ) ){
-				fprintf(stderr, "Error resampling tiles\n");
-				return success;
-			}
 
-		}
+	if(resample>1){
+		desired_resolution=smallest_res*resample;
 	}
 
+	// Don't resize large 1 deg tiles in large multi-degree plots as it gets messy
+	if(tiles[0].width != 3600){
+	
+		for (size_t i = 0; i< fc; i++) {
+			float rescale = tiles[i].resolution / (float)desired_resolution;
+			if(debug)
+				fprintf(stderr,"res %.5f desired_res %.5f\n",tiles[i].resolution,(float)desired_resolution);
+			if (rescale != 1){
+				if( (success = tile_rescale(&tiles[i], rescale) != 0 ) ){
+					fprintf(stderr, "Error resampling tiles\n");
+					return success;
+				}
+			}
+		}
+
+	}
 	/* Now we work out the size of the giant lidar tile. */
 	double total_width = max_west - min_west >= 0 ? max_west - min_west : max_west + (360 - min_west);
 	double total_height = max_north - min_north;
 	if (debug) {
-		fprintf(stderr,"totalh: %.7f - %.7f = %.7f totalw: %.7f - %.7f = %.7f\n", max_north, min_north, total_height, max_west, min_west, total_width);
+		fprintf(stderr,"totalh: %.7f - %.7f = %.7f totalw: %.7f - %.7f = %.7f fc: %d\n", max_north, min_north, total_height, max_west, min_west, total_width,fc);
 		fprintf(stderr,"mw:%lf Mnw:%lf\n", max_west, min_west);
 	}
 	//detect problematic layouts eg. vertical rectangles
-	//create a synthetic empty tile on the fly with dimensions to make a square
-	/*if(total_height > total_width*1.5){
+	// 1x2
+	if(fc == 2 && desired_resolution < 28 && total_height > total_width){
 		tiles[fc].max_north=max_north;
 		tiles[fc].min_north=min_north;
-		westoffset=westoffset-total_width; // WGS84 for stdout only
-		max_west=max_west+total_width; // Positive westing
+		westoffset=westoffset-(total_height-total_width); // WGS84 for stdout only
+		max_west=max_west+(total_height-total_width); // Positive westing
 		tiles[fc].max_west=max_west; // Positive westing
 		tiles[fc].min_west=max_west;
 		tiles[fc].ppdy=tiles[fc-1].ppdy;
 		tiles[fc].ppdy=tiles[fc-1].ppdx;
-		tiles[fc].width=tiles[fc-1].width;
+		tiles[fc].width=(total_height-total_width);
 		tiles[fc].height=total_height;
 		tiles[fc].data=tiles[fc-1].data;
 		fc++;
-	}*/
 
-	/* This is how we should _theoretically_ work this out, but due to
-	 * the nature of floating point arithmetic and rounding errors, 499 usd to gbpwe need to
-	 * crunch the numbers the hard way */
-	// size_t new_width = total_width * pix_per_deg;
-	// size_t new_height = total_height * pix_per_deg;
+		//calculate deficit
+
+		if (debug) {
+			fprintf(stderr,"deficit: %.4f cellsize: %.9f tiles needed to square: %.1f\n", total_width-total_height,avgCellsize,(total_width-total_height)/avgCellsize);
+		}	
+	}
+	// 2x1
+	if(fc == 2 && desired_resolution < 28 && total_width > total_height){
+		tiles[fc].max_north=max_north+(total_width-total_height);
+		tiles[fc].min_north=max_north;
+		tiles[fc].max_west=max_west; // Positive westing
+		max_north=max_north+(total_width-total_height); // Positive westing
+		tiles[fc].min_west=min_west;
+		tiles[fc].ppdy=tiles[fc-1].ppdy;
+		tiles[fc].ppdy=tiles[fc-1].ppdx;
+		tiles[fc].width=total_width; 
+		tiles[fc].height=(total_width-total_height);
+		tiles[fc].data=tiles[fc-1].data;
+		fc++;
+
+		//calculate deficit
+
+		if (debug) {
+			fprintf(stderr,"deficit: %.4f cellsize: %.9f tiles needed to square: %.1f\n", total_width-total_height,avgCellsize,(total_width-total_height)/avgCellsize);
+		}	
+	
+	}
 
 	size_t new_height = 0;
 	size_t new_width = 0;
@@ -287,12 +342,12 @@ int loadLIDAR(char *filenames, int resample)
 		if ( north_pixel_offset + tiles[i].height > new_height )
 			new_height = north_pixel_offset + tiles[i].height;
 		if (debug)
-			fprintf(stderr,"north_pixel_offset %d west_pixel_offset %d\n", north_pixel_offset, west_pixel_offset);
+			fprintf(stderr,"north_pixel_offset %d west_pixel_offset %d, %d x %d\n", north_pixel_offset, west_pixel_offset,new_height,new_width);
 	}
-	
-	size_t new_tile_alloc = new_width * new_height;
 
+	size_t new_tile_alloc = new_width * new_height;
 	short * new_tile = (short*) calloc( new_tile_alloc, sizeof(short) );
+
 	if ( new_tile == NULL ){
 		free(tiles);
 		return ENOMEM;
@@ -309,12 +364,6 @@ int loadLIDAR(char *filenames, int resample)
 		double west_offset = max_west - tiles[i].max_west >= 0 ? max_west - tiles[i].max_west : max_west + (360 - tiles[i].max_west);
 		size_t north_pixel_offset = north_offset * tiles[i].ppdy;
 		size_t west_pixel_offset = west_offset * tiles[i].ppdx; 
-		/*if(i>0){
-			if(tiles[i].ppdx>tiles[i-1].ppdx+10 || tiles[i].ppdx<tiles[i-1].ppdx-10)
-				west_pixel_offset=prevPixelOffsetW; //PROBLEM WITH DIFF SIZED TILES eg. 4096, 4069..
-			if(tiles[i].ppdy>tiles[i-1].ppdy+10 || tiles[i].ppdy<tiles[i-1].ppdy-10)
-				north_pixel_offset=prevPixelOffsetW; //PROBLEM WITH DIFF SIZED TILES eg. 4096, 4069..
-		}*/
 		prevPixelOffsetW=west_pixel_offset;
 		prevPixelOffsetN=north_pixel_offset;
 
@@ -341,13 +390,18 @@ int loadLIDAR(char *filenames, int resample)
 	// SUPER tile
 	MAXPAGES = 1;
 	IPPD = MAX(new_width,new_height);
+	ippd=IPPD;
+
+	ARRAYSIZE = (MAXPAGES * IPPD) + 10;
+	do_allocs();
+
+	height = new_height;
+	width = new_width;
 
 	if(debug){
-		fprintf(stderr,"Setting IPPD to %d\n",IPPD);
+		fprintf(stderr,"Setting IPPD to %d height %d width %d\n",IPPD,height,width);
 		fflush(stderr);
 	}
-	ARRAYSIZE = (MAXPAGES * IPPD) + 50;
-	do_allocs();
 
 	/* Load the data into the global dem array */
 	dem[0].max_north = max_north;
@@ -357,9 +411,9 @@ int loadLIDAR(char *filenames, int resample)
 	dem[0].max_el = max_elevation;
 	dem[0].min_el = min_elevation;
 
-	/* 
-	 * Copy the lidar tile data into the dem array. The dem array is rotated
-	 * 90 degrees (christ knows why...)
+	/*
+	 * Copy the lidar tile data into the dem array. The dem array is then rotated
+	 * 90 degrees(!)...it's a legacy thing.
 	 */
 	int y = new_height-1;
 	for (size_t h = 0; h < new_height; h++, y--) {
@@ -371,9 +425,17 @@ int loadLIDAR(char *filenames, int resample)
 		}
 	}
 
-	ippd=IPPD;
-	height = new_height;
-	width = new_width;
+	//Polyfilla for warped tiles
+	y = new_height-2;
+	for (size_t h = 0; h < new_height-2; h++, y--) {
+		int x = new_width-2;
+		for (size_t w = 0; w < new_width-2; w++, x--) {
+
+			if(dem[0].data[y][x]<=0){
+				dem[0].data[y][x] = averageHeight(new_height,new_width,x,y);
+			}
+		}
+	}	
 
 	if (debug)
 		fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height);
@@ -1225,7 +1287,7 @@ int LoadLossColors(struct site xmtr)
 
 	/* Default values */
 
-	region.level[0] = 80;
+/*	region.level[0] = 80;
 	region.color[0][0] = 255;
 	region.color[0][1] = 0;
 	region.color[0][2] = 0;
@@ -1304,8 +1366,15 @@ int LoadLossColors(struct site xmtr)
 	region.color[15][0] = 255;
 	region.color[15][1] = 194;
 	region.color[15][2] = 204;
+*/
 
-	region.levels = 16;
+	region.levels = 120; // 240dB max PL
+	for(int i=0; i<region.levels;i++){
+		region.level[i] = i*2;
+		region.color[i][0] = i*2;
+		region.color[i][1] = i*2;
+		region.color[i][2] = i*2;
+	}
 
 	if ( (fd = fopen(filename, "r")) == NULL && xmtr.filename[0] == '\0' )
 	//if ( xmtr.filename[0] == '\0' && (fd = fopen(filename, "r")) == NULL )
diff --git a/inputs.hh b/inputs.hh
index 6d7d353..3ab8599 100644
--- a/inputs.hh
+++ b/inputs.hh
@@ -17,7 +17,7 @@ int LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat);
 int LoadUDT(char *filename);
 int loadLIDAR(char *filename, int resample);
 int loadClutter(char *filename, double radius, struct site tx);
-
+int averageHeight(int h, int w, int x, int y);
 static const char AZ_FILE_SUFFIX[] = ".az";
 static const char EL_FILE_SUFFIX[] = ".el"; 
 
diff --git a/main.cc b/main.cc
index 6e2be69..316645f 100644
--- a/main.cc
+++ b/main.cc
@@ -1,4 +1,4 @@
-double version = 3.07;
+double version = 3.08;
 /****************************************************************************\
 *  Signal Server: Radio propagation simulator by Alex Farrant QCVS, 2E0TDW   *
 ******************************************************************************
@@ -50,10 +50,10 @@ double earthradius, max_range = 0.0, forced_erp, dpp, ppd, yppd,
     min_north = 90, max_north = -90, min_west = 360, max_west = -1, westoffset=180, eastoffset=-180, delta=0, rxGain=0,
     cropLat=-70, cropLon=0;
 
-int ippd, mpi,
+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, resample = 0;
+    0, MAXRAD, hottest = 0, height, width, resample = 0;
 
 unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0;
 
@@ -1097,7 +1097,7 @@ int main(int argc, char *argv[])
 		fprintf(stdout,	"     -tercon Terrain conductivity 0.01-0.0001 (optional)\n");
 		fprintf(stdout, "     -cl Climate code 1-6 (optional)\n");
 		fprintf(stdout, "     -rel Reliability for ITM model 50 to 99 (optional)\n");
-		fprintf(stdout, "     -resample Resample Lidar input to specified resolution in meters (optional)\n");
+		fprintf(stdout, "     -resample Reduce Lidar resolution by specified factor (2 = 50%)\n");
 		fprintf(stdout, "Output:\n");
 		fprintf(stdout,	"     -dbm Plot Rxd signal power instead of field strength\n");
 		fprintf(stdout, "     -rt Rx Threshold (dB / dBm / dBuV/m)\n");
@@ -1619,9 +1619,9 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	if (contour_threshold < -200 || contour_threshold > 200) {
+	if (contour_threshold < -200 || contour_threshold > 240) {
 		fprintf(stderr,
-			"ERROR: Receiver threshold out of range (-200 / +200)");
+			"ERROR: Receiver threshold out of range (-200 / +240)");
 		exit(EINVAL);
 	}
 	if (propmodel > 2 && propmodel < 7 && LR.frq_mhz < 150) {
@@ -1636,6 +1636,11 @@ int main(int argc, char *argv[])
 		exit(EINVAL);
 	}
 
+	if(resample > 10){
+		fprintf(stderr,
+			"ERROR: Cannot resample higher than a factor of 10");
+		exit(EINVAL);	
+	}
 	if (metric) {
 		altitudeLR /= METERS_PER_FOOT;	/* 10ft * 0.3 = 3.3m */
 		max_range /= KM_PER_MILE;	/* 10 / 1.6 = 7.5 */
@@ -1708,14 +1713,14 @@ int main(int argc, char *argv[])
 			exit(result);
 		}
 
-		if(debug){
-			fprintf(stderr,"%.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) (height / (max_north-min_north));
+		yppd=(double) (width / (max_west-min_west));
+
+		if(debug){
+			fprintf(stderr,"ppd %lf, yppd %lf, %.4f,%.4f,%.4f,%.4f,%d x %d\n",ppd,yppd,max_north,min_west,min_north,max_west,width,height);
+		}
+
 
-		//ppd=(double)ippd;
-		//yppd=ppd; 
 
 
 		if(delta>0){
@@ -1853,7 +1858,7 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	if(max_range>100){
+	if(max_range>100 || LR.frq_mhz==446.446){
 		cropping=false;
 	}
 	if (ppa == 0) {
@@ -1869,21 +1874,16 @@ int main(int argc, char *argv[])
                         if(debug)
                         	fprintf(stderr,"Finished PlotPropagation()\n");
 
-			if(!lidar){
-				if (LR.erp == 0.0)
-					hottest=9; // 9dB nearfield
 
-				// nearfield bugfix
-				for (lat = tx_site[0].lat - 0.0005;
-					 lat <= tx_site[0].lat + 0.0005;
-					 lat = lat + 0.0005) {
-					for (lon = tx_site[0].lon - 0.0005;
-						 lon <= tx_site[0].lon + 0.0005;
-						 lon = lon + 0.0005) {
-						PutSignal(lat, lon, hottest);
+				// nearfield void
+
+				for (float x=-0.001; x<0.001;x=x+0.0001){
+					for (float y=-0.001; y<0.001;y=y+0.0001){
+						if(GetSignal(tx_site[0].lat+y, tx_site[0].lon+x)<=0){
+							PutSignal(tx_site[0].lat+y, tx_site[0].lon+x, hottest);
+						}
 					}
 				}
-			}
 
 			if(cropping){
 				// CROPPING. croplat assigned in propPathLoss()
diff --git a/models/sui.cc b/models/sui.cc
index 223db65..a1b969b 100644
--- a/models/sui.cc
+++ b/models/sui.cc
@@ -15,13 +15,12 @@ double SUIpathLoss(double f, double TxH, double RxH, double d, int mode)
 	   http://www.cl.cam.ac.uk/research/dtg/lce-pub/public/vsa23/VTC05_Empirical.pdf
 	 */
 	d = d * 1000.0;		// km to m
-	RxH = RxH * 1000.0; // Correction factor for CPE units.
 
 	// Urban (A1) is default
 	double a = 4.6;
 	double b = 0.0075;
 	double c = 12.6;
-	double s = 0.0; // Optional fading value. Max 10.6dB
+	double s = 10.6; // Optional fading value. Max 10.6dB
 	double XhCF = -10.8;
 
 	if (mode == 2) { // Suburban
diff --git a/outputs.cc b/outputs.cc
index 5509223..6fdce53 100644
--- a/outputs.cc
+++ b/outputs.cc
@@ -46,7 +46,7 @@ void DoPathLoss(char *filename, unsigned char geo, unsigned char kml,
 
 	if( (success = LoadLossColors(xmtr[0])) != 0 ){
 		fprintf(stderr,"Error loading loss colors\n");
-		exit(success);
+		//exit(success);
 	}
 
 	if( filename != NULL ) {
@@ -289,7 +289,7 @@ int DoSigStr(char *filename, unsigned char geo, unsigned char kml,
 
 	if( (success = LoadSignalColors(xmtr[0])) != 0 ){
 		fprintf(stderr,"Error loading signal colors\n");
-		return success;
+		//exit(success);
 	}
 
 	if( filename != NULL ) {
@@ -540,7 +540,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml,
 
 	if( (success = LoadDBMColors(xmtr[0])) != 0 ){
 		fprintf(stderr,"Error loading DBM colors\n");
-		exit(success);
+		//exit(success);
 	}
 
 	if( filename != NULL ) {
@@ -1070,7 +1070,6 @@ void PathReport(struct site source, struct site destination, char *name,
 	angle1 = ElevationAngle(source, destination);
 	angle2 = ElevationAngle2(source, destination, earthradius);
 
-	if (got_azimuth_pattern || got_elevation_pattern) {
 		x = (int)rint(10.0 * (10.0 - angle2));
 
 		if (x >= 0 && x <= 1000)
@@ -1078,7 +1077,6 @@ void PathReport(struct site source, struct site destination, char *name,
 			    (double)LR.antenna_pattern[(int)rint(azimuth)][x];
 
 		patterndB = 20.0 * log10(pattern);
-	}
 
 	if (metric)
 		fprintf(fd2, "Distance to %s: %.2f kilometers\n",
@@ -1307,9 +1305,7 @@ void PathReport(struct site source, struct site destination, char *name,
 		fprintf(fd2, "\nSummary for the link between %s and %s:\n\n",
 			source.name, destination.name);
 
-		if (patterndB != 0.0)
-			fprintf(fd2,
-				"%s antenna pattern towards %s: %.3f (%.2f dB)\n",
+		fprintf(fd2, "%s antenna pattern towards %s: %.3f (%.2f dB)\n",
 				source.name, destination.name, pattern,
 				patterndB);
 
@@ -1508,8 +1504,11 @@ void PathReport(struct site source, struct site destination, char *name,
 				pattern =
 				    (double)LR.antenna_pattern[(int)azimuth][x];
 
-				if (pattern != 0.0)
+				if (pattern != 0.0){
 					patterndB = 20.0 * log10(pattern);
+				}else{
+					patterndB = 0.0;
+				}
 			}
 
 			else
@@ -1548,9 +1547,7 @@ void PathReport(struct site source, struct site destination, char *name,
 				"Attenuation due to terrain shielding: %.2f dB\n",
 				loss - free_space_loss);
 
-		if (patterndB != 0.0)
-			fprintf(fd2,
-				"Total path loss including %s antenna pattern: %.2f dB\n",
+		fprintf(fd2,"Total path loss including %s antenna pattern: %.2f dB\n",
 				source.name, total_loss);
 
 		if (LR.erp != 0.0) {
diff --git a/tiles.cc b/tiles.cc
index d946a18..186794d 100644
--- a/tiles.cc
+++ b/tiles.cc
@@ -6,7 +6,7 @@
 #include "tiles.hh"
 #include "common.h"
 
-#define MAX_LINE 25000
+#define MAX_LINE 50000
 
 /* Computes the distance between two long/lat points */
 double haversine_formula(double th1, double ph1, double th2, double ph2)
@@ -16,7 +16,6 @@ double haversine_formula(double th1, double ph1, double th2, double ph2)
 	double dx, dy, dz;
 	ph1 -= ph2;
 	ph1 *= TO_RAD, th1 *= TO_RAD, th2 *= TO_RAD;
- 
 	dz = sin(th1) - sin(th2);
 	dx = cos(ph1) * cos(th1) - cos(th2);
 	dy = sin(ph1) * cos(th1);
@@ -120,7 +119,7 @@ int tile_load_lidar(tile_t *tile, char *filename){
 	tile->precise_resolution = (current_res_km/MAX(tile->width,tile->height)*1000);
 
 	// Round to nearest 0.5
-	tile->resolution = tile->precise_resolution < 0.5f ? 0.5f : floor((tile->precise_resolution * 2)+0.5) / 2;
+	tile->resolution = tile->precise_resolution < 0.5f ? 0.5f : ceil((tile->precise_resolution * 2)+0.5) / 2;
 
 	tile->width_deg = tile->max_west - tile->min_west >= 0 ? tile->max_west - tile->min_west : tile->max_west + (360 - tile->min_west);
 	tile->height_deg = tile->max_north - tile->min_north;
@@ -235,7 +234,7 @@ int tile_resize(tile_t* tile, int resolution){
 	int current_res = (int) ceil((current_res_km/IPPD)*1000);
 	float scaling_factor = resolution / current_res;
 	if (debug)
-		fprintf(stderr, "Resampling: Current %dm Desired %dm Scale %d\n", current_res, resolution, scaling_factor);
+		fprintf(stderr, "Resampling: Current %dm Desired %dm Scale %.1f\n", current_res, resolution, scaling_factor);
 	return tile_rescale(tile, scaling_factor);
 }