2.8 LIDAR improvements

This commit is contained in:
alex
2016-06-14 20:36:47 +01:00
parent 7845afd2f0
commit adac44f2f3
13 changed files with 372 additions and 5831 deletions

360
inputs.cc
View File

@@ -6,29 +6,30 @@
#include "common.h"
#include "main.hh"
void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx, double n, double e, double s, double west)
void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx,
double n, double e, double s, double west)
{
int x=0,y=0,reads=0;
int x = 0, y = 0, reads = 0;
char line[150000];
char * pch;
char *pch;
// use offsets to determine middle lat/lon for 5 x 10k tiles
// TALL
if(voffset==0 && h==10000){
if (voffset == 0 && h == 10000) {
s = (s+n)/2;
h=5000;
h = 5000;
}
if(voffset==5000 && h==10000){
if (voffset == 5000 && h == 10000) {
n = (s+n)/2;
}
// WIDE
if(hoffset==0 && w==10000){
if (hoffset == 0 && w == 10000) {
e = (e+west)/2;
w=5000;
w = 5000;
}
if(hoffset==5000 && w==10000){
if (hoffset == 5000 && w == 10000) {
west = (e+west)/2;
w=5000;
w = 5000;
}
dem[indx].max_north=n;
@@ -36,270 +37,271 @@ void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx, doubl
dem[indx].min_north=s;
dem[indx].max_west=west;
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)
if (max_west == -1) {
max_west = dem[indx].max_west;
else {
} else {
if (abs(dem[indx].max_west - max_west) < 180) {
if (dem[indx].max_west > max_west)
max_west = dem[indx].max_west;
}
else {
} else {
if (dem[indx].max_west < max_west)
max_west = dem[indx].max_west;
}
}
if (min_west == 360)
if (min_west == 360) {
min_west = dem[indx].min_west;
else {
} 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 {
} 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, 150000, fd) != NULL) {
// do nothing until n rows have passed
if(y<voffset || voffset==0){
pch = strtok (line, " ");
x = w-1;
if (fgets(line, 150000, fd) != NULL) {
// do nothing until n rows have passed
if (y < voffset || voffset == 0) {
pch = strtok(line, " ");
//dummy reads until we reach offset
// for 5000 offset, width must be 10e3
for(n=0;n<hoffset;n++){
pch = strtok(NULL, " ");
}
//dummy reads until we reach offset
// for 5000 offset, width must be 10e3
for (n = 0; n < hoffset; n++)
pch = strtok(NULL, " ");
while(pch != NULL && x > -1){
if(atoi(pch)<-999){
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
} //voffset
}else{
fprintf(stdout,"LIDAR error @ x %d y %d indx %d\n",x,y,indx);
}//if
}//for
while (pch != NULL && x > -1) {
if (atoi(pch) < -999)
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
}//voffset
} else {
fprintf(stdout, "LIDAR error @ x %d y %d indx %d\n",
x, y, indx);
}//if
}//for
}
int loadLIDAR(char *filenames)
{
/* This function reads either 9 LIDAR tiles of n rows and n columns in ASCII grid format OR a single super tile composed of 2 or more tiles.
The tile must have WGS84 bounds in the header in the order: WEST,SOUTH,EAST,NORTH
ncols 5000
nrows 5000
xllcorner -2.291359
yllcorner 51.788295
xurcorner -2.146674
yurcorner 51.878474
cellsize 2
NODATA_value -9999
Tiles must be entered in the format -lid tile1.asc,tile2.asc,tile3.asc
*/
/*
* This function reads either 9 LIDAR tiles of n rows and n columns
* in ASCII grid format OR a single super tile composed of 2 or more
* tiles. The tile must have WGS84 bounds in the header in the order:
* WEST,SOUTH,EAST,NORTH
* ncols 5000
* nrows 5000
* xllcorner -2.291359
* yllcorner 51.788295
* xurcorner -2.146674
* yurcorner 51.878474
* cellsize 2
* NODATA_value -9999
*
* Tiles must be entered in the format
*
* -lid tile1.asc,tile2.asc,tile3.asc
*/
char *filename;
char *files[4]; // 4 tiles
int x, y, cellsize,indx=0,fc=0,hoffset=0,voffset=0,pos;
double xll, yll, xur, yur;
int x, y, indx = 0, fc = 0, hoffset = 0, voffset = 0, pos,
dem_alloced = 0;
double xll, yll, xur, yur, cellsize;
char found, free_page = 0, jline[20], lid_file[255],
path_plus_name[255], *junk = NULL;
path_plus_name[255], *junk = NULL;
char line[50000];
char * pch;
FILE *fd;
// test for multiple files
filename = strtok(filenames, " ,");
while (filename != NULL)
{
while (filename != NULL) {
files[fc] = filename;
filename = strtok(NULL, " ,");
fc++;
}
while (indx<fc) {
while (indx < fc) {
fd = fopen(files[indx], "rb");
if (fd != NULL) {
if (debug) {
fprintf(stdout,"Loading \"%s\" into page %d...\n",files[indx], indx);
fprintf(stdout, "Loading \"%s\" into page %d...\n", files[indx], indx);
fflush(stdout);
}
if (fgets(line, 19, fd) != NULL) {
pch = strtok (line," ");
pch = strtok (NULL, " ");
width=atoi(pch);
pch = strtok (line," ");
pch = strtok (NULL, " ");
width = atoi(pch);
if (!dem_alloced) {
IPPD = width;
ARRAYSIZE = (MAXPAGES * IPPD) + 10;
do_allocs();
dem_alloced = 1;
}
}
if (fgets(line, 19, fd) != NULL) {
height=atoi(pch);
}
fgets(line, 24, fd); //
if (fgets(line, 19, fd) != NULL)
height = atoi(pch);
fgets(line, 24, fd);
if (fgets(line, 24, fd) != NULL) {
//xll=atof(pch);
sscanf(pch, "%lf", &xll);
//xll=atof(pch);
sscanf(pch, "%lf", &xll);
}
fgets(line, 24, fd); //
fgets(line, 24, fd);
if (fgets(line, 24, fd) != NULL) {
//yll=atof(pch);
sscanf(pch, "%lf", &yll);
//yll=atof(pch);
sscanf(pch, "%lf", &yll);
}
fgets(line, 24, fd); //
fgets(line, 24, fd);
if (fgets(line, 24, fd) != NULL) {
//xur=atof(pch);
sscanf(pch, "%lf", &xur);
//xur=atof(pch);
sscanf(pch, "%lf", &xur);
}
fgets(line, 24, fd); //
fgets(line, 24, fd);
if (fgets(line, 24, fd) != NULL) {
//yur=atof(pch);
sscanf(pch, "%lf", &yur);
}
fgets(line, 15, fd); //
if (fgets(line, 15, fd) != NULL) {
cellsize=atoi(pch);
//yur=atof(pch);
sscanf(pch, "%lf", &yur);
}
fgets(line, 15, fd);
if (fgets(line, 15, fd) != NULL)
cellsize = strtod(pch, NULL);
// LIDAR 10m @ 10800 PPD
if (cellsize == 10.0)
MAXRAD = 30;
// LIDAR 2m @ 54000 PPD
if(cellsize==2){
ippd=5000;
MAXRAD=15;
}
if (cellsize == 2.0)
MAXRAD = 15;
// LIDAR 1m @ 108000 PPD!
if(cellsize==1){
ippd=10000;
MAXRAD=10;
}
if (cellsize == 1.0)
MAXRAD = 10;
if(xur<eastoffset)
eastoffset=xur;
if(xll>westoffset)
westoffset=xll;
if (xur < eastoffset)
eastoffset = xur;
if (xll > westoffset)
westoffset = xll;
if (debug)
fprintf(stdout, "PRE yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta);
if(debug){
fprintf(stdout,"PRE yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n",yll,yur,xur,xll,delta);
}
// Greenwich straddling hack
if(xll < 0 && xur > 0){
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{
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 (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(stdout, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta);
}
if(debug){
fprintf(stdout,"POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n",yll,yur,xur,xll,delta);
}
if (yll < min_north)
min_north=yll;
min_north = yll;
if (yur > max_north)
max_north=yur;
max_north = yur;
fgets(line, 30, fd); // NODATA
pos=ftell(fd);
pos = ftell(fd);
// tile 0 [x| ]
if(debug){
fprintf(stdout,"readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n",0,0,height,width,indx,yur,xur,yll,xll);
}
readLIDAR(fd,0,0,height,width,indx,yur,xur,yll,xll);
if (debug)
fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", 0, 0, height, width, indx, yur, xur, yll, xll);
readLIDAR(fd, 0, 0, height, width, indx, yur, xur, yll,
xll);
//rewind
fseek(fd,pos,SEEK_SET);
fseek(fd, pos, SEEK_SET);
// tile 1 [ |x]
if(width==10000){
if (width == 2000 ||
width == 10000 ||
width == 10082) {
indx++;
if(debug){
fprintf(stdout,"readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n",5000,0,height,width,indx,yur,xur,yll,xll);
}
readLIDAR(fd,5000,0,height,width,indx,yur,xur,yll,xll);
}
if (debug)
fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", width / 2, 0, height, width, indx, yur, xur, yll, xll);
readLIDAR(fd, width / 2, 0, height, width,
indx, yur, xur, yll, xll);
}
//rewind
fseek(fd,pos,SEEK_SET);
fseek(fd, pos, SEEK_SET);
// tile 2 [x | ]
if(height==10000){
if (height == 2000 ||
height == 10000 ||
height == 10082) {
indx++;
if(debug){
fprintf(stdout,"readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n",0,5000,height,width,indx,yur,xur,yll,xll);
}
readLIDAR(fd,0,5000,height,width,indx,yur,xur,yll,xll);
}
if (debug)
fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", 0, height / 2, height, width, indx, yur, xur, yll, xll);
readLIDAR(fd, 0, height / 2, height, width,
indx, yur, xur, yll, xll);
}
//rewind
fseek(fd,pos,SEEK_SET);
fseek(fd, pos, SEEK_SET);
// tile 3 [ |x]
if(width==10000 && height==10000){
if ((width == 2000 && height == 2000) ||
(width == 10000 && height == 10000) ||
(width == 10082 && height == 10082)) {
indx++;
if(debug){
fprintf(stdout,"readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n",5000,5000,height,width,indx,yur,xur,yll,xll);
}
readLIDAR(fd,5000,5000,height,width,indx,yur,xur,yll,xll);
if (debug)
fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", width / 2, height / 2, height, width, indx, yur, xur, yll, xll);
readLIDAR(fd, width / 2, height / 2, height,
width, indx, yur, xur, yll,
xll);
}
fclose(fd);
fprintf(stdout,"LIDAR LOADED %d x %d\n",width,height);
} // if (fd != NULL)
else
if (debug)
fprintf(stdout, "LIDAR LOADED %d x %d\n", width, height);
} else {
return -1;
indx++;
}
indx++;
} // filename(s)
return 0;
}
int LoadSDF_SDF(char *name)