Skip to content

Commit 15b3aec

Browse files
hfattahiHeresh Fattahi
authored and
GitHub Enterprise
committed
GSLC interpolation consistent with resample (#805)
* GSLC interpolation consistent with resample * remove sin length from input arguments * clang-fromat * half pixel shift to get to the center of each geo pixel * fix documentation Co-authored-by: Heresh Fattahi <[email protected]>
1 parent 6d4d9c9 commit 15b3aec

File tree

3 files changed

+90
-76
lines changed

3 files changed

+90
-76
lines changed

‎cxx/isce3/geocode/geocodeSlc.cpp‎

Lines changed: 8 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -78,10 +78,6 @@ void isce3::geocode::geocodeSlc(
7878
std::valarray<double> radarX(blockSize);
7979
std::valarray<double> radarY(blockSize);
8080

81-
// container for the sum of the carrier phase (Doppler) to be added back
82-
// and the geometrical phase to be removed for flattening the SLC phase.
83-
std::valarray<std::complex<double>> geometricalPhase(blockSize);
84-
8581
int localAzimuthFirstLine = radarGrid.length() - 1;
8682
int localAzimuthLastLine = 0;
8783
int localRangeFirstPixel = radarGrid.width() - 1;
@@ -102,11 +98,14 @@ void isce3::geocode::geocodeSlc(
10298
// Global line index
10399
const size_t line = lineStart + blockLine;
104100

105-
// y coordinate in the out put grid
106-
double y = geoGrid.startY() + geoGrid.spacingY() * line;
101+
// y coordinate in the out put grid
102+
// Assuming geoGrid.startY() and geoGrid.startX() represent the top-left
103+
// corner of the first pixel, then 0.5 pixel shift is needed to get
104+
// to the center of each pixel
105+
double y = geoGrid.startY() + geoGrid.spacingY() * (line + 0.5);
107106

108107
// x in the output geocoded Grid
109-
double x = geoGrid.startX() + geoGrid.spacingX() * pixel;
108+
double x = geoGrid.startX() + geoGrid.spacingX() * (pixel + 0.5);
110109

111110
// compute the azimuth time and slant range for the
112111
// x,y coordinates in the output grid
@@ -158,20 +157,6 @@ void isce3::geocode::geocodeSlc(
158157
radarX[blockLine * geoGrid.width() + pixel] = rdrX;
159158
radarY[blockLine * geoGrid.width() + pixel] = rdrY;
160159

161-
// doppler to be added back after interpolation
162-
double phase =
163-
nativeDoppler.eval(aztime, srange) * 2 * M_PI * aztime;
164-
//
165-
166-
if (flatten) {
167-
phase += (4.0 * (M_PI / radarGrid.wavelength())) * srange;
168-
}
169-
170-
const std::complex<double> cpxPhase(std::cos(phase),
171-
std::sin(phase));
172-
173-
geometricalPhase[blockLine * geoGrid.width() + pixel] = cpxPhase;
174-
175160
} // end loops over lines and pixel of output grid
176161

177162
// Get min and max swath extents from among all threads
@@ -208,28 +193,10 @@ void isce3::geocode::geocodeSlc(
208193
azimuthFirstLine, rdrBlockWidth,
209194
rdrBlockLength, band + 1);
210195

211-
// baseband the SLC in the radar grid
212-
const double blockStartingRange =
213-
radarGrid.startingRange() +
214-
rangeFirstPixel * radarGrid.rangePixelSpacing();
215-
const double blockSensingStart = radarGrid.sensingStart() +
216-
azimuthFirstLine / radarGrid.prf();
217-
218-
isce3::geocode::baseband(rdrDataBlock, blockStartingRange,
219-
blockSensingStart,
220-
radarGrid.rangePixelSpacing(),
221-
radarGrid.prf(), nativeDoppler);
222-
223196
// interpolate the data in radar grid to the geocoded grid.
224-
// Also the geometrical phase, which is the phase of the carrier
225-
// to be added back and the geometrical phase to be removed is
226-
// applied.
227-
std::cout << "resample " << std::endl;
228197
isce3::geocode::interpolate(rdrDataBlock, geoDataBlock, radarX,
229-
radarY, geometricalPhase, rdrBlockWidth,
230-
rdrBlockLength, azimuthFirstLine,
231-
rangeFirstPixel, interp.get());
232-
198+
radarY, azimuthFirstLine, rangeFirstPixel, interp.get(),
199+
radarGrid, nativeDoppler, flatten);
233200
// set output
234201
std::cout << "set output " << std::endl;
235202
outputRaster.setBlock(geoDataBlock.data(), 0, lineStart,

‎cxx/isce3/geocode/interpolate.cpp‎

Lines changed: 67 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -4,20 +4,21 @@ void isce3::geocode::interpolate(
44
const isce3::core::Matrix<std::complex<float>>& rdrDataBlock,
55
isce3::core::Matrix<std::complex<float>>& geoDataBlock,
66
const std::valarray<double>& radarX,
7-
const std::valarray<double>& radarY,
8-
const std::valarray<std::complex<double>>& geometricalPhase,
9-
const int radarBlockWidth, const int radarBlockLength,
10-
const int azimuthFirstLine, const int rangeFirstPixel,
11-
const isce3::core::Interpolator<std::complex<float>>* interp)
7+
const std::valarray<double>& radarY, const int azimuthFirstLine,
8+
const int rangeFirstPixel,
9+
const isce3::core::Interpolator<std::complex<float>>* interp,
10+
const isce3::product::RadarGridParameters& radarGrid,
11+
const isce3::core::LUT2d<double>& dopplerLUT, const bool& flatten)
1212
{
13-
14-
size_t length = geoDataBlock.length();
15-
size_t width = geoDataBlock.width();
16-
int extraMargin = isce3::core::SINC_HALF;
13+
const int chipSize = isce3::core::SINC_ONE;
14+
const int width = geoDataBlock.width();
15+
const int length = geoDataBlock.length();
16+
const int inWidth = rdrDataBlock.width();
17+
const int inLength = rdrDataBlock.length();
18+
const int chipHalf = chipSize / 2;
1719

1820
#pragma omp parallel for
1921
for (size_t kk = 0; kk < length * width; ++kk) {
20-
2122
size_t i = kk / width;
2223
size_t j = kk % width;
2324

@@ -26,22 +27,65 @@ void isce3::geocode::interpolate(
2627
double rdrY = radarY[i * width + j] - azimuthFirstLine;
2728
double rdrX = radarX[i * width + j] - rangeFirstPixel;
2829

29-
if (rdrX < extraMargin || rdrY < extraMargin ||
30-
rdrX >= (radarBlockWidth - extraMargin) ||
31-
rdrY >= (radarBlockLength - extraMargin)) {
30+
const int intX = static_cast<int>(rdrX);
31+
const int intY = static_cast<int>(rdrY);
32+
const double fracX = rdrX - intX;
33+
const double fracY = rdrY - intY;
34+
35+
if ((intX < chipHalf) || (intX >= (inWidth - chipHalf)))
36+
continue;
37+
if ((intY < chipHalf) || (intY >= (inLength - chipHalf)))
38+
continue;
39+
40+
// Slant Range at the current output pixel
41+
const double rng =
42+
radarGrid.startingRange() +
43+
radarX[i * width + j] * radarGrid.rangePixelSpacing();
44+
45+
// Azimuth time at the current output pixel
46+
const double az = radarGrid.sensingStart() +
47+
radarY[i * width + j] / radarGrid.prf();
3248

33-
geoDataBlock(i,j) = std::complex<float> (0.0, 0.0);
49+
if (not dopplerLUT.contains(az, rng))
50+
continue;
3451

35-
} else {
52+
// Evaluate Doppler at current range and azimuth time
53+
const double dop =
54+
dopplerLUT.eval(az, rng) * 2 * M_PI / radarGrid.prf();
3655

37-
// Interpolate chip
38-
const std::complex<double> cval =
39-
interp->interpolate(rdrX, rdrY, rdrDataBlock);
56+
// Doppler to be added back. Simultaneously evaluate carrier
57+
// that needs to be added back after interpolation
58+
double carrier_phase = (dop * fracY); // +
59+
//_rgCarrier.eval(rdrX, rdrY) +
60+
//_azCarrier.eval(rdrX, rdrY);
4061

41-
// geometricalPhase is the sum of carrier (Doppler) phase to be added
42-
// back and the geometrical phase to be removed: exp(1J* (carrier
43-
// - 4.0*PI*slantRange/wavelength))
44-
geoDataBlock(i, j) = cval * geometricalPhase[i * width + j];
62+
if (flatten) {
63+
carrier_phase += (4.0 * (M_PI / radarGrid.wavelength())) * rng;
4564
}
46-
} // end for
65+
66+
isce3::core::Matrix<std::complex<float>> chip(chipSize, chipSize);
67+
// Read data chip without the carrier phases
68+
for (int ii = 0; ii < chipSize; ++ii) {
69+
// Row to read from
70+
const int chipRow = intY + ii - chipHalf;
71+
72+
// Carrier phase
73+
const double phase = dop * (ii - chipHalf);
74+
const std::complex<float> cval(std::cos(phase), -std::sin(phase));
75+
76+
// Set the data values after removing doppler in azimuth
77+
for (int jj = 0; jj < chipSize; ++jj) {
78+
// Column to read from
79+
const int chipCol = intX + jj - chipHalf;
80+
chip(ii, jj) = rdrDataBlock(chipRow, chipCol) * cval;
81+
}
82+
}
83+
// Interpolate chip
84+
const std::complex<float> cval =
85+
interp->interpolate(isce3::core::SINC_HALF + fracX,
86+
isce3::core::SINC_HALF + fracY, chip);
87+
88+
geoDataBlock(i, j) = cval * std::complex<float>(std::cos(carrier_phase),
89+
std::sin(carrier_phase));
90+
}
4791
}

‎cxx/isce3/geocode/interpolate.h‎

Lines changed: 15 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,8 @@
66

77
#include <isce3/core/Interpolator.h>
88
#include <isce3/core/Matrix.h>
9+
#include <isce3/product/Product.h>
10+
#include <isce3/product/RadarGridParameters.h>
911

1012
namespace isce3 { namespace geocode {
1113

@@ -14,20 +16,21 @@ namespace isce3 { namespace geocode {
1416
* @param[out] geoDataBlock a block of data in geo coordinates
1517
* @param[in] radarX the radar-coordinates x-index of the pixels in geo-grid
1618
* @param[in] radarY the radar-coordinates y-index of the pixels in geo-grid
17-
* @param[in] geometricalPhase the geometrical phase of each pixel in geo-grid
18-
* @param[in] radarBlockWidth width of the data block in radar coordinates
19-
* @param[in] radarBlockLength length of the data block in radar coordinates
20-
* @param[in] azimuthFirstLine azimuth time of the first sample
21-
* @param[in] rangeFirstPixel range of the first sample
19+
* @param[in] azimuthFirstLine line index of the first sample of the block
20+
* @param[in] rangeFirstPixel pixel index of the first sample of the block
2221
* @param[in] interp interpolator object
22+
* @param[in] radarGrid RadarGridParameters
23+
* @param[in] dopplerLUT Doppler LUT
24+
* @param[in] flatten flag to determine if the geocoded SLC will be flattened or
25+
* not
2326
*/
2427
void interpolate(const isce3::core::Matrix<std::complex<float>>& rdrDataBlock,
25-
isce3::core::Matrix<std::complex<float>>& geoDataBlock,
26-
const std::valarray<double>& radarX,
27-
const std::valarray<double>& radarY,
28-
const std::valarray<std::complex<double>>& geometricalPhase,
29-
const int radarBlockWidth, const int radarBlockLength,
30-
const int azimuthFirstLine, const int rangeFirstPixel,
31-
const isce3::core::Interpolator<std::complex<float>>* interp);
28+
isce3::core::Matrix<std::complex<float>>& geoDataBlock,
29+
const std::valarray<double>& radarX,
30+
const std::valarray<double>& radarY, const int azimuthFirstLine,
31+
const int rangeFirstPixel,
32+
const isce3::core::Interpolator<std::complex<float>>* interp,
33+
const isce3::product::RadarGridParameters& radarGrid,
34+
const isce3::core::LUT2d<double>& dopplerLUT, const bool& flatten);
3235

3336
}} // namespace isce3::geocode

0 commit comments

Comments
 (0)