point2dem: Add optio --auto-proj-center

NEWS.rst

@@ -40,6 +40,10 @@ parallel_stereo (:numref:`parallel_stereo`):
   * For the ``asp_sgm`` and ``asp_mgm`` algorithms allow ``cost-mode`` to
     have the value 3 or 4 only, as other values produce bad results. 
 
+point2dem (:numref:`point2dem`):
+  * Added the option ``--auto-proj-center``, to automatically compute the 
+    projection center for stereographic and other projections. 
+
 stereo_gui (:numref:`stereo_gui`):
   * Can show scattered data with a colorbar and axes 
     (:numref:`scattered_points_colorbar`).

docs/next_steps.rst

@@ -315,17 +315,16 @@ The ``parallel_stereo`` command can also take multiple input images,
 performing multi-view stereo (:numref:`multiview`), though this
 approach is rather discouraged as better results can be obtained with
 bundle adjustment followed by pairwise stereo and merging of DEMs with
-``dem_mosaic``.
+``dem_mosaic`` (:numref:`dem_mosaic`).
 
 Running the GUI frontend
 ~~~~~~~~~~~~~~~~~~~~~~~~
 
-The ``stereo_gui`` program is a GUI frontend to
-``parallel_stereo``. It is invoked with the same options as
-``parallel_stereo`` (except for the more specialized ones such as
-``--job-size-h``, etc.). It displays the input images, and makes it
-possible to zoom in and select smaller regions to run stereo on. The
-GUI is described in :numref:`stereo_gui`.
+The ``stereo_gui`` program (:numref:`stereo_gui`) is a GUI frontend to
+``parallel_stereo``. It is invoked with the same options as ``parallel_stereo``
+(except for the more specialized ones such as ``--job-size-h``, etc.). It
+displays the input images, and makes it possible to zoom in and select smaller
+regions to run stereo on. 
 
 .. _cmdline:
 
@@ -424,16 +423,15 @@ additional 190 meter vertical offset (such as the dataset
 ``molaMarsPlanetaryRadius0001.cub`` shipped with ISIS data), which can
 be taken care of with ``image_calc`` (:numref:`image_calc`).
 
-Alternatively, a low-resolution smooth DEM can be obtained by running
-ASP itself as described in previous sections. In such a run, subpixel
-mode may be set to parabola (``subpixel-mode 1``) for speed. To make it
-sufficiently coarse and smooth, the resolution can be set to about 40
-times coarser than either the default ``point2dem`` resolution or the
-resolution of the input images. If the resulting DEM turns out to be
-noisy or have holes, one could change in ``point2dem`` the search radius
-factor, use hole-filling, invoke more aggressive outlier removal, and
-erode pixels at the boundary (those tend to be less reliable).
-Alternatively, holes can be filled with ``dem_mosaic``.
+Alternatively, a low-resolution smooth DEM can be obtained by running ASP itself
+as described in previous sections. In such a run, subpixel mode may be set to
+parabola (``subpixel-mode 1``) for speed. To make it sufficiently coarse and
+smooth, the resolution can be set to about 40 times coarser than either the
+default ``point2dem`` (:numref:`point2dem`) resolution or the resolution of the
+input images. If the resulting DEM turns out to be noisy or have holes, one
+could change in ``point2dem`` the search radius factor, use hole-filling, invoke
+more aggressive outlier removal, and erode pixels at the boundary (those tend to
+be less reliable). Alternatively, holes can be filled with ``dem_mosaic``.
 
 .. _conv_to_ellipsoid:
 
@@ -542,13 +540,15 @@ We used ``--search-radius-factor 5`` to expand the DEM a
 bit, to counteract future erosion at image boundary in stereo due to
 the correlation kernel size. This is optional.
 
-If this terrain is close to the poles, say within 25 degrees of
-latitude, it is advised to use a stereographic projection, centered
-either at the nearest pole, or close to the center of the current DEM.
-Its center's longitude and latitude can be found with
-``gdalinfo -stats``, which can then be passed to ``point2dem`` such as::
+If this terrain is close to the poles, say within 25 degrees of latitude, it is
+advised to use a stereographic projection, centered either at the nearest pole,
+or close to the center of the current DEM. 
 
-     point2dem --stereographic --proj-lon <lon_ctr> --proj-lat <lat_ctr> ...
+The center of this projection can be auto-computed by ``point2dem`` if invoked 
+with the option ``--auto-proj-center``. Or, it can be found in advance
+with ``gdalinfo -stats``, which can then be passed to ``point2dem`` such as::
+
+     point2dem --stereographic --proj-lon <lon_ctr> --proj-lat <lat_ctr> 
 
 By calling ``gdalinfo -proj4``, the PROJ.4 string of the obtained DEM
 can be found, which can be used in mapprojection later, and with the
