Skip to content

Commit

Permalink
Aztimuth and slant range corrections for CUDA geocode (#1423)
Browse files Browse the repository at this point in the history 
* az and srg corrections for CUDA geocode
* corrections incorporated into masking geogrid radar indices
* update geocode function call
* workflow CUDA geocode object init consolidated inside gpu_geocode_rasters
* more comprehensive unit tests
* output to temporary directory
* check number of invalids in tests
* removed redundant lines_per_block
  • Loading branch information
Liang Yu authored and GitHub Enterprise committed Sep 25, 2023
1 parent f685460 commit 4cd3505
Show file tree
Hide file tree
Showing 7 changed files with 514 additions and 154 deletions.
150 changes: 104 additions & 46 deletions cxx/isce3/cuda/geocode/Geocode.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#include <isce3/core/Projections.h>
#include <isce3/cuda/container/RadarGeometry.h>
#include <isce3/cuda/core/OrbitView.h>
#include <isce3/cuda/core/gpuLUT2d.h>
#include <isce3/cuda/core/gpuProjections.h>
#include <isce3/cuda/except/Error.h>
#include <isce3/cuda/geometry/gpuDEMInterpolator.h>
Expand Down Expand Up @@ -45,21 +44,24 @@ namespace isce3::cuda::geocode {
/** Coverts a batch of rows from input geogrid into radar coordinates. Outputs
* the pixel-space coordinates (x,y) of each resulting (range, azimuth) pair
* with respect to specified radargrid, as well as mask invalid pixels (out
* of bounds or failed to converge in geo2rdr).
* of bounds, failed to converge in geo2rdr, or not contained in
* azTimeCorrections, sRangeCorrections, or nativeDoppler).
*
* \param[out] rdr_x pointer to device_vector of computed radar grid
* x / slant range indices
* \param[out] rdr_y pointer to device_vector of computed radar grid
* y / azimuth time indices
* \param[out] mask pointer to device_vector mask of valid
* pixels. Mask follows numpy mask_array
* convention where True is masked
* \param[out] mask pointer to device_vector mask of valid pixels
* Mask follows numpy mask_array convention where
* True is masked.
* \param[in] ellipsoid Ellipsoid based on output geogrid coordinate
* system
* \param[in] orbit Orbit associated with radar data
* \param[in] dem DEM interpolator. Maybe be of different
* coordinate system than output geogrid.
* \param[in] doppler doppler
* \param[in] nativeDoppler doppler centroid of data associated with radar
* grid, in Hz, as a function of azimuth and range
* \param[in] imageDoppler image grid doppler
* \param[in] wvl wavelength
* \param[in] side look side
* \param[in] geo2rdr_params geo2rdr params
Expand All @@ -69,15 +71,23 @@ namespace isce3::cuda::geocode {
* \param[in] block_size Number of elements in a block
* \param[in] proj Projection used to covert geogrid XYZ to LLH
* of output coordinate system
* \param[in] azTimeCorrection geo2rdr azimuth additive correction, in
* seconds, as a function of azimuth and range
* \param[in] sRangeCorrection geo2rdr slant range additive correction, in
* seconds, as a function of azimuth and range
*/
__global__ void geoToRdrIndices(double* rdr_x, double* rdr_y, bool* mask,
const isce3::core::Ellipsoid ellipsoid, const DeviceOrbitView orbit,
isce3::cuda::geometry::gpuDEMInterpolator dem,
const DeviceLUT2d<double> doppler, const double wvl,
const DeviceLUT2d<double> nativeDoppler,
const DeviceLUT2d<double> imageDoppler,
const double wvl,
const isce3::core::LookSide side, const Geo2RdrParams geo2rdr_params,
const isce3::product::GeoGridParameters geogrid,
const RadarGridParamsLite radargrid, const size_t line_start,
const size_t block_size, isce3::cuda::core::ProjectionBase** proj)
const size_t block_size, isce3::cuda::core::ProjectionBase** proj,
const DeviceLUT2d<double> azTimeCorrection,
const DeviceLUT2d<double> sRangeCorrection)
{
// thread index (1d grid of 1d blocks)
const auto tid = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
Expand All @@ -97,25 +107,47 @@ __global__ void geoToRdrIndices(double* rdr_x, double* rdr_y, bool* mask,
(*proj)->inverse(xyz, llh);
llh[2] = dem.interpolateLonLat(llh[0], llh[1]);

// geo2rdr slant range
double r;
// geo2rdr azimuth time initial value
double t = radargrid.sensing_mid;
// slant range to be computed by geo2rdr
double srange;
// initial guess for azimuth time to be updated by geo2rdr
double aztime = radargrid.sensing_mid;

// returns 0 if geo2rdr converges else 1
int converged = isce3::cuda::geometry::geo2rdr(llh, ellipsoid, orbit,
doppler, &t, &r, wvl, side, geo2rdr_params.threshold,
geo2rdr_params.maxiter, geo2rdr_params.delta_range);
imageDoppler, &aztime, &srange, wvl, side,
geo2rdr_params.threshold, geo2rdr_params.maxiter,
geo2rdr_params.delta_range);

// Bool to track if aztime and slant range is not contained in any of the
// following LUT2d's: azTimeCorrection, sRangeCorrection, and nativeDoppler.
// When set to true, this geogrid point will be marked as invalid
bool not_in_lut = false;

// Apply timing corrections
if (azTimeCorrection.contains(aztime, srange)) {
const auto aztimeCor = azTimeCorrection.eval(aztime, srange);
aztime += aztimeCor;
} else
not_in_lut = true;

if (sRangeCorrection.contains(aztime, srange)) {
const auto srangeCor = sRangeCorrection.eval(aztime, srange);
srange += srangeCor;
} else
not_in_lut = true;

if (!nativeDoppler.contains(aztime, srange))
not_in_lut = true;

// convert aztime and range to indices
double y = (t - radargrid.sensing_start) * radargrid.prf;
double x = (r - radargrid.starting_range) / radargrid.range_pixel_spacing;
double y = (aztime - radargrid.sensing_start) * radargrid.prf;
double x = (srange - radargrid.starting_range) / radargrid.range_pixel_spacing;

// check if indinces in bounds and set accordingly
// check if indinces in bounds of radar grid
const bool not_in_rdr_grid =
y < 0 || y >= radargrid.length || x < 0 || x >= radargrid.width;

const bool invalid_index = not_in_rdr_grid || converged == 0;
const bool invalid_index = not_in_rdr_grid || converged == 0 || not_in_lut;
rdr_y[tid] = invalid_index ? 0.0 : y;
rdr_x[tid] = invalid_index ? 0.0 : x;
mask[tid] = invalid_index;
Expand Down Expand Up @@ -287,8 +319,7 @@ void invalidValueCheck(const int dtype, const double invalid_value)
__host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid,
const isce3::container::RadarGeometry& rdr_geom,
const size_t lines_per_block) :
_lines_per_block(lines_per_block),
_n_blocks((geogrid.length() + _lines_per_block - 1) / _lines_per_block),
_n_blocks((geogrid.length() + lines_per_block - 1) / lines_per_block),
_geogrid(geogrid),
_rdr_geom(rdr_geom),
_ellipsoid(isce3::core::makeProjection(_geogrid.epsg())->ellipsoid()),
Expand All @@ -298,23 +329,23 @@ __host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid,
// correctly
_az_first_line(_rdr_geom.radarGrid().length() - 1),
_az_last_line(0),
_range_first_pixel(_rdr_geom.radarGrid().width() - 1),
_range_first_pixel(rdr_geom.radarGrid().width() - 1),
_range_last_pixel(0),
_geo_block_length(_lines_per_block),
_geo_block_length(lines_per_block),
_proj_handle(geogrid.epsg())
{
// init light weight radar grid
_radar_grid.sensing_start = _rdr_geom.radarGrid().sensingStart();
_radar_grid.sensing_mid = _rdr_geom.radarGrid().sensingMid();
_radar_grid.prf = _rdr_geom.radarGrid().prf();
_radar_grid.starting_range = _rdr_geom.radarGrid().startingRange();
_radar_grid.range_pixel_spacing = _rdr_geom.radarGrid().rangePixelSpacing();
_radar_grid.length = _rdr_geom.gridLength();
_radar_grid.width = _rdr_geom.gridWidth();
_radar_grid.sensing_start = rdr_geom.radarGrid().sensingStart();
_radar_grid.sensing_mid = rdr_geom.radarGrid().sensingMid();
_radar_grid.prf = rdr_geom.radarGrid().prf();
_radar_grid.starting_range = rdr_geom.radarGrid().startingRange();
_radar_grid.range_pixel_spacing = rdr_geom.radarGrid().rangePixelSpacing();
_radar_grid.length = rdr_geom.gridLength();
_radar_grid.width = rdr_geom.gridWidth();

// Determine max number of elements per block. Last block to be processed
// may not contain the max number of elements.
auto n_elem = _lines_per_block * _geogrid.width();
auto n_elem = _geo_block_length * _geogrid.width();

// Resize all device vectors to max block size.
_radar_x.resize(n_elem);
Expand All @@ -324,8 +355,12 @@ __host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid,

void Geocode::setBlockRdrCoordGrid(const size_t block_number,
Raster& dem_raster,
const isce3::core::dataInterpMethod _dem_interp_method,
const isce3::geometry::detail::Geo2RdrParams _geo2rdr_params)
const isce3::core::dataInterpMethod dem_interp_method,
const isce3::cuda::container::RadarGeometry& dev_rdr_geom,
const isce3::geometry::detail::Geo2RdrParams geo2rdr_params,
const DeviceLUT2d<double>& nativeDoppler,
const DeviceLUT2d<double>& azTimeCorrection,
const DeviceLUT2d<double>& sRangeCorrection)
{
// make sure block index does not exceed actual number of blocks
if (block_number >= _n_blocks) {
Expand All @@ -334,11 +369,10 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number,
}

// Get block extents (of the geocoded grid)
_line_start = block_number * _lines_per_block;
_line_start = block_number * _geo_block_length;

// Set block sizes for everything but last block
size_t block_size = _geo_block_length * _geogrid.width();
_geo_block_length = _lines_per_block;

// Adjust for last block sized assuming it is sized differently than others
if (block_number == (_n_blocks - 1)) {
Expand All @@ -356,12 +390,9 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number,
isce3::geometry::DEMInterpolator host_dem_interp =
isce3::geometry::DEMRasterToInterpolator(dem_raster, _geogrid,
_line_start, _geo_block_length, _geogrid.width(),
dem_margin_in_pixels, _dem_interp_method);
dem_margin_in_pixels, dem_interp_method);
isce3::cuda::geometry::gpuDEMInterpolator dev_dem_interp(host_dem_interp);

// copy RadarGeometry to device
isce3::cuda::container::RadarGeometry dev_rdr_geom(_rdr_geom);

// Create geogrid on device
{
const unsigned threads_per_block = 256;
Expand All @@ -370,16 +401,20 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number,

geoToRdrIndices<<<n_blocks, threads_per_block>>>(_radar_x.data().get(),
_radar_y.data().get(), _mask.data().get(), _ellipsoid,
dev_rdr_geom.orbit(), dev_dem_interp, dev_rdr_geom.doppler(),
dev_rdr_geom.orbit(), dev_dem_interp,
nativeDoppler, dev_rdr_geom.doppler(),
dev_rdr_geom.wavelength(), dev_rdr_geom.lookSide(),
_geo2rdr_params, _geogrid, _radar_grid, _line_start, block_size,
_proj_handle.get_proj());
geo2rdr_params, _geogrid, _radar_grid, _line_start, block_size,
_proj_handle.get_proj(), azTimeCorrection, sRangeCorrection);

checkCudaErrors(cudaPeekAtLastError());
checkCudaErrors(cudaDeviceSynchronize());
}

// find index of min and max values in x/range
// Mask is utilized to _radar_x and _radar_y when determining respective
// min and max values which are then used to determine corresponding radar
// data to read from radar raster
// Find index of unmasked min and max values in x/range
const auto [rdr_x_min, rdr_x_max] = masked_minmax(_radar_x, _mask);

_range_first_pixel = std::min(_rdr_geom.radarGrid().width() - 1,
Expand All @@ -388,7 +423,7 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number,
_range_last_pixel = std::max(static_cast<size_t>(0),
static_cast<size_t>(std::ceil(rdr_x_max) - 1));

// find index of min and max values in y/azimuth
// Find index of unmasked min and max values in y/azimuth
const auto [rdr_y_min, rdr_y_max] = masked_minmax(_radar_y, _mask);

_az_first_line = std::min(_rdr_geom.radarGrid().length() - 1,
Expand Down Expand Up @@ -586,6 +621,9 @@ void Geocode::geocodeRasters(
const std::vector<GDALDataType>& raster_datatypes,
const std::vector<double>& invalid_values_double,
Raster& dem_raster,
const isce3::core::LUT2d<double>& hostNativeDoppler,
const isce3::core::LUT2d<double>& hostAzTimeCorrection,
const isce3::core::LUT2d<double>& hostSRangeCorrection,
const isce3::core::dataInterpMethod dem_interp_method,
const double threshold, const int maxiter, const double dr)
{
Expand All @@ -603,7 +641,18 @@ void Geocode::geocodeRasters(
std::vector<std::any> _invalid_values;

// geo2rdr params used in radar index calculation
isce3::geometry::detail::Geo2RdrParams _geo2rdr_params{threshold, maxiter, dr};
isce3::geometry::detail::Geo2RdrParams geo2rdr_params{threshold, maxiter, dr};

// 2D LUT Doppler of the SLC image
auto nativeDoppler = DeviceLUT2d<double>(hostNativeDoppler);

// geo2rdr azimuth additive correction, in seconds, as a function of
// azimuth and range
auto azTimeCorrection = DeviceLUT2d<double>(hostAzTimeCorrection);;

// geo2rdr slant range additive correction, in seconds, as a function of
// azimuth and range
auto sRangeCorrection = DeviceLUT2d<double>(hostSRangeCorrection);

// Prepare interpolator handler and invalid value for each raster
for (size_t i_raster = 0; i_raster < n_raster_pairs; ++i_raster) {
Expand Down Expand Up @@ -680,12 +729,21 @@ void Geocode::geocodeRasters(
}
}

// copy RadarGeometry to device
const auto dev_rdr_geom = isce3::cuda::container::RadarGeometry(_rdr_geom);

// iterate over blocks
for (size_t i_block = 0; i_block < _n_blocks; ++i_block) {

// compute then set radar coords for each geocode obj for curret block
setBlockRdrCoordGrid(i_block, dem_raster, dem_interp_method,
_geo2rdr_params);
setBlockRdrCoordGrid(i_block,
dem_raster,
dem_interp_method,
dev_rdr_geom,
geo2rdr_params,
nativeDoppler,
azTimeCorrection,
sRangeCorrection);

// iterate over rasters and geocode with associated datatype and
// interpolation method
Expand Down
49 changes: 40 additions & 9 deletions cxx/isce3/cuda/geocode/Geocode.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,16 +11,21 @@
#include <thrust/device_vector.h>

#include <isce3/container/RadarGeometry.h>
#include <isce3/cuda/container/forward.h>
#include <isce3/cuda/core/InterpolatorHandle.h>
#include <isce3/cuda/core/ProjectionBaseHandle.h>
#include <isce3/cuda/core/gpuInterpolator.h>
#include <isce3/cuda/core/gpuLUT2d.h>
#include <isce3/cuda/core/gpuProjections.h>
#include <isce3/geometry/detail/Geo2Rdr.h>
#include <isce3/io/Raster.h>
#include <isce3/product/GeoGridParameters.h>

namespace isce3::cuda::geocode {

template<typename T>
using DeviceLUT2d = isce3::cuda::core::gpuLUT2d<T>;

/* light weight radar grid container */
struct RadarGridParamsLite {
double sensing_start;
Expand Down Expand Up @@ -83,6 +88,15 @@ class Geocode {
* initialize templated interpolators.
* \param[in] invalid_values Invalid values for each raster
* \param[in] dem_raster DEM used to calculate radar grid indices
* \param[in] nativeDoppler Doppler centroid of data in Hz associated
* radar grid, as a function azimuth and
* range
* \param[in] azTimeCorrection geo2rdr azimuth additive correction, in
* seconds, as a function of azimuth and
* range
* \param[in] sRangeCorrection geo2rdr slant range additive
* correction, in seconds, as a function
* of azimuth and range
* \param[in] dem_interp_method DEMinterpolation method
* \param[in] threshold Convergence threshold for geo2rdr
* \param[in] maxiter Maximum iterations for geo2rdr
Expand All @@ -96,15 +110,16 @@ class Geocode {
const std::vector<GDALDataType>& raster_datatypes,
const std::vector<double>& invalid_values,
Raster& dem_raster,
const isce3::core::LUT2d<double>& hostNativeDoppler = {},
const isce3::core::LUT2d<double>& hostAzTimeCorrection = {},
const isce3::core::LUT2d<double>& hostSRangeCorrection = {},
const isce3::core::dataInterpMethod dem_interp_method =
isce3::core::BIQUINTIC_METHOD,
const double threshold = 1e-8, const int maxiter = 50,
const double threshold = 1e-8,
const int maxiter = 50,
const double dr = 10);

private:
// number of lines to be processed in a block
size_t _lines_per_block;

// total number of blocks necessary to geocoding a provided geogrid
size_t _n_blocks;

Expand Down Expand Up @@ -156,7 +171,8 @@ class Geocode {
* \param[in] raster_datatypes GDAL type of each raster represented by
* equivalent int. Necessary to correctly
* initialize templated interpolators.
* \param[in] invalid_values Invalid values for each raster
* \param[in] invalid_values_double Invalid values for each raster as double, to
be converted later to the correct type
*/
void ensureRasterConsistency(
const std::vector<std::reference_wrapper<isce3::io::Raster>>& output_rasters,
Expand All @@ -173,14 +189,29 @@ class Geocode {
* \param[in] block_number Index of block where radar grid coordinates
* coordinates will be calculated and set
* \param[in] dem_raster DEM used to calculate radar grid indices
* \param[in] _dem_interp_method DEMinterpolation method
* \param[in] _geo2rdr_params Parameters used by geo2rdr to compute
* \param[in] dem_interp_method DEMinterpolation method
* \param[in] dev_rdr_geom Radar geometry describing input raster
* stored on-device
* \param[in] geo2rdr_params Parameters used by geo2rdr to compute
* geogrids radar grid coordinates
* \param[in] nativeDoppler Doppler centroid of data in Hz associated
* radar grid, as a function azimuth and
* range
* \param[in] azTimeCorrection geo2rdr azimuth additive correction, in
* seconds, as a function of azimuth and
* range
* \param[in] sRangeCorrection geo2rdr slant range additive
* correction, in seconds, as a function
* of azimuth and range
*/
void setBlockRdrCoordGrid(const size_t block_number,
Raster& dem_raster,
const isce3::core::dataInterpMethod _dem_interp_method,
const isce3::geometry::detail::Geo2RdrParams _geo2rdr_params);
const isce3::core::dataInterpMethod dem_interp_method,
const isce3::cuda::container::RadarGeometry& dev_rdr_geom,
const isce3::geometry::detail::Geo2RdrParams geo2rdr_params,
const DeviceLUT2d<double>& nativeDoppler,
const DeviceLUT2d<double>& azTimeCorrection,
const DeviceLUT2d<double>& sRangeCorrection);

/** Geocode a block of raster according to grid last set in
* setBlockRdrCoordGrid. Only operates on 1st raster band of input/output
Expand Down
Loading

0 comments on commit 4cd3505

Please sign in to comment.