fly/terragear
1
0
Fork 0
Code Activity

Fixed height rectification along boundaries and proper clipping of

Browse source 
open contours.
This commit is contained in:
James.Hester 2018-12-17 23:50:28 +11:00
parent 34b3bd89ff
commit a7af59b61f
7 changed files with 141 additions and 33 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/BuildTiles/Main/tgconstruct_poly.cxx
View file

@ -70,7 +70,7 @@ int TGConstruct::LoadLandclassPolys( void ) {
string lastext = p.extension();
if ((lext == "arr") || (lext == "arr.gz") || (lext == "btg.gz") ||
(lext == "fit") || (lext == "fit.gz") || (lext == "ind") ||
(lastext == "cliffs"))
(lastext == "cliffs") || (lext == "arr.rectified.gz") || (lext == "arr.new.gz"))
{
// skipped!
} else {

11
src/Lib/Array/CMakeLists.txt
View file

@ -15,6 +15,8 @@ install( FILES array.hxx DESTINATION include/tg )
add_executable(test_array testarray.cxx)
add_executable(rectify_height rectify_height.cxx)
target_link_libraries(test_array
Array
terragear
@ -22,7 +24,14 @@ target_link_libraries(test_array
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
install(TARGETS test_array RUNTIME DESTINATION bin)
target_link_libraries(rectify_height
Array
terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
install(TARGETS test_array rectify_height RUNTIME DESTINATION bin)
if (MSVC)
set_target_properties( test_array PROPERTIES DEBUG_POSTFIX d )
endif ()

107
src/Lib/Array/array.cxx
View file

@ -57,16 +57,27 @@ TGArray::TGArray( const string &file ):
// open an Array file (and fitted file if it exists)
// Also open cliffs file if it exists
// Also open cliffs file if it exists. By default a
// rectified file is searched for, if it doesn't exist
// we load the unrectified version.
bool TGArray::open( const string& file_base ) {
// open array data file
string array_name = file_base + ".arr.gz";
string array_name = file_base + ".arr.rectified.gz";
rectified = true;
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;
}
}
SG_LOG(SG_GENERAL,SG_DEBUG,"Loaded height array " << array_name);
// open fitted data file
string fitted_name = file_base + ".fit.gz";
fitted_in = new sg_gzifstream( fitted_name );
@ -206,10 +217,7 @@ TGArray::parse( SGBucket& b ) {
}
}
// Fix heights near cliffs
if(cliffs_list.size()>0) rectify_heights();
return true;
return true;
}
void TGArray::parse_bin()
@ -240,6 +248,40 @@ void TGArray::parse_bin()
sgReadShort(array_in, cols * rows, in_data);
}
// Write out an array. If rectified is true, the heights have been adjusted
// for discontinuities.
void TGArray::write_bin(const string root_dir, bool rectified, SGBucket& b) {
// generate output file name
string base = b.gen_base_path();
string path = root_dir + "/" + base;
string extension = ".arr.new.gz";
if (rectified) extension = ".arr.rectified.gz";
SGPath sgp( path );
sgp.append( "dummy" );
sgp.create_dir( 0755 );
string array_file = path + "/" + b.gen_index_str() + extension;
SG_LOG(SG_GENERAL, SG_DEBUG, "array_file = " << array_file );
// write the file
gzFile fp;
if ( (fp = gzopen( array_file.c_str(), "wb9" )) == NULL ) {
SG_LOG(SG_GENERAL, SG_ALERT, "ERROR: cannot open " << array_file << " for writing!" );
return;
}
int32_t header = 0x54474152; //'TGAR'
sgWriteLong(fp,header);
sgWriteInt(fp,originx);
sgWriteInt(fp,originy);
sgWriteInt(fp,cols);
sgWriteInt(fp,col_step);
sgWriteInt(fp,rows);
sgWriteInt(fp,row_step);
sgWriteShort(fp, rows*cols, in_data);
gzclose(fp);
}
// write an Array file
bool TGArray::write( const string root_dir, SGBucket& b ) {
// generate output file name
@ -273,7 +315,6 @@ bool TGArray::write( const string root_dir, SGBucket& b ) {
return true;
}
// do our best to remove voids by picking data from the nearest neighbor.
void TGArray::remove_voids( ) {
// need two passes to ensure that all voids are removed (unless entire
@ -384,7 +425,7 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
}
}
std::vector<int> TGArray::collect_bad_points() {
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
@ -392,7 +433,7 @@ std::vector<int> TGArray::collect_bad_points() {
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)) {
if(is_near_cliff(lon,lat,bad_zone)) {
bad_points.push_back(horiz+vert*cols);
}
}
@ -412,10 +453,10 @@ bool TGArray::is_bad_point(const int xgrid, const int ygrid, const std::vector<i
//This may collide with other threads, but as they will both be writing
//the correct height, this is harmless.
void TGArray::rectify_heights() {
void TGArray::rectify_heights(const double bad_zone) {
double new_ht;
std::vector<int> rectified,bad_points;
bad_points = collect_bad_points();
bad_points = collect_bad_points(bad_zone);
while(1) {
for (auto pt : bad_points) {
int ygrid = pt/cols;
@ -433,7 +474,10 @@ void TGArray::rectify_heights() {
}
rectified.clear();
} else {
break;
if(bad_points.size() > 0) {
std::cout << "Failed to rectify " << bad_points.size() << " points" << std::endl;
}
break; // Cant do any more
}
}
}
@ -465,8 +509,10 @@ double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vecto
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
@ -479,16 +525,36 @@ double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vecto
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)) break;
}
if (pt == pt_cnt) return -9999; // not found
!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 are in a bay, 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);
double height = horiz + (vert - corner);
if (corner < -9000) corner = 0;
if (horiz < -9000) horiz = 0;
if (vert < -9000) vert = 0;
height = horiz + (vert - corner);
}
std::cout << xgrid << "," << ygrid << ": was " << original_height << " , now " << height << std::endl;
return height;
}
@ -768,13 +834,12 @@ bool TGArray::check_points (const double lon1, const double lat1, const double l
//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 {
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);
double mindist = 30; // 1 arcsec
for (int i=0;i<cliffs_list.size();i++) {
double dist = cliffs_list[i].MinDist(pt1);
if (dist < mindist) return true;
if (dist < bad_zone) return true;
}
return false;
}

