From 9fa61bc68d8fa742abb77b11141f441f3ee24788 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 1 Jun 2017 22:03:50 +0100 Subject: [PATCH] Add resampling code --- inputs.cc | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- inputs.hh | 3 +++ main.cc | 25 ++++++++++++++++--- tiles.hh | 2 +- 4 files changed, 98 insertions(+), 6 deletions(-) diff --git a/inputs.cc b/inputs.cc index e7ef294..c67b871 100644 --- a/inputs.cc +++ b/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); } diff --git a/inputs.hh b/inputs.hh index c946286..30c1e4b 100644 --- a/inputs.hh +++ b/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); diff --git a/main.cc b/main.cc index 7634848..26b5026 100644 --- a/main.cc +++ b/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; diff --git a/tiles.hh b/tiles.hh index 5755e10..2341faa 100644 --- a/tiles.hh +++ b/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 \ No newline at end of file