From a7af59b61f47fa6d8c57a64c1e244a735b64e4c3 Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Mon, 17 Dec 2018 23:50:28 +1100
Subject: [PATCH] Fixed height rectification along boundaries and proper
clipping of open contours.
---
src/BuildTiles/Main/tgconstruct_poly.cxx | 2 +-
src/Lib/Array/CMakeLists.txt | 11 ++-
src/Lib/Array/array.cxx | 107 ++++++++++++++++++-----
src/Lib/Array/array.hxx | 13 ++-
src/Lib/terragear/tg_chopper.cxx | 15 +++-
src/Lib/terragear/tg_contour.cxx | 20 +++--
src/Lib/terragear/tg_polygon_clip.cxx | 6 +-
7 files changed, 141 insertions(+), 33 deletions(-)
diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx
index 5ebec27a..bfb6e76d 100644
--- a/src/BuildTiles/Main/tgconstruct_poly.cxx
+++ b/src/BuildTiles/Main/tgconstruct_poly.cxx
@@ -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 {
diff --git a/src/Lib/Array/CMakeLists.txt b/src/Lib/Array/CMakeLists.txt
index 989241aa..2d4f870e 100644
--- a/src/Lib/Array/CMakeLists.txt
+++ b/src/Lib/Array/CMakeLists.txt
@@ -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 ()
diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx
index d5b1f1a5..9c159d96 100644
--- a/src/Lib/Array/array.cxx
+++ b/src/Lib/Array/array.cxx
@@ -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;
}
diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx
index cfa2b1bb..9de84269 100644
--- a/src/Lib/Array/array.hxx
+++ b/src/Lib/Array/array.hxx
@@ -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;
diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx
index f9a52e9f..25cc71e4 100644
--- a/src/Lib/terragear/tg_chopper.cxx
+++ b/src/Lib/terragear/tg_chopper.cxx
@@ -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 );
diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx
index b4aea891..160f92f3 100644
--- a/src/Lib/terragear/tg_contour.cxx
+++ b/src/Lib/terragear/tg_contour.cxx
@@ -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
diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx
index ec1c8436..c73b7436 100644
--- a/src/Lib/terragear/tg_polygon_clip.cxx
+++ b/src/Lib/terragear/tg_polygon_clip.cxx
@@ -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;
}