Bug fixes and improvements.

src/Prep/E00Lines/main.cxx

@@ -60,6 +60,10 @@ SG_USING_STD(vector);
 // Utility stuff.
 ////////////////////////////////////////////////////////////////////////
 
+
+/**
+ * Inline function to clamp an angle between 0 and 360 degrees.
+ */
 static inline double 
 ANGLE (double a)
 {
@@ -71,7 +75,41 @@ ANGLE (double a)
 }
 
 
+/**
+ * Calculate the intersection of two lines.
+ *
+ * @param p0 First point on the first line.
+ * @param p1 A second point on the first line.
+ * @param p2 First point on the second line.
+ * @param p3 A second point on the second line.
+ * @param intersection A variable to hold the calculated intersection.
+ * @return true if there was an intersection, false if the lines
+ *         are parallel or coincident.
+ */
+static bool
+getIntersection (const Point3D &p0, const Point3D &p1,
+		 const Point3D &p2, const Point3D &p3,
+		 Point3D &intersection)
+{
+  double u_num =
+    ((p3.x()-p2.x())*(p0.y()-p2.y()))-((p3.y()-p2.y())*(p0.x()-p2.x()));
+  double u_den =
+    ((p3.y()-p2.y())*(p1.x()-p0.x()))-((p3.x()-p2.x())*(p1.y()-p0.y()));
 
+  if (u_den == 0) {
+    if (u_num == 0)
+      SG_LOG(SG_GENERAL, SG_ALERT, "Intersection: coincident lines");
+    else
+      SG_LOG(SG_GENERAL, SG_ALERT, "Intersection: parallel lines");
+    return false;
+  } else {
+    double u = u_num/u_den;
+    intersection = Point3D((p0.x()+u*(p1.x()-p0.x())),
+			   (p0.y()+u*(p1.y()-p0.y())),
+			   0);
+    return true;
+  }
+}
 
 
 
@@ -203,17 +241,23 @@ checkAttribute (const E00 &data, int index, const Attribute &att)
 
 /**
  * Create polygons out of points.
+ *
+ * Note that simple geometry doesn't work here, because the scale is
+ * not even -- the points on the x-axis (longitude) become closer and
+ * closer as the y-axis (latitude) approaches the poles, meeting in
+ * a single point at y=90 and y=-90.  As a result, this function
+ * uses the WGS80 functions, rather than simple Pythagorean stuff.
  */
 static void
 processPoints (const E00 &data, const Rectangle &bounds,
 	       AreaType areaType, const string &workDir, int width)
 {
-  FGPolygon shape;
   double x, y, az;
 
   int nPoints = data.nPoints();
   cout << "Processing " << nPoints << " points" << endl;
   for (int i = 1; i <= nPoints; i++) {
+    FGPolygon shape;
     const E00::LAB &lab = data.getLAB(i);
     double lon = lab.coord.x;
     double lat = lab.coord.y;
@@ -243,17 +287,22 @@ processPoints (const E00 &data, const Rectangle &bounds,
 
 /**
  * Create polygons out of all loose line segments.
+ *
+ * Note that simple geometry doesn't work here, because the scale is
+ * not even -- the points on the x-axis (longitude) become closer and
+ * closer as the y-axis (latitude) approaches the poles, meeting in
+ * a single point at y=90 and y=-90.  As a result, this function
+ * uses the WGS80 functions, rather than simple Pythagorean stuff.
  */
 static void
 processLines (const E00 &data, const Rectangle &bounds,
 	      AreaType areaType, const string &workDir, int width,
 	      const vector<Attribute> &aat_list)
 {
-  FGPolygon shape;
-
   int nLines = data.nLines();
   cout << "Processing " << nLines << " lines." << endl;
   for (int i = 1; i <= nLines; i++) {
+    FGPolygon shape;
     const E00::ARC &arc = data.getARC(i);
     Rectangle arcBounds = makeBounds(arc);
     if (!bounds.isOverlapping(arcBounds)) {
@@ -279,51 +328,88 @@ processLines (const E00 &data, const Rectangle &bounds,
       }
     }
 
-				// Make the line into a polygon.
+				// Put the rectangles for the segments
-    cout << "Line has " << arc.numberOfCoordinates << " coordinates" << endl;
+				// into a list
-    for (int j = 0; j < arc.numberOfCoordinates - 1; j++) {
+    vector<FGPolygon> segment_list;
+
+    int j;
+    for (j = 0; j < arc.numberOfCoordinates - 1; j++) {
       double lon1 = arc.coordinates[j].x;
       double lat1 = arc.coordinates[j].y;
       double lon2 = arc.coordinates[j+1].x;
       double lat2 = arc.coordinates[j+1].y;
       double angle1, angle2, dist;
       geo_inverse_wgs_84(0, lat1, lon1, lat2, lon2, &angle1, &angle2, &dist);
-      cout << "angle1 = " << angle1 << endl;
-      cout << "angle2 = " << angle2 << endl;
-      cout << "dist = " << dist << endl;
       shape.erase();
 
       double x, y, az;
       
+				// Wind each rectangle counterclockwise
+
 				// Corner 1
       geo_direct_wgs_84(0, lat1, lon1, ANGLE(angle1+90), width/2, &y, &x, &az);
-      cout << x << '\t' << y << endl;
       shape.add_node(0, Point3D(x, y, 0));
 
 				// Corner 2
-      geo_direct_wgs_84(0, lat1, lon1, ANGLE(angle1-90), width/2, &y, &x, &az);
+      geo_direct_wgs_84(0, lat2, lon2, ANGLE(angle1+90), width/2, &y, &x, &az);
-      cout << x << '\t' << y << endl;
       shape.add_node(0, Point3D(x, y, 0));
 
 				// Corner 3
-      geo_direct_wgs_84(0, lat2, lon2, ANGLE(angle2+90), width/2, &y, &x, &az);
+      geo_direct_wgs_84(0, lat2, lon2, ANGLE(angle1-90), width/2, &y, &x, &az);
-      cout << x << '\t' << y << endl;
       shape.add_node(0, Point3D(x, y, 0));
 
 				// Corner 4
-      geo_direct_wgs_84(0, lat2, lon2, ANGLE(angle2-90), width/2, &y, &x, &az);
+      geo_direct_wgs_84(0, lat1, lon1, ANGLE(angle1-90), width/2, &y, &x, &az);
-      cout << x << '\t' << y << endl;
       shape.add_node(0, Point3D(x, y, 0));
 
-				// Corner 1, again
+				// Save this rectangle
-      geo_direct_wgs_84(0, lat1, lon1, ANGLE(angle1+90), width/2, &y, &x, &az);
+      segment_list.push_back(shape);
-      cout << x << '\t' << y << endl;
-      shape.add_node(0, Point3D(x, y, 0));
-
-				// Split into tiles
-      split_polygon(workDir, areaType, shape);
     }
+
+    // Build one big polygon out of all the rectangles by intersecting
+    // the lines running through the bottom and top sides
+
+    shape.erase();
+
+				// Connect the bottom part.
+    int nSegments = segment_list.size();
+    Point3D intersection;
+    shape.add_node(0, segment_list[0].get_pt(0, 0));
+    for (j = 0; j < nSegments - 1; j++) {
+      if (getIntersection(segment_list[j].get_pt(0, 0),
+			  segment_list[j].get_pt(0, 1),
+			  segment_list[j+1].get_pt(0, 0),
+			  segment_list[j+1].get_pt(0, 1),
+			  intersection))
+	shape.add_node(0, intersection);
+      else
+	shape.add_node(0, segment_list[j].get_pt(0, 1));
+    }
+    shape.add_node(0, segment_list[nSegments-1].get_pt(0, 1));
+
+				// Connect the top part
+    shape.add_node(0, segment_list[nSegments-1].get_pt(0, 2));
+    for (j = nSegments - 1; j > 0; j--) {
+      if (getIntersection(segment_list[j].get_pt(0, 2),
+			  segment_list[j].get_pt(0, 3),
+			  segment_list[j-1].get_pt(0, 2),
+			  segment_list[j-1].get_pt(0, 3),
+			  intersection))
+	shape.add_node(0, intersection);
+      else
+	shape.add_node(0, segment_list[j].get_pt(0, 3));
+    }
+    shape.add_node(0, segment_list[0].get_pt(0, 3));
+    
+    				// Split into tiles
+    cout << "Splitting polygon..." << endl;
+    cout << "  Total size: " << shape.total_size() << endl;
+    cout << "  Minimum angle: "
+	 << (shape.minangle_contour(0) * SGD_RADIANS_TO_DEGREES) << endl;
+    
+    split_polygon(workDir, areaType, shape);
   }
+  cout << "Done lines" << endl;
 }
 
 
@@ -335,13 +421,11 @@ processPolygons (const E00 &data, const Rectangle &bounds,
 		 AreaType areaType, const string &workDir,
 		 const vector<Attribute> pat_list)
 {
-  FGPolygon shape;
-
   int nPolygons = data.nPolygons();
   cout << "Processing " << nPolygons << " polygons" << endl;
 
   for (int i = 2; i <= nPolygons; i++) {
-
+    FGPolygon shape;
 				// Test whether the polygon matches
 				// at least one of the attributes
 				// provided.
@@ -490,9 +574,9 @@ main (int argc, const char **argv)
     }
 
     else if (arg.find("--lines=") == 0) {
-      if (arg.substr(9) == "yes")
+      if (arg.substr(8) == "yes")
 	useLines = true;
-      else if (arg.substr(9) == "no")
+      else if (arg.substr(8) == "no")
 	useLines = false;
       else {
 	cerr << "--lines option needs 'yes' or 'no'" << endl;
@@ -502,9 +586,9 @@ main (int argc, const char **argv)
     }
 
     else if (arg.find("--polygons=") == 0) {
-      if (arg.substr(9) == "yes")
+      if (arg.substr(11) == "yes")
 	usePolygons = true;
-      else if (arg.substr(9) == "no")
+      else if (arg.substr(11) == "no")
 	usePolygons = false;
       else {
 	cerr << "--polygons option needs 'yes' or 'no'" << endl;
@@ -676,4 +760,4 @@ main (int argc, const char **argv)
   return 0;
 }
 
-
+// end of main.cxx

src/Prep/TGVPF/tgvpf.cxx

@@ -55,6 +55,14 @@ SG_USING_STD(vector);
 #endif
 
 
+
+////////////////////////////////////////////////////////////////////////
+// Program-wide variables.
+////////////////////////////////////////////////////////////////////////
+
+static const char * progname;
+
+
 
 ////////////////////////////////////////////////////////////////////////
 // Utility stuff.
@@ -326,18 +334,20 @@ makePolygon (const VpfPolygon &polygon)
  * Print the command-line usage and exit.
  */
 static void
-usage (const char * prog)
+usage ()
 {
-  cerr << "Usage: " << prog << " [opts] <db> <library> <coverage> <feature>"
+  cerr << "Usage: "
+       << progname
+       << " [opts] <db> <library> <coverage> <feature>"
        << endl;
   cerr << "Options:" << endl;
+  cerr << "--chunk=<chunk> (default: none)" << endl;
   cerr << "--min-lon=<longitude> (default: -180.0)" << endl;
   cerr << "--min-lat=<latitude> (default: -90.0)" << endl;
   cerr << "--max-lon=<longitude> (default: 180.0)" << endl;
   cerr << "--max-lat=<latitude> (default: 90.0)" << endl;
   cerr << "--area=<area_type> (default: Default)" << endl;
-  cerr << "--point-width=<meters> (default: 500)" << endl;
+  cerr << "--width=<meters> (default: 50 line, 500 point)" << endl;
-  cerr << "--line-width=<meters> (default: 50)" << endl;
   cerr << "--work-dir=<dir> (default: .)" << endl;
   cerr << "--att=<item>:<value> (may be repeated)" << endl;
   cerr << "--att=!<item>:<value> (may be repeated)" << endl;
@@ -345,11 +355,54 @@ usage (const char * prog)
 }
 
 
+/**
+ * Parse a 10x10 degree chunk name.
+ */
+static VpfRectangle
+parseChunk (string chunk)
+{
+  VpfRectangle bounds;
+  int x_factor;
+  int y_factor;
+
+  if (chunk.size() != 7) {
+    cerr << "Bad length for chunk specifier " << chunk << endl;
+    usage();
+  }
+
+  if (chunk[0] == 'w')
+    x_factor = -1;
+  else if (chunk[0] == 'e')
+    x_factor = 1;
+  else {
+    cerr << "Chunk specifier must begin with 'e' or 'w'" << endl;
+    usage();
+  }
+
+  if (chunk[4] == 's')
+    y_factor = -1;
+  else if (chunk[4] == 'n')
+    y_factor = 1;
+  else {
+    cerr << "Second part of chunk specifier must begin with 's' or 'n'"
+	 << endl;
+    usage();
+  }
+
+  bounds.minX = atoi(chunk.substr(1,3).c_str()) * x_factor;
+  bounds.minY = atoi(chunk.substr(5).c_str()) * y_factor;
+  bounds.maxX = bounds.minX + 10;
+  bounds.maxY = bounds.minY + 10;
+
+  return bounds;
+}
+
+
 /**
  * Parse an attribute value specification from the command line.
  */
 static const Attribute
-parseAttribute (const char * prog, string arg)
+parseAttribute (string arg)
 {
   Attribute att;
 
@@ -363,7 +416,7 @@ parseAttribute (const char * prog, string arg)
   int pos = arg.find(':');
   if (pos == -1) {
     cerr << "Bad attribute specification: " << arg << endl;
-    usage(prog);
+    usage();
   }
 
   att.name = arg.substr(0, pos);
@@ -379,6 +432,10 @@ parseAttribute (const char * prog, string arg)
 int
 main (int argc, const char **argv)
 {
+				// Store the program name for future
+				// reference.
+  progname = argv[0];
+
   vector<Attribute> attributes;
   VpfRectangle bounds;
 
@@ -404,7 +461,12 @@ main (int argc, const char **argv)
   while (argPos < argc) {
     string arg = argv[argPos];
 
-    if (arg.find("--min-lon=") == 0) {
+    if (arg.find("--chunk=") == 0) {
+      bounds = parseChunk(arg.substr(8));
+      argPos++;
+    }
+
+    else if (arg.find("--min-lon=") == 0) {
       bounds.minX = strtod(arg.substr(10).c_str(), 0);
       argPos++;
     }
@@ -445,7 +507,7 @@ main (int argc, const char **argv)
     }
 
     else if (arg.find("--att=") == 0) {
-      attributes.push_back(parseAttribute(argv[0], arg.substr(6)));
+      attributes.push_back(parseAttribute(arg.substr(6)));
       argPos++;
     }
 
@@ -456,7 +518,7 @@ main (int argc, const char **argv)
 
     else if (arg.find("-") == 0) {
       cerr << "Unrecognized option: " << arg << endl;
-      usage(argv[0]);
+      usage();
     }
 
     else {
@@ -464,13 +526,40 @@ main (int argc, const char **argv)
     }
   }
 
+  //
+  // Sanity check on bounds.
+  //
+  if (bounds.minX < -180) {
+    cerr << "Minimum longitude out of range (-180:180): "
+	 << bounds.minX << endl;
+    usage();
+  } else if (bounds.maxX > 180) {
+    cerr << "Maximum longitude out of range (-180:180): "
+	 << bounds.maxX << endl;
+    usage();
+  } else if (bounds.minY < -90) {
+    cerr << "Minimum latitude out of range (-90:90): "
+	 << bounds.minY << endl;
+    usage();
+  } else if (bounds.maxY > 90) {
+    cerr << "Maximum latitude out of range (-90:90): "
+	 << bounds.maxY << endl;
+    usage();
+  } else if (bounds.minX >= bounds.maxX) {
+    cerr << "Minimum longitude less than maximum longitude" << endl;
+    usage();
+  } else if (bounds.minY >= bounds.maxY) {
+    cerr << "Minimum latitude less than maximum latitude" << endl;
+    usage();
+  }
+
 
   //
   // Process command-line arguments.
   //
 
   if (argPos != (argc - 4))
-    usage(argv[0]);
+    usage();
 
   const char * database_name = argv[argPos++];
   const char * library_name = argv[argPos++];