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

Add support for working out IPPD based on LIDAR tiles

Browse Source
This commit is contained in:
Gareth Evans
2017-06-01 18:54:49 +01:00
parent ee67a86a9f
commit 83cdc915d2
3 changed files with 277 additions and 104 deletions
Download Patch File Download Diff File Expand all files Collapse all files

194
inputs.cc
View File

@@ -7,6 +7,7 @@
#include <limits.h> #include <limits.h>
#include "common.h" #include "common.h"
#include "main.hh" #include "main.hh"
#include "tiles.hh"
int loadClutter(char *filename, double radius, struct site tx) int loadClutter(char *filename, double radius, struct site tx)
@@ -215,15 +216,15 @@ int loadLIDAR(char *filenames)
{ {
char *filename; char *filename;
char *files[100]; // 10x10 tiles char *files[100]; // 10x10 tiles
int x, y, indx = 0, fc = 0, hoffset = 0, voffset = 0, pos, int indx = 0, fc = 0, hoffset = 0, voffset = 0, pos, success;
dem_alloced = 0; double xll, yll, xur, yur, cellsize, avgCellsize = 0, smCellsize = 0;
double xll, yll, xur, yur, cellsize, avgCellsize;
char found, free_page = 0, jline[20], lid_file[255], 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 line[25000];
char * pch; char * pch;
double TO_DEG = (180 / PI); double TO_DEG = (180 / PI);
FILE *fd; FILE *fd;
tile_t *tiles;
// test for multiple files // test for multiple files
filename = strtok(filenames, " ,"); filename = strtok(filenames, " ,");
@@ -233,110 +234,141 @@ int loadLIDAR(char *filenames)
fc++; fc++;
} }
while (indx < fc) { /* Allocate the tile array */
if( (fd = fopen(files[indx], "rb")) == NULL ) tiles = (tile_t*)calloc(fc,sizeof(tile_t));
return errno;
if (fgets(line, 255, fd) != NULL) { /* Load each tile in turn */
pch = strtok (line," "); for (indx = 0; indx < fc; indx++) {
pch = strtok (NULL, " ");
width = atoi(pch); // ncols /* 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) { if (debug) {
fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, width); fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, tiles[indx].width);
fflush(stderr); fflush(stderr);
} }
if (fgets(line, 255, fd) != NULL) // Increase the average cell size
height = atoi(pch); // nrows avgCellsize += tiles[indx].cellsize;
// Update the smallest cell size
if (smCellsize == 0 || tiles[indx].cellsize < smCellsize) {
smCellsize = tiles[indx].cellsize;
}
if (!dem_alloced) { }
//Reduce MAXPAGES to increase speed
MAXPAGES=fc; /* Iterate through all tiles to find the largest x/y dimension
if(width>height){ to use as the global IPPD value. We do this here so that we are
IPPD = width; sure to allocate enough memory for the remainder of the calculations */
}else{ int dimension_max = -1;
IPPD = height; 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 // add fudge as reprojected tiles sometimes vary by a pixel or ten
IPPD+=50; IPPD += 50;
ARRAYSIZE = (MAXPAGES * IPPD)+50; ARRAYSIZE = (MAXPAGES * IPPD) + 50;
do_allocs(); do_allocs();
dem_alloced = 1; // reset the IPPD after allocations
} IPPD -= 50;
} /* Load the data into the global dem array */
if (fgets(line, 255, fd) != NULL) { // for (size_t indx = 0; indx < fc; indx++) {
sscanf(pch, "%lf", &xll); // xll // 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 (fgets(line, 255, fd) != NULL) { // //Not sure why we set these here...
sscanf(pch, "%lf", &yll); // yll // 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 (fgets(line, 255, fd) != NULL) // if (max_north == -90)
sscanf(pch, "%lf", &cellsize); // max_north = dem[indx].max_north;
avgCellsize=avgCellsize+cellsize; // else if (dem[indx].max_north > max_north)
// max_north = dem[indx].max_north;
/*if(cellsize>=0.5){ // 50cm LIDAR? // if (min_north == 90)
// compute xur and yur with inverse haversine if cellsize in *metres* // min_north = dem[indx].min_north;
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) // else if (dem[indx].min_north < min_north)
eastoffset = xur; // min_north = dem[indx].min_north;
if (xll < westoffset)
westoffset = xll;
if (debug) // if (dem[indx].max_west > max_west)
fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur); // 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;
// }
// }
// Greenwich straddling hack // if (min_west == 360) {
if (xll <= 0 && xur > 0) { // min_west = dem[indx].min_west;
xll = (xur - xll); // full width // } else {
xur = 0.0; // budge it along so it's west of greenwich // if (fabs(dem[indx].min_west - min_west) < 180.0) {
delta = eastoffset; // add to Tx longitude later // if (dem[indx].min_west < min_west)
} else { // min_west = dem[indx].min_west;
// Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;) // } else {
if (xll >= 0) // if (dem[indx].min_west > min_west)
xll = 360-xll; // min_west = dem[indx].min_west;
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);
// // /* 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;
// // }
// // }
fgets(line, 255, fd); // NODATA // }
pos = ftell(fd);
// tile 0 [x| ] 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) 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", 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; ippd=IPPD;
height = (unsigned)((max_north-min_north) / cellsize); height = (unsigned)((max_north-min_north) / smCellsize);
width = (unsigned)((max_west-min_west) / cellsize); width = (unsigned)((max_west-min_west) / smCellsize);
fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height);
if (debug) 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); 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);

113
tiles.cc Normal file
View File

@@ -0,0 +1,113 @@
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <string.h>
#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;
}

28
tiles.hh Normal file
View File

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