From 563568561d1393a1e8e0ebae21cc29a596b3d799 Mon Sep 17 00:00:00 2001 From: Lijun Zhu Date: Sat, 13 Feb 2021 12:56:58 -0800 Subject: [PATCH] change runDenseOffset(s) in Alos2Proc/TopsProc in accordance with updated PyCuAmpcor --- .../isceobj/Alos2Proc/runDenseOffset.py | 90 +++++++++---- .../isceobj/TopsProc/runDenseOffsets.py | 120 ++++++++++++------ contrib/PyCuAmpcor/README.md | 5 +- contrib/PyCuAmpcor/src/PyCuAmpcor.pyx | 9 +- contrib/PyCuAmpcor/src/cuAmpcorChunk.cu | 2 +- contrib/PyCuAmpcor/src/cuAmpcorController.cu | 7 +- contrib/PyCuAmpcor/src/cudaUtil.h | 4 +- 7 files changed, 158 insertions(+), 79 deletions(-) diff --git a/components/isceobj/Alos2Proc/runDenseOffset.py b/components/isceobj/Alos2Proc/runDenseOffset.py index 1c87b9c..e3d3853 100644 --- a/components/isceobj/Alos2Proc/runDenseOffset.py +++ b/components/isceobj/Alos2Proc/runDenseOffset.py @@ -70,7 +70,7 @@ def runDenseOffset(self): #get water body mask wbd0=np.memmap('wbd.rdr', dtype=np.int8, mode='r', shape=(length0, width0)) - wbd0=wbd0[0+self._insar.offsetImageTopoffset:length0:self.offsetSkipHeight, + wbd0=wbd0[0+self._insar.offsetImageTopoffset:length0:self.offsetSkipHeight, 0+self._insar.offsetImageLeftoffset:width0:self.offsetSkipWidth] wbd = np.zeros((length+100, width+100), dtype=np.int8) wbd[0:wbd0.shape[0], 0:wbd0.shape[1]]=wbd0 @@ -240,15 +240,11 @@ def runDenseOffsetGPU(self): print('dense offset skip width: %d' % (self.offsetSkipWidth)) print('dense offset skip hight: %d' % (self.offsetSkipHeight)) print('dense offset covariance surface oversample factor: %d' % (self.offsetCovarianceOversamplingFactor)) - print('dense offset covariance surface oversample window size: %d\n' % (self.offsetCovarianceOversamplingWindowsize)) objOffset = PyCuAmpcor.PyCuAmpcor() objOffset.algorithm = 0 - objOffset.deviceID = -1 - objOffset.nStreams = 2 - #original ampcor program in roi_pac uses phase gradient to deramp - objOffset.derampMethod = 2 + objOffset.derampMethod = 1 # 1=linear phase ramp, 0=take mag, 2=skip objOffset.referenceImageName = self._insar.referenceSlc objOffset.referenceImageHeight = m.length objOffset.referenceImageWidth = m.width @@ -256,56 +252,79 @@ def runDenseOffsetGPU(self): objOffset.secondaryImageHeight = s.length objOffset.secondaryImageWidth = s.width objOffset.offsetImageName = self._insar.denseOffset + objOffset.grossOffsetImageName = self._insar.denseOffset + ".gross" objOffset.snrImageName = self._insar.denseOffsetSnr + objOffset.covImageName = self._insar.denseOffsetCov objOffset.windowSizeWidth = self.offsetWindowWidth objOffset.windowSizeHeight = self.offsetWindowHeight - #objOffset.halfSearchRangeAcross = int(self.offsetSearchWindowWidth / 2 + 0.5) - #objOffset.halfSearchRangeDown = int(self.offsetSearchWindowHeight / 2 + 0.5) + objOffset.halfSearchRangeAcross = self.offsetSearchWindowWidth objOffset.halfSearchRangeDown = self.offsetSearchWindowHeight + objOffset.skipSampleDown = self.offsetSkipHeight objOffset.skipSampleAcross = self.offsetSkipWidth + #Oversampling method for correlation surface(0=fft,1=sinc) - objOffset.corrSufaceOverSamplingMethod = 0 + objOffset.corrSurfaceOverSamplingMethod = 0 objOffset.corrSurfaceOverSamplingFactor = self.offsetCovarianceOversamplingFactor - objOffset.corrSurfaceZoomInWindow = self.offsetCovarianceOversamplingWindowsize + + # set gross offset objOffset.grossOffsetAcrossStatic = 0 objOffset.grossOffsetDownStatic = 0 + # set the margin + margin = 0 - objOffset.referenceStartPixelDownStatic = self.offsetWindowHeight//2 - objOffset.referenceStartPixelAcrossStatic = self.offsetWindowWidth//2 + # adjust the margin + margin = max(margin, abs(objOffset.grossOffsetAcrossStatic), abs(objOffset.grossOffsetDownStatic)) - objOffset.numberWindowDown = (m.length - 2*self.offsetSearchWindowHeight - self.offsetWindowHeight) // self.offsetSkipHeight - objOffset.numberWindowAcross = (m.width - 2*self.offsetSearchWindowWidth - self.offsetWindowWidth) // self.offsetSkipWidth + # set the starting pixel of the first reference window + objOffset.referenceStartPixelDownStatic = margin + self.offsetSearchWindowHeight + objOffset.referenceStartPixelAcrossStatic = margin + self.offsetSearchWindowWidth - # generic control - objOffset.numberWindowDownInChunk = 8 - objOffset.numberWindowAcrossInChunk = 8 + # find out the total number of windows + objOffset.numberWindowDown = (m.length - 2*margin - 2*self.offsetSearchWindowHeight - self.offsetWindowHeight) // self.offsetSkipHeight + objOffset.numberWindowAcross = (m.width - 2*margin - 2*self.offsetSearchWindowWidth - self.offsetWindowWidth) // self.offsetSkipWidth + + # gpu job control + objOffset.deviceID = 0 + objOffset.nStreams = 2 + objOffset.numberWindowDownInChunk = 1 + objOffset.numberWindowAcrossInChunk = 64 objOffset.mmapSize = 16 + # pass/adjust the parameters objOffset.setupParams() - objOffset.setConstantGrossOffset(0, 0) + # set up the starting pixels for each window, based on the gross offset + objOffset.setConstantGrossOffset(objOffset.grossOffsetAcrossStatic, objOffset.grossOffsetDownStatic) + # check whether all pixels are in image range (optional) objOffset.checkPixelInImageRange() + print('\n======================================') + print('Running PyCuAmpcor...') + print('======================================\n') objOffset.runAmpcor() ### Store params for later - self._insar.offsetImageTopoffset = objOffset.halfSearchRangeDown - self._insar.offsetImageLeftoffset = objOffset.halfSearchRangeAcross + # location of the center of the first reference window + self._insar.offsetImageTopoffset = objOffset.referenceStartPixelDownStatic + (objOffset.windowSizeHeight-1)//2 + self._insar.offsetImageLeftoffset = objOffset.referenceStartPixelAcrossStatic +(objOffset.windowSizeWidth-1)//2 - + # offset image dimension, the number of windows width = objOffset.numberWindowAcross length = objOffset.numberWindowDown - offsetBIP = np.fromfile(objOffset.offsetImageName.decode('utf-8'), dtype=np.float32).reshape(length, width*2) + + # convert the offset image from BIP to BIL + offsetBIP = np.fromfile(objOffset.offsetImageName, dtype=np.float32).reshape(length, width*2) offsetBIL = np.zeros((length*2, width), dtype=np.float32) offsetBIL[0:length*2:2, :] = offsetBIP[:, 1:width*2:2] offsetBIL[1:length*2:2, :] = offsetBIP[:, 0:width*2:2] - os.remove(objOffset.offsetImageName.decode('utf-8')) - offsetBIL.astype(np.float32).tofile(objOffset.offsetImageName.decode('utf-8')) + os.remove(objOffset.offsetImageName) + offsetBIL.astype(np.float32).tofile(objOffset.offsetImageName) + # generate offset image description files outImg = isceobj.createImage() outImg.setDataType('FLOAT') - outImg.setFilename(objOffset.offsetImageName.decode('utf-8')) + outImg.setFilename(objOffset.offsetImageName) outImg.setBands(2) outImg.scheme = 'BIL' outImg.setWidth(objOffset.numberWindowAcross) @@ -314,8 +333,11 @@ def runDenseOffsetGPU(self): outImg.setAccessMode('read') outImg.renderHdr() + # gross offset image is not needed, since all zeros + + # generate snr image description files snrImg = isceobj.createImage() - snrImg.setFilename( objOffset.snrImageName.decode('utf8')) + snrImg.setFilename( objOffset.snrImageName) snrImg.setDataType('FLOAT') snrImg.setBands(1) snrImg.setWidth(objOffset.numberWindowAcross) @@ -323,4 +345,20 @@ def runDenseOffsetGPU(self): snrImg.setAccessMode('read') snrImg.renderHdr() + # generate cov image description files + # covariance of azimuth/range offsets. + # 1st band: cov(az, az), 2nd band: cov(rg, rg), 3rd band: cov(az, rg) + covImg = isceobj.createImage() + covImg.setFilename(objOffset.covImageName) + covImg.setDataType('FLOAT') + covImg.setBands(3) + covImg.scheme = 'BIP' + covImg.setWidth(objOffset.numberWindowAcross) + covImg.setLength(objOffset.numberWindowDown) + outImg.addDescription('covariance of azimuth/range offsets') + covImg.setAccessMode('read') + covImg.renderHdr() + return (objOffset.numberWindowAcross, objOffset.numberWindowDown) + +# end of file diff --git a/components/isceobj/TopsProc/runDenseOffsets.py b/components/isceobj/TopsProc/runDenseOffsets.py index 75023cf..a1ffb29 100644 --- a/components/isceobj/TopsProc/runDenseOffsets.py +++ b/components/isceobj/TopsProc/runDenseOffsets.py @@ -156,11 +156,11 @@ def runDenseOffsetsGPU(self): secondary = os.path.join(self._insar.mergedDirname, sf) ####For this module currently, we need to create an actual file on disk - for infile in [reference,secondary]: if os.path.isfile(infile): continue + print('Creating actual file {} ...\n'.format(infile)) cmd = 'gdal_translate -of ENVI {0}.vrt {0}'.format(infile) status = os.system(cmd) if status: @@ -170,21 +170,23 @@ def runDenseOffsetsGPU(self): m = isceobj.createSlcImage() m.load(reference + '.xml') m.setAccessMode('READ') -# m.createImage() ### Load the secondary object s = isceobj.createSlcImage() s.load(secondary + '.xml') s.setAccessMode('READ') -# s.createImage() + # get the dimension width = m.getWidth() length = m.getLength() + ### create the GPU processor objOffset = PyCuAmpcor.PyCuAmpcor() + + ### Set parameters + # cross-correlation method, 0=Frequency domain, 1= Time domain objOffset.algorithm = 0 - objOffset.deviceID = -1 - objOffset.nStreams = 2 + # deramping method: 0 to take magnitude (fixed for Tops) objOffset.derampMethod = 0 objOffset.referenceImageName = reference + '.vrt' objOffset.referenceImageHeight = length @@ -193,36 +195,56 @@ def runDenseOffsetsGPU(self): objOffset.secondaryImageHeight = length objOffset.secondaryImageWidth = width - objOffset.numberWindowDown = (length-100-self.winhgt)//self.skiphgt - objOffset.numberWindowAcross = (width-100-self.winwidth)//self.skipwidth + # adjust the margin + margin = max(self.margin, abs(self.azshift), abs(self.rgshift)) + # set the start pixel in the reference image + objOffset.referenceStartPixelDownStatic = margin + self.srchgt + objOffset.referenceStartPixelAcrossStatic = margin + self.srcwidth + + # compute the number of windows + objOffset.numberWindowDown = (length-2*margin-2*self.srchgt-self.winhgt)//self.skiphgt + objOffset.numberWindowAcross = (width-2*margin-2*self.srcwidth-self.winwidth)//self.skipwidth + + # set the template window size objOffset.windowSizeHeight = self.winhgt objOffset.windowSizeWidth = self.winwidth + # set the (half) search range objOffset.halfSearchRangeDown = self.srchgt objOffset.halfSearchRangeAcross = self.srcwidth - objOffset.referenceStartPixelDownStatic = 50 - objOffset.referenceStartPixelAcrossStatic = 50 - + # set the skip distance between windows objOffset.skipSampleDown = self.skiphgt objOffset.skipSampleAcross = self.skipwidth - objOffset.corrSufaceOverSamplingMethod = 0 + # correlation surface oversampling method, # 0=FFT, 1=Sinc + objOffset.corrSurfaceOverSamplingMethod = 0 + # oversampling factor objOffset.corrSurfaceOverSamplingFactor = self.oversample - # generic control - objOffset.numberWindowDownInChunk = 10 - objOffset.numberWindowAcrossInChunk = 10 + ### gpu control + objOffset.deviceID = 0 + objOffset.nStreams = 2 + # number of windows in a chunk/batch + objOffset.numberWindowDownInChunk = 1 + objOffset.numberWindowAcrossInChunk = 64 + # memory map cache size in GB objOffset.mmapSize = 16 + # Modify BIL in filename to BIP if needed and store for future use + prefix, ext = os.path.splitext(self._insar.offsetfile) + if ext == '.bil': + ext = '.bip' + self._insar.offsetfile = prefix + ext - objOffset.setupParams() + # set the output file name + objOffset.offsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile) + objOffset.grossOffsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile + ".gross") + objOffset.snrImageName = os.path.join(self._insar.mergedDirname, self._insar.snrfile) + objOffset.covImageName = os.path.join(self._insar.mergedDirname, self._insar.covfile) - objOffset.setConstantGrossOffset(self.azshift,self.rgshift) - -# objOffset.numberThreads = 1 - ### Configure dense Ampcor object + ### print the settings print('\nReference frame: %s' % (mf)) print('Secondary frame: %s' % (sf)) print('Main window size width: %d' % (self.winwidth)) @@ -231,41 +253,41 @@ def runDenseOffsetsGPU(self): print('Search window size height: %d' % (self.srchgt)) print('Skip sample across: %d' % (self.skipwidth)) print('Skip sample down: %d' % (self.skiphgt)) - print('Field margin: %d' % (self.margin)) + print('Field margin: %d' % (margin)) print('Oversampling factor: %d' % (self.oversample)) print('Gross offset across: %d' % (self.rgshift)) print('Gross offset down: %d\n' % (self.azshift)) - - #Modify BIL in filename to BIP if needed and store for future use - prefix, ext = os.path.splitext(self._insar.offsetfile) - if ext == '.bil': - ext = '.bip' - self._insar.offsetfile = prefix + ext - - objOffset.offsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile) - objOffset.snrImageName = os.path.join(self._insar.mergedDirname, self._insar.snrfile) - print('Output dense offsets file name: %s' % (objOffset.offsetImageName)) + print('Output gross offsets file name: %s' % (objOffset.grossOffsetImageName)) print('Output SNR file name: %s' % (objOffset.snrImageName)) + print('Output COV file name: %s' % (objOffset.covImageName)) + + # pass the parameters to C++ programs + objOffset.setupParams() + # set the (static) gross offset + objOffset.setConstantGrossOffset(self.azshift,self.rgshift) + # make sure all pixels are in range + objOffset.checkPixelInImageRange() + print('\n======================================') - print('Running dense ampcor...') + print('Running PyCuAmpcor...') print('======================================\n') - - objOffset.checkPixelInImageRange() + # run ampcor objOffset.runAmpcor() - #objOffset.denseampcor(m, s) ### Where the magic happens... - ### Store params for later + # offset width x length, also number of windows self._insar.offset_width = objOffset.numberWindowAcross self._insar.offset_length = objOffset.numberWindowDown - self._insar.offset_top = 50 - self._insar.offset_left = 50 + # the center of the first reference window + self._insar.offset_top = objOffset.referenceStartPixelDownStatic + (objOffset.windowSizeHeight-1)//2 + self._insar.offset_left = objOffset.referenceStartPixelAcrossStatic + (objOffset.windowSizeWidth-1)//2 + # generate description files for output images outImg = isceobj.createImage() outImg.setDataType('FLOAT') - outImg.setFilename(objOffset.offsetImageName.decode('utf-8')) + outImg.setFilename(objOffset.offsetImageName) outImg.setBands(2) outImg.scheme = 'BIP' outImg.setWidth(objOffset.numberWindowAcross) @@ -273,8 +295,19 @@ def runDenseOffsetsGPU(self): outImg.setAccessMode('read') outImg.renderHdr() + # gross offset + goutImg = isceobj.createImage() + goutImg.setDataType('FLOAT') + goutImg.setFilename(objOffset.grossOffsetImageName) + goutImg.setBands(2) + goutImg.scheme = 'BIP' + goutImg.setWidth(objOffset.numberWindowAcross) + goutImg.setLength(objOffset.numberWindowDown) + goutImg.setAccessMode('read') + goutImg.renderHdr() + snrImg = isceobj.createImage() - snrImg.setFilename( objOffset.snrImageName.decode('utf-8')) + snrImg.setFilename(objOffset.snrImageName) snrImg.setDataType('FLOAT') snrImg.setBands(1) snrImg.setWidth(objOffset.numberWindowAcross) @@ -282,6 +315,15 @@ def runDenseOffsetsGPU(self): snrImg.setAccessMode('read') snrImg.renderHdr() + covImg = isceobj.createImage() + covImg.setFilename(objOffset.covImageName) + covImg.setDataType('FLOAT') + covImg.setBands(3) + covImg.scheme = 'BIP' + covImg.setWidth(objOffset.numberWindowAcross) + covImg.setLength(objOffset.numberWindowDown) + covImg.setAccessMode('read') + covImg.renderHdr() if __name__ == '__main__' : diff --git a/contrib/PyCuAmpcor/README.md b/contrib/PyCuAmpcor/README.md index 6632ea1..62ce22e 100644 --- a/contrib/PyCuAmpcor/README.md +++ b/contrib/PyCuAmpcor/README.md @@ -91,7 +91,7 @@ mm=0 # margin to be neglected gross=0 # whether to use a varying gross offset azshift=0 # constant gross offset along height/azimuth rgshift=0 # constant gross offset along width/range -deramp=0 # 0 for mag (TOPS), 1 for complex +deramp=0 # 0 for mag (TOPS), 1 for complex linear ramp, 2 for complex no deramping oo=32 # correlation surface oversampling factor outprefix=./merged/20151120_20151214/offset # output prefix outsuffix=_ww64_wh64 # output suffix @@ -368,8 +368,7 @@ RIOPAC extracts the secondary window centering at the correlation surface peak. | PyCuAmpcor | CUDA variable | ampcor.F equivalent | Notes | | :--- | :--- | :---- | :--- | | rawDataOversamplingFactor | rawDataOversamplingFactor | i_ovs=2 | the oversampling factor for reference and secondary windows, use 2 for InSAR SLCs. | -| derampMethod | derampMethod | 1 or no effect on TOPS | 0=mag for TOPS, 1=deramping (default), else=skip deramping. - +| derampMethod | derampMethod | 1 or no effect on TOPS | Only for complex: 0=take mag (TOPS), 1=linear deramp (default), else=skip deramp. **Difference to ROIPAC** diff --git a/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx b/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx index 9db1c8f..7fc9cc5 100644 --- a/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx +++ b/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx @@ -259,12 +259,7 @@ cdef class PyCuAmpcor(object): @secondaryImageName.setter def secondaryImageName(self, str a): self.c_cuAmpcor.param.secondaryImageName = a.encode() - @property - def referenceImageName(self): - return self.c_cuAmpcor.param.referenceImageName - @referenceImageName.setter - def referenceImageName(self, str a): - self.c_cuAmpcor.param.referenceImageName = a.encode() + @property def referenceImageHeight(self): return self.c_cuAmpcor.param.referenceImageHeight @@ -449,4 +444,4 @@ cdef class PyCuAmpcor(object): self.c_cuAmpcor.runAmpcor() -# end of file \ No newline at end of file +# end of file diff --git a/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu b/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu index 286b099..c7c894d 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu +++ b/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu @@ -3,7 +3,7 @@ /** * Run ampcor process for a batch of images (a chunk) - * @param[in] idxDown_ index oIDIVUP(i,j) ((i+j-1)/j)f the chunk along Down/Azimuth direction + * @param[in] idxDown_ index of the chunk along Down/Azimuth direction * @param[in] idxAcross_ index of the chunk along Across/Range direction */ void cuAmpcorChunk::run(int idxDown_, int idxAcross_) diff --git a/contrib/PyCuAmpcor/src/cuAmpcorController.cu b/contrib/PyCuAmpcor/src/cuAmpcorController.cu index fc9a086..eee12fe 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorController.cu +++ b/contrib/PyCuAmpcor/src/cuAmpcorController.cu @@ -41,7 +41,9 @@ void cuAmpcorController::runAmpcor() GDALAllRegister(); // reference and secondary images; use band=1 as default // TODO: selecting band + std::cout << "Opening reference image " << param->referenceImageName << "...\n"; GDALImage *referenceImage = new GDALImage(param->referenceImageName, 1, param->mmapSizeInGB); + std::cout << "Opening secondary image " << param->secondaryImageName << "...\n"; GDALImage *secondaryImage = new GDALImage(param->secondaryImageName, 1, param->mmapSizeInGB); cuArrays *offsetImage, *offsetImageRun; @@ -100,9 +102,12 @@ void cuAmpcorController::runAmpcor() << nChunksDown << " x " << nChunksAcross << std::endl; // iterative over chunks down + int message_interval = nChunksDown/10; for(int i = 0; inStreams) { diff --git a/contrib/PyCuAmpcor/src/cudaUtil.h b/contrib/PyCuAmpcor/src/cudaUtil.h index 5a2629c..cc174e0 100644 --- a/contrib/PyCuAmpcor/src/cudaUtil.h +++ b/contrib/PyCuAmpcor/src/cudaUtil.h @@ -81,7 +81,7 @@ inline int gpuDeviceInit(int devID) } checkCudaErrors(cudaSetDevice(devID)); - printf("gpuDeviceInit() Using CUDA Device %d ...\n", devID); + printf("Using CUDA Device %d ...\n", devID); return devID; } @@ -120,4 +120,4 @@ inline void gpuDeviceList() } #endif //__CUDAUTIL_H -//end of file \ No newline at end of file +//end of file