Merge pull request isce-framework#6 from jungkyoj/develop

merged the most recent develop branch into split_main_band

VERSION.txt

@@ -1 +1 @@
-0.4.1-dev
+0.5.0-dev

cxx/isce3/Headers.cmake

@@ -83,7 +83,6 @@ focus/Presum.icc
 focus/RangeComp.h
 geocode/baseband.h
 geocode/geocodeSlc.h
-geocode/interpolate.h
 geocode/loadDem.h
 geometry/DEMInterpolator.h
 geometry/forward.h

cxx/isce3/Sources.cmake

@@ -43,7 +43,6 @@ focus/Presum.cpp
 focus/RangeComp.cpp
 geocode/baseband.cpp
 geocode/geocodeSlc.cpp
-geocode/interpolate.cpp
 geocode/loadDem.cpp
 geometry/DEMInterpolator.cpp
 geometry/Geo2rdr.cpp

cxx/isce3/cuda/geocode/Geocode.cu

@@ -5,6 +5,7 @@
 #include <thrust/functional.h>
 #include <thrust/host_vector.h>
 
+#include <gdal_priv.h>
 #include <pyre/journal.h>
 
 #include <isce3/core/Ellipsoid.h>
@@ -414,6 +415,57 @@ void Geocode::geocodeRasterBlock(Raster& output_raster, Raster& input_raster)
  _geogrid.width(), _geo_block_length, 1);
 }
 
+void Geocode::geocodeRasters(
+ std::vector<std::reference_wrapper<isce3::io::Raster>> output_rasters,
+ std::vector<std::reference_wrapper<isce3::io::Raster>> input_rasters)
+{
+ // check if vectors are of same length
+ if (output_rasters.size() != input_rasters.size()) {
+ throw isce3::except::LengthError(
+ ISCE_SRCINFO(), "number of input and output rasters not equal");
+ }
+
+ const auto n_raster_pairs = output_rasters.size();
+
+ // iterate over blocks
+ for (size_t i_block = 0; i_block < _n_blocks; ++i_block) {
+
+ // set radar coords for each geocode obj for curret block
+ setBlockRdrCoordGrid(i_block);
+
+ for (size_t i_raster = 0; i_raster < n_raster_pairs; ++i_raster)
+ {
+ const int dtype = input_rasters[i_raster].get().dtype();
+ switch (dtype) {
+ case GDT_Float32: {
+ geocodeRasterBlock<float>(
+ output_rasters[i_raster], input_rasters[i_raster]);
+ break; }
+ case GDT_CFloat32: {
+ geocodeRasterBlock<thrust::complex<float>>(
+ output_rasters[i_raster], input_rasters[i_raster]);
+ break;}
+ case GDT_Float64: {
+ geocodeRasterBlock<double>(
+ output_rasters[i_raster], input_rasters[i_raster]);
+ break; }
+ case GDT_CFloat64: {
+ geocodeRasterBlock<thrust::complex<double>>(
+ output_rasters[i_raster], input_rasters[i_raster]);
+ break;}
+ case GDT_Byte: {
+ geocodeRasterBlock<unsigned char>(
+ output_rasters[i_raster], input_rasters[i_raster]);
+ break;}
+ default: {
+ throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+ "unsupported datatype");
+ }
+ }
+ }
+ }
+}
+
 #define EXPLICIT_INSTATIATION(T) \
  template void Geocode::geocodeRasterBlock<T>( \
  Raster & output_raster, Raster & input_raster);

cxx/isce3/cuda/geocode/Geocode.h

@@ -1,5 +1,8 @@
 #pragma once
 
+#include <functional>
+#include <vector>
+
 #include <isce3/core/forward.h>
 #include <isce3/geometry/forward.h>
 
@@ -104,6 +107,16 @@ class Geocode {
  void geocodeRasterBlock(
  isce3::io::Raster& output_raster, isce3::io::Raster& input_raster);
 
+ /** Geocode rasters with a shared geogrid. Block processing handled
+ * internally in function.
+ *
+ * \param[in] output_rasters Geocoded rasters
+ * \param[in] input_rasters Rasters to be geocoded
+ */
+ void geocodeRasters(
+ std::vector<std::reference_wrapper<isce3::io::Raster>> output_rasters,
+ std::vector<std::reference_wrapper<isce3::io::Raster>> input_rasters);
+
  size_t numBlocks() const { return _n_blocks; }
  size_t linesPerBlock() const { return _lines_per_block; }
 

cxx/isce3/cuda/geometry/gpuGeometry.cu

@@ -74,10 +74,14 @@ int geo2rdr(const Vec3& inputLLH, const isce3::core::Ellipsoid& ellipsoid,
  // Use mid-orbit epoch as initial guess
  aztime = orbit.midTime();
 
+ // Newton step, init to zero.
+ double dt = 0.0;
+
  // Begin iterations
  int converged = 0;
- double slantRange_old = 0.0;
  for (int i = 0; i < maxIter; ++i) {
+ // Apply Newton step from previous iteration.
+ aztime -= dt;
 
  // Interpolate the orbit to current estimate of azimuth time
  Vec3 pos, vel;
@@ -95,16 +99,6 @@ int geo2rdr(const Vec3& inputLLH, const isce3::core::Ellipsoid& ellipsoid,
  return converged;
  }
 
- // Check convergence
- if (std::abs(slantRange - slantRange_old) < threshold) {
- converged = 1;
- *slantRange_result = slantRange;
- *aztime_result = aztime;
- return converged;
- } else {
- slantRange_old = slantRange;
- }
-
  // Compute doppler
  const double dopfact = dr.dot(vel);
  const double fdop = doppler.eval(slantRange) * dopscale;
@@ -119,8 +113,17 @@ int geo2rdr(const Vec3& inputLLH, const isce3::core::Ellipsoid& ellipsoid,
  const double c2 = (fdop / slantRange) + fdopder;
  const double fnprime = c1 + c2 * dopfact;
 
- // Update guess for azimuth time
- aztime -= fn / fnprime;
+ // Compute Newton step. Don't apply until start of next iteration so
+ // that returned (slantRange, aztime) are always consistent.
+ dt = fn / fnprime;
+
+ // Check convergence
+ if (std::abs(dt) < threshold) {
+ converged = 1;
+ *slantRange_result = slantRange;
+ *aztime_result = aztime;
+ return converged;
+ }
  }
 
  // If we reach this point, no convergence for specified threshold

cxx/isce3/geocode/GeocodePolygon.h

@@ -33,7 +33,7 @@ class GeocodePolygon {
  * @param[in] orbit Orbit
  * @param[in] input_dop Doppler LUT associated with the radar grid
  * @param[in] dem_raster Input DEM raster
- * @param[in] threshold Distance threshold for convergence
+ * @param[in] threshold Azimuth time threshold for convergence (s)
  * @param[in] num_iter Maximum number of Newton-Raphson iterations
  * @param[in] delta_range Step size used for computing Doppler derivative
  */