13
src/Lib/Array/array.hxx
View file

@ -47,6 +47,8 @@ private:
// number of columns and rows
int cols, rows;
// Whether or not the input data have been rectified
bool rectified;
// Distance between column and row data points (in arc seconds)
double col_step, row_step;
@ -62,11 +64,10 @@ private:
void parse_bin();
// Routines for height rectification
void rectify_heights();
std::vector<int> collect_bad_points();
std::vector<int> collect_bad_points(const double bad_zone);
bool is_bad_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const;
double rectify_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const;
bool is_near_cliff(const double lon1,const double lon2) const;
bool is_near_cliff(const double lon1,const double lon2, const double bad_zone) const;
public:
// Constructor
@ -94,10 +95,16 @@ public:
// write an Array file
bool write( const std::string root_dir, SGBucket& b );
// write an Array file in binary format. If ht_rect is true,
// the file will have extension 'arr.rectified.gz'
void write_bin(const std::string root_dir, bool ht_rect, SGBucket& b);
// do our best to remove voids by picking data from the nearest
// neighbor.
void remove_voids();
void rectify_heights(const double bad_zone);
// Return the elevation of the closest non-void grid point to lon, lat
double closest_nonvoid_elev( double lon, double lat ) const;

15
src/Lib/terragear/tg_chopper.cxx
View file

@ -21,8 +21,18 @@ 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);
result = tgPolygon::Intersect( subject, base);
// Debug: See if numbers of nodes have changed
if ( result.Contours() > 0 ) {
/*if (result.ContourSize(0) != subject.ContourSize(0)) {
tgRectangle rr = subject.GetBoundingBox();
SG_LOG(SG_GENERAL,SG_INFO,"---- Bucket ID " << b.gen_index() << " ------- " );
SG_LOG(SG_GENERAL,SG_INFO,"A contour has changed size");
SG_LOG(SG_GENERAL,SG_INFO,"Bounding box " << rr.getMin() << " , " << rr.getMax() );
SG_LOG(SG_GENERAL,SG_INFO,"Old size " << subject.ContourSize(0) << " New size " << result.ContourSize(0));
SG_LOG(SG_GENERAL,SG_INFO,"Old: First node " << subject.GetNode(0,0) << " New: First node " << result.GetNode(0,0));
SG_LOG(SG_GENERAL,SG_INFO,"Clipping rectangle: " << base.GetContour(0));
}*/
if ( subject.GetPreserve3D() ) {
result.InheritElevations( subject );
result.SetPreserve3D( true );
@ -33,6 +43,9 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
result.SetTexMethod( TG_TEX_BY_GEODE, b.get_center_lat() );
}
result.SetFlag(type);
if (!subject.IsClosed()) {
result.SetOpen();
}
lock.lock();
bp_map[b.gen_index()].push_back( result );

20
src/Lib/terragear/tg_contour.cxx
View file

@ -83,8 +83,8 @@ bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) cons
double y1 = firstpt.getLongitudeDeg();
double y2 = secondpt.getLongitudeDeg();
//Store differences for later
double xdif = x1-x2;
double ydif = y1-y2;
double xdif = x2-x1;
double ydif = y2-y1;
/*We describe a line parametrically:
x1 (x2-x1)
@ -94,6 +94,16 @@ bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) cons
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).
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)
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;
@ -104,11 +114,11 @@ Then the line segments intersect if 0 <= t,u <= 1.
double ny1 = node_list[i].getLongitudeDeg();
double nx2 = node_list[i+1].getLatitudeDeg();
double ny2 = node_list[i+1].getLongitudeDeg();
double nydif = ny1-ny2;
double nxdif = nx1-nx2;
double nydif = ny2-ny1;
double nxdif = nx2-nx1;
double denom = xdif*nydif - ydif*nxdif;
if (denom != 0) { //Not parallel
double crossx = x1-nx1; double crossy = y1-ny1;
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 have

6
src/Lib/terragear/tg_polygon_clip.cxx
View file

@ -178,7 +178,11 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
result.SetTexParams( subject.GetTexParams() );
result.SetId( subject.GetId() );
result.SetPreserve3D( subject.GetPreserve3D() );
if(subject.IsClosed()) {
result.SetClosed();
} else {
result.SetOpen();
}
return result;
}