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 *);