@@ -919,10 +919,10 @@ tune the parameters to improve the results.
    significantly different absolute ranges of disparity.
 
 Whenever ``parallel_stereo``, ``point2dem``, and other executables are run, they
-create log files in given tool's results directory, containing a copy of
-the configuration file, the command that was run, your system settings,
-and tool's console output. This will help track what was performed so
-that others in the future can recreate your work.
+create log files in given tool's results directory, containing a copy of the
+configuration file, the command that was run, your system settings, and tool's
+console output. This will help track what was performed so that others in the
+future can recreate your work.
 
 Another handy debugging tool is the ``disparitydebug`` program
 (:numref:`disparitydebug`), which allows you to generate viewable
@@ -1167,6 +1167,10 @@ Elevation Model (DEM) from the point cloud file.
 The resulting TIFF file is mapprojected and will contain georeference
 information stored as a GeoTIFF header.
 
+The default projection is geographic (latitude and longitude), which 
+is not great at the poles. See :numref:`point2dem` for how to change
+the projection and auto-compute its center, if desired.
+
 The tool will infer the datum and projection from the input images, if
 present. You can explicitly specify a coordinate system (e.g., mercator,
 sinusoidal) and a reference spheroid (i.e., calculated for the Moon,

docs/tools/bundle_adjust.rst

@@ -77,12 +77,15 @@ Here we assumed that the cameras point towards some planet's surface and
 used the ``nadirpinhole`` session. If this assumption is not true, one
 should use the ``pinhole`` session or the ``--no-datum`` option.
 
-Using uniformly distributed interest points
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Use cases
+~~~~~~~~~
+
+Uniformly distributed interest points
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 When different parts of the image have different properties, such as rock vs snow,
 additional work may be needed to ensure interest points are created somewhat
-uniformly. Then, one can use the option ``--matches-per-tile``::
+uniformly. For that, use the option ``--matches-per-tile``::
 
     bundle_adjust image1.tif image2.tif       \
         image1.tsai image2.tsai               \
@@ -94,7 +97,7 @@ uniformly. Then, one can use the option ``--matches-per-tile``::
         -o run_ba/run 
 
 For very large images, the number of interest points and matches per tile (whose
-size is 1024 pixels on the side) should decrease from the above. 
+size is 1024 pixels on the side) should be decreased from the above. 
 
 Controlling where interest points are placed
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -162,6 +165,9 @@ can then be passed directly to stereo::
      stereo file1.JPG file2.JPG run_ba/run-file1.tsai \
        run_ba/run-file2.tsai run_stereo/run
 
+When cameras are of CSM type (:numref:`csm`), self-contained optimized
+cameras will be written to disk (:numref:`csm_state`).
+
 The ``bundle_adjust`` program can read camera adjustments from a
 previous run, via ``--input-adjustments-prefix string``. It can also
 apply to the input cameras a transform as output by ``pc_align``, via

docs/tools/point2dem.rst

@@ -28,15 +28,28 @@ The obtained DEMs can be colorized or hillshaded
 Examples
 ~~~~~~~~
 
-Create a DEM only::
+When creating an DEM it is very important to pick a projection 
+that is appropriate for the region of interest. The default
+projection is geographic (longitude and latitude), which is not 
+good for regions close to the poles. In such cases, it is best
+to pick a stereographic projection.
 
-    point2dem run/run-PC.tif
+Auto-guess projection center and datum
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+::
+
+    point2dem --stereographic --auto-proj-center run/run-PC.tif
 
 This creates ``run/run-DEM.tif``, which is a GeoTIFF file, with each
 32-bit floating point pixel value being the height above the datum
 (ellipsoid). The datum is saved in the geoheader and can be seen with
 ``gdalinfo`` (:numref:`gdal_tools`).
 
+In this case the stereographic projection was used, and its center was
+auto-guessed as the median longitude and latitude for the 
+points in the cloud. The grid size was also auto-guessed.
+
 ASP normally auto-guesses the datum, otherwise the option ``-r`` can
 be used. If desired to change the output no-data value (which can also
 be inspected with ``gdalinfo``), use the options ``--nodata-value``.
@@ -45,15 +58,20 @@ If desired to change the range of longitudes from [0, 360] to [-180,
 180], or vice-versa, post-process obtained DEM with ``image_calc``
 (:numref:`image_calc`).
 
-Create a DEM, orthoimage, and intersection error image::
+Orthoimage and error image
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+::
 
     point2dem run/run-PC.tif -r moon --errorimage \
         --orthoimage run/run-L.tif
 
-This produced the DEM, and also takes the left input image and
-orthographically projects it onto the DEM. The resulting
-``run/run-DRG.tif`` file will be saved as a GeoTIFF image with the
-same geoheader as the DEM.
+This produced the DEM, in the default geographic projection (longitude and
+latitude, which sometimes is problematic).
+
+Then, the left aligned image was used to create an orthoimage, by
+orthographically projecting it onto the DEM. The resulting ``run/run-DRG.tif``
+file will be saved as a GeoTIFF image with the same geoheader as the DEM.
 
 In addition, the file ``run/run-IntersectionErr.tif`` is created,
 based on the 4th band of the ``PC.tif`` file, having the gridded
@@ -66,7 +84,10 @@ Here we have explicitly specified the spheroid (``-r moon``), rather
 than have it inferred automatically. The Moon spheroid will have a
 radius of 1737.4 km.
 
-Example with setting the grid size::
+Custom grid size with geographic projection
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+::
 
     point2dem --tr 0.0001 run/run-PC.tif
 
@@ -79,12 +100,18 @@ automatically, so not specifying ``--tr`` at all, or otherwise use a
 multiple of the automatically determined grid size
 (:numref:`post-spacing`).
 
-Example with stereographic projection (for data close to poles)::
+Polar stereographic projection
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+::
 
      point2dem --stereographic --proj-lon 0 --proj-lat -90 \
        run/run-PC.tif
 
-Example with multiple input clouds::
+Multiple clouds
+^^^^^^^^^^^^^^^
+
+::
 
      point2dem in1.las in2.csv run/run-PC.tif -o combined \
        --dem-spacing 0.001 --nodata-value -32768
@@ -115,7 +142,7 @@ of 3,396,195 m, in the model returned with the ``-r mars`` option, that
 pixel would just be 5 m.
 
 You may want to compare the output to MOLA data. MOLA data is released
-in three ‘flavors,’ namely: Topography, Radius, and Areoid. The MOLA
+in three 'flavors', namely: Topography, Radius, and Areoid. The MOLA
 Topography data product that most people use is just the MOLA Radius
 product with the MOLA Areoid product subtracted. Additionally, it is
 important to note that all of these data products have a reference value
@@ -184,7 +211,7 @@ If you attempt to derive science results from an ASP-produced terrain
 model with the default DEM spacing, expect serious questions from
 reviewers.
 
-Using with LAS or CSV Clouds
+Using with LAS or CSV clouds
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The ``point2dem`` program can take as inputs point clouds in LAS and CSV
@@ -351,6 +378,11 @@ Command-line options for point2dem
 --proj-lon <float>
     The center of projection longitude (if applicable).
 
+--auto-proj-center
+    Automatically compute the projection center, when the projection is
+    stereographic, etc. This overrides the values of ``--proj-lat`` and
+    ``--proj-lon``. 
+
 --proj-scale <float>
     The projection scale (if applicable).
 

src/asp/Camera/BundleAdjustCamera.cc

@@ -758,7 +758,8 @@ void asp::init_camera_using_gcp(boost::shared_ptr<vw::ba::ControlNetwork> const&
 
   // Sanity check
   if (camera_models.size() != 1) 
-    vw::vw_throw(vw::ArgumentErr() << "Cannot initialize more than a camera at a time using GCP. "
+    vw::vw_throw(vw::ArgumentErr() 
+                 << "Cannot initialize more than a camera at a time using GCP. "
                  << "Consider using --transform-cameras-with-shared-gcp or "
                  << "--transform-cameras-using-gcp.\n");
 

src/asp/Camera/BundleAdjustCamera.h

@@ -629,8 +629,7 @@ void saveConvergenceAngles(std::string const& conv_angles_file,
                            std::vector<asp::MatchPairStats> const& convAngles,
                            std::vector<std::string> const& imageFiles);
 
-// Save stats of horizontal and vertical errors propagated from cameras
-// to triangulation
+// Save propagated horizontal and vertical errors
 void saveHorizVertErrors(std::string const& horiz_vert_errors_file,
                          std::vector<asp::HorizVertErrorStats> const& horizVertErrors,
                          std::vector<std::string> const& imageFiles);