diff --git a/NEWS.rst b/NEWS.rst index 7fceefcb3..ebcec7542 100644 --- a/NEWS.rst +++ b/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`). diff --git a/docs/next_steps.rst b/docs/next_steps.rst index a3940a3c7..e1e7464e4 100644 --- a/docs/next_steps.rst +++ b/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 --proj-lat ... +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 --proj-lat 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, diff --git a/docs/tools/bundle_adjust.rst b/docs/tools/bundle_adjust.rst index 48cd40ad4..169873a41 100644 --- a/docs/tools/bundle_adjust.rst +++ b/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 diff --git a/docs/tools/point2dem.rst b/docs/tools/point2dem.rst index ea0dab6e5..d2da61d36 100644 --- a/docs/tools/point2dem.rst +++ b/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 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 The projection scale (if applicable). diff --git a/src/asp/Camera/BundleAdjustCamera.cc b/src/asp/Camera/BundleAdjustCamera.cc index d8d46cf1c..fecdac1c0 100644 --- a/src/asp/Camera/BundleAdjustCamera.cc +++ b/src/asp/Camera/BundleAdjustCamera.cc @@ -758,7 +758,8 @@ void asp::init_camera_using_gcp(boost::shared_ptr 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"); diff --git a/src/asp/Camera/BundleAdjustCamera.h b/src/asp/Camera/BundleAdjustCamera.h index 163048b06..296035bef 100644 --- a/src/asp/Camera/BundleAdjustCamera.h +++ b/src/asp/Camera/BundleAdjustCamera.h @@ -629,8 +629,7 @@ void saveConvergenceAngles(std::string const& conv_angles_file, std::vector const& convAngles, std::vector 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 const& horizVertErrors, std::vector const& imageFiles); diff --git a/src/asp/Core/PointUtils.cc b/src/asp/Core/PointUtils.cc index c9be40f57..1ec5ac927 100644 --- a/src/asp/Core/PointUtils.cc +++ b/src/asp/Core/PointUtils.cc @@ -213,7 +213,6 @@ namespace asp{ if (!valid) vw_throw(ArgumentErr() << "Fatal error reading PCD file: " << m_pcd_file); } - PcdReader::PcdReader(std::string const & pcd_file) : m_pcd_file(pcd_file), m_has_valid_point(false){ @@ -292,7 +291,7 @@ namespace asp{ /// important that the writer invoking this image be single-threaded, /// as we read from the las file sequentially. template - class LasOrCsvToTif_Class : public ImageViewBase< LasOrCsvToTif_Class > { + class LasOrCsvToTif_Class: public ImageViewBase< LasOrCsvToTif_Class> { typedef typename ImageT::pixel_type PixelT; @@ -373,7 +372,7 @@ namespace asp{ }; // End class LasOrCsvToTif_Class -} // namespace asp +} // end namespace asp //------------------------------------------------------------------------------------------ // Class CsvConv functions @@ -1283,20 +1282,37 @@ vw::BBox3 asp::pointcloud_bbox(vw::ImageViewRef const& point_image, return result; } -// Find the average longitude for a given point image with lon, lat, height values -double asp::find_avg_lon(ImageViewRef const& point_image){ +// Find the average longitude for a given point image with lon, lat, height values. +// The computation is done in cartesian coordinates. Then look at the produced +// x value. That correlates with a hemisphere. +double asp::find_avg_lon(ImageViewRef const& point_image) { Stopwatch sw; sw.start(); int32 subsample_amt = int32(norm_2(Vector2(point_image.cols(), point_image.rows()))/32.0); - if (subsample_amt < 1 ) + if (subsample_amt < 1) subsample_amt = 1; - PixelAccumulator > mean_accum; - for_each_pixel( subsample(point_image, subsample_amt), - mean_accum, - TerminalProgressCallback("asp","Statistics: ") ); - Vector3 avg_location = mean_accum.value(); + + ImageViewRef sub_image = subsample(point_image, subsample_amt); + + // Accumulate valid values + Vector3 avg_location(0, 0, 0); + double num = 0; + for (int32 col = 0; col < sub_image.cols(); col++) { + for (int32 row = 0; row < sub_image.rows(); row++) { + auto pix = sub_image(col, row); + if (pix == Vector3()) + continue; + avg_location += pix; + num += 1.0; + } + } + + // TODO(oalexan1): What it do if we hit no valid value? + if (num > 0) + avg_location /= num; + double avg_lon = avg_location.x() >= 0 ? 0 : 180; sw.stop(); vw_out(DebugMessage,"asp") << "Statistics time: " << sw.elapsed_seconds() << std::endl; @@ -1304,6 +1320,62 @@ double asp::find_avg_lon(ImageViewRef const& point_image){ return avg_lon; } +// Find the median longitude and latitude for a subset of the point cloud +void asp::median_lon_lat(vw::ImageViewRef const& point_image, + vw::cartography::GeoReference const& georef, + double & lon, double & lat) { + + vw_out() << "Estimating the median longitude and latitude for the cloud.\n"; + + // Initialize the outputs + lon = 0.0; lat = 0.0; + + // Try to subsample with these amounts + std::vector sub = {32.0, 16.0}; + + Stopwatch sw; + sw.start(); + + // Iterate over sub + bool success = false; + for (size_t s = 0; s < sub.size(); s++) { + + int32 subsample_amt = int32(norm_2(Vector2(point_image.cols(), + point_image.rows()))/sub[s]); + if (subsample_amt < 1) + subsample_amt = 1; + + ImageViewRef sub_image = subsample(point_image, subsample_amt); + + // Accumulate valid values + std::vector lons, lats; + for (int32 col = 0; col < sub_image.cols(); col++) { + for (int32 row = 0; row < sub_image.rows(); row++) { + Vector3 ecef = sub_image(col, row); + if (ecef == Vector3()) + continue; + Vector3 llh = georef.datum().cartesian_to_geodetic(ecef); + lons.push_back(llh[0]); + lats.push_back(llh[1]); + } + } + if (lons.empty() || lats.empty()) + continue; + + lon = vw::math::destructive_median(lons); + lat = vw::math::destructive_median(lats); + success = true; + break; + } + + if (!success) + vw_throw(ArgumentErr() << "Could not find a valid median longitude and latitude. " + << "Check if your cloud is empty.\n"); + + sw.stop(); + vw_out(DebugMessage,"asp") << "Median longitude and latitude elapsed time: " + << sw.elapsed_seconds() << std::endl; +} /// Analyze a file name to determine the file type std::string asp::get_cloud_type(std::string const& file_name){ diff --git a/src/asp/Core/PointUtils.h b/src/asp/Core/PointUtils.h index d64c71cc4..296b13a7e 100644 --- a/src/asp/Core/PointUtils.h +++ b/src/asp/Core/PointUtils.h @@ -435,6 +435,11 @@ vw::ImageViewRef form_point_cloud_composite(std::vector con // Find the average longitude for a given point image with lon, lat, height values double find_avg_lon(vw::ImageViewRef const& point_image); +// Find the median longitude and latitude for a subset of the point cloud +void median_lon_lat(vw::ImageViewRef const& point_image, + vw::cartography::GeoReference const& georef, + double & lon, double & lat); + std::string get_cloud_type(std::string const& file_name); // Find the number of channels in the point clouds. diff --git a/src/asp/Tools/point2dem.cc b/src/asp/Tools/point2dem.cc index e09426d6a..f43d5affc 100644 --- a/src/asp/Tools/point2dem.cc +++ b/src/asp/Tools/point2dem.cc @@ -15,9 +15,8 @@ // limitations under the License. // __END_LICENSE__ - -/// \file point2dem.cc -/// +// \file point2dem.cc +// #include #include @@ -87,7 +86,7 @@ struct Options : vw::GdalWriteOptions { std::string csv_format_str, csv_proj4_str, filter; double search_radius_factor, sigma_factor, default_grid_size_multiplier; bool use_surface_sampling; - bool has_las_or_csv_or_pcd; + bool has_las_or_csv_or_pcd, auto_proj_center; Vector2i max_output_size; bool input_is_projected; @@ -103,7 +102,8 @@ struct Options : vw::GdalWriteOptions { max_valid_triangulation_error(0), erode_len(0), search_radius_factor(0), sigma_factor(0), default_grid_size_multiplier(1.0), use_surface_sampling(false), - has_las_or_csv_or_pcd(false), max_output_size(9999999, 9999999), input_is_projected(false){} + has_las_or_csv_or_pcd(false), max_output_size(9999999, 9999999), + auto_proj_center(false), input_is_projected(false) {} }; void parse_input_clouds_textures(std::vector const& files, @@ -191,11 +191,11 @@ void parse_input_clouds_textures(std::vector const& files, } -/// Convert any LAS or CSV files to ASP tif files. We do some binning -/// to make the spatial data more localized, to improve performance. -/// - We will later wipe these temporary tifs. +// Convert any LAS or CSV files to ASP tif files. We do some binning +// to make the spatial data more localized, to improve performance. +// - We will later wipe these temporary tifs. void las_or_csv_or_pcd_to_tifs(Options& opt, cartography::Datum const& datum, - std::vector & tmp_tifs){ + std::vector & tmp_tifs) { if (!opt.has_las_or_csv_or_pcd) return; @@ -247,7 +247,7 @@ void las_or_csv_or_pcd_to_tifs(Options& opt, cartography::Datum const& datum, } // No tif files exist. Find a reasonable value for the number of rows. - if (num_rows == 0){ + if (num_rows == 0) { std::int64_t max_num_pts = 0; for (int i = 0; i < num_files; i++){ std::string file = opt.pointcloud_files[i]; @@ -310,7 +310,7 @@ void las_or_csv_or_pcd_to_tifs(Options& opt, cartography::Datum const& datum, } // TODO: Move this somewhere? -/// Parses a string containing a list of numbers +// Parses a string containing a list of numbers void split_number_string(const std::string &input, std::vector &output) { // Get a space delimited string @@ -365,6 +365,9 @@ void handle_arguments(int argc, char *argv[], Options& opt) { ("utm", po::value(&opt.utm_zone), "Save using a UTM projection with the given zone.") ("proj-lat", po::value(&opt.proj_lat)->default_value(0), "The center of projection latitude (if applicable).") ("proj-lon", po::value(&opt.proj_lon)->default_value(0), "The center of projection longitude (if applicable).") + ("auto-proj-center", po::bool_switch(&opt.auto_proj_center)->default_value(false), + "Automatically compute the projection center, when the projection is stereographic, " + "etc. This overrides the values of --proj-lat and --proj-lon.") ("proj-scale", po::value(&opt.proj_scale)->default_value(1), "The projection scale (if applicable).") ("false-easting", po::value(&opt.false_easting)->default_value(0), "The projection false easting (if applicable).") ("false-northing", po::value(&opt.false_northing)->default_value(0), "The projection false northing (if applicable)."); @@ -546,7 +549,10 @@ void handle_arguments(int argc, char *argv[], Options& opt) { << "Invalid values were provided for outlier removal params.\n"); } - if (opt.max_valid_triangulation_error > 0){ + if (opt.input_is_projected && opt.auto_proj_center) + vw_throw(ArgumentErr() << "Cannot use both --input-is-projected and --auto-proj-center.\n"); + + if (opt.max_valid_triangulation_error > 0) { // Since the user passed in a threshold, will use that to rm // outliers, instead of using the percentage. opt.remove_outliers_with_pct = false; @@ -594,8 +600,13 @@ void handle_arguments(int argc, char *argv[], Options& opt) { else if (vm.count("gnomonic")) opt.projection = GNOMONIC; else if (vm.count("lambert-azimuthal")) opt.projection = LAMBERTAZIMUTHAL; else if (vm.count("utm")) opt.projection = UTM; - else opt.projection = PLATECARREE; // Default projection -} + else opt.projection = PLATECARREE; // Default + + // Sanity check + if (opt.auto_proj_center && opt.projection == PLATECARREE) + vw::vw_throw(ArgumentErr() << "No projection was set. Cannot use --auto-proj-center.\n"); + +} // end function handle_arguments() template ImageViewRef< PixelGray> @@ -628,7 +639,7 @@ generate_fsaa_raster(ImageViewBase const& rasterizer, Options const& opt return rasterizer_fsaa; } -namespace asp{ +namespace asp { // If the third component of a vector is NaN, mask that vector as invalid template @@ -691,12 +702,11 @@ namespace asp{ ErrorToNED(georef)); } - /// Write an image to disk while handling some common options. + // Write an image to disk while handling some common options. template void save_image(Options& opt, ImageT img, GeoReference const& georef, int hole_fill_len, std::string const& imgName){ - // When hole-filling is used, we need to look hole_fill_len beyond // the current block. If the block size is 256, and hole fill len // is big, like 512 or 1024, we end up processing a huge block @@ -729,8 +739,7 @@ namespace asp{ vw::cartography::write_gdal_image(output_file, img, georef, opt, tpc); } // End function save_image - - /// A class for combining the three channels of errors and finding their absolute values. + // A class for combining the three channels of errors and finding their absolute values. template class CombinedView : public ImageViewBase > { @@ -772,7 +781,6 @@ namespace asp{ return Vector3f(std::abs(error[0]), std::abs(error[1]), std::abs(error[2])); } - /// \cond INTERNAL typedef CombinedView prerasterize_type; inline prerasterize_type prerasterize(BBox2i const& bbox) const { @@ -785,7 +793,6 @@ namespace asp{ template inline void rasterize(DestT const& dest, BBox2i const& bbox) const { vw::rasterize(prerasterize(bbox), dest, bbox); } - /// \endcond }; template CombinedView combine_channels(double nodata_value, @@ -801,15 +808,15 @@ namespace asp{ return CombinedView(nodata_value, image1.impl(), image2.impl(), image3.impl()); } - /// Round pixels in given image to multiple of given scale. - /// Don't round nodata values. + // Round pixels in given image to multiple of given scale. + // Don't round nodata values. template struct RoundImagePixelsSkipNoData: public vw::ReturnFixedType { double m_scale, m_nodata; - RoundImagePixelsSkipNoData(double scale, double nodata) : m_scale(scale), m_nodata(nodata){ - } + RoundImagePixelsSkipNoData(double scale, double nodata): + m_scale(scale), m_nodata(nodata) {} PixelT operator() (PixelT const& pt) const { @@ -1189,12 +1196,37 @@ void do_software_rasterization_multi_spacing(const ImageViewRef& proj_p opt.out_prefix = base_out_prefix; // Restore the original value } +// Set the projection +void set_projection(Options const& opt, cartography::GeoReference & output_georef) { + + // This is an error + if (!opt.target_srs_string.empty()) + vw::vw_throw(ArgumentErr() + << "The --t_srs option must not be used when setting a projection.\n"); + + switch (opt.projection) { + case SINUSOIDAL: output_georef.set_sinusoidal (opt.proj_lon, opt.false_easting, opt.false_northing); break; + case MERCATOR: output_georef.set_mercator (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; + case TRANSVERSEMERCATOR: output_georef.set_transverse_mercator (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; + case ORTHOGRAPHIC: output_georef.set_orthographic (opt.proj_lat, opt.proj_lon, opt.false_easting, opt.false_northing); break; + case STEREOGRAPHIC: output_georef.set_stereographic (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; + case OSTEREOGRAPHIC: output_georef.set_oblique_stereographic(opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; + case GNOMONIC: output_georef.set_gnomonic (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; + case LAMBERTAZIMUTHAL: output_georef.set_lambert_azimuthal (opt.proj_lat, opt.proj_lon, opt.false_easting, opt.false_northing); break; + case UTM: output_georef.set_UTM(opt.utm_zone); break; + default: // Handles plate carree + break; + } + + return; +} + int main(int argc, char *argv[]) { Options opt; try { handle_arguments(argc, argv, opt); - // Set up the georeferencing information. We specify everything + // Set up the georeferencing information. We specify everything // here except for the affine transform, which is defined later once // we know the bounds of the orthorasterizer view. However, we can // still reproject the points in the point image without the affine @@ -1221,20 +1253,9 @@ int main(int argc, char *argv[]) { if (have_user_datum) output_georef.set_datum(user_datum); - - switch (opt.projection) { - case SINUSOIDAL: output_georef.set_sinusoidal (opt.proj_lon, opt.false_easting, opt.false_northing); break; - case MERCATOR: output_georef.set_mercator (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; - case TRANSVERSEMERCATOR: output_georef.set_transverse_mercator (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; - case ORTHOGRAPHIC: output_georef.set_orthographic (opt.proj_lat, opt.proj_lon, opt.false_easting, opt.false_northing); break; - case STEREOGRAPHIC: output_georef.set_stereographic (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; - case OSTEREOGRAPHIC: output_georef.set_oblique_stereographic(opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; - case GNOMONIC: output_georef.set_gnomonic (opt.proj_lat, opt.proj_lon, opt.proj_scale, opt.false_easting, opt.false_northing); break; - case LAMBERTAZIMUTHAL: output_georef.set_lambert_azimuthal (opt.proj_lat, opt.proj_lon, opt.false_easting, opt.false_northing); break; - case UTM: output_georef.set_UTM(opt.utm_zone); break; - default: // Handles plate carree - break; - } + + set_projection(opt, output_georef); + } else { // The user specified the target srs_string // Set the srs string into georef. @@ -1245,8 +1266,8 @@ int main(int argc, char *argv[]) { // If the input PROJ.4 string is empty, use the output one. if (!opt.csv_format_str.empty() && opt.csv_proj4_str.empty()) { opt.csv_proj4_str = output_georef.overall_proj4_str(); - std::cout << "The option --csv-proj4 was not specified. Using the output projection " - << "when interpreting csv files.\n"; + vw_out() << "The option --csv-proj4 was not specified. Using the output projection " + << "when interpreting csv files.\n"; } // Convert any input LAS or CSV files to ASP's point cloud tif format @@ -1299,7 +1320,13 @@ int main(int argc, char *argv[]) { // Forcing the georef object outside its comfort zone is not safe for all projections! if (output_georef.overall_proj4_str().find("+proj=aea") == std::string::npos) output_georef.set_lon_center(avg_lon < 100); - + + if (opt.auto_proj_center) { + // Find the median lon lat and reapply this to the georef + asp::median_lon_lat(point_image, output_georef, opt.proj_lon, opt.proj_lat); + set_projection(opt, output_georef); + } + // Convert xyz points to projected points // - The cartesian_to_geodetic call converts invalid (0,0,0,0) points to NaN, // which is checked for in the OrthoRasterizer class.