Numerous small fixes following code review:
(1) Whitespace use matched better to surrounding code (2) While(1) in TGArray::rectify_heights replaced with do/while loop (3) Incorrect interpolation for 2 points fixed in altitude_from_grid (4) exit(1) replaced with return EXIT_FAILURE in cliff-decode.cxx
This commit is contained in:
parent
dfe81ce9fa
commit
3421ce2018
5 changed files with 404 additions and 352 deletions
|
@ -67,14 +67,14 @@ bool TGArray::open( const string& file_base ) {
|
|||
|
||||
array_in = gzopen( array_name.c_str(), "rb" );
|
||||
if (array_in == NULL) {
|
||||
// try unrectified
|
||||
array_name = file_base + ".arr.gz";
|
||||
array_in = gzopen(array_name.c_str(), "rb");
|
||||
if (array_in == NULL) {
|
||||
return false;
|
||||
} else {
|
||||
rectified = false;
|
||||
}
|
||||
// try unrectified
|
||||
array_name = file_base + ".arr.gz";
|
||||
array_in = gzopen(array_name.c_str(), "rb");
|
||||
if (array_in == NULL) {
|
||||
return false;
|
||||
} else {
|
||||
rectified = false;
|
||||
}
|
||||
}
|
||||
|
||||
SG_LOG(SG_GENERAL,SG_DEBUG,"Loaded height array " << array_name);
|
||||
|
@ -121,34 +121,34 @@ TGArray::close() {
|
|||
|
||||
void TGArray::load_cliffs(const string & height_base)
|
||||
{
|
||||
//Get the directory so we can list the children
|
||||
tgPolygon poly; //actually a contour but whatever...
|
||||
int total_contours_read = 0;
|
||||
SGPath b(height_base);
|
||||
simgear::Dir d(b.dir());
|
||||
simgear::PathList files = d.children(simgear::Dir::TYPE_FILE);
|
||||
BOOST_FOREACH(const SGPath& p, files) {
|
||||
if (p.file_base() != b.file_base()) {
|
||||
continue;
|
||||
}
|
||||
|
||||
string lext = p.lower_extension();
|
||||
if (lext == "cliffs") {
|
||||
gzFile fp = gzopen( p.c_str(), "rb" );
|
||||
unsigned int count;
|
||||
sgReadUInt( fp, &count );
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " Load " << count << " contours from " << p.realpath() );
|
||||
|
||||
for ( unsigned int i=0; i<count; i++ ) {
|
||||
poly.LoadFromGzFile( fp );
|
||||
if ( poly.Contours()==1 ) { //should always have one contour
|
||||
cliffs_list.push_back(poly.GetContour(0));
|
||||
} else {
|
||||
SG_LOG( SG_GENERAL, SG_WARN, " Found " << poly.Contours() << " contours in " << p.realpath() );
|
||||
//Get the directory so we can list the children
|
||||
tgPolygon poly; //actually a contour but whatever...
|
||||
int total_contours_read = 0;
|
||||
SGPath b(height_base);
|
||||
simgear::Dir d(b.dir());
|
||||
simgear::PathList files = d.children(simgear::Dir::TYPE_FILE);
|
||||
for (const SGPath& p: files) {
|
||||
if (p.file_base() != b.file_base()) {
|
||||
continue;
|
||||
}
|
||||
|
||||
string lext = p.lower_extension();
|
||||
if (lext == "cliffs") {
|
||||
gzFile fp = gzopen( p.c_str(), "rb" );
|
||||
unsigned int count;
|
||||
sgReadUInt( fp, &count );
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " Load " << count << " contours from " << p.realpath() );
|
||||
|
||||
for ( unsigned int i=0; i<count; i++ ) {
|
||||
poly.LoadFromGzFile( fp );
|
||||
if ( poly.Contours()==1 ) { //should always have one contour
|
||||
cliffs_list.push_back(poly.GetContour(0));
|
||||
} else {
|
||||
SG_LOG( SG_GENERAL, SG_WARN, " Found " << poly.Contours() << " contours in " << p.realpath() );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
@ -408,7 +408,7 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
|
|||
|
||||
for ( int row = 0; row < rows; row++ ) {
|
||||
for ( int col = 0; col < cols; col++ ) {
|
||||
SGGeod p1 = SGGeod::fromDeg( (originx + col * col_step)/3600.0, (originy + row * row_step)/3600.0 );
|
||||
SGGeod p1 = SGGeod::fromDeg( (originx + col * col_step)/3600.0, (originy + row * row_step)/3600.0 );
|
||||
double dist = SGGeodesy::distanceM( p0, p1 );
|
||||
double elev = get_array_elev(col, row);
|
||||
if ( dist < mindist && elev > -9000 ) {
|
||||
|
@ -425,61 +425,68 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
|
|||
}
|
||||
}
|
||||
|
||||
//Find and remember all points that are bad because they are
|
||||
//too close to a cliff
|
||||
std::vector<int> TGArray::collect_bad_points(const double bad_zone) {
|
||||
//Find and remember all points that are bad because they are
|
||||
//too close to a cliff
|
||||
std::vector<int> bad_points; //local to avoid multi-thread issues
|
||||
for(int horiz=0;horiz<cols;horiz++) {
|
||||
double lon = (originx + col_step*horiz)/3600;
|
||||
for(int vert=0;vert<rows;vert++) {
|
||||
double lat = (originy + row_step*vert)/3600;
|
||||
if(is_near_cliff(lon,lat,bad_zone)) {
|
||||
bad_points.push_back(horiz+vert*cols);
|
||||
|
||||
std::vector<int> bad_points; //local to avoid multi-thread issues
|
||||
|
||||
for( int horiz=0;horiz<cols;horiz++ ) {
|
||||
double lon = (originx + col_step*horiz)/3600;
|
||||
for( int vert=0;vert<rows;vert++ ) {
|
||||
double lat = (originy + row_step*vert)/3600;
|
||||
if( is_near_cliff(lon,lat,bad_zone) ) {
|
||||
bad_points.push_back(horiz+vert*cols);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return bad_points;
|
||||
|
||||
return bad_points;
|
||||
}
|
||||
|
||||
// Check to see if the specified grid point is bad
|
||||
bool TGArray::is_bad_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const {
|
||||
int grididx;
|
||||
grididx = xgrid+ygrid*cols;
|
||||
auto result = std::find(std::begin(bad_points),std::end(bad_points),grididx);
|
||||
if (result != std::end(bad_points)) return true;
|
||||
return false;
|
||||
}
|
||||
int grididx;
|
||||
grididx = xgrid+ygrid*cols;
|
||||
auto result = std::find( std::begin(bad_points),std::end(bad_points),grididx );
|
||||
if ( result != std::end(bad_points) ) return true;
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
//This may collide with other threads, but as they will both be writing
|
||||
//the correct height, this is harmless.
|
||||
void TGArray::rectify_heights(const double bad_zone) {
|
||||
double new_ht;
|
||||
std::vector<int> rectified,bad_points;
|
||||
bad_points = collect_bad_points(bad_zone);
|
||||
while(1) {
|
||||
for (auto pt : bad_points) {
|
||||
int ygrid = pt/cols;
|
||||
int xgrid = pt - ygrid*cols;
|
||||
new_ht = rectify_point(xgrid,ygrid,bad_points);
|
||||
if (new_ht > -9999) {
|
||||
rectified.push_back(pt);
|
||||
set_array_elev(xgrid,ygrid,(int) new_ht);
|
||||
}
|
||||
}
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Rectified " << rectified.size() << " points ");
|
||||
if(rectified.size()>0) {
|
||||
for(auto r : rectified) {
|
||||
bad_points.erase(std::remove(std::begin(bad_points),std::end(bad_points),r));
|
||||
}
|
||||
rectified.clear();
|
||||
} else {
|
||||
if(bad_points.size() > 0) {
|
||||
void TGArray::rectify_heights( const double bad_zone ) {
|
||||
double new_ht;
|
||||
std::vector<int> rectified,bad_points;
|
||||
int total_rectified;
|
||||
bad_points = collect_bad_points( bad_zone );
|
||||
|
||||
do {
|
||||
for ( auto pt : bad_points ) {
|
||||
int ygrid = pt/cols;
|
||||
int xgrid = pt - ygrid*cols;
|
||||
new_ht = rectify_point( xgrid,ygrid,bad_points );
|
||||
if (new_ht > -9999) {
|
||||
rectified.push_back(pt);
|
||||
set_array_elev( xgrid,ygrid,(int) new_ht );
|
||||
}
|
||||
}
|
||||
total_rectified = rectified.size();
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Rectified " << total_rectified << " points ");
|
||||
|
||||
if( total_rectified > 0 ) {
|
||||
for( auto r : rectified ) {
|
||||
bad_points.erase( std::remove( std::begin(bad_points), std::end(bad_points),r) );
|
||||
}
|
||||
rectified.clear();
|
||||
}
|
||||
} while ( total_rectified > 0 );
|
||||
|
||||
if( bad_points.size() > 0 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Failed to rectify " << bad_points.size() << " points");
|
||||
}
|
||||
break; // Cant do any more
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* If we have cliffs, it is possible that a grid point will be too close
|
||||
|
@ -496,67 +503,82 @@ through the three known points.
|
|||
* * *
|
||||
|
||||
TODO: Handle points on the boundaries. */
|
||||
|
||||
double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const {
|
||||
//xgrid: grid units horizontally
|
||||
//ygrid: grid units vertically
|
||||
//Loop over corner points, if no points available, give up
|
||||
int corners[4][2]; //possible corners
|
||||
int final_pts[3][2]; // rectangle corners
|
||||
int pt_cnt = 0;
|
||||
double centre_long, centre_lat;
|
||||
double cliff_error = col_step; //Assume row step, col step the same
|
||||
int original_height = get_array_elev(xgrid,ygrid);
|
||||
centre_long = (originx + col_step*xgrid)/3600;
|
||||
centre_lat = (originy + row_step*ygrid)/3600;
|
||||
for (int horiz = -1; horiz <= 1; horiz+=2) {
|
||||
if (xgrid + horiz >= cols || xgrid + horiz < 0) continue; //edge of bucket
|
||||
double test_long = centre_long + (col_step*horiz)/3600;
|
||||
for (int vert = -1; vert <= 1; vert+=2) {
|
||||
if (ygrid + vert >= rows || ygrid + vert < 0) continue; //edge of bucket
|
||||
double test_lat = centre_lat + (row_step*vert)/3600;
|
||||
if (!is_bad_point(xgrid+horiz,ygrid+vert,bad_points) && //can trust height
|
||||
check_points(test_long,test_lat,centre_long,centre_lat)) { //same side
|
||||
corners[pt_cnt][0] = horiz;
|
||||
corners[pt_cnt][1] = vert;
|
||||
pt_cnt++;
|
||||
//xgrid: grid units horizontally
|
||||
//ygrid: grid units vertically
|
||||
//Loop over corner points, if no points available, give up
|
||||
int corners[4][2]; //possible corners
|
||||
int final_pts[3][2]; // rectangle corners
|
||||
int pt_cnt = 0;
|
||||
double centre_long, centre_lat;
|
||||
double cliff_error = col_step; //Assume row step, col step the same
|
||||
int original_height = get_array_elev(xgrid,ygrid);
|
||||
centre_long = (originx + col_step*xgrid)/3600;
|
||||
centre_lat = (originy + row_step*ygrid)/3600;
|
||||
|
||||
for ( int horiz = -1; horiz <= 1; horiz+=2 ) {
|
||||
if (xgrid + horiz >= cols || xgrid + horiz < 0) continue; //edge of bucket
|
||||
|
||||
double test_long = centre_long + (col_step*horiz)/3600;
|
||||
for ( int vert = -1; vert <= 1; vert+=2 ) {
|
||||
if (ygrid + vert >= rows || ygrid + vert < 0) continue; //edge of bucket
|
||||
|
||||
double test_lat = centre_lat + (row_step*vert)/3600;
|
||||
if ( !is_bad_point( xgrid+horiz,ygrid+vert,bad_points ) && //can trust height
|
||||
check_points( test_long,test_lat,centre_long,centre_lat ) ) { //same side
|
||||
|
||||
corners[pt_cnt][0] = horiz;
|
||||
corners[pt_cnt][1] = vert;
|
||||
pt_cnt++;
|
||||
}
|
||||
}
|
||||
}
|
||||
} // end of search for corners
|
||||
if (pt_cnt == 0) return -9999; // no corners found
|
||||
// Find two points that form a rectangle with a corner
|
||||
int pt;
|
||||
double height = 0;
|
||||
for (pt = 0; pt < pt_cnt; pt++) {
|
||||
if (!is_bad_point(xgrid+corners[pt][0],ygrid,bad_points) &&
|
||||
!is_bad_point(xgrid, ygrid+corners[pt][1],bad_points)) {
|
||||
double test_horiz = centre_long + corners[pt][0]*col_step/3600;
|
||||
double test_vert = centre_lat + corners[pt][1]*row_step/3600;
|
||||
if (check_points(test_horiz,centre_lat,centre_long,centre_lat) &&
|
||||
check_points(centre_long,test_vert,centre_long,centre_lat)) break;
|
||||
} // end of search for corners
|
||||
|
||||
if (pt_cnt == 0) return -9999; // no corners found
|
||||
|
||||
// Find two points that form a rectangle with a corner
|
||||
int pt;
|
||||
double height = 0;
|
||||
for ( pt = 0; pt < pt_cnt; pt++ ) {
|
||||
|
||||
if ( !is_bad_point( xgrid+corners[pt][0],ygrid,bad_points ) &&
|
||||
!is_bad_point( xgrid, ygrid+corners[pt][1],bad_points ) ) {
|
||||
|
||||
double test_horiz = centre_long + corners[pt][0]*col_step/3600;
|
||||
double test_vert = centre_lat + corners[pt][1]*row_step/3600;
|
||||
|
||||
if ( check_points( test_horiz,centre_lat,centre_long,centre_lat ) &&
|
||||
check_points( centre_long,test_vert,centre_long,centre_lat ) ) break;
|
||||
}
|
||||
}
|
||||
|
||||
if (pt == pt_cnt) { // perhaps we have a concave cliff, just take the
|
||||
// average of the known points
|
||||
double totht = 0;
|
||||
for(int pti = 0; pti <pt_cnt; pti++) {
|
||||
totht = totht + get_array_elev(xgrid+corners[pti][0],ygrid+corners[pti][1]);
|
||||
}
|
||||
height = totht/pt_cnt;
|
||||
} else {
|
||||
if (pt == pt_cnt) {
|
||||
// perhaps we have a concave cliff, just take the
|
||||
// average of the known points
|
||||
double totht = 0;
|
||||
for( int pti = 0; pti <pt_cnt; pti++ ) {
|
||||
totht = totht + get_array_elev( xgrid+corners[pti][0],ygrid+corners[pti][1] );
|
||||
}
|
||||
|
||||
height = totht/pt_cnt;
|
||||
|
||||
} else {
|
||||
|
||||
// We have three points, calculate the height
|
||||
// Set anything very negative to zero
|
||||
double corner = get_array_elev(xgrid+corners[pt][0],ygrid+corners[pt][1]);
|
||||
double horiz = get_array_elev(xgrid,ygrid+corners[pt][1]);
|
||||
double vert = get_array_elev(xgrid+corners[pt][0],ygrid);
|
||||
if (corner < -9000) corner = 0;
|
||||
if (horiz < -9000) horiz = 0;
|
||||
if (vert < -9000) vert = 0;
|
||||
height = horiz + (vert - corner);
|
||||
}
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, xgrid << "," << ygrid << ": was " << original_height << " , now " << height);
|
||||
return height;
|
||||
// We have three points, calculate the height
|
||||
// Set anything very negative to zero
|
||||
double corner = get_array_elev( xgrid+corners[pt][0],ygrid+corners[pt][1] );
|
||||
double horiz = get_array_elev( xgrid,ygrid+corners[pt][1] );
|
||||
double vert = get_array_elev( xgrid+corners[pt][0],ygrid );
|
||||
if ( corner < -9000 ) corner = 0;
|
||||
if ( horiz < -9000 ) horiz = 0;
|
||||
if ( vert < -9000 ) vert = 0;
|
||||
height = horiz + ( vert - corner );
|
||||
}
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, xgrid << "," << ygrid << ": was " << original_height << " , now " << height);
|
||||
|
||||
return height;
|
||||
}
|
||||
|
||||
// return the current altitude based on grid data.
|
||||
|
@ -578,7 +600,7 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
|
|||
------
|
||||
|
||||
then calculate our end points
|
||||
*/
|
||||
*/
|
||||
|
||||
// Store in degrees for later
|
||||
double londeg = lon/3600;
|
||||
|
@ -589,8 +611,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
|
|||
xindex = (int)(xlocal);
|
||||
yindex = (int)(ylocal);
|
||||
|
||||
// printf("xindex = %d yindex = %d\n", xindex, yindex);
|
||||
|
||||
if ( xindex + 1 == cols ) {
|
||||
xindex--;
|
||||
}
|
||||
|
@ -601,7 +621,9 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
|
|||
|
||||
if ( (xindex < 0) || (xindex + 1 >= cols) ||
|
||||
(yindex < 0) || (yindex + 1 >= rows) ) {
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "WARNING: Attempt to interpolate value outside of array!!!" );
|
||||
|
||||
return -9999;
|
||||
}
|
||||
|
||||
|
@ -613,95 +635,101 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
|
|||
int corners[4][2];
|
||||
int ccnt = 0;
|
||||
int missing = -1; //the missing point when 3 available
|
||||
|
||||
double lon1 = (originx+(xindex*col_step))/3600;
|
||||
double lat1 = (originy+(yindex*row_step))/3600;
|
||||
double lon2 = lon1 + col_step/3600;
|
||||
double lat2 = lat1 + row_step/3600;
|
||||
if (check_points(lon1,lat1,londeg,latdeg)) {
|
||||
corners[ccnt][0] = xindex;
|
||||
corners[ccnt][1] = yindex;
|
||||
ccnt++;
|
||||
} else missing = 0;
|
||||
if (check_points(lon1,lat2,londeg,latdeg)) {
|
||||
corners[ccnt][0] = xindex;
|
||||
corners[ccnt][1] = yindex+1;
|
||||
ccnt++;
|
||||
} else missing = 1;
|
||||
if (check_points(lon2,lat2,londeg,latdeg)) {
|
||||
corners[ccnt][0] = xindex+1;
|
||||
corners[ccnt][1] = yindex+1;
|
||||
ccnt++;
|
||||
} else missing = 2;
|
||||
if (check_points(lon2,lat1,londeg,latdeg)) {
|
||||
corners[ccnt][0] = xindex+1;
|
||||
corners[ccnt][1] = yindex;
|
||||
ccnt++;
|
||||
} else missing = 3;
|
||||
|
||||
switch (ccnt) {
|
||||
case 3: //3 points are corners of a rectangle
|
||||
if ( check_points(lon1,lat1,londeg,latdeg) ) {
|
||||
corners[ccnt][0] = xindex;
|
||||
corners[ccnt][1] = yindex;
|
||||
ccnt++;
|
||||
} else missing = 0;
|
||||
if ( check_points(lon1,lat2,londeg,latdeg) ) {
|
||||
corners[ccnt][0] = xindex;
|
||||
corners[ccnt][1] = yindex+1;
|
||||
ccnt++;
|
||||
} else missing = 1;
|
||||
if ( check_points(lon2,lat2,londeg,latdeg) ) {
|
||||
corners[ccnt][0] = xindex+1;
|
||||
corners[ccnt][1] = yindex+1;
|
||||
ccnt++;
|
||||
} else missing = 2;
|
||||
if ( check_points(lon2,lat1,londeg,latdeg) ) {
|
||||
corners[ccnt][0] = xindex+1;
|
||||
corners[ccnt][1] = yindex;
|
||||
ccnt++;
|
||||
} else missing = 3;
|
||||
|
||||
switch (ccnt) {
|
||||
case 3: //3 points are corners of a rectangle
|
||||
// choose the points so that x2 is the right angle
|
||||
// and x1-x2 is the x arm of the triangle
|
||||
// dx,dy are the (positive) distances from the x1 corner
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "3 points, missing #" << missing);
|
||||
|
||||
dx = xlocal -xindex;
|
||||
dy = ylocal -yindex;
|
||||
switch (missing) {
|
||||
|
||||
switch ( missing ) {
|
||||
case 0: //SW corner missing
|
||||
x1 = corners[0][0];
|
||||
y1 = corners[0][1];
|
||||
x1 = corners[0][0];
|
||||
y1 = corners[0][1];
|
||||
|
||||
x2 = corners[1][0];
|
||||
y2 = corners[1][1];
|
||||
x2 = corners[1][0];
|
||||
y2 = corners[1][1];
|
||||
|
||||
x3 = corners[2][0];
|
||||
y3 = corners[2][1];
|
||||
x3 = corners[2][0];
|
||||
y3 = corners[2][1];
|
||||
|
||||
dy = 1 - dy;
|
||||
break;
|
||||
dy = 1 - dy;
|
||||
break;
|
||||
case 1: //NW corner missing
|
||||
x1 = corners[0][0];
|
||||
y1 = corners[0][1];
|
||||
x1 = corners[0][0];
|
||||
y1 = corners[0][1];
|
||||
|
||||
x2 = corners[2][0];
|
||||
y2 = corners[2][1];
|
||||
x2 = corners[2][0];
|
||||
y2 = corners[2][1];
|
||||
|
||||
x3 = corners[1][0];
|
||||
y3 = corners[1][1];
|
||||
x3 = corners[1][0];
|
||||
y3 = corners[1][1];
|
||||
|
||||
break;
|
||||
break;
|
||||
case 2: //NE corner missing
|
||||
x1 = corners[2][0];
|
||||
y1 = corners[2][1];
|
||||
x1 = corners[2][0];
|
||||
y1 = corners[2][1];
|
||||
|
||||
x2 = corners[0][0];
|
||||
y2 = corners[0][1];
|
||||
x2 = corners[0][0];
|
||||
y2 = corners[0][1];
|
||||
|
||||
x3 = corners[1][0];
|
||||
y3 = corners[1][1];
|
||||
x3 = corners[1][0];
|
||||
y3 = corners[1][1];
|
||||
|
||||
dx = 1 - dx; //x1 is SE corner
|
||||
break;
|
||||
dx = 1 - dx; //x1 is SE corner
|
||||
break;
|
||||
|
||||
case 3: //SE corner missing
|
||||
x1 = corners[2][0];
|
||||
y1 = corners[2][1];
|
||||
x1 = corners[2][0];
|
||||
y1 = corners[2][1];
|
||||
|
||||
x2 = corners[1][0];
|
||||
y2 = corners[1][1];
|
||||
x2 = corners[1][0];
|
||||
y2 = corners[1][1];
|
||||
|
||||
x3 = corners[0][0];
|
||||
y3 = corners[0][1];
|
||||
x3 = corners[0][0];
|
||||
y3 = corners[0][1];
|
||||
|
||||
dx = 1 - dx; //x1 is NE corner
|
||||
dy = 1 - dy;
|
||||
break;
|
||||
dx = 1 - dx; //x1 is NE corner
|
||||
dy = 1 - dy;
|
||||
break;
|
||||
|
||||
}
|
||||
// Now do the calcs on the triangle
|
||||
// We interpolate on height along x1-x2 and
|
||||
// x1 - x3. Then interpolate between these
|
||||
// two points along y.
|
||||
|
||||
z1 = get_array_elev(x1,y1);
|
||||
z2 = get_array_elev(x2,y2);
|
||||
z3 = get_array_elev(x3,y3);
|
||||
|
@ -709,13 +737,14 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
|
|||
zB = dx * (z3 - z1) + z1;
|
||||
|
||||
if ( dx > SG_EPSILON ) {
|
||||
elev = dy * (zB - zA) / dx + zA;
|
||||
elev = dy * (zB - zA) / dx + zA;
|
||||
} else {
|
||||
elev = zA;
|
||||
elev = zA;
|
||||
}
|
||||
|
||||
break;
|
||||
case 2: //project onto line connecting two points
|
||||
|
||||
case 2: //project onto line connecting two points
|
||||
x1 = corners[0][0];
|
||||
y1 = corners[0][1];
|
||||
z1 = get_array_elev(x1,y1);
|
||||
|
@ -729,119 +758,135 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
|
|||
dx = xlocal - x1;
|
||||
dy = ylocal - y1;
|
||||
if (x1==x2) {
|
||||
elev = z1+dy*(z2-z1);
|
||||
elev = z1+dy*(z2-z1)/(y2-y1);
|
||||
}
|
||||
else if (y1==y2) {
|
||||
elev = z1+dx*(z2-z1);
|
||||
elev = z1+dx*(z2-z1)/(x2-x1);
|
||||
}
|
||||
else { //diagonal: project onto 45 degree line
|
||||
int comp1 = x2-x1;
|
||||
int comp2 = y2-y1;
|
||||
double dotprod = (dx*comp1 + dy*comp2)/sqrt(2);
|
||||
double projlen = sqrt(dx*dx+dy*dy)*dotprod;
|
||||
elev = (z2-z1)*projlen/sqrt(2);
|
||||
int comp1 = x2-x1;
|
||||
int comp2 = y2-y1;
|
||||
double dotprod = (dx*comp1 + dy*comp2)/sqrt(2);
|
||||
double projlen = sqrt(dx*dx+dy*dy)*dotprod;
|
||||
elev = (z2-z1)*projlen/sqrt(2);
|
||||
}
|
||||
break;
|
||||
case 1: //only one point found
|
||||
elev = get_array_elev(corners[0][0],corners[0][1]);
|
||||
break;
|
||||
case 0: // all points on wrong side, fall through to normal calc
|
||||
|
||||
case 1: //only one point found
|
||||
elev = get_array_elev( corners[0][0],corners[0][1] );
|
||||
break;
|
||||
|
||||
case 0: // all points on wrong side, fall through to normal calc
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_WARN, "All elevation grid points on wrong side of cliff for " << std::setprecision(10) << londeg << "," << latdeg );
|
||||
SG_LOG(SG_GENERAL, SG_WARN, "Grid points ("<< std::setprecision(9) << lon1 << "," << lat1 << "),("<<lon2<<","<<lat2<<")");
|
||||
default: // all corners
|
||||
|
||||
default: // all corners
|
||||
dx = xlocal - xindex;
|
||||
dy = ylocal - yindex;
|
||||
|
||||
if ( dx > dy ) {
|
||||
// lower triangle
|
||||
// lower triangle
|
||||
|
||||
x1 = xindex;
|
||||
y1 = yindex;
|
||||
z1 = get_array_elev(x1, y1);
|
||||
x1 = xindex;
|
||||
y1 = yindex;
|
||||
z1 = get_array_elev(x1, y1);
|
||||
|
||||
x2 = xindex + 1;
|
||||
y2 = yindex;
|
||||
z2 = get_array_elev(x2, y2);
|
||||
x2 = xindex + 1;
|
||||
y2 = yindex;
|
||||
z2 = get_array_elev(x2, y2);
|
||||
|
||||
x3 = xindex + 1;
|
||||
y3 = yindex + 1;
|
||||
z3 = get_array_elev(x3, y3);
|
||||
x3 = xindex + 1;
|
||||
y3 = yindex + 1;
|
||||
z3 = get_array_elev(x3, y3);
|
||||
|
||||
if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) {
|
||||
// don't interpolate off a void
|
||||
return closest_nonvoid_elev( lon, lat );
|
||||
}
|
||||
if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) {
|
||||
// don't interpolate off a void
|
||||
return closest_nonvoid_elev( lon, lat );
|
||||
}
|
||||
|
||||
zA = dx * (z2 - z1) + z1;
|
||||
zB = dx * (z3 - z1) + z1;
|
||||
zA = dx * (z2 - z1) + z1;
|
||||
zB = dx * (z3 - z1) + z1;
|
||||
|
||||
if ( dx > SG_EPSILON ) {
|
||||
elev = dy * (zB - zA) / dx + zA;
|
||||
} else {
|
||||
elev = zA;
|
||||
}
|
||||
if ( dx > SG_EPSILON ) {
|
||||
elev = dy * (zB - zA) / dx + zA;
|
||||
} else {
|
||||
elev = zA;
|
||||
}
|
||||
} else {
|
||||
// upper triangle
|
||||
// upper triangle
|
||||
|
||||
x1 = xindex;
|
||||
y1 = yindex;
|
||||
z1 = get_array_elev(x1, y1);
|
||||
x1 = xindex;
|
||||
y1 = yindex;
|
||||
z1 = get_array_elev(x1, y1);
|
||||
|
||||
x2 = xindex;
|
||||
y2 = yindex + 1;
|
||||
z2 = get_array_elev(x2, y2);
|
||||
x2 = xindex;
|
||||
y2 = yindex + 1;
|
||||
z2 = get_array_elev(x2, y2);
|
||||
|
||||
x3 = xindex + 1;
|
||||
y3 = yindex + 1;
|
||||
z3 = get_array_elev(x3, y3);
|
||||
x3 = xindex + 1;
|
||||
y3 = yindex + 1;
|
||||
z3 = get_array_elev(x3, y3);
|
||||
|
||||
if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) {
|
||||
// don't interpolate off a void
|
||||
return closest_nonvoid_elev( lon, lat );
|
||||
}
|
||||
if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) {
|
||||
// don't interpolate off a void
|
||||
return closest_nonvoid_elev( lon, lat );
|
||||
}
|
||||
|
||||
zA = dy * (z2 - z1) + z1;
|
||||
zB = dy * (z3 - z1) + z1;
|
||||
zA = dy * (z2 - z1) + z1;
|
||||
zB = dy * (z3 - z1) + z1;
|
||||
|
||||
if ( dy > SG_EPSILON ) {
|
||||
elev = dx * (zB - zA) / dy + zA;
|
||||
} else {
|
||||
elev = zA;
|
||||
}
|
||||
if ( dy > SG_EPSILON ) {
|
||||
elev = dx * (zB - zA) / dy + zA;
|
||||
} else {
|
||||
elev = zA;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return elev;
|
||||
}
|
||||
|
||||
// Check that two points are on the same side of all cliff contours
|
||||
// Could speed up by checking bounding box first
|
||||
bool TGArray::check_points (const double lon1, const double lat1, const double lon2, const double lat2) const {
|
||||
if (cliffs_list.size()==0) return true;
|
||||
if (fabs(lon1-lon2)<SG_EPSILON && fabs(lat1-lat2)<SG_EPSILON) return true;
|
||||
SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
|
||||
SGGeod pt2 = SGGeod::fromDeg(lon2,lat2);
|
||||
bool same_side = true;
|
||||
for (int i=0;i<cliffs_list.size();i++) {
|
||||
bool check_result = cliffs_list[i].AreSameSide(pt1,pt2);
|
||||
if(!check_result) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Cliff " << i <<":" <<pt1 << " and " << pt2 << " on opposite sides");
|
||||
same_side = false;
|
||||
break;
|
||||
bool TGArray::check_points( const double lon1, const double lat1, const double lon2, const double lat2 ) const {
|
||||
|
||||
if ( cliffs_list.size()==0 ) return true;
|
||||
|
||||
if ( fabs(lon1-lon2)<SG_EPSILON && fabs(lat1-lat2)<SG_EPSILON ) return true;
|
||||
|
||||
SGGeod pt1 = SGGeod::fromDeg( lon1,lat1 );
|
||||
SGGeod pt2 = SGGeod::fromDeg( lon2,lat2 );
|
||||
bool same_side = true;
|
||||
|
||||
for ( int i=0;i<cliffs_list.size();i++ ) {
|
||||
bool check_result = cliffs_list[i].AreSameSide( pt1,pt2 );
|
||||
|
||||
if(!check_result) {
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Cliff " << i <<":" <<pt1 << " and " << pt2 << " on opposite sides");
|
||||
|
||||
same_side = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
return same_side;
|
||||
|
||||
return same_side;
|
||||
}
|
||||
|
||||
//Check that a point is more than given distance from any cliff
|
||||
//Could speed up by checking bounding box
|
||||
bool TGArray::is_near_cliff(const double lon1, const double lat1, const double bad_zone) const {
|
||||
if (cliffs_list.size()==0) return false;
|
||||
SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
|
||||
for (int i=0;i<cliffs_list.size();i++) {
|
||||
double dist = cliffs_list[i].MinDist(pt1);
|
||||
if (dist < bad_zone) return true;
|
||||
}
|
||||
return false;
|
||||
bool TGArray::is_near_cliff( const double lon1, const double lat1, const double bad_zone ) const {
|
||||
|
||||
if (cliffs_list.size()==0) return false;
|
||||
|
||||
SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
|
||||
|
||||
for ( int i=0;i<cliffs_list.size();i++ ) {
|
||||
double dist = cliffs_list[i].MinDist(pt1);
|
||||
if (dist < bad_zone) return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
TGArray::~TGArray( void )
|
||||
|
@ -875,10 +920,10 @@ void TGArray::set_array_elev( int col, int row, int val )
|
|||
|
||||
bool TGArray::is_open() const
|
||||
{
|
||||
if ( array_in != NULL ) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
if ( array_in != NULL ) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
@ -11,7 +11,7 @@
|
|||
#include "tg_misc.hxx"
|
||||
|
||||
tgPolygon tgChopper::Clip( const tgPolygon& subject,
|
||||
const std::string& type, SGBucket& b)
|
||||
const std::string& type, SGBucket& b )
|
||||
{
|
||||
tgPolygon base, result;
|
||||
|
||||
|
@ -21,8 +21,7 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
|
|||
base.AddNode( 0, b.get_corner( SG_BUCKET_NE ) );
|
||||
base.AddNode( 0, b.get_corner( SG_BUCKET_NW ) );
|
||||
|
||||
result = tgPolygon::Intersect( subject, base);
|
||||
// Debug: See if numbers of nodes have changed
|
||||
result = tgPolygon::Intersect( subject, base );
|
||||
if ( result.Contours() > 0 ) {
|
||||
|
||||
if ( subject.GetPreserve3D() ) {
|
||||
|
@ -68,7 +67,7 @@ void tgChopper::ClipRow( const tgPolygon& subject, const double& center_lat, con
|
|||
}
|
||||
|
||||
|
||||
void tgChopper::Add( const tgPolygon& subject, const std::string& type)
|
||||
void tgChopper::Add( const tgPolygon& subject, const std::string& type )
|
||||
{
|
||||
// bail out immediately if polygon is empty
|
||||
if ( subject.Contours() == 0 )
|
||||
|
@ -95,7 +94,7 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type)
|
|||
// We just have a single row - no need to intersect first
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " UN_CLIPPED row - center lat is " << b_min.get_center_lat() );
|
||||
|
||||
ClipRow( subject, b_min.get_center_lat(), type);
|
||||
ClipRow( subject, b_min.get_center_lat(), type );
|
||||
}
|
||||
else
|
||||
{
|
||||
|
@ -122,7 +121,7 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type)
|
|||
clip_row.AddNode( 0, SGGeod::fromDeg( 180.0, clip_top) );
|
||||
clip_row.AddNode( 0, SGGeod::fromDeg(-180.0, clip_top) );
|
||||
|
||||
clipped = tgPolygon::Intersect( subject, clip_row);
|
||||
clipped = tgPolygon::Intersect( subject, clip_row );
|
||||
if ( clipped.TotalNodes() > 0 ) {
|
||||
|
||||
if ( subject.GetPreserve3D() ) {
|
||||
|
|
|
@ -61,7 +61,7 @@ double tgContour::GetArea( void ) const
|
|||
SGVec2d a, b;
|
||||
unsigned int i, j;
|
||||
|
||||
if (node_list.size() ) {
|
||||
if ( node_list.size() ) {
|
||||
j = node_list.size() - 1;
|
||||
for (i=0; i<node_list.size(); i++) {
|
||||
a = SGGeod_ToSGVec2d( node_list[i] );
|
||||
|
@ -77,80 +77,88 @@ double tgContour::GetArea( void ) const
|
|||
// Check that the two supplied points are on the same side of the contour
|
||||
bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) const
|
||||
{
|
||||
//Find equation of line segment joining the points
|
||||
double x1 = firstpt.getLatitudeDeg();
|
||||
double x2 = secondpt.getLatitudeDeg();
|
||||
double y1 = firstpt.getLongitudeDeg();
|
||||
double y2 = secondpt.getLongitudeDeg();
|
||||
//Store differences for later
|
||||
double xdif = x2-x1;
|
||||
double ydif = y2-y1;
|
||||
/*We describe a line parametrically:
|
||||
//Find equation of line segment joining the points
|
||||
double x1 = firstpt.getLatitudeDeg();
|
||||
double x2 = secondpt.getLatitudeDeg();
|
||||
double y1 = firstpt.getLongitudeDeg();
|
||||
double y2 = secondpt.getLongitudeDeg();
|
||||
|
||||
//Store differences for later
|
||||
double xdif = x2-x1;
|
||||
double ydif = y2-y1;
|
||||
|
||||
/*We describe a line parametrically:
|
||||
|
||||
x1 (x2-x1)
|
||||
L = + t
|
||||
y1 (y2-y1)
|
||||
x1 (x2-x1)
|
||||
L = + t
|
||||
y1 (y2-y1)
|
||||
|
||||
with u the parametric coefficient for the second line.
|
||||
Then the line segments intersect if 0 <= t,u <= 1.
|
||||
with u the parametric coefficient for the second line.
|
||||
Then the line segments intersect if 0 <= t,u <= 1.
|
||||
|
||||
To determine t and u we use the approach of Goldman ("Graphics
|
||||
Gems" as described in Stack Overflow question 563198).
|
||||
To determine t and u we use the approach of Goldman ("Graphics
|
||||
Gems" as described in Stack Overflow question 563198).
|
||||
|
||||
if r x s = r_x * s_y - r_y * s_x, then
|
||||
if r x s = r_x * s_y - r_y * s_x, then
|
||||
|
||||
t = (q - p) x s / (r x s)
|
||||
and
|
||||
u = (q - p) x r / (r x s)
|
||||
t = (q - p) x s / (r x s)
|
||||
and
|
||||
u = (q - p) x r / (r x s)
|
||||
|
||||
for line 1 = p + t r, line 2 = q + u s
|
||||
*/
|
||||
//Now cycle over all nodes and count how many times we intersect
|
||||
int intersect_ct = 0;
|
||||
if (node_list.size()) {
|
||||
int j = node_list.size() - 1;
|
||||
for (int i=0;i<node_list.size()-1;i++) {
|
||||
double nx1 = node_list[i].getLatitudeDeg();
|
||||
double ny1 = node_list[i].getLongitudeDeg();
|
||||
double nx2 = node_list[i+1].getLatitudeDeg();
|
||||
double ny2 = node_list[i+1].getLongitudeDeg();
|
||||
double nydif = ny2-ny1;
|
||||
double nxdif = nx2-nx1;
|
||||
double denom = xdif*nydif - ydif*nxdif;
|
||||
if (denom != 0) { //Not parallel
|
||||
double crossx = nx1-x1; double crossy = ny1-y1;
|
||||
double t = (crossx*nydif - crossy*nxdif)/denom;
|
||||
double u = -1*(xdif*crossy - ydif*crossx)/denom;
|
||||
// We consider that an intersection at the edge of the line has
|
||||
// crossed
|
||||
// over, that is, they lie on opposite sides. This way we capture
|
||||
// places where the chopper has clipped a cliff on the tile edge
|
||||
if (t > -0.0001 && t < 1.0001 && u > -0.0001 && u < 1.0001) intersect_ct++;
|
||||
}
|
||||
for line 1 = p + t r, line 2 = q + u s
|
||||
*/
|
||||
|
||||
//Now cycle over all nodes and count how many times we intersect
|
||||
int intersect_ct = 0;
|
||||
if (node_list.size()) {
|
||||
int j = node_list.size() - 1;
|
||||
for (int i=0;i<node_list.size()-1;i++) {
|
||||
double nx1 = node_list[i].getLatitudeDeg();
|
||||
double ny1 = node_list[i].getLongitudeDeg();
|
||||
double nx2 = node_list[i+1].getLatitudeDeg();
|
||||
double ny2 = node_list[i+1].getLongitudeDeg();
|
||||
double nydif = ny2-ny1;
|
||||
double nxdif = nx2-nx1;
|
||||
double denom = xdif*nydif - ydif*nxdif;
|
||||
|
||||
if (denom != 0) { //Not parallel
|
||||
double crossx = nx1-x1; double crossy = ny1-y1;
|
||||
double t = (crossx*nydif - crossy*nxdif)/denom;
|
||||
double u = -1*(xdif*crossy - ydif*crossx)/denom;
|
||||
// We consider that an intersection at the edge of the line has
|
||||
// crossed
|
||||
// over, that is, they lie on opposite sides. This way we capture
|
||||
// places where the chopper has clipped a cliff on the tile edge
|
||||
if (t > -0.0001 && t < 1.0001 && u > -0.0001 && u < 1.0001) intersect_ct++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
bool isinter = (intersect_ct%2 == 0);
|
||||
return isinter;
|
||||
|
||||
bool isinter = (intersect_ct%2 == 0);
|
||||
return isinter;
|
||||
}
|
||||
|
||||
double tgContour::MinDist(const SGGeod& probe) const
|
||||
{
|
||||
SGVec3d probexyz;
|
||||
SGGeodesy::SGGeodToCart(probe,probexyz);
|
||||
double mindist = 100000.0;
|
||||
double dist;
|
||||
if (node_list.size()) {
|
||||
int j = node_list.size() - 1;
|
||||
for (int i=0;i<j;i++) {
|
||||
SGVec3d start,end;
|
||||
SGGeodesy::SGGeodToCart(node_list[i],start);
|
||||
SGGeodesy::SGGeodToCart(node_list[i+1],end);
|
||||
SGLineSegment<double> piece = SGLineSegment<double>(start,end);
|
||||
dist = distSqr(piece,probexyz);
|
||||
if (dist < mindist) mindist = dist;
|
||||
double tgContour::MinDist(const SGGeod& probe) const {
|
||||
SGVec3d probexyz;
|
||||
SGGeodesy::SGGeodToCart( probe,probexyz );
|
||||
double mindist = 100000.0;
|
||||
double dist;
|
||||
|
||||
if ( node_list.size() ) {
|
||||
|
||||
int j = node_list.size() - 1;
|
||||
|
||||
for (int i=0;i<j;i++) {
|
||||
SGVec3d start,end;
|
||||
SGGeodesy::SGGeodToCart( node_list[i],start );
|
||||
SGGeodesy::SGGeodToCart( node_list[i+1],end );
|
||||
SGLineSegment<double> piece = SGLineSegment<double>(start,end);
|
||||
dist = distSqr( piece,probexyz );
|
||||
if (dist < mindist) mindist = dist;
|
||||
}
|
||||
}
|
||||
}
|
||||
return sqrt(mindist);
|
||||
|
||||
return sqrt(mindist);
|
||||
}
|
||||
|
||||
bool tgContour::IsInside( const tgContour& inside, const tgContour& outside )
|
||||
|
|
|
@ -137,7 +137,7 @@ tgPolygon tgPolygon::Diff( const tgPolygon& subject, tgPolygon& clip )
|
|||
}
|
||||
|
||||
//Intersect will keep open paths open
|
||||
tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
|
||||
tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip )
|
||||
{
|
||||
tgPolygon result;
|
||||
UniqueSGGeodSet all_nodes;
|
||||
|
@ -162,14 +162,14 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
|
|||
|
||||
ClipperLib::Clipper c;
|
||||
c.Clear();
|
||||
c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed());
|
||||
c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed() );
|
||||
c.AddPaths(clipper_clip, ClipperLib::PolyType::Clip, true);
|
||||
if(subject.IsClosed()) {
|
||||
c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
|
||||
c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd );
|
||||
result = tgPolygon::FromClipper( clipper_result );
|
||||
}
|
||||
else {
|
||||
c.Execute(ClipperLib::ClipType::Intersection, clipper_tree_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
|
||||
c.Execute(ClipperLib::ClipType::Intersection, clipper_tree_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd );
|
||||
result = tgPolygon::FromClipper( clipper_tree_result );
|
||||
}
|
||||
result = tgPolygon::AddColinearNodes( result, all_nodes );
|
||||
|
|
|
@ -391,7 +391,7 @@ int main( int argc, char **argv ) {
|
|||
if( poDS == NULL )
|
||||
{
|
||||
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening datasource " << datasource );
|
||||
exit( 1 );
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
|
||||
SG_LOG( SG_GENERAL, SG_ALERT, "Processing datasource " << datasource );
|
||||
|
@ -404,7 +404,7 @@ int main( int argc, char **argv ) {
|
|||
if (poLayer == NULL )
|
||||
{
|
||||
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource );
|
||||
exit( 1 );
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
processLayer(poLayer, results );
|
||||
}
|
||||
|
|
Loading…
Reference in a new issue