From ee67a86a9f9855341c2d3a9797697d46883f7237 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 1 Jun 2017 18:54:15 +0100 Subject: [PATCH 01/16] Improve error message on image output code --- outputs.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/outputs.cc b/outputs.cc index 5d79961..d132849 100644 --- a/outputs.cc +++ b/outputs.cc @@ -35,7 +35,7 @@ void DoPathLoss(char *filename, unsigned char geo, unsigned char kml, int success; if( (success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0 ){ - fprintf(stderr,"Error initializing image\n"); + fprintf(stderr,"Error initializing image: %s\n", strerror(success)); exit(success); } @@ -278,7 +278,7 @@ int DoSigStr(char *filename, unsigned char geo, unsigned char kml, int success; if((success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0){ - fprintf(stderr,"Error initializing image\n"); + fprintf(stderr,"Error initializing image: %s\n", strerror(success)); exit(success); } @@ -529,7 +529,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, int success; if( (success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0 ){ - fprintf(stderr,"Error initializing image\n"); + fprintf(stderr,"Error initializing image: %s\n", strerror(success)); exit(success); } @@ -775,7 +775,7 @@ void DoLOS(char *filename, unsigned char geo, unsigned char kml, int success; if((success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0){ - fprintf(stderr,"Error initializing image\n"); + fprintf(stderr,"Error initializing image: %s\n", strerror(success)); exit(success); } From 83cdc915d254ae01449cbb18d8c42ed0fd26bfcb Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 1 Jun 2017 18:54:49 +0100 Subject: [PATCH 02/16] Add support for working out IPPD based on LIDAR tiles --- inputs.cc | 240 +++++++++++++++++++++++++++++++----------------------- tiles.cc | 113 +++++++++++++++++++++++++ tiles.hh | 28 +++++++ 3 files changed, 277 insertions(+), 104 deletions(-) create mode 100644 tiles.cc create mode 100644 tiles.hh diff --git a/inputs.cc b/inputs.cc index 1f45304..e52f94a 100644 --- a/inputs.cc +++ b/inputs.cc @@ -7,6 +7,7 @@ #include #include "common.h" #include "main.hh" +#include "tiles.hh" int loadClutter(char *filename, double radius, struct site tx) @@ -215,15 +216,15 @@ int loadLIDAR(char *filenames) { char *filename; char *files[100]; // 10x10 tiles - int x, y, indx = 0, fc = 0, hoffset = 0, voffset = 0, pos, - dem_alloced = 0; - double xll, yll, xur, yur, cellsize, avgCellsize; + int indx = 0, fc = 0, hoffset = 0, voffset = 0, pos, success; + double xll, yll, xur, yur, cellsize, avgCellsize = 0, smCellsize = 0; char found, free_page = 0, jline[20], lid_file[255], - path_plus_name[255], *junk = NULL; + path_plus_name[255], *junk = NULL; char line[25000]; char * pch; - double TO_DEG = (180 / PI); + double TO_DEG = (180 / PI); FILE *fd; + tile_t *tiles; // test for multiple files filename = strtok(filenames, " ,"); @@ -233,110 +234,141 @@ int loadLIDAR(char *filenames) fc++; } - while (indx < fc) { - if( (fd = fopen(files[indx], "rb")) == NULL ) + /* Allocate the tile array */ + tiles = (tile_t*)calloc(fc,sizeof(tile_t)); + + /* Load each tile in turn */ + for (indx = 0; indx < fc; indx++) { + + /* Grab the tile metadata */ + if( (success = tile_load_lidar(&tiles[indx], files[indx])) != 0 ){ + fprintf(stderr,"Failed to load LIDAR tile %s\n",files[indx]); + fflush(stderr); + free(tiles); + return success; + } + + if (debug) { + fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, tiles[indx].width); + fflush(stderr); + } + + // Increase the average cell size + avgCellsize += tiles[indx].cellsize; + // Update the smallest cell size + if (smCellsize == 0 || tiles[indx].cellsize < smCellsize) { + smCellsize = tiles[indx].cellsize; + } + + } + + /* Iterate through all tiles to find the largest x/y dimension + to use as the global IPPD value. We do this here so that we are + sure to allocate enough memory for the remainder of the calculations */ + int dimension_max = -1; + for (size_t i = 0; i < fc; i++) { + if( tiles[i].width > tiles[i].height && tiles[i].width > dimension_max ) + dimension_max = tiles[i].width; + else if( tiles[i].height > dimension_max ) + dimension_max = tiles[i].height; + } + + MAXPAGES = fc; + IPPD = dimension_max; + if(debug){ + fprintf(stderr,"Setting IPPD to %d\n",IPPD); + fflush(stderr); + } + // add fudge as reprojected tiles sometimes vary by a pixel or ten + IPPD += 50; + ARRAYSIZE = (MAXPAGES * IPPD) + 50; + do_allocs(); + // reset the IPPD after allocations + IPPD -= 50; + + /* Load the data into the global dem array */ + // for (size_t indx = 0; indx < fc; indx++) { + // dem[indx].max_north = tiles[indx].yur; + // dem[indx].min_west = tiles[indx].xur; + // dem[indx].min_north = tiles[indx].yll; + // dem[indx].max_west = tiles[indx].xll; + // dem[indx].max_el = tiles[indx].max_el; + // dem[indx].min_el = tiles[indx].min_el; + + // //Not sure why we set these here... + // if (tiles[indx].max_el > max_elevation) + // max_elevation = tiles[indx].max_el; + // if (tiles[indx].min_el < min_elevation) + // min_elevation = tiles[indx].min_el; + + // 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) + // min_west = dem[indx].min_west; + + // if (max_west == -1) { + // max_west = dem[indx].max_west; + // } else { + // if (abs(dem[indx].max_west - max_west) < 180) { + // if (dem[indx].max_west > max_west) + // max_west = dem[indx].max_west; + // } else { + // if (dem[indx].max_west < max_west) + // max_west = dem[indx].max_west; + // } + // } + + // if (min_west == 360) { + // min_west = dem[indx].min_west; + // } else { + // if (fabs(dem[indx].min_west - min_west) < 180.0) { + // if (dem[indx].min_west < min_west) + // min_west = dem[indx].min_west; + // } else { + // if (dem[indx].min_west > min_west) + // min_west = dem[indx].min_west; + // } + // } + + // // /* Copy the lidar tile data into the dem array */ + // // int x = tiles[indx].width-1; + // // int y = tiles[indx].height-1; + // // for (size_t h = 0; h < tiles[indx].height; h++, y--) { + // // for (size_t w = 0; w < tiles[indx].width; w++, x--) { + // // dem[indx].data[y][x] = tiles[indx].data[h*tiles[indx].width + w]; + // // dem[indx].signal[h][w] = 0; + // // dem[indx].mask[h][w] = 0; + // // } + // // } + + // } + + for ( int indx = 0; indx < fc; indx++ ) { + if ( (fd = fopen(files[indx],"rb")) == NULL) return errno; - - if (fgets(line, 255, fd) != NULL) { - pch = strtok (line," "); - pch = strtok (NULL, " "); - width = atoi(pch); // ncols - - if (debug) { - fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, width); - fflush(stderr); - } - - if (fgets(line, 255, fd) != NULL) - height = atoi(pch); // nrows - - if (!dem_alloced) { - //Reduce MAXPAGES to increase speed - MAXPAGES=fc; - if(width>height){ - IPPD = width; - }else{ - IPPD = height; - } - // add fudge as reprojected tiles sometimes vary by a pixel or ten - IPPD+=50; - ARRAYSIZE = (MAXPAGES * IPPD)+50; - do_allocs(); - dem_alloced = 1; - } - - } - if (fgets(line, 255, fd) != NULL) { - sscanf(pch, "%lf", &xll); // xll - } - - if (fgets(line, 255, fd) != NULL) { - sscanf(pch, "%lf", &yll); // yll - } - - if (fgets(line, 255, fd) != NULL) - sscanf(pch, "%lf", &cellsize); - - avgCellsize=avgCellsize+cellsize; - - /*if(cellsize>=0.5){ // 50cm LIDAR? - // compute xur and yur with inverse haversine if cellsize in *metres* - 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); - //} - - if (xur > eastoffset) - eastoffset = xur; - if (xll < westoffset) - westoffset = xll; - - if (debug) - fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur); - - - // Greenwich straddling hack - if (xll <= 0 && xur > 0) { - xll = (xur - xll); // full width - xur = 0.0; // budge it along so it's west of greenwich - delta = eastoffset; // add to Tx longitude later - } else { - // Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;) - if (xll >= 0) - xll = 360-xll; - if(xur >= 0) - xur = 360-xur; - if(xll < 0) - xll = xll * -1; - if(xur < 0) - xur = xur * -1; - } - if (debug) - fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta); - - - fgets(line, 255, fd); // NODATA - pos = ftell(fd); - - // tile 0 [x| ] + 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); + readLIDAR(fd, tiles[indx].height, tiles[indx].width, indx, tiles[indx].yur, tiles[indx].xur, tiles[indx].yll, tiles[indx].xll); + } - readLIDAR(fd, height, width, indx, yur, xur, yll, xll); - - fclose(fd); - if (debug) - fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); - indx++; - } // filename(s) - IPPD=width; ippd=IPPD; - height = (unsigned)((max_north-min_north) / cellsize); - width = (unsigned)((max_west-min_west) / cellsize); + height = (unsigned)((max_north-min_north) / smCellsize); + width = (unsigned)((max_west-min_west) / smCellsize); + + fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); if (debug) fprintf(stderr, "fc %d WIDTH %d HEIGHT %d ippd %d minN %.5f maxN %.5f minW %.5f maxW %.5f avgCellsize %.5f\n", fc, width, height, ippd,min_north,max_north,min_west,max_west,avgCellsize); diff --git a/tiles.cc b/tiles.cc new file mode 100644 index 0000000..edcdd77 --- /dev/null +++ b/tiles.cc @@ -0,0 +1,113 @@ +#include +#include +#include +#include +#include "tiles.hh" +#include "common.h" + +#define MAX_LINE 25000 + +int tile_load_lidar(tile_t *tile, char *filename){ + FILE *fd; + char line[MAX_LINE]; + int nextval; + char *pch; + int w, x, y; + + /* Clear the tile data */ + memset(tile,0x00,sizeof(tile_t)); + + /* Open the file handle and return on error */ + if ( (fd = fopen(filename,"r")) == NULL ) + return errno; + + /* This is where we read the header data */ + /* The string is split for readability but is parsed as a block */ + if( fscanf(fd,"%*s %d\n" "%*s %d\n" "%*s %lf\n" "%*s %lf\n" "%*s %lf\n" "%*s %d\n",&tile->width,&tile->height,&tile->xll,&tile->yll,&tile->cellsize,&tile->nodata) != 6 ){ + fclose(fd); + return -1; + } + + tile->datastart = ftell(fd); + + if(debug){ + fprintf(stderr,"w:%d h:%d s:%lf\n", tile->width, tile->height, tile->cellsize); + fflush(stderr); + } + + /* Set the filename */ + tile->filename = strdup(filename); + + /* Perform xur calcs */ + // Degrees with GDAL option: -co "FORCE_CELLSIZE=YES" + tile->xur = tile->xll+(tile->cellsize*tile->width); + tile->yur = tile->yll+(tile->cellsize*tile->height); + + if (tile->xur > eastoffset) + eastoffset = tile->xur; + if (tile->xll < westoffset) + westoffset = tile->xll; + + // if (debug) + // fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur); + + // Greenwich straddling hack + if (tile->xll <= 0 && tile->xur > 0) { + tile->xll = (tile->xur - tile->xll); // full width + tile->xur = 0.0; // budge it along so it's west of greenwich + delta = eastoffset; // add to Tx longitude later + } else { + // Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;) + if (tile->xll >= 0) + tile->xll = 360-tile->xll; + if(tile->xur >= 0) + tile->xur = 360-tile->xur; + if(tile->xll < 0) + tile->xll = tile->xll * -1; + if(tile->xur < 0) + tile->xur = tile->xur * -1; + } + + if (debug) + fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", tile->yll, tile->yur, tile->xur, tile->xll, delta); + +#ifdef _LOAD_TILEDATA + /* Read the actual tile data */ + /* Allocate the array for the lidar data */ + if ( (tile->data = (int*) calloc(tile->width * tile->height, sizeof(int))) == NULL ) { + fclose(fd); + free(tile->filename); + return ENOMEM; + } + + size_t loaded = 0; + for (size_t h = 0; h < tile->height; h++) { + if (fgets(line, MAX_LINE, fd) != NULL) { + pch = strtok(line, " "); // split line into values + for (w = 0; w < tile->width && pch != NULL; w++) { + /* If the data is less than a *magic* minimum, normalize it to zero */ + nextval = atoi(pch); + if (nextval <= -9999) + nextval = 0; + tile->data[h*tile->width + w] = nextval; + loaded++; + if ( nextval > tile->max_el ) + tile->max_el = nextval; + if ( nextval < tile->min_el ) + tile->min_el = nextval; + pch = strtok(NULL, " "); + }//while + } else { + fprintf(stderr, "LIDAR error @ x %d y %d file %s\n", x, y, filename); + }//if + } + + if (debug) + fprintf(stderr,"Pixels loaded: %zu/%zu\n",loaded,tile->width*tile->height); +#endif + + /* All done, close the LIDAR file */ + fclose(fd); + + return 0; +} \ No newline at end of file diff --git a/tiles.hh b/tiles.hh new file mode 100644 index 0000000..5755e10 --- /dev/null +++ b/tiles.hh @@ -0,0 +1,28 @@ +#ifndef _TILES_HH_ +#define _TILES_HH_ + +typedef struct _tile_t{ + char *filename; + union{ + int cols; + int width; + }; + union{ + int rows; + int height; + }; + double xll; + double yll; + double xur; + double yur; + double cellsize; + long long datastart; + int nodata; + int max_el; + int min_el; + int *data; +} tile_t, *ptile_t; + +int tile_load_lidar(tile_t*, char *, int); + +#endif \ No newline at end of file From d2919ece4d0f9e8a9619ad6a6092863ebbd0cc27 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 1 Jun 2017 18:56:01 +0100 Subject: [PATCH 03/16] Clear out dead tile loading code --- inputs.cc | 69 ------------------------------------------------------- 1 file changed, 69 deletions(-) diff --git a/inputs.cc b/inputs.cc index e52f94a..e7ef294 100644 --- a/inputs.cc +++ b/inputs.cc @@ -286,75 +286,6 @@ int loadLIDAR(char *filenames) // reset the IPPD after allocations IPPD -= 50; - /* Load the data into the global dem array */ - // for (size_t indx = 0; indx < fc; indx++) { - // dem[indx].max_north = tiles[indx].yur; - // dem[indx].min_west = tiles[indx].xur; - // dem[indx].min_north = tiles[indx].yll; - // dem[indx].max_west = tiles[indx].xll; - // dem[indx].max_el = tiles[indx].max_el; - // dem[indx].min_el = tiles[indx].min_el; - - // //Not sure why we set these here... - // if (tiles[indx].max_el > max_elevation) - // max_elevation = tiles[indx].max_el; - // if (tiles[indx].min_el < min_elevation) - // min_elevation = tiles[indx].min_el; - - // 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) - // min_west = dem[indx].min_west; - - // if (max_west == -1) { - // max_west = dem[indx].max_west; - // } else { - // if (abs(dem[indx].max_west - max_west) < 180) { - // if (dem[indx].max_west > max_west) - // max_west = dem[indx].max_west; - // } else { - // if (dem[indx].max_west < max_west) - // max_west = dem[indx].max_west; - // } - // } - - // if (min_west == 360) { - // min_west = dem[indx].min_west; - // } else { - // if (fabs(dem[indx].min_west - min_west) < 180.0) { - // if (dem[indx].min_west < min_west) - // min_west = dem[indx].min_west; - // } else { - // if (dem[indx].min_west > min_west) - // min_west = dem[indx].min_west; - // } - // } - - // // /* Copy the lidar tile data into the dem array */ - // // int x = tiles[indx].width-1; - // // int y = tiles[indx].height-1; - // // for (size_t h = 0; h < tiles[indx].height; h++, y--) { - // // for (size_t w = 0; w < tiles[indx].width; w++, x--) { - // // dem[indx].data[y][x] = tiles[indx].data[h*tiles[indx].width + w]; - // // dem[indx].signal[h][w] = 0; - // // dem[indx].mask[h][w] = 0; - // // } - // // } - - // } - for ( int indx = 0; indx < fc; indx++ ) { if ( (fd = fopen(files[indx],"rb")) == NULL) return errno; From 9fa61bc68d8fa742abb77b11141f441f3ee24788 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 1 Jun 2017 22:03:50 +0100 Subject: [PATCH 04/16] 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 From 67c145dd7d2290c7f0c7038aae99bddc02a7d7ef Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 1 Jun 2017 22:58:06 +0100 Subject: [PATCH 05/16] Fix height width and makefule issues --- Makefile | 4 +++- inputs.cc | 2 ++ main.cc | 38 +++++++++++++++++++------------------- 3 files changed, 24 insertions(+), 20 deletions(-) diff --git a/Makefile b/Makefile index 23ef504..56926bb 100644 --- a/Makefile +++ b/Makefile @@ -8,7 +8,7 @@ LIBS = -lm -lpthread -ldl VPATH = models objects = main.o cost.o ecc33.o ericsson.o fspl.o hata.o itwom3.0.o \ - los.o sui.o pel.o egli.o inputs.o outputs.o image.o image-ppm.o + los.o sui.o pel.o egli.o inputs.o outputs.o image.o image-ppm.o tiles.o GCC_MAJOR := $(shell $(CXX) -dumpversion 2>&1 | cut -d . -f 1) GCC_MINOR := $(shell $(CXX) -dumpversion 2>&1 | cut -d . -f 2) @@ -52,6 +52,8 @@ image-ppm.o: image-ppm.cc los.o: los.cc common.h main.hh cost.hh ecc33.hh ericsson.hh fspl.hh hata.hh \ itwom3.0.hh sui.hh pel.hh egli.hh +tiles.o: tiles.cc tiles.hh common.h + .PHONY: clean clean: rm -f $(objects) signalserver signalserverHD signalserverLIDAR diff --git a/inputs.cc b/inputs.cc index c67b871..12364a9 100644 --- a/inputs.cc +++ b/inputs.cc @@ -76,6 +76,8 @@ int resample_data(int scaling_factor){ 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 */ + height /= scaling_factor; + width /= scaling_factor; IPPD = new_ippd; ippd = IPPD; ARRAYSIZE = (MAXPAGES * IPPD) + 50; diff --git a/main.cc b/main.cc index 26b5026..4a0fbfb 100644 --- a/main.cc +++ b/main.cc @@ -1884,26 +1884,26 @@ int main(int argc, char *argv[]) } } - if(max_range<=100){ - // CROPPING. croplat assigned in propPathLoss() - max_north=cropLat; // MAX(path.lat[y]) - // Edge case #1 - EAST/WEST - if(cropLon>357 && tx_site[0].lon < 3) - cropLon=tx_site[0].lon+3; - // Edge case #2 - EAST/EAST - if(cropLon>359.5 && tx_site[0].lon > 359.5) - cropLon=362; - max_west=cropLon; // MAX(path.lon[y]) - cropLat-=tx_site[0].lat; // angle from tx to edge - cropLon-=tx_site[0].lon; - width=(int)((cropLon*ppd)*2); - height=(int)((cropLat*ppd)*2); + // if(max_range<=100){ + // // CROPPING. croplat assigned in propPathLoss() + // max_north=cropLat; // MAX(path.lat[y]) + // // Edge case #1 - EAST/WEST + // if(cropLon>357 && tx_site[0].lon < 3) + // cropLon=tx_site[0].lon+3; + // // Edge case #2 - EAST/EAST + // if(cropLon>359.5 && tx_site[0].lon > 359.5) + // cropLon=362; + // max_west=cropLon; // MAX(path.lon[y]) + // cropLat-=tx_site[0].lat; // angle from tx to edge + // cropLon-=tx_site[0].lon; + // width=(int)((cropLon*ppd)*2); + // height=(int)((cropLat*ppd)*2); - if(width>3600*10){ - fprintf(stderr,"FATAL BOUNDS! max_west: %.4f cropLat: %.4f cropLon: %.4f longitude: %.5f\n",max_west,cropLat,cropLon,tx_site[0].lon); - return 0; - } - } + // if(width>3600*10){ + // fprintf(stderr,"FATAL BOUNDS! max_west: %.4f cropLat: %.4f cropLon: %.4f longitude: %.5f\n",max_west,cropLat,cropLon,tx_site[0].lon); + // return 0; + // } + // } // Write bitmap if (LR.erp == 0.0) From 90bd2a3d3208c7038625626879bb3ebc987db165 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Fri, 2 Jun 2017 21:18:31 +0100 Subject: [PATCH 06/16] 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; } From fbeffe2e9806d85872b3591008020679ab7f872c Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Mon, 5 Jun 2017 23:09:29 +0100 Subject: [PATCH 07/16] Add debug switch --- inputs.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/inputs.cc b/inputs.cc index a4653b7..6a26f16 100644 --- a/inputs.cc +++ b/inputs.cc @@ -88,7 +88,8 @@ int resample_data(int scaling_factor){ } /* Report to the user */ - fprintf(stderr, "Resampling IPPD %d->%d min/max el %d/%d\n", IPPD, new_ippd, min_elevation, max_elevation); + if (debug) + 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 */ height /= scaling_factor; @@ -403,7 +404,8 @@ int loadLIDAR(char *filenames) height = (unsigned)((max_north-min_north) / smCellsize); width = (unsigned)((max_west-min_west) / smCellsize); - fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); + if (debug) + fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); if (debug) fprintf(stderr, "fc %d WIDTH %d HEIGHT %d ippd %d minN %.5f maxN %.5f minW %.5f maxW %.5f avgCellsize %.5f\n", fc, width, height, ippd,min_north,max_north,min_west,max_west,avgCellsize); From 580fc804e2c116bb0715d7491f0ac271da7740d6 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Tue, 6 Jun 2017 20:54:39 +0100 Subject: [PATCH 08/16] Move tile loading into tiles.cc --- inputs.cc | 162 +++++++++++++++++++++++------------------------------- tiles.cc | 9 +-- 2 files changed, 72 insertions(+), 99 deletions(-) diff --git a/inputs.cc b/inputs.cc index 6a26f16..1d8652b 100644 --- a/inputs.cc +++ b/inputs.cc @@ -231,92 +231,6 @@ int loadClutter(char *filename, double radius, struct site tx) return 0; } -void readLIDAR(FILE *fd, int h, int w, int indx,double n, double e, double s, double west) -{ - int x = 0, y = 0, reads = 0, a=0, b=0, avg=0, tWidth = 0, tHeight = 0; - char line[25000]; - char *pch; - - 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) - min_west = dem[indx].min_west; - - if (max_west == -1) { - max_west = dem[indx].max_west; - } else { - if (abs(dem[indx].max_west - max_west) < 180) { - if (dem[indx].max_west > max_west) - max_west = dem[indx].max_west; - } else { - if (dem[indx].max_west < max_west) - max_west = dem[indx].max_west; - } - } - - if (min_west == 360) { - min_west = dem[indx].min_west; - } else { - if (fabs(dem[indx].min_west - min_west) < 180.0) { - if (dem[indx].min_west < min_west) - min_west = dem[indx].min_west; - } else { - if (dem[indx].min_west > min_west) - min_west = dem[indx].min_west; - } - } - - for (y = h-1; y > -1; y--) { - x = w-1; - - if (fgets(line, 25000, fd) != NULL) { - pch = strtok(line, " "); // split line into values - - while (pch != NULL && x > -1) { - if (atoi(pch) <= -9999) - pch = "0"; - dem[indx].data[y][x] = atoi(pch); - dem[indx].signal[x][y] = 0; - dem[indx].mask[x][y] = 0; - if (atoi(pch) > dem[indx].max_el) { - dem[indx].max_el = atoi(pch); - max_elevation = atoi(pch); - } - if (atoi(pch) < dem[indx].min_el) { - dem[indx].min_el = atoi(pch); - min_elevation = dem[indx].min_el; - } - - x--; - pch = strtok(NULL, " "); - }//while - - - } else { - fprintf(stderr, "LIDAR error @ x %d y %d indx %d\n", - x, y, indx); - }//if - }//for - -} - int loadLIDAR(char *filenames) { char *filename; @@ -391,13 +305,75 @@ int loadLIDAR(char *filenames) // reset the IPPD after allocations IPPD -= 50; - for ( int indx = 0; indx < fc; indx++ ) { - if ( (fd = fopen(files[indx],"rb")) == NULL) - return errno; - fseek(fd, tiles[indx].datastart, SEEK_SET); - if (debug) - 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); + /* Load the data into the global dem array */ + for (size_t indx = 0; indx < fc; indx++) { + dem[indx].max_north = tiles[indx].yur; + dem[indx].min_west = tiles[indx].xur; + dem[indx].min_north = tiles[indx].yll; + dem[indx].max_west = tiles[indx].xll; + dem[indx].max_el = tiles[indx].max_el; + dem[indx].min_el = tiles[indx].min_el; + + if (tiles[indx].max_el > max_elevation) + max_elevation = tiles[indx].max_el; + if (tiles[indx].min_el < min_elevation) + min_elevation = tiles[indx].min_el; + + 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) + min_west = dem[indx].min_west; + + if (max_west == -1) { + max_west = dem[indx].max_west; + } else { + if (abs(dem[indx].max_west - max_west) < 180) { + if (dem[indx].max_west > max_west) + max_west = dem[indx].max_west; + } else { + if (dem[indx].max_west < max_west) + max_west = dem[indx].max_west; + } + } + + if (min_west == 360) { + min_west = dem[indx].min_west; + } else { + if (fabs(dem[indx].min_west - min_west) < 180.0) { + if (dem[indx].min_west < min_west) + min_west = dem[indx].min_west; + } else { + if (dem[indx].min_west > min_west) + min_west = dem[indx].min_west; + } + } + + /* + * Copy the lidar tile data into the dem array. The dem array is rotated + * 90 degrees (christ knows why...) + */ + int y = tiles[indx].height-1; + for (size_t h = 0; h < tiles[indx].height; h++, y--) { + int x = tiles[indx].width-1; + for (size_t w = 0; w < tiles[indx].width; w++, x--) { + dem[indx].data[y][x] = tiles[indx].data[h*tiles[indx].width + w]; + dem[indx].signal[y][x] = 0; + dem[indx].mask[y][x] = 0; + } + } + } ippd=IPPD; diff --git a/tiles.cc b/tiles.cc index edcdd77..e5fb66d 100644 --- a/tiles.cc +++ b/tiles.cc @@ -12,7 +12,6 @@ int tile_load_lidar(tile_t *tile, char *filename){ char line[MAX_LINE]; int nextval; char *pch; - int w, x, y; /* Clear the tile data */ memset(tile,0x00,sizeof(tile_t)); @@ -71,7 +70,6 @@ int tile_load_lidar(tile_t *tile, char *filename){ if (debug) fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", tile->yll, tile->yur, tile->xur, tile->xll, delta); -#ifdef _LOAD_TILEDATA /* Read the actual tile data */ /* Allocate the array for the lidar data */ if ( (tile->data = (int*) calloc(tile->width * tile->height, sizeof(int))) == NULL ) { @@ -84,7 +82,7 @@ int tile_load_lidar(tile_t *tile, char *filename){ for (size_t h = 0; h < tile->height; h++) { if (fgets(line, MAX_LINE, fd) != NULL) { pch = strtok(line, " "); // split line into values - for (w = 0; w < tile->width && pch != NULL; w++) { + for (size_t w = 0; w < tile->width && pch != NULL; w++) { /* If the data is less than a *magic* minimum, normalize it to zero */ nextval = atoi(pch); if (nextval <= -9999) @@ -98,13 +96,12 @@ int tile_load_lidar(tile_t *tile, char *filename){ pch = strtok(NULL, " "); }//while } else { - fprintf(stderr, "LIDAR error @ x %d y %d file %s\n", x, y, filename); + fprintf(stderr, "LIDAR error @ h %zu file %s\n", h, filename); }//if } if (debug) - fprintf(stderr,"Pixels loaded: %zu/%zu\n",loaded,tile->width*tile->height); -#endif + fprintf(stderr,"Pixels loaded: %zu/%d\n",loaded,tile->width*tile->height); /* All done, close the LIDAR file */ fclose(fd); From c6873f8dead0d03b606134f6d52549004968a77a Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Tue, 6 Jun 2017 21:47:05 +0100 Subject: [PATCH 09/16] Move tile post processing --- inputs.cc | 89 ++++++++++++++++++++++++++----------------------------- tiles.hh | 20 ++++++++++--- 2 files changed, 58 insertions(+), 51 deletions(-) diff --git a/inputs.cc b/inputs.cc index 1d8652b..1d94c19 100644 --- a/inputs.cc +++ b/inputs.cc @@ -254,7 +254,7 @@ int loadLIDAR(char *filenames) } /* Allocate the tile array */ - tiles = (tile_t*)calloc(fc,sizeof(tile_t)); + tiles = (tile_t*) calloc(fc, sizeof(tile_t)); /* Load each tile in turn */ for (indx = 0; indx < fc; indx++) { @@ -279,6 +279,47 @@ int loadLIDAR(char *filenames) smCellsize = tiles[indx].cellsize; } + // Update a bunch of globals + if (tiles[indx].max_el > max_elevation) + max_elevation = tiles[indx].max_el; + if (tiles[indx].min_el < min_elevation) + min_elevation = tiles[indx].min_el; + + if (max_north == -90 || tiles[indx].max_north > max_north) + max_north = tiles[indx].max_north; + + if (min_north == 90 || tiles[indx].min_north < min_north) + min_north = tiles[indx].min_north; + + if (tiles[indx].max_west > max_west) + max_west = tiles[indx].max_west; + if (tiles[indx].min_west < min_west) + min_west = tiles[indx].min_west; + + if (max_west == -1) { + max_west = tiles[indx].max_west; + } else { + if (abs(tiles[indx].max_west - max_west) < 180) { + if (tiles[indx].max_west > max_west) + max_west = tiles[indx].max_west; + } else { + if (tiles[indx].max_west < max_west) + max_west = tiles[indx].max_west; + } + } + + if (min_west == 360) { + min_west = tiles[indx].min_west; + } else { + if (fabs(tiles[indx].min_west - min_west) < 180.0) { + if (tiles[indx].min_west < min_west) + min_west = tiles[indx].min_west; + } else { + if (tiles[indx].min_west > min_west) + min_west = tiles[indx].min_west; + } + } + } /* Iterate through all tiles to find the largest x/y dimension @@ -314,52 +355,6 @@ int loadLIDAR(char *filenames) dem[indx].max_el = tiles[indx].max_el; dem[indx].min_el = tiles[indx].min_el; - if (tiles[indx].max_el > max_elevation) - max_elevation = tiles[indx].max_el; - if (tiles[indx].min_el < min_elevation) - min_elevation = tiles[indx].min_el; - - 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) - min_west = dem[indx].min_west; - - if (max_west == -1) { - max_west = dem[indx].max_west; - } else { - if (abs(dem[indx].max_west - max_west) < 180) { - if (dem[indx].max_west > max_west) - max_west = dem[indx].max_west; - } else { - if (dem[indx].max_west < max_west) - max_west = dem[indx].max_west; - } - } - - if (min_west == 360) { - min_west = dem[indx].min_west; - } else { - if (fabs(dem[indx].min_west - min_west) < 180.0) { - if (dem[indx].min_west < min_west) - min_west = dem[indx].min_west; - } else { - if (dem[indx].min_west > min_west) - min_west = dem[indx].min_west; - } - } - /* * Copy the lidar tile data into the dem array. The dem array is rotated * 90 degrees (christ knows why...) diff --git a/tiles.hh b/tiles.hh index 2341faa..e76756f 100644 --- a/tiles.hh +++ b/tiles.hh @@ -11,10 +11,22 @@ typedef struct _tile_t{ int rows; int height; }; - double xll; - double yll; - double xur; - double yur; + union{ + double xll; + double max_west; + }; + union{ + double yll; + double min_north; + }; + union{ + double xur; + double min_west; + }; + union{ + double yur; + double max_north; + }; double cellsize; long long datastart; int nodata; From 545b7997e36ee5789515a02cc11b5cc8d3989afa Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 8 Jun 2017 20:22:04 +0100 Subject: [PATCH 10/16] Add feature to merge lidar tiles into single dem entry --- common.h | 2 + inputs.cc | 133 ++++++++++++++++++++++++++++++++++++++++-------------- tiles.cc | 96 ++++++++++++++++++++++++++++++++++++++- tiles.hh | 2 + 4 files changed, 198 insertions(+), 35 deletions(-) diff --git a/common.h b/common.h index 7abd1ad..73568a2 100644 --- a/common.h +++ b/common.h @@ -22,6 +22,8 @@ #define KM_PER_MILE 1.609344 #define FOUR_THIRDS 1.3333333333333 +#define MAX(x,y)((x)>(y)?(x):(y)) + struct dem { float min_north; float max_north; diff --git a/inputs.cc b/inputs.cc index 1d94c19..2956fdc 100644 --- a/inputs.cc +++ b/inputs.cc @@ -322,58 +322,123 @@ int loadLIDAR(char *filenames) } - /* Iterate through all tiles to find the largest x/y dimension - to use as the global IPPD value. We do this here so that we are - sure to allocate enough memory for the remainder of the calculations */ - int dimension_max = -1; + /* Iterate through all of the tiles to find the smallest resolution. We will + * need to rescale every tile from here on out to this value */ + int smallest_res = 0; + int pix_per_deg = 0; for (size_t i = 0; i < fc; i++) { - if( tiles[i].width > tiles[i].height && tiles[i].width > dimension_max ) - dimension_max = tiles[i].width; - else if( tiles[i].height > dimension_max ) - dimension_max = tiles[i].height; + if ( smallest_res == 0 || tiles[i].resolution < smallest_res ){ + smallest_res = tiles[i].resolution; + pix_per_deg = MAX(tiles[i].width,tiles[i].height) / MAX(tiles[i].max_north - tiles[i].min_north, tiles[i].max_west - tiles[i].min_west); + } } - MAXPAGES = fc; - IPPD = dimension_max; + /* Now we need to rescale all tiles the the lowest resolution. ie if we have + * one 1m lidar and one 2m lidar, resize the 2m to fake 1m */ + for (size_t i = 0; i< fc; i++) { + float rescale = tiles[i].resolution / smallest_res; + if (tiles[i].resolution != smallest_res) + tile_rescale(&tiles[i],rescale); + } + + /* Now we work out the size of the giant lidar tile. */ + double total_width = max_west - min_west; + double total_height = max_north - min_north; + /* This is how we should _theoretically_ work this out, but due to + * the nature of floating point arithmetic and rounding errors, we 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; + + size_t new_height = 0; + size_t new_width = 0; + for ( size_t i = 0; i < fc; i++ ) { + double north_offset = max_north - tiles[i].max_north; + double west_offset = max_west - tiles[i].max_west; + size_t north_pixel_offset = north_offset * pix_per_deg; + size_t west_pixel_offset = west_offset * pix_per_deg; + + if ( west_pixel_offset + tiles[i].width > new_width ) + new_width = west_pixel_offset + tiles[i].width; + if ( north_pixel_offset + tiles[i].height > new_height ) + new_height = north_pixel_offset + tiles[i].height; + } + size_t new_tile_alloc = new_width * new_height; + + int * new_tile = (int*) calloc( new_tile_alloc, sizeof(int) ); + if (debug) + fprintf(stderr,"Lidar tile dimensions w:%lf(%zu) h:%lf(%zu)\n", total_width, new_width, total_height, new_height); + + /* ...If we wanted a value other than sea level here, we would need to initialize the array... */ + + /* Fill out the array one tile at a time */ + for (size_t i = 0; i< fc; i++) { + double north_offset = max_north - tiles[i].max_north; + double west_offset = max_west - tiles[i].max_west; + size_t north_pixel_offset = north_offset * pix_per_deg; + size_t west_pixel_offset = west_offset * pix_per_deg; + + if (debug) { + fprintf(stderr,"mn: %lf mw:%lf globals: %lf %lf\n", tiles[i].max_north, tiles[i].max_west, max_north, max_west); + fprintf(stderr,"Offset n:%zu w:%zu\n", north_pixel_offset, west_pixel_offset); + fprintf(stderr,"Height: %d\n", tiles[i].height); + } + + /* Copy it row-by-row from the tile */ + for (size_t h = 0; h < tiles[i].height; h++) { + register int *dest_addr = &new_tile[ (north_pixel_offset+h)*new_width + west_pixel_offset]; + register int *src_addr = &tiles[i].data[h*tiles[i].width]; + // Check if we might overflow + if ( dest_addr + tiles[i].width > new_tile + new_tile_alloc ){ + if (debug) + fprintf(stderr, "Overflow %zu\n",i); + continue; + } + //fprintf(stderr,"dest:%p src:%p\n", dest_addr, src_addr); + memcpy( dest_addr, src_addr, tiles[i].width * sizeof(int) ); + } + } + + MAXPAGES = 1; + IPPD = MAX(new_width,new_height); if(debug){ fprintf(stderr,"Setting IPPD to %d\n",IPPD); fflush(stderr); } // add fudge as reprojected tiles sometimes vary by a pixel or ten - IPPD += 50; + // IPPD += 50; ARRAYSIZE = (MAXPAGES * IPPD) + 50; do_allocs(); // reset the IPPD after allocations - IPPD -= 50; + // IPPD -= 50; /* Load the data into the global dem array */ - for (size_t indx = 0; indx < fc; indx++) { - dem[indx].max_north = tiles[indx].yur; - dem[indx].min_west = tiles[indx].xur; - dem[indx].min_north = tiles[indx].yll; - dem[indx].max_west = tiles[indx].xll; - dem[indx].max_el = tiles[indx].max_el; - dem[indx].min_el = tiles[indx].min_el; + dem[0].max_north = max_north; + dem[0].min_west = min_west; + dem[0].min_north = min_north; + dem[0].max_west = max_west; + 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...) - */ - int y = tiles[indx].height-1; - for (size_t h = 0; h < tiles[indx].height; h++, y--) { - int x = tiles[indx].width-1; - for (size_t w = 0; w < tiles[indx].width; w++, x--) { - dem[indx].data[y][x] = tiles[indx].data[h*tiles[indx].width + w]; - dem[indx].signal[y][x] = 0; - dem[indx].mask[y][x] = 0; - } + /* + * Copy the lidar tile data into the dem array. The dem array is rotated + * 90 degrees (christ knows why...) + */ + int y = new_height-1; + for (size_t h = 0; h < new_height; h++, y--) { + int x = new_width-1; + for (size_t w = 0; w < new_width; w++, x--) { + dem[0].data[y][x] = new_tile[h*new_width + w]; + dem[0].signal[y][x] = 0; + dem[0].mask[y][x] = 0; } - } ippd=IPPD; - height = (unsigned)((max_north-min_north) / smCellsize); - width = (unsigned)((max_west-min_west) / smCellsize); + // height = (unsigned)((max_north-min_north) / smCellsize); + // width = (unsigned)((max_west-min_west) / smCellsize); + height = (unsigned)((max_north-min_north) * pix_per_deg); + width = (unsigned)((max_west-min_west) * pix_per_deg); if (debug) fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); diff --git a/tiles.cc b/tiles.cc index e5fb66d..8d18dea 100644 --- a/tiles.cc +++ b/tiles.cc @@ -2,11 +2,27 @@ #include #include #include +#include #include "tiles.hh" #include "common.h" #define MAX_LINE 25000 +/* Computes the distance between two long/lat points */ +double haversine_formulaz(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; +} + int tile_load_lidar(tile_t *tile, char *filename){ FILE *fd; char line[MAX_LINE]; @@ -100,6 +116,9 @@ int tile_load_lidar(tile_t *tile, char *filename){ }//if } + double current_res_km = haversine_formulaz(tile->max_north, tile->max_west, tile->max_north, tile->min_west); + tile->resolution = (int) ceil((current_res_km/MAX(tile->width,tile->height))*1000); + if (debug) fprintf(stderr,"Pixels loaded: %zu/%d\n",loaded,tile->width*tile->height); @@ -107,4 +126,79 @@ int tile_load_lidar(tile_t *tile, char *filename){ fclose(fd); return 0; -} \ No newline at end of file +} + +/* + * A positive scale will _increase_ the size of the data + */ +int tile_rescale(tile_t *tile, float scale){ + int *new_data; + size_t skip_count = 1; + size_t copy_count = 1; + + size_t new_height = tile->height * scale; + size_t new_width = tile->width * scale; + + /* Allocate the array for the lidar data */ + if ( (new_data = (int*) calloc(new_height * new_width, sizeof(int))) == NULL ) { + return ENOMEM; + } + + tile->max_el = -32768; + tile->min_el = 32768; + + /* Making the tile data smaller */ + if (scale < 0) { + skip_count = 1 / scale; + } else { + copy_count = (size_t) scale; + } + + fprintf(stderr,"Skip: %zu Copy: %zu\n", skip_count, copy_count); + + + + /* Nearest neighbour normalization. For each subsample of the original, simply + * assign the value in the top left to the new pixel */ + if (scale < 0){ + for (size_t x = 0, i = 0; i < new_width; x += skip_count, i++) { + for (size_t y = 0, j = 0; j < new_height; y += skip_count, j++){ + new_data[i*new_width+j] = tile->data[x*tile->width+y]; + /* Update local min / max values */ + if (tile->data[x * tile->width + y] > tile->max_el) + tile->max_el = tile->data[x * tile->width + y]; + if (tile->data[x * tile->width + y] < tile->min_el) + tile->min_el = tile->data[x * tile->width + y]; + } + } + }else{ + for (size_t x = 0; x < tile->width; x++) { + for (size_t y = 0; y < tile->height; y++){ + /* These are for scaling up the data */ + for (size_t copy_x = 0; copy_x < copy_count; copy_x++) { + for (size_t copy_y = 0; copy_y < copy_count; copy_y++) { + size_t new_x = (x * skip_count) + copy_x; + size_t new_y = (y * skip_count) + copy_y; + new_data[ new_x * new_width + new_y ] = tile->data[x * tile->width + y]; + // new_data[(x + copy_x) * new_width + (y + copy_y)] = tile->data[x * tile->width + y]; + } + } + /* Update local min / max values */ + if (tile->data[x * tile->width + y] > tile->max_el) + tile->max_el = tile->data[x * tile->width + y]; + if (tile->data[x * tile->width + y] < tile->min_el) + tile->min_el = tile->data[x * tile->width + y]; + } + } + } + + /* Update the date in the tile */ + free(tile->data); + tile->data = new_data; + + /* Update the height and width values */ + tile->height = new_height; + tile->width = new_width; + + return 0; +} diff --git a/tiles.hh b/tiles.hh index e76756f..357c081 100644 --- a/tiles.hh +++ b/tiles.hh @@ -33,8 +33,10 @@ typedef struct _tile_t{ int max_el; int min_el; int *data; + int resolution; } tile_t, *ptile_t; int tile_load_lidar(tile_t*, char *); +int tile_rescale(tile_t *, float); #endif \ No newline at end of file From c8f044380fa9f615a9c7e66e1c10301752a2388a Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 8 Jun 2017 20:32:58 +0100 Subject: [PATCH 11/16] Fix cropping switches --- main.cc | 59 +++++++++++++++++++++++++++++---------------------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/main.cc b/main.cc index 5fc0baa..24e7a42 100644 --- a/main.cc +++ b/main.cc @@ -57,7 +57,7 @@ int ippd, mpi, unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0; -bool to_stdout = false; +bool to_stdout = false, cropping = false; __thread double *elev; __thread struct path path; @@ -1885,26 +1885,29 @@ int main(int argc, char *argv[]) } } - // if(max_range<=100){ - // // CROPPING. croplat assigned in propPathLoss() - // max_north=cropLat; // MAX(path.lat[y]) - // // Edge case #1 - EAST/WEST - // if(cropLon>357 && tx_site[0].lon < 3) - // cropLon=tx_site[0].lon+3; - // // Edge case #2 - EAST/EAST - // if(cropLon>359.5 && tx_site[0].lon > 359.5) - // cropLon=362; - // max_west=cropLon; // MAX(path.lon[y]) - // cropLat-=tx_site[0].lat; // angle from tx to edge - // cropLon-=tx_site[0].lon; - // width=(int)((cropLon*ppd)*2); - // height=(int)((cropLat*ppd)*2); +#ifdef CROPPING + if(max_range<=100){ + // CROPPING. croplat assigned in propPathLoss() + max_north=cropLat; // MAX(path.lat[y]) + // Edge case #1 - EAST/WEST + if(cropLon>357 && tx_site[0].lon < 3) + cropLon=tx_site[0].lon+3; + // Edge case #2 - EAST/EAST + if(cropLon>359.5 && tx_site[0].lon > 359.5) + cropLon=362; + max_west=cropLon; // MAX(path.lon[y]) + cropLat-=tx_site[0].lat; // angle from tx to edge + cropLon-=tx_site[0].lon; + width=(int)((cropLon*ppd)*2); + height=(int)((cropLat*ppd)*2); + cropping = true; - // if(width>3600*10){ - // fprintf(stderr,"FATAL BOUNDS! max_west: %.4f cropLat: %.4f cropLon: %.4f longitude: %.5f\n",max_west,cropLat,cropLon,tx_site[0].lon); - // return 0; - // } - // } + if(width>3600*10){ + fprintf(stderr,"FATAL BOUNDS! max_west: %.4f cropLat: %.4f cropLon: %.4f longitude: %.5f\n",max_west,cropLat,cropLon,tx_site[0].lon); + return 0; + } + } +#endif // Write bitmap if (LR.erp == 0.0) @@ -1926,21 +1929,19 @@ int main(int argc, char *argv[]) tx_site[0].lon *= -1; } if (tx_site[0].lon < -180.0){ - tx_site[0].lon += 360; + tx_site[0].lon += 360; } - if (propmodel == 2 || max_range > 100) { - // No croppping because this is LOS - fprintf(stderr, "|%.6f", max_north); - fprintf(stderr, "|%.6f", east); - fprintf(stderr, "|%.6f", min_north); - fprintf(stderr, "|%.6f|",west); - }else{ - // Cropped EPSG4326 coordinates + if (cropping) { fprintf(stderr, "|%.6f", tx_site[0].lat+cropLat); fprintf(stderr, "|%.6f", tx_site[0].lon+cropLon); fprintf(stderr, "|%.6f", tx_site[0].lat-cropLat); fprintf(stderr, "|%.6f|",tx_site[0].lon-cropLon); + }else{ + fprintf(stderr, "|%.6f", max_north); + fprintf(stderr, "|%.6f", east); + fprintf(stderr, "|%.6f", min_north); + fprintf(stderr, "|%.6f|",west); } fprintf(stderr, "\n"); From 5ea0e1f1c57301822d62b8ba36ef761d942d71dc Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Thu, 8 Jun 2017 23:49:04 +0100 Subject: [PATCH 12/16] Fixup some arithmetic. Still no grenwich --- inputs.cc | 27 ++++++++++++++++++--------- tiles.cc | 12 ++++++------ 2 files changed, 24 insertions(+), 15 deletions(-) diff --git a/inputs.cc b/inputs.cc index 2956fdc..94debcb 100644 --- a/inputs.cc +++ b/inputs.cc @@ -329,7 +329,7 @@ int loadLIDAR(char *filenames) for (size_t i = 0; i < fc; i++) { if ( smallest_res == 0 || tiles[i].resolution < smallest_res ){ smallest_res = tiles[i].resolution; - pix_per_deg = MAX(tiles[i].width,tiles[i].height) / MAX(tiles[i].max_north - tiles[i].min_north, tiles[i].max_west - tiles[i].min_west); + pix_per_deg = MAX(tiles[i].width,tiles[i].height) / MAX(tiles[i].max_north - tiles[i].min_north, tiles[i].max_west - tiles[i].min_west >= 0 ? tiles[i].max_west - tiles[i].min_west : tiles[i].max_west + (360 - tiles[i].min_west)); } } @@ -342,8 +342,12 @@ int loadLIDAR(char *filenames) } /* Now we work out the size of the giant lidar tile. */ - double total_width = max_west - min_west; + 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, "totalwidth: %.7f - %.7f = %.7f\n", max_west, min_west, total_width); + fprintf(stderr,"mw:%lf Mnw:%lf\n", max_west, min_west); + } /* This is how we should _theoretically_ work this out, but due to * the nature of floating point arithmetic and rounding errors, we need to * crunch the numbers the hard way */ @@ -354,7 +358,7 @@ int loadLIDAR(char *filenames) size_t new_width = 0; for ( size_t i = 0; i < fc; i++ ) { double north_offset = max_north - tiles[i].max_north; - double west_offset = max_west - tiles[i].max_west; + 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 * pix_per_deg; size_t west_pixel_offset = west_offset * pix_per_deg; @@ -366,6 +370,10 @@ int loadLIDAR(char *filenames) size_t new_tile_alloc = new_width * new_height; int * new_tile = (int*) calloc( new_tile_alloc, sizeof(int) ); + if ( new_tile == NULL ){ + free(tiles); + return ENOMEM; + } if (debug) fprintf(stderr,"Lidar tile dimensions w:%lf(%zu) h:%lf(%zu)\n", total_width, new_width, total_height, new_height); @@ -374,13 +382,14 @@ int loadLIDAR(char *filenames) /* Fill out the array one tile at a time */ for (size_t i = 0; i< fc; i++) { double north_offset = max_north - tiles[i].max_north; - double west_offset = max_west - tiles[i].max_west; + 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 * pix_per_deg; size_t west_pixel_offset = west_offset * pix_per_deg; + if (debug) { fprintf(stderr,"mn: %lf mw:%lf globals: %lf %lf\n", tiles[i].max_north, tiles[i].max_west, max_north, max_west); - fprintf(stderr,"Offset n:%zu w:%zu\n", north_pixel_offset, west_pixel_offset); + fprintf(stderr,"Offset n:%zu(%lf) w:%zu(%lf)\n", north_pixel_offset, north_offset, west_pixel_offset, west_offset); fprintf(stderr,"Height: %d\n", tiles[i].height); } @@ -389,12 +398,12 @@ int loadLIDAR(char *filenames) register int *dest_addr = &new_tile[ (north_pixel_offset+h)*new_width + west_pixel_offset]; register int *src_addr = &tiles[i].data[h*tiles[i].width]; // Check if we might overflow - if ( dest_addr + tiles[i].width > new_tile + new_tile_alloc ){ + if ( dest_addr + tiles[i].width > new_tile + new_tile_alloc || dest_addr < new_tile ){ if (debug) fprintf(stderr, "Overflow %zu\n",i); continue; } - //fprintf(stderr,"dest:%p src:%p\n", dest_addr, src_addr); + // fprintf(stderr,"dest:%p src:%p\n", dest_addr, src_addr); memcpy( dest_addr, src_addr, tiles[i].width * sizeof(int) ); } } @@ -437,8 +446,8 @@ int loadLIDAR(char *filenames) ippd=IPPD; // height = (unsigned)((max_north-min_north) / smCellsize); // width = (unsigned)((max_west-min_west) / smCellsize); - height = (unsigned)((max_north-min_north) * pix_per_deg); - width = (unsigned)((max_west-min_west) * pix_per_deg); + height = (unsigned)((total_height) * pix_per_deg); + width = (unsigned)((total_width) * pix_per_deg); if (debug) fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); diff --git a/tiles.cc b/tiles.cc index 8d18dea..43cffeb 100644 --- a/tiles.cc +++ b/tiles.cc @@ -67,11 +67,11 @@ int tile_load_lidar(tile_t *tile, char *filename){ // fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur); // Greenwich straddling hack - if (tile->xll <= 0 && tile->xur > 0) { - tile->xll = (tile->xur - tile->xll); // full width - tile->xur = 0.0; // budge it along so it's west of greenwich - delta = eastoffset; // add to Tx longitude later - } else { + // if (tile->xll <= 0 && tile->xur > 0) { + // tile->xll = (tile->xur - tile->xll); // full width + // tile->xur = 0.0; // budge it along so it's west of greenwich + // delta = eastoffset; // add to Tx longitude later + // } else { // Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;) if (tile->xll >= 0) tile->xll = 360-tile->xll; @@ -81,7 +81,7 @@ int tile_load_lidar(tile_t *tile, char *filename){ tile->xll = tile->xll * -1; if(tile->xur < 0) tile->xur = tile->xur * -1; - } + // } if (debug) fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", tile->yll, tile->yur, tile->xur, tile->xll, delta); From cb2939b244a2f4ca5859dfdc6b0ddf7db097c612 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Wed, 14 Jun 2017 18:44:29 +0100 Subject: [PATCH 13/16] Fix LIDAR tile resampling --- inputs.cc | 147 ++++++++---------------------------------------------- inputs.hh | 2 +- main.cc | 20 +++----- tiles.cc | 101 ++++++++++++++++++++++++------------- tiles.hh | 4 ++ 5 files changed, 100 insertions(+), 174 deletions(-) diff --git a/inputs.cc b/inputs.cc index 94debcb..395f0e5 100644 --- a/inputs.cc +++ b/inputs.cc @@ -9,112 +9,6 @@ #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 - * 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 */ - if (debug) - 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 */ - height /= scaling_factor; - width /= scaling_factor; - IPPD = new_ippd; - ippd = IPPD; - ARRAYSIZE = (MAXPAGES * IPPD) + 50; - 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. @@ -231,7 +125,7 @@ int loadClutter(char *filename, double radius, struct site tx) return 0; } -int loadLIDAR(char *filenames) +int loadLIDAR(char *filenames, int resample) { char *filename; char *files[100]; // 10x10 tiles @@ -254,7 +148,8 @@ int loadLIDAR(char *filenames) } /* Allocate the tile array */ - tiles = (tile_t*) calloc(fc, sizeof(tile_t)); + if( (tiles = (tile_t*) calloc(fc, sizeof(tile_t))) == NULL ) + return ENOMEM; /* Load each tile in turn */ for (indx = 0; indx < fc; indx++) { @@ -325,20 +220,26 @@ int loadLIDAR(char *filenames) /* Iterate through all of the tiles to find the smallest resolution. We will * need to rescale every tile from here on out to this value */ int smallest_res = 0; - int pix_per_deg = 0; for (size_t i = 0; i < fc; i++) { if ( smallest_res == 0 || tiles[i].resolution < smallest_res ){ smallest_res = tiles[i].resolution; - pix_per_deg = MAX(tiles[i].width,tiles[i].height) / MAX(tiles[i].max_north - tiles[i].min_north, tiles[i].max_west - tiles[i].min_west >= 0 ? tiles[i].max_west - tiles[i].min_west : tiles[i].max_west + (360 - tiles[i].min_west)); } } - /* Now we need to rescale all tiles the the lowest resolution. ie if we have + /* 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 */ + int desired_resolution = smallest_res < resample ? resample : smallest_res; + if (desired_resolution > resample && debug ) + fprintf(stderr, "Warning: Unable to rescale to requested resolution\n"); for (size_t i = 0; i< fc; i++) { - float rescale = tiles[i].resolution / smallest_res; - if (tiles[i].resolution != smallest_res) - tile_rescale(&tiles[i],rescale); + float rescale = 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. */ @@ -359,8 +260,8 @@ int loadLIDAR(char *filenames) for ( size_t i = 0; i < fc; i++ ) { double north_offset = max_north - tiles[i].max_north; 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 * pix_per_deg; - size_t west_pixel_offset = west_offset * pix_per_deg; + size_t north_pixel_offset = north_offset * tiles[i].ppdy; + size_t west_pixel_offset = west_offset * tiles[i].ppdx; if ( west_pixel_offset + tiles[i].width > new_width ) new_width = west_pixel_offset + tiles[i].width; @@ -383,8 +284,8 @@ int loadLIDAR(char *filenames) for (size_t i = 0; i< fc; i++) { double north_offset = max_north - tiles[i].max_north; 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 * pix_per_deg; - size_t west_pixel_offset = west_offset * pix_per_deg; + size_t north_pixel_offset = north_offset * tiles[i].ppdy; + size_t west_pixel_offset = west_offset * tiles[i].ppdx; if (debug) { @@ -414,12 +315,8 @@ int loadLIDAR(char *filenames) fprintf(stderr,"Setting IPPD to %d\n",IPPD); fflush(stderr); } - // add fudge as reprojected tiles sometimes vary by a pixel or ten - // IPPD += 50; ARRAYSIZE = (MAXPAGES * IPPD) + 50; do_allocs(); - // reset the IPPD after allocations - // IPPD -= 50; /* Load the data into the global dem array */ dem[0].max_north = max_north; @@ -444,10 +341,8 @@ int loadLIDAR(char *filenames) } ippd=IPPD; - // height = (unsigned)((max_north-min_north) / smCellsize); - // width = (unsigned)((max_west-min_west) / smCellsize); - height = (unsigned)((total_height) * pix_per_deg); - width = (unsigned)((total_width) * pix_per_deg); + height = new_height; + width = new_width; if (debug) fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height); diff --git a/inputs.hh b/inputs.hh index 2d79a52..6d7d353 100644 --- a/inputs.hh +++ b/inputs.hh @@ -15,7 +15,7 @@ int LoadLossColors(struct site xmtr); int LoadDBMColors(struct site xmtr); int LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat); int LoadUDT(char *filename); -int loadLIDAR(char *filename); +int loadLIDAR(char *filename, int resample); int loadClutter(char *filename, double radius, struct site tx); static const char AZ_FILE_SUFFIX[] = ".az"; diff --git a/main.cc b/main.cc index 24e7a42..b0d2549 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, resample; + 0, MAXRAD, hottest = 10, height, width, resample = 0; unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0; @@ -1327,8 +1327,11 @@ 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"); + if(!lidar){ + fprintf(stderr, "Error, this should only be used with LIDAR tiles.\n"); + return -1; + } + sscanf(argv[z], "%d", &resample); } @@ -1697,21 +1700,14 @@ int main(int argc, char *argv[]) /* Load the required tiles */ if(lidar){ - if( (result = loadLIDAR(lidar_tiles)) != 0 ){ + if( (result = loadLIDAR(lidar_tiles, resample)) != 0 ){ fprintf(stderr, "Couldn't find one or more of the " "lidar files. Please ensure their paths are " "correct and try again.\n"); + fprintf(stderr, "Error %d: %s\n", result, strerror(result)); exit(result); } - /* If we have been asked to resample the input data; do it now. */ - if (resample != -1 ){ - if ((result = resize_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); } diff --git a/tiles.cc b/tiles.cc index 43cffeb..c4b47b4 100644 --- a/tiles.cc +++ b/tiles.cc @@ -9,7 +9,7 @@ #define MAX_LINE 25000 /* Computes the distance between two long/lat points */ -double haversine_formulaz(double th1, double ph1, double th2, double ph2) +double haversine_formula(double th1, double ph1, double th2, double ph2) { #define TO_RAD (3.1415926536 / 180) int R = 6371; @@ -116,11 +116,17 @@ int tile_load_lidar(tile_t *tile, char *filename){ }//if } - double current_res_km = haversine_formulaz(tile->max_north, tile->max_west, tile->max_north, tile->min_west); + double current_res_km = haversine_formula(tile->max_north, tile->max_west, tile->max_north, tile->min_west); tile->resolution = (int) ceil((current_res_km/MAX(tile->width,tile->height))*1000); + 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; + + tile->ppdx = tile->width / tile->width_deg; + tile->ppdy = tile->height / tile->height_deg; + if (debug) - fprintf(stderr,"Pixels loaded: %zu/%d\n",loaded,tile->width*tile->height); + fprintf(stderr,"Pixels loaded: %zu/%d (PPD %dx%d)\n", loaded, tile->width*tile->height, tile->ppdx, tile->ppdy); /* All done, close the LIDAR file */ fclose(fd); @@ -129,13 +135,23 @@ int tile_load_lidar(tile_t *tile, char *filename){ } /* - * A positive scale will _increase_ the size of the data + * tile_rescale + * 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 tile_rescale(tile_t *tile, float scale){ int *new_data; size_t skip_count = 1; size_t copy_count = 1; + if (scale == 1) { + return 0; + } + size_t new_height = tile->height * scale; size_t new_width = tile->width * scale; @@ -148,47 +164,42 @@ int tile_rescale(tile_t *tile, float scale){ tile->min_el = 32768; /* Making the tile data smaller */ - if (scale < 0) { + if (scale < 1) { skip_count = 1 / scale; } else { copy_count = (size_t) scale; } - fprintf(stderr,"Skip: %zu Copy: %zu\n", skip_count, copy_count); - - + if (debug) + fprintf(stderr,"Resampling tile:\n\tOld %zux%zu. New %zux%zu\n\tScale %f Skip %zu Copy %zu\n", tile->width, tile->height, new_width, new_height, scale, skip_count, copy_count); /* Nearest neighbour normalization. For each subsample of the original, simply - * assign the value in the top left to the new pixel */ - if (scale < 0){ - for (size_t x = 0, i = 0; i < new_width; x += skip_count, i++) { - for (size_t y = 0, j = 0; j < new_height; y += skip_count, j++){ - new_data[i*new_width+j] = tile->data[x*tile->width+y]; - /* Update local min / max values */ - if (tile->data[x * tile->width + y] > tile->max_el) - tile->max_el = tile->data[x * tile->width + y]; - if (tile->data[x * tile->width + y] < tile->min_el) - tile->min_el = tile->data[x * tile->width + y]; - } - } - }else{ - for (size_t x = 0; x < tile->width; x++) { - for (size_t y = 0; y < tile->height; y++){ - /* These are for scaling up the data */ + * assign the value in the top left to the new pixel + * SOURCE: X / Y + * DEST: I / J */ + + for (size_t y = 0, j = 0; + y < tile->height && j < new_height; + y += skip_count, j += copy_count){ + + for (size_t x = 0, i = 0; + x < tile->width && i < new_width; + x += skip_count, i += copy_count) { + + /* These are for scaling up the data */ + for (size_t copy_y = 0; copy_y < copy_count; copy_y++) { for (size_t copy_x = 0; copy_x < copy_count; copy_x++) { - for (size_t copy_y = 0; copy_y < copy_count; copy_y++) { - size_t new_x = (x * skip_count) + copy_x; - size_t new_y = (y * skip_count) + copy_y; - new_data[ new_x * new_width + new_y ] = tile->data[x * tile->width + y]; - // new_data[(x + copy_x) * new_width + (y + copy_y)] = tile->data[x * tile->width + y]; - } + size_t new_j = j + copy_y; + size_t new_i = i + copy_x; + /* Do the copy */ + new_data[ new_j * new_width + new_i ] = tile->data[y * tile->width + x]; } - /* Update local min / max values */ - if (tile->data[x * tile->width + y] > tile->max_el) - tile->max_el = tile->data[x * tile->width + y]; - if (tile->data[x * tile->width + y] < tile->min_el) - tile->min_el = tile->data[x * tile->width + y]; } + /* Update local min / max values */ + if (tile->data[y * tile->width + x] > tile->max_el) + tile->max_el = tile->data[y * tile->width + x]; + if (tile->data[y * tile->width + x] < tile->min_el) + tile->min_el = tile->data[y * tile->width + x]; } } @@ -199,6 +210,26 @@ int tile_rescale(tile_t *tile, float scale){ /* Update the height and width values */ tile->height = new_height; tile->width = new_width; + tile->resolution *= scale; + tile->ppdy *= scale; + tile->ppdx *= scale; + tile->width_deg *= scale; + tile->height_deg *= scale; return 0; } + +/* + * tile_resize + * 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 tile_resize(tile_t* tile, int resolution){ + double current_res_km = haversine_formula(tile->max_north, tile->max_west, tile->max_north, tile->min_west); + 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); + return tile_rescale(tile, scaling_factor); +} diff --git a/tiles.hh b/tiles.hh index 357c081..899ef92 100644 --- a/tiles.hh +++ b/tiles.hh @@ -34,6 +34,10 @@ typedef struct _tile_t{ int min_el; int *data; int resolution; + double width_deg; + double height_deg; + int ppdx; + int ppdy; } tile_t, *ptile_t; int tile_load_lidar(tile_t*, char *); From 02bd67f11bbd2c2d2a77b6f5d888d778bd5e4fb8 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Wed, 14 Jun 2017 20:33:32 +0100 Subject: [PATCH 14/16] Change data to short to restrain memory usage --- inputs.cc | 10 +++++----- main.cc | 2 +- tiles.cc | 10 +++++----- tiles.hh | 8 ++++---- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/inputs.cc b/inputs.cc index 395f0e5..995e3c4 100644 --- a/inputs.cc +++ b/inputs.cc @@ -228,7 +228,7 @@ 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 */ - int desired_resolution = smallest_res < resample ? resample : smallest_res; + int desired_resolution = resample != 0 && smallest_res < resample ? resample : smallest_res; if (desired_resolution > resample && debug ) fprintf(stderr, "Warning: Unable to rescale to requested resolution\n"); for (size_t i = 0; i< fc; i++) { @@ -270,7 +270,7 @@ int loadLIDAR(char *filenames, int resample) } size_t new_tile_alloc = new_width * new_height; - int * new_tile = (int*) calloc( new_tile_alloc, sizeof(int) ); + short * new_tile = (short*) calloc( new_tile_alloc, sizeof(int) ); if ( new_tile == NULL ){ free(tiles); return ENOMEM; @@ -296,8 +296,8 @@ int loadLIDAR(char *filenames, int resample) /* Copy it row-by-row from the tile */ for (size_t h = 0; h < tiles[i].height; h++) { - register int *dest_addr = &new_tile[ (north_pixel_offset+h)*new_width + west_pixel_offset]; - register int *src_addr = &tiles[i].data[h*tiles[i].width]; + register short *dest_addr = &new_tile[ (north_pixel_offset+h)*new_width + west_pixel_offset]; + register short *src_addr = &tiles[i].data[h*tiles[i].width]; // Check if we might overflow if ( dest_addr + tiles[i].width > new_tile + new_tile_alloc || dest_addr < new_tile ){ if (debug) @@ -305,7 +305,7 @@ int loadLIDAR(char *filenames, int resample) continue; } // fprintf(stderr,"dest:%p src:%p\n", dest_addr, src_addr); - memcpy( dest_addr, src_addr, tiles[i].width * sizeof(int) ); + memcpy( dest_addr, src_addr, tiles[i].width * sizeof(short) ); } } diff --git a/main.cc b/main.cc index b0d2549..cc19762 100644 --- a/main.cc +++ b/main.cc @@ -1146,7 +1146,7 @@ int main(int argc, char *argv[]) max_txsites = 30; fzone_clearance = 0.6; contour_threshold = 0; - resample = -1; + resample = 0; ano_filename[0] = 0; earthradius = EARTHRADIUS; diff --git a/tiles.cc b/tiles.cc index c4b47b4..05912af 100644 --- a/tiles.cc +++ b/tiles.cc @@ -26,11 +26,11 @@ double haversine_formula(double th1, double ph1, double th2, double ph2) int tile_load_lidar(tile_t *tile, char *filename){ FILE *fd; char line[MAX_LINE]; - int nextval; + short nextval; char *pch; /* Clear the tile data */ - memset(tile,0x00,sizeof(tile_t)); + memset(tile, 0x00, sizeof(tile_t)); /* Open the file handle and return on error */ if ( (fd = fopen(filename,"r")) == NULL ) @@ -88,7 +88,7 @@ int tile_load_lidar(tile_t *tile, char *filename){ /* Read the actual tile data */ /* Allocate the array for the lidar data */ - if ( (tile->data = (int*) calloc(tile->width * tile->height, sizeof(int))) == NULL ) { + if ( (tile->data = (short*) calloc(tile->width * tile->height, sizeof(short))) == NULL ) { fclose(fd); free(tile->filename); return ENOMEM; @@ -144,7 +144,7 @@ int tile_load_lidar(tile_t *tile, char *filename){ * (ie 2m LIDAR can be 4/6/8/... and 20m can be 40/60) */ int tile_rescale(tile_t *tile, float scale){ - int *new_data; + short *new_data; size_t skip_count = 1; size_t copy_count = 1; @@ -156,7 +156,7 @@ int tile_rescale(tile_t *tile, float scale){ size_t new_width = tile->width * scale; /* Allocate the array for the lidar data */ - if ( (new_data = (int*) calloc(new_height * new_width, sizeof(int))) == NULL ) { + if ( (new_data = (short*) calloc(new_height * new_width, sizeof(short))) == NULL ) { return ENOMEM; } diff --git a/tiles.hh b/tiles.hh index 899ef92..8bb0caf 100644 --- a/tiles.hh +++ b/tiles.hh @@ -29,10 +29,10 @@ typedef struct _tile_t{ }; double cellsize; long long datastart; - int nodata; - int max_el; - int min_el; - int *data; + short nodata; + short max_el; + short min_el; + short *data; int resolution; double width_deg; double height_deg; From ff3bc0c164f17b09ecc9e705e843bb81f309a2a9 Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Wed, 14 Jun 2017 20:40:52 +0100 Subject: [PATCH 15/16] Add cleanup code to LIDAR resampling --- inputs.cc | 11 ++++++++++- tiles.cc | 10 ++++++++++ tiles.hh | 1 + 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/inputs.cc b/inputs.cc index 995e3c4..9017348 100644 --- a/inputs.cc +++ b/inputs.cc @@ -304,7 +304,6 @@ int loadLIDAR(char *filenames, int resample) fprintf(stderr, "Overflow %zu\n",i); continue; } - // fprintf(stderr,"dest:%p src:%p\n", dest_addr, src_addr); memcpy( dest_addr, src_addr, tiles[i].width * sizeof(short) ); } } @@ -349,6 +348,16 @@ int loadLIDAR(char *filenames, int resample) if (debug) fprintf(stderr, "fc %d WIDTH %d HEIGHT %d ippd %d minN %.5f maxN %.5f minW %.5f maxW %.5f avgCellsize %.5f\n", fc, width, height, ippd,min_north,max_north,min_west,max_west,avgCellsize); + +cleanup: + + if ( tiles != NULL ) { + for (size_t i = 0; i < fc; i++) { + tile_destroy(&tiles[i]); + } + } + free(tiles); + return 0; } diff --git a/tiles.cc b/tiles.cc index 05912af..f81bcb7 100644 --- a/tiles.cc +++ b/tiles.cc @@ -233,3 +233,13 @@ int tile_resize(tile_t* tile, int resolution){ fprintf(stderr, "Resampling: Current %dm Desired %dm Scale %d\n", current_res, resolution, scaling_factor); return tile_rescale(tile, scaling_factor); } + +/* + * tile_destroy + * This function simply destroys any data associated with a tile + */ +void tile_destroy(tile_t* tile){ + if (tile->data != NULL) + free(tile->data); +} + diff --git a/tiles.hh b/tiles.hh index 8bb0caf..a2d25f4 100644 --- a/tiles.hh +++ b/tiles.hh @@ -42,5 +42,6 @@ typedef struct _tile_t{ int tile_load_lidar(tile_t*, char *); int tile_rescale(tile_t *, float); +void tile_destroy(tile_t *); #endif \ No newline at end of file From 8f13dae8b4c2703ab5d6213f53d2407d5e0d1cca Mon Sep 17 00:00:00 2001 From: Gareth Evans Date: Wed, 14 Jun 2017 21:32:51 +0100 Subject: [PATCH 16/16] Minor short type fix --- inputs.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inputs.cc b/inputs.cc index 9017348..00ebceb 100644 --- a/inputs.cc +++ b/inputs.cc @@ -167,7 +167,7 @@ int loadLIDAR(char *filenames, int resample) fflush(stderr); } - // Increase the average cell size + // Increase the "average" cell size avgCellsize += tiles[indx].cellsize; // Update the smallest cell size if (smCellsize == 0 || tiles[indx].cellsize < smCellsize) { @@ -270,7 +270,7 @@ int loadLIDAR(char *filenames, int resample) } size_t new_tile_alloc = new_width * new_height; - short * new_tile = (short*) calloc( new_tile_alloc, sizeof(int) ); + short * new_tile = (short*) calloc( new_tile_alloc, sizeof(short) ); if ( new_tile == NULL ){ free(tiles); return ENOMEM;