point2dem: Fix support of nonspherical datums

Apparently we were not handling datums like earth correctly. This
required fixing in Vision Workbench.

src/asp/Core/OrthoRasterizer.h

@@ -42,9 +42,14 @@ namespace cartography {
     typedef std::pair<BBox3, BBox2i> BBoxPair;
     std::vector<BBoxPair > m_point_image_boundaries;
 
+    // This is growing a bbox of points in point projection and Z
+    // values which are altitude.
     struct GrowBBoxAccumulator {
       BBox3 bbox;
-      void operator()( Vector3 const& v ) { if (v != Vector3()) bbox.grow(v); }
+      void operator()( Vector3 const& v ) {
+        if ( !boost::math::isnan(v.z()) )
+          bbox.grow(v);
+      }
     };
 
     struct RemoveSoftInvalid : ReturnFixedType<PixelT> {
@@ -199,7 +204,7 @@ namespace cartography {
           for ( int32 row = 0; row < point_copy.rows()-1; ++row ) {
             PointAcc point_ul = row_acc;
             for ( int32 col = 0; col < point_copy.cols()-1; ++col ) {
-              if ( *point_ul != Vector3() ) {
+              if ( !boost::math::isnan((*point_ul).z()) ) {
 
                 PointAcc point_ur = point_ul; point_ur.next_col();
                 PointAcc point_ll = point_ul; point_ll.next_row();
@@ -208,7 +213,7 @@ namespace cartography {
                 // Verify the point UL and LR are in the rasterable
                 // area. Then just verify that at least one vertice is
                 // placed in the viewable area.
-                if ( *point_lr == Vector3() ||
+                if ( boost::math::isnan((*point_lr).z()) ||
                      !buf_local_3d_bbox.contains(*point_ul) ||
                      !buf_local_3d_bbox.contains(*point_lr) ||
                      !(local_3d_bbox.contains(*point_ul) ||
@@ -237,12 +242,12 @@ namespace cartography {
                 intensities[4] = texture_copy(col+1,row);
 
                 // triangle 1 is: LL, LR, UL
-                if ( *point_ll != Vector3() &&
+                if ( !boost::math::isnan((*point_ll).z()) &&
                      buf_local_3d_bbox.contains(*point_lr) )
                   renderer.DrawPolygon(0, 3);
 
                 // triangle 2 is: UL, LR, UR
-                if ( *point_ur != Vector3() &&
+                if ( !boost::math::isnan((*point_ur).z()) &&
                      buf_local_3d_bbox.contains(*point_ur) )
                   renderer.DrawPolygon(2, 3);
               }

src/asp/Tools/point2dem.cc

@@ -74,6 +74,28 @@ inline point_image_offset( ImageViewBase<ImageT> const& image, Vector3 const& of
   return UnaryPerPixelView<ImageT,PointOffsetFunc>( image.impl(), PointOffsetFunc(offset) );
 }
 
+// Center Longitudes
+class CenterLongitudeFunc : public UnaryReturnSameType {
+  double center;
+public:
+  CenterLongitudeFunc(double c = 0) : center(c) {}
+
+  Vector3 operator()( Vector3 const& v ) const {
+    if ( v[0] < center - 180 )
+      return (*this)(v + Vector3(360,0,0));
+    else if ( v[0] > center + 180 )
+      return (*this)(v - Vector3(360,0,0));
+    return v;
+  }
+};
+
+template <class ImageT>
+UnaryPerPixelView<ImageT, CenterLongitudeFunc>
+inline recenter_longitude( ImageViewBase<ImageT> const& image, double center ) {
+  return UnaryPerPixelView<ImageT, CenterLongitudeFunc>(image.impl(),
+                                                        CenterLongitudeFunc(center));
+}
+
 // Allows FileIO to correctly read/write these pixel types
 namespace vw {
   template<> struct PixelFormatID<Vector3>   { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
@@ -284,15 +306,14 @@ int main( int argc, char *argv[] ) {
       vw_out() << "\t--> Applying offset: " << opt.x_offset
                << " " << opt.y_offset << " " << opt.z_offset << "\n";
       point_image =
-        project_point_image(
-          point_image_offset(
-            xyz_to_lon_lat_radius(point_image, true, avg_location[0] >= 0 ),
-            Vector3(opt.x_offset,opt.y_offset,opt.z_offset)), georef);
+        geodetic_to_point(point_image_offset(recenter_longitude(cartesian_to_geodetic(point_image,georef),
+                                                                avg_location.x() >= 0 ? 0 : 180),
+                                             Vector3(opt.x_offset,
+                                                     opt.y_offset,
+                                                     opt.z_offset)),georef);
     } else {
       point_image =
-        project_point_image(
-          xyz_to_lon_lat_radius(point_image, true, avg_location[0] >= 0 ),
-          georef);
+        geodetic_to_point(recenter_longitude(cartesian_to_geodetic(point_image,georef), avg_location.x() >= 0 ? 0 : 180),georef);
     }
 
     // Rasterize the results to a temporary file on disk so as to speed