diff --git a/cxx/isce3/cuda/signal/gpuCrossMul.cu b/cxx/isce3/cuda/signal/gpuCrossMul.cu index e262d51cd..49c3690bf 100644 --- a/cxx/isce3/cuda/signal/gpuCrossMul.cu +++ b/cxx/isce3/cuda/signal/gpuCrossMul.cu @@ -119,8 +119,7 @@ __global__ void flatten_g(thrust::complex *ifg, /** Set number of range looks */ -void -gpuCrossmul::rangeLooks(int rngLks) { +void gpuCrossmul::rangeLooks(int rngLks) { if (rngLks < 1) { std::string error_msg = "ERROR CUDA crossmul range multilook < 1"; throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg); @@ -134,8 +133,7 @@ gpuCrossmul::rangeLooks(int rngLks) { } /** Set number of azimuth looks */ -void -gpuCrossmul::azimuthLooks(int azLks) { +void gpuCrossmul::azimuthLooks(int azLks) { if (azLks < 1) { std::string error_msg = "ERROR CUDA crossmul azimuth multilook < 1"; throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg); @@ -148,9 +146,8 @@ gpuCrossmul::azimuthLooks(int azLks) { _multiLookEnabled = false; } -void -gpuCrossmul::doppler(isce3::core::LUT1d refDoppler, - isce3::core::LUT1d secDoppler) +void gpuCrossmul::doppler(isce3::core::LUT2d refDoppler, + isce3::core::LUT2d secDoppler) { _refDoppler = refDoppler; _secDoppler = secDoppler; @@ -214,7 +211,8 @@ gpuCrossmul::crossmul(isce3::io::Raster& refSlcRaster, isce3::io::Raster& secSlcRaster, isce3::io::Raster& ifgRaster, isce3::io::Raster& coherenceRaster, - isce3::io::Raster* rngOffsetRaster) const + isce3::io::Raster* rngOffsetRaster, + isce3::io::Raster* aziOffsetRaster) const { // set flatten flag based range offset raster ptr value bool flatten = rngOffsetRaster ? true : false; diff --git a/cxx/isce3/cuda/signal/gpuCrossMul.h b/cxx/isce3/cuda/signal/gpuCrossMul.h index 4e96feaff..199d797c1 100644 --- a/cxx/isce3/cuda/signal/gpuCrossMul.h +++ b/cxx/isce3/cuda/signal/gpuCrossMul.h @@ -2,9 +2,11 @@ #include "forward.h" #include // Raster +#include #include -#include +#include + #include namespace isce3::cuda::signal { @@ -24,28 +26,32 @@ class gpuCrossmul { * \param[out] coherenceRaster output coherence raster * \param[in] rngOffsetRaster optional pointer to range offset raster * if provided, interferogram will be flattened + * \param[in] aziOffsetRaster optional pointer to azimuth offset raster + * it is a placeholder here, and will implement the + * azimuth common band filter based on it */ void crossmul(isce3::io::Raster& refSlcRaster, isce3::io::Raster& secSlcRaster, isce3::io::Raster& ifgRaster, isce3::io::Raster& coherenceRaster, - isce3::io::Raster* rngOffsetRaster = nullptr) const; + isce3::io::Raster* rngOffsetRaster = nullptr, + isce3::io::Raster* aziOffsetRaster = nullptr) const; /** Set doppler LUTs for reference and secondary SLCs*/ - void doppler(isce3::core::LUT1d refDoppler, - isce3::core::LUT1d secDoppler); + void doppler(isce3::core::LUT2d refDoppler, + isce3::core::LUT2d secDoppler); /** Set reference doppler */ - inline void refDoppler(isce3::core::LUT1d refDopp) {_refDoppler = refDopp;}; + inline void refDoppler(isce3::core::LUT2d refDopp) {_refDoppler = refDopp;}; /** Get reference doppler */ - inline const isce3::core::LUT1d & refDoppler() const {return _refDoppler;}; + inline const isce3::core::LUT2d & refDoppler() const {return _refDoppler;}; /** Set secondary doppler */ - inline void secDoppler(isce3::core::LUT1d secDopp) {_secDoppler = secDopp;}; + inline void secDoppler(isce3::core::LUT2d secDopp) {_secDoppler = secDopp;}; /** Get secondary doppler */ - inline const isce3::core::LUT1d & secDoppler() const {return _secDoppler;}; + inline const isce3::core::LUT2d & secDoppler() const {return _secDoppler;}; /** Set reference and secondary starting range shift */ inline void startingRangeShift(double rng_shift) { _offsetStartingRangeShift = rng_shift; } @@ -93,11 +99,11 @@ class gpuCrossmul { inline bool multiLookEnabled() const { return _multiLookEnabled; } private: - //Doppler LUT for the reference SLC - isce3::core::LUT1d _refDoppler; + //Doppler LUT for the refernce SLC + isce3::core::LUT2d _refDoppler; //Doppler LUT for the secondary SLC - isce3::core::LUT1d _secDoppler; + isce3::core::LUT2d _secDoppler; // starting range shifts between the secondary and reference RSLC in meters double _offsetStartingRangeShift = 0.0; diff --git a/cxx/isce3/geocode/GeocodeCov.cpp b/cxx/isce3/geocode/GeocodeCov.cpp index 09f9e6ce4..25d9566ab 100644 --- a/cxx/isce3/geocode/GeocodeCov.cpp +++ b/cxx/isce3/geocode/GeocodeCov.cpp @@ -491,7 +491,7 @@ void Geocode::geocodeInterp( std::make_unique( vsimem_ref, radar_grid.width(), radar_grid.length(), 1, GDT_Float32, "ENVI"); - rtc_sigma0_raster = + rtc_sigma0_raster = rtc_raster_sigma0_unique_ptr.get(); } else { rtc_sigma0_raster = output_rtc_sigma; @@ -542,7 +542,7 @@ void Geocode::geocodeInterp( block_length = geogrid.length(); } else { - if (geocode_memory_mode == + if (geocode_memory_mode == isce3::core::GeocodeMemoryMode::BlocksGeogridAndRadarGrid) { warning << "WARNING the geocode memory mode" << " BlocksGeogridAndRadarGrid is not available" @@ -554,7 +554,7 @@ void Geocode::geocodeInterp( isce3::core::getBlockProcessingParametersY( geogrid.length(), geogrid.width(), nbands, sizeof(T), &info, &block_length, &nBlocks, min_block_size, max_block_size); - } + } info << "number of blocks: " << nBlocks << pyre::journal::newline; info << "block length: " << block_length << pyre::journal::newline; @@ -1012,7 +1012,6 @@ inline void Geocode::_interpolate( } if (sub_swaths != nullptr) { - // read the sub-swath value at the center to save the output mask uint8_t sample_sub_swath_center = sub_swaths->getSampleSubSwath( rdr_y_rslc, rdr_x_rslc); @@ -1047,10 +1046,10 @@ inline void Geocode::_interpolate( } } - /* + /* check within the interpolation kernel (approximated by `interp_margin`) if any of the samples is marked as shadow or layover-and-shadow - in which case we skip to the next position, i.e., we "break" the + in which case we skip to the next position, i.e., we "break" the 2 inner for-loop bellow (vars: yy and xx) and "continue" from the parent for-loop (var: kk) above. */ @@ -2079,7 +2078,7 @@ void Geocode::geocodeAreaProj( output_rtc_sigma == nullptr) { std::string vsimem_ref = ( "/vsimem/" + getTempString("geocode_cov_areaproj_rtc_sigma0")); - rtc_raster_sigma0_unique_ptr = + rtc_raster_sigma0_unique_ptr = std::make_unique( vsimem_ref, radar_grid_cropped.width(), radar_grid_cropped.length(), 1, GDT_Float32, "ENVI"); @@ -2259,7 +2258,7 @@ void Geocode::geocodeAreaProj( } else { isce3::core::getBlockProcessingParametersXY( imax, jmax, nbands + nbands_off_diag_terms, sizeof(T_out), - &info, &block_size_with_upsampling_y, &nblocks_y, + &info, &block_size_with_upsampling_y, &nblocks_y, &block_size_with_upsampling_x, &nblocks_x, min_block_size, max_block_size, geogrid_upsampling); block_size_x = block_size_with_upsampling_x / geogrid_upsampling; @@ -2509,7 +2508,7 @@ void Geocode::_runBlock( isce3::io::Raster& dem_raster, isce3::io::Raster* out_off_diag_terms, isce3::io::Raster* out_geo_rdr, isce3::io::Raster* out_geo_dem, isce3::io::Raster* out_geo_nlooks, isce3::io::Raster* out_geo_rtc, - isce3::io::Raster* out_geo_rtc_gamma0_to_sigma0, + isce3::io::Raster* out_geo_rtc_gamma0_to_sigma0, isce3::core::ProjectionBase* proj, bool flag_apply_rtc, bool flag_rtc_raster_is_in_memory, bool flag_rtc_sigma0_raster_is_in_memory, isce3::io::Raster* rtc_raster, isce3::io::Raster* rtc_sigma0_raster, @@ -2905,7 +2904,7 @@ void Geocode::_runBlock( this_block_size_y, this_block_size_x)); nan_t_out *= std::numeric_limits::quiet_NaN(); - + for (int band = 0; band < nbands; ++band) geoDataBlock[band]->fill(nan_t_out); @@ -3286,7 +3285,7 @@ void Geocode::_runBlock( also need to add `offset_x` and `offset_y` that represent the offsets in X- and Y- directions over the radar-grid coordinates. - in which case we skip to the next position, i.e., we "break" the + in which case we skip to the next position, i.e., we "break" the 2 inner for-loop bellow (vars: yy and xx) and "continue" from the parent for-loop (var: kk) above. */ @@ -3647,7 +3646,7 @@ std::string _get_geocode_memory_mode_str( } template -void Geocode::_print_parameters(pyre::journal::info_t& channel, +void Geocode::_print_parameters(pyre::journal::info_t& channel, isce3::core::GeocodeMemoryMode& geocode_memory_mode, const long long min_block_size, const long long max_block_size) { diff --git a/cxx/isce3/math/Sinc.h b/cxx/isce3/math/Sinc.h index a1ef33fd3..ee4161fe9 100644 --- a/cxx/isce3/math/Sinc.h +++ b/cxx/isce3/math/Sinc.h @@ -1,10 +1,15 @@ #pragma once #include +#include namespace isce3 { namespace math { /** sinc function defined as \f$ \frac{\sin(\pi x)}{\pi x} \f$ */ +template +CUDA_HOSTDEV +std::complex sinc(std::complex t); + template CUDA_HOSTDEV T sinc(T t); diff --git a/cxx/isce3/math/Sinc.icc b/cxx/isce3/math/Sinc.icc index c2a4e67a2..1864c3a5e 100644 --- a/cxx/isce3/math/Sinc.icc +++ b/cxx/isce3/math/Sinc.icc @@ -90,4 +90,29 @@ T sinc(T t) return std::sin(x) / x; } +template +CUDA_HOSTDEV +inline +std::complex sinc(std::complex t) +{ + static_assert(std::is_floating_point::value); + + using U = typename std::remove_cv::type; + static constexpr auto eps1 = detail::root_epsilon(); + static constexpr auto eps2 = detail::fourth_root_epsilon(); + + auto u = T(M_PI) * t; + auto x = std::abs(u); + + // approximate by Taylor series expansion near zero + if (x < eps2) { + std::complex out = static_cast(1); + if (x > eps1) { + out -= u * u / static_cast(6); + } + return out; + } + return std::sin(u) / u; +} + }} diff --git a/cxx/isce3/product/RadarGridProduct.h b/cxx/isce3/product/RadarGridProduct.h index 1802b6bda..860cc9597 100644 --- a/cxx/isce3/product/RadarGridProduct.h +++ b/cxx/isce3/product/RadarGridProduct.h @@ -27,7 +27,7 @@ namespace isce3 { namespace product { /** * Return the path to each child group of `group` that ends with the substring * `group_name`. - * + * * \param[in] group Parent group * \param[in] group_name Search string * \returns List of child group paths @@ -39,13 +39,13 @@ std::vector findGroupPath( * Return grids or swaths group paths within the base_group. * Start by assigning an empty string to image_group_str in case * grids and swaths group are not found. - * + * * \param[in] file * \param[in] base_dir Path to `base_group` object (e.g. '/science/') * \param[in] base_group Base group * \param[in] key_vector Vector containing possible image groups * (e.g., 'swaths', 'grids', or both) to look for - * \param[out] image_group_str Path to first image group found containing + * \param[out] image_group_str Path to first image group found containing * one of the `key_vector` keys (e.g., '/science/LSAR/RSLC/swaths') * \param[in] metadata_group_str Path to first metadata group found by * substituting `key` with `metadata` in `image_group_str` @@ -60,10 +60,10 @@ void setImageMetadataGroupStr( std::string &metadata_group_str); /** RadarGridProduct class declaration - * + * * The Produt attribute Swaths map, i.e. _swaths, associates the - * frequency (key) with the Swath object (value). The RadarGridProduct object - * is usually initiated with an empty map and the serialization of + * frequency (key) with the Swath object (value). The RadarGridProduct object + * is usually initiated with an empty map and the serialization of * the SAR product is responsible for populating the Swath map * from the product's metadata. * @@ -75,7 +75,7 @@ class RadarGridProduct { RadarGridProduct(isce3::io::IH5File &); /** Constructor with Metadata and Swath map. */ inline RadarGridProduct(const Metadata &, const std::map &); - + /** Get a read-only reference to the metadata */ inline const Metadata & metadata() const { return _metadata; } /** Get a reference to the metadata. */ diff --git a/cxx/isce3/product/Swath.h b/cxx/isce3/product/Swath.h index 180dc825b..58780e288 100644 --- a/cxx/isce3/product/Swath.h +++ b/cxx/isce3/product/Swath.h @@ -118,7 +118,7 @@ class isce3::product::Swath { inline isce3::product::SubSwaths& subSwaths() { return _subSwaths; } - + inline const isce3::product::SubSwaths& subSwaths() const { return _subSwaths; } diff --git a/cxx/isce3/signal/Crossmul.cpp b/cxx/isce3/signal/Crossmul.cpp index d4a0ef15b..68a1acbf9 100644 --- a/cxx/isce3/signal/Crossmul.cpp +++ b/cxx/isce3/signal/Crossmul.cpp @@ -4,6 +4,9 @@ #include "Looks.h" #include "Signal.h" +#include +#include + /** * Compute the frequency response due to a subpixel shift introduced by * upsampling and downsampling @@ -32,7 +35,7 @@ void lookdownShiftImpact(size_t oversample, size_t fft_size, size_t blockRows, // As an example for a signal with length of 5 and : // original sample locations: 0 1 2 3 4 // upsampled sample locations: 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 - // Looked dow sample locations: 0.25 1.25 2.25 3.25 4.25 + // Looked down sample locations: 0.25 1.25 2.25 3.25 4.25 // Obviously the signal after looking down would be shifted by 0.25 pixel in // range comared to the original signal. Since a shift in time domain introduces // a linear phase in frequency domain, we compute the impact in frequency domain. @@ -65,16 +68,314 @@ size_t omp_thread_count() { return n; } +/** +* @param[in, out] refSlc a block of the reference SLC to be filtered +* @param[in, out] secSlc a block of second SLC to be filtered +* @param[in] geometryIfgram a simulated interferogram that contains the geometrical phase due to baseline separation +* @param[in] geometryIfgramConj conjugate of geometryIfgram +* @param[in, out] refSpectrum spectrum of geometryIfgramConj in range direction +* @param[in, out] secSpectrum spectrum of geometryIfgram in range direction +* @param[in] rangeFrequencies frequencies in range direction +* @param[in] rngFilter a filter object +* @param[in] blockLength number of rows +* @param[in] ncols number of columns +* @param[in] maxRangeFilterKernelSize maximum range filter kernel size +* @returns the common bandwidth +*/ +double isce3::signal::Crossmul:: +rangeCommonBandFilter(std::valarray> &refSlc, + std::valarray> &secSlc, + const std::valarray> &geometryIfgram, + const std::valarray> &geometryIfgramConj, + std::valarray> &refSpectrum, + std::valarray> &secSpectrum, + std::valarray &rangeFrequencies, + isce3::signal::Filter &rngFilter, + size_t blockLength, + size_t ncols, + const size_t maxRangeFilterKernelSize) +{ + pyre::journal::debug_t debug("isce.signal.Crossmul.rangeCommonBandFilter"); + + // size of the arrays + size_t vectorLength = refSlc.size(); + + // Aligning the spectrum of the two SLCs + // Shifting the range spectrum of each image according to the local + // (slope-dependent) wavenumber. This shift in frequency domain is + // achieved by removing/adding the geometrical phase (representing topography) + // from/to reference and secondary SLCs in time domain. + #pragma omp parallel for + for (size_t i = 0; i < refSlc.size(); i++) { + refSlc[i] *= geometryIfgramConj[i]; + secSlc[i] *= geometryIfgram[i]; + } + + + // determine the frequency shift, in Hz based on the power spectral density of + // the geometrical interferometric phase using an empirical approach + double frequencyShift = computeRangeFrequencyShift(refSpectrum, + secSpectrum, + rangeFrequencies, + blockLength, + ncols); + + debug << "rangeFrequencyShift (MHz): "<< frequencyShift/1e6 << pyre::journal::endl; + debug << "range bandwidth (MHz): " << _rangeBandwidth/1e6 << pyre::journal::endl; + + // Since the spectrum of the ref and sec SLCs are already aligned, + // we design the low-pass filter as a band-pass at zero frequency with + // bandwidth of (range bandwidth - frequency shift) + const double filterCenterFrequency = 0.0; + const double filterBandwidth = _rangeBandwidth - fabs(frequencyShift); + + if (_windowType == "kaiser") { + auto [n, _] = rngFilter._kaiserord(_ripple, + _transitionWidth * filterBandwidth/_rangeSamplingFrequency); + if (n > maxRangeFilterKernelSize) { + pyre::journal::error_t err( + "isce.signal.Crossmul.rangeCommonBandFilter"); + err << "Max filter length exceeded due to insufficient " + << "spectral overlap (perpendicular baseline is too large)" + << pyre::journal::endl; + throw isce3::except::LengthError(ISCE_SRCINFO(), + "Max filter length exceeded due to insufficient spectral overlap (perpendicular baseline is too large)"); + } + } + + // Contruct the low pass filter for this block. This filter is + // common for both SLCs + rngFilter.constructRangeCommonBandFilter(_rangeSamplingFrequency, + filterCenterFrequency, + filterBandwidth, + ncols, + blockLength, + _windowType, + _windowParameter, + maxRangeFilterKernelSize); + + // low pass filter the ref slc + rngFilter.initiateRangeFilter(refSlc,refSpectrum, ncols, blockLength); + rngFilter.filter(refSlc, refSpectrum); + + // low pass filter the sec slc + rngFilter.initiateRangeFilter(secSlc,secSpectrum, ncols, blockLength); + rngFilter.filter(secSlc, secSpectrum); + + // restore the original phase without the geometry phase + // in case other steps will use the original phase + #pragma omp parallel for + for (size_t i = 0; i < refSlc.size(); i++) { + refSlc[i] *= geometryIfgram[i]; + secSlc[i] *= geometryIfgramConj[i]; + } + + return filterBandwidth; +} + +/** +* @param[in, out] refSlc a block of the reference SLC to be filtered +* @param[in, out] secSlc a block of second SLC to be filtered +* @param[in] refDoppCentroids reference doppler centroid +* @param[in] secDoppCentroids secondary doppler centroid +* @param[in, out] refAzimuthSpectrum spectrum of the reference after filtering +* @param[in, out] secAzimuthSpectrum spectrum of the secondary after filtering +* @param[in] azimuthFilter a filter object +* @param[in] blockLength number of rows +* @param[in] ncols number of columns +* @returns the common azimuth bandwidth +*/ +double isce3::signal::Crossmul:: +azimuthCommonBandFilter(std::valarray> &refSlc, + std::valarray> &secSlc, + const std::valarray &refDoppCentroids, + const std::valarray &secDoppCentroids, + std::valarray> &refAzimuthSpectrum, + std::valarray> &secAzimuthSpectrum, + isce3::signal::Filter &azimuthFilter, + size_t blockRows, + size_t ncols) +{ + // Construct azimuth common bandpass filter for both reference and secondary + double processedAzimuthBandwidth = azimuthFilter.constructAzimuthCommonBandFilter( + refDoppCentroids, secDoppCentroids, + _azimuthBandwidth, + _prf, _windowParameter, ncols, + blockRows, _windowType); + + // Filter a block of data of the reference + azimuthFilter.initiateAzimuthFilter(refSlc, refAzimuthSpectrum, ncols, blockRows); + azimuthFilter.filter(refSlc, refAzimuthSpectrum); + + // Filter a block of data of the secondary + azimuthFilter.initiateAzimuthFilter(secSlc, secAzimuthSpectrum, ncols, blockRows); + azimuthFilter.filter(secSlc, secAzimuthSpectrum); + + return processedAzimuthBandwidth; +} + +/** +* @param[in] refDopplerCentroids doppler centroid frequency for reference, w.r.t slant range +* @param[in] secDopplerCentroids doppler centroid frequency for secondary, w.r.t slant range +* @param[in] numOfDopplerCentroids number of valid doppler centroids +* @param[in] bandwidth intput SLCs bandwidth +* @param[in] prf pulse repeat frequency +* @param[in] beta kaiser window parameter +* @param[in] aziFilter azimuth filter +* @returns the maximum azimuth filter kernel size +*/ +int +isce3::signal::Crossmul::_computeMaxAzimuthFilterKernelSize(const std::valarray &refDopplerCentroids, + const std::valarray &secDopplerCentroids, + const std::valarray &numOfValidDopplerCentroids, + const double bandwidth, + const double prf, + isce3::signal::Filter &aziFilter) +{ + int max_kernel_size = 0; + + // Loop over columns to get the maximum kernel size + #pragma omp parallel for reduction(max:max_kernel_size) + for (size_t j = 0; j < refDopplerCentroids.size(); ++j) { + if (numOfValidDopplerCentroids[j] > 0) { + double refFreq = refDopplerCentroids[j]; + double secFreq = secDopplerCentroids[j]; + double fshift = std::abs(refFreq - secFreq); + + // The bandwidth should be less than frequency shift + if (bandwidth < fshift) { + std::string err_msg = "Bandwidth " + + std::to_string(bandwidth) + + " is less than frequency shift " + + std::to_string(fshift); + pyre::journal::error_t err( + "isce.signal.Crossmul._computeMaxAzimuthFilterKernelSize"); + err << err_msg << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), err_msg); + } + + // Normalized bandwdith + const double bw = (bandwidth - fshift)/ prf; + + // Transition width is specified in terms of output bandwidth, so scale to + // get width at sample rate of filter. + const double tw = _transitionWidth * bw; + if ((bw + tw / 2.0) > 1.0) { + pyre::journal::error_t err( + "isce.signal.Crossmul._maximum_kernel_size"); + err << "Passband + transition cannot exceed Nyquist" + << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Passband + transition cannot exceed Nyquist"); + } + + auto [n, _] = aziFilter._kaiserord(_ripple, tw); + if (n > max_kernel_size) max_kernel_size = n; + } + } + + return max_kernel_size; +} + +/** +* @param[in] refDoppler 2d LUT doppler centroid frequency for reference +* @param[in] secDoppler 2d LUT doppler centroid frequency for secondary +* @param[in] rngOffsetRaster range offset product pointer +* @param[out] refDopplerCentroids azimuth mean reference doppler centroids +* @param[out] secDopplerCentroids azimuth mean secondary doppler centroids +* @param[out] numOfValidDopplerCentroids number of valid doppler centroids + +*/ +void +isce3::signal::Crossmul::_computeDoppCentroids(const isce3::core::LUT2d & refDoppler, + const isce3::core::LUT2d & secDoppler, + isce3::io::Raster* rngOffsetRaster, + isce3::io::Raster* aziOffsetRaster, + std::valarray &refDopplerCentroids, + std::valarray &secDopplerCentroids, + std::valarray &numOfValidDopplerCentroids) +{ + const size_t ncols = rngOffsetRaster->width(); + const size_t nrows = rngOffsetRaster->length(); + + if (refDopplerCentroids.size() <= 0) { + refDopplerCentroids.resize(ncols); + refDopplerCentroids = 0.0; + } + if (secDopplerCentroids.size() <= 0) { + secDopplerCentroids.resize(ncols); + secDopplerCentroids = 0.0; + } + if (numOfValidDopplerCentroids.size() <= 0) { + numOfValidDopplerCentroids.resize(ncols); + numOfValidDopplerCentroids = 0; + } + + std::valarray rangeOffsets(ncols); + std::valarray azimuthOffsets(ncols); + + // Compute the doppler centroids for the reference and secondary images + for (size_t row = 0; row < nrows; row++) { + + rngOffsetRaster->getLine(rangeOffsets, row); + aziOffsetRaster->getLine(azimuthOffsets, row); + + #pragma omp parallel for + for (size_t col = 0; col < ncols; col++) { + // Convert the line/pixel to range doppler coordinates + double refX = col * rangePixelSpacing() + refStartRange(); + double refY = row / prf() + refStartAzimuthTime(); + + double secX = (col + rangeOffsets[col]) * rangePixelSpacing() + secStartRange(); + double secY = (row + azimuthOffsets[col]) / prf() + secStartAzimuthTime(); + + // Interpolate the doppler centroids + if (refDoppler.contains(refY, refX) && secDoppler.contains(secY, secX)) { + refDopplerCentroids[col] += refDoppler.eval(refY, refX); + secDopplerCentroids[col] += secDoppler.eval(secY, secX); + numOfValidDopplerCentroids[col]++; + } + } + } + + // Using the avergae doppler centroid for each column to avoid the artifacts when + // process block by block + #pragma omp parallel for + for (size_t col = 0; col < ncols; col++) { + if (numOfValidDopplerCentroids[col] > 0) { + refDopplerCentroids[col] /= numOfValidDopplerCentroids[col]; + secDopplerCentroids[col] /= numOfValidDopplerCentroids[col]; + } + } +} + void isce3::signal::Crossmul:: crossmul(isce3::io::Raster& refSlcRaster, isce3::io::Raster& secSlcRaster, isce3::io::Raster& ifgRaster, isce3::io::Raster& coherenceRaster, - isce3::io::Raster* rngOffsetRaster) const + isce3::io::Raster* rngOffsetRaster, + isce3::io::Raster* aziOffsetRaster) { + pyre::journal::info_t info("isce.signal.Crossmul.crossmul"); + pyre::journal::debug_t debug("isce.signal.Crossmul.crossmul"); + + // number of threads + size_t nthreads = omp_thread_count(); + // setting local lines per block to avoid modifying class member size_t linesPerBlock = _linesPerBlock; + //filter objects which will be used for azimuth and range common band filtering + isce3::signal::Filter azimuthFilter; + isce3::signal::Filter rangeFilter; + + //signal object for refSlc + isce3::signal::Signal refSignal(nthreads); + + //signal object for secSlc + isce3::signal::Signal secSignal(nthreads); + // check consistency of input/output raster shapes size_t nrows = refSlcRaster.length(); size_t ncols = refSlcRaster.width(); @@ -87,12 +388,87 @@ crossmul(isce3::io::Raster& refSlcRaster, throw isce3::except::LengthError(ISCE_SRCINFO(), "interferogram and coherence rasters width do not match"); + // maximum range filter kernel size + size_t maxRangeFilterKernelSize = 0; + + // Compute the range sampling frequency in Hz using the range pixel spacing + _rangeSamplingFrequency = 1.0 / (_rangePixelSpacing*2.0/isce3::core::speed_of_light); + if (_doCommonRangeBandFilter) { + // Determine the max range filter kernel size + auto [n, _] = rangeFilter._kaiserord(_ripple, + _minRangeSpectrumOverlapFraction * + rangeBandwidth()/_rangeSamplingFrequency * _transitionWidth); + maxRangeFilterKernelSize = n; + debug << "max range filter kernel size: " << maxRangeFilterKernelSize << pyre::journal::endl; + } + + // Compute FFT size (power of 2) and set up the maximum range window size + size_t fft_size = isce3::fft::nextFastPower(ncols + maxRangeFilterKernelSize); + + if (fft_size > INT_MAX) + throw isce3::except::LengthError(ISCE_SRCINFO(), "fft_size > INT_MAX"); + if (_oversampleFactor * fft_size > INT_MAX) + throw isce3::except::LengthError(ISCE_SRCINFO(), "_oversampleFactor * fft_size > INT_MAX"); + + // force the multilook to be true if azimuth or range looks > 1 + _multiLookEnabled = ((_rangeLooks > 1) || (_azimuthLooks > 1)); + + // Declare valarray for range and azimuth doppler centroids used by filter + std::valarray refDoppCentroids; + std::valarray secDoppCentroids; + std::valarray numOfValidDoppCentroids; + int maxAzimuthFilterKernelSize = 0; + int overlapSize = 0; + int halfOverlapSize = 0; + + if (_doCommonAzimuthBandFilter) { + // Compute the mean doppler centroid of each slant range for + // the reference and secondary images for each slant range + refDoppCentroids.resize(fft_size); refDoppCentroids = 0.0; + secDoppCentroids.resize(fft_size); secDoppCentroids = 0.0; + numOfValidDoppCentroids.resize(fft_size); numOfValidDoppCentroids = 0; + + info << "determine the max azimuth kernel size" << pyre::journal::newline; + _computeDoppCentroids(_refDoppler, _secDoppler, + rngOffsetRaster, + aziOffsetRaster, + refDoppCentroids, + secDoppCentroids, + numOfValidDoppCentroids); + + maxAzimuthFilterKernelSize = _computeMaxAzimuthFilterKernelSize(refDoppCentroids, + secDoppCentroids, + numOfValidDoppCentroids, + _azimuthBandwidth, + _prf, + azimuthFilter); + + debug << "max azimuth window kernel size: " << maxAzimuthFilterKernelSize + << pyre::journal::endl; + // to ensure the overlap size is an integer multiple number of azimuth looks. + int k = (maxAzimuthFilterKernelSize + _azimuthLooks - 1) / _azimuthLooks; + + // Force the kernel size to be even + k = (k % 2 == 0) ? k : k + 1; + + // The overlap size will be even + overlapSize = k * _azimuthLooks; + halfOverlapSize = overlapSize / 2; + + // Compute the lines per block to account for the overlaps between two blocks + linesPerBlock = ((_linesPerBlock + _azimuthLooks - 1)/ _azimuthLooks) * _azimuthLooks; + _linesPerBlock = linesPerBlock + overlapSize; + + // The lines per block will be multiple times of the azimuth looks + linesPerBlock = _linesPerBlock; + } + const auto output_rows = ifgRaster.length(); const auto output_cols = ifgRaster.width(); if (_multiLookEnabled) { // Making sure that the number of rows in each block (linesPerBlock) // to be an integer multiple of the number of azimuth looks. - linesPerBlock = (_linesPerBlock / _azimuthLooks) * _azimuthLooks; + linesPerBlock = ((_linesPerBlock + _azimuthLooks - 1)/ _azimuthLooks) * _azimuthLooks; // checking only multilook interferogram shape is sufficient // interferogram and coherence shapes checked to match above @@ -115,16 +491,18 @@ crossmul(isce3::io::Raster& refSlcRaster, "full resolution input/output raster widths do not match"); } - size_t nthreads = omp_thread_count(); - // Set flatten flag based range offset raster ptr value - bool flatten = rngOffsetRaster ? true : false; + //signal object for geometryIfgram + isce3::signal::Signal geometryIfgramSignal(nthreads); - //signal object for refSlc - isce3::signal::Signal refSignal(nthreads); + //signal object for geometryIfgramConj + isce3::signal::Signal geometryIfgramConjSignal(nthreads); - //signal object for secSlc - isce3::signal::Signal secSignal(nthreads); + //signal object for refSlc windowing effects revert + isce3::signal::Signal refWindowSignal(nthreads); + + //signal object for secSlc windowing effects revert + isce3::signal::Signal secWindowSignal(nthreads); // instantiate Looks used for multi-looking the interferogram isce3::signal::Looks looksObj; @@ -138,23 +516,20 @@ crossmul(isce3::io::Raster& refSlcRaster, looksObj.nrowsLooked(linesPerBlockMLooked); looksObj.ncolsLooked(ncolsMultiLooked); - // Compute FFT size (power of 2) - size_t fft_size; - refSignal.nextPowerOfTwo(ncols, fft_size); - - if (fft_size > INT_MAX) - throw isce3::except::LengthError(ISCE_SRCINFO(), "fft_size > INT_MAX"); - if (_oversampleFactor * fft_size > INT_MAX) - throw isce3::except::LengthError(ISCE_SRCINFO(), "_oversampleFactor * fft_size > INT_MAX"); - // number of blocks to process - size_t nblocks = nrows / linesPerBlock; + size_t nblocks = nrows / (linesPerBlock - overlapSize); if (nblocks == 0) { nblocks = 1; - } else if (nrows % (nblocks * linesPerBlock) != 0) { + } else if (nrows % (nblocks * (linesPerBlock - overlapSize)) != 0) { nblocks += 1; } + // set up the processed azimuth and range bandwidth + // if common band filters are enabled, they will be the mean bandwidth of all data block + // otherwise, they are the original SLC bandwidths + _processedAzimuthBandwidth = _doCommonAzimuthBandFilter ? 0.0 : _azimuthBandwidth; + _processedRangeBandwidth = _doCommonRangeBandFilter ? 0.0 : _rangeBandwidth; + // size of not-unsampled valarray const auto spectrumSize = fft_size * linesPerBlock; @@ -168,12 +543,14 @@ crossmul(isce3::io::Raster& refSlcRaster, std::valarray> secSlc(spectrumSize); // storage for a block of range offsets - std::valarray rngOffset(ncols*linesPerBlock); + std::valarray rngOffset(spectrumSize); + + // InSAR phase due to topography + std::valarray> geometryIfgram; - // storage for a simulated interferogram which its phase is the - // interferometric phase due to the imaging geometry: - // phase = (4*PI/wavelength)*(rangePixelSpacing)*(rngOffset) // complex conjugate of geometryIfgram + // both the flattening and range common band filtering will use + // this variable, so creating a buffer here std::valarray> geometryIfgramConj(spectrumSize); // upsampled interferogram @@ -224,10 +601,13 @@ crossmul(isce3::io::Raster& refSlcRaster, // upsampled block of secondary SLC std::valarray> secSlcUpsampled; - // only resize valarrays and init FFT when oversampling - if (_oversampleFactor > 1) { + if ((_oversampleFactor > 1) || _doCommonRangeBandFilter) { refSpectrum.resize(spectrumSize); secSpectrum.resize(spectrumSize); + } + + // only resize valarrays and init FFT when oversampling + if (_oversampleFactor > 1) { refSpectrumUpsampled.resize(spectrumUpsampleSize); secSpectrumUpsampled.resize(spectrumUpsampleSize); @@ -252,13 +632,73 @@ crossmul(isce3::io::Raster& refSlcRaster, lookdownShiftImpact(_oversampleFactor, fft_size, linesPerBlock, shiftImpact); - // loop over all blocks - std::cout << "nblocks : " << nblocks << std::endl; + // Declare valarray for range frequencies used by filter + std::valarray rangeFrequencies; + + // Declare valarray for azimuth spectrum used by filter + std::valarray> refAzimuthSpectrum; + std::valarray> secAzimuthSpectrum; + + // Retrieve the original filter to revert the range windowing effects + std::valarray> originalRangeFilter; + + if (_doCommonAzimuthBandFilter) { + // Allocate storage for azimuth spectrum + refAzimuthSpectrum.resize(spectrumSize); + secAzimuthSpectrum.resize(spectrumSize); + } + + if (_doCommonRangeBandFilter) { + geometryIfgram.resize(spectrumSize); + rangeFrequencies.resize(fft_size); + + // Compute the range frequency for each pixel + fftfreq(1.0/_rangeSamplingFrequency, rangeFrequencies); + + // For NISAR, use a standard Kaiser window in frequency domain to + // compensate the windowing effects in range direction + // NOTE: This function is sensor dependent + if (_sensorType == "NISAR") { + originalRangeFilter.resize(fft_size); + std::valarray subBandCenterFrequencies{0.0}; + std::valarray subBandBandwidths{_rangeSamplingFrequency}; + rangeFilter.constructRangeBandpassKaiser(subBandCenterFrequencies, + subBandBandwidths, + 1.0/_rangeSamplingFrequency, + fft_size, rangeFrequencies, 1.6, + originalRangeFilter); + } + + // Construct the FFTW plan for both geometryIfgram and its conjugation + geometryIfgramSignal.forwardRangeFFT(geometryIfgram, refSpectrum, + fft_size, linesPerBlock); + geometryIfgramSignal.inverseRangeFFT(refSpectrum, geometryIfgram, + fft_size, linesPerBlock); + + geometryIfgramConjSignal.forwardRangeFFT(geometryIfgramConj, secSpectrum, + fft_size, linesPerBlock); + geometryIfgramConjSignal.inverseRangeFFT(secSpectrum, geometryIfgramConj, + fft_size, linesPerBlock); + + refWindowSignal.forwardRangeFFT(refSlc, refSpectrum, + fft_size, linesPerBlock); + refWindowSignal.inverseRangeFFT(refSpectrum, refSlc, + fft_size, linesPerBlock); + + secWindowSignal.forwardRangeFFT(secSlc, secSpectrum, + fft_size, linesPerBlock); + secWindowSignal.inverseRangeFFT(secSpectrum, secSlc, + fft_size, linesPerBlock); + } + // loop over all blocks + std::cout << "nblocks: " << nblocks << std::endl; + info << "nblocks: " << nblocks << pyre::journal::newline; for (size_t block = 0; block < nblocks; ++block) { std::cout << "block: " << block << std::endl; + info << "block: " << block << pyre::journal::newline; // start row for this block - const auto rowStart = block * linesPerBlock; + const auto rowStart = block * (linesPerBlock - overlapSize); //number of lines of data in this block. blockRowsData<= linesPerBlock //Note that linesPerBlock is fixed number of lines @@ -272,6 +712,9 @@ crossmul(isce3::io::Raster& refSlcRaster, secSlc = 0; ifgramUpsampled = 0; ifgram = 0; + geometryIfgramConj = 0; + geometryIfgram = 0; + rngOffset = 0; // get a block of reference and secondary SLC data // and a block of range offsets @@ -286,6 +729,71 @@ crossmul(isce3::io::Raster& refSlcRaster, secSlc[std::slice(line*fft_size, ncols, 1)] = dataLine; } + // load the range offsets that are required by flattening, + // common range filters + if (_doFlatten || _doCommonRangeBandFilter) { + std::valarray offsetLine(ncols); + for (size_t line = 0; line < blockRowsData; ++line) { + rngOffsetRaster->getLine(offsetLine, rowStart + line); + rngOffset[std::slice(line*fft_size, ncols, 1)] = offsetLine; + } + } + + // common range band-pass filtering + if (_doCommonRangeBandFilter) { + // Some diagnostic messages to make sure everything has been configured + debug << "range pixel spacing (m): " << _rangePixelSpacing << pyre::journal::endl; + debug << "wavelength (m): " << _wavelength << pyre::journal::endl; + + refWindowSignal.forward(refSlc, refSpectrum); + secWindowSignal.forward(secSlc, secSpectrum); + + // Convert range offset from meters to complex one-way geometric phase + #pragma omp parallel for + for (size_t line = 0; line < blockRowsData; ++line) { + for (size_t col = 0; col < ncols; ++col) { + + // one-way geometric phase to shift the spectrum + double phase = 2.0 * M_PI + * _rangePixelSpacing*rngOffset[line*fft_size+col] + / _wavelength; + + geometryIfgram[line*fft_size + col] = + std::complex (std::cos(phase), std::sin(phase)); + geometryIfgramConj[line*fft_size + col] = std::conj(geometryIfgram[line*fft_size + col]); + + // revert the original windowing effects along the slant range direction + if (originalRangeFilter.size() > 0) { + refSpectrum[line*fft_size + col] /= originalRangeFilter[col]; + secSpectrum[line*fft_size + col] /= originalRangeFilter[col]; + } + } + } + + refWindowSignal.inverse(refSpectrum, refSlc); + secWindowSignal.inverse(secSpectrum, secSlc); + + // Forward FFT to compute geometry-dependent spectrum to determine the + // frequency shift + geometryIfgramConjSignal.forward(geometryIfgramConj, refSpectrum); + geometryIfgramSignal.forward(geometryIfgram, secSpectrum); + + // Apply range common band filter to ref and sec SLC + // and the resultant refSlc and secSlc have no topo phase removal + _processedRangeBandwidth += rangeCommonBandFilter(refSlc, secSlc, geometryIfgram, + geometryIfgramConj, refSpectrum, secSpectrum, + rangeFrequencies, rangeFilter, linesPerBlock, fft_size, + maxRangeFilterKernelSize); + } + + // Apply the azimuth common band-pass filter to the reference and secondary SLCs + if (_doCommonAzimuthBandFilter) { + _processedAzimuthBandwidth += azimuthCommonBandFilter(refSlc, secSlc, + refDoppCentroids, secDoppCentroids, + refAzimuthSpectrum, secAzimuthSpectrum, + azimuthFilter, linesPerBlock, fft_size); + } + // upsample the reference and secondary SLCs if (_oversampleFactor == 1) { refSlcUpsampled = refSlc; @@ -307,7 +815,7 @@ crossmul(isce3::io::Raster& refSlcRaster, } } - if (flatten) { + if (_doFlatten) { // Read range offsets std::valarray offsetLine(ncols); for (size_t line = 0; line < blockRowsData; ++line) { @@ -318,10 +826,9 @@ crossmul(isce3::io::Raster& refSlcRaster, #pragma omp parallel for for (size_t line = 0; line < blockRowsData; ++line) { for (size_t col = 0; col < ncols; ++col) { - double phase = 4.0*M_PI*_rangePixelSpacing*rngOffset[line*ncols+col]/_wavelength; + double phase = 4.0*M_PI*_rangePixelSpacing*rngOffset[line*fft_size+col]/_wavelength; geometryIfgramConj[line*fft_size + col] = std::complex (std::cos(phase), -1.0*std::sin(phase)); - } } } @@ -336,20 +843,42 @@ crossmul(isce3::io::Raster& refSlcRaster, sum += ifgramUpsampled[line*(ncols*_oversampleFactor) + j + col*_oversampleFactor]; ifgram[line*ncols + col] = sum/ov; - if (flatten) + if (_doFlatten) { ifgram[line*ncols + col] *= geometryIfgramConj[line*fft_size + col]; + } } } // Take looks down (summing columns) if (_multiLookEnabled) { - // mulitlook interferogram and set raster looksObj.ncols(ncols); looksObj.colsLooks(_rangeLooks); looksObj.multilook(ifgram, ifgramMultiLooked); - ifgRaster.setBlock(ifgramMultiLooked, 0, rowStart/_azimuthLooks, - ncols/_rangeLooks, blockRowsData/_azimuthLooks); + + size_t rowNewStart = 0; + // if there is only one block, the nrowsMultiLooked will be equal to blockRowsData + size_t nrowsMultiLooked = (nblocks == 1) ? blockRowsData : blockRowsData - halfOverlapSize; + nrowsMultiLooked /= _azimuthLooks; + + // The first block + std::slice dataSlice = std::slice(0, ncolsMultiLooked * nrowsMultiLooked, 1); + if (block > 0) { + rowNewStart = (rowStart + halfOverlapSize)/_azimuthLooks; + if (block != (nblocks - 1)) { + nrowsMultiLooked = (blockRowsData - overlapSize) / _azimuthLooks; + dataSlice = std::slice(halfOverlapSize/_azimuthLooks * ncolsMultiLooked, + ncolsMultiLooked * nrowsMultiLooked, 1); + } else { // The last block + nrowsMultiLooked = (blockRowsData - halfOverlapSize) / _azimuthLooks; + dataSlice = std::slice(halfOverlapSize/_azimuthLooks * ncolsMultiLooked, + ncolsMultiLooked * nrowsMultiLooked, 1); + } + } + // set the block of interferogram + std::valarray> ifgramSlice= ifgramMultiLooked[dataSlice]; + ifgRaster.setBlock(ifgramSlice, 0, rowNewStart, + ncolsMultiLooked, nrowsMultiLooked); // multilook SLC to power for coherence computation // refPowerLooked = average(abs(refSlc)^2) @@ -373,18 +902,41 @@ crossmul(isce3::io::Raster& refSlcRaster, } // set coherence raster - coherenceRaster.setBlock(coherence, 0, rowStart/_azimuthLooks, - ncols/_rangeLooks, blockRowsData/_azimuthLooks); + std::valarray coherenceSlice = coherence[dataSlice]; + coherenceRaster.setBlock(coherenceSlice, 0, rowNewStart, + ncolsMultiLooked,nrowsMultiLooked); } else { + // The first block + size_t rowNewStart = 0; + // if there is only one block, the nrowsValid will be equal to blockRowsData + size_t nrowsValid = (nblocks == 1) ? blockRowsData : blockRowsData - halfOverlapSize; + std::slice dataSlice = std::slice(0, nrowsValid * ncols, 1); + if (block > 0) { + rowNewStart = rowStart + halfOverlapSize; + if (block != (nblocks - 1)) { + nrowsValid = blockRowsData - overlapSize; + dataSlice = std::slice(halfOverlapSize * ncols, nrowsValid, 1); + } else { // the last block + nrowsValid = blockRowsData - halfOverlapSize; + dataSlice = std::slice(halfOverlapSize * ncols,nrowsValid, 1); + } + } + + std::valarray> ifgramSlice= ifgram[dataSlice]; // set the block of interferogram - ifgRaster.setBlock(ifgram, 0, rowStart, ncols, blockRowsData); + ifgRaster.setBlock(ifgramSlice, 0, rowNewStart, ncols, nrowsValid); // fill coherence with ones (no need to compute result) coherence = 1.0; - // set the block of coherence - coherenceRaster.setBlock(coherence, 0, rowStart, ncols, - blockRowsData); + std::valarray coherenceSlice = coherence[dataSlice]; + coherenceRaster.setBlock(coherenceSlice, 0, rowNewStart, ncols, + nrowsValid); } } + + // update the azimuth and range bandwidth after common band filtering + // using the mean bandwidth + if (_doCommonRangeBandFilter) _processedRangeBandwidth /= nblocks; + if (_doCommonAzimuthBandFilter) _processedAzimuthBandwidth /= nblocks; } diff --git a/cxx/isce3/signal/Crossmul.h b/cxx/isce3/signal/Crossmul.h index c7b269888..6ae73b49e 100644 --- a/cxx/isce3/signal/Crossmul.h +++ b/cxx/isce3/signal/Crossmul.h @@ -3,7 +3,9 @@ #include "forward.h" #include -#include +#include +#include +#include #include /** \brief Intereferogram generation by cross-multiplication of reference and secondary SLCs. @@ -24,30 +26,32 @@ class isce3::signal::Crossmul { * \param[in] secSlcRaster input raster of secondary SLC * \param[out] ifgRaster output interferogram raster * \param[out] coherenceRaster output coherence raster - * \param[in] rngOffsetRaster optional pointer to range offset raster - * if provided, interferogram will be flattened + * \param[in] rngOffsetRaster optional pointer to range offset raster, in pixels + * \param[in] aziOffsetRaster optional pointer to azimuth offset raster, in pixels + * */ void crossmul(isce3::io::Raster& refSlcRaster, isce3::io::Raster& secSlcRaster, isce3::io::Raster& ifgRaster, isce3::io::Raster& coherence, - isce3::io::Raster* rngOffsetRaster = nullptr) const; + isce3::io::Raster* rngOffsetRaster = nullptr, + isce3::io::Raster* aziOffsetRaster = nullptr); /** Set doppler LUTs for reference and secondary SLCs*/ - inline void doppler(isce3::core::LUT1d, - isce3::core::LUT1d); + inline void doppler(isce3::core::LUT2d, + isce3::core::LUT2d); /** Set dopplers LUT for reference SLC */ - inline void refDoppler(isce3::core::LUT1d refDopp) { _refDoppler = refDopp; } + inline void refDoppler(isce3::core::LUT2d refDopp) { _refDoppler = refDopp; } /** Get doppler LUT for reference SLC */ - inline const isce3::core::LUT1d & refDoppler() const { return _refDoppler; } + inline const isce3::core::LUT2d & refDoppler() const { return _refDoppler; } /** Set dopplers LUT for secondary SLC */ - inline void secDoppler(isce3::core::LUT1d secDopp) { _secDoppler = secDopp; } + inline void secDoppler(isce3::core::LUT2d secDopp) { _secDoppler = secDopp; } /** Get doppler LUT for secondary SLC */ - inline const isce3::core::LUT1d & secDoppler() const { return _secDoppler; } + inline const isce3::core::LUT2d & secDoppler() const { return _secDoppler; } /** Set reference and seconary starting range shift */ inline void startingRangeShift(double rng_shift) { _offsetStartingRangeShift = rng_shift; } @@ -55,13 +59,29 @@ class isce3::signal::Crossmul { /** Get reference and secondary starting range shift */ inline double startingRangeShift() const { return _offsetStartingRangeShift; } - /** Set range pixel spacing */ + /** Set range pixel spacing, in meters */ inline void rangePixelSpacing(double rgPxlSpacing) { _rangePixelSpacing = rgPxlSpacing; } - /** Get range pixel spacing */ + /** Get range pixel spacing, in meters */ inline double rangePixelSpacing() const { return _rangePixelSpacing; } - /** Set Wavelength*/ + /** Set start range for reference and secondary, in meters */ + inline void refStartRange(double startRng) { _refStartRange = startRng; } + inline void secStartRange(double startRng) { _secStartRange = startRng; } + + /** Get start range for reference and secondary , in meters */ + inline double refStartRange() const { return _refStartRange; } + inline double secStartRange() const { return _secStartRange; } + + /** Set start azimuth time for reference and secondary, in seconds */ + inline void refStartAzimuthTime(double startAziTime) { _refStartAzimuthTime = startAziTime; } + inline void secStartAzimuthTime(double startAziTime) { _secStartAzimuthTime = startAziTime; } + + /** Get start azimuth time for reference and secondary, in seconds */ + inline double refStartAzimuthTime() const { return _refStartAzimuthTime; } + inline double secStartAzimuthTime() const { return _secStartAzimuthTime; } + + /** Set Wavelength, in meters*/ inline void wavelength(double wvl) { _wavelength = wvl; } /** Get Wavelength*/ @@ -79,6 +99,85 @@ class isce3::signal::Crossmul { /** Get number of azimuth looks */ inline int azimuthLooks() const { return _azimuthLooks; } + /** Set common azimuth band filtering flag */ + inline void doCommonAzimuthBandFilter(bool doAzBandFilter) { + _doCommonAzimuthBandFilter = doAzBandFilter; } + + /** Get common azimuth band filtering flag */ + inline bool doCommonAzimuthBandFilter() const { + return _doCommonAzimuthBandFilter; } + + /** Set azimuth bandwidth, in Hz*/ + inline void azimuthBandwidth(double azBandwidth) { + _azimuthBandwidth = azBandwidth; } + + /** Get azimuth bandwidth, in Hz */ + inline double azimuthBandwidth() const { + return _azimuthBandwidth; } + + /** Get processed azimuth bandwidth, in Hz. + * This is the average azimuth bandwidth after + * common-band filtering from among all blocks. + */ + inline double processedAzimuthBandwidth() const { + return _processedAzimuthBandwidth; } + + /** Set window parameter for the common band filter + The meaning of this parameter depends on the `window_type`. + For a raised-cosine window, it is the pedestal height of the window. + For a Kaiser window, it is the beta parameter.*/ + inline void windowParameter(double windowParameter) { _windowParameter = windowParameter; } + + /** Get window parameter for the common band filter */ + inline double windowParameter() const { return _windowParameter; } + + /** Get window type for the common band filter */ + inline std::string windowType() const { return _windowType; } + + /** Set the window type */ + inline void windowType(std::string windowType) { _windowType = windowType; } + + /** Get sensor type*/ + inline std::string sensorType() const { return _sensorType; } + + /** Set the sensor type */ + inline void sensorType(std::string sensorType) { _sensorType = sensorType;} + + /** Set common range band filtering flag */ + inline void doCommonRangeBandFilter(bool doRgBandFilter) { + _doCommonRangeBandFilter = doRgBandFilter; } + + /** Get common range band filtering flag */ + inline bool doCommonRangeBandFilter() const { + return _doCommonRangeBandFilter; } + + /** Set flatten flag */ + inline void doFlatten(bool doFlatten) { + _doFlatten = doFlatten; } + + /** Get flatten flag */ + inline bool doFlatten() const { + return _doFlatten; } + + /** Set pulse repetition frequency (PRF), in Hz */ + inline void prf(double prf) { _prf = prf; } + + /** Get pulse repetition frequency (PRF), in Hz */ + inline double prf() const { return _prf; } + + /** Set the range bandwidth in Hz*/ + inline void rangeBandwidth(double rngBandwidth) { _rangeBandwidth = rngBandwidth; } + + /** Get the range bandwidth, in Hz */ + inline double rangeBandwidth() const {return _rangeBandwidth; } + + /** Get processed range bandwidth after common band filter + * This is the average range bandwidth after + * common-band filtering from among all blocks. + */ + inline double processedRangeBandwidth() const { + return _processedRangeBandwidth; } + /** Set oversample factor */ inline void oversampleFactor(size_t oversamp) { _oversampleFactor = oversamp; } @@ -95,32 +194,86 @@ class isce3::signal::Crossmul { inline bool multiLookEnabled() const { return _multiLookEnabled; } /** Compute the avergae frequency shift in range direction between two SLCs*/ - inline void rangeFrequencyShift(std::valarray> &refAvgSpectrum, + inline double computeRangeFrequencyShift(std::valarray> &refAvgSpectrum, std::valarray> &secAvgSpectrum, std::valarray &rangeFrequencies, size_t linesPerBlockData, - size_t fft_size, - double &frequencyShift); + size_t fft_size); /** estimate the index of the maximum of a vector of data */ inline void getPeakIndex(std::valarray data, size_t &peakIndex); + /** Range common band filtering block by block @insar2007product */ + double rangeCommonBandFilter(std::valarray> &refSlc, + std::valarray> &secSlc, + const std::valarray> &geometryIfgram, + const std::valarray> &geometryIfgramConj, + std::valarray> &refSpectrum, + std::valarray> &secSpectrum, + std::valarray &rangeFrequencies, + isce3::signal::Filter &rngFilter, + size_t blockRows, + size_t ncols, + const size_t maxRangeFilterKernelSize = 256); + + /** Azimuth common band filtering block by block @insar2007product + //TODO: since there is no windowing is applied in azimuth, no need to revert + //the windowing effects, and the antenna pattern coefficients are not stored in the + //SLC yet, will wait for the implementation and revert the attena pattern then + */ + double azimuthCommonBandFilter(std::valarray> &refSlc, + std::valarray> &secSlc, + const std::valarray &refDoppCentroids, + const std::valarray &secDoppCentroids, + std::valarray> &refAzimuthSpectrum, + std::valarray> &secAzimuthSpectrum, + isce3::signal::Filter &azimuthFilter, + size_t blockRows, + size_t ncols); + private: + void _computeDoppCentroids(const isce3::core::LUT2d & refDoppler, + const isce3::core::LUT2d & secDoppler, + isce3::io::Raster* rngOffsetRaster, + isce3::io::Raster* aziOffsetRaster, + std::valarray &refDopplerCentroids, + std::valarray &secDopplerCentroids, + std::valarray &numOfValidDopplerCentroids); + + int _computeMaxAzimuthFilterKernelSize(const std::valarray &refDopplerCentroids, + const std::valarray &secDopplerCentroids, + const std::valarray &numOfValidDopplerCentroids, + const double bandwidth, + const double prf, + isce3::signal::Filter &aziFilter); + //Doppler LUT for the refernce SLC - isce3::core::LUT1d _refDoppler; + isce3::core::LUT2d _refDoppler; //Doppler LUT for the secondary SLC - isce3::core::LUT1d _secDoppler; + isce3::core::LUT2d _secDoppler; + + // range pixel spacing in meters + double _rangePixelSpacing = 0.0; // starting range shifts between the secondary and reference RSLC in meters double _offsetStartingRangeShift = 0.0; - // range pixel spacing - double _rangePixelSpacing; + // reference starting range + double _refStartRange = 0.0; + + // reference starting azimuth time + double _refStartAzimuthTime = 0.0; + + // secondary starting range + double _secStartRange = 0.0; - // radar wavelength - double _wavelength; + // secondary starting azimuth time + double _secStartAzimuthTime = 0.0; + + // radar wavelength in meters + double _wavelength = 0.0; // number of range looks int _rangeLooks = 1; @@ -130,13 +283,58 @@ class isce3::signal::Crossmul { bool _multiLookEnabled = false; + // Flag for topo phase removal + bool _doFlatten = false; + + // Flag for common azimuth band filtering + bool _doCommonAzimuthBandFilter = false; + + // Azimuth bandwidth in Hz + double _azimuthBandwidth = 0.0; + + // Processed azimuth bandwidth after the common band filtering in Hz + double _processedAzimuthBandwidth = 0.0; + + // Sensor type + std::string _sensorType = ""; + + // Window type + std::string _windowType = "kaiser"; + + // Window parameter for constructing common band filter + double _windowParameter = 1.6; + + // Flag for common range band filtering + bool _doCommonRangeBandFilter = false; + + //pulse repetition frequency, in Hz + double _prf = 0.0; + + // range samping frequency, in Hz + double _rangeSamplingFrequency = 0.0; + + // range signal bandwidth in Hz + double _rangeBandwidth = 0.0; + + // Processed range bandwidth after the common band filtering, in Hz + double _processedRangeBandwidth = 0.0; + // number of lines per block size_t _linesPerBlock = 1024; // upsampling factor size_t _oversampleFactor = 1; + // minimum range spectrum overlap fraction between reference and secondary RSLC + double _minRangeSpectrumOverlapFraction = 0.2; + + // maximum stop band ripple amplitude, + // in dB below the pass band gain, for the FIR filter + double _ripple = 40.0; + // transition width of the common band filter(s), + // relative to the pass band width (total width of both transition bands) + double _transitionWidth = 0.15; }; // Get inline implementations for Crossmul diff --git a/cxx/isce3/signal/Crossmul.icc b/cxx/isce3/signal/Crossmul.icc index 742d63e2c..f092036ae 100644 --- a/cxx/isce3/signal/Crossmul.icc +++ b/cxx/isce3/signal/Crossmul.icc @@ -17,7 +17,7 @@ * @param[in] secSlcDoppler 2D Doppler polynomial for secondary SLC */ void isce3::signal::Crossmul:: -doppler(isce3::core::LUT1d refSlcDoppler, isce3::core::LUT1d secSlcDoppler) +doppler(isce3::core::LUT2d refSlcDoppler, isce3::core::LUT2d secSlcDoppler) { _refDoppler = refSlcDoppler; @@ -65,19 +65,18 @@ azimuthLooks(int azimuthLooks) } /** @param[in] refSpectrum the spectrum of a block of a complex data -@param[in] secSpectrum the spectrum of a block of complex data -@param[in] rangeFrequencies the frequencies in range direction -@param[in] linesPerBlock number of rows in refSpectrum and secSpectrum -@param[in] fft_size number of columns in refSpectrum and secSpectrum -@param[out] estimated frequency shift between refSpectrum and secSpectrum + @param[in] secSpectrum the spectrum of a block of complex data + @param[in] rangeFrequencies the frequencies in range direction + @param[in] linesPerBlock number of rows in refSpectrum and secSpectrum + @param[in] fft_size number of columns in refSpectrum and secSpectrum + @return estimated frequency shift between refSpectrum and secSpectrum */ -void isce3::signal::Crossmul:: -rangeFrequencyShift(std::valarray> &refSpectrum, +double isce3::signal::Crossmul:: +computeRangeFrequencyShift(std::valarray> &refSpectrum, std::valarray> &secSpectrum, std::valarray &rangeFrequencies, size_t linesPerBlock, - size_t fft_size, - double &frequencyShift) + size_t fft_size) { std::valarray refAvgSpectrum(fft_size); @@ -85,6 +84,7 @@ rangeFrequencyShift(std::valarray> &refSpectrum, std::valarray secAvgSpectrum(fft_size); secAvgSpectrum = 0.0; + #pragma omp parallel for for (size_t line = 0; line < linesPerBlock; ++line){ for (size_t col = 0; col < fft_size; ++col){ refAvgSpectrum[col] += std::abs(refSpectrum[line*fft_size + col]); @@ -98,17 +98,23 @@ rangeFrequencyShift(std::valarray> &refSpectrum, // index of the max of secAvgSpectrum size_t idx2; + // the index for the frequency shift + size_t idx; + // get the index of the max of refAvgSpectrum getPeakIndex(refAvgSpectrum, idx1); // get the index of the max of secAvgSpectrum getPeakIndex(secAvgSpectrum, idx2); - std::cout << "idx1, idx2: "<< idx1 << " , " << idx2 << std::endl; + // use the smallest idx to determine the frequency shift + // this is because their peaks can lie at two sides if + // the fftshift operation is not performed. + idx = std::min(idx1, idx2); // frequency shift is the difference of the two peaks - frequencyShift = rangeFrequencies[idx1] - rangeFrequencies[idx2]; - + // i.e., two times of the offset of the one peak + return 2.0 * rangeFrequencies[idx]; } /**@param[in] data a vector to be evaluated for the index of its maximum diff --git a/cxx/isce3/signal/Filter.cpp b/cxx/isce3/signal/Filter.cpp index c414f7c69..df8df22c1 100644 --- a/cxx/isce3/signal/Filter.cpp +++ b/cxx/isce3/signal/Filter.cpp @@ -6,6 +6,9 @@ // #include "Filter.h" +#include +#include +using namespace isce3::math::complex_operations; /** * @param[in] signal a block of data to filter @@ -24,7 +27,7 @@ initiateRangeFilter(std::valarray> &signal, size_t nrows) { _signal.forwardRangeFFT(signal, spectrum, ncols, nrows); - _signal.inverseRangeFFT(spectrum, signal, ncols, nrows); + _signal.inverseRangeFFT(spectrum, signal, ncols, nrows); } /** @@ -47,9 +50,9 @@ initiateAzimuthFilter(std::valarray> &signal, } /** - * @param[in] rangeSamplingFrequency range sampling frequency - * @param[in] subBandCenterFrequencies a vector of center frequencies for each band - * @param[in] subBandBandwidths a vector of bandwidths for each band + * @param[in] rangeSamplingFrequency range sampling frequency, in Hz + * @param[in] subBandCenterFrequencies a vector of center frequencies for each band, in Hz + * @param[in] subBandBandwidths a vector of bandwidths for each band, in Hz * @param[in] signal a block of data to filter * @param[in] spectrum a block of spectrum, which is internally used for FFT computations * @param[in] ncols number of columns of the block of data @@ -68,16 +71,15 @@ constructRangeBandpassFilter(double rangeSamplingFrequency, size_t nrows, std::string filterType) { + _signal.forwardRangeFFT(signal, spectrum, ncols, nrows); + _signal.inverseRangeFFT(spectrum, signal, ncols, nrows); + constructRangeBandpassFilter(rangeSamplingFrequency, subBandCenterFrequencies, subBandBandwidths, ncols, nrows, filterType); - - _signal.forwardRangeFFT(signal, spectrum, ncols, nrows); - _signal.inverseRangeFFT(spectrum, signal, ncols, nrows); - } template @@ -117,9 +119,22 @@ constructRangeBandpassFilter(double rangeSamplingFrequency, frequency, beta, _filter1D); - + } else if (filterType=="kaiser"){ + double beta = 1.6; + constructRangeBandpassKaiser(subBandCenterFrequencies, + subBandBandwidths, + dt, + fft_size, + frequency, + beta, + _filter1D); } else { - std::cout << filterType << " filter has not been implemented" << std::endl; + std::string err_str = filterType + " filter has not been implemented"; + pyre::journal::error_t err( + "isce.signal.Filter.constructRangeBandpassFilter"); + err << err_str << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + err_str); } //construct a block of the filter with normalization @@ -129,16 +144,70 @@ constructRangeBandpassFilter(double rangeSamplingFrequency, _filter[line*fft_size+col] = _filter1D[col] / norm; } } - } +/** + * @param[in] rangeSamplingFrequency range sampling frequency, in Hz + * @param[in] subBandCenterFrequency center frequency, in Hz + * @param[in] subBandBandwidth bandwidth, in Hz + * @param[in, out] signal a block of data to filter + * @param[in, out] spectrum a block of spectrum, which is internally used for FFT computations + * @param[in] ncols number of columns of the block of data + * @param[in] nrows number of rows of the block of data + * @param[in] filterType type of the band-pass filter + * @param[in] windowParameter filter parameter + * @param[in] maxFilterKernelSize maximum kernel size of the FIR filter in time domain + + */ +template +void +isce3::signal::Filter::constructRangeCommonBandFilter(const double rangeSamplingFrequency, + const double subBandCenterFrequency, + const double subBandBandwidth, + size_t ncols, + size_t nrows, + const std::string& filterType, + const double windowParameter, + const int maxFilterKernelSize) +{ + + int fft_size = ncols; + + _filter.resize(fft_size*nrows); + std::valarray> _filter1D(fft_size); // + _filter1D = std::complex(0.0,0.0); + + if (filterType=="kaiser"){ + constructRangeCommonBandKaiserFilter(subBandCenterFrequency, + subBandBandwidth, + rangeSamplingFrequency, + fft_size, + windowParameter, + _filter1D, + maxFilterKernelSize); + } else { + pyre::journal::error_t err( + "isce.signal.Filter.constructRangeCommonBandFilter"); + err << filterType <<" filter has not been implemented" + << pyre::journal::endl; + } + + // construct a block of the filter, and normalize the transform to recovery the + // input signal + #pragma omp parallel for + for (size_t line = 0; line < nrows; line++ ){ + for (size_t col = 0; col < fft_size; col++ ){ + _filter[line*fft_size+col] = _filter1D[col] / static_cast(fft_size); + } + } +} /** * @param[in] subBandCenterFrequencies a vector of center frequencies for each band * @param[in] subBandBandwidths a vector of bandwidths for each band - * @param[in] dt samplig rate of the signal + * @param[in] dt time interval of the signal * @param[in] fft_size length of the spectrum - * @param[out] _filter1D one dimensional boxcar bandpass filter in frequency domain + * @param[out] _filter1D one dimensional boxcar bandpass filter in frequency domain */ template void @@ -149,8 +218,8 @@ constructRangeBandpassBoxcar(std::valarray subBandCenterFrequencies, int fft_size, std::valarray>& _filter1D) { - // construct a boxcar bandpass filter in frequency domian - // which may have several bands defined by centerferquencies and + // construct a boxcar bandpass filter in frequency domain + // which may have several bands defined by centerferquencies and // subBandBandwidths for (size_t i = 0; i subBandCenterFrequencies, double fH = subBandCenterFrequencies[i] + subBandBandwidths[i]/2; //index of frequencies for fL and fH - int indL; - indexOfFrequency(dt, fft_size, fL, indL); + int indL; + indexOfFrequency(dt, fft_size, fL, indL); int indH; indexOfFrequency(dt, fft_size, fH, indH); std::cout << "fL: "<< fL << " , fH: " << fH << " indL: " << indL << " , indH: " << indH << std::endl; @@ -186,7 +255,7 @@ constructRangeBandpassBoxcar(std::valarray subBandCenterFrequencies, /** * @param[in] subBandCenterFrequencies a vector of center frequencies for each band * @param[in] subBandBandwidths a vector of bandwidths for each band - * @param[in] dt samplig rate of the signal + * @param[in] dt time interval of the signal * @param[in] frequency a vector of frequencies * @param[in] beta parameter for the raised cosine filter (0 <= beta <= 1) * @param[out] _filter1D one dimensional boxcar bandpass filter in frequency domain @@ -202,8 +271,8 @@ constructRangeBandpassCosine(std::valarray subBandCenterFrequencies, std::valarray>& _filter1D) { - const double norm = 1.0; - + const double norm = 1.0; + for (size_t i = 0; i subBandCenterFrequencies, } /** -* @param[in] refDoppler Doppler LUT1d of the reference SLC -* @param[in] secDoppler Doppler LUT1d of the secondary SLC -* @param[in] bandwidth common bandwidth in azimuth -* @param[in] prf pulse repetition frequency + * @param[in] subBandCenterFrequency the center frequency, in Hz + * @param[in] subBandBandwidth the bandwidth, in Hz + * @param[in] rangeSamplingFrequency sampling rate of the signal + * @param[in] fft_size length of the spectrum + * @param[in] beta parameter for the kaiser filter + * @param[out] filter1D one dimensional kaiser filter in frequency domain + * @param[in] maxFilterKernelSize maximum kernel size of the FIR in the time domain + */ +template +void +isce3::signal::Filter:: +constructRangeCommonBandKaiserFilter(const double subBandCenterFrequency, + const double subBandBandwidth, + const double rangeSamplingFrequency, + const int fft_size, + const double beta, + std::valarray>& filter1D, + const int maxFilterKernelSize) +{ + if (filter1D.size() <= 0) filter1D.resize(fft_size); + filter1D = std::complex(0.0, 0.0); + + isce3::signal::Signal signal; + signal.forwardRangeFFT(filter1D, filter1D, fft_size, 1); + + std::valarray> kaiser; + _design_shaped_lowpass_filter(subBandBandwidth, rangeSamplingFrequency, beta, kaiser); + + const int sizeOfKaiser = kaiser.size(); + int halfSizeOfKaiser = (sizeOfKaiser - 1)/2; + + if (maxFilterKernelSize < sizeOfKaiser) { + pyre::journal::error_t err( + "isce.signal.Filter.constructRangeCommonBandKaiserFilter"); + err << "kaiser kernel size " << sizeOfKaiser <<" is greater han maximum kernel size" + << maxFilterKernelSize + << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "kaiser kernel size is greater than maximum kernel size"); + } + + // Zero padding the filter in the middle + for (size_t ind = halfSizeOfKaiser; ind < sizeOfKaiser; ind ++) { + filter1D[ind - halfSizeOfKaiser] = kaiser[ind]; + } + for (size_t ind = 0; ind < halfSizeOfKaiser; ind ++) { + filter1D[fft_size - halfSizeOfKaiser + ind] = kaiser[ind]; + } + + // Transform to frequency domain + signal.forward(filter1D, filter1D); +} + +/** + * @param[in] subBandCenterFrequencies a vector of center frequencies for each band + * @param[in] subBandBandwidths a vector of bandwidths for each band + * @param[in] dt time interval of the signal + * @param[in] fft_size length of the spectrum + * @param[in] frequency a vector of frequencies + * @param[in] beta parameter for the kaiser filter + * @param[out] _filter1D one dimensional kaiser filter in frequency domain + */ +template +void +isce3::signal::Filter:: +constructRangeBandpassKaiser(std::valarray subBandCenterFrequencies, + std::valarray subBandBandwidths, + double dt, + int fft_size, + std::valarray& frequency, + double beta, + std::valarray>& _filter1D) +{ + // construct a kaiser bandpass filter in frequency domian + // which may have several bands defined by centerferquencies and + // subBandBandwidths + for (size_t i = 0; i=0){ + for (size_t ind = indL; ind < fft_size; ++ind){ + double fre = fL + (ind - indL) / (dt * fft_size); + double tmp = 2.0 * fre * dt; + double kaiserCoefficient = isce3::math::bessel_i0(beta * sqrt(1.0 - tmp * tmp)); + _filter1D[ind] = std::complex(kaiserCoefficient/bessel_i0_beta, 0.0); + } + for (size_t ind = 0; ind < indH; ++ind){ + double fre = ind / (dt * fft_size); + double tmp = 2.0 * fre * dt; + double kaiserCoefficient = isce3::math::bessel_i0(beta * sqrt(1.0 - tmp * tmp)); + _filter1D[ind] = std::complex(kaiserCoefficient/bessel_i0_beta, 0.0); + } + + }else{ + for (size_t ind = indL; ind < indH; ++ind){ + double fre = (ind - 0.5 * (indL + indH)) / (dt * fft_size); + double tmp = 2.0 * fre * dt; + double kaiserCoefficient = isce3::math::bessel_i0(beta * sqrt(1.0 - tmp * tmp)); + _filter1D[ind] = std::complex(kaiserCoefficient/bessel_i0_beta, 0.0); + } + } + } +} + +/** +* @param[in] refDoppler Doppler Centroids, in Hz, of the reference SLC w.r.t slant range axis +* @param[in] secDoppler Doppler Centroids, in Hz, of the secondary SLC w.r.t slant range axis +* @param[in] bandwidth input bandwidth in azimuth of the pair of SLCs, in Hz +* @param[in] prf pulse repetition frequency, in Hz +* @param[in] windowParameter window parameter of the filter +* @param[in] signal a block of data to filter +* @param[in] spectrum of the block of data +* @param[in] ncols number of columns of the block of data +* @param[in] nrows number of rows of the block of data +* @return new bandwidth +*/ +template +double +isce3::signal::Filter:: +constructAzimuthCommonBandFilter(const std::valarray & refDoppler, + const std::valarray & secDoppler, + double bandwidth, + double prf, + double windowParameter, + size_t ncols, + size_t nrows, + std::string& filterType) +{ + if (filterType=="cosine"){ + // Cosine filter is constructed in the frequency domain + return constructAzimuthCommonBandCosineFilter( + refDoppler, + secDoppler, + bandwidth, + prf, windowParameter, + ncols, nrows); + + } + if (filterType=="kaiser"){ + // Kaiser filter is constructed in the time domain + return constructAzimuthCommonBandKaiserFilter( + refDoppler, + secDoppler, + bandwidth, + prf, windowParameter, + ncols, nrows); + } + + std::string err_str = filterType + " filter has not been implemented"; + pyre::journal::error_t err( + "isce.signal.Filter.constructAzimuthCommonBandFilter"); + err << err_str << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + err_str); +} + +/** +* @param[in] refDoppler Doppler Centroids, in Hz, of the reference SLC w.r.t slant range axis +* @param[in] secDoppler Doppler Centroids, in Hz, of the secondary SLC w.r.t slant range axis +* @param[in] bandwidth input bandwidth in azimuth of the pair of SLCs, in Hz +* @param[in] prf pulse repetition frequency, in Hz * @param[in] beta parameter for raised cosine filter * @param[in] signal a block of data to filter * @param[in] spectrum of the block of data @@ -241,15 +482,13 @@ constructRangeBandpassCosine(std::valarray subBandCenterFrequencies, * @param[in] nrows number of rows of the block of data */ template -void +double isce3::signal::Filter:: -constructAzimuthCommonbandFilter(const isce3::core::LUT1d & refDoppler, - const isce3::core::LUT1d & secDoppler, +constructAzimuthCommonBandCosineFilter(const std::valarray & refDoppler, + const std::valarray & secDoppler, double bandwidth, double prf, double beta, - std::valarray> &signal, - std::valarray> &spectrum, size_t ncols, size_t nrows) { @@ -268,17 +507,39 @@ constructAzimuthCommonbandFilter(const isce3::core::LUT1d & refDoppler, // Construct vector of frequencies std::valarray frequency(fft_size); fftfreq(1.0/prf, frequency); - + + // mean doppler center frequency and + // dopper frequency shift between Reference and Secondary + double meanDopCenterFreq = 0.0; + double meanDopCenterFreqShift = 0.0; + // Loop over range bins for (int j = 0; j < ncols; ++j) { - // Compute center frequency of common band - const double fmid = 0.5 * (refDoppler.eval(j) + secDoppler.eval(j)); + double refFreq = refDoppler[j]; + double secFreq = secDoppler[j]; + + double fmid = 0.5 * (refFreq + secFreq); + double fshift = std::abs(refFreq - secFreq); + + if (fshift > bandwidth) { + pyre::journal::error_t err( + "isce.signal.Filter.constructAzimuthCommonBandCosineFilter"); + err << "Bandwith is less than frequency shift" + << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Bandwith is less than frequency shift"); + } + meanDopCenterFreqShift += fshift; + meanDopCenterFreq += fmid; // Compute filter for (size_t i = 0; i < frequency.size(); ++i) { // Get the absolute value of shifted frequency - const double freq = std::abs(frequency[i] - fmid); + double freq = std::abs(frequency[i] - fmid); + + // Wrap the frequency within the RPF + freq = freq - int(freq/prf) * prf; // Passband if (freq <= (0.5 * bandwidth - df)) { @@ -286,8 +547,8 @@ constructAzimuthCommonbandFilter(const isce3::core::LUT1d & refDoppler, // Transition region } else if (freq > (0.5 * bandwidth - df) && freq <= (0.5 * bandwidth + df)) { - _filter[i*ncols+j] = std::complex(norm * 0.5 * - (1.0 + std::cos(M_PI / (bandwidth*beta) * + _filter[i*ncols+j] = std::complex(norm * 0.5 * + (1.0 + std::cos(M_PI / (bandwidth*beta) * (freq - 0.5 * (1.0 - beta) * bandwidth))), 0.0); // Stop band @@ -305,13 +566,125 @@ constructAzimuthCommonbandFilter(const isce3::core::LUT1d & refDoppler, } } - _signal.forwardAzimuthFFT(signal, spectrum, ncols, nrows); - _signal.inverseAzimuthFFT(spectrum, signal, ncols, nrows); + meanDopCenterFreq /= ncols; + meanDopCenterFreqShift /= ncols; + std::cout << " - mean doppler center freq:" << meanDopCenterFreq << std::endl; + std::cout << " - mean doppler center freq shift:" << meanDopCenterFreqShift << std::endl; + + return (bandwidth - meanDopCenterFreqShift); +} + +/** +* @param[in] refDoppler Doppler Centroids, in Hz, of the reference SLC w.r.t slant range axis +* @param[in] secDoppler Doppler Centroids, in Hz, of the secondary SLC w.r.t slant range axis +* @param[in] bandwidth input bandwidth in azimuth of the pair of SLCs, in Hz +* @param[in] prf pulse repetition frequency, in Hz +* @param[in] beta parameter for kaiser filter +* @param[in] signal a block of data to filter +* @param[in] spectrum of the block of data +* @param[in] ncols number of columns of the block of data +* @param[in] nrows number of rows of the block of data +* @returns common azimuth bandwidth +*/ +template +double +isce3::signal::Filter:: +constructAzimuthCommonBandKaiserFilter(const std::valarray & refDoppler, + const std::valarray & secDoppler, + double bandwidth, + double prf, + double beta, + size_t ncols, + size_t nrows) +{ + pyre::journal::debug_t debug("isce.signal.Filter.constructAzimuthCommonBandKaiserFilter"); + + _filter.resize(ncols*nrows); + _filter = std::complex(0.0, 0.0); + int fft_size = nrows; + + std::valarray> filter1D(fft_size); + isce3::signal::Signal doppSignal; + doppSignal.forwardRangeFFT(filter1D, filter1D, fft_size, 1); + + // azimuth mean doppler centroid frequency, in Hz + double meanDopCenterFreq = 0.0; + + // azimuth mean doppler center frequency shift, in Hz + double meanDopCenterFreqShift = 0.0; + + // Loop over range bins + for (int j = 0; j < ncols; ++j) { + // Compute the mean doppler center frequency + double refFreq = refDoppler[j]; + double secFreq = secDoppler[j]; + + double fmid = 0.5 * (refFreq + secFreq); + double fshift = std::abs(refFreq - secFreq); + + if (bandwidth < fshift) { + pyre::journal::error_t err( + "isce.signal.Filter.constructAzimuthCommonBandKaiserFilter"); + err << "Bandwith is less than frequency shift" + << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Bandwith is less than frequency shift"); + } + + std::valarray> kaiser; + _design_shaped_lowpass_filter(bandwidth - fshift, prf, beta, kaiser); + + // Turn into the bandpass filter + _lowpass2bandpass(kaiser, fmid/prf, kaiser); + + const int sizeOfKaiser = kaiser.size(); + const int halfSizeOfKaiser = (sizeOfKaiser - 1)/2; + if (fft_size < sizeOfKaiser) { + pyre::journal::error_t err( + "isce.signal.Filter.constructAzimuthCommonBandKaiserFilter"); + err << "FFT size " << fft_size <<" is less than the filter kernel size " + << sizeOfKaiser + << pyre::journal::endl; + throw isce3::except::LengthError(ISCE_SRCINFO(), + "FFT size is less than the filter kernel size"); + } + + // Zero padding the filter in the middle + filter1D = std::complex(0.0, 0.0); + for (size_t ind = halfSizeOfKaiser; ind < sizeOfKaiser; ind ++) { + filter1D[ind - halfSizeOfKaiser] = kaiser[ind]; + } + for (size_t ind = 0; ind < halfSizeOfKaiser; ind ++) { + filter1D[fft_size - halfSizeOfKaiser + ind] = kaiser[ind]; + } + + // Transform to frequency domain + doppSignal.forward(filter1D, filter1D); + + meanDopCenterFreqShift += fshift; + meanDopCenterFreq += fmid; + + // Copy the the filter centered at fmid to a block filter + // and normalize the transform to recovery the input signal via + // dividing by fft_size + for (size_t i = 0; i < filter1D.size(); ++i) { + _filter[i*ncols+j] = filter1D[i] / static_cast(fft_size); + } + } + + meanDopCenterFreq /= ncols; + meanDopCenterFreqShift /= ncols; + + debug << " - mean doppler center freq (Hz):" << meanDopCenterFreq << pyre::journal::endl; + debug << " - mean doppler center freq shift (Hz):" << meanDopCenterFreqShift << pyre::journal::endl; + + return (bandwidth - meanDopCenterFreqShift); } + /** * @param[in] signal a block of data to filter. -* @param[in] spectrum of the block of the data +* @param[in] spectrum of the block of the data */ template void @@ -321,13 +694,13 @@ filter(std::valarray> &signal, { _signal.forward(signal, spectrum); spectrum = spectrum*_filter; - _signal.inverse(spectrum, signal); + _signal.inverse(spectrum, signal); } /** * @param[in] N length of the signal * @param[in] dt sampling interval of the signal - * @param[out] freq output vector of the frequencies + * @param[out] freq output vector of the frequencies */ void isce3::signal:: @@ -362,7 +735,7 @@ void isce3::signal::Filter:: indexOfFrequency(double dt, int N, double f, int &n) // deterrmine the index (n) of a given frequency f -// dt: sampling rate, +// dt: sampling rate, // N: length of a signal // f: frequency of interest // Assumption: for indices 0 to (N-1)/2, frequency is positive @@ -386,6 +759,157 @@ writeFilter(size_t ncols, size_t nrows) } +template +std::tuple +isce3::signal::Filter::_kaiserord(const double ripple, const double transition_width) +{ + const double A = std::abs(ripple); // in case somebody is confused as to what's meant + if (A < 8) { + // Formula for N is not valid in this range. + pyre::journal::error_t err( + "isce.signal.Filter._kaiserord"); + err << "Requested maximum ripple attentuation " << A + << " is too small for the Kaiser formula." + << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Requested maximum ripple attentuation is too small for the Kaiser formula."); + } + + // Kaiser's formula (as given in Oppenheim and Schafer) is for the filter + // order, so we have to add 1 to get the number of taps. + const double numtaps = (A - 7.95) / 2.285 / (M_PI * transition_width) + 1; + + return std::make_tuple(int(std::ceil(numtaps)), _kaiser_beta(A)); +} + +template +double +isce3::signal::Filter::_kaiser_beta(const double ripple) +{ + if (ripple > 50) { + return 0.1102 * (ripple - 8.7); + } + else if (ripple > 21){ + return 0.5842 * std::pow(ripple - 21.0, 0.4) + 0.07886 * (ripple - 21); + } + else { + return 0.0; + } +} + +template +void +isce3::signal::Filter::_kaiser_design(const double stopatt, + const double transition_width, + const bool force_odd_len, + int &n, + double &beta, + std::valarray &t) +{ + std::tie(n, beta) = _kaiserord(stopatt, transition_width); + if (force_odd_len && (n % 2 == 0)) n += 1; + + if (t.size() <= 0) t.resize(n); + for (size_t i = 0; i < t.size(); i++) t[i] = i - (n - 1) / 2.0; +} + +template +void +isce3::signal::Filter::_kaiser_irf(const std::valarray &t, + const double beta, + std::valarray> &irf) +{ + const double alpha = beta / M_PI; + const double beta0 = isce3::math::bessel_i0(beta); + + if (irf.size() <= 0) irf.resize(t.size()); + for (size_t i = 0; i < t.size(); i++) { + auto val = std::complex(t[i] * t[i] - alpha * alpha, 0.0); + auto sincVal = isce3::math::sinc(std::sqrt(val)) / beta0; + irf[i] = std::complex(sincVal.real(), 0.0); + } +} + +template +void +isce3::signal::Filter::_kaiser(const int n, + const double beta, + std::valarray> &coeffs) +{ + const double beta0 = isce3::math::bessel_i0(beta); + + if (coeffs.size() <= 0) coeffs.resize(n); + + if (n == 1) { + coeffs[0] = std::complex(1.0,0.0); + return; + } + const double alpha = (n - 1) / 2.0; + for (size_t i = 0; i < n; i++) { + const double t = (i - alpha)/alpha; + coeffs[i] = isce3::math::bessel_i0(beta * std::sqrt(1.0 - t * t)) / beta0; + } +} + +template +void +isce3::signal::Filter::_lowpass2bandpass(const std::valarray> &low_pass_filter, + const double fc, + std::valarray> &band_pass_filter) +{ + if (band_pass_filter.size() <= 0) band_pass_filter.resize(low_pass_filter.size()); + const size_t n = low_pass_filter.size(); + for (size_t i = 0; i < n; i++) { + const double t = (i - (n - 1) / 2.0); + const double phase = 2.0 * M_PI * fc * t; + const std::complex ramp_phase(std::cos(phase), std::sin(phase)); + band_pass_filter[i] = low_pass_filter[i] * ramp_phase; + } +} + +template +void +isce3::signal::Filter::_design_shaped_lowpass_filter(const double bandwidth, + const double fs, + const double window_shape, + std::valarray> &coeffs, + const double stopatt, + const double transition_width, + const bool force_odd_len) +{ + // Normalized bandwdith + const double bw = bandwidth / fs; + + // Transition width is specified in terms of output bandwidth, so scale to + // get width at sample rate of filter. + const double tw = transition_width * bw; + + if ((bw + tw / 2.0) > 1.0) { + pyre::journal::error_t err( + "isce.signal.Filter._design_shaped_lowpass_filter"); + err << "Passband + transition cannot exceed Nyquist" + << pyre::journal::endl; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Passband + transition cannot exceed Nyquist"); + } + + int n = 0; + double beta = 0.0; + std::valarray t; + std::valarray> irf; + std::valarray> kaiser; + + _kaiser_design(stopatt, tw, force_odd_len, n, beta, t); + _kaiser_irf(t * bw, window_shape, irf); + _kaiser(n, beta, kaiser); + + if (coeffs.size() <= 0) coeffs.resize(n); + + for (size_t i = 0; i < n; i++) { + coeffs[i] = irf[i] * kaiser[i] * bw; + } +} + template class isce3::signal::Filter; template class isce3::signal::Filter; diff --git a/cxx/isce3/signal/Filter.h b/cxx/isce3/signal/Filter.h index e321b48fa..80714f0ce 100644 --- a/cxx/isce3/signal/Filter.h +++ b/cxx/isce3/signal/Filter.h @@ -9,12 +9,19 @@ #include "forward.h" +#include #include +#include #include #include #include #include +#include +#include + +#include +#include #include "Signal.h" // Declaration @@ -31,7 +38,7 @@ class isce3::signal::Filter { ~Filter() {}; - /** constructs forward abd backward FFT plans for filtering a block of data in range direction. */ + /** constructs forward abd backward FFT plans for filtering a block of data in range direction. */ void initiateRangeFilter(std::valarray> &signal, std::valarray> &spectrum, size_t ncols, @@ -70,6 +77,7 @@ class isce3::signal::Filter { int fft_size, std::valarray> &_filter1D); + /** Construct a cosine range band-pass filter for multiple bands*/ void constructRangeBandpassCosine(std::valarray subBandCenterFrequencies, std::valarray subBandBandwidths, double dt, @@ -77,16 +85,67 @@ class isce3::signal::Filter { double beta, std::valarray>& _filter1D); - //T constructRangeCommonbandFilter(); + /** Construct a kaiser range band-pass filter for multiple bands*/ + void constructRangeBandpassKaiser(std::valarray subBandCenterFrequencies, + std::valarray subBandBandwidths, + double dt, + int fft_size, + std::valarray& frequency, + double beta, + std::valarray>& _filter1D); + + /** Construct the range common band filter*/ + void constructRangeCommonBandFilter(const double rangeSamplingFrequency, + const double subBandCenterFrequency, + const double subBandBandwidth, + size_t ncols, + size_t nrows, + const std::string& filterType, + const double windowParameter, + const int maxFilterKernelSize = 256); + + /** Construct a kaiser range band-pass filter for one band + * First constructs a time-domain FIR filter, then transforms to the frequency + * domain. + * Returns the frequency-domain filter coefficients. + */ + void constructRangeCommonBandKaiserFilter(const double subBandCenterFrequency, + const double subBandBandwidth, + const double rangeSamplingFrequency, + const int fft_size, + const double beta, + std::valarray>& filter1D, + const int maxFilterKernelSize = 256); /** Construct azimuth common band filter*/ - void constructAzimuthCommonbandFilter(const isce3::core::LUT1d & refDoppler, - const isce3::core::LUT1d & secDoppler, + double constructAzimuthCommonBandFilter(const std::valarray & refDoppler, + const std::valarray & secDoppler, + double bandwidth, + double prf, + double windowParameter, + size_t ncols, + size_t nrows, + std::string& filterType); + + /** Construct azimuth common band cosine filter with the doppler centroid compensation*/ + double constructAzimuthCommonBandCosineFilter(const std::valarray & refDoppler, + const std::valarray & secDoppler, + double bandwidth, + double prf, + double beta, + size_t ncols, + size_t nrows); + + /** Construct a kaiser range band-pass filter for one band + * First constructs a time-domain FIR filter, then transforms to the frequency + * domain. + * Returns the frequency-domain filter coefficients. + */ + double constructAzimuthCommonBandKaiserFilter(const std::valarray & refDoppler, + const std::valarray & secDoppler, double bandwidth, double prf, double beta, - std::valarray> &signal, - std::valarray> &spectrum, size_t ncols, size_t nrows); @@ -99,6 +158,78 @@ class isce3::signal::Filter { void writeFilter(size_t ncols, size_t nrows); + public: + /** Determine the filter window parameters for the Kaiser window method + * @param[in] ripple Upper bound for the deviation (in dB) of the magnitude of the filter's frequency response from that of the desired filter (not including frequencies in any transition intervals). + * @param[in] transition_width Width of transition region, normalized so that 1 corresponds to pi radians / sample. + * @returns the length and the beta of the Kaiser window. + */ + std::tuple _kaiserord(const double ripple, const double transition_width); + + /** Compute the Kaiser parameter `beta`, given the attenuation 'ripple` + * @param[in] ripple The desired attenuation in the stopband and maximum ripple in the passband, in dB. + * @return beta + */ + double _kaiser_beta(const double ripple); + + /** Return length, shape, and time samples for Kaiser filter design method + * @param[in] stopatt Upper bound for the deviation (in dB) of the magnitude of the filter's frequency response from that of the desired filter (not including frequencies in any transition intervals). + * @param[in] transition_width Width of transition region, normalized so that 1 corresponds to pi radians / sample. + * @param[in] force_odd_len Force to be odd length + * @param[out] n the length of the Kaiser window. + * @param[out] beta the beta parameter for the Kaiser window + * @param[out] t time samples for Kaiser filter design method + */ + void _kaiser_design(const double stopatt, + const double transition_width, + const bool force_odd_len, + int &n, + double &beta, + std::valarray &t); + + /** Impulse response (Fourier transform) of Kaiser window + * @param[in] t time samples for Kaiser filter design method + * @param[in] beta the beta parameter for the Kaiser window + * @param[out] irf time samples for Kaiser filter design method + */ + void _kaiser_irf(const std::valarray &t, + const double beta, + std::valarray> &irf); + + /** Kaiser window with length n + * @param[in] n the length of the Kaiser window. + * @param[in] beta the beta parameter for the Kaiser window + * @param[out] coeffs the Kaiser filter coefficients + */ + void _kaiser(const int n, + const double beta, + std::valarray> &coeffs); + + /** Turn a low pass filter into a band pass filter by applying a phase ramp. + * @param[in] low_pass_filter low pass filter + * @param[in] fc the center frequency, in Hz divded by sampling rate, in Hz + * @param[out] band_pass_filter band pass filter + */ + void _lowpass2bandpass(const std::valarray> &low_pass_filter, + const double fc, + std::valarray> &band_pass_filter); + + /** Design a low pass filter having a passband shaped like a window using the Kaiser method + * @param[in] bandwidth the signal bandwidth + * @param[in] fs the sampling frequency + * @param[in] beta the Kaiser window beta parameter. + * @param[in] stopatt Upper bound for the deviation (in dB) of the magnitude of the filter's frequency response from that of the desired filter (not including frequencies in any transition intervals). + * @param[in] transition_width transition width [0-1] + * @param[in] force_odd_len Force to be odd length + * @param[out] kaiser_window low bandpass kaiser window + */ + void _design_shaped_lowpass_filter(const double bandwidth, + const double fs, + const double window_shape, + std::valarray> &coeffs, + const double stopatt = 40.0, + const double transition_width = 0.15, + const bool force_odd_len = false); private: isce3::signal::Signal _signal; std::valarray> _filter; diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 3e253c187..b5ddd8f22 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -37,6 +37,14 @@ @article{envisat2007product url={https://earth.esa.int/pub/ESA_DOC/ENVISAT/ASAR/asar.ProductHandbook.2_2.pdf}, } +@article{insar2007product, + title={InSAR processing:a practical approach}, + author={European Space Agency}, + journal={European Space}, + year={2007}, + url={https://www.esa.int/esapub/tm/tm19/TM-19_ptB.pdf}, +} + @article{small2008guide, title={Guide to ASAR geocoding}, author={Small, David and Schubert, Adrian}, @@ -155,8 +163,8 @@ @article{breit2010 @ARTICLE{small2011, author={D. {Small}}, - journal={IEEE Transactions on Geoscience and Remote Sensing}, - title={Flattening Gamma: Radiometric Terrain Correction for SAR Imagery}, + journal={IEEE Transactions on Geoscience and Remote Sensing}, + title={Flattening Gamma: Radiometric Terrain Correction for SAR Imagery}, year={2011}, volume={49}, number={8}, diff --git a/python/extensions/pybind_isce3/cuda/signal/Crossmul.cpp b/python/extensions/pybind_isce3/cuda/signal/Crossmul.cpp index 12c5f81ff..7f58b9845 100644 --- a/python/extensions/pybind_isce3/cuda/signal/Crossmul.cpp +++ b/python/extensions/pybind_isce3/cuda/signal/Crossmul.cpp @@ -30,7 +30,8 @@ void addbinding(py::class_ & pyCrossmul) py::arg("sec_slc"), py::arg("interferogram"), py::arg("coherence"), - py::arg("range_offset") = nullptr, R"( + py::arg("range_offset") = nullptr, + py::arg("azimuth_offset") = nullptr, R"( Crossmultiply reference and secondary SLCs to generate interferogram and coherence products. Parameters @@ -44,17 +45,18 @@ void addbinding(py::class_ & pyCrossmul) coherence: Raster Output coherence raster interferogram: Raster - Optional range offset raster usef for flattening + Optional range offset raster usef for flattening and common band filtering + Optional azimuth offset raster usef for azimuth common band filtering )") .def("set_dopplers", &gpuCrossmul::doppler, py::arg("ref_doppler"), py::arg("sec_doppler")) .def_property("ref_doppler", py::overload_cast<>(&gpuCrossmul::refDoppler, py::const_), - py::overload_cast>(&gpuCrossmul::refDoppler)) + py::overload_cast>(&gpuCrossmul::refDoppler)) .def_property("sec_doppler", py::overload_cast<>(&gpuCrossmul::secDoppler, py::const_), - py::overload_cast>(&gpuCrossmul::secDoppler)) + py::overload_cast>(&gpuCrossmul::secDoppler)) .def_property("ref_sec_offset_starting_range_shift", py::overload_cast<>(&gpuCrossmul::startingRangeShift, py::const_), py::overload_cast(&gpuCrossmul::startingRangeShift)) diff --git a/python/extensions/pybind_isce3/signal/Crossmul.cpp b/python/extensions/pybind_isce3/signal/Crossmul.cpp index 613d33f47..7ba31cfa6 100644 --- a/python/extensions/pybind_isce3/signal/Crossmul.cpp +++ b/python/extensions/pybind_isce3/signal/Crossmul.cpp @@ -30,7 +30,8 @@ void addbinding(py::class_ & pyCrossmul) py::arg("sec_slc"), py::arg("interferogram"), py::arg("coherence"), - py::arg("range_offset") = nullptr, R"( + py::arg("range_offset") = nullptr, + py::arg("azimuth_offset") = nullptr, R"( Crossmultiply reference and secondary SLCs to generate interferogram and coherence products. Parameters @@ -43,24 +44,69 @@ void addbinding(py::class_ & pyCrossmul) Output interferogram raster coherence: Raster Output coherence raster - interferogram: Raster - Optional range offset raster usef for flattening + range_offset: Raster + Optional range offset raster usef for flattening and common band filter + azimuth_offset: Raster + Optional azimuth offset raster usef for azimuth common band filter + range_bandwidth: float + range bandwidth of the reference SLC + azimuth_bandwidth: float + azimuth bandwidth of the reference SLC )") .def("set_dopplers", &Crossmul::doppler, py::arg("ref_doppler"), py::arg("sec_doppler")) .def_property("ref_doppler", py::overload_cast<>(&Crossmul::refDoppler, py::const_), - py::overload_cast>(&Crossmul::refDoppler)) + py::overload_cast>(&Crossmul::refDoppler)) .def_property("sec_doppler", py::overload_cast<>(&Crossmul::secDoppler, py::const_), - py::overload_cast>(&Crossmul::secDoppler)) + py::overload_cast>(&Crossmul::secDoppler)) .def_property("ref_sec_offset_starting_range_shift", py::overload_cast<>(&Crossmul::startingRangeShift, py::const_), py::overload_cast(&Crossmul::startingRangeShift)) .def_property("range_pixel_spacing", py::overload_cast<>(&Crossmul::rangePixelSpacing, py::const_), py::overload_cast(&Crossmul::rangePixelSpacing)) + .def_property("window_parameter", + py::overload_cast<>(&Crossmul::windowParameter, py::const_), + py::overload_cast(&Crossmul::windowParameter)) + .def_property("sensor_type", + py::overload_cast<>(&Crossmul::sensorType, py::const_), + py::overload_cast(&Crossmul::sensorType)) + .def_property("window_type", + py::overload_cast<>(&Crossmul::windowType, py::const_), + py::overload_cast(&Crossmul::windowType)) + .def_property("do_common_range_band_filter", + py::overload_cast<>(&Crossmul::doCommonRangeBandFilter, py::const_), + py::overload_cast(&Crossmul::doCommonRangeBandFilter)) + .def_property("do_common_azimuth_band_filter", + py::overload_cast<>(&Crossmul::doCommonAzimuthBandFilter, py::const_), + py::overload_cast(&Crossmul::doCommonAzimuthBandFilter)) + .def_property("do_flatten", + py::overload_cast<>(&Crossmul::doFlatten, py::const_), + py::overload_cast(&Crossmul::doFlatten)) + .def_property("range_bandwidth", + py::overload_cast<>(&Crossmul::rangeBandwidth, py::const_), + py::overload_cast(&Crossmul::rangeBandwidth)) + .def_property("ref_start_range", + py::overload_cast<>(&Crossmul::refStartRange, py::const_), + py::overload_cast(&Crossmul::refStartRange)) + .def_property("ref_start_azimuth_time", + py::overload_cast<>(&Crossmul::refStartAzimuthTime, py::const_), + py::overload_cast(&Crossmul::refStartAzimuthTime)) + .def_property("sec_start_range", + py::overload_cast<>(&Crossmul::secStartRange, py::const_), + py::overload_cast(&Crossmul::secStartRange)) + .def_property("sec_start_azimuth_time", + py::overload_cast<>(&Crossmul::secStartAzimuthTime, py::const_), + py::overload_cast(&Crossmul::secStartAzimuthTime)) + .def_property("azimuth_bandwidth", + py::overload_cast<>(&Crossmul::azimuthBandwidth, py::const_), + py::overload_cast(&Crossmul::azimuthBandwidth)) + .def_property("prf", + py::overload_cast<>(&Crossmul::prf, py::const_), + py::overload_cast(&Crossmul::prf)) .def_property("wavelength", py::overload_cast<>(&Crossmul::wavelength, py::const_), py::overload_cast(&Crossmul::wavelength)) @@ -77,5 +123,7 @@ void addbinding(py::class_ & pyCrossmul) py::overload_cast<>(&Crossmul::linesPerBlock, py::const_), py::overload_cast(&Crossmul::linesPerBlock)) .def_property_readonly("multilook_enabled", &Crossmul::multiLookEnabled) + .def_property_readonly("processed_range_bandwidth", &Crossmul::processedRangeBandwidth) + .def_property_readonly("processed_azimuth_bandwidth", &Crossmul::processedAzimuthBandwidth) ; } diff --git a/python/packages/nisar/workflows/crossmul.py b/python/packages/nisar/workflows/crossmul.py index 4b79fb65a..73768dd5a 100644 --- a/python/packages/nisar/workflows/crossmul.py +++ b/python/packages/nisar/workflows/crossmul.py @@ -14,6 +14,7 @@ from isce3.io import HDF5OptimizedReader from nisar.products.readers import SLC + from nisar.workflows import prepare_insar_hdf5 from nisar.workflows.compute_stats import compute_stats_real_data from nisar.workflows.crossmul_runconfig import CrossmulRunConfig @@ -34,6 +35,9 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', crossmul_params = cfg['processing']['crossmul'] scratch_path = pathlib.Path(cfg['product_path_group']['scratch_path']) flatten = crossmul_params['flatten'] + do_common_range_band_filter = crossmul_params['common_band_range_filter'] + do_common_azimuth_band_filter = crossmul_params['common_band_azimuth_filter'] + lines_per_block = crossmul_params['lines_per_block'] if rg_looks == None: @@ -41,7 +45,7 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', if az_looks == None: az_looks = crossmul_params['azimuth_looks'] - if flatten: + if flatten or do_common_range_band_filter or do_common_azimuth_band_filter: flatten_path = crossmul_params['flatten_path'] if output_hdf5 is None: @@ -63,8 +67,22 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', device = isce3.cuda.core.Device(cfg['worker']['gpu_id']) isce3.cuda.core.set_device(device) crossmul = isce3.cuda.signal.Crossmul() + + if do_common_range_band_filter or do_common_azimuth_band_filter: + raise NotImplementedError("Common band filters have not been implemented for GPU") else: crossmul = isce3.signal.Crossmul() + # do common range band filter + crossmul.do_common_range_band_filter = \ + do_common_range_band_filter + # do common azimuth band filter + crossmul.do_common_azimuth_band_filter = \ + do_common_azimuth_band_filter + # do the flatten + crossmul.do_flatten = flatten + # sensor type + crossmul.sensor_type = \ + ref_slc.identification.missionId crossmul.range_looks = rg_looks crossmul.az_looks = az_looks @@ -87,10 +105,9 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', crossmul_dir = scratch_path / f'crossmul/freq{freq}' crossmul_dir.mkdir(parents=True, exist_ok=True) # get 2d doppler, discard azimuth dependency, and set crossmul dopplers - ref_dopp = isce3.core.avg_lut2d_to_lut1d( - ref_slc.getDopplerCentroid(frequency=freq)) - sec_dopp = isce3.core.avg_lut2d_to_lut1d( - sec_slc.getDopplerCentroid(frequency=freq)) + ref_dopp = ref_slc.getDopplerCentroid(frequency=freq) + sec_dopp = sec_slc.getDopplerCentroid(frequency=freq) + crossmul.set_dopplers(ref_dopp, sec_dopp) freq_group_path = f'{RIFGGroupsPaths().SwathsPath}/frequency{freq}' @@ -100,12 +117,38 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', crossmul.range_pixel_spacing = ref_radar_grid.range_pixel_spacing crossmul.wavelength = ref_radar_grid.wavelength + # CPU version + # TODO: Add parameters to GPU implementation + if not use_gpu: + + sec_rdr_grid = sec_slc.getRadarGrid(freq) + + # range bandwidth + crossmul.range_bandwidth = \ + ref_slc.getSwathMetadata(freq).processed_range_bandwidth + # azimuth band width and PRF + crossmul.azimuth_bandwidth = \ + ref_slc.getSwathMetadata(freq).processed_azimuth_bandwidth + crossmul.prf = ref_radar_grid.prf + + # start range and azimuth time for reference and secondary images + crossmul.ref_start_range = ref_radar_grid.starting_range + crossmul.sec_start_range = sec_rdr_grid.starting_range + crossmul.ref_start_azimuth_time = ref_radar_grid.sensing_start + crossmul.sec_start_azimuth_time = sec_rdr_grid.sensing_start + # enable/disable flatten accordingly - if flatten: + if flatten or do_common_range_band_filter or do_common_azimuth_band_filter: # set frequency dependent range offset raster - flatten_raster = isce3.io.Raster( + range_offsets_raster = isce3.io.Raster( f'{flatten_path}/geo2rdr/freq{freq}/range.off') + if do_common_azimuth_band_filter: + azimuth_offsets_raster = isce3.io.Raster( + f'{flatten_path}/geo2rdr/freq{freq}/azimuth.off') + else: + azimuth_offsets_raster = None + # Calculate the starting range shift between reference and secondary in meters sec_radar_grid = sec_slc.getRadarGrid(freq) rng_shift = (sec_radar_grid.starting_range - @@ -114,7 +157,8 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', crossmul.ref_sec_offset_starting_range_shift\ = rng_shift else: - flatten_raster = None + range_offsets_raster = None + azimuth_offsets_raster = None for pol in pol_list: output_dir = crossmul_dir / f'{pol}' @@ -166,7 +210,25 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse', # Compute multilooked interferogram and coherence raster crossmul.crossmul(ref_slc_raster, sec_slc_raster, ifg_raster, - coh_raster, flatten_raster) + coh_raster, + range_offsets_raster, + azimuth_offsets_raster) + + # populate the new bandwidth along azimuth and range after the common band filter + # if there is no common band filter applied, the bandwith will remain the same with + # the orignal SLC bandwidth. + # NOTE: Those bandwidths have already been in the dataset. + # TODO: GPU + if not use_gpu: + processing_info_path = \ + RIFGGroupsPaths().ProcessingInformationPath + ifgram_processing_parameter = \ + f'{processing_info_path}/parameters/interferogram/frequency{freq}' + # Update the bandwidth + dst_h5[f'{ifgram_processing_parameter}/azimuthBandwidth'][...] = \ + crossmul.processed_azimuth_bandwidth + dst_h5[f'{ifgram_processing_parameter}/rangeBandwidth'][...] = \ + crossmul.processed_range_bandwidth del ifg_raster diff --git a/python/packages/nisar/workflows/crossmul_runconfig.py b/python/packages/nisar/workflows/crossmul_runconfig.py index a3d5fa060..53e571c13 100644 --- a/python/packages/nisar/workflows/crossmul_runconfig.py +++ b/python/packages/nisar/workflows/crossmul_runconfig.py @@ -62,7 +62,11 @@ def yaml_check(self, resample_type): # flatten defaults to bool True # flatten_path defaults to scratch_path flatten = crossmul_cfg['flatten'] - if flatten: + do_common_band_range_filter = crossmul_cfg['common_band_range_filter'] + do_common_band_azimuth_filter = crossmul_cfg['common_band_azimuth_filter'] + do_common_band_filter = \ + do_common_band_range_filter or do_common_band_azimuth_filter + if flatten or do_common_band_filter: if 'flatten_path' in crossmul_cfg: flatten_path = crossmul_cfg['flatten_path'] else: diff --git a/python/packages/nisar/workflows/troposphere.py b/python/packages/nisar/workflows/troposphere.py index b69cb535f..fa93fd1b9 100644 --- a/python/packages/nisar/workflows/troposphere.py +++ b/python/packages/nisar/workflows/troposphere.py @@ -162,7 +162,6 @@ def _convert_HRES_to_raider_NetCDF(weather_model_file, weather_model_output_dir): ''' Internal convenience function to convert the ECMWF NetCDF to RAiDER NetCDF - Parameters ---------- weather_model_file: str diff --git a/tests/cxx/isce3/cuda/signal/gpuCrossMul.cpp b/tests/cxx/isce3/cuda/signal/gpuCrossMul.cpp index 2e5d85530..be9490326 100644 --- a/tests/cxx/isce3/cuda/signal/gpuCrossMul.cpp +++ b/tests/cxx/isce3/cuda/signal/gpuCrossMul.cpp @@ -48,18 +48,18 @@ TEST(gpuCrossmul, Crossmul) // Create a product isce3::product::RadarGridProduct product(file); - // get the Doppler polynomial for reference SLC - const isce3::core::LUT1d dop1 = - avgLUT2dToLUT1d(product.metadata().procInfo().dopplerCentroid('A')); + // get the Doppler for refernce SLC + const isce3::core::LUT2d dop1 = + product.metadata().procInfo().dopplerCentroid('A'); // Since this test careates an interferogram between the reference SLC and itself, // the second Doppler is the same as the first - isce3::core::LUT1d dop2 = dop1; + isce3::core::LUT2d dop2 = dop1; //instantiate the Crossmul class isce3::cuda::signal::gpuCrossmul crsmul; - // set Doppler polynomials for reference and secondary SLCs + // set Doppler for refernce and secondary SLCs crsmul.doppler(dop1, dop2); // set number of interferogram looks in range @@ -124,18 +124,18 @@ TEST(gpuCrossmul, MultilookCrossmul) // Create a product isce3::product::RadarGridProduct product(file); - // get the Doppler polynomial for reference SLC - isce3::core::LUT1d dop1 = - avgLUT2dToLUT1d(product.metadata().procInfo().dopplerCentroid('A')); + // get the Doppler polynomial for refernce SLC + isce3::core::LUT2d dop1 = + product.metadata().procInfo().dopplerCentroid('A'); // Since this test careates an interferogram between the reference SLC and itself, // the second Doppler is the same as the first - isce3::core::LUT1d dop2 = dop1; + isce3::core::LUT2d dop2 = dop1; //instantiate the Crossmul class isce3::cuda::signal::gpuCrossmul crsmul; - // set Doppler polynomials for reference and secondary SLCs + // set Doppler for refernce and secondary SLCs crsmul.doppler(dop1, dop2); // set number of interferogram looks in range diff --git a/tests/cxx/isce3/signal/crossmul.cpp b/tests/cxx/isce3/signal/crossmul.cpp index 664c2b1ee..5ee1dbb99 100644 --- a/tests/cxx/isce3/signal/crossmul.cpp +++ b/tests/cxx/isce3/signal/crossmul.cpp @@ -26,6 +26,10 @@ TEST(Crossmul, RunCrossmul) //a raster object for the reference SLC isce3::io::Raster referenceSlc(TESTDATA_DIR "warped_envisat.slc.vrt"); + //a raster object for the range and azimuth offsets + isce3::io::Raster aziOffsets(TESTDATA_DIR "envisat_offsets.tif"); + isce3::io::Raster rngOffsets(TESTDATA_DIR "envisat_offsets.tif"); + // get the length and width of the SLC int width = referenceSlc.width(); int length = referenceSlc.length(); @@ -47,18 +51,19 @@ TEST(Crossmul, RunCrossmul) // Create a product isce3::product::RadarGridProduct product(file); - // get the Doppler polynomial for reference SLC - isce3::core::LUT1d dop1 = - avgLUT2dToLUT1d(product.metadata().procInfo().dopplerCentroid('A')); + // get the Doppler for refernce SLC + const char freq = 'A'; + isce3::core::LUT2d dop1 = product.metadata().procInfo().dopplerCentroid(freq); + auto swath = product.swath(freq); // Since this test careates an interferogram between the reference SLC and itself, // the second Doppler is the same as the first - isce3::core::LUT1d dop2 = dop1; + isce3::core::LUT2d dop2 = dop1; //instantiate the Crossmul class isce3::signal::Crossmul crsmul; - // set Doppler polynomials for reference and secondary SLCs + // set Doppler for refernce and secondary SLCs crsmul.doppler(dop1, dop2); // set number of interferogram looks in range @@ -67,8 +72,34 @@ TEST(Crossmul, RunCrossmul) // set number of interferogram looks in azimuth crsmul.azimuthLooks(1); + // set the product information for the common band filtering and flattenning + const double wavelength = swath.processedWavelength(); + const double azimuthBandwidth = swath.processedAzimuthBandwidth(); + const double rangeBandwidth = swath.processedRangeBandwidth(); + const double prf = swath.nominalAcquisitionPRF(); + const double rngPixelSampling = swath.rangePixelSpacing(); + const double startAzimuthTime = swath.zeroDopplerTime()[0]; + const double startSlantRange = swath.slantRange()[0]; + + crsmul.wavelength(wavelength); + crsmul.prf(prf); + crsmul.rangePixelSpacing(rngPixelSampling); + crsmul.azimuthBandwidth(azimuthBandwidth); + crsmul.rangeBandwidth(rangeBandwidth); + crsmul.refStartRange(startSlantRange); + crsmul.secStartRange(startSlantRange); + crsmul.refStartAzimuthTime(startAzimuthTime); + crsmul.secStartAzimuthTime(startAzimuthTime); + + // Enable the flattening and common band filtering + crsmul.doFlatten(true); + crsmul.doCommonRangeBandFilter(true); + crsmul.doCommonAzimuthBandFilter(true); + // running crossmul - crsmul.crossmul(referenceSlc, referenceSlc, interferogram, coherence); + crsmul.crossmul(referenceSlc, referenceSlc, + interferogram, coherence, + &rngOffsets, &aziOffsets); // an array for the computed interferogram std::valarray> data(width*length); @@ -97,6 +128,10 @@ TEST(Crossmul, RunCrossmulMLook) //a raster object for the reference SLC isce3::io::Raster referenceSlc(TESTDATA_DIR "warped_envisat.slc.vrt"); + //a raster object for the range and azimuth offsets + isce3::io::Raster aziOffsets(TESTDATA_DIR "envisat_offsets.tif"); + isce3::io::Raster rngOffsets(TESTDATA_DIR "envisat_offsets.tif"); + // define looks const int rngLooks = 3; const int azLooks = 13; @@ -122,13 +157,14 @@ TEST(Crossmul, RunCrossmulMLook) // Create a product isce3::product::RadarGridProduct product(file); - // get the Doppler polynomial for reference SLC - isce3::core::LUT1d dop1 = - avgLUT2dToLUT1d(product.metadata().procInfo().dopplerCentroid('A')); + // get the Doppler polynomial for refernce SLC + const char freq = 'A'; + isce3::core::LUT2d dop1 = product.metadata().procInfo().dopplerCentroid(freq); + auto swath = product.swath(freq); // Since this test careates an interferogram between the reference SLC and itself, // the second Doppler is the same as the first - isce3::core::LUT1d dop2 = dop1; + isce3::core::LUT2d dop2 = dop1; //instantiate the Crossmul class isce3::signal::Crossmul crsmul; @@ -142,8 +178,34 @@ TEST(Crossmul, RunCrossmulMLook) // set number of interferogram looks in azimuth crsmul.azimuthLooks(azLooks); + // set the product information for the common band filtering and flattenning + const double wavelength = swath.processedWavelength(); + const double azimuthBandwidth = swath.processedAzimuthBandwidth(); + const double rangeBandwidth = swath.processedRangeBandwidth(); + const double prf = swath.nominalAcquisitionPRF(); + const double rngPixelSampling = swath.rangePixelSpacing(); + const double startAzimuthTime = swath.zeroDopplerTime()[0]; + const double startSlantRange = swath.slantRange()[0]; + + crsmul.wavelength(wavelength); + crsmul.prf(prf); + crsmul.rangePixelSpacing(rngPixelSampling); + crsmul.azimuthBandwidth(azimuthBandwidth); + crsmul.rangeBandwidth(rangeBandwidth); + crsmul.refStartRange(startSlantRange); + crsmul.secStartRange(startSlantRange); + crsmul.refStartAzimuthTime(startAzimuthTime); + crsmul.secStartAzimuthTime(startAzimuthTime); + + // Enable the flattening and common band filtering + crsmul.doFlatten(true); + crsmul.doCommonRangeBandFilter(true); + crsmul.doCommonAzimuthBandFilter(true); + // running crossmul - crsmul.crossmul(referenceSlc, referenceSlc, interferogram, coherence); + crsmul.crossmul(referenceSlc, referenceSlc, + interferogram, coherence, + &rngOffsets, &aziOffsets); // an array for the computed interferogram std::valarray> data(width*length); diff --git a/tests/cxx/isce3/signal/filter.cpp b/tests/cxx/isce3/signal/filter.cpp index 55ccdf070..312344127 100644 --- a/tests/cxx/isce3/signal/filter.cpp +++ b/tests/cxx/isce3/signal/filter.cpp @@ -14,7 +14,7 @@ #include #include -TEST(Filter, constructAzimuthCommonbandFilter) +TEST(Filter, constructAzimuthCommonBandFilter) { //This test constructs a common azimuth band filter. @@ -35,9 +35,9 @@ TEST(Filter, constructAzimuthCommonbandFilter) isce3::product::RadarGridProduct product(file); const isce3::product::Swath & swath = product.swath('A'); - // Get the Doppler polynomial and use it for both reference and secondary SLCs - auto dop1 = isce3::core::avgLUT2dToLUT1d(product.metadata().procInfo().dopplerCentroid('A')); - auto dop2 = dop1; + // Get the Doppler polynomial and use it for both refernce and secondary SLCs + std::valarray dop1(nfft); + std::valarray dop2(nfft); // get pulase repetition frequency (prf) double prf = swath.nominalAcquisitionPRF(); @@ -50,15 +50,13 @@ TEST(Filter, constructAzimuthCommonbandFilter) double commonAzimuthBandwidth = 1000.0; isce3::signal::Filter filter; - filter.constructAzimuthCommonbandFilter(dop1, - dop2, - commonAzimuthBandwidth, - prf, - beta, - refSlc, refSpectrum, - ncols, blockRows); + filter.constructAzimuthCommonBandCosineFilter(dop1, + dop2, + commonAzimuthBandwidth, + prf, + beta, + ncols, blockRows); filter.writeFilter(ncols, blockRows); - } TEST(Filter, constructBoxcarRangeBandpassFilter) @@ -80,10 +78,10 @@ TEST(Filter, constructBoxcarRangeBandpassFilter) // get the range bandwidth double BW = swath.processedRangeBandwidth(); - + //The bands are specified by two vectors: // 1) a vector of center frequencies for each sub-band - std::valarray subBandCenterFrequencies{-3.0e6, 0.0, 3e6}; + std::valarray subBandCenterFrequencies{-3.0e6, 0.0, 3e6}; // 2) a vector of bandwidth of each sub-band std::valarray subBandBandwidths{2.0e6, 2.0e6, 2.0e6}; @@ -101,7 +99,7 @@ TEST(Filter, constructBoxcarRangeBandpassFilter) blockRows, filterType); - //filter.writeFilter(ncols, blockRows); + filter.writeFilter(ncols, blockRows); // change the filter type to cosine filterType = "cosine"; @@ -114,8 +112,31 @@ TEST(Filter, constructBoxcarRangeBandpassFilter) blockRows, filterType); - //filter.writeFilter(ncols, blockRows); - + filter.writeFilter(ncols, blockRows); +} + +TEST(Filter, constructRangeCommonBandKaiserFilter) +{ + //This test constructs a kaiser common band range band-pass FIR filter. + int fft_size = 256; + std::valarray> kaiser; + + double subBandCenterFrequency = 0.0; + double subBandBandwidth = 2.0e6; + // Assume range sampling frequency equals 1.2 times bandwidth for this test + double rangeSamplingFrequency = 1.2 * subBandBandwidth; + + isce3::signal::Filter filter; + filter.constructRangeCommonBandKaiserFilter(subBandCenterFrequency, + subBandBandwidth, + rangeSamplingFrequency, + fft_size, + 1.6, + kaiser); + + ASSERT_LT(std::abs(std::abs(kaiser[0]) - 0.9997906684875488), 1.0e-6); + ASSERT_LT(std::abs(std::abs(kaiser[127]) - 0.003600120544433594), 1.0e-6); + ASSERT_LT(std::abs(std::abs(kaiser[255]) - 0.99970543384552), 1.0e-6); } int main(int argc, char * argv[]) { diff --git a/tests/data/README.md b/tests/data/README.md index f101f795f..5a21c41d4 100644 --- a/tests/data/README.md +++ b/tests/data/README.md @@ -27,6 +27,17 @@ $ md5sum NISAR_ANC_L_PR_FRP_20210114T023357_20200107T200000_20200109T040000.xml To reduce the size of the test data, this file was trimmed to the first ten records to produce `attitude.xml`. +## Crossmul + +- **envisat.h5** + + Simulated NISAR SLC using the envisat data with the image dimenson of 500 by 500. + +- **envisat_offsets.tif** + + The range and azimuth offests of the envisat.h5, and all values are 0. + + ## Orbit Orbit sample data in NISAR format (see JPL D-102253) provided by Paul Ries diff --git a/tests/data/envisat_offsets.tif b/tests/data/envisat_offsets.tif new file mode 100644 index 000000000..bf86bbee7 Binary files /dev/null and b/tests/data/envisat_offsets.tif differ diff --git a/tests/python/extensions/pybind/cuda/signal/crossmul.py b/tests/python/extensions/pybind/cuda/signal/crossmul.py index aa83e172d..dbc702405 100644 --- a/tests/python/extensions/pybind/cuda/signal/crossmul.py +++ b/tests/python/extensions/pybind/cuda/signal/crossmul.py @@ -21,7 +21,7 @@ def common_crossmul_obj(): ''' # make SLC object and extract parameters slc_obj = SLC(hdf5file=os.path.join(iscetest.data, 'envisat.h5')) - dopp = isce3.core.avg_lut2d_to_lut1d(slc_obj.getDopplerCentroid()) + dopp = slc_obj.getDopplerCentroid() prf = slc_obj.getRadarGrid().prf crossmul = isce3.cuda.signal.Crossmul() diff --git a/tests/python/extensions/pybind/signal/crossmul.py b/tests/python/extensions/pybind/signal/crossmul.py index d5b5b6920..18ec878c9 100644 --- a/tests/python/extensions/pybind/signal/crossmul.py +++ b/tests/python/extensions/pybind/signal/crossmul.py @@ -21,7 +21,7 @@ def common_crossmul_obj(): ''' # make SLC object and extract parameters slc_obj = SLC(hdf5file=os.path.join(iscetest.data, 'envisat.h5')) - dopp = isce3.core.avg_lut2d_to_lut1d(slc_obj.getDopplerCentroid()) + dopp = slc_obj.getDopplerCentroid() prf = slc_obj.getRadarGrid().prf crossmul = isce3.signal.Crossmul()