HFNOC/Signal-Server
2
0
Fork 0
forked from ExternalVendorCode/Signal-Server
Code Pull Requests Activity

Merge pull request #8 from kryc/feature-lidar

Browse Source 
Improved LIDAR support with resampling and multi-resolution, mixed tile size handling.
This commit is contained in:
CloudRF
2017-06-18 21:50:42 +01:00
committed by GitHub
parent 9038f32c3f 8f13dae8b4
commit 5442357bea
8 changed files with 532 additions and 199 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
Makefile
View File

@@ -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

2
common.h
View File

@@ -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;

374
inputs.cc
View File

@@ -7,7 +7,7 @@
#include <limits.h>
#include "common.h"
#include "main.hh"
#include "tiles.hh"
int loadClutter(char *filename, double radius, struct site tx)
{
@@ -125,105 +125,19 @@ 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)
int loadLIDAR(char *filenames, int resample)
{
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,113 +147,217 @@ int loadLIDAR(char *filenames)
fc++;
}
while (indx < fc) {
if( (fd = fopen(files[indx], "rb")) == NULL )
return errno;
/* Allocate the tile array */
if( (tiles = (tile_t*) calloc(fc, sizeof(tile_t))) == NULL )
return ENOMEM;
if (fgets(line, 255, fd) != NULL) {
pch = strtok (line," ");
pch = strtok (NULL, " ");
width = atoi(pch); // ncols
/* Load each tile in turn */
for (indx = 0; indx < fc; indx++) {
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
/* 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 (fgets(line, 255, fd) != NULL) {
sscanf(pch, "%lf", &yll); // yll
if (debug) {
fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, tiles[indx].width);
fflush(stderr);
}
if (fgets(line, 255, fd) != NULL)
sscanf(pch, "%lf", &cellsize);
// 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;
}
avgCellsize=avgCellsize+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(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 (max_north == -90 || tiles[indx].max_north > max_north)
max_north = tiles[indx].max_north;
if (xur > eastoffset)
eastoffset = xur;
if (xll < westoffset)
westoffset = xll;
if (min_north == 90 || tiles[indx].min_north < min_north)
min_north = tiles[indx].min_north;
if (debug)
fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur);
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;
// 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
if (max_west == -1) {
max_west = tiles[indx].max_west;
} 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 (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 (debug)
fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta);
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 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;
for (size_t i = 0; i < fc; i++) {
if ( smallest_res == 0 || tiles[i].resolution < smallest_res ){
smallest_res = tiles[i].resolution;
}
}
/* 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 = 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++) {
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. */
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 */
// 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 >= 0 ? max_west - tiles[i].max_west : max_west + (360 - tiles[i].max_west);
size_t north_pixel_offset = north_offset * tiles[i].ppdy;
size_t west_pixel_offset = west_offset * tiles[i].ppdx;
if ( 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;
short * new_tile = (short*) calloc( new_tile_alloc, sizeof(short) );
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);
/* ...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 >= 0 ? max_west - tiles[i].max_west : max_west + (360 - tiles[i].max_west);
size_t north_pixel_offset = north_offset * tiles[i].ppdy;
size_t west_pixel_offset = west_offset * tiles[i].ppdx;
fgets(line, 255, fd); // NODATA
pos = ftell(fd);
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(%lf) w:%zu(%lf)\n", north_pixel_offset, north_offset, west_pixel_offset, west_offset);
fprintf(stderr,"Height: %d\n", tiles[i].height);
}
// tile 0 [x| ]
if (debug)
fprintf(stderr, "readLIDAR(fd,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", height, width, indx, yur, xur, yll, xll);
/* Copy it row-by-row from the tile */
for (size_t h = 0; h < tiles[i].height; h++) {
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)
fprintf(stderr, "Overflow %zu\n",i);
continue;
}
memcpy( dest_addr, src_addr, tiles[i].width * sizeof(short) );
}
}
readLIDAR(fd, height, width, indx, yur, xur, yll, xll);
MAXPAGES = 1;
IPPD = MAX(new_width,new_height);
if(debug){
fprintf(stderr,"Setting IPPD to %d\n",IPPD);
fflush(stderr);
}
ARRAYSIZE = (MAXPAGES * IPPD) + 50;
do_allocs();
/* Load the data into the global dem array */
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 = 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;
}
}
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 = new_height;
width = new_width;
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);
cleanup:
if ( tiles != NULL ) {
for (size_t i = 0; i < fc; i++) {
tile_destroy(&tiles[i]);
}
}
free(tiles);
return 0;
}

6
inputs.hh
View File

@@ -3,6 +3,10 @@
#include "common.h"
/* 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);
int LoadPAT(char *az_filename, char *el_filename);
@@ -11,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";

45
main.cc
View File

@@ -53,11 +53,11 @@ 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 = 0;
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;
@@ -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");
@@ -1145,6 +1146,7 @@ int main(int argc, char *argv[])
max_txsites = 30;
fzone_clearance = 0.6;
contour_threshold = 0;
resample = 0;
ano_filename[0] = 0;
earthradius = EARTHRADIUS;
@@ -1322,6 +1324,17 @@ int main(int argc, char *argv[])
}
}
if (strcmp(argv[x], "-resample") == 0) {
z = x + 1;
if(!lidar){
fprintf(stderr, "Error, this should only be used with LIDAR tiles.\n");
return -1;
}
sscanf(argv[z], "%d", &resample);
}
if (strcmp(argv[x], "-lat") == 0) {
z = x + 1;
@@ -1687,10 +1700,11 @@ 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);
}
@@ -1710,9 +1724,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;
@@ -1867,6 +1881,7 @@ int main(int argc, char *argv[])
}
}
#ifdef CROPPING
if(max_range<=100){
// CROPPING. croplat assigned in propPathLoss()
max_north=cropLat; // MAX(path.lat[y])
@@ -1881,12 +1896,14 @@ int main(int argc, char *argv[])
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;
}
}
#endif
// Write bitmap
if (LR.erp == 0.0)
@@ -1908,21 +1925,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");

8
outputs.cc
View File

@@ -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);
}

245
tiles.cc Normal file
View File

@@ -0,0 +1,245 @@
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <math.h>
#include "tiles.hh"
#include "common.h"
#define MAX_LINE 25000
/* 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;
}
int tile_load_lidar(tile_t *tile, char *filename){
FILE *fd;
char line[MAX_LINE];
short nextval;
char *pch;
/* 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);
/* Read the actual tile data */
/* Allocate the array for the lidar data */
if ( (tile->data = (short*) calloc(tile->width * tile->height, sizeof(short))) == 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 (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)
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 @ h %zu file %s\n", h, filename);
}//if
}
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 (PPD %dx%d)\n", loaded, tile->width*tile->height, tile->ppdx, tile->ppdy);
/* All done, close the LIDAR file */
fclose(fd);
return 0;
}
/*
* 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){
short *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;
/* Allocate the array for the lidar data */
if ( (new_data = (short*) calloc(new_height * new_width, sizeof(short))) == NULL ) {
return ENOMEM;
}
tile->max_el = -32768;
tile->min_el = 32768;
/* Making the tile data smaller */
if (scale < 1) {
skip_count = 1 / scale;
} else {
copy_count = (size_t) scale;
}
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
* 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++) {
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[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];
}
}
/* 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;
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);
}
/*
* 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);
}

47
tiles.hh Normal file
View File

@@ -0,0 +1,47 @@
#ifndef _TILES_HH_
#define _TILES_HH_
typedef struct _tile_t{
char *filename;
union{
int cols;
int width;
};
union{
int rows;
int height;
};
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;
short nodata;
short max_el;
short min_el;
short *data;
int resolution;
double width_deg;
double height_deg;
int ppdx;
int ppdy;
} tile_t, *ptile_t;
int tile_load_lidar(tile_t*, char *);
int tile_rescale(tile_t *, float);
void tile_destroy(tile_t *);
#endif