cuDenseOffsets: add prepareGeometry()
PyCuAmpcor/examples/cuDenseOffsets: + add prepareGeometry() to generate the cooresponding geometry file + update file path in xml file to support dir change + fix str.decode bug with try/except statement + printout used time + check --redo before initiate PyCuAmpcor to avoid the long GPU memory error message PyCuAmpcor/README.md: fix typos for ROIPACLT1AB
Zhang Yunjun
Zhang Yunjun
parent
d6ba085ef5
commit
d81a98094e
|
|
@ -24,7 +24,7 @@ In practice, we
|
|||
|
||||
A detailed formulation can be found, e.g., by J. P. Lewis with [the frequency domain approach](http://scribblethink.org/Work/nvisionInterface/nip.html).
|
||||
|
||||
PyCuAmpcor follows the same procedure as the FORTRAN code, ampcor.F, in RIOPAC. In order to optimize the performance on GPU, some implementations are slightly different. In the [list the procedures](#5-list-of-procedures), we show the detailed steps of PyCuAmpcor, as well as their differences.
|
||||
PyCuAmpcor follows the same procedure as the FORTRAN code, ampcor.F, in ROIPAC. In order to optimize the performance on GPU, some implementations are slightly different. In the [list the procedures](#5-list-of-procedures), we show the detailed steps of PyCuAmpcor, as well as their differences.
|
||||
|
||||
## 2. Installation
|
||||
|
||||
|
|
@ -338,8 +338,8 @@ This step provides an initial estimate of the offset, usually with a large searc
|
|||
|
||||
**Difference to ROIPAC**
|
||||
|
||||
* RIOPAC only offers the time-domain algorithm. The frequency-domain algorithm is faster and is set as default in PyCuAmpcor.
|
||||
* RIOPAC proceeds from here only for windows with *good* match, or with high coherence. To maintain parallelism, PyCuAmpcor proceeds anyway while leaving the *filtering* to users in post processing.
|
||||
* ROIPAC only offers the time-domain algorithm. The frequency-domain algorithm is faster and is set as default in PyCuAmpcor.
|
||||
* ROIPAC proceeds from here only for windows with *good* match, or with high coherence. To maintain parallelism, PyCuAmpcor proceeds anyway while leaving the *filtering* to users in post processing.
|
||||
|
||||
|
||||
### 5.3 Extract a smaller window from the secondary window for oversampling
|
||||
|
|
@ -354,12 +354,12 @@ This step provides an initial estimate of the offset, usually with a large searc
|
|||
|
||||
**Difference to ROIPAC**
|
||||
|
||||
RIOPAC extracts the secondary window centering at the correlation surface peak. If the peak locates near the edge, zeros are padded if the extraction zone exceeds the window range. In PyCuAmpcor, the extraction center may be shifted away from peak to warrant all pixels being in the range of the original window.
|
||||
ROIPAC extracts the secondary window centering at the correlation surface peak. If the peak locates near the edge, zeros are padded if the extraction zone exceeds the window range. In PyCuAmpcor, the extraction center may be shifted away from peak to warrant all pixels being in the range of the original window.
|
||||
|
||||
|
||||
### 5.4 Oversampling reference and (extracted) secondary windows
|
||||
|
||||
* Oversample both the reference and the (extracted) secondary windows by a factor of 2, which is to avoid aliasing in the complex multiplication of the SAR images. The oversampling is performed with FFT (zero padding), same as in RIOPAC.
|
||||
* Oversample both the reference and the (extracted) secondary windows by a factor of 2, which is to avoid aliasing in the complex multiplication of the SAR images. The oversampling is performed with FFT (zero padding), same as in ROIPAC.
|
||||
* A deramping procedure is in general required for complex signals before oversampling, to shift the band center to 0. The procedure is only designed to remove a linear phase ramp. It doesn't work for TOPSAR, whose ramp goes quadratic. Instead, the amplitudes are taken before oversampling.
|
||||
* the amplitudes (real) are then taken for each pixel of the complex signals in reference and secondary windows.
|
||||
|
||||
|
|
@ -372,9 +372,9 @@ RIOPAC extracts the secondary window centering at the correlation surface peak.
|
|||
|
||||
**Difference to ROIPAC**
|
||||
|
||||
RIOPAC enlarges both windows to a size which is a power of 2; ideal for FFT. PyCuAmpcor uses their original sizes for FFT.
|
||||
ROIPAC enlarges both windows to a size which is a power of 2; ideal for FFT. PyCuAmpcor uses their original sizes for FFT.
|
||||
|
||||
RIOPAC always performs deramping with Method 1, to obtain the ramp by averaging the phase difference between neighboring pixels. For TOPS mode, users need to specify 'mag' as the image *datatype* such that the amplitudes are taken before oversampling. Therefore, deramping has no effect. In PyCuAmpcor, derampMethod=0 is equivalent to *datatype='mag'*, taking amplitudes but skipping deramping. derampMethod=1 always performs deramping, no matter the 'complex' or 'real' image datatypes.
|
||||
ROIPAC always performs deramping with Method 1, to obtain the ramp by averaging the phase difference between neighboring pixels. For TOPS mode, users need to specify 'mag' as the image *datatype* such that the amplitudes are taken before oversampling. Therefore, deramping has no effect. In PyCuAmpcor, derampMethod=0 is equivalent to *datatype='mag'*, taking amplitudes but skipping deramping. derampMethod=1 always performs deramping, no matter the 'complex' or 'real' image datatypes.
|
||||
|
||||
### 5.5 Cross-Correlate the oversampled reference and secondary windows
|
||||
|
||||
|
|
@ -386,11 +386,11 @@ RIOPAC always performs deramping with Method 1, to obtain the ramp by averaging
|
|||
|
||||
| PyCuAmpcor | CUDA variable | ampcor.F equivalent | Notes |
|
||||
| :--- | :--- | :---- | :--- |
|
||||
| corrSurfaceZoomInWindow | zoomWindowSize | i_cw | The size of correlation surface of the (anti-aliasing) oversampled reference/secondary windows, also used to set halfZoomWindowSizeRaw. Set it to 16 to be consistent with RIOPAC. |
|
||||
| corrSurfaceZoomInWindow | zoomWindowSize | i_cw | The size of correlation surface of the (anti-aliasing) oversampled reference/secondary windows, also used to set halfZoomWindowSizeRaw. Set it to 16 to be consistent with ROIPAC. |
|
||||
|
||||
**Difference to ROIPAC**
|
||||
|
||||
In RIOPAC, an extra resizing step is performed on the correlation surface, from (2\*halfZoomWindowSizeRaw\*rawDataOversamplingFactor+1, 2\*halfZoomWindowSizeRaw\*rawDataOversamplingFactor+1) to (i_cw, i_cw), centered at the peak (in RIOPAC, the peak seeking is incorporated in the correlation module while is seperate in PyCuAmpcor). i_cw is a user configurable variable; it could be smaller or bigger than 2\*i_srchp\*i_ovs+1=17 (fixed), leading to extraction or enlargement by padding 0s. This procedure is not performed in PyCuAmpcor, as it makes little difference in the next oversampling procedure.
|
||||
In ROIPAC, an extra resizing step is performed on the correlation surface, from (2\*halfZoomWindowSizeRaw\*rawDataOversamplingFactor+1, 2\*halfZoomWindowSizeRaw\*rawDataOversamplingFactor+1) to (i_cw, i_cw), centered at the peak (in ROIPAC, the peak seeking is incorporated in the correlation module while is seperate in PyCuAmpcor). i_cw is a user configurable variable; it could be smaller or bigger than 2\*i_srchp\*i_ovs+1=17 (fixed), leading to extraction or enlargement by padding 0s. This procedure is not performed in PyCuAmpcor, as it makes little difference in the next oversampling procedure.
|
||||
|
||||
### 5.6 Oversample the correlation surface and find the peak position
|
||||
|
||||
|
|
@ -412,4 +412,4 @@ Note that this offset does not include the pre-defined gross offset. Users need
|
|||
|
||||
**Difference to ROIPAC**
|
||||
|
||||
RIOPAC by default uses the sinc interpolator (the FFT method is included but one needs to change the FORTRAN code to switch). For since interpolator, there is no difference in implementations. For FFT, RIOPAC always enlarges the window to a size in power of 2.
|
||||
ROIPAC by default uses the sinc interpolator (the FFT method is included but one needs to change the FORTRAN code to switch). For since interpolator, there is no difference in implementations. For FFT, ROIPAC always enlarges the window to a size in power of 2.
|
||||
|
|
|
|||
|
|
@ -4,19 +4,25 @@
|
|||
|
||||
|
||||
import os
|
||||
import sys
|
||||
import time
|
||||
import argparse
|
||||
import numpy as np
|
||||
from osgeo import gdal
|
||||
|
||||
import isce
|
||||
import isceobj
|
||||
from isceobj.Util.decorators import use_api
|
||||
from isceobj.Util.ImageUtil import ImageLib as IML
|
||||
from contrib.PyCuAmpcor.PyCuAmpcor import PyCuAmpcor
|
||||
|
||||
|
||||
EXAMPLE = '''example
|
||||
cuDenseOffsets.py -r ./merged/SLC/20151120/20151120.slc.full -s ./merged/SLC/20151214/20151214.slc.full
|
||||
--outprefix ./merged/offsets/20151120_20151214/offset
|
||||
--ww 256 --wh 256 --oo 32 --kw 300 --kh 100 --nwac 100 --nwdc 1 --sw 8 --sh 8 --gpuid 2
|
||||
cuDenseOffsets.py -r ./SLC/20151120/20151120.slc.full -s ./SLC/20151214/20151214.slc.full
|
||||
cuDenseOffsets.py -r ./SLC/20151120/20151120.slc.full -s ./SLC/20151214/20151214.slc.full --outprefix ./offsets/20151120_20151214/offset --ww 256 --wh 256 --sw 8 --sh 8 --oo 32 --kw 300 --kh 100 --nwac 100 --nwdc 1 --gpuid 2
|
||||
|
||||
# offset and its geometry
|
||||
cuDenseOffsets.py -r ./SLC/20151120/20151120.slc.full -s ./SLC/20151214/20151214.slc.full --outprefix ./offsets/20151120_20151214/offset --ww 256 --wh 256 --sw 8 --sh 8 --oo 32 --kw 300 --kh 100 --nwac 100 --nwdc 1 --gpuid 2 --full-geom ./geom_reference --out-geom ./offset/geom_reference
|
||||
'''
|
||||
|
||||
|
||||
|
|
@ -25,8 +31,7 @@ def createParser():
|
|||
Command line parser.
|
||||
'''
|
||||
|
||||
|
||||
parser = argparse.ArgumentParser(description='Generate offset field between two Sentinel slc',
|
||||
parser = argparse.ArgumentParser(description='Generate offset field between two SLCs',
|
||||
formatter_class=argparse.RawTextHelpFormatter,
|
||||
epilog=EXAMPLE)
|
||||
|
||||
|
|
@ -38,9 +43,9 @@ def createParser():
|
|||
|
||||
parser.add_argument('--op','--outprefix','--output-prefix', type=str, dest='outprefix',
|
||||
default='offset', required=True,
|
||||
help='Output prefix, default: offset.')
|
||||
help='Output prefix (default: %(default)s).')
|
||||
parser.add_argument('--os','--outsuffix', type=str, dest='outsuffix', default='',
|
||||
help='Output suffix, default:.')
|
||||
help='Output suffix (default: %(default)s).')
|
||||
|
||||
# window size settings
|
||||
parser.add_argument('--ww', type=int, dest='winwidth', default=64,
|
||||
|
|
@ -66,9 +71,11 @@ def createParser():
|
|||
parser.add_argument('--nwd', type=int, dest='numWinDown', default=-1,
|
||||
help='Number of window down (default: %(default)s).')
|
||||
parser.add_argument('--startpixelac', dest='startpixelac', type=int, default=-1,
|
||||
help='Starting Pixel across of the reference image(default: %(default)s to be determined by margin and search range).')
|
||||
help='Starting Pixel across of the reference image.' +
|
||||
'Default: %(default)s to be determined by margin and search range.')
|
||||
parser.add_argument('--startpixeldw', dest='startpixeldw', type=int, default=-1,
|
||||
help='Starting Pixel down of the reference image (default: %(default)s).')
|
||||
help='Starting Pixel down of the reference image.' +
|
||||
'Default: %(default)s to be determined by margin and search range.')
|
||||
|
||||
# cross-correlation algorithm
|
||||
parser.add_argument('--alg', '--algorithm', dest='algorithm', type=int, default=0,
|
||||
|
|
@ -103,6 +110,12 @@ def createParser():
|
|||
corr.add_argument('--corr-osm', '--corr-over-samp-method', type=int, dest='corr_oversamplemethod', default=0,
|
||||
help = 'Oversampling method for the correlation surface 0=fft, 1=sinc (default: %(default)s).')
|
||||
|
||||
geom = parser.add_argument_group('Geometry', 'generate corresponding geometry datasets ')
|
||||
geom.add_argument('--full-geom', dest='full_geometry_dir', type=str,
|
||||
help='(Input) Directory of geometry files in full resolution.')
|
||||
geom.add_argument('--out-geom', dest='out_geometry_dir', type=str,
|
||||
help='(Output) Directory of geometry files corresponding to the offset field.')
|
||||
|
||||
# gpu settings
|
||||
proc = parser.add_argument_group('Processing parameters')
|
||||
proc.add_argument('--gpuid', '--gid', '--gpu-id', dest='gpuid', type=int, default=0,
|
||||
|
|
@ -126,7 +139,7 @@ def createParser():
|
|||
|
||||
def cmdLineParse(iargs = None):
|
||||
parser = createParser()
|
||||
inps = parser.parse_args(args=iargs)
|
||||
inps = parser.parse_args(args=iargs)
|
||||
|
||||
# check oversampled window size
|
||||
if (inps.winwidth + 2 * inps.srcwidth ) * inps.raw_oversample > 1024:
|
||||
|
|
@ -144,6 +157,21 @@ def cmdLineParse(iargs = None):
|
|||
@use_api
|
||||
def estimateOffsetField(reference, secondary, inps=None):
|
||||
|
||||
# check redo
|
||||
print('redo: ', inps.redo)
|
||||
if not inps.redo:
|
||||
offsetImageName = '{}{}.bip'.format(inps.outprefix, inps.outsuffix)
|
||||
if os.path.exists(offsetImageName):
|
||||
print('offset field file: {} exists and w/o redo, skip re-estimation.'.format(offsetImageName))
|
||||
return 0
|
||||
|
||||
# update file path in xml file
|
||||
for fname in [reference, secondary]:
|
||||
fname = os.path.abspath(fname)
|
||||
img = IML.loadImage(fname)[0]
|
||||
img.filename = fname
|
||||
img.setAccessMode('READ')
|
||||
img.renderHdr()
|
||||
|
||||
###Loading the secondary image object
|
||||
sim = isceobj.createSlcImage()
|
||||
|
|
@ -170,7 +198,6 @@ def estimateOffsetField(reference, secondary, inps=None):
|
|||
objOffset.derampMethod = inps.deramp
|
||||
print('deramp method (0 for magnitude, 1 for complex): ', objOffset.derampMethod)
|
||||
|
||||
|
||||
objOffset.referenceImageName = reference+'.vrt'
|
||||
objOffset.referenceImageHeight = length
|
||||
objOffset.referenceImageWidth = width
|
||||
|
|
@ -231,10 +258,11 @@ def estimateOffsetField(reference, secondary, inps=None):
|
|||
print('correlation surface oversampling factor:', inps.corr_oversample)
|
||||
|
||||
# output filenames
|
||||
objOffset.offsetImageName = str(inps.outprefix) + str(inps.outsuffix) + '.bip'
|
||||
objOffset.grossOffsetImageName = str(inps.outprefix) + str(inps.outsuffix) + '_gross.bip'
|
||||
objOffset.snrImageName = str(inps.outprefix) + str(inps.outsuffix) + '_snr.bip'
|
||||
objOffset.covImageName = str(inps.outprefix) + str(inps.outsuffix) + '_cov.bip'
|
||||
fbase = '{}{}'.format(inps.outprefix, inps.outsuffix)
|
||||
objOffset.offsetImageName = fbase + '.bip'
|
||||
objOffset.grossOffsetImageName = fbase + '_gross.bip'
|
||||
objOffset.snrImageName = fbase + '_snr.bip'
|
||||
objOffset.covImageName = fbase + '_cov.bip'
|
||||
print("offsetfield: ",objOffset.offsetImageName)
|
||||
print("gross offsetfield: ",objOffset.grossOffsetImageName)
|
||||
print("snr: ",objOffset.snrImageName)
|
||||
|
|
@ -243,14 +271,16 @@ def estimateOffsetField(reference, secondary, inps=None):
|
|||
# whether to include the gross offset in offsetImage
|
||||
objOffset.mergeGrossOffset = inps.merge_gross_offset
|
||||
|
||||
offsetImageName = objOffset.offsetImageName
|
||||
grossOffsetImageName = objOffset.grossOffsetImageName
|
||||
snrImageName = objOffset.snrImageName
|
||||
covImageName = objOffset.covImageName
|
||||
|
||||
if os.path.exists(offsetImageName) and not inps.redo:
|
||||
print('offsetfield file {} exists while the redo flag is {}.'.format(offsetImageName, inps.redo))
|
||||
return 0
|
||||
try:
|
||||
offsetImageName = objOffset.offsetImageName.decode('utf8')
|
||||
grossOffsetImageName = objOffset.grossOffsetImageName.decode('utf8')
|
||||
snrImageName = objOffset.snrImageName.decode('utf8')
|
||||
covImageName = objOffset.covImageName.decode('utf8')
|
||||
except:
|
||||
offsetImageName = objOffset.offsetImageName
|
||||
grossOffsetImageName = objOffset.grossOffsetImageName
|
||||
snrImageName = objOffset.snrImageName
|
||||
covImageName = objOffset.covImageName
|
||||
|
||||
# generic control
|
||||
objOffset.numberWindowDownInChunk = inps.numWinDownInChunk
|
||||
|
|
@ -271,7 +301,8 @@ def estimateOffsetField(reference, secondary, inps=None):
|
|||
grossOffset = np.fromfile(inps.gross_offset_file, dtype=np.int32)
|
||||
numberWindows = objOffset.numberWindowDown*objOffset.numberWindowAcross
|
||||
if grossOffset.size != 2*numberWindows :
|
||||
print('The input gross offsets do not match the number of windows {} by {} in int32 type'.format(objOffset.numberWindowDown, objOffset.numberWindowAcross))
|
||||
print(('WARNING: The input gross offsets do not match the number of windows:'
|
||||
' {} by {} in int32 type').format(objOffset.numberWindowDown, objOffset.numberWindowAcross))
|
||||
return 0;
|
||||
grossOffset = grossOffset.reshape(numberWindows, 2)
|
||||
grossAzimuthOffset = grossOffset[:, 0]
|
||||
|
|
@ -338,22 +369,86 @@ def estimateOffsetField(reference, secondary, inps=None):
|
|||
covImg.setAccessMode('read')
|
||||
covImg.renderHdr()
|
||||
|
||||
return objOffset
|
||||
|
||||
|
||||
def prepareGeometry(full_dir, out_dir, x_start, y_start, x_step, y_step, num_win_x, num_win_y,
|
||||
fbases=['hgt','lat','lon','los','shadowMask','waterMask']):
|
||||
"""Generate multilooked geometry datasets in the same grid as the estimated offset field
|
||||
from the full resolution geometry datasets.
|
||||
Parameters: full_dir - str, path of input geometry directory in full resolution
|
||||
out_dir - str, path of output geometry directory
|
||||
x/y_start - int, starting column/row number
|
||||
x/y_step - int, output pixel step in column/row direction
|
||||
num_win_x/y - int, number of columns/rows
|
||||
"""
|
||||
print('-'*50)
|
||||
print('generate the corresponding multi-looked geometry datasets using gdal ...')
|
||||
in_files = [os.path.join(full_dir, '{}.rdr.full'.format(i)) for i in fbases]
|
||||
in_files = [i for i in in_files if os.path.isfile(i)]
|
||||
if len(in_files) == 0:
|
||||
raise ValueError('No full resolution geometry file found in: {}'.format(full_dir))
|
||||
|
||||
fbases = [os.path.basename(i).split('.')[0] for i in in_files]
|
||||
out_files = [os.path.join(out_dir, '{}.rdr'.format(i)) for i in fbases]
|
||||
os.makedirs(out_dir, exist_ok=True)
|
||||
|
||||
for i in range(len(in_files)):
|
||||
in_file = in_files[i]
|
||||
out_file = out_files[i]
|
||||
|
||||
# input file size
|
||||
ds = gdal.Open(in_file, gdal.GA_ReadOnly)
|
||||
in_wid = ds.RasterXSize
|
||||
in_len = ds.RasterYSize
|
||||
|
||||
# starting column/row and column/row number
|
||||
src_win = [x_start, y_start, num_win_x * x_step, num_win_y * y_step]
|
||||
print('read {} from file: {}'.format(src_win, in_file))
|
||||
|
||||
# write binary data file
|
||||
print('write file: {}'.format(out_file))
|
||||
opts = gdal.TranslateOptions(format='ENVI',
|
||||
width=num_win_x,
|
||||
height=num_win_y,
|
||||
srcWin=src_win,
|
||||
noData=0)
|
||||
gdal.Translate(out_file, ds, options=opts)
|
||||
ds = None
|
||||
|
||||
# write VRT file
|
||||
print('write file: {}'.format(out_file+'.vrt'))
|
||||
ds = gdal.Open(out_file, gdal.GA_ReadOnly)
|
||||
gdal.Translate(out_file+'.vrt', ds, options=gdal.TranslateOptions(format='VRT'))
|
||||
ds = None
|
||||
|
||||
return
|
||||
|
||||
|
||||
def main(iargs=None):
|
||||
|
||||
inps = cmdLineParse(iargs)
|
||||
outDir = os.path.dirname(inps.outprefix)
|
||||
print(inps.outprefix)
|
||||
start_time = time.time()
|
||||
|
||||
print(inps.outprefix)
|
||||
outDir = os.path.dirname(inps.outprefix)
|
||||
os.makedirs(outDir, exist_ok=True)
|
||||
|
||||
estimateOffsetField(inps.reference, inps.secondary, inps)
|
||||
# estimate offset
|
||||
objOffset = estimateOffsetField(inps.reference, inps.secondary, inps)
|
||||
|
||||
# generate geometry
|
||||
if inps.full_geometry_dir and inps.out_geometry_dir:
|
||||
prepareGeometry(inps.full_geometry_dir, inps.out_geometry_dir,
|
||||
x_start=objOffset.referenceStartPixelAcrossStatic,
|
||||
y_start=objOffset.referenceStartPixelDownStatic,
|
||||
x_step=objOffset.skipSampleAcross,
|
||||
y_step=objOffset.skipSampleDown,
|
||||
num_win_x=objOffset.numberWindowAcross,
|
||||
num_win_y=objOffset.numberWindowDown)
|
||||
|
||||
m, s = divmod(time.time() - start_time, 60)
|
||||
print('time used: {:02.0f} mins {:02.1f} secs.\n'.format(m, s))
|
||||
return
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
main()
|
||||
main(sys.argv[1:])
|
||||
|
|
|
|||
Loading…
Reference in New Issue