From 90bd2a3d3208c7038625626879bb3ebc987db165 Mon Sep 17 00:00:00 2001
From: Gareth Evans <gareth.evans@contextis.co.uk>
Date: Fri, 2 Jun 2017 21:18:31 +0100
Subject: [PATCH] Add resize data to specify resolution

---
 inputs.cc | 30 ++++++++++++++++++++++++++++++
 inputs.hh |  1 +
 main.cc   |  3 ++-
 3 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/inputs.cc b/inputs.cc
index 12364a9..a4653b7 100644
--- a/inputs.cc
+++ b/inputs.cc
@@ -9,6 +9,21 @@
 #include "main.hh"
 #include "tiles.hh"
 
+/* Computes the distance between two long/lat points */
+double haversine_formula(double th1, double ph1, double th2, double ph2)
+{
+	#define TO_RAD (3.1415926536 / 180)
+	int R = 6371;
+	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);
+	return asin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * R;
+}
+
 /*
  * resample_data
  * This is used to resample tile data. It is particularly designed for
@@ -84,6 +99,21 @@ int resample_data(int scaling_factor){
 	return 0;
 }
 
+/*
+ * resize_data
+ * This function works in conjuntion with resample_data. It takes a
+ * resolution value in meters as its argument. It then calculates the
+ * nearest (via averaging) resample value and calls resample_data
+ */
+int resize_data(int resolution){
+	double current_res_km = haversine_formula(dem[0].max_north, dem[0].max_west, dem[0].max_north, dem[0].min_west);
+	int current_res = (int) ceil((current_res_km/IPPD)*1000);
+	int scaling_factor = resolution / current_res;
+	if (debug)
+		fprintf(stderr, "Resampling: Current %dm Desired %dm Scale %d\n", current_res, resolution, scaling_factor);
+	return resample_data(scaling_factor);
+}
+
 int loadClutter(char *filename, double radius, struct site tx)
 {
 	/* This function reads a MODIS 17-class clutter file in ASCII Grid format.
diff --git a/inputs.hh b/inputs.hh
index 30c1e4b..2d79a52 100644
--- a/inputs.hh
+++ b/inputs.hh
@@ -5,6 +5,7 @@
 
 /* Resample input tiles to new resolution */
 int resample_data(int scaling_factor);
+int resize_data(int resolution);
 
 int LoadSDF_SDF(char *name, int winfiles);
 int LoadSDF(char *name, int winfiles);
diff --git a/main.cc b/main.cc
index 4a0fbfb..5fc0baa 100644
--- a/main.cc
+++ b/main.cc
@@ -1097,6 +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, "Output:\n");
 		fprintf(stdout,	"     -dbm Plot Rxd signal power instead of field strength\n");
 		fprintf(stdout, "     -rt Rx Threshold (dB / dBm / dBuV/m)\n");
@@ -1705,7 +1706,7 @@ int main(int argc, char *argv[])
 
 		/* If we have been asked to resample the input data; do it now. */
 		if (resample != -1 ){
-			if ((result = resample_data(resample)) != 0) {
+			if ((result = resize_data(resample)) != 0) {
 				fprintf(stderr, "Error resampling data\n");
 				return result;
 			}