From 804223271862fe393aaed43e8d7dd1189f578530 Mon Sep 17 00:00:00 2001 From: Joseph H Kennedy Date: Mon, 28 Oct 2019 14:51:58 -0800 Subject: [PATCH 01/27] Property is not callable --- components/isceobj/Orbit/Orbit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/isceobj/Orbit/Orbit.py b/components/isceobj/Orbit/Orbit.py index 548ade6..e3cdc1f 100644 --- a/components/isceobj/Orbit/Orbit.py +++ b/components/isceobj/Orbit/Orbit.py @@ -1061,7 +1061,7 @@ class Orbit(Component): ###This wont break the old interface but could cause ###issues at midnight crossing if reference is None: - reference = self.minTime() + reference = self.minTime refEpoch = reference.replace(hour=0, minute=0, second=0, microsecond=0) From 431f857ba7e37f29aca24ffb726d2157f81c8689 Mon Sep 17 00:00:00 2001 From: piyushrpt Date: Mon, 11 Nov 2019 10:42:11 -0800 Subject: [PATCH 02/27] Look for autorift --- contrib/SConscript | 3 +++ 1 file changed, 3 insertions(+) diff --git a/contrib/SConscript b/contrib/SConscript index fb91670..268c0cb 100644 --- a/contrib/SConscript +++ b/contrib/SConscript @@ -78,3 +78,6 @@ SConscript(rfi) SConscript('PyCuAmpcor/SConscript') SConscript('splitSpectrum/SConscript') SConscript('alos2proc/SConscript') + +if os.path.exists('geoAutorift'): + SConscript('geoAutorift/SConscript') From 959a278e477c42fc8dfd138a7b86aaca24ea1615 Mon Sep 17 00:00:00 2001 From: piyushrpt Date: Mon, 11 Nov 2019 10:42:11 -0800 Subject: [PATCH 03/27] Look for autorift and install if found --- contrib/SConscript | 3 +++ 1 file changed, 3 insertions(+) diff --git a/contrib/SConscript b/contrib/SConscript index fb91670..ecc0baf 100644 --- a/contrib/SConscript +++ b/contrib/SConscript @@ -78,3 +78,6 @@ SConscript(rfi) SConscript('PyCuAmpcor/SConscript') SConscript('splitSpectrum/SConscript') SConscript('alos2proc/SConscript') + +if os.path.exists('geo_autoRIFT'): + SConscript('geo_autoRIFT/SConscript') From 5bd6e15a435db60f72c2eb5353ee051f29bce836 Mon Sep 17 00:00:00 2001 From: Lijun Zhu Date: Tue, 19 Nov 2019 16:40:56 -0800 Subject: [PATCH 04/27] setup/setup.py: update modules to work with newest versions of python3 --- setup/setup.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/setup/setup.py b/setup/setup.py index e25db37..3e01849 100755 --- a/setup/setup.py +++ b/setup/setup.py @@ -12,7 +12,7 @@ from __future__ import print_function import sys import os -import urllib2 +import urllib import getopt import re import shutil @@ -57,7 +57,7 @@ def print2log(msg, withtime=True, cmd=False): if withtime: now = datetime.datetime.today() msg = "%s >> %s" % (now.isoformat(), msg) - LOGFILE.write(msg + '\n') + LOGFILE.write((msg + '\n').encode('utf-8')) LOGFILE.flush() os.fsync(LOGFILE) @@ -157,9 +157,9 @@ def downloadfile(url, fname, repeat=1): counter = 0 while counter < repeat: try: - response = urllib2.urlopen(url) + response = urllib.request.urlopen(url) break - except urllib2.URLError, e: + except urllib.request.URLError as e: counter += 1 if hasattr(e, 'reason'): print2log("Failed to reach server. Reason: %s" % e.reason) @@ -851,7 +851,7 @@ class ISCEDeps(object): f = open(self.config, 'rb') lines = f.readlines() for line in lines: - m = re.match("([^#].*?)=([^#]+?)$", line.strip()) + m = re.match("([^#].*?)=([^#]+?)$", line.strip().decode('utf-8')) if m: var = m.group(1).strip() val = m.group(2).strip() @@ -867,7 +867,7 @@ def readSetupConfig(setup_config): f = open(setup_config, 'rb') lines = f.readlines() for line in lines: - m = re.match("([^#].*?)=([^#]+?)$", line.strip()) + m = re.match("([^#].*?)=([^#]+?)$", line.strip().decode('utf-8')) if m: var = m.group(1).strip() val = m.group(2).strip().replace('"', '') @@ -885,7 +885,7 @@ def checkArgs(args): """ try: opts, args = getopt.getopt(args, "h", ["help", "prefix=", "ping=", "config=", "uname=", "download=", "unpack=", "install=", "gcc=", "gpp=", "verbose"]) - except getopt.GetoptError, err: + except getopt.GetoptError as err: print2log("ProgError: %s" % str(err)) usage() sys.exit(2) From d1d951689077675684bb470b47f5e0bac7662188 Mon Sep 17 00:00:00 2001 From: Lijun Zhu Date: Tue, 19 Nov 2019 16:51:06 -0800 Subject: [PATCH 05/27] cuda: remove arch specific flag; add /lib64 to LIBPATH for libcuda.so for most Linux systems --- components/zerodop/GPUampcor/cuda/SConscript | 10 +++++----- components/zerodop/GPUtopozero/cuda/compilation | 2 +- scons_tools/cuda.py | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/components/zerodop/GPUampcor/cuda/SConscript b/components/zerodop/GPUampcor/cuda/SConscript index 7d76208..c0ec346 100644 --- a/components/zerodop/GPUampcor/cuda/SConscript +++ b/components/zerodop/GPUampcor/cuda/SConscript @@ -2,19 +2,19 @@ #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Copyright 2010 California Institute of Technology. ALL RIGHTS RESERVED. -# +# # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at -# +# # http://www.apache.org/licenses/LICENSE-2.0 -# +# # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -# +# # United States Government Sponsorship acknowledged. This software is subject to # U.S. export control laws and regulations and has been classified as 'EAR99 NLR' # (No [Export] License Required except when exporting to an embargoed country, @@ -49,7 +49,7 @@ if envGPUampcor['GPU_ACC_ENABLED']: build_base += "-ccbin " + envGPUampcor['NVCC_CCBIN'] + " " else: print('Assuming default system compiler for nvcc.') - build_base += "-arch=sm_35 -shared -Xcompiler -fPIC -O3 " + build_base += "-shared -Xcompiler -fPIC -O3 " build_cmd = build_base + "-dc -m64 -o $TARGET -c $SOURCE" built_path = os.path.join(build, 'gpu-ampcor.o') linked_path = os.path.join(build, 'gpu-ampcor-linked.o') diff --git a/components/zerodop/GPUtopozero/cuda/compilation b/components/zerodop/GPUtopozero/cuda/compilation index ff97c2a..246d366 100644 --- a/components/zerodop/GPUtopozero/cuda/compilation +++ b/components/zerodop/GPUtopozero/cuda/compilation @@ -1,2 +1,2 @@ -nvcc -arch=sm_35 -Xcompiler -fPIC -o gpu-topo.o -c Topo.cu +nvcc -Xcompiler -fPIC -o gpu-topo.o -c Topo.cu cp -f gpu-topo.o .. diff --git a/scons_tools/cuda.py b/scons_tools/cuda.py index fe2f6d0..9dce35b 100644 --- a/scons_tools/cuda.py +++ b/scons_tools/cuda.py @@ -52,7 +52,7 @@ def generate(env): # default flags for the NVCC compiler env['STATICNVCCFLAGS'] = '' env['SHAREDNVCCFLAGS'] = '' - env['ENABLESHAREDNVCCFLAG'] = '-arch=sm_35 -shared -Xcompiler -fPIC' + env['ENABLESHAREDNVCCFLAG'] = '-shared -Xcompiler -fPIC' # default NVCC commands env['STATICNVCCCMD'] = '$NVCC -o $TARGET -c $NVCCFLAGS $STATICNVCCFLAGS $SOURCES' @@ -153,7 +153,7 @@ def generate(env): #env.Append(LIBPATH=[cudaSDKPath + '/lib', cudaSDKPath + '/common/lib' + cudaSDKSubLibDir, cudaToolkitPath + '/lib']) env.Append(CUDACPPPATH=[cudaToolkitPath + '/include']) - env.Append(CUDALIBPATH=[cudaToolkitPath + '/lib', cudaToolkitPath + '/lib64']) + env.Append(CUDALIBPATH=[cudaToolkitPath + '/lib', cudaToolkitPath + '/lib64', '/lib64']) env.Append(CUDALIBS=['cudart']) def exists(env): From df4db88ffae9ac6f2019297e01bcca12ea07f05b Mon Sep 17 00:00:00 2001 From: Lijun Zhu Date: Tue, 19 Nov 2019 16:59:49 -0800 Subject: [PATCH 06/27] PyCuAmpcor: updated to the most recent version with gdal input --- contrib/PyCuAmpcor/SConscript | 4 +- contrib/PyCuAmpcor/examples/GeoTiffSample.py | 63 +++ ...dGrossOffsetSample.py => glacierSample.py} | 21 +- .../examples/varyGrossOffsetSample.py | 26 +- contrib/PyCuAmpcor/src/GDALImage.cu | 154 +++++++ contrib/PyCuAmpcor/src/GDALImage.h | 79 ++++ contrib/PyCuAmpcor/src/Makefile | 25 +- contrib/PyCuAmpcor/src/PyCuAmpcor.pyx | 190 ++++---- contrib/PyCuAmpcor/src/SConscript | 2 +- contrib/PyCuAmpcor/src/cuAmpcorChunk.cu | 408 ++++++++++++------ contrib/PyCuAmpcor/src/cuAmpcorChunk.h | 81 ++-- contrib/PyCuAmpcor/src/cuAmpcorController.cu | 131 +++--- contrib/PyCuAmpcor/src/cuAmpcorParameter.cu | 164 +++---- contrib/PyCuAmpcor/src/cuAmpcorParameter.h | 110 ++--- contrib/PyCuAmpcor/src/cuAmpcorUtil.h | 39 +- contrib/PyCuAmpcor/src/cuArrays.cu | 13 + contrib/PyCuAmpcor/src/cuArraysCopy.cu | 316 +++++++++----- contrib/PyCuAmpcor/src/cuCorrNormalization.cu | 8 +- contrib/PyCuAmpcor/src/cuEstimateStats.cu | 81 +++- contrib/PyCuAmpcor/src/setup.py | 13 +- 20 files changed, 1333 insertions(+), 595 deletions(-) create mode 100644 contrib/PyCuAmpcor/examples/GeoTiffSample.py rename contrib/PyCuAmpcor/examples/{fixedGrossOffsetSample.py => glacierSample.py} (78%) create mode 100644 contrib/PyCuAmpcor/src/GDALImage.cu create mode 100644 contrib/PyCuAmpcor/src/GDALImage.h diff --git a/contrib/PyCuAmpcor/SConscript b/contrib/PyCuAmpcor/SConscript index 5ef82ea..227a34f 100644 --- a/contrib/PyCuAmpcor/SConscript +++ b/contrib/PyCuAmpcor/SConscript @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 import os @@ -28,7 +28,7 @@ if envPyCuAmpcor['GPU_ACC_ENABLED']: if not os.path.exists(initFile): with open(initFile, 'w') as fout: - fout.write("#!/usr/bin/env python") + fout.write("#!/usr/bin/env python3") listFiles = [initFile] envPyCuAmpcor.Install(install, listFiles) diff --git a/contrib/PyCuAmpcor/examples/GeoTiffSample.py b/contrib/PyCuAmpcor/examples/GeoTiffSample.py new file mode 100644 index 0000000..a700524 --- /dev/null +++ b/contrib/PyCuAmpcor/examples/GeoTiffSample.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +# +# Test program to run ampcor with GPU +# For two GeoTiff images +# + +import argparse +import numpy as np +from PyCuAmpcor import PyCuAmpcor + + +def main(): + ''' + main program + ''' + + objOffset = PyCuAmpcor() # create the processor + + objOffset.algorithm = 0 # cross-correlation method 0=freq 1=time + objOffset.deviceID = 0 # GPU device id to be used + objOffset.nStreams = 2 # cudaStreams; multiple streams to overlap data transfer with gpu calculations + objOffset.masterImageName = "master.tif" + objOffset.masterImageHeight = 16480 # RasterYSize + objOffset.masterImageWidth = 17000 # RasterXSize + objOffset.slaveImageName = "slave.tif" + objOffset.slaveImageHeight = 16480 + objOffset.slaveImageWidth = 17000 + objOffset.windowSizeWidth = 64 # template window size + objOffset.windowSizeHeight = 64 + objOffset.halfSearchRangeDown = 20 # search range + objOffset.halfSearchRangeAcross = 20 + objOffset.derampMethod = 1 # deramping for complex signal, set to 1 for real images + + objOffset.skipSampleDown = 128 # strides between windows + objOffset.skipSampleAcross = 64 + # gpu processes several windows in one batch/Chunk + # total windows in Chunk = numberWindowDownInChunk*numberWindowAcrossInChunk + # the max number of windows depending on gpu memory and type + objOffset.numberWindowDownInChunk = 1 + objOffset.numberWindowAcrossInChunk = 10 + objOffset.corrSurfaceOverSamplingFactor = 8 # oversampling factor for correlation surface + objOffset.corrSurfaceZoomInWindow = 16 # area in correlation surface to be oversampled + objOffset.corrSufaceOverSamplingMethod = 1 # fft or sinc oversampler + objOffset.useMmap = 1 # default using memory map as buffer, if having troubles, set to 0 + objOffset.mmapSize = 1 # mmap or buffer size used for transferring data from file to gpu, in GB + + objOffset.numberWindowDown = 40 # number of windows to be processed + objOffset.numberWindowAcross = 100 + # if to process the whole image; some math needs to be done + # margin = 0 # margins to be neglected + #objOffset.numberWindowDown = (objOffset.slaveImageHeight - 2*margin - 2*objOffset.halfSearchRangeDown - objOffset.windowSizeHeight) // objOffset.skipSampleDown + #objOffset.numberWindowAcross = (objOffset.slaveImageWidth - 2*margin - 2*objOffset.halfSearchRangeAcross - objOffset.windowSizeWidth) // objOffset.skipSampleAcross + + objOffset.setupParams() + objOffset.masterStartPixelDownStatic = objOffset.halfSearchRangeDown # starting pixel offset + objOffset.masterStartPixelAcrossStatic = objOffset.halfSearchRangeDown + objOffset.setConstantGrossOffset(0, 0) # gross offset between master and slave images + objOffset.checkPixelInImageRange() # check whether there is something wrong with + objOffset.runAmpcor() + + +if __name__ == '__main__': + diff --git a/contrib/PyCuAmpcor/examples/fixedGrossOffsetSample.py b/contrib/PyCuAmpcor/examples/glacierSample.py similarity index 78% rename from contrib/PyCuAmpcor/examples/fixedGrossOffsetSample.py rename to contrib/PyCuAmpcor/examples/glacierSample.py index a81c897..8946512 100644 --- a/contrib/PyCuAmpcor/examples/fixedGrossOffsetSample.py +++ b/contrib/PyCuAmpcor/examples/glacierSample.py @@ -1,14 +1,14 @@ #!/usr/bin/env python3 -# +# # test_cuAmpcor.py # Test program to run ampcor with GPU # -# +# import argparse import numpy as np -#from PyCuAmpcor import PyCuAmpcor -from isce.components.contrib.PyCuAmpcor import PyCuAmpcor +from PyCuAmpcor import PyCuAmpcor + def main(): ''' @@ -20,10 +20,10 @@ def main(): objOffset.algorithm = 0 objOffset.deviceID = 0 # -1:let system find the best GPU objOffset.nStreams = 2 #cudaStreams - objOffset.masterImageName = "master.slc" + objOffset.masterImageName = "20131213.slc.vrt" objOffset.masterImageHeight = 43008 objOffset.masterImageWidth = 24320 - objOffset.slaveImageName = "slave.slc" + objOffset.slaveImageName = "20131221.slc.vrt" objOffset.slaveImageHeight = 43008 objOffset.slaveImageWidth = 24320 objOffset.windowSizeWidth = 64 @@ -38,8 +38,9 @@ def main(): objOffset.numberWindowDownInChunk = 10 objOffset.numberWindowAcrossInChunk = 10 objOffset.corrSurfaceOverSamplingFactor = 8 - objOffset.corrSurfaceZoomInWindow = 16 - objOffset.corrSufaceOverSamplingMethod = 1 + objOffset.corrSurfaceZoomInWindow = 16 + objOffset.corrSufaceOverSamplingMethod = 1 + objOffset.useMmap = 1 objOffset.mmapSize = 8 objOffset.setupParams() @@ -48,8 +49,8 @@ def main(): objOffset.setConstantGrossOffset(642, -30) objOffset.checkPixelInImageRange() objOffset.runAmpcor() - + if __name__ == '__main__': - + main() diff --git a/contrib/PyCuAmpcor/examples/varyGrossOffsetSample.py b/contrib/PyCuAmpcor/examples/varyGrossOffsetSample.py index e4a56f2..6feeb89 100644 --- a/contrib/PyCuAmpcor/examples/varyGrossOffsetSample.py +++ b/contrib/PyCuAmpcor/examples/varyGrossOffsetSample.py @@ -1,27 +1,27 @@ #!/usr/bin/env python3 -# +# from PyCuAmpcor import PyCuAmpcor import numpy as np -def main(): +def main(): ''' Set parameters manually and run ampcor ''' objOffset = PyCuAmpcor() - + #step 1 set constant parameters - objOffset.masterImageName = "master.slc" + objOffset.masterImageName = "master.slc.vrt" objOffset.masterImageHeight = 128 objOffset.masterImageWidth = 128 - objOffset.slaveImageName = "slave.slc" + objOffset.slaveImageName = "slave.slc.vrt" objOffset.masterImageHeight = 128 - objOffset.masterImageWidth = 128 + objOffset.masterImageWidth = 128 objOffset.skipSampleDown = 2 objOffset.skipSampleAcross = 2 objOffset.windowSizeHeight = 16 objOffset.windowSizeWidth = 16 - objOffset.halfSearchRangeDown = 20 + objOffset.halfSearchRangeDown = 20 objOffset.halfSearchRangeAcross = 20 objOffset.numberWindowDown = 2 objOffset.numberWindowAcross = 2 @@ -29,19 +29,19 @@ def main(): objOffset.numberWindowAcrossInChunk = 2 # 2 set other dependent parameters and allocate aray parameters objOffset.setupParams() - + #3 set gross offsets: constant or varying - objOffset.masterStartPixelDownStatic = objOffset.halfSearchRangeDown + objOffset.masterStartPixelDownStatic = objOffset.halfSearchRangeDown objOffset.masterStartPixelAcrossStatic = objOffset.halfSearchRangeAcross - vD = np.random.randint(0, 10, size =objOffset.numberWindows, dtype=np.int32) + vD = np.random.randint(0, 10, size =objOffset.numberWindows, dtype=np.int32) vA = np.random.randint(0, 1, size = objOffset.numberWindows, dtype=np.int32) objOffset.setVaryingGrossOffset(vD, vA) - + objOffset.checkPixelInImageRange() #4 run ampcor objOffset.runAmpcor() - - + + if __name__ == '__main__': main() diff --git a/contrib/PyCuAmpcor/src/GDALImage.cu b/contrib/PyCuAmpcor/src/GDALImage.cu new file mode 100644 index 0000000..2c42cbb --- /dev/null +++ b/contrib/PyCuAmpcor/src/GDALImage.cu @@ -0,0 +1,154 @@ +#include "GDALImage.h" +#include + +#include +#include +#include +#include +#include +#include +#include "cudaError.h" +#include +#include + + +/** + * \brief Constructor + * + * @param filename a std::string with the raster image file name + */ + +GDALImage::GDALImage(std::string filename, int band, int cacheSizeInGB, int useMmap) + : _useMmap(useMmap) +{ + // open the file as dataset + _poDataset = (GDALDataset *) GDALOpen(filename.c_str(), GA_ReadOnly ); + // if something is wrong, throw an exception + // GDAL reports the error message + if(!_poDataset) + throw; + + // check the band info + int count = _poDataset->GetRasterCount(); + if(band > count) + { + std::cout << "The desired band " << band << " is greated than " << count << " bands available"; + throw; + } + + // get the desired band + _poBand = _poDataset->GetRasterBand(band); + if(!_poBand) + throw; + + // get the width(x), and height(y) + _width = _poBand->GetXSize(); + _height = _poBand->GetYSize(); + + _dataType = _poBand->GetRasterDataType(); + // determine the image type + _isComplex = GDALDataTypeIsComplex(_dataType); + // determine the pixel size in bytes + _pixelSize = GDALGetDataTypeSize(_dataType); + + _bufferSize = 1024*1024*cacheSizeInGB; + + // checking whether using memory map + if(_useMmap) { + + char **papszOptions = NULL; + // if cacheSizeInGB = 0, use default + // else set the option + if(cacheSizeInGB > 0) + papszOptions = CSLSetNameValue( papszOptions, + "CACHE_SIZE", + std::to_string(_bufferSize).c_str()); + + // space between two lines + GIntBig pnLineSpace; + // set up the virtual mem buffer + _poBandVirtualMem = GDALGetVirtualMemAuto( + static_cast(_poBand), + GF_Read, + &_pixelSize, + &pnLineSpace, + papszOptions); + + // check it + if(!_poBandVirtualMem) + throw; + + // get the starting pointer + _memPtr = CPLVirtualMemGetAddr(_poBandVirtualMem); + } + else { // use a buffer + checkCudaErrors(cudaMallocHost((void **)&_memPtr, _bufferSize)); + } + + // make sure memPtr is not Null + if (!_memPtr) + throw; + + // all done +} + + +/// load a tile of data h_tile x w_tile from CPU (mmap) to GPU +/// @param dArray pointer for array in device memory +/// @param h_offset Down/Height offset +/// @param w_offset Across/Width offset +/// @param h_tile Down/Height tile size +/// @param w_tile Across/Width tile size +/// @param stream CUDA stream for copying +void GDALImage::loadToDevice(void *dArray, size_t h_offset, size_t w_offset, size_t h_tile, size_t w_tile, cudaStream_t stream) +{ + size_t tileStartOffset = (h_offset*_width + w_offset)*_pixelSize; + + char * startPtr = (char *)_memPtr ; + startPtr += tileStartOffset; + + // @note + // We assume down/across directions as rows/cols. Therefore, SLC mmap and device array are both row major. + // cuBlas assumes both source and target arrays are column major. + // To use cublasSetMatrix, we need to switch w_tile/h_tile for rows/cols + // checkCudaErrors(cublasSetMatrixAsync(w_tile, h_tile, sizeof(float2), startPtr, width, dArray, w_tile, stream)); + if (_useMmap) + checkCudaErrors(cudaMemcpy2DAsync(dArray, w_tile*_pixelSize, startPtr, _width*_pixelSize, + w_tile*_pixelSize, h_tile, cudaMemcpyHostToDevice,stream)); + else { + // get the total tile size in bytes + size_t tileSize = h_tile*w_tile*_pixelSize; + // if the size is bigger than existing buffer, reallocate + if (tileSize > _bufferSize) { + // maybe we need to make it to fit the pagesize + _bufferSize = tileSize; + checkCudaErrors(cudaFree(_memPtr)); + checkCudaErrors(cudaMallocHost((void **)&_memPtr, _bufferSize)); + } + // copy from file to buffer + CPLErr err = _poBand->RasterIO(GF_Read, //eRWFlag + w_offset, h_offset, //nXOff, nYOff + w_tile, h_tile, // nXSize, nYSize + _memPtr, // pData + w_tile*h_tile, 1, // nBufXSize, nBufYSize + _dataType, //eBufType + 0, 0, //nPixelSpace, nLineSpace in pData + NULL //psExtraArg extra resampling callback + ); + + if(err != CE_None) + throw; + // copy from buffer to gpu + checkCudaErrors(cudaMemcpyAsync(dArray, _memPtr, tileSize, cudaMemcpyHostToDevice, stream)); + } +} + +GDALImage::~GDALImage() +{ + // free the virtual memory + CPLVirtualMemFree(_poBandVirtualMem), + // free the GDAL Dataset, close the file + delete _poDataset; +} + +// end of file diff --git a/contrib/PyCuAmpcor/src/GDALImage.h b/contrib/PyCuAmpcor/src/GDALImage.h new file mode 100644 index 0000000..3322a2a --- /dev/null +++ b/contrib/PyCuAmpcor/src/GDALImage.h @@ -0,0 +1,79 @@ +// -*- c++ -*- +/** + * \brief Class for an image described GDAL vrt + * + * only complex (pixelOffset=8) or real(pixelOffset=4) images are supported, such as SLC and single-precision TIFF + */ + +#ifndef __GDALIMAGE_H +#define __GDALIMAGE_H + +#include +#include +#include +#include + +class GDALImage{ + +public: + using size_t = std::size_t; + +private: + size_t _fileSize; + int _height; + int _width; + + // buffer pointer + void * _memPtr = NULL; + + int _pixelSize; //in bytes + + int _isComplex; + + size_t _bufferSize; + int _useMmap; + + GDALDataType _dataType; + CPLVirtualMem * _poBandVirtualMem = NULL; + GDALDataset * _poDataset = NULL; + GDALRasterBand * _poBand = NULL; + +public: + GDALImage() = delete; + GDALImage(std::string fn, int band=1, int cacheSizeInGB=0, int useMmap=1); + + void * getmemPtr() + { + return(_memPtr); + } + + size_t getFileSize() + { + return (_fileSize); + } + + size_t getHeight() { + return (_height); + } + + size_t getWidth() + { + return (_width); + } + + int getPixelSize() + { + return _pixelSize; + } + + bool isComplex() + { + return _isComplex; + } + + void loadToDevice(void *dArray, size_t h_offset, size_t w_offset, size_t h_tile, size_t w_tile, cudaStream_t stream); + ~GDALImage(); + +}; + +#endif //__GDALIMAGE_H diff --git a/contrib/PyCuAmpcor/src/Makefile b/contrib/PyCuAmpcor/src/Makefile index 2138093..b789a59 100644 --- a/contrib/PyCuAmpcor/src/Makefile +++ b/contrib/PyCuAmpcor/src/Makefile @@ -3,23 +3,24 @@ PROJECT = CUAMPCOR LDFLAGS = -lcuda -lcudart -lcufft -lcublas CXXFLAGS = -std=c++11 -fpermissive -fPIC -shared NVCCFLAGS = -ccbin g++ -m64 \ - -gencode arch=compute_35,code=sm_35 \ + -gencode arch=compute_35,code=sm_35 \ + -gencode arch=compute_60,code=sm_60 \ -Xcompiler -fPIC -shared -Wno-deprecated-gpu-targets \ -ftz=false -prec-div=true -prec-sqrt=true CXX=g++ NVCC=nvcc -DEPS = cudaUtil.h cudaError.h cuArrays.h SlcImage.h cuAmpcorParameter.h -OBJS = SlcImage.o cuArrays.o cuArraysCopy.o cuArraysPadding.o cuOverSampler.o \ +DEPS = cudaUtil.h cudaError.h cuArrays.h GDALImage.h cuAmpcorParameter.h +OBJS = GDALImage.o cuArrays.o cuArraysCopy.o cuArraysPadding.o cuOverSampler.o \ cuSincOverSampler.o cuDeramp.o cuOffset.o \ cuCorrNormalization.o cuAmpcorParameter.o cuCorrTimeDomain.o cuCorrFrequency.o \ cuAmpcorChunk.o cuAmpcorController.o cuEstimateStats.o -all: cuampcor +all: pyampcor -SlcImage.o: SlcImage.cu $(DEPS) - $(NVCC) $(NVCCFLAGS) -c -o $@ SlcImage.cu +GDALImage.o: GDALImage.cu $(DEPS) + $(NVCC) $(NVCCFLAGS) -c -o $@ GDALImage.cu cuArrays.o: cuArrays.cu $(DEPS) $(NVCC) $(NVCCFLAGS) -c -o $@ cuArrays.cu @@ -45,8 +46,8 @@ cuOffset.o: cuOffset.cu $(DEPS) cuCorrNormalization.o: cuCorrNormalization.cu $(DEPS) $(NVCC) $(NVCCFLAGS) -c -o $@ cuCorrNormalization.cu -cuAmpcorParameter.o: cuAmpcorParameter.cu - $(NVCC) $(NVCCFLAGS) -c -o $@ cuAmpcorParameter.cu +cuAmpcorParameter.o: cuAmpcorParameter.cu + $(NVCC) $(NVCCFLAGS) -c -o $@ cuAmpcorParameter.cu cuCorrTimeDomain.o: cuCorrTimeDomain.cu $(DEPS) $(NVCC) $(NVCCFLAGS) -c -o $@ cuCorrTimeDomain.cu @@ -54,8 +55,8 @@ cuCorrTimeDomain.o: cuCorrTimeDomain.cu $(DEPS) cuCorrFrequency.o: cuCorrFrequency.cu $(DEPS) cuCorrFrequency.h $(NVCC) $(NVCCFLAGS) -c -o $@ cuCorrFrequency.cu -cuAmpcorChunk.o: cuAmpcorChunk.cu cuAmpcorUtil.h $(DEPS) - $(NVCC) $(NVCCFLAGS) -c -o $@ cuAmpcorChunk.cu +cuAmpcorChunk.o: cuAmpcorChunk.cu cuAmpcorUtil.h $(DEPS) + $(NVCC) $(NVCCFLAGS) -c -o $@ cuAmpcorChunk.cu cuAmpcorController.o: cuAmpcorController.cu $(NVCC) $(NVCCFLAGS) -c -o $@ cuAmpcorController.cu @@ -64,8 +65,8 @@ cuEstimateStats.o: cuEstimateStats.cu $(NVCC) $(NVCCFLAGS) -c -o $@ cuEstimateStats.cu -cuampcor: $(OBJS) +pyampcor: $(OBJS) rm -f PyCuAmpcor.cpp && python3 setup.py build_ext --inplace clean: - rm -rf *.o *so build *~ PyCuAmpcor.cpp ctest *.dat + rm -rf *.o *so build *~ PyCuAmpcor.cpp ctest *.dat diff --git a/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx b/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx index d5f6c8d..857966d 100644 --- a/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx +++ b/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx @@ -1,6 +1,6 @@ -# +# # PYX file to control Python module interface to underlying CUDA-Ampcor code -# +# from libcpp.string cimport string import numpy as np cimport numpy as np @@ -9,13 +9,13 @@ cimport numpy as np cdef extern from "cudaUtil.h": int gpuDeviceInit(int) void gpuDeviceList() - int gpuGetMaxGflopsDeviceId() + int gpuGetMaxGflopsDeviceId() def listGPU(): gpuDeviceList() def findGPU(): - return gpuGetMaxGflopsDeviceId() + return gpuGetMaxGflopsDeviceId() def setGPU(int id): return gpuDeviceInit(id) @@ -24,90 +24,92 @@ def setGPU(int id): cdef extern from "cuAmpcorParameter.h": cdef cppclass cuAmpcorParameter: cuAmpcorParameter() except + - int algorithm ## Cross-correlation algorithm: 0=freq domain 1=time domain - int deviceID ## Targeted GPU device ID: use -1 to auto select - int nStreams ## Number of streams to asynchonize data transfers and compute kernels + int algorithm ## Cross-correlation algorithm: 0=freq domain 1=time domain + int deviceID ## Targeted GPU device ID: use -1 to auto select + int nStreams ## Number of streams to asynchonize data transfers and compute kernels int derampMethod ## Method for deramping 0=None, 1=average, 2=phase gradient - + ## chip or window size for raw data int windowSizeHeightRaw ## Template window height (original size) - int windowSizeWidthRaw ## Template window width (original size) - int searchWindowSizeHeightRaw ## Search window height (original size) + int windowSizeWidthRaw ## Template window width (original size) + int searchWindowSizeHeightRaw ## Search window height (original size) int searchWindowSizeWidthRaw ## Search window width (orignal size) - int halfSearchRangeDownRaw ##(searchWindowSizeHeightRaw-windowSizeHeightRaw)/2 + int halfSearchRangeDownRaw ##(searchWindowSizeHeightRaw-windowSizeHeightRaw)/2 int halfSearchRangeAcrossRaw ##(searchWindowSizeWidthRaw-windowSizeWidthRaw)/2 ## chip or window size after oversampling int rawDataOversamplingFactor ## Raw data overampling factor (from original size to oversampled size) - - ## strides between chips/windows + + ## strides between chips/windows int skipSampleDownRaw ## Skip size between neighboring windows in Down direction (original size) int skipSampleAcrossRaw ## Skip size between neighboring windows in across direction (original size) - + ## Zoom in region near location of max correlation - int zoomWindowSize ## Zoom-in window size in correlation surface (same for down and across directions) + int zoomWindowSize ## Zoom-in window size in correlation surface (same for down and across directions) int oversamplingFactor ## Oversampling factor for interpolating correlation surface - int oversamplingMethod - - float thresholdSNR ## Threshold of Signal noise ratio to remove noisy data - + int oversamplingMethod + + float thresholdSNR ## Threshold of Signal noise ratio to remove noisy data + ##master image string masterImageName ## master SLC image name int imageDataType1 ## master image data type, 2=cfloat=complex=float2 1=float - int masterImageHeight ## master image height + int masterImageHeight ## master image height int masterImageWidth ## master image width - + ##slave image string slaveImageName ## slave SLC image name int imageDataType2 ## slave image data type, 2=cfloat=complex=float2 1=float - int slaveImageHeight ## slave image height + int slaveImageHeight ## slave image height int slaveImageWidth ## slave image width - - int mmapSizeInGB ## mmap buffer size in unit of Gigabytes - + + int useMmap ## whether to use mmap + int mmapSizeInGB ## mmap buffer size in unit of Gigabytes (if not mmmap, the buffer size) + ## total number of chips/windows int numberWindowDown ## number of total windows (down) int numberWindowAcross ## number of total windows (across) int numberWindows ## numberWindowDown*numberWindowAcross - + ## number of chips/windows in a batch/chunk int numberWindowDownInChunk ## number of windows processed in a chunk (down) int numberWindowAcrossInChunk ## number of windows processed in a chunk (across) int numberWindowsInChunk ## numberWindowDownInChunk*numberWindowAcrossInChunk int numberChunkDown ## number of chunks (down) int numberChunkAcross ## number of chunks (across) - int numberChunks + int numberChunks - int *masterStartPixelDown ## master starting pixels for each window (down) + int *masterStartPixelDown ## master starting pixels for each window (down) int *masterStartPixelAcross ## master starting pixels for each window (across) - int *slaveStartPixelDown ## slave starting pixels for each window (down) - int *slaveStartPixelAcross ## slave starting pixels for each window (across) + int *slaveStartPixelDown ## slave starting pixels for each window (down) + int *slaveStartPixelAcross ## slave starting pixels for each window (across) int *grossOffsetDown ## Gross offsets between master and slave windows (down) : slaveStartPixel - masterStartPixel - int *grossOffsetAcross ## Gross offsets between master and slave windows (across) + int *grossOffsetAcross ## Gross offsets between master and slave windows (across) int grossOffsetDown0 ## constant gross offset (down) int grossOffsetAcross0 ## constant gross offset (across) - int masterStartPixelDown0 ## the first pixel of master image (down), be adjusted with margins and gross offset + int masterStartPixelDown0 ## the first pixel of master image (down), be adjusted with margins and gross offset int masterStartPixelAcross0 ## the first pixel of master image (across) int *masterChunkStartPixelDown ## array of starting pixels for all master chunks (down) int *masterChunkStartPixelAcross ## array of starting pixels for all master chunks (across) int *slaveChunkStartPixelDown ## array of starting pixels for all slave chunks (down) int *slaveChunkStartPixelAcross ## array of starting pixels for all slave chunks (across) - int *masterChunkHeight ## array of heights of all master chunks, required when loading chunk to GPU + int *masterChunkHeight ## array of heights of all master chunks, required when loading chunk to GPU int *masterChunkWidth ## array of width of all master chunks int *slaveChunkHeight ## array of width of all master chunks int *slaveChunkWidth ## array of width of all slave chunks - int maxMasterChunkHeight ## max height for all master/slave chunks, determine the size of reading cache in GPU - int maxMasterChunkWidth ## max width for all master chunks, determine the size of reading cache in GPU + int maxMasterChunkHeight ## max height for all master/slave chunks, determine the size of reading cache in GPU + int maxMasterChunkWidth ## max width for all master chunks, determine the size of reading cache in GPU int maxSlaveChunkHeight int maxSlaveChunkWidth - - string grossOffsetImageName - string offsetImageName ## Output Offset fields filename + + string grossOffsetImageName + string offsetImageName ## Output Offset fields filename string snrImageName ## Output SNR filename - void setStartPixels(int*, int*, int*, int*) - void setStartPixels(int, int, int*, int*) - void setStartPixels(int, int, int, int) - void checkPixelInImageRange() ## check whether - + string covImageName ## Output COV filename + void setStartPixels(int*, int*, int*, int*) + void setStartPixels(int, int, int*, int*) + void setStartPixels(int, int, int, int) + void checkPixelInImageRange() ## check whether + void setupParameters() ## Process other parameters after Python Inpu cdef extern from "cuAmpcorController.h": @@ -115,34 +117,40 @@ cdef extern from "cuAmpcorController.h": cuAmpcorController() except + cuAmpcorParameter *param void runAmpcor() - + cdef class PyCuAmpcor(object): ''' - Python interface for cuda Ampcor + Python interface for cuda Ampcor ''' cdef cuAmpcorController c_cuAmpcor def __cinit__(self): - return - + return + @property def algorithm(self): - return self.c_cuAmpcor.param.algorithm + return self.c_cuAmpcor.param.algorithm @algorithm.setter def algorithm(self, int a): self.c_cuAmpcor.param.algorithm = a @property def deviceID(self): - return self.c_cuAmpcor.param.deviceID + return self.c_cuAmpcor.param.deviceID @deviceID.setter def deviceID(self, int a): self.c_cuAmpcor.param.deviceID = a @property def nStreams(self): - return self.c_cuAmpcor.param.nStreams + return self.c_cuAmpcor.param.nStreams @nStreams.setter def nStreams(self, int a): self.c_cuAmpcor.param.nStreams = a - @property + @property + def useMmap(self): + return self.c_cuAmpcor.param.useMmap + @useMmap.setter + def useMmap(self, int a): + self.c_cuAmpcor.param.useMmap = a + @property def mmapSize(self): return self.c_cuAmpcor.param.mmapSizeInGB @mmapSize.setter @@ -150,19 +158,19 @@ cdef class PyCuAmpcor(object): self.c_cuAmpcor.param.mmapSizeInGB = a @property def derampMethod(self): - return self.c_cuAmpcor.param.derampMethod + return self.c_cuAmpcor.param.derampMethod @derampMethod.setter def derampMethod(self, int a): self.c_cuAmpcor.param.derampMethod = a @property def windowSizeHeight(self): - return self.c_cuAmpcor.param.windowSizeHeightRaw + return self.c_cuAmpcor.param.windowSizeHeightRaw @windowSizeHeight.setter def windowSizeHeight(self, int a): self.c_cuAmpcor.param.windowSizeHeightRaw = a @property def windowSizeWidth(self): - return self.c_cuAmpcor.param.windowSizeWidthRaw + return self.c_cuAmpcor.param.windowSizeWidthRaw @windowSizeWidth.setter def windowSizeWidth(self, int a): self.c_cuAmpcor.param.windowSizeWidthRaw = a @@ -200,7 +208,7 @@ cdef class PyCuAmpcor(object): @skipSampleAcross.setter def skipSampleAcross(self, int a): self.c_cuAmpcor.param.skipSampleAcrossRaw = a - + @property def rawDataOversamplingFactor(self): """anti-aliasing oversampling factor""" @@ -229,7 +237,7 @@ cdef class PyCuAmpcor(object): @corrSufaceOverSamplingMethod.setter def corrSufaceOverSamplingMethod(self, int a): self.c_cuAmpcor.param.oversamplingMethod = a - @property + @property def masterImageName(self): return self.c_cuAmpcor.param.masterImageName @masterImageName.setter @@ -241,12 +249,12 @@ cdef class PyCuAmpcor(object): @slaveImageName.setter def slaveImageName(self, str a): self.c_cuAmpcor.param.slaveImageName = a.encode() - @property + @property def masterImageName(self): return self.c_cuAmpcor.param.masterImageName @masterImageName.setter def masterImageName(self, str a): - self.c_cuAmpcor.param.masterImageName = a.encode() + self.c_cuAmpcor.param.masterImageName = a.encode() @property def masterImageHeight(self): return self.c_cuAmpcor.param.masterImageHeight @@ -258,7 +266,7 @@ cdef class PyCuAmpcor(object): return self.c_cuAmpcor.param.masterImageWidth @masterImageWidth.setter def masterImageWidth(self, int a): - self.c_cuAmpcor.param.masterImageWidth=a + self.c_cuAmpcor.param.masterImageWidth=a @property def slaveImageHeight(self): return self.c_cuAmpcor.param.slaveImageHeight @@ -270,8 +278,8 @@ cdef class PyCuAmpcor(object): return self.c_cuAmpcor.param.slaveImageWidth @slaveImageWidth.setter def slaveImageWidth(self, int a): - self.c_cuAmpcor.param.slaveImageWidth=a - + self.c_cuAmpcor.param.slaveImageWidth=a + @property def numberWindowDown(self): return self.c_cuAmpcor.param.numberWindowDown @@ -283,11 +291,11 @@ cdef class PyCuAmpcor(object): return self.c_cuAmpcor.param.numberWindowAcross @numberWindowAcross.setter def numberWindowAcross(self, int a): - self.c_cuAmpcor.param.numberWindowAcross = a + self.c_cuAmpcor.param.numberWindowAcross = a @property def numberWindows(self): return self.c_cuAmpcor.param.numberWindows - + @property def numberWindowDownInChunk(self): return self.c_cuAmpcor.param.numberWindowDownInChunk @@ -299,7 +307,7 @@ cdef class PyCuAmpcor(object): return self.c_cuAmpcor.param.numberWindowAcrossInChunk @numberWindowAcrossInChunk.setter def numberWindowAcrossInChunk(self, int a): - self.c_cuAmpcor.param.numberWindowAcrossInChunk = a + self.c_cuAmpcor.param.numberWindowAcrossInChunk = a @property def numberChunkDown(self): return self.c_cuAmpcor.param.numberChunkDown @@ -309,9 +317,9 @@ cdef class PyCuAmpcor(object): @property def numberChunks(self): return self.c_cuAmpcor.param.numberChunks - - - ## gross offets + + + ## gross offets @property def grossOffsetImageName(self): return self.c_cuAmpcor.param.grossOffsetImageName @@ -324,13 +332,21 @@ cdef class PyCuAmpcor(object): @offsetImageName.setter def offsetImageName(self, str a): self.c_cuAmpcor.param.offsetImageName = a.encode() + @property def snrImageName(self): return self.c_cuAmpcor.param.snrImageName @snrImageName.setter def snrImageName(self, str a): self.c_cuAmpcor.param.snrImageName = a.encode() - + + @property + def covImageName(self): + return self.c_cuAmpcor.param.covImageName + @covImageName.setter + def covImageName(self, str a): + self.c_cuAmpcor.param.covImageName = a.encode() + @property def masterStartPixelDownStatic(self): return self.c_cuAmpcor.param.masterStartPixelDown0 @@ -342,20 +358,20 @@ cdef class PyCuAmpcor(object): return self.c_cuAmpcor.param.masterStartPixelAcross0 @masterStartPixelAcrossStatic.setter def masterStartPixelAcrossStatic(self, int a): - self.c_cuAmpcor.param.masterStartPixelAcross0 = a + self.c_cuAmpcor.param.masterStartPixelAcross0 = a @property def grossOffsetDownStatic(self): return self.c_cuAmpcor.param.grossOffsetDown0 @grossOffsetDownStatic.setter def grossOffsetDownStatic(self, int a): - self.c_cuAmpcor.param.grossOffsetDown0 =a + self.c_cuAmpcor.param.grossOffsetDown0 =a @property def grossOffsetAcrossStatic(self): return self.c_cuAmpcor.param.grossOffsetAcross0 @grossOffsetAcrossStatic.setter def grossOffsetAcrossStatic(self, int a): - self.c_cuAmpcor.param.grossOffsetAcross0 =a - + self.c_cuAmpcor.param.grossOffsetAcross0 =a + @property def grossOffsetDownDynamic(self): cdef int *c_data @@ -366,12 +382,12 @@ cdef class PyCuAmpcor(object): return p_data @grossOffsetDownDynamic.setter def grossOffsetDownDynamic (self, np.ndarray[np.int32_t,ndim=1,mode="c"] pa): - cdef int *c_data + cdef int *c_data cdef int *p_data c_data = self.c_cuAmpcor.param.grossOffsetDown p_data = pa.data for i in range (self.numberWindows): - c_data[i] = p_data[i] + c_data[i] = p_data[i] @property def grossOffsetAcrossDynamic(self): cdef int *c_data @@ -382,23 +398,23 @@ cdef class PyCuAmpcor(object): return p_data @grossOffsetAcrossDynamic.setter def grossOffsetAcrossDynamic (self, np.ndarray[np.int32_t,ndim=1,mode="c"] pa): - cdef int *c_data + cdef int *c_data cdef int *p_data c_data = self.c_cuAmpcor.param.grossOffsetAcross p_data = pa.data for i in range (self.numberWindows): - c_data[i] = p_data[i] + c_data[i] = p_data[i] return - + def setConstantGrossOffset(self, int goDown, int goAcross): - """ + """ constant gross offsets param goDown gross offset in azimuth direction param goAcross gross offset in range direction """ self.c_cuAmpcor.param.setStartPixels(self.masterStartPixelDownStatic, self.masterStartPixelAcrossStatic, goDown, goAcross) - + def setVaryingGrossOffset(self, np.ndarray[np.int32_t,ndim=1,mode="c"] vD, np.ndarray[np.int32_t,ndim=1,mode="c"] vA): """ varying gross offsets for each window @@ -411,21 +427,21 @@ cdef class PyCuAmpcor(object): def checkPixelInImageRange(self): """ check whether each window is with image range """ self.c_cuAmpcor.param.checkPixelInImageRange() - + def setupParams(self): """ set up constant parameters and allocate array parameters (offsets) should be called after number of windows is set and before setting varying gross offsets """ - self.c_cuAmpcor.param.setupParameters() + self.c_cuAmpcor.param.setupParameters() def runAmpcor(self): """ main procedure to run ampcor """ self.c_cuAmpcor.runAmpcor() - - - - - - + + + + + + diff --git a/contrib/PyCuAmpcor/src/SConscript b/contrib/PyCuAmpcor/src/SConscript index 1bb8d59..5de06c7 100644 --- a/contrib/PyCuAmpcor/src/SConscript +++ b/contrib/PyCuAmpcor/src/SConscript @@ -6,7 +6,7 @@ package = envPyCuAmpcor['PACKAGE'] project = envPyCuAmpcor['PROJECT'] build = envPyCuAmpcor['PRJ_LIB_DIR'] install = envPyCuAmpcor['PRJ_SCONS_INSTALL'] + '/' + package + '/' + project -listFiles = ['SlcImage.cu', 'cuArrays.cu', 'cuArraysCopy.cu', +listFiles = ['GDALImage.cu', 'cuArrays.cu', 'cuArraysCopy.cu', 'cuArraysPadding.cu', 'cuOverSampler.cu', 'cuSincOverSampler.cu', 'cuDeramp.cu', 'cuOffset.cu', 'cuCorrNormalization.cu', diff --git a/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu b/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu index 4f0b43d..ef911e7 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu +++ b/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu @@ -2,58 +2,74 @@ #include "cuAmpcorUtil.h" /** - * Run ampcor process for a batch of images (a chunk) + * 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] idxAcross_ index of the chunk along Across/Range direction - */ + */ void cuAmpcorChunk::run(int idxDown_, int idxAcross_) { // set chunk index setIndex(idxDown_, idxAcross_); - - // load master image chunk - loadMasterChunk(); - + + // load master image chunk + loadMasterChunk(); + //std::cout << "load master chunk ok\n"; - + cuArraysAbs(c_masterBatchRaw, r_masterBatchRaw, stream); cuArraysSubtractMean(r_masterBatchRaw, stream); // load slave image chunk loadSlaveChunk(); cuArraysAbs(c_slaveBatchRaw, r_slaveBatchRaw, stream); - + //std::cout << "load slave chunk ok\n"; - - + + //cross correlation for none-oversampled data if(param->algorithm == 0) { cuCorrFreqDomain->execute(r_masterBatchRaw, r_slaveBatchRaw, r_corrBatchRaw); } else { cuCorrTimeDomain(r_masterBatchRaw, r_slaveBatchRaw, r_corrBatchRaw, stream); //time domain cross correlation - } + } cuCorrNormalize(r_masterBatchRaw, r_slaveBatchRaw, r_corrBatchRaw, stream); - //find the maximum location of none-oversampled correlation - cuArraysMaxloc2D(r_corrBatchRaw, offsetInit, stream); -// Estimate SNR (Minyan Zhong) - //std::cout<< "flag stats 1" <outputToFile("offsetInit1", stream); - //std::cout<< "flag stats 3" <outputToFile("r_maxval",stream); + r_corrBatchRaw->outputToFile("r_corrBatchRaw",stream); + r_corrBatchRawZoomIn->outputToFile("r_corrBatchRawZoomIn",stream); + i_corrBatchZoomInValid->outputToFile("i_corrBatchZoomInValid",stream); + + // Summation of correlation and data point values + cuArraysSumCorr(r_corrBatchRawZoomIn, i_corrBatchZoomInValid, r_corrBatchSum, i_corrBatchValidCount, stream); + + // SNR + cuEstimateSnr(r_corrBatchSum, i_corrBatchValidCount, r_maxval, r_snrValue, stream); + + // Variance + // cuEstimateVariance(r_corrBatchRaw, offsetInit, r_maxval, r_covValue, stream); + + // Using the approximate estimation to adjust slave image (half search window size becomes only 4 pixels) //offsetInit->debuginfo(stream); // determine the starting pixel to extract slave images around the max location - cuDetermineSlaveExtractOffset(offsetInit, + cuDetermineSlaveExtractOffset(offsetInit, param->halfSearchRangeDownRaw, // old range - param->halfSearchRangeAcrossRaw, + param->halfSearchRangeAcrossRaw, param->halfZoomWindowSizeRaw, // new range param->halfZoomWindowSizeRaw, stream); @@ -63,58 +79,67 @@ void cuAmpcorChunk::run(int idxDown_, int idxAcross_) masterBatchOverSampler->execute(c_masterBatchRaw, c_masterBatchOverSampled, param->derampMethod); cuArraysAbs(c_masterBatchOverSampled, r_masterBatchOverSampled, stream); cuArraysSubtractMean(r_masterBatchOverSampled, stream); - + // extract slave and oversample cuArraysCopyExtract(c_slaveBatchRaw, c_slaveBatchZoomIn, offsetInit, stream); slaveBatchOverSampler->execute(c_slaveBatchZoomIn, c_slaveBatchOverSampled, param->derampMethod); cuArraysAbs(c_slaveBatchOverSampled, r_slaveBatchOverSampled, stream); - + // correlate oversampled images if(param->algorithm == 0) { cuCorrFreqDomain_OverSampled->execute(r_masterBatchOverSampled, r_slaveBatchOverSampled, r_corrBatchZoomIn); } else { - cuCorrTimeDomain(r_masterBatchOverSampled, r_slaveBatchOverSampled, r_corrBatchZoomIn, stream); - } + cuCorrTimeDomain(r_masterBatchOverSampled, r_slaveBatchOverSampled, r_corrBatchZoomIn, stream); + } cuCorrNormalize(r_masterBatchOverSampled, r_slaveBatchOverSampled, r_corrBatchZoomIn, stream); - + //std::cout << "debug correlation oversample\n"; //std::cout << r_masterBatchOverSampled->height << " " << r_masterBatchOverSampled->width << "\n"; //std::cout << r_slaveBatchOverSampled->height << " " << r_slaveBatchOverSampled->width << "\n"; //std::cout << r_corrBatchZoomIn->height << " " << r_corrBatchZoomIn->width << "\n"; - - // oversample the correlation surface + + // oversample the correlation surface cuArraysCopyExtract(r_corrBatchZoomIn, r_corrBatchZoomInAdjust, make_int2(0,0), stream); - + //std::cout << "debug oversampling " << r_corrBatchZoomInAdjust << " " << r_corrBatchZoomInOverSampled << "\n"; - + if(param->oversamplingMethod) { corrSincOverSampler->execute(r_corrBatchZoomInAdjust, r_corrBatchZoomInOverSampled); } else { - corrOverSampler->execute(r_corrBatchZoomInAdjust, r_corrBatchZoomInOverSampled); + corrOverSampler->execute(r_corrBatchZoomInAdjust, r_corrBatchZoomInOverSampled); } - + //find the max again - + cuArraysMaxloc2D(r_corrBatchZoomInOverSampled, offsetZoomIn, corrMaxValue, stream); - - // determine the final offset from non-oversampled (pixel) and oversampled (sub-pixel) - cuSubPixelOffset(offsetInit, offsetZoomIn, offsetFinal, - param->oversamplingFactor, param->rawDataOversamplingFactor, + + // determine the final offset from non-oversampled (pixel) and oversampled (sub-pixel) + cuSubPixelOffset(offsetInit, offsetZoomIn, offsetFinal, + param->oversamplingFactor, param->rawDataOversamplingFactor, param->halfSearchRangeDownRaw, param->halfSearchRangeAcrossRaw, param->halfZoomWindowSizeRaw, param->halfZoomWindowSizeRaw, stream); //offsetInit->debuginfo(stream); //offsetZoomIn->debuginfo(stream); - //offsetFinal->debuginfo(stream); - + //offsetFinal->debuginfo(stream); + + // Do insertion. + // Offsetfields. cuArraysCopyInsert(offsetFinal, offsetImage, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); - // Minyan Zhong - //cuArraysCopyInsert(corrMaxValue, snrImage, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); - //cuArraysCopyInsert(r_snrValue, snrImage, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); + // Debugging matrix. + cuArraysCopyInsert(r_corrBatchSum, floatImage1, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); + cuArraysCopyInsert(i_corrBatchValidCount, intImage1, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); + // Old: save max correlation coefficients. + //cuArraysCopyInsert(corrMaxValue, snrImage, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); + // New: save SNR + cuArraysCopyInsert(r_snrValue, snrImage, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); + + // Variance. + cuArraysCopyInsert(r_covValue, covImage, idxDown_*param->numberWindowDownInChunk, idxAcross_*param->numberWindowAcrossInChunk,stream); } void cuAmpcorChunk::setIndex(int idxDown_, int idxAcross_) @@ -122,14 +147,14 @@ void cuAmpcorChunk::setIndex(int idxDown_, int idxAcross_) idxChunkDown = idxDown_; idxChunkAcross = idxAcross_; idxChunk = idxChunkAcross + idxChunkDown*param->numberChunkAcross; - + if(idxChunkDown == param->numberChunkDown -1) { nWindowsDown = param->numberWindowDown - param->numberWindowDownInChunk*(param->numberChunkDown -1); } else { nWindowsDown = param->numberWindowDownInChunk; } - + if(idxChunkAcross == param->numberChunkAcross -1) { nWindowsAcross = param->numberWindowAcross - param->numberWindowAcrossInChunk*(param->numberChunkAcross -1); } @@ -137,20 +162,20 @@ void cuAmpcorChunk::setIndex(int idxDown_, int idxAcross_) nWindowsAcross = param->numberWindowAcrossInChunk; } //std::cout << "DEBUG setIndex" << idxChunk << " " << nWindowsDown << " " << nWindowsAcross << "\n"; - + } /// obtain the starting pixels for each chip -/// @param[in] oStartPixel +/// @param[in] oStartPixel /// void cuAmpcorChunk::getRelativeOffset(int *rStartPixel, const int *oStartPixel, int diff) { for(int i=0; inumberWindowDownInChunk; ++i) { int iDown = i; - if(i>=nWindowsDown) iDown = nWindowsDown-1; + if(i>=nWindowsDown) iDown = nWindowsDown-1; for(int j=0; jnumberWindowAcrossInChunk; ++j){ int iAcross = j; - if(j>=nWindowsAcross) iAcross = nWindowsAcross-1; + if(j>=nWindowsAcross) iAcross = nWindowsAcross-1; int idxInChunk = iDown*param->numberWindowAcrossInChunk+iAcross; int idxInAll = (iDown+idxChunkDown*param->numberWindowDownInChunk)*param->numberWindowAcross + idxChunkAcross*param->numberWindowAcrossInChunk+iAcross; @@ -158,108 +183,179 @@ void cuAmpcorChunk::getRelativeOffset(int *rStartPixel, const int *oStartPixel, //fprintf(stderr, "relative offset %d %d %d %d\n", i, j, rStartPixel[idxInChunk], diff); } } -} +} void cuAmpcorChunk::loadMasterChunk() { - //load a chunk from mmap to gpu - int startD = param->masterChunkStartPixelDown[idxChunk]; - int startA = param->masterChunkStartPixelAcross[idxChunk]; - int height = param->masterChunkHeight[idxChunk]; - int width = param->masterChunkWidth[idxChunk]; - masterImage->loadToDevice(c_masterChunkRaw->devData, startD, startA, height, width, stream); - std::cout << "debug load master: " << startD << " " << startA << " " << height << " " << width << "\n"; - //copy the chunk to a batch of images format (nImages, height, width) - //use cpu for some simple math + + // we first load the whole chunk of image from cpu to a gpu buffer c(r)_masterChunkRaw + // then copy to a batch of windows with (nImages, height, width) (leading dimension on the right) + + // get the chunk size to be loaded to gpu + int startD = param->masterChunkStartPixelDown[idxChunk]; //start pixel down (along height) + int startA = param->masterChunkStartPixelAcross[idxChunk]; // start pixel across (along width) + int height = param->masterChunkHeight[idxChunk]; // number of pixels along height + int width = param->masterChunkWidth[idxChunk]; // number of pixels along width + + //use cpu to compute the starting positions for each window getRelativeOffset(ChunkOffsetDown->hostData, param->masterStartPixelDown, param->masterChunkStartPixelDown[idxChunk]); + // copy the positions to gpu ChunkOffsetDown->copyToDevice(stream); + // same for the across direction getRelativeOffset(ChunkOffsetAcross->hostData, param->masterStartPixelAcross, param->masterChunkStartPixelAcross[idxChunk]); ChunkOffsetAcross->copyToDevice(stream); - // if derampMethod = 0 (no deramp), take amplitudes; otherwise, copy complex data - if(param->derampMethod == 0) { - cuArraysCopyToBatchAbsWithOffset(c_masterChunkRaw, param->masterChunkWidth[idxChunk], - c_masterBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + + // check whether the image is complex (e.g., SLC) or real( e.g. TIFF) + if(masterImage->isComplex()) + { + // allocate a gpu buffer to load data from cpu/file + // try allocate/deallocate the buffer on the fly to save gpu memory 07/09/19 + c_masterChunkRaw = new cuArrays (param->maxMasterChunkHeight, param->maxMasterChunkWidth); + c_masterChunkRaw->allocate(); + + // load the data from cpu + masterImage->loadToDevice((void *)c_masterChunkRaw->devData, startD, startA, height, width, stream); + //std::cout << "debug load master: " << startD << " " << startA << " " << height << " " << width << "\n"; + + //copy the chunk to a batch format (nImages, height, width) + // if derampMethod = 0 (no deramp), take amplitudes; otherwise, copy complex data + if(param->derampMethod == 0) { + cuArraysCopyToBatchAbsWithOffset(c_masterChunkRaw, param->masterChunkWidth[idxChunk], + c_masterBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + } + else { + cuArraysCopyToBatchWithOffset(c_masterChunkRaw, param->masterChunkWidth[idxChunk], + c_masterBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + } + // deallocate the gpu buffer + delete c_masterChunkRaw; } + // if the image is real else { - cuArraysCopyToBatchWithOffset(c_masterChunkRaw, param->masterChunkWidth[idxChunk], - c_masterBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + r_masterChunkRaw = new cuArrays (param->maxMasterChunkHeight, param->maxMasterChunkWidth); + r_masterChunkRaw->allocate(); + + // load the data from cpu + masterImage->loadToDevice((void *)r_masterChunkRaw->devData, startD, startA, height, width, stream); + + // copy the chunk (real) to a batch format (complex) + cuArraysCopyToBatchWithOffsetR2C(r_masterChunkRaw, param->masterChunkWidth[idxChunk], + c_masterBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + // deallocate the gpu buffer + delete r_masterChunkRaw; } + + } void cuAmpcorChunk::loadSlaveChunk() { - //load a chunk from mmap to gpu - slaveImage->loadToDevice(c_slaveChunkRaw->devData, - param->slaveChunkStartPixelDown[idxChunk], - param->slaveChunkStartPixelAcross[idxChunk], - param->slaveChunkHeight[idxChunk], - param->slaveChunkWidth[idxChunk], - stream); + //copy to a batch format (nImages, height, width) getRelativeOffset(ChunkOffsetDown->hostData, param->slaveStartPixelDown, param->slaveChunkStartPixelDown[idxChunk]); ChunkOffsetDown->copyToDevice(stream); getRelativeOffset(ChunkOffsetAcross->hostData, param->slaveStartPixelAcross, param->slaveChunkStartPixelAcross[idxChunk]); ChunkOffsetAcross->copyToDevice(stream); - if(param->derampMethod == 0) { - cuArraysCopyToBatchAbsWithOffset(c_slaveChunkRaw, param->slaveChunkWidth[idxChunk], - c_slaveBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); - } - else + + if(slaveImage->isComplex()) { - cuArraysCopyToBatchWithOffset(c_slaveChunkRaw, param->slaveChunkWidth[idxChunk], - c_slaveBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); - } + c_slaveChunkRaw = new cuArrays (param->maxSlaveChunkHeight, param->maxSlaveChunkWidth); + c_slaveChunkRaw->allocate(); + + //load a chunk from mmap to gpu + slaveImage->loadToDevice(c_slaveChunkRaw->devData, + param->slaveChunkStartPixelDown[idxChunk], + param->slaveChunkStartPixelAcross[idxChunk], + param->slaveChunkHeight[idxChunk], + param->slaveChunkWidth[idxChunk], + stream); + + if(param->derampMethod == 0) { + cuArraysCopyToBatchAbsWithOffset(c_slaveChunkRaw, param->slaveChunkWidth[idxChunk], + c_slaveBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + } + else { + cuArraysCopyToBatchWithOffset(c_slaveChunkRaw, param->slaveChunkWidth[idxChunk], + c_slaveBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + } + delete c_slaveChunkRaw; + } + else { //real image + //allocate the gpu buffer + r_slaveChunkRaw = new cuArrays (param->maxSlaveChunkHeight, param->maxSlaveChunkWidth); + r_slaveChunkRaw->allocate(); + + //load a chunk from mmap to gpu + slaveImage->loadToDevice(r_slaveChunkRaw->devData, + param->slaveChunkStartPixelDown[idxChunk], + param->slaveChunkStartPixelAcross[idxChunk], + param->slaveChunkHeight[idxChunk], + param->slaveChunkWidth[idxChunk], + stream); + + // convert to the batch format + cuArraysCopyToBatchWithOffsetR2C(r_slaveChunkRaw, param->slaveChunkWidth[idxChunk], + c_slaveBatchRaw, ChunkOffsetDown->devData, ChunkOffsetAcross->devData, stream); + delete r_slaveChunkRaw; + } } -cuAmpcorChunk::cuAmpcorChunk(cuAmpcorParameter *param_, SlcImage *master_, SlcImage *slave_, - cuArrays *offsetImage_, cuArrays *snrImage_, cudaStream_t stream_) +cuAmpcorChunk::cuAmpcorChunk(cuAmpcorParameter *param_, GDALImage *master_, GDALImage *slave_, + cuArrays *offsetImage_, cuArrays *snrImage_, cuArrays *covImage_, cuArrays *intImage1_, cuArrays *floatImage1_, cudaStream_t stream_) + { param = param_; masterImage = master_; - slaveImage = slave_; + slaveImage = slave_; offsetImage = offsetImage_; snrImage = snrImage_; + covImage = covImage_; + + intImage1 = intImage1_; + floatImage1 = floatImage1_; + stream = stream_; - - std::cout << "debug Chunk creator " << param->maxMasterChunkHeight << " " << param->maxMasterChunkWidth << "\n"; - c_masterChunkRaw = new cuArrays (param->maxMasterChunkHeight, param->maxMasterChunkWidth); - c_masterChunkRaw->allocate(); - - c_slaveChunkRaw = new cuArrays (param->maxSlaveChunkHeight, param->maxSlaveChunkWidth); - c_slaveChunkRaw->allocate(); - + + // std::cout << "debug Chunk creator " << param->maxMasterChunkHeight << " " << param->maxMasterChunkWidth << "\n"; + // try allocate/deallocate on the fly to save gpu memory 07/09/19 + // c_masterChunkRaw = new cuArrays (param->maxMasterChunkHeight, param->maxMasterChunkWidth); + // c_masterChunkRaw->allocate(); + + // c_slaveChunkRaw = new cuArrays (param->maxSlaveChunkHeight, param->maxSlaveChunkWidth); + // c_slaveChunkRaw->allocate(); + ChunkOffsetDown = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); ChunkOffsetDown->allocate(); ChunkOffsetDown->allocateHost(); ChunkOffsetAcross = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); ChunkOffsetAcross->allocate(); ChunkOffsetAcross->allocateHost(); - + c_masterBatchRaw = new cuArrays ( param->windowSizeHeightRaw, param->windowSizeWidthRaw, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); c_masterBatchRaw->allocate(); - + c_slaveBatchRaw = new cuArrays ( param->searchWindowSizeHeightRaw, param->searchWindowSizeWidthRaw, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); c_slaveBatchRaw->allocate(); - + r_masterBatchRaw = new cuArrays ( param->windowSizeHeightRaw, param->windowSizeWidthRaw, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_masterBatchRaw->allocate(); - + r_slaveBatchRaw = new cuArrays ( param->searchWindowSizeHeightRaw, param->searchWindowSizeWidthRaw, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_slaveBatchRaw->allocate(); - + c_slaveBatchZoomIn = new cuArrays ( param->searchWindowSizeHeightRawZoomIn, param->searchWindowSizeWidthRawZoomIn, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); c_slaveBatchZoomIn->allocate(); - + c_masterBatchOverSampled = new cuArrays ( param->windowSizeHeight, param->windowSizeWidth, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); @@ -269,7 +365,7 @@ cuAmpcorChunk::cuAmpcorChunk(cuAmpcorParameter *param_, SlcImage *master_, SlcIm param->searchWindowSizeHeight, param->searchWindowSizeWidth, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); c_slaveBatchOverSampled->allocate(); - + r_masterBatchOverSampled = new cuArrays ( param->windowSizeHeight, param->windowSizeWidth, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); @@ -279,66 +375,114 @@ cuAmpcorChunk::cuAmpcorChunk(cuAmpcorParameter *param_, SlcImage *master_, SlcIm param->searchWindowSizeHeight, param->searchWindowSizeWidth, param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_slaveBatchOverSampled->allocate(); - + masterBatchOverSampler = new cuOverSamplerC2C( c_masterBatchRaw->height, c_masterBatchRaw->width, //orignal size - c_masterBatchOverSampled->height, c_masterBatchOverSampled->width, //oversampled size + c_masterBatchOverSampled->height, c_masterBatchOverSampled->width, //oversampled size c_masterBatchRaw->count, stream); - - slaveBatchOverSampler = new cuOverSamplerC2C(c_slaveBatchZoomIn->height, c_slaveBatchZoomIn->width, + + slaveBatchOverSampler = new cuOverSamplerC2C(c_slaveBatchZoomIn->height, c_slaveBatchZoomIn->width, c_slaveBatchOverSampled->height, c_slaveBatchOverSampled->width, c_slaveBatchRaw->count, stream); - + r_corrBatchRaw = new cuArrays ( - param->searchWindowSizeHeightRaw-param->windowSizeHeightRaw+1, - param->searchWindowSizeWidthRaw-param->windowSizeWidthRaw+1, - param->numberWindowDownInChunk, + param->searchWindowSizeHeightRaw-param->windowSizeHeightRaw+1, + param->searchWindowSizeWidthRaw-param->windowSizeWidthRaw+1, + param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_corrBatchRaw->allocate(); - + r_corrBatchZoomIn = new cuArrays ( - param->searchWindowSizeHeight - param->windowSizeHeight+1, - param->searchWindowSizeWidth - param->windowSizeWidth+1, - param->numberWindowDownInChunk, + param->searchWindowSizeHeight - param->windowSizeHeight+1, + param->searchWindowSizeWidth - param->windowSizeWidth+1, + param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_corrBatchZoomIn->allocate(); - + r_corrBatchZoomInAdjust = new cuArrays ( - param->searchWindowSizeHeight - param->windowSizeHeight, - param->searchWindowSizeWidth - param->windowSizeWidth, - param->numberWindowDownInChunk, + param->searchWindowSizeHeight - param->windowSizeHeight, + param->searchWindowSizeWidth - param->windowSizeWidth, + param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_corrBatchZoomInAdjust->allocate(); - - + + r_corrBatchZoomInOverSampled = new cuArrays ( - param->zoomWindowSize * param->oversamplingFactor, param->zoomWindowSize * param->oversamplingFactor, - param->numberWindowDownInChunk, + param->zoomWindowSize * param->oversamplingFactor, + param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); r_corrBatchZoomInOverSampled->allocate(); - + offsetInit = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); offsetInit->allocate(); - + offsetZoomIn = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); offsetZoomIn->allocate(); - + offsetFinal = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); offsetFinal->allocate(); - + corrMaxValue = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); corrMaxValue->allocate(); - + + + // new arrays due to snr estimation + std::cout<< "corrRawZoomInHeight: " << param->corrRawZoomInHeight << "\n"; + std::cout<< "corrRawZoomInWidth: " << param->corrRawZoomInWidth << "\n"; + + r_corrBatchRawZoomIn = new cuArrays ( + param->corrRawZoomInHeight, + param->corrRawZoomInWidth, + param->numberWindowDownInChunk, + param->numberWindowAcrossInChunk); + r_corrBatchRawZoomIn->allocate(); + + i_corrBatchZoomInValid = new cuArrays ( + param->corrRawZoomInHeight, + param->corrRawZoomInWidth, + param->numberWindowDownInChunk, + param->numberWindowAcrossInChunk); + i_corrBatchZoomInValid->allocate(); + + + r_corrBatchSum = new cuArrays ( + param->numberWindowDownInChunk, + param->numberWindowAcrossInChunk); + r_corrBatchSum->allocate(); + + i_corrBatchValidCount = new cuArrays ( + param->numberWindowDownInChunk, + param->numberWindowAcrossInChunk); + i_corrBatchValidCount->allocate(); + + i_maxloc = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); + + i_maxloc->allocate(); + + r_maxval = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); + + r_maxval->allocate(); + + r_snrValue = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); + + r_snrValue->allocate(); + + r_covValue = new cuArrays (param->numberWindowDownInChunk, param->numberWindowAcrossInChunk); + + r_covValue->allocate(); + + // end of new arrays + if(param->oversamplingMethod) { corrSincOverSampler = new cuSincOverSamplerR2R(param->zoomWindowSize, param->oversamplingFactor, stream); } - else { + else { corrOverSampler= new cuOverSamplerR2R(param->zoomWindowSize, param->zoomWindowSize, - (param->zoomWindowSize)*param->oversamplingFactor, + (param->zoomWindowSize)*param->oversamplingFactor, (param->zoomWindowSize)*param->oversamplingFactor, - param->numberWindowDownInChunk*param->numberWindowAcrossInChunk, - stream); - } + param->numberWindowDownInChunk*param->numberWindowAcrossInChunk, + stream); + } if(param->algorithm == 0) { cuCorrFreqDomain = new cuFreqCorrelator( param->searchWindowSizeHeightRaw, param->searchWindowSizeWidthRaw, @@ -347,10 +491,10 @@ cuAmpcorChunk::cuAmpcorChunk(cuAmpcorParameter *param_, SlcImage *master_, SlcIm cuCorrFreqDomain_OverSampled = new cuFreqCorrelator( param->searchWindowSizeHeight, param->searchWindowSizeWidth, param->numberWindowDownInChunk*param->numberWindowAcrossInChunk, - stream); + stream); } - + debugmsg("all objects in chunk are created ...\n"); diff --git a/contrib/PyCuAmpcor/src/cuAmpcorChunk.h b/contrib/PyCuAmpcor/src/cuAmpcorChunk.h index 0de0aed..9ca2766 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorChunk.h +++ b/contrib/PyCuAmpcor/src/cuAmpcorChunk.h @@ -1,4 +1,4 @@ -/* +/* * cuAmpcorChunk.h * Purpose: a group of chips processed at the same time */ @@ -6,7 +6,7 @@ #ifndef __CUAMPCORCHUNK_H #define __CUAMPCORCHUNK_H -#include "SlcImage.h" +#include "GDALImage.h" #include "cuArrays.h" #include "cuAmpcorParameter.h" #include "cuOverSampler.h" @@ -22,64 +22,81 @@ private: int nWindowsAcross; int devId; - cudaStream_t stream; - - SlcImage *masterImage; - SlcImage *slaveImage; + cudaStream_t stream; + + GDALImage *masterImage; + GDALImage *slaveImage; cuAmpcorParameter *param; cuArrays *offsetImage; cuArrays *snrImage; - - cuArrays * c_masterChunkRaw, * c_slaveChunkRaw; + cuArrays *covImage; + + // added for test + cuArrays *intImage1; + cuArrays *floatImage1; + + // gpu buffer + cuArrays * c_masterChunkRaw, * c_slaveChunkRaw; + cuArrays * r_masterChunkRaw, * r_slaveChunkRaw; + + // gpu windows raw data cuArrays * c_masterBatchRaw, * c_slaveBatchRaw, * c_slaveBatchZoomIn; cuArrays * r_masterBatchRaw, * r_slaveBatchRaw; - cuArrays * c_masterBatchOverSampled, * c_slaveBatchOverSampled; + + // gpu windows oversampled data + cuArrays * c_masterBatchOverSampled, * c_slaveBatchOverSampled; cuArrays * r_masterBatchOverSampled, * r_slaveBatchOverSampled; cuArrays * r_corrBatchRaw, * r_corrBatchZoomIn, * r_corrBatchZoomInOverSampled, * r_corrBatchZoomInAdjust; - + cuArrays *ChunkOffsetDown, *ChunkOffsetAcross; - + cuOverSamplerC2C *masterBatchOverSampler, *slaveBatchOverSampler; - + cuOverSamplerR2R *corrOverSampler; - cuSincOverSamplerR2R *corrSincOverSampler; - + cuSincOverSamplerR2R *corrSincOverSampler; + //for frequency domain cuFreqCorrelator *cuCorrFreqDomain, *cuCorrFreqDomain_OverSampled; - + cuArrays *offsetInit; cuArrays *offsetZoomIn; cuArrays *offsetFinal; + cuArrays *corrMaxValue; + + + //SNR estimation + + cuArrays *r_corrBatchRawZoomIn; + cuArrays *r_corrBatchSum; + cuArrays *i_corrBatchZoomInValid, *i_corrBatchValidCount; + + cuArrays *r_snrValue; - //corr statistics cuArrays *i_maxloc; cuArrays *r_maxval; - - cuArrays *r_corrBatchSum; - cuArrays *i_corrBatchZoomInValid, *i_corrBatchValidCount; - - cuArrays *corrMaxValue; - cuArrays *r_snrValue; - + + // Varince estimation. + cuArrays *r_covValue; + public: cuAmpcorChunk() {} //cuAmpcorChunk(cuAmpcorParameter *param_, SlcImage *master_, SlcImage *slave_); - + void setIndex(int idxDown_, int idxAcross_); + cuAmpcorChunk(cuAmpcorParameter *param_, GDALImage *master_, GDALImage *slave_, cuArrays *offsetImage_, + cuArrays *snrImage_, cuArrays *covImage_, cuArrays *intImage1_, cuArrays *floatImage1_, cudaStream_t stream_); + - cuAmpcorChunk(cuAmpcorParameter *param_, SlcImage *master_, SlcImage *slave_, cuArrays *offsetImage_, - cuArrays *snrImage_, cudaStream_t stream_); - void loadMasterChunk(); void loadSlaveChunk(); void getRelativeOffset(int *rStartPixel, const int *oStartPixel, int diff); - - ~cuAmpcorChunk(); - - void run(int, int); + + ~cuAmpcorChunk(); + + void run(int, int); }; -#endif +#endif diff --git a/contrib/PyCuAmpcor/src/cuAmpcorController.cu b/contrib/PyCuAmpcor/src/cuAmpcorController.cu index 4798716..e9a6c33 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorController.cu +++ b/contrib/PyCuAmpcor/src/cuAmpcorController.cu @@ -1,113 +1,142 @@ // Implementation of cuAmpcorController #include "cuAmpcorController.h" -#include "SlcImage.h" +#include "GDALImage.h" #include "cuArrays.h" #include "cudaUtil.h" #include "cuAmpcorChunk.h" #include "cuAmpcorUtil.h" #include -cuAmpcorController::cuAmpcorController() { param = new cuAmpcorParameter();} -cuAmpcorController::~cuAmpcorController() { delete param; } +cuAmpcorController::cuAmpcorController() { param = new cuAmpcorParameter();} +cuAmpcorController::~cuAmpcorController() { delete param; } -void cuAmpcorController::runAmpcor() { - +void cuAmpcorController::runAmpcor() { + + // set the gpu id param->deviceID = gpuDeviceInit(param->deviceID); - SlcImage *masterImage; - SlcImage *slaveImage; - + // initialize the gdal driver + GDALAllRegister(); + // master and slave images; use band=1 as default + // TODO: selecting band + GDALImage *masterImage = new GDALImage(param->masterImageName, 1, param->mmapSizeInGB); + GDALImage *slaveImage = new GDALImage(param->slaveImageName, 1, param->mmapSizeInGB); + cuArrays *offsetImage, *offsetImageRun; cuArrays *snrImage, *snrImageRun; - - -// cuArrays *floatImage; -// cuArrays *intImage; + cuArrays *covImage, *covImageRun; + + // For debugging. + cuArrays *intImage1; + cuArrays *floatImage1; + + int nWindowsDownRun = param->numberChunkDown * param->numberWindowDownInChunk; + int nWindowsAcrossRun = param->numberChunkAcross * param->numberWindowAcrossInChunk; - masterImage = new SlcImage(param->masterImageName, param->masterImageHeight, param->masterImageWidth, param->mmapSizeInGB); - slaveImage = new SlcImage(param->slaveImageName, param->slaveImageHeight, param->slaveImageWidth, param->mmapSizeInGB); - - int nWindowsDownRun = param->numberChunkDown*param->numberWindowDownInChunk; - int nWindowsAcrossRun = param->numberChunkAcross*param->numberWindowAcrossInChunk; - std::cout << "Debug " << nWindowsDownRun << " " << param->numberWindowDown << "\n"; - + offsetImageRun = new cuArrays(nWindowsDownRun, nWindowsAcrossRun); - snrImageRun = new cuArrays(nWindowsDownRun, nWindowsAcrossRun); offsetImageRun->allocate(); + + snrImageRun = new cuArrays(nWindowsDownRun, nWindowsAcrossRun); snrImageRun->allocate(); - + + covImageRun = new cuArrays(nWindowsDownRun, nWindowsAcrossRun); + covImageRun->allocate(); + + // intImage 1 and floatImage 1 are added for debugging issues + + intImage1 = new cuArrays(nWindowsDownRun, nWindowsAcrossRun); + intImage1->allocate(); + + floatImage1 = new cuArrays(nWindowsDownRun, nWindowsAcrossRun); + floatImage1->allocate(); + + // Offsetfields. offsetImage = new cuArrays(param->numberWindowDown, param->numberWindowAcross); - snrImage = new cuArrays(param->numberWindowDown, param->numberWindowAcross); offsetImage->allocate(); + + // SNR. + snrImage = new cuArrays(param->numberWindowDown, param->numberWindowAcross); snrImage->allocate(); -// Minyan Zhong -// floatImage = new cuArrays(param->numberWindowDown, param->numberWindowAcross); -// intImage = new cuArrays(param->numberWindowDown, param->numberWindowAcross); + // Variance. + covImage = new cuArrays(param->numberWindowDown, param->numberWindowAcross); + covImage->allocate(); -// floatImage->allocate(); -// intImage->allocate(); -// cudaStream_t streams[param->nStreams]; cuAmpcorChunk *chunk[param->nStreams]; - for(int ist=0; istnStreams; ist++) + for(int ist=0; istnStreams; ist++) { cudaStreamCreate(&streams[ist]); - chunk[ist]= new cuAmpcorChunk(param, masterImage, slaveImage, offsetImageRun, snrImageRun, streams[ist]); + chunk[ist]= new cuAmpcorChunk(param, masterImage, slaveImage, offsetImageRun, snrImageRun, covImageRun, intImage1, floatImage1, streams[ist]); + } - + int nChunksDown = param->numberChunkDown; - int nChunksAcross = param->numberChunkAcross; - + int nChunksAcross = param->numberChunkAcross; + std::cout << "Total number of windows (azimuth x range): " <numberWindowDown << " x " << param->numberWindowAcross << std::endl; std::cout << "to be processed in the number of chunks: " <nStreams) { //std::cout << "Processing chunk(" << i <<", " << j <<")" << std::endl; for(int ist = 0; istnStreams; ist++) - { + { if(j+ist < nChunksAcross) { - + chunk[ist]->run(i, j+ist); } - } + } } } - + cudaDeviceSynchronize(); - + + // Do extraction. cuArraysCopyExtract(offsetImageRun, offsetImage, make_int2(0,0), streams[0]); - cuArraysCopyExtract(snrImageRun, snrImage, make_int2(0,0), streams[0]); - + cuArraysCopyExtract(snrImageRun, snrImage, make_int2(0,0), streams[0]); + cuArraysCopyExtract(covImageRun, covImage, make_int2(0,0), streams[0]); + offsetImage->outputToFile(param->offsetImageName, streams[0]); snrImage->outputToFile(param->snrImageName, streams[0]); + covImage->outputToFile(param->covImageName, streams[0]); -// Minyan Zhong -// floatImage->allocate(); -// intImage->allocate(); -// + // Output debugging arrays. + intImage1->outputToFile("intImage1", streams[0]); + floatImage1->outputToFile("floatImage1", streams[0]); outputGrossOffsets(); + + // Delete arrays. delete offsetImage; delete snrImage; + delete covImage; + + delete intImage1; + delete floatImage1; + delete offsetImageRun; delete snrImageRun; + delete covImageRun; + for (int ist=0; istnStreams; ist++) delete chunk[ist]; + delete masterImage; - delete slaveImage; -} + delete slaveImage; + +} void cuAmpcorController::outputGrossOffsets() { cuArrays *grossOffsets = new cuArrays(param->numberWindowDown, param->numberWindowAcross); grossOffsets->allocateHost(); - + for(int i=0; i< param->numberWindows; i++) grossOffsets->hostData[i] = make_float2(param->grossOffsetDown[i], param->grossOffsetAcross[i]); grossOffsets->outputHostToFile(param->grossOffsetImageName); @@ -176,7 +205,7 @@ void cuAmpcorController::setGrossOffsets(int *in, int size) { param->grossOffsets = (int *)malloc(size*sizeof(int)); mempcpy(param->grossOffsets, in, size*sizeof(int)); fprintf(stderr, "copy grossOffsets %d\n", size); -} +} void cuAmpcorController::setOffsetImageName(std::string s) { param->offsetImageName = s; } void cuAmpcorController::setSNRImageName(std::string s) { param->snrImageName = s; } //void cuAmpcorController::setMargin(int n) { param->margin = n; } diff --git a/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu b/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu index 1f42340..d50f5d3 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu +++ b/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu @@ -1,6 +1,6 @@ /** * cuAmpcorParameter.cu - * Input parameters for ampcor + * Input parameters for ampcor */ #include "cuAmpcorParameter.h" @@ -11,17 +11,19 @@ #endif /// -/// Constructor for cuAmpcorParameter class +/// Constructor for cuAmpcorParameter class /// also sets the default/initial values of various parameters /// - + cuAmpcorParameter::cuAmpcorParameter() { - algorithm = 0; //0 freq; 1 time - deviceID = 0; - nStreams = 1; + // default settings + // will be changed if they are set by python scripts + algorithm = 0; //0 freq; 1 time + deviceID = 0; + nStreams = 1; derampMethod = 1; - + windowSizeWidthRaw = 64; windowSizeHeightRaw = 64; halfSearchRangeDownRaw = 20; @@ -31,9 +33,9 @@ cuAmpcorParameter::cuAmpcorParameter() skipSampleDownRaw = 64; rawDataOversamplingFactor = 2; zoomWindowSize = 8; - oversamplingFactor = 16; - oversamplingMethod = 0; - + oversamplingFactor = 16; + oversamplingMethod = 0; + masterImageName = "master.slc"; masterImageWidth = 1000; masterImageHeight = 1000; @@ -43,50 +45,58 @@ cuAmpcorParameter::cuAmpcorParameter() offsetImageName = "DenseOffset.off"; grossOffsetImageName = "GrossOffset.off"; snrImageName = "snr.snr"; + covImageName = "cov.cov"; numberWindowDown = 1; - numberWindowAcross = 1; + numberWindowAcross = 1; numberWindowDownInChunk = 1; - numberWindowAcrossInChunk = 1 ; - + numberWindowAcrossInChunk = 1 ; + masterStartPixelDown0 = 0; masterStartPixelAcross0 = 0; + + corrRawZoomInHeight = 17; // 8*2+1 + corrRawZoomInWidth = 17; + + useMmap = 1; // use mmap + mmapSizeInGB = 1; + } /** - * To determine other process parameters after reading essential parameters from python - */ + * To determine other process parameters after reading essential parameters from python + */ void cuAmpcorParameter::setupParameters() -{ +{ zoomWindowSize *= rawDataOversamplingFactor; //8 * 2 - halfZoomWindowSizeRaw = zoomWindowSize/(2*rawDataOversamplingFactor); // 8*2/(2*2) = 4 - + halfZoomWindowSizeRaw = zoomWindowSize/(2*rawDataOversamplingFactor); // 8*2/(2*2) = 4 + windowSizeWidth = windowSizeWidthRaw*rawDataOversamplingFactor; // windowSizeHeight = windowSizeHeightRaw*rawDataOversamplingFactor; - searchWindowSizeWidthRaw = windowSizeWidthRaw + 2*halfSearchRangeDownRaw; + searchWindowSizeWidthRaw = windowSizeWidthRaw + 2*halfSearchRangeDownRaw; searchWindowSizeHeightRaw = windowSizeHeightRaw + 2*halfSearchRangeAcrossRaw; - + searchWindowSizeWidthRawZoomIn = windowSizeWidthRaw + 2*halfZoomWindowSizeRaw; searchWindowSizeHeightRawZoomIn = windowSizeHeightRaw + 2*halfZoomWindowSizeRaw; - + searchWindowSizeWidth = searchWindowSizeWidthRawZoomIn*rawDataOversamplingFactor; searchWindowSizeHeight = searchWindowSizeHeightRawZoomIn*rawDataOversamplingFactor; - + numberWindows = numberWindowDown*numberWindowAcross; if(numberWindows <=0) { fprintf(stderr, "Incorrect number of windows! (%d, %d)\n", numberWindowDown, numberWindowAcross); exit(EXIT_FAILURE); - } - + } + // modified 02/12/2018 to include one more chunk // e.g. numberWindowDownInChunk=102, numberWindowDown=10, results in numberChunkDown=11 - // the last chunk will include 2 windows, numberWindowDownInChunkRun = 2. - + // the last chunk will include 2 windows, numberWindowDownInChunkRun = 2. + numberChunkDown = IDIVUP(numberWindowDown, numberWindowDownInChunk); numberChunkAcross = IDIVUP(numberWindowAcross, numberWindowAcrossInChunk); numberChunks = numberChunkDown*numberChunkAcross; - allocateArrays(); + allocateArrays(); } @@ -99,7 +109,7 @@ void cuAmpcorParameter::allocateArrays() masterStartPixelAcross = (int *)malloc(arraySize); slaveStartPixelDown = (int *)malloc(arraySize); slaveStartPixelAcross = (int *)malloc(arraySize); - + int arraySizeChunk = numberChunks*sizeof(int); masterChunkStartPixelDown = (int *)malloc(arraySizeChunk); masterChunkStartPixelAcross = (int *)malloc(arraySizeChunk); @@ -130,18 +140,18 @@ void cuAmpcorParameter::deallocateArrays() } -/// Set starting pixels for master and slave windows from arrays +/// Set starting pixels for master and slave windows from arrays /// set also gross offsets between master and slave windows -/// +/// void cuAmpcorParameter::setStartPixels(int *mStartD, int *mStartA, int *gOffsetD, int *gOffsetA) { for(int i=0; i vpixel) mChunkSD = vpixel; + if(mChunkSD > vpixel) mChunkSD = vpixel; if(mChunkED < vpixel) mChunkED = vpixel; vpixel = masterStartPixelAcross[idxWindow]; - if(mChunkSA > vpixel) mChunkSA = vpixel; + if(mChunkSA > vpixel) mChunkSA = vpixel; if(mChunkEA < vpixel) mChunkEA = vpixel; vpixel = slaveStartPixelDown[idxWindow]; - if(sChunkSD > vpixel) sChunkSD = vpixel; + if(sChunkSD > vpixel) sChunkSD = vpixel; if(sChunkED < vpixel) sChunkED = vpixel; vpixel = slaveStartPixelAcross[idxWindow]; - if(sChunkSA > vpixel) sChunkSA = vpixel; + if(sChunkSA > vpixel) sChunkSA = vpixel; if(sChunkEA < vpixel) sChunkEA = vpixel; } } @@ -261,58 +271,58 @@ void cuAmpcorParameter::checkPixelInImageRange() for(int col = 0; col < numberWindowAcross; col++) { int i = row*numberWindowAcross + col; - if(masterStartPixelDown[i] <0) + if(masterStartPixelDown[i] <0) { fprintf(stderr, "Master Window start pixel out ot range in Down, window (%d,%d), pixel %d\n", row, col, masterStartPixelDown[i]); exit(EXIT_FAILURE); //or raise range error - } - if(masterStartPixelAcross[i] <0) + } + if(masterStartPixelAcross[i] <0) { fprintf(stderr, "Master Window start pixel out ot range in Across, window (%d,%d), pixel %d\n", row, col, masterStartPixelAcross[i]); exit(EXIT_FAILURE); - } - endPixel = masterStartPixelDown[i] + windowSizeHeightRaw; - if(endPixel >= masterImageHeight) + } + endPixel = masterStartPixelDown[i] + windowSizeHeightRaw; + if(endPixel >= masterImageHeight) { fprintf(stderr, "Master Window end pixel out ot range in Down, window (%d,%d), pixel %d\n", row, col, endPixel); exit(EXIT_FAILURE); - } - endPixel = masterStartPixelAcross[i] + windowSizeWidthRaw; - if(endPixel >= masterImageWidth) + } + endPixel = masterStartPixelAcross[i] + windowSizeWidthRaw; + if(endPixel >= masterImageWidth) { fprintf(stderr, "Master Window end pixel out ot range in Across, window (%d,%d), pixel %d\n", row, col, endPixel); exit(EXIT_FAILURE); - } + } //slave - if(slaveStartPixelDown[i] <0) + if(slaveStartPixelDown[i] <0) { fprintf(stderr, "Slave Window start pixel out ot range in Down, window (%d,%d), pixel %d\n", row, col, slaveStartPixelDown[i]); - exit(EXIT_FAILURE); - } - if(slaveStartPixelAcross[i] <0) + exit(EXIT_FAILURE); + } + if(slaveStartPixelAcross[i] <0) { fprintf(stderr, "Slave Window start pixel out ot range in Across, window (%d,%d), pixel %d\n", row, col, slaveStartPixelAcross[i]); exit(EXIT_FAILURE); - } - endPixel = slaveStartPixelDown[i] + searchWindowSizeHeightRaw; - if(endPixel >= slaveImageHeight) + } + endPixel = slaveStartPixelDown[i] + searchWindowSizeHeightRaw; + if(endPixel >= slaveImageHeight) { fprintf(stderr, "Slave Window end pixel out ot range in Down, window (%d,%d), pixel %d\n", row, col, endPixel); exit(EXIT_FAILURE); - } - endPixel = slaveStartPixelAcross[i] + searchWindowSizeWidthRaw; - if(endPixel >= slaveImageWidth) + } + endPixel = slaveStartPixelAcross[i] + searchWindowSizeWidthRaw; + if(endPixel >= slaveImageWidth) { fprintf(stderr, "Slave Window end pixel out ot range in Across, window (%d,%d), pixel %d\n", row, col, endPixel); exit(EXIT_FAILURE); - } - - } + } + + } } } -cuAmpcorParameter::~cuAmpcorParameter() +cuAmpcorParameter::~cuAmpcorParameter() { deallocateArrays(); } diff --git a/contrib/PyCuAmpcor/src/cuAmpcorParameter.h b/contrib/PyCuAmpcor/src/cuAmpcorParameter.h index 6c258b7..eed1027 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorParameter.h +++ b/contrib/PyCuAmpcor/src/cuAmpcorParameter.h @@ -1,7 +1,7 @@ /** * cuAmpcorParameter.h * Header file for Ampcor Parameter Class - * + * * Author: Lijun Zhu @ Seismo Lab, Caltech * March 2017 */ @@ -12,13 +12,13 @@ #include /// Class container for all parameters -/// +/// /// @note -/// The dimension/direction names used are: +/// The dimension/direction names used are: /// The inner-most dimension: x, row, height, down, azimuth, along the track. /// The outer-most dimension: y, column, width, across, range, along the sight. -/// C/C++/Python use row-major indexing: a[i][j] -> a[i*WIDTH+j] -/// FORTRAN/BLAS/CUBLAS use column-major indexing: a[i][j]->a[i+j*LENGTH] +/// C/C++/Python use row-major indexing: a[i][j] -> a[i*WIDTH+j] +/// FORTRAN/BLAS/CUBLAS use column-major indexing: a[i][j]->a[i+j*LENGTH] /// @note /// Common procedures to use cuAmpcorParameter @@ -27,72 +27,74 @@ /// 3. Call setupParameters() to determine related parameters and allocate starting pixels for each window: param->setupParameters() /// 4. Provide/set Master window starting pixel(s), and gross offset(s): param->setStartPixels(masterStartDown, masterStartAcross, grossOffsetDown, grossOffsetAcross) /// 4a. Optionally, check the range of windows is within the SLC image range: param->checkPixelInImageRange() -/// Steps 1, 3, 4 are mandatory. If step 2 is missing, default values will be used +/// Steps 1, 3, 4 are mandatory. If step 2 is missing, default values will be used class cuAmpcorParameter{ public: - int algorithm; /// Cross-correlation algorithm: 0=freq domain (default) 1=time domain - int deviceID; /// Targeted GPU device ID: use -1 to auto select - int nStreams; /// Number of streams to asynchonize data transfers and compute kernels + int algorithm; /// Cross-correlation algorithm: 0=freq domain (default) 1=time domain + int deviceID; /// Targeted GPU device ID: use -1 to auto select + int nStreams; /// Number of streams to asynchonize data transfers and compute kernels int derampMethod; /// Method for deramping 0=None, 1=average, 2=phase gradient - + // chip or window size for raw data int windowSizeHeightRaw; /// Template window height (original size) - int windowSizeWidthRaw; /// Template window width (original size) - int searchWindowSizeHeightRaw; /// Search window height (original size) + int windowSizeWidthRaw; /// Template window width (original size) + int searchWindowSizeHeightRaw; /// Search window height (original size) int searchWindowSizeWidthRaw; /// Search window width (orignal size) - - int halfSearchRangeDownRaw; ///(searchWindowSizeHeightRaw-windowSizeHeightRaw)/2 + + int halfSearchRangeDownRaw; ///(searchWindowSizeHeightRaw-windowSizeHeightRaw)/2 int halfSearchRangeAcrossRaw; ///(searchWindowSizeWidthRaw-windowSizeWidthRaw)/2 // search range is (-halfSearchRangeRaw, halfSearchRangeRaw) - + int searchWindowSizeHeightRawZoomIn; int searchWindowSizeWidthRawZoomIn; - - + + int corrRawZoomInHeight; // window to estimate snr + int corrRawZoomInWidth; + // chip or window size after oversampling int rawDataOversamplingFactor; /// Raw data overampling factor (from original size to oversampled size) int windowSizeHeight; /// Template window length (oversampled size) int windowSizeWidth; /// Template window width (original size) - int searchWindowSizeHeight; /// Search window height (oversampled size) - int searchWindowSizeWidth; /// Search window width (oversampled size) - - // strides between chips/windows + int searchWindowSizeHeight; /// Search window height (oversampled size) + int searchWindowSizeWidth; /// Search window width (oversampled size) + + // strides between chips/windows int skipSampleDownRaw; /// Skip size between neighboring windows in Down direction (original size) int skipSampleAcrossRaw; /// Skip size between neighboring windows in across direction (original size) //int skipSampleDown; /// Skip size between neighboring windows in Down direction (oversampled size) //int skipSampleAcross; /// Skip size between neighboring windows in Across direction (oversampled size) - + // Zoom in region near location of max correlation - int zoomWindowSize; /// Zoom-in window size in correlation surface (same for down and across directions) - int halfZoomWindowSizeRaw; /// = half of zoomWindowSize/rawDataOversamplingFactor - - - + int zoomWindowSize; /// Zoom-in window size in correlation surface (same for down and across directions) + int halfZoomWindowSizeRaw; /// = half of zoomWindowSize/rawDataOversamplingFactor + + + int oversamplingFactor; /// Oversampling factor for interpolating correlation surface - int oversamplingMethod; /// 0 = fft (default) 1 = sinc - - - float thresholdSNR; /// Threshold of Signal noise ratio to remove noisy data - + int oversamplingMethod; /// 0 = fft (default) 1 = sinc + + + float thresholdSNR; /// Threshold of Signal noise ratio to remove noisy data + //master image std::string masterImageName; /// master SLC image name int imageDataType1; /// master image data type, 2=cfloat=complex=float2 1=float - int masterImageHeight; /// master image height + int masterImageHeight; /// master image height int masterImageWidth; /// master image width - + //slave image std::string slaveImageName; /// slave SLC image name int imageDataType2; /// slave image data type, 2=cfloat=complex=float2 1=float - int slaveImageHeight; /// slave image height + int slaveImageHeight; /// slave image height int slaveImageWidth; /// slave image width - + // total number of chips/windows int numberWindowDown; /// number of total windows (down) int numberWindowAcross; /// number of total windows (across) int numberWindows; /// numberWindowDown*numberWindowAcross - + // number of chips/windows in a batch/chunk int numberWindowDownInChunk; /// number of windows processed in a chunk (down) int numberWindowAcrossInChunk; /// number of windows processed in a chunk (across) @@ -100,20 +102,21 @@ public: int numberChunkDown; /// number of chunks (down) int numberChunkAcross; /// number of chunks (across) int numberChunks; - - int mmapSizeInGB; - + + int useMmap; /// whether to use mmap 0=not 1=yes (default = 0) + int mmapSizeInGB; /// size for mmap buffer(useMmap=1) or a cpu memory buffer (useMmap=0) + int masterStartPixelDown0; int masterStartPixelAcross0; - int *masterStartPixelDown; /// master starting pixels for each window (down) + int *masterStartPixelDown; /// master starting pixels for each window (down) int *masterStartPixelAcross;/// master starting pixels for each window (across) - int *slaveStartPixelDown; /// slave starting pixels for each window (down) - int *slaveStartPixelAcross; /// slave starting pixels for each window (across) + int *slaveStartPixelDown; /// slave starting pixels for each window (down) + int *slaveStartPixelAcross; /// slave starting pixels for each window (across) int grossOffsetDown0; int grossOffsetAcross0; int *grossOffsetDown; /// Gross offsets between master and slave windows (down) : slaveStartPixel - masterStartPixel - int *grossOffsetAcross; /// Gross offsets between master and slave windows (across) - + int *grossOffsetAcross; /// Gross offsets between master and slave windows (across) + int *masterChunkStartPixelDown; int *masterChunkStartPixelAcross; int *slaveChunkStartPixelDown; @@ -124,18 +127,19 @@ public: int *slaveChunkWidth; int maxMasterChunkHeight, maxMasterChunkWidth; int maxSlaveChunkHeight, maxSlaveChunkWidth; - - std::string grossOffsetImageName; - std::string offsetImageName; /// Output Offset fields filename + + std::string grossOffsetImageName; + std::string offsetImageName; /// Output Offset fields filename std::string snrImageName; /// Output SNR filename - + std::string covImageName; + cuAmpcorParameter(); /// Class constructor and default parameters setter - ~cuAmpcorParameter(); /// Class descontructor - + ~cuAmpcorParameter(); /// Class descontructor + void allocateArrays(); /// Allocate various arrays after the number of Windows is given void deallocateArrays(); /// Deallocate arrays on exit - + /// Three methods to set master/slave starting pixels and gross offsets from input master start pixel(s) and gross offset(s) /// 1 (int *, int *, int *, int *): varying master start pixels and gross offsets /// 2 (int, int, int *, int *): fixed master start pixel (first window) and varying gross offsets @@ -144,7 +148,7 @@ public: void setStartPixels(int, int, int*, int*); void setStartPixels(int, int, int, int); void setChunkStartPixels(); - void checkPixelInImageRange(); /// check whether + void checkPixelInImageRange(); /// check whether void setupParameters(); /// Process other parameters after Python Input }; diff --git a/contrib/PyCuAmpcor/src/cuAmpcorUtil.h b/contrib/PyCuAmpcor/src/cuAmpcorUtil.h index a8c8a2e..fe20b19 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorUtil.h +++ b/contrib/PyCuAmpcor/src/cuAmpcorUtil.h @@ -1,10 +1,10 @@ -/* +/* * cuAmpcorUtil.h * header file to include the various routines for ampcor - * serves as an index + * serves as an index */ - - + + #ifndef __CUAMPCORUTIL_H #define __CUMAPCORUTIL_H @@ -18,20 +18,27 @@ //in cuArraysCopy.cu: various utitlies for copy images file in gpu memory void cuArraysCopyToBatch(cuArrays *image1, cuArrays *image2, int strideH, int strideW, cudaStream_t stream); -void cuArraysCopyToBatchWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, +void cuArraysCopyToBatchWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, const int *offsetH, const int* offsetW, cudaStream_t stream); -void cuArraysCopyToBatchAbsWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, +void cuArraysCopyToBatchAbsWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, + const int *offsetH, const int* offsetW, cudaStream_t stream); +void cuArraysCopyToBatchWithOffsetR2C(cuArrays *image1, const int lda1, cuArrays *image2, const int *offsetH, const int* offsetW, cudaStream_t stream); void cuArraysCopyC2R(cuArrays *image1, cuArrays *image2, int strideH, int strideW, cudaStream_t stream); +// same routine name overloaded for different data type void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, cuArrays *offset, cudaStream_t stream); void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, int2 offset, cudaStream_t stream); void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, int2 offset, cudaStream_t stream); void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, int2 offset, cudaStream_t stream); void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, cuArrays *offsets, cudaStream_t stream); +void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, int2 offset, cudaStream_t stream); void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int offsetX, int offersetY, cudaStream_t stream); +void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int offsetX, int offersetY, cudaStream_t stream); void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int offsetX, int offsetY, cudaStream_t stream); +void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int offsetX, int offersetY, cudaStream_t stream); + void cuArraysCopyInversePadded(cuArrays *imageIn, cuArrays *imageOut,cudaStream_t stream); void cuArraysCopyPadded(cuArrays *imageIn, cuArrays *imageOut,cudaStream_t stream); @@ -46,8 +53,8 @@ void cuDerampMethod2(cuArrays *images, cudaStream_t stream); void cpuDerampMethod3(cuArrays *images, cudaStream_t stream); //in cuArraysPadding.cu: various utilities for oversampling padding -void cuArraysPadding(cuArrays *image1, cuArrays *image2, cudaStream_t stream); -void cuArraysPaddingMany(cuArrays *image1, cuArrays *image2, cudaStream_t stream); +void cuArraysPadding(cuArrays *image1, cuArrays *image2, cudaStream_t stream); +void cuArraysPaddingMany(cuArrays *image1, cuArrays *image2, cudaStream_t stream); void cuArraysR2C(cuArrays *image1, cuArrays *image2, cudaStream_t stream); void cuArraysC2R(cuArrays *image1, cuArrays *image2, cudaStream_t stream); void cuArraysAbs(cuArrays *image1, cuArrays *image2, cudaStream_t stream); @@ -57,21 +64,21 @@ void cuArraysSubtractMean(cuArrays *images, cudaStream_t stream); void cuCorrNormalize(cuArrays *templates, cuArrays *images, cuArrays *results, cudaStream_t stream); //in cuOffset.cu: utitilies for determining the max locaiton of cross correlations or the offset -void cuArraysMaxloc2D(cuArrays *images, cuArrays *maxloc, cuArrays *maxval, cudaStream_t stream); -void cuArraysMaxloc2D(cuArrays *images, cuArrays *maxloc, cudaStream_t stream); +void cuArraysMaxloc2D(cuArrays *images, cuArrays *maxloc, cuArrays *maxval, cudaStream_t stream); +void cuArraysMaxloc2D(cuArrays *images, cuArrays *maxloc, cudaStream_t stream); void cuSubPixelOffset(cuArrays *offsetInit, cuArrays *offsetZoomIn, cuArrays *offsetFinal, int OverSampleRatioZoomin, int OverSampleRatioRaw, int xHalfRangeInit, int yHalfRangeInit, int xHalfRangeZoomIn, int yHalfRangeZoomIn, cudaStream_t stream); - + void cuDetermineInterpZone(cuArrays *maxloc, cuArrays *zoomInOffset, cuArrays *corrOrig, cuArrays *corrZoomIn, cudaStream_t stream); void cuDetermineSlaveExtractOffset(cuArrays *maxLoc, int xOldRange, int yOldRange, int xNewRange, int yNewRange, cudaStream_t stream); //in cuCorrTimeDomain.cu: cross correlation in time domain -void cuCorrTimeDomain(cuArrays *templates, cuArrays *images, cuArrays *results, cudaStream_t stream); +void cuCorrTimeDomain(cuArrays *templates, cuArrays *images, cuArrays *results, cudaStream_t stream); //in cuCorrFrequency.cu: cross correlation in freq domain, also include fft correlatior class -void cuArraysElementMultiply(cuArrays *image1, cuArrays *image2, cudaStream_t stream); +void cuArraysElementMultiply(cuArrays *image1, cuArrays *image2, cudaStream_t stream); void cuArraysElementMultiplyConjugate(cuArrays *image1, cuArrays *image2, float coef, cudaStream_t stream); @@ -80,7 +87,11 @@ void cuArraysElementMultiplyConjugate(cuArrays *image1, cuArrays void cuArraysCopyExtractCorr(cuArrays *imagesIn, cuArrays *imagesOut, cuArrays *imagesValid, cuArrays *maxloc, cudaStream_t stream); // implemented in cuCorrNormalization.cu void cuArraysSumCorr(cuArrays *images, cuArrays *imagesValid, cuArrays *imagesSum, cuArrays *imagesValidCount, cudaStream_t stream); + // implemented in cuEstimateStats.cu void cuEstimateSnr(cuArrays *corrSum, cuArrays *corrValidCount, cuArrays *maxval, cuArrays *snrValue, cudaStream_t stream); -#endif +// implemented in cuEstimateStats.cu +void cuEstimateVariance(cuArrays *corrBatchRaw, cuArrays *maxloc, cuArrays *maxval, cuArrays *covValue, cudaStream_t stream); + +#endif diff --git a/contrib/PyCuAmpcor/src/cuArrays.cu b/contrib/PyCuAmpcor/src/cuArrays.cu index 7417e3c..82dfe51 100644 --- a/contrib/PyCuAmpcor/src/cuArrays.cu +++ b/contrib/PyCuAmpcor/src/cuArrays.cu @@ -154,8 +154,21 @@ file.write((char *)data, size*count*sizeof(float2)); file.close(); } + + template<> + void cuArrays::outputToFile(std::string fn, cudaStream_t stream) + { + float *data; + data = (float *)malloc(size*count*sizeof(float3)); + checkCudaErrors(cudaMemcpyAsync(data, devData, size*count*sizeof(float3), cudaMemcpyDeviceToHost, stream)); + std::ofstream file; + file.open(fn.c_str(), std::ios_base::binary); + file.write((char *)data, size*count*sizeof(float3)); + file.close(); + } template class cuArrays; template class cuArrays; + template class cuArrays; template class cuArrays; template class cuArrays; diff --git a/contrib/PyCuAmpcor/src/cuArraysCopy.cu b/contrib/PyCuAmpcor/src/cuArraysCopy.cu index dd7a42f..822b3d2 100644 --- a/contrib/PyCuAmpcor/src/cuArraysCopy.cu +++ b/contrib/PyCuAmpcor/src/cuArraysCopy.cu @@ -4,7 +4,7 @@ * Lijun Zhu @ Seismo Lab, Caltech * v1.0 Jan 2017 */ - + #include "cuArrays.h" #include "cudaUtil.h" #include "cudaError.h" @@ -16,14 +16,14 @@ inline __device__ float cuAbs(float2 a) return sqrtf(a.x*a.x+a.y*a.y); }*/ -//copy a chunk into a series of chips -__global__ void cuArraysCopyToBatch_kernel(const float2 *imageIn, const int inNX, const int inNY, - float2 *imageOut, const int outNX, const int outNY, - const int nImagesX, const int nImagesY, +// copy a chunk into a batch of chips for a given stride +__global__ void cuArraysCopyToBatch_kernel(const float2 *imageIn, const int inNX, const int inNY, + float2 *imageOut, const int outNX, const int outNY, + const int nImagesX, const int nImagesY, const int strideX, const int strideY) { int idxImage = blockIdx.z; - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; if(idxImage >=nImagesX*nImagesY|| outx >= outNX || outy >= outNY) return; int idxOut = idxImage*outNX*outNY + outx*outNY + outy; @@ -33,8 +33,7 @@ __global__ void cuArraysCopyToBatch_kernel(const float2 *imageIn, const int inNX imageOut[idxOut] = imageIn[idxIn]; } -//tested -void cuArraysCopyToBatch(cuArrays *image1, cuArrays *image2, +void cuArraysCopyToBatch(cuArrays *image1, cuArrays *image2, int strideH, int strideW, cudaStream_t stream) { const int nthreads = NTHREADS2D; @@ -48,12 +47,14 @@ void cuArraysCopyToBatch(cuArrays *image1, cuArrays *image2, getLastCudaError("cuArraysCopyToBatch_kernel"); } -__global__ void cuArraysCopyToBatchWithOffset_kernel(const float2 *imageIn, const int inNY, - float2 *imageOut, const int outNX, const int outNY, const int nImages, + +// copy a chunk into a batch of chips for a set of offsets (varying strides), from complex to complex +__global__ void cuArraysCopyToBatchWithOffset_kernel(const float2 *imageIn, const int inNY, + float2 *imageOut, const int outNX, const int outNY, const int nImages, const int *offsetX, const int *offsetY) { int idxImage = blockIdx.z; - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; if(idxImage>=nImages || outx >= outNX || outy >= outNY) return; int idxOut = idxImage*outNX*outNY + outx*outNY + outy; @@ -61,11 +62,8 @@ __global__ void cuArraysCopyToBatchWithOffset_kernel(const float2 *imageIn, cons imageOut[idxOut] = imageIn[idxIn]; } -/// @param[in] image1 input image in a large chunk -/// @param[in] lda1 width of image 1 -/// @param[out] image2 output image with a batch of small windows - -void cuArraysCopyToBatchWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, +// lda1 (inNY) is the leading dimension of image1, usually, its width +void cuArraysCopyToBatchWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, const int *offsetH, const int* offsetW, cudaStream_t stream) { const int nthreads = 16; @@ -79,12 +77,13 @@ void cuArraysCopyToBatchWithOffset(cuArrays *image1, const int lda1, cuA getLastCudaError("cuArraysCopyToBatchAbsWithOffset_kernel"); } -__global__ void cuArraysCopyToBatchAbsWithOffset_kernel(const float2 *imageIn, const int inNY, - float2 *imageOut, const int outNX, const int outNY, const int nImages, +// copy a chunk into a batch of chips for a set of offsets (varying strides), from complex to real(take amplitudes) +__global__ void cuArraysCopyToBatchAbsWithOffset_kernel(const float2 *imageIn, const int inNY, + float2 *imageOut, const int outNX, const int outNY, const int nImages, const int *offsetX, const int *offsetY) { int idxImage = blockIdx.z; - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; if(idxImage>=nImages || outx >= outNX || outy >= outNY) return; int idxOut = idxImage*outNX*outNY + outx*outNY + outy; @@ -92,7 +91,7 @@ __global__ void cuArraysCopyToBatchAbsWithOffset_kernel(const float2 *imageIn, c imageOut[idxOut] = make_float2(complexAbs(imageIn[idxIn]), 0.0); } -void cuArraysCopyToBatchAbsWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, +void cuArraysCopyToBatchAbsWithOffset(cuArrays *image1, const int lda1, cuArrays *image2, const int *offsetH, const int* offsetW, cudaStream_t stream) { const int nthreads = 16; @@ -106,14 +105,42 @@ void cuArraysCopyToBatchAbsWithOffset(cuArrays *image1, const int lda1, getLastCudaError("cuArraysCopyToBatchAbsWithOffset_kernel"); } +// copy a chunk into a batch of chips for a set of offsets (varying strides), from real to complex(to real part) +__global__ void cuArraysCopyToBatchWithOffsetR2C_kernel(const float *imageIn, const int inNY, + float2 *imageOut, const int outNX, const int outNY, const int nImages, + const int *offsetX, const int *offsetY) +{ + int idxImage = blockIdx.z; + int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outy = threadIdx.y + blockDim.y*blockIdx.y; + if(idxImage>=nImages || outx >= outNX || outy >= outNY) return; + int idxOut = idxImage*outNX*outNY + outx*outNY + outy; + int idxIn = (offsetX[idxImage]+outx)*inNY + offsetY[idxImage] + outy; + imageOut[idxOut] = make_float2(imageIn[idxIn], 0.0f); +} + +void cuArraysCopyToBatchWithOffsetR2C(cuArrays *image1, const int lda1, cuArrays *image2, + const int *offsetH, const int* offsetW, cudaStream_t stream) +{ + const int nthreads = 16; + dim3 blockSize(nthreads, nthreads, 1); + dim3 gridSize(IDIVUP(image2->height,nthreads), IDIVUP(image2->width,nthreads), image2->count); + //fprintf(stderr, "copy tile to batch, %d %d\n", lda1, image2->count); + cuArraysCopyToBatchWithOffsetR2C_kernel<<>> ( + image1->devData, lda1, + image2->devData, image2->height, image2->width, image2->count, + offsetH, offsetW); + getLastCudaError("cuArraysCopyToBatchWithOffsetR2C_kernel"); +} + //copy a chunk into a series of chips -__global__ void cuArraysCopyC2R_kernel(const float2 *imageIn, const int inNX, const int inNY, - float *imageOut, const int outNX, const int outNY, - const int nImagesX, const int nImagesY, +__global__ void cuArraysCopyC2R_kernel(const float2 *imageIn, const int inNX, const int inNY, + float *imageOut, const int outNX, const int outNY, + const int nImagesX, const int nImagesY, const int strideX, const int strideY, const float factor) { int idxImage = blockIdx.z; - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; if(idxImage >=nImagesX*nImagesY|| outx >= outNX || outy >= outNY) return; int idxOut = idxImage*outNX*outNY + outx*outNY + outy; @@ -121,17 +148,17 @@ __global__ void cuArraysCopyC2R_kernel(const float2 *imageIn, const int inNX, co int idxImageY = idxImage%nImagesY; int idxIn = (idxImageX*strideX+outx)*inNY + idxImageY*strideY+outy; imageOut[idxOut] = complexAbs(imageIn[idxIn])*factor; - //printf( "%d\n", idxOut); + //printf( "%d\n", idxOut); } //tested -void cuArraysCopyC2R(cuArrays *image1, cuArrays *image2, +void cuArraysCopyC2R(cuArrays *image1, cuArrays *image2, int strideH, int strideW, cudaStream_t stream) { const int nthreads = 16; dim3 blockSize(nthreads, nthreads, 1); dim3 gridSize(IDIVUP(image2->height,nthreads), IDIVUP(image2->width,nthreads), image2->count); - float factor = 1.0f/image1->size; //the FFT factor + float factor = 1.0f/image1->size; //the FFT factor cuArraysCopyC2R_kernel<<>> ( image1->devData, image1->height, image1->width, image2->devData, image2->height, image2->width, @@ -141,22 +168,22 @@ void cuArraysCopyC2R(cuArrays *image1, cuArrays *image2, } __global__ void cuArraysCopyExtractVaryingOffset(const float *imageIn, const int inNX, const int inNY, - float *imageOut, const int outNX, const int outNY, const int nImages, + float *imageOut, const int outNX, const int outNY, const int nImages, const int2 *offsets) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) { - int idxImage = blockIdx.z; + int idxImage = blockIdx.z; int idxOut = (blockIdx.z * outNX + outx)*outNY+outy; - int idxIn = (blockIdx.z*inNX + outx + offsets[idxImage].x)*inNY + outy + offsets[idxImage].y; + int idxIn = (blockIdx.z*inNX + outx + offsets[idxImage].x)*inNY + outy + offsets[idxImage].y; imageOut[idxOut] = imageIn[idxIn]; } } -/* copy a tile of images to another image, with starting pixels offsets +/* copy a tile of images to another image, with starting pixels offsets * param[in] imageIn inut images * param[out] imageOut output images of dimension nImages*outNX*outNY */ @@ -166,24 +193,24 @@ void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, const int nthreads = 16; dim3 threadsperblock(nthreads, nthreads,1); dim3 blockspergrid(IDIVUP(imagesOut->height,nthreads), IDIVUP(imagesOut->width,nthreads), imagesOut->count); - cuArraysCopyExtractVaryingOffset<<>>(imagesIn->devData, imagesIn->height, imagesIn->width, + cuArraysCopyExtractVaryingOffset<<>>(imagesIn->devData, imagesIn->height, imagesIn->width, imagesOut->devData, imagesOut->height, imagesOut->width, imagesOut->count, offsets->devData); getLastCudaError("cuArraysCopyExtract error"); } __global__ void cuArraysCopyExtractVaryingOffset_C2C(const float2 *imageIn, const int inNX, const int inNY, - float2 *imageOut, const int outNX, const int outNY, const int nImages, + float2 *imageOut, const int outNX, const int outNY, const int nImages, const int2 *offsets) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) { - int idxImage = blockIdx.z; + int idxImage = blockIdx.z; int idxOut = (blockIdx.z * outNX + outx)*outNY+outy; - int idxIn = (blockIdx.z*inNX + outx + offsets[idxImage].x)*inNY + outy + offsets[idxImage].y; + int idxIn = (blockIdx.z*inNX + outx + offsets[idxImage].x)*inNY + outy + offsets[idxImage].y; imageOut[idxOut] = imageIn[idxIn]; } } @@ -194,7 +221,7 @@ void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut const int nthreads = 16; dim3 threadsperblock(nthreads, nthreads,1); dim3 blockspergrid(IDIVUP(imagesOut->height,nthreads), IDIVUP(imagesOut->width,nthreads), imagesOut->count); - cuArraysCopyExtractVaryingOffset_C2C<<>>(imagesIn->devData, imagesIn->height, imagesIn->width, + cuArraysCopyExtractVaryingOffset_C2C<<>>(imagesIn->devData, imagesIn->height, imagesIn->width, imagesOut->devData, imagesOut->height, imagesOut->width, imagesOut->count, offsets->devData); getLastCudaError("cuArraysCopyExtractC2C error"); @@ -202,26 +229,29 @@ void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut // correlation surface extraction (Minyan Zhong) __global__ void cuArraysCopyExtractVaryingOffsetCorr(const float *imageIn, const int inNX, const int inNY, - float *imageOut, const int outNX, const int outNY, int *imageValid, const int nImages, + float *imageOut, const int outNX, const int outNY, int *imageValid, const int nImages, const int2 *maxloc) { int idxImage = blockIdx.z; - int outx = threadIdx.x + blockDim.x*blockIdx.x; + // One thread per out point. Find the coordinates within the current image. + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; + // Find the correponding input. int inx = outx + maxloc[idxImage].x - outNX/2; int iny = outy + maxloc[idxImage].y - outNY/2; - - if (outx < outNX && outy < outNY) + + if (outx < outNX && outy < outNY) { + // Find the location in full array. int idxOut = ( blockIdx.z * outNX + outx ) * outNY + outy; int idxIn = ( blockIdx.z * inNX + inx ) * inNY + iny; if (inx>=0 && iny>=0 && inx *imagesIn, cuArrays *imagesO __global__ void cuArraysCopyExtractFixedOffset(const float *imageIn, const int inNX, const int inNY, - float *imageOut, const int outNX, const int outNY, const int nImages, + float *imageOut, const int outNX, const int outNY, const int nImages, const int offsetX, const int offsetY) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) - { + { int idxOut = (blockIdx.z * outNX + outx)*outNY+outy; - int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; + int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; imageOut[idxOut] = imageIn[idxIn]; } } -/* copy a tile of images to another image, with starting pixels offsets +/* copy a tile of images to another image, with starting pixels offsets * param[in] imageIn inut images * param[out] imageOut output images of dimension nImages*outNX*outNY */ @@ -279,23 +309,24 @@ void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, const int nthreads = 16; dim3 threadsperblock(nthreads, nthreads,1); dim3 blockspergrid(IDIVUP(imagesOut->height,nthreads), IDIVUP(imagesOut->width,nthreads), imagesOut->count); - cuArraysCopyExtractFixedOffset<<>>(imagesIn->devData, imagesIn->height, imagesIn->width, + cuArraysCopyExtractFixedOffset<<>>(imagesIn->devData, imagesIn->height, imagesIn->width, imagesOut->devData, imagesOut->height, imagesOut->width, imagesOut->count, offset.x, offset.y); getLastCudaError("cuArraysCopyExtract error"); } +// __global__ void cuArraysCopyExtract_C2C_FixedOffset(const float2 *imageIn, const int inNX, const int inNY, - float2 *imageOut, const int outNX, const int outNY, const int nImages, + float2 *imageOut, const int outNX, const int outNY, const int nImages, const int offsetX, const int offsetY) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) - { + { int idxOut = (blockIdx.z * outNX + outx)*outNY+outy; - int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; + int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; imageOut[idxOut] = imageIn[idxIn]; } } @@ -311,27 +342,64 @@ void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut //imagesIn->debuginfo(stream); //imagesOut->debuginfo(stream); cuArraysCopyExtract_C2C_FixedOffset<<>> - (imagesIn->devData, imagesIn->height, imagesIn->width, + (imagesIn->devData, imagesIn->height, imagesIn->width, imagesOut->devData, imagesOut->height, imagesOut->width, imagesOut->count, offset.x, offset.y); getLastCudaError("cuArraysCopyExtractC2C error"); } +// -__global__ void cuArraysCopyExtract_C2R_FixedOffset(const float2 *imageIn, const int inNX, const int inNY, - float *imageOut, const int outNX, const int outNY, const int nImages, +// float3 +__global__ void cuArraysCopyExtract_C2C_FixedOffset(const float3 *imageIn, const int inNX, const int inNY, + float3 *imageOut, const int outNX, const int outNY, const int nImages, const int offsetX, const int offsetY) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) - { + { int idxOut = (blockIdx.z * outNX + outx)*outNY+outy; - int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; + int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; + imageOut[idxOut] = imageIn[idxIn]; + } +} + + +void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, int2 offset, cudaStream_t stream) +{ + //assert(imagesIn->height >= imagesOut && inNY >= outNY); + const int nthreads = NTHREADS2D; + dim3 threadsperblock(nthreads, nthreads,1); + dim3 blockspergrid(IDIVUP(imagesOut->height,nthreads), IDIVUP(imagesOut->width,nthreads), imagesOut->count); + //std::cout << "debug copyExtract" << imagesOut->width << imagesOut->height << "\n"; + //imagesIn->debuginfo(stream); + //imagesOut->debuginfo(stream); + cuArraysCopyExtract_C2C_FixedOffset<<>> + (imagesIn->devData, imagesIn->height, imagesIn->width, + imagesOut->devData, imagesOut->height, imagesOut->width, imagesOut->count, offset.x, offset.y); + getLastCudaError("cuArraysCopyExtractFloat3 error"); +} + +// + + +__global__ void cuArraysCopyExtract_C2R_FixedOffset(const float2 *imageIn, const int inNX, const int inNY, + float *imageOut, const int outNX, const int outNY, const int nImages, + const int offsetX, const int offsetY) +{ + int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outy = threadIdx.y + blockDim.y*blockIdx.y; + + if(outx < outNX && outy < outNY) + { + int idxOut = (blockIdx.z * outNX + outx)*outNY+outy; + int idxIn = (blockIdx.z*inNX + outx + offsetX)*inNY + outy + offsetY; imageOut[idxOut] = imageIn[idxIn].x; } } + void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, int2 offset, cudaStream_t stream) { //assert(imagesIn->height >= imagesOut && inNY >= outNY); @@ -339,16 +407,16 @@ void cuArraysCopyExtract(cuArrays *imagesIn, cuArrays *imagesOut, dim3 threadsperblock(nthreads, nthreads,1); dim3 blockspergrid(IDIVUP(imagesOut->height,nthreads), IDIVUP(imagesOut->width,nthreads), imagesOut->count); cuArraysCopyExtract_C2R_FixedOffset<<>> - (imagesIn->devData, imagesIn->height, imagesIn->width, + (imagesIn->devData, imagesIn->height, imagesIn->width, imagesOut->devData, imagesOut->height, imagesOut->width, imagesOut->count, offset.x, offset.y); getLastCudaError("cuArraysCopyExtractC2C error"); } - +// __global__ void cuArraysCopyInsert_kernel(const float2* imageIn, const int inNX, const int inNY, float2* imageOut, const int outNY, const int offsetX, const int offsetY) { - int inx = threadIdx.x + blockDim.x*blockIdx.x; + int inx = threadIdx.x + blockDim.x*blockIdx.x; int iny = threadIdx.y + blockDim.y*blockIdx.y; if(inx < inNX && iny < inNY) { int idxOut = IDX2R(inx+offsetX, iny+offsetY, outNY); @@ -363,16 +431,40 @@ void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, i const int nthreads = 16; dim3 threadsperblock(nthreads, nthreads); dim3 blockspergrid(IDIVUP(imageIn->height,nthreads), IDIVUP(imageIn->width,nthreads)); - cuArraysCopyInsert_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, + cuArraysCopyInsert_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, + imageOut->devData, imageOut->width, offsetX, offsetY); + getLastCudaError("cuArraysCopyInsert error"); +} +// +// float3 +__global__ void cuArraysCopyInsert_kernel(const float3* imageIn, const int inNX, const int inNY, + float3* imageOut, const int outNY, const int offsetX, const int offsetY) +{ + int inx = threadIdx.x + blockDim.x*blockIdx.x; + int iny = threadIdx.y + blockDim.y*blockIdx.y; + if(inx < inNX && iny < inNY) { + int idxOut = IDX2R(inx+offsetX, iny+offsetY, outNY); + int idxIn = IDX2R(inx, iny, inNY); + imageOut[idxOut] = make_float3(imageIn[idxIn].x, imageIn[idxIn].y, imageIn[idxIn].z); + } +} + +void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int offsetX, int offsetY, cudaStream_t stream) +{ + const int nthreads = 16; + dim3 threadsperblock(nthreads, nthreads); + dim3 blockspergrid(IDIVUP(imageIn->height,nthreads), IDIVUP(imageIn->width,nthreads)); + cuArraysCopyInsert_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, imageOut->devData, imageOut->width, offsetX, offsetY); getLastCudaError("cuArraysCopyInsert error"); } +// __global__ void cuArraysCopyInsert_kernel(const float* imageIn, const int inNX, const int inNY, float* imageOut, const int outNY, const int offsetX, const int offsetY) { - int inx = threadIdx.x + blockDim.x*blockIdx.x; + int inx = threadIdx.x + blockDim.x*blockIdx.x; int iny = threadIdx.y + blockDim.y*blockIdx.y; if(inx < inNX && iny < inNY) { int idxOut = IDX2R(inx+offsetX, iny+offsetY, outNY); @@ -387,18 +479,44 @@ void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int const int nthreads = 16; dim3 threadsperblock(nthreads, nthreads); dim3 blockspergrid(IDIVUP(imageIn->height,nthreads), IDIVUP(imageIn->width,nthreads)); - cuArraysCopyInsert_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, + cuArraysCopyInsert_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, imageOut->devData, imageOut->width, offsetX, offsetY); getLastCudaError("cuArraysCopyInsert Float error"); } +// + +__global__ void cuArraysCopyInsert_kernel(const int* imageIn, const int inNX, const int inNY, + int* imageOut, const int outNY, const int offsetX, const int offsetY) +{ + int inx = threadIdx.x + blockDim.x*blockIdx.x; + int iny = threadIdx.y + blockDim.y*blockIdx.y; + if(inx < inNX && iny < inNY) { + int idxOut = IDX2R(inx+offsetX, iny+offsetY, outNY); + int idxIn = IDX2R(inx, iny, inNY); + imageOut[idxOut] = imageIn[idxIn]; + } +} + + +void cuArraysCopyInsert(cuArrays *imageIn, cuArrays *imageOut, int offsetX, int offsetY, cudaStream_t stream) +{ + const int nthreads = 16; + dim3 threadsperblock(nthreads, nthreads); + dim3 blockspergrid(IDIVUP(imageIn->height,nthreads), IDIVUP(imageIn->width,nthreads)); + cuArraysCopyInsert_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, + imageOut->devData, imageOut->width, offsetX, offsetY); + getLastCudaError("cuArraysCopyInsert Integer error"); +} +// + __global__ void cuArraysCopyInversePadded_kernel(float *imageIn, int inNX, int inNY, int sizeIn, float *imageOut, int outNX, int outNY, int sizeOut, int nImages) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) { int idxImage = blockIdx.z; @@ -409,27 +527,27 @@ __global__ void cuArraysCopyInversePadded_kernel(float *imageIn, int inNX, int i } else { imageOut[idxOut] = 0.0f; } - } + } } void cuArraysCopyInversePadded(cuArrays *imageIn, cuArrays *imageOut,cudaStream_t stream) { const int nthreads = 16; - int nImages = imageIn->count; + int nImages = imageIn->count; dim3 blockSize(nthreads, nthreads,1); dim3 gridSize(IDIVUP(imageOut->height,nthreads), IDIVUP(imageOut->width,nthreads), nImages); cuArraysCopyInversePadded_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, imageIn->size, imageOut->devData, imageOut->height, imageOut->width, imageOut->size, nImages); - getLastCudaError("cuArraysCopyInversePadded error"); + getLastCudaError("cuArraysCopyInversePadded error"); } __global__ void cuArraysCopyPadded_R2R_kernel(float *imageIn, int inNX, int inNY, int sizeIn, float *imageOut, int outNX, int outNY, int sizeOut, int nImages) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) { int idxImage = blockIdx.z; @@ -440,26 +558,26 @@ __global__ void cuArraysCopyPadded_R2R_kernel(float *imageIn, int inNX, int inNY } else { imageOut[idxOut] = 0.0f; } - } + } } void cuArraysCopyPadded(cuArrays *imageIn, cuArrays *imageOut,cudaStream_t stream) { const int nthreads = 16; - int nImages = imageIn->count; + int nImages = imageIn->count; dim3 blockSize(nthreads, nthreads,1); dim3 gridSize(IDIVUP(imageOut->height,nthreads), IDIVUP(imageOut->width,nthreads), nImages); cuArraysCopyPadded_R2R_kernel<<>>(imageIn->devData, imageIn->height, imageIn->width, imageIn->size, imageOut->devData, imageOut->height, imageOut->width, imageOut->size, nImages); - getLastCudaError("cuArraysCopyPaddedR2R error"); + getLastCudaError("cuArraysCopyPaddedR2R error"); } __global__ void cuArraysCopyPadded_C2C_kernel(float2 *imageIn, int inNX, int inNY, int sizeIn, float2 *imageOut, int outNX, int outNY, int sizeOut, int nImages) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) { int idxImage = blockIdx.z; @@ -468,31 +586,31 @@ __global__ void cuArraysCopyPadded_C2C_kernel(float2 *imageIn, int inNX, int inN int idxIn = IDX2R(outx, outy, inNY)+idxImage*sizeIn; imageOut[idxOut] = imageIn[idxIn]; } - else{ - imageOut[idxOut] = make_float2(0.0f, 0.0f); + else{ + imageOut[idxOut] = make_float2(0.0f, 0.0f); } - } + } } void cuArraysCopyPadded(cuArrays *imageIn, cuArrays *imageOut,cudaStream_t stream) { const int nthreads = NTHREADS2D; - int nImages = imageIn->count; + int nImages = imageIn->count; dim3 blockSize(nthreads, nthreads,1); dim3 gridSize(IDIVUP(imageOut->height,nthreads), IDIVUP(imageOut->width,nthreads), nImages); cuArraysCopyPadded_C2C_kernel<<>> (imageIn->devData, imageIn->height, imageIn->width, imageIn->size, imageOut->devData, imageOut->height, imageOut->width, imageOut->size, nImages); - getLastCudaError("cuArraysCopyInversePadded error"); + getLastCudaError("cuArraysCopyInversePadded error"); } __global__ void cuArraysCopyPadded_R2C_kernel(float *imageIn, int inNX, int inNY, int sizeIn, float2 *imageOut, int outNX, int outNY, int sizeOut, int nImages) { - int outx = threadIdx.x + blockDim.x*blockIdx.x; + int outx = threadIdx.x + blockDim.x*blockIdx.x; int outy = threadIdx.y + blockDim.y*blockIdx.y; - + if(outx < outNX && outy < outNY) { int idxImage = blockIdx.z; @@ -501,42 +619,42 @@ __global__ void cuArraysCopyPadded_R2C_kernel(float *imageIn, int inNX, int inNY int idxIn = IDX2R(outx, outy, inNY)+idxImage*sizeIn; imageOut[idxOut] = make_float2(imageIn[idxIn], 0.0f); } - else{ - imageOut[idxOut] = make_float2(0.0f, 0.0f); + else{ + imageOut[idxOut] = make_float2(0.0f, 0.0f); } - } + } } void cuArraysCopyPadded(cuArrays *imageIn, cuArrays *imageOut,cudaStream_t stream) { const int nthreads = NTHREADS2D; - int nImages = imageIn->count; + int nImages = imageIn->count; dim3 blockSize(nthreads, nthreads,1); dim3 gridSize(IDIVUP(imageOut->height,nthreads), IDIVUP(imageOut->width,nthreads), nImages); cuArraysCopyPadded_R2C_kernel<<>> (imageIn->devData, imageIn->height, imageIn->width, imageIn->size, imageOut->devData, imageOut->height, imageOut->width, imageOut->size, nImages); - getLastCudaError("cuArraysCopyPadded error"); + getLastCudaError("cuArraysCopyPadded error"); } __global__ void cuArraysSetConstant_kernel(float *image, int size, float value) { - int idx = threadIdx.x + blockDim.x*blockIdx.x; - + int idx = threadIdx.x + blockDim.x*blockIdx.x; + if(idx < size) { image[idx] = value; - } + } } void cuArraysSetConstant(cuArrays *imageIn, float value, cudaStream_t stream) { const int nthreads = 256; - int size = imageIn->getSize(); + int size = imageIn->getSize(); cuArraysSetConstant_kernel<<>> (imageIn->devData, imageIn->size, value); - getLastCudaError("cuArraysCopyPadded error"); + getLastCudaError("cuArraysCopyPadded error"); } diff --git a/contrib/PyCuAmpcor/src/cuCorrNormalization.cu b/contrib/PyCuAmpcor/src/cuCorrNormalization.cu index 83dc1cf..7c9c339 100644 --- a/contrib/PyCuAmpcor/src/cuCorrNormalization.cu +++ b/contrib/PyCuAmpcor/src/cuCorrNormalization.cu @@ -195,7 +195,6 @@ __device__ float2 partialSums(const float v, volatile float* shmem, const int st return make_float2(Sum, Sum2); } -__forceinline__ __device__ int __mul(const int a, const int b) { return a*b; } template __global__ void cuCorrNormalize_kernel( @@ -232,7 +231,7 @@ __global__ void cuCorrNormalize_kernel( templateSum += templateD[i]; } templateSum = sumReduceBlock(templateSum, shmem); - + __syncthreads(); float templateSum2 = 0.0f; for (int i = tid; i < templateSize; i += Nthreads) @@ -241,11 +240,12 @@ __global__ void cuCorrNormalize_kernel( templateSum2 += t*t; } templateSum2 = sumReduceBlock(templateSum2, shmem); + __syncthreads(); //if(tid ==0) printf("template sum %d %g %g \n", imageIdx, templateSum, templateSum2); /*********/ - shmem[tid] = shmem[tid + Nthreads] = 0.0f; + shmem[tid] = shmem[tid + Nthreads] = shmem[tid + 2*Nthreads] = 0.0f; __syncthreads(); float imageSum = 0.0f; @@ -281,7 +281,7 @@ __global__ void cuCorrNormalize_kernel( if (tid < resultNY) { const int ix = iaddr/imageNY; - const int addr = __mul(ix-templateNX, resultNY); + const int addr = (ix-templateNX)*resultNY; //printf("test norm %d %d %d %d %f\n", tid, ix, addr, addr+tid, resultD[addr + tid]); diff --git a/contrib/PyCuAmpcor/src/cuEstimateStats.cu b/contrib/PyCuAmpcor/src/cuEstimateStats.cu index 8fa7e45..49c1d71 100644 --- a/contrib/PyCuAmpcor/src/cuEstimateStats.cu +++ b/contrib/PyCuAmpcor/src/cuEstimateStats.cu @@ -25,7 +25,7 @@ __global__ void cudaKernel_estimateSnr(const float* corrSum, const int* corrVali float mean = (corrSum[idx] - maxval[idx] * maxval[idx]) / (corrValidCount[idx] - 1); - snrValue[idx] = maxval[idx] / mean; + snrValue[idx] = maxval[idx] * maxval[idx] / mean; } void cuEstimateSnr(cuArrays *corrSum, cuArrays *corrValidCount, cuArrays *maxval, cuArrays *snrValue, cudaStream_t stream) @@ -55,7 +55,7 @@ void cuEstimateSnr(cuArrays *corrSum, cuArrays *corrValidCount, cuAr //for (int i=0; ihostData[i]<hostData[i]< *corrSum, cuArrays *corrValidCount, cuAr getLastCudaError("cuda kernel estimate stats error\n"); } + + +template // number of threads per block. +__global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX, const int NY, const int2* maxloc, const float* maxval, float3* covValue, const int size) +{ + + // Find image id. + int idxImage = threadIdx.x + blockDim.x*blockIdx.x; + + if (idxImage >= size) return; + + // Preparation. + int px = maxloc[idxImage].x; + int py = maxloc[idxImage].y; + float peak = maxval[idxImage]; + + // Check if maxval is on the margin. + if (px-1 < 0 || py-1 <0 || px + 1 >=NX || py+1 >=NY) { + + covValue[idxImage] = make_float3(99.0, 99.0, 99.0); + + } + else { + int offset = NX * NY * idxImage; + int idx00 = offset + (px - 1) * NY + py - 1; + int idx01 = offset + (px - 1) * NY + py ; + int idx02 = offset + (px - 1) * NY + py + 1; + int idx10 = offset + (px ) * NY + py - 1; + int idx11 = offset + (px ) * NY + py ; + int idx12 = offset + (px ) * NY + py + 1; + int idx20 = offset + (px + 1) * NY + py - 1; + int idx21 = offset + (px + 1) * NY + py ; + int idx22 = offset + (px + 1) * NY + py + 1; + + float dxx = - ( corrBatchRaw[idx21] + corrBatchRaw[idx01] - 2*corrBatchRaw[idx11] ) * 0.5; + float dyy = - ( corrBatchRaw[idx12] + corrBatchRaw[idx10] - 2*corrBatchRaw[idx11] ) * 0.5; + float dxy = - ( corrBatchRaw[idx22] + corrBatchRaw[idx00] - corrBatchRaw[idx20] - corrBatchRaw[idx02] ) *0.25; + + float n2 = fmaxf(1 - peak, 0.0); + + int winSize = NX*NY; + + dxx = dxx * winSize; + dyy = dyy * winSize; + dxy = dxy * winSize; + + float n4 = n2*n2; + n2 = n2 * 2; + n4 = n4 * 0.5 * winSize; + + float u = dxy * dxy - dxx * dyy; + float u2 = u*u; + + if (fabsf(u) < 1e-2) { + + covValue[idxImage] = make_float3(99.0, 99.0, 99.0); + + } + else { + float cov_xx = (- n2 * u * dyy + n4 * ( dyy*dyy + dxy*dxy) ) / u2; + float cov_yy = (- n2 * u * dxx + n4 * ( dxx*dxx + dxy*dxy) ) / u2; + float cov_xy = ( n2 * u * dxy - n4 * ( dxx + dyy ) * dxy ) / u2; + covValue[idxImage] = make_float3(cov_xx, cov_yy, cov_xy); + } + } +} + +void cuEstimateVariance(cuArrays *corrBatchRaw, cuArrays *maxloc, cuArrays *maxval, cuArrays *covValue, cudaStream_t stream) +{ + + int size = corrBatchRaw->count; + + // One dimensional launching parameters to loop over every correlation surface. + cudaKernel_estimateVar<<< IDIVUP(size, NTHREADS), NTHREADS, 0, stream>>> + (corrBatchRaw->devData, corrBatchRaw->height, corrBatchRaw->width, maxloc->devData, maxval->devData, covValue->devData, size); + getLastCudaError("cudaKernel_estimateVar error\n"); +} diff --git a/contrib/PyCuAmpcor/src/setup.py b/contrib/PyCuAmpcor/src/setup.py index 6384c91..aa34441 100644 --- a/contrib/PyCuAmpcor/src/setup.py +++ b/contrib/PyCuAmpcor/src/setup.py @@ -7,20 +7,21 @@ from distutils.core import setup from distutils.extension import Extension from Cython.Build import cythonize -import os -os.environ["CC"] = "g++" +import numpy setup( name = 'PyCuAmpcor', ext_modules = cythonize(Extension( "PyCuAmpcor", sources=['PyCuAmpcor.pyx'], - include_dirs=['/usr/local/cuda/include'], # REPLACE WITH YOUR PATH TO YOUR CUDA LIBRARY HEADERS + include_dirs=['/usr/local/cuda/include', numpy.get_include()], # REPLACE WITH YOUR PATH TO YOUR CUDA LIBRARY HEADERS extra_compile_args=['-fPIC','-fpermissive'], - extra_objects=['SlcImage.o','cuAmpcorChunk.o','cuAmpcorParameter.o','cuCorrFrequency.o', + extra_objects=['GDALImage.o','cuAmpcorChunk.o','cuAmpcorParameter.o','cuCorrFrequency.o', 'cuCorrNormalization.o','cuCorrTimeDomain.o','cuArraysCopy.o', 'cuArrays.o','cuArraysPadding.o','cuOffset.o','cuOverSampler.o', - 'cuSincOverSampler.o', 'cuDeramp.o','cuAmpcorController.o'], - extra_link_args=['-L/usr/local/cuda/lib64','-lcuda','-lcudart','-lcufft','-lcublas'], # REPLACE FIRST PATH WITH YOUR PATH TO YOUR CUDA LIBRARIES + 'cuSincOverSampler.o', 'cuDeramp.o','cuAmpcorController.o','cuEstimateStats.o'], + extra_link_args=['-L/usr/local/cuda/lib64', + '-L/usr/lib64/nvidia', + '-lcuda','-lcudart','-lcufft','-lcublas','-lgdal'], # REPLACE FIRST PATH WITH YOUR PATH TO YOUR CUDA LIBRARIES language='c++' ))) From 6e36a73fb0a0e34cc2e30bb06ca5995bd14deb3a Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Mon, 18 Mar 2019 19:00:18 -0400 Subject: [PATCH 07/27] stripmapStack/unwrap: enable corrlooks for SNAPHU Stack: + add alks and rlks while writing config file for unwrap unwrap: + calculate the az/rg spacing/resolution and the equivalent number of independent looks in extractInfoFromPickle() + uncomment corrlooks setting in runUnwrap() --- contrib/stack/stripmapStack/Stack.py | 2 ++ contrib/stack/stripmapStack/unwrap.py | 16 +++++++++++++--- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/contrib/stack/stripmapStack/Stack.py b/contrib/stack/stripmapStack/Stack.py index 74d0df4..b0a116e 100755 --- a/contrib/stack/stripmapStack/Stack.py +++ b/contrib/stack/stripmapStack/Stack.py @@ -197,6 +197,8 @@ class config(object): self.f.write('nomcf : ' + self.noMCF + '\n') self.f.write('master : ' + self.master + '\n') self.f.write('defomax : ' + self.defoMax + '\n') + self.f.write('alks : ' + self.alks + '\n') + self.f.write('rlks : ' + self.rlks + '\n') self.f.write('method : ' + self.unwMethod + '\n') self.f.write('##########################'+'\n') diff --git a/contrib/stack/stripmapStack/unwrap.py b/contrib/stack/stripmapStack/unwrap.py index bae4893..f1de7f9 100755 --- a/contrib/stack/stripmapStack/unwrap.py +++ b/contrib/stack/stripmapStack/unwrap.py @@ -113,9 +113,19 @@ def extractInfoFromPickle(pckfile, inps): data['earthRadius'] = elp.local_radius_of_curvature(llh.lat, hdg) #azspacing = burst.azimuthTimeInterval * sv.getScalarVelocity() - azres = 20.0 + #azres = 20.0 + azspacing = sv.getScalarVelocity() / burst.PRF + azres = burst.platform.antennaLength / 2.0 + azfact = azres / azspacing + + burst.getInstrument() + rgBandwidth = burst.instrument.pulseLength * burst.instrument.chirpSlope + rgres = abs(SPEED_OF_LIGHT / (2.0 * rgBandwidth)) + rgspacing = burst.instrument.rangePixelSize + rgfact = rgres / rgspacing #data['corrlooks'] = inps.rglooks * inps.azlooks * azspacing / azres + data['corrlooks'] = inps.rglooks * inps.azlooks / (azfact * rgfact) data['rglooks'] = inps.rglooks data['azlooks'] = inps.azlooks @@ -149,7 +159,7 @@ def runUnwrap(infile, outfile, corfile, config, costMode = None,initMethod = Non altitude = config['altitude'] rangeLooks = config['rglooks'] azimuthLooks = config['azlooks'] - #corrLooks = config['corrlooks'] + corrLooks = config['corrlooks'] maxComponents = 20 snp = Snaphu() @@ -163,7 +173,7 @@ def runUnwrap(infile, outfile, corfile, config, costMode = None,initMethod = Non snp.setAltitude(altitude) snp.setCorrfile(corfile) snp.setInitMethod(initMethod) - # snp.setCorrLooks(corrLooks) + snp.setCorrLooks(corrLooks) snp.setMaxComponents(maxComponents) snp.setDefoMaxCycles(defomax) snp.setRangeLooks(rangeLooks) From 6e4fcadce863e37a997c395c21aa5b9631ba1bdd Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Thu, 21 Nov 2019 13:18:25 -0800 Subject: [PATCH 08/27] fix bug while calculating corrLooks for snaphu(_mcf) Currently, the number of looks in azimuth direction is used twice instead of once while calculating the `corrLooks` value. Below is the relevant document of `NCORRLOOKS` in SNPAHU: https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/snaphu.conf.full --- components/isceobj/Unwrap/snaphu.py | 2 +- components/isceobj/Unwrap/snaphu_mcf.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/components/isceobj/Unwrap/snaphu.py b/components/isceobj/Unwrap/snaphu.py index 88073a6..cbb51b4 100755 --- a/components/isceobj/Unwrap/snaphu.py +++ b/components/isceobj/Unwrap/snaphu.py @@ -54,7 +54,7 @@ class snaphu(Component): self.azimuthLooks = obj.insar.topo.numberAzimuthLooks azres = obj.insar.masterFrame.platform.antennaLength/2.0 - azfact = obj.insar.topo.numberAzimuthLooks *azres / obj.insar.topo.azimuthSpacing + azfact = azres / obj.insar.topo.azimuthSpacing rBW = obj.insar.masterFrame.instrument.pulseLength * obj.insar.masterFrame.instrument.chirpSlope rgres = abs(SPEED_OF_LIGHT / (2.0 * rBW)) diff --git a/components/isceobj/Unwrap/snaphu_mcf.py b/components/isceobj/Unwrap/snaphu_mcf.py index 0705d8a..f6cd5fd 100755 --- a/components/isceobj/Unwrap/snaphu_mcf.py +++ b/components/isceobj/Unwrap/snaphu_mcf.py @@ -54,7 +54,7 @@ class snaphu_mcf(Component): self.azimuthLooks = obj.insar.topo.numberAzimuthLooks azres = obj.insar.masterFrame.platform.antennaLength/2.0 - azfact = obj.insar.topo.numberAzimuthLooks *azres / obj.insar.topo.azimuthSpacing + azfact = azres / obj.insar.topo.azimuthSpacing rBW = obj.insar.masterFrame.instrument.pulseLength * obj.insar.masterFrame.instrument.chirpSlope rgres = abs(SPEED_OF_LIGHT / (2.0 * rBW)) From 43aee1dcb327df8c6e2d2f4d8d53f8aa82e78244 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Thu, 21 Nov 2019 13:29:34 -0800 Subject: [PATCH 09/27] stripmapStack/unwrap: adjust indentation --- contrib/stack/stripmapStack/unwrap.py | 90 ++++++++++++++------------- 1 file changed, 47 insertions(+), 43 deletions(-) diff --git a/contrib/stack/stripmapStack/unwrap.py b/contrib/stack/stripmapStack/unwrap.py index f1de7f9..77080ce 100755 --- a/contrib/stack/stripmapStack/unwrap.py +++ b/contrib/stack/stripmapStack/unwrap.py @@ -258,33 +258,34 @@ def runUnwrapIcu(infile, outfile): unwImage.finalizeImage() unwImage.renderHdr() -def runUnwrap2Stage(unwrappedIntFilename,connectedComponentsFilename,unwrapped2StageFilename, unwrapper_2stage_name=None, solver_2stage=None): +def runUnwrap2Stage(unwrappedIntFilename,connectedComponentsFilename,unwrapped2StageFilename, + unwrapper_2stage_name=None, solver_2stage=None): - if unwrapper_2stage_name is None: - unwrapper_2stage_name = 'REDARC0' - - if solver_2stage is None: - # If unwrapper_2state_name is MCF then solver is ignored - # and relaxIV MCF solver is used by default - solver_2stage = 'pulp' - - print('Unwrap 2 Stage Settings:') - print('Name: %s'%unwrapper_2stage_name) - print('Solver: %s'%solver_2stage) - - inpFile = unwrappedIntFilename - ccFile = connectedComponentsFilename - outFile = unwrapped2StageFilename - - # Hand over to 2Stage unwrap - unw = UnwrapComponents() - unw.setInpFile(inpFile) - unw.setConnCompFile(ccFile) - unw.setOutFile(outFile) - unw.setSolver(solver_2stage) - unw.setRedArcs(unwrapper_2stage_name) - unw.unwrapComponents() - return + if unwrapper_2stage_name is None: + unwrapper_2stage_name = 'REDARC0' + + if solver_2stage is None: + # If unwrapper_2state_name is MCF then solver is ignored + # and relaxIV MCF solver is used by default + solver_2stage = 'pulp' + + print('Unwrap 2 Stage Settings:') + print('Name: %s'%unwrapper_2stage_name) + print('Solver: %s'%solver_2stage) + + inpFile = unwrappedIntFilename + ccFile = connectedComponentsFilename + outFile = unwrapped2StageFilename + + # Hand over to 2Stage unwrap + unw = UnwrapComponents() + unw.setInpFile(inpFile) + unw.setConnCompFile(ccFile) + unw.setOutFile(outFile) + unw.setSolver(solver_2stage) + unw.setRedArcs(unwrapper_2stage_name) + unw.unwrapComponents() + return def main(iargs=None): @@ -303,24 +304,26 @@ def main(iargs=None): if inps.method != 'icu': - masterShelveDir = os.path.join(interferogramDir , 'masterShelve') - if not os.path.exists(masterShelveDir): - os.makedirs(masterShelveDir) + masterShelveDir = os.path.join(interferogramDir , 'masterShelve') + if not os.path.exists(masterShelveDir): + os.makedirs(masterShelveDir) + + inps.master = os.path.dirname(inps.master) + cpCmd='cp ' + os.path.join(inps.master, 'data*') +' '+masterShelveDir + os.system(cpCmd) + pckfile = os.path.join(masterShelveDir,'data') + print(pckfile) + metadata = extractInfoFromPickle(pckfile, inps) - inps.master = os.path.dirname(inps.master) - cpCmd='cp ' + os.path.join(inps.master, 'data*') +' '+masterShelveDir - os.system(cpCmd) - pckfile = os.path.join(masterShelveDir,'data') - print(pckfile) - metadata = extractInfoFromPickle(pckfile, inps) ######## print ('unwrapping method : ' , inps.method) if inps.method == 'snaphu': - if inps.nomcf: - fncall = runUnwrap - else: - fncall = runUnwrapMcf - fncall(inps.intfile, inps.unwprefix + '_snaphu.unw', inps.cohfile, metadata, defomax=inps.defomax) + if inps.nomcf: + fncall = runUnwrap + else: + fncall = runUnwrapMcf + fncall(inps.intfile, inps.unwprefix + '_snaphu.unw', inps.cohfile, metadata, defomax=inps.defomax) + elif inps.method == 'snaphu2stage': if inps.nomcf: fncall = runUnwrap @@ -329,11 +332,12 @@ def main(iargs=None): fncall(inps.intfile, inps.unwprefix + '_snaphu.unw', inps.cohfile, metadata, defomax=inps.defomax) # adding in the two-stage - runUnwrap2Stage(inps.unwprefix + '_snaphu.unw', inps.unwprefix + '_snaphu.unw.conncomp',inps.unwprefix + '_snaphu2stage.unw') - + runUnwrap2Stage(inps.unwprefix + '_snaphu.unw', + inps.unwprefix + '_snaphu.unw.conncomp', + inps.unwprefix + '_snaphu2stage.unw') elif inps.method == 'icu': - runUnwrapIcu(inps.intfile, inps.unwprefix + '_icu.unw') + runUnwrapIcu(inps.intfile, inps.unwprefix + '_icu.unw') if __name__ == '__main__': From 21675df72476af2dfa81d5aba557ab07527ce535 Mon Sep 17 00:00:00 2001 From: piyushrpt Date: Tue, 10 Dec 2019 09:54:36 -0800 Subject: [PATCH 10/27] Bug fixes for rubbersheeting --- components/isceobj/StripmapProc/runResampleSlc.py | 1 + examples/input_files/stripmapApp.xml | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/components/isceobj/StripmapProc/runResampleSlc.py b/components/isceobj/StripmapProc/runResampleSlc.py index c108c98..b1dfa72 100644 --- a/components/isceobj/StripmapProc/runResampleSlc.py +++ b/components/isceobj/StripmapProc/runResampleSlc.py @@ -75,6 +75,7 @@ def runResampleSlc(self, kind='coarse'): if kind in ['coarse', 'refined']: azname = os.path.join(offsetsDir, self.insar.azimuthOffsetFilename) rgname = os.path.join(offsetsDir, self.insar.rangeOffsetFilename) + flatten = True else: azname = os.path.join(offsetsDir, self.insar.azimuthRubbersheetFilename) if self.doRubbersheetingRange: diff --git a/examples/input_files/stripmapApp.xml b/examples/input_files/stripmapApp.xml index f4eff23..87c1365 100644 --- a/examples/input_files/stripmapApp.xml +++ b/examples/input_files/stripmapApp.xml @@ -12,7 +12,8 @@ /Users/fattahi/process/test_roiApp/Alos_Maule_T116/demLat_S39_S35_Lon_W074_W071.dem.wgs84 True True From 80d25ca7279b447cf27d49593e293eecc165f61c Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Sat, 14 Dec 2019 15:30:51 -0800 Subject: [PATCH 11/27] Simplify python version info --- SConstruct | 35 ++--------------------------------- __init__.py | 10 +++------- 2 files changed, 5 insertions(+), 40 deletions(-) diff --git a/SConstruct b/SConstruct index 8e653fe..a979b3c 100644 --- a/SConstruct +++ b/SConstruct @@ -216,43 +216,12 @@ else: ### End of GPU branch-specific modifications -file = '__init__.py' -if not os.path.exists(file): - fout = open(file,"w") - fout.write("#!/usr/bin/env python3") - fout.close() - -env.Install(inst,file) -try: - from subprocess import check_output - svn_revision = check_output('svnversion').strip() or 'Unknown' - if sys.version_info[0] == 3: - svn_revision = svn_revision.decode('utf-8') -except ImportError: - try: - import popen2 - stdout, stdin, stderr = popen2.popen3('svnversion') - svn_revision = stdout.read().strip() - if stderr.read(): - raise Exception - except Exception: - svn_revision = 'Unknown' -except OSError: - svn_revision = 'Unknown' +env.Install(inst, '__init__.py') +env.Install(inst, 'release_history.py') if not os.path.exists(inst): os.makedirs(inst) -fvers = open(os.path.join(inst,'version.py'),'w') - -from release_history import release_version, release_svn_revision, release_date -fvers_lines = ["release_version = '"+release_version+"'\n", - "release_svn_revision = '"+release_svn_revision+"'\n", - "release_date = '"+release_date+"'\n", - "svn_revision = '"+svn_revision+"'\n\n"] - -fvers.write(''.join(fvers_lines)) -fvers.close() v = 0 if isrerun == 'no': cmd = 'scons -Q install --isrerun=yes' diff --git a/__init__.py b/__init__.py index 138b442..a6e2c5b 100755 --- a/__init__.py +++ b/__init__.py @@ -25,13 +25,9 @@ # Author: Giangi Sacco #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - - -from __future__ import print_function - -from .version import release_version, release_svn_revision, release_date -from .version import svn_revision +from .release_history import release_version, release_svn_revision, release_date +svn_revision = release_svn_revision +version = release_history # compatibility alias __version__ = release_version From 537bae03d958fd7b1aee216bd7700f6f07f62da6 Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Sat, 14 Dec 2019 12:20:41 -0800 Subject: [PATCH 12/27] Hide rubbersheeting scipy/astropy imports Recent rubbersheeting changes added scipy and astropy imports which are not otherwise needed. Due to isce2's eager import behavior this causes errors when using unrelated functionality. We can fix this by tucking these imports away inside their function definitions. To reproduce: `python3 -c "import isce; from stripmapApp import Insar"` (without scipy/astropy installed) --- .../isceobj/StripmapProc/runRubbersheet.py | 6 +++++- .../StripmapProc/runRubbersheetAzimuth.py | 15 ++++++++++---- .../StripmapProc/runRubbersheetRange.py | 20 ++++++++++++++----- 3 files changed, 31 insertions(+), 10 deletions(-) diff --git a/components/isceobj/StripmapProc/runRubbersheet.py b/components/isceobj/StripmapProc/runRubbersheet.py index b16a58e..ea417bd 100644 --- a/components/isceobj/StripmapProc/runRubbersheet.py +++ b/components/isceobj/StripmapProc/runRubbersheet.py @@ -6,7 +6,6 @@ import isce import isceobj from osgeo import gdal -from scipy import ndimage import numpy as np import os @@ -24,6 +23,9 @@ def fill(data, invalid=None): Output: Return a filled array. """ + + from scipy import ndimage + if invalid is None: invalid = np.isnan(data) ind = ndimage.distance_transform_edt(invalid, @@ -35,6 +37,8 @@ def fill(data, invalid=None): def mask_filter(denseOffsetFile, snrFile, band, snrThreshold, filterSize, outName): #masking and Filtering + from scipy import ndimage + ##Read in the offset file ds = gdal.Open(denseOffsetFile + '.vrt', gdal.GA_ReadOnly) Offset = ds.GetRasterBand(1).ReadAsArray() diff --git a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py index 75583f1..43ee4a8 100755 --- a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py +++ b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py @@ -9,14 +9,14 @@ import isce import isceobj from osgeo import gdal -from scipy import ndimage -from astropy.convolution import convolve import numpy as np import os def mask_filterNoSNR(denseOffsetFile,filterSize,outName): # Masking the offsets with a data-based approach - + + from scipy import ndimage + # Open the offsets ds = gdal.Open(denseOffsetFile+'.vrt',gdal.GA_ReadOnly) off_az = ds.GetRasterBand(1).ReadAsArray() @@ -79,6 +79,9 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): def off_masking(off,filterSize,thre=2): + + from scipy import ndimage + # Define the mask to fill the offsets vram = ndimage.median_filter(off.real, filterSize) vazm = ndimage.median_filter(off.imag, filterSize) @@ -112,6 +115,8 @@ def fill(data, invalid=None): def mask_filter(denseOffsetFile, snrFile, band, snrThreshold, filterSize, outName): #masking and Filtering + from scipy import ndimage + ##Read in the offset file ds = gdal.Open(denseOffsetFile + '.vrt', gdal.GA_ReadOnly) Offset = ds.GetRasterBand(band).ReadAsArray() @@ -154,7 +159,9 @@ def mask_filter(denseOffsetFile, snrFile, band, snrThreshold, filterSize, outNam return None def fill_with_smoothed(off,filterSize): - + + from astropy.convolution import convolve + off_2filt=np.copy(off) kernel = np.ones((filterSize,filterSize),np.float32)/(filterSize*filterSize) loop = 0 diff --git a/components/isceobj/StripmapProc/runRubbersheetRange.py b/components/isceobj/StripmapProc/runRubbersheetRange.py index 2f1ab3f..04fc2a2 100755 --- a/components/isceobj/StripmapProc/runRubbersheetRange.py +++ b/components/isceobj/StripmapProc/runRubbersheetRange.py @@ -9,15 +9,14 @@ import isce import isceobj from osgeo import gdal -from scipy import ndimage import numpy as np import os -from astropy.convolution import convolve - def mask_filterNoSNR(denseOffsetFile,filterSize,outName): # Masking the offsets with a data-based approach - + + from scipy import ndimage + # Open the offsets ds = gdal.Open(denseOffsetFile+'.vrt',gdal.GA_ReadOnly) off_az = ds.GetRasterBand(1).ReadAsArray() @@ -78,6 +77,9 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): return def off_masking(off,filterSize,thre=2): + + from scipy import ndimage + vram = ndimage.median_filter(off.real, filterSize) vazm = ndimage.median_filter(off.imag, filterSize) @@ -100,6 +102,8 @@ def fill(data, invalid=None): Output: Return a filled array. """ + from scipy import ndimage + if invalid is None: invalid = np.isnan(data) ind = ndimage.distance_transform_edt(invalid, @@ -108,7 +112,9 @@ def fill(data, invalid=None): return data[tuple(ind)] def fill_with_smoothed(off,filterSize): - + + from astropy.convolution import convolve + off_2filt=np.copy(off) kernel = np.ones((filterSize,filterSize),np.float32)/(filterSize*filterSize) loop = 0 @@ -131,6 +137,8 @@ def fill_with_smoothed(off,filterSize): def mask_filter(denseOffsetFile, snrFile, band, snrThreshold, filterSize, outName): #masking and Filtering + from scipy import ndimage + ##Read in the offset file ds = gdal.Open(denseOffsetFile + '.vrt', gdal.GA_ReadOnly) Offset = ds.GetRasterBand(band).ReadAsArray() @@ -236,6 +244,8 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): def runRubbersheetRange(self): + from scipy import ndimage + if not self.doRubbersheetingRange: print('Rubber sheeting in azimuth not requested ... skipping') return From e2a81bbd6aca1b121e6612e07e21b17e37bab0f9 Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Mon, 16 Dec 2019 10:22:01 -0800 Subject: [PATCH 13/27] Hide scipy imports for stripmapapp --- components/isceobj/StripmapProc/runDispersive.py | 3 ++- components/isceobj/TopsProc/runIon.py | 12 +++++++----- components/isceobj/TopsProc/runOffsetFilter.py | 3 ++- components/isceobj/TopsProc/runPrepESD.py | 2 +- 4 files changed, 12 insertions(+), 8 deletions(-) diff --git a/components/isceobj/StripmapProc/runDispersive.py b/components/isceobj/StripmapProc/runDispersive.py index 305022b..724c189 100644 --- a/components/isceobj/StripmapProc/runDispersive.py +++ b/components/isceobj/StripmapProc/runDispersive.py @@ -9,7 +9,6 @@ from isceobj.Constants import SPEED_OF_LIGHT import numpy as np import gdal -from scipy import ndimage try: import cv2 except ImportError: @@ -296,6 +295,8 @@ def fill(data, invalid=None): Output: Return a filled array. """ + from scipy import ndimage + if invalid is None: invalid = np.isnan(data) ind = ndimage.distance_transform_edt(invalid, diff --git a/components/isceobj/TopsProc/runIon.py b/components/isceobj/TopsProc/runIon.py index 8407287..d0c4af6 100644 --- a/components/isceobj/TopsProc/runIon.py +++ b/components/isceobj/TopsProc/runIon.py @@ -9,9 +9,6 @@ import shutil import datetime import numpy as np import numpy.matlib -import scipy.signal as ss -from scipy import interpolate -from scipy.interpolate import interp1d import isceobj import logging @@ -638,6 +635,7 @@ def cal_coherence(inf, win=5, edge=0): 4: keep all samples ''' + import scipy.signal as ss if win % 2 != 1: raise Exception('window size must be odd!') @@ -1682,6 +1680,9 @@ def computeDopplerOffset(burst, firstline, lastline, firstcolumn, lastcolumn, nr output: first lines > 0, last lines < 0 ''' + from scipy import interpolate + from scipy.interpolate import interp1d + Vs = np.linalg.norm(burst.orbit.interpolateOrbit(burst.sensingMid, method='hermite').getVelocity()) Ks = 2 * Vs * burst.azimuthSteeringRate / burst.radarWavelength @@ -1830,6 +1831,7 @@ def adaptive_gaussian(ionos, wgt, size_max, size_min): size_max: maximum window size size_min: minimum window size ''' + import scipy.signal as ss length = (ionos.shape)[0] width = (ionos.shape)[1] @@ -1892,6 +1894,8 @@ def filt_gaussian(self, ionParam): currently not implemented. a less accurate method is to use ionsphere without any projection ''' + from scipy import interpolate + from scipy.interpolate import interp1d ################################################# #SET PARAMETERS HERE @@ -2659,5 +2663,3 @@ def runIon(self): #esd_noion(self, ionParam) return - - \ No newline at end of file diff --git a/components/isceobj/TopsProc/runOffsetFilter.py b/components/isceobj/TopsProc/runOffsetFilter.py index 51735d9..e9989d9 100644 --- a/components/isceobj/TopsProc/runOffsetFilter.py +++ b/components/isceobj/TopsProc/runOffsetFilter.py @@ -3,7 +3,6 @@ # Copyright 2016 # -from scipy.ndimage.filters import median_filter import numpy as np import isce import isceobj @@ -20,6 +19,8 @@ def runOffsetFilter(self): if not self.doDenseOffsets: return + from scipy.ndimage.filters import median_filter + offsetfile = os.path.join(self._insar.mergedDirname, self._insar.offsetfile) snrfile = os.path.join(self._insar.mergedDirname, self._insar.snrfile) print('\n======================================') diff --git a/components/isceobj/TopsProc/runPrepESD.py b/components/isceobj/TopsProc/runPrepESD.py index 5d65b78..0cf9c3c 100755 --- a/components/isceobj/TopsProc/runPrepESD.py +++ b/components/isceobj/TopsProc/runPrepESD.py @@ -8,7 +8,6 @@ import numpy as np import os import isceobj import logging -import scipy.signal as SS from isceobj.Util.ImageUtil import ImageLib as IML import datetime import pprint @@ -177,6 +176,7 @@ def createCoherence(intfile, win=5): ''' Compute coherence using scipy convolve 2D. ''' + import scipy.signal as SS corfile = os.path.splitext(intfile)[0] + '.cor' filt = np.ones((win,win))/ (1.0*win*win) From 6b42742b005d3635af2f2ae11d4cf3edb1de812c Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Tue, 17 Dec 2019 12:42:06 -0800 Subject: [PATCH 14/27] Remove use of malloc.h This header is deprecated and Linux-specific. The standard include for malloc should be stdlib.h. --- components/stdproc/alosreformat/ALOS_lib/src/cfft1d_fftpack.c | 1 - 1 file changed, 1 deletion(-) diff --git a/components/stdproc/alosreformat/ALOS_lib/src/cfft1d_fftpack.c b/components/stdproc/alosreformat/ALOS_lib/src/cfft1d_fftpack.c index 773be84..56ddc80 100644 --- a/components/stdproc/alosreformat/ALOS_lib/src/cfft1d_fftpack.c +++ b/components/stdproc/alosreformat/ALOS_lib/src/cfft1d_fftpack.c @@ -1,7 +1,6 @@ #include #include #include -#include /************************************************************************ * cfft1d is a subroutine used to call and initialize perflib Fortran FFT * * routines. * From 7df2b44a0cdc96016b56a4432710e7422c5a805c Mon Sep 17 00:00:00 2001 From: HBaldwin3 <49163135+HBaldwin3@users.noreply.github.com> Date: Wed, 18 Dec 2019 14:44:34 -0600 Subject: [PATCH 15/27] Change doRubbersheeting to doRubbersheetingAzimuth. --- .../isceobj/StripmapProc/runDenseOffsets.py | 10 +++++----- .../isceobj/StripmapProc/runRubbersheet.py | 16 +++++++--------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/components/isceobj/StripmapProc/runDenseOffsets.py b/components/isceobj/StripmapProc/runDenseOffsets.py index eb27366..0071b1e 100644 --- a/components/isceobj/StripmapProc/runDenseOffsets.py +++ b/components/isceobj/StripmapProc/runDenseOffsets.py @@ -62,7 +62,7 @@ def estimateOffsetField(master, slave, denseOffsetFileName, else: objOffset.setImageDataType2('real') - + objOffset.offsetImageName = denseOffsetFileName + '.bil' objOffset.snrImageName = denseOffsetFileName +'_snr.bil' objOffset.covImageName = denseOffsetFileName +'_cov.bil' @@ -75,11 +75,11 @@ def estimateOffsetField(master, slave, denseOffsetFileName, def runDenseOffsets(self): - if self.doDenseOffsets or self.doRubbersheeting: + if self.doDenseOffsets or self.doRubbersheetingAzimuth: if self.doDenseOffsets: print('Dense offsets explicitly requested') - if self.doRubbersheeting: + if self.doRubbersheetingAzimuth: print('Generating offsets as rubber sheeting requested') else: return @@ -96,7 +96,7 @@ def runDenseOffsets(self): os.makedirs(dirname) denseOffsetFilename = os.path.join(dirname , self.insar.denseOffsetFilename) - + field = estimateOffsetField(masterSlc, slaveSlc, denseOffsetFilename, ww = self.denseWindowWidth, wh = self.denseWindowHeight, @@ -107,5 +107,5 @@ def runDenseOffsets(self): self._insar.offset_top = field[0] self._insar.offset_left = field[1] - + return None diff --git a/components/isceobj/StripmapProc/runRubbersheet.py b/components/isceobj/StripmapProc/runRubbersheet.py index ea417bd..6f6afc8 100644 --- a/components/isceobj/StripmapProc/runRubbersheet.py +++ b/components/isceobj/StripmapProc/runRubbersheet.py @@ -13,13 +13,13 @@ def fill(data, invalid=None): """ Replace the value of invalid 'data' cells (indicated by 'invalid') by the value of the nearest valid data cell - + Input: data: numpy array of any dimension invalid: a binary array of same shape as 'data'. data value are replaced where invalid is True If None (default), use: invalid = np.isnan(data) - + Output: Return a filled array. """ @@ -97,7 +97,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): ###Currently making the assumption that top left of dense offsets and interfeorgrams are the same. ###This is not true for now. We need to update DenseOffsets to have the ability to have same top left ###As the input images. Once that is implemente, the math here should all be consistent. - ###However, this is not too far off since the skip for doing dense offsets is generally large. + ###However, this is not too far off since the skip for doing dense offsets is generally large. ###The offset is not too large to worry about right now. If the skip is decreased, this could be an issue. print('oversampling the filtered and masked offsets to the width and length:', width, ' ', length ) @@ -125,7 +125,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): for ll in range(length): val = geomoff.bands[0][ll,:] + osoff.bands[0][ll,:] val.tofile(fid) - + fid.close() img = isceobj.createImage() @@ -138,13 +138,13 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): img.scheme = 'BIP' img.renderHdr() - + return None def runRubbersheet(self): - if not self.doRubbersheeting: + if not self.doRubbersheetingAzimuth: print('Rubber sheeting not requested ... skipping') return @@ -168,11 +168,9 @@ def runRubbersheet(self): sheetOffset = os.path.join(offsetsDir, self.insar.azimuthRubbersheetFilename) # oversampling the filtAzOffsetFile to the same size of geometryAzimuthOffset - # and then update the geometryAzimuthOffset by adding the oversampled + # and then update the geometryAzimuthOffset by adding the oversampled # filtAzOffsetFile to it. resampleOffset(filtAzOffsetFile, geometryAzimuthOffset, sheetOffset) print("I'm here") return None - - From 6061697fb92854199f681955373a4ee6d83faba6 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Wed, 18 Dec 2019 17:00:38 -0800 Subject: [PATCH 16/27] stripmapStack: generate multilooked geometry files Stack: + add alks and rlks options in config.topo() topo: if number of looks is larger than 1, generate multilooked geometry files + generate `geom_master` folder in the root directory, same level as `Igrams` + import mroipac.looks.Looks module to multilook the full resolution geometry files from `./merged/geom_master` + copy over the full resolution xml/vrt files from `./merged/geome` to `./geom_master` as `*.full.xml/vrt` files, to support the number of looks extraction for mintpy + bug fixed in extractInfo(). This was not a problem before because the default values are 1 and the custom values from argument was not inputed in the config file. --- contrib/stack/stripmapStack/Stack.py | 2 + contrib/stack/stripmapStack/topo.py | 55 +++++++++++++++++++++++++++- 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/contrib/stack/stripmapStack/Stack.py b/contrib/stack/stripmapStack/Stack.py index b0a116e..818393e 100755 --- a/contrib/stack/stripmapStack/Stack.py +++ b/contrib/stack/stripmapStack/Stack.py @@ -65,6 +65,8 @@ class config(object): self.f.write('master : ' + self.slcDir +'\n') self.f.write('dem : ' + self.dem +'\n') self.f.write('output : ' + self.geometryDir +'\n') + self.f.write('alks : ' + self.alks +'\n') + self.f.write('rlks : ' + self.rlks +'\n') if self.nativeDoppler: self.f.write('native : True\n') if self.useGPU: diff --git a/contrib/stack/stripmapStack/topo.py b/contrib/stack/stripmapStack/topo.py index b8d7281..509a09b 100755 --- a/contrib/stack/stripmapStack/topo.py +++ b/contrib/stack/stripmapStack/topo.py @@ -6,8 +6,11 @@ import numpy as np import shelve import os import datetime +import shutil from isceobj.Constants import SPEED_OF_LIGHT from isceobj.Util.Poly2D import Poly2D +from mroipac.looks.Looks import Looks + def createParser(): ''' @@ -354,6 +357,47 @@ def runSimamp(outdir, hname='z.rdr'): simImage.renderHdr() hgtImage.finalizeImage() simImage.finalizeImage() + return + + +def runMultilook(in_dir, out_dir, alks, rlks): + print('generate multilooked geometry files with alks={} and rlks={}'.format(alks, rlks)) + from iscesys.Parsers.FileParserFactory import createFileParser + FP = createFileParser('xml') + + if not os.path.isdir(out_dir): + os.makedirs(out_dir) + print('create directory: {}'.format(out_dir)) + + for fbase in ['hgt', 'incLocal', 'lat', 'lon', 'los', 'shadowMask', 'waterMask']: + fname = '{}.rdr'.format(fbase) + in_file = os.path.join(in_dir, fname) + out_file = os.path.join(out_dir, fname) + + if os.path.isfile(in_file): + xmlProp = FP.parse(in_file+'.xml')[0] + if('image_type' in xmlProp and xmlProp['image_type'] == 'dem'): + inImage = isceobj.createDemImage() + else: + inImage = isceobj.createImage() + + inImage.load(in_file+'.xml') + inImage.filename = in_file + + lkObj = Looks() + lkObj.setDownLooks(alks) + lkObj.setAcrossLooks(rlks) + lkObj.setInputImage(inImage) + lkObj.setOutputFilename(out_file) + lkObj.looks() + + # copy the full resolution xml/vrt file from ./merged/geom_master to ./geom_master + # to facilitate the number of looks extraction + # the file path inside .xml file is not, but should, updated + shutil.copy(in_file+'.xml', out_file+'.full.xml') + shutil.copy(in_file+'.vrt', out_file+'.full.vrt') + + return out_dir def extractInfo(frame, inps): @@ -369,8 +413,8 @@ def extractInfo(frame, inps): info.lookSide = frame.instrument.platform.pointingDirection info.rangeFirstSample = frame.startingRange - info.numberRangeLooks = inps.rlks - info.numberAzimuthLooks = inps.alks + info.numberRangeLooks = 1 #inps.rlks + info.numberAzimuthLooks = 1 #inps.alks fsamp = frame.rangeSamplingRate @@ -443,6 +487,13 @@ def main(iargs=None): runTopo(info,demImage,dop=doppler,nativedop=inps.nativedop, legendre=inps.legendre) runSimamp(os.path.dirname(info.heightFilename),os.path.basename(info.heightFilename)) + # write multilooked geometry files in "geom_master" directory, same level as "Igrams" + if inps.rlks * inps.rlks > 1: + out_dir = os.path.join(os.path.dirname(os.path.dirname(info.outdir)), 'geom_master') + runMultilook(in_dir=info.outdir, out_dir=out_dir, alks=inps.alks, rlks=inps.rlks) + + return + if __name__ == '__main__': ''' From 851858b2284440f98d1188e23b711ec40038f7a7 Mon Sep 17 00:00:00 2001 From: HBaldwin3 <49163135+HBaldwin3@users.noreply.github.com> Date: Thu, 19 Dec 2019 19:13:43 -0600 Subject: [PATCH 17/27] Update runRubbersheetAzimuth.py Add "from scipy import ndimage" --- .../StripmapProc/runRubbersheetAzimuth.py | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py index 43ee4a8..fc07e25 100755 --- a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py +++ b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py @@ -22,39 +22,39 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): off_az = ds.GetRasterBand(1).ReadAsArray() off_rg = ds.GetRasterBand(2).ReadAsArray() ds = None - - # Remove missing values from ampcor + + # Remove missing values from ampcor off_rg[np.where(off_rg < -9999)]=0 off_az[np.where(off_az < -9999)]=0 - - + + # Store the offsets in a complex variable off = off_rg + 1j*off_az # Mask the azimuth offsets based on the MAD - mask = off_masking(off,filterSize,thre=3) - + mask = off_masking(off,filterSize,thre=3) + xoff_masked = np.ma.array(off.real,mask=mask) yoff_masked = np.ma.array(off.imag,mask=mask) - + # Delete unused variables mask = None off = None - + # Remove residual noisy spots with a median filter on the azimuth offmap yoff_masked.mask = yoff_masked.mask | \ (ndimage.median_filter(xoff_masked.filled(fill_value=0),3) == 0) | \ (ndimage.median_filter(yoff_masked.filled(fill_value=0),3) == 0) - - # Fill the data by iteratively using smoothed values + + # Fill the data by iteratively using smoothed values data = yoff_masked.data data[yoff_masked.mask]=np.nan - + off_az_filled = fill_with_smoothed(data,filterSize) - + # Apply median filter to smooth the azimuth offset map off_az_filled = ndimage.median_filter(off_az_filled,filterSize) - + # Save the filtered offsets length, width = off_az_filled.shape @@ -74,8 +74,8 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): img.scheme = 'BIP' img.renderHdr() - - return + + return def off_masking(off,filterSize,thre=2): @@ -94,16 +94,18 @@ def fill(data, invalid=None): """ Replace the value of invalid 'data' cells (indicated by 'invalid') by the value of the nearest valid data cell - + Input: data: numpy array of any dimension invalid: a binary array of same shape as 'data'. data value are replaced where invalid is True If None (default), use: invalid = np.isnan(data) - + Output: Return a filled array. """ + from scipy import ndimage + if invalid is None: invalid = np.isnan(data) ind = ndimage.distance_transform_edt(invalid, @@ -166,7 +168,7 @@ def fill_with_smoothed(off,filterSize): kernel = np.ones((filterSize,filterSize),np.float32)/(filterSize*filterSize) loop = 0 cnt2=1 - + while (cnt2!=0 & loop<100): loop += 1 idx2= np.isnan(off_2filt) @@ -178,9 +180,9 @@ def fill_with_smoothed(off,filterSize): idx3 = np.where(off_filt == 0) off_2filt[idx3]=np.nan off_filt=None - + return off_2filt - + def resampleOffset(maskedFiltOffset, geometryOffset, outName): ''' Oversample offset and add. @@ -198,7 +200,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): ###Currently making the assumption that top left of dense offsets and interfeorgrams are the same. ###This is not true for now. We need to update DenseOffsets to have the ability to have same top left ###As the input images. Once that is implemente, the math here should all be consistent. - ###However, this is not too far off since the skip for doing dense offsets is generally large. + ###However, this is not too far off since the skip for doing dense offsets is generally large. ###The offset is not too large to worry about right now. If the skip is decreased, this could be an issue. print('oversampling the filtered and masked offsets to the width and length:', width, ' ', length ) @@ -226,7 +228,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): for ll in range(length): val = geomoff.bands[0][ll,:] + osoff.bands[0][ll,:] val.tofile(fid) - + fid.close() img = isceobj.createImage() @@ -239,7 +241,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): img.scheme = 'BIP' img.renderHdr() - + return None @@ -264,7 +266,7 @@ def runRubbersheetAzimuth(self): if not self.doRubbersheetingRange: print('Rubber sheeting in range is off, filtering the offsets with a SNR-based mask') mask_filter(denseOffsetFile, snrFile, band[0], snrThreshold, filterSize, filtAzOffsetFile) - else: + else: print('Rubber sheeting in range is on, filtering the offsets with data-based mask') mask_filterNoSNR(denseOffsetFile, filterSize, filtAzOffsetFile) @@ -274,10 +276,8 @@ def runRubbersheetAzimuth(self): sheetOffset = os.path.join(offsetsDir, self.insar.azimuthRubbersheetFilename) # oversampling the filtAzOffsetFile to the same size of geometryAzimuthOffset - # and then update the geometryAzimuthOffset by adding the oversampled + # and then update the geometryAzimuthOffset by adding the oversampled # filtAzOffsetFile to it. resampleOffset(filtAzOffsetFile, geometryAzimuthOffset, sheetOffset) return None - - From 18bd583ec20e68704b28bdc97d7c009bba89510e Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Tue, 24 Dec 2019 16:12:34 -0800 Subject: [PATCH 18/27] topo: multilook geometry file with gdal_translate --- contrib/stack/stripmapStack/topo.py | 42 ++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/contrib/stack/stripmapStack/topo.py b/contrib/stack/stripmapStack/topo.py index 509a09b..b730af5 100755 --- a/contrib/stack/stripmapStack/topo.py +++ b/contrib/stack/stripmapStack/topo.py @@ -400,6 +400,45 @@ def runMultilook(in_dir, out_dir, alks, rlks): return out_dir +def runMultilookGdal(in_dir, out_dir, alks, rlks): + print('generate multilooked geometry files with alks={} and rlks={}'.format(alks, rlks)) + import gdal + + # create 'geom_master' directory + if not os.path.isdir(out_dir): + os.makedirs(out_dir) + print('create directory: {}'.format(out_dir)) + + # multilook files one by one + for fbase in ['hgt', 'incLocal', 'lat', 'lon', 'los', 'shadowMask', 'waterMask']: + fname = '{}.rdr'.format(fbase) + in_file = os.path.join(in_dir, fname) + out_file = os.path.join(out_dir, fname) + + if os.path.isfile(in_file): + ds = gdal.Open(in_file, gdal.GA_ReadOnly) + in_wid = ds.RasterXSize + in_len = ds.RasterYSize + + out_wid = int(in_wid / rlks) + out_len = int(in_len / alks) + src_wid = out_wid * rlks + src_len = out_len * alks + + cmd = 'gdal_translate -of ENVI -a_nodata 0 -outsize {ox} {oy} '.format(ox=out_wid, oy=out_len) + cmd += ' -srcwin 0 0 {sx} {sy} {fi} {fo} '.format(sx=src_wid, sy=src_len, fi=in_file, fo=out_file) + print(cmd) + os.system(cmd) + + # copy the full resolution xml/vrt file from ./merged/geom_master to ./geom_master + # to facilitate the number of looks extraction + # the file path inside .xml file is not, but should, updated + shutil.copy(in_file+'.xml', out_file+'.full.xml') + shutil.copy(in_file+'.vrt', out_file+'.full.vrt') + + return out_dir + + def extractInfo(frame, inps): ''' Extract relevant information only. @@ -490,7 +529,8 @@ def main(iargs=None): # write multilooked geometry files in "geom_master" directory, same level as "Igrams" if inps.rlks * inps.rlks > 1: out_dir = os.path.join(os.path.dirname(os.path.dirname(info.outdir)), 'geom_master') - runMultilook(in_dir=info.outdir, out_dir=out_dir, alks=inps.alks, rlks=inps.rlks) + runMultilookGdal(in_dir=info.outdir, out_dir=out_dir, alks=inps.alks, rlks=inps.rlks) + #runMultilook(in_dir=info.outdir, out_dir=out_dir, alks=inps.alks, rlks=inps.rlks) return From 0b3d49744d5b918cfcc3dea50417ebdf24f5ff1f Mon Sep 17 00:00:00 2001 From: Heresh Fattahi Date: Thu, 26 Dec 2019 23:02:35 -0800 Subject: [PATCH 19/27] fixing bug in naming interferogram without rubbersheeting --- applications/stripmapApp.py | 2 +- .../isceobj/StripmapProc/runInterferogram.py | 16 ++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/applications/stripmapApp.py b/applications/stripmapApp.py index 938cafb..6f76b66 100755 --- a/applications/stripmapApp.py +++ b/applications/stripmapApp.py @@ -265,7 +265,7 @@ RUBBERSHEET_SNR_THRESHOLD = Application.Parameter('rubberSheetSNRThreshold', RUBBERSHEET_FILTER_SIZE = Application.Parameter('rubberSheetFilterSize', public_name='rubber sheet filter size', - default = 8, + default = 9, type = int, mandatory = False, doc = '') diff --git a/components/isceobj/StripmapProc/runInterferogram.py b/components/isceobj/StripmapProc/runInterferogram.py index fbbde4f..3d43136 100644 --- a/components/isceobj/StripmapProc/runInterferogram.py +++ b/components/isceobj/StripmapProc/runInterferogram.py @@ -56,7 +56,7 @@ def compute_FlatEarth(self,ifgFilename,width,length,radarWavelength): # Open the interferogram #ifgFilename= os.path.join(self.insar.ifgDirname, self.insar.ifgFilename) - intf = np.memmap(ifgFilename+'.full',dtype=np.complex64,mode='r+',shape=(length,width)) + intf = np.memmap(ifgFilename,dtype=np.complex64,mode='r+',shape=(length,width)) for ll in range(length): intf[ll,:] *= np.exp(cJ*fact*rng2[ll,:]) @@ -155,10 +155,14 @@ def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks,radarW else: resampAmp += '.amp' - resampInt = resampName + if not self.doRubbersheetingRange: + resampInt = resampName + else: + resampInt = resampName + ".full" + #resampAmp = resampAmp + ".full" objInt = isceobj.createIntImage() - objInt.setFilename(resampInt+'.full') + objInt.setFilename(resampInt) objInt.setWidth(intWidth) imageInt = isceobj.createIntImage() IU.copyAttributes(objInt, imageInt) @@ -166,7 +170,7 @@ def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks,radarW objInt.createImage() objAmp = isceobj.createAmpImage() - objAmp.setFilename(resampAmp+'.full') + objAmp.setFilename(resampAmp) objAmp.setWidth(intWidth) imageAmp = isceobj.createAmpImage() IU.copyAttributes(objAmp, imageAmp) @@ -196,8 +200,8 @@ def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks,radarW compute_FlatEarth(self,resampInt,intWidth,lines,radarWavelength) # Perform Multilook - multilook(resampInt+'.full', outname=resampInt, alks=azLooks, rlks=rgLooks) #takeLooks(objAmp,azLooks,rgLooks) - multilook(resampAmp+'.full', outname=resampAmp, alks=azLooks, rlks=rgLooks) #takeLooks(objInt,azLooks,rgLooks) + multilook(resampInt, outname=resampName, alks=azLooks, rlks=rgLooks) #takeLooks(objAmp,azLooks,rgLooks) + multilook(resampAmp, outname=resampAmp.replace(".full",""), alks=azLooks, rlks=rgLooks) #takeLooks(objInt,azLooks,rgLooks) #os.system('rm ' + resampInt+'.full* ' + resampAmp + '.full* ') # End of modification From 0704e98ac4f87246587bc656d10cf3f554ae7b49 Mon Sep 17 00:00:00 2001 From: Heresh Fattahi Date: Thu, 26 Dec 2019 23:08:17 -0800 Subject: [PATCH 20/27] remove a commented out line --- components/isceobj/StripmapProc/runInterferogram.py | 1 - 1 file changed, 1 deletion(-) diff --git a/components/isceobj/StripmapProc/runInterferogram.py b/components/isceobj/StripmapProc/runInterferogram.py index 3d43136..da42508 100644 --- a/components/isceobj/StripmapProc/runInterferogram.py +++ b/components/isceobj/StripmapProc/runInterferogram.py @@ -159,7 +159,6 @@ def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks,radarW resampInt = resampName else: resampInt = resampName + ".full" - #resampAmp = resampAmp + ".full" objInt = isceobj.createIntImage() objInt.setFilename(resampInt) From b2f320459531e07e233450a15f0e9981e4c4f302 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Mon, 18 Mar 2019 19:20:18 -0400 Subject: [PATCH 21/27] stripmapStack/createWaterMask: add --dem_file opt add -d / --dem_file option to input DEM file as an alternative for --bbox option. add dem2bbox() to grab bbox from existing DEM file adjust download_waterMask() that: if DEM file is input, put downloaded water body mask file to the same directory as the DEM file. adjust the default water mask file in radar coordinates that: 1. in the same directory as the input lat/lon files 2. named as "waterMask.rdr" instead of "swbdLat_*.wbd.rdr" add example docstring and update parser description. --- .../stack/stripmapStack/createWaterMask.py | 89 ++++++++++++++----- 1 file changed, 66 insertions(+), 23 deletions(-) diff --git a/contrib/stack/stripmapStack/createWaterMask.py b/contrib/stack/stripmapStack/createWaterMask.py index 8ffa9fa..97f9875 100755 --- a/contrib/stack/stripmapStack/createWaterMask.py +++ b/contrib/stack/stripmapStack/createWaterMask.py @@ -2,28 +2,42 @@ #Author: Heresh Fattahi -import isce -import isceobj -from contrib.demUtils.SWBDStitcher import SWBDStitcher -from iscesys.DataManager import createManager +import os import argparse import configparser -from numpy import round +import numpy as np +import isce +import isceobj +from iscesys.DataManager import createManager +from contrib.demUtils.SWBDStitcher import SWBDStitcher + + +EXAMPLE = """example: + createWaterMask.py -b 31 33 130 132 + createWaterMask.py -b 31 33 130 132 -l lat.rdr -L lon.rdr -o waterMask.rdr + createWaterMask.py -d ../DEM/demLat_N31_N33_Lon_E130_E132.dem.wgs84 -l lat.rdr -L lon.rdr -o waterMask.rdr +""" def createParser(): ''' Create command line parser. ''' - parser = argparse.ArgumentParser( description='extracts the overlap geometry between master bursts') - # parser.add_argument('-b', '--bbox', dest='bbox', type=str, default=None, - # help='Lat/Lon Bounding SNWE') - parser.add_argument('-b', '--bbox', type = int, default = None, nargs = '+', dest = 'bbox', help = 'Defines the spatial region in the format south north west east.\ - The values should be integers from (-90,90) for latitudes and (0,360) or (-180,180) for longitudes.') + parser = argparse.ArgumentParser(description='Create water body mask in geo and/or radar coordinates', + formatter_class=argparse.RawTextHelpFormatter, + epilog=EXAMPLE) + parser.add_argument('-b', '--bbox', dest='bbox', type=int, default=None, nargs=4, metavar=('S','N','W','E'), + help = 'Defines the spatial region in the format south north west east.\n' + 'The values should be integers from (-90,90) for latitudes ' + 'and (0,360) or (-180,180) for longitudes.') + parser.add_argument('-d','--dem_file', dest='demName', type=str, default=None, + help='DEM file in geo coordinates, i.e. demLat*.dem.wgs84.') parser.add_argument('-l', '--lat_file', dest='latName', type=str, default=None, help='pixel by pixel lat file in radar coordinate') parser.add_argument('-L', '--lon_file', dest='lonName', type=str, default=None, help='pixel by pixel lat file in radar coordinate') + parser.add_argument('-o', '--output', dest='outfile', type=str, + help='output filename of water mask in radar coordinates') return parser def cmdLineParse(iargs = None): @@ -33,26 +47,54 @@ def cmdLineParse(iargs = None): parser = createParser() inps = parser.parse_args(args=iargs) - #inps.bbox = [int(round(val)) for val in inps.bbox.split()] + + if not inps.bbox and not inps.demName: + parser.print_usage() + raise SystemExit('ERROR: no --bbox/--dem_file input, at least one is required.') + + if not inps.outfile and (inps.latName and inps.lonName): + inps.outfile = os.path.join(os.path.dirname(inps.latName), 'waterMask.rdr') + return inps -def download_waterMask(inps): +def dem2bbox(dem_file): + """Grab bbox from DEM file in geo coordinates""" + demImage = isceobj.createDemImage() + demImage.load(dem_file + '.xml') + demImage.setAccessMode('read') + N = demImage.getFirstLatitude() + W = demImage.getFirstLongitude() + S = N + demImage.getDeltaLatitude() * demImage.getLength() + E = W + demImage.getDeltaLongitude() * demImage.getWidth() + bbox = [np.floor(S).astype(int), np.ceil(N).astype(int), + np.floor(W).astype(int), np.ceil(E).astype(int)] + return bbox + + +def download_waterMask(bbox, dem_file): + out_dir = os.getcwd() + # update out_dir and/or bbox if dem_file is input + if dem_file: + out_dir = os.path.dirname(dem_file) + if not bbox: + bbox = dem2bbox(dem_file) sw = createManager('wbd') sw.configure() - inps.waterBodyGeo = sw.defaultName(inps.bbox) + #inps.waterBodyGeo = sw.defaultName(inps.bbox) + sw.outputFile = os.path.join(out_dir, sw.defaultName(bbox)) sw._noFilling = False #sw._fillingValue = -1.0 sw._fillingValue = 0.0 - sw.stitch(inps.bbox[0:2],inps.bbox[2:]) + sw.stitch(bbox[0:2], bbox[2:]) + return sw.outputFile - return inps - -def geo2radar(inps): - inps.waterBodyRadar = inps.waterBodyGeo + '.rdr' +def geo2radar(geo_file, rdr_file, lat_file, lon_file): + #inps.waterBodyRadar = inps.waterBodyGeo + '.rdr' sw = SWBDStitcher() - sw.toRadar(inps.waterBodyGeo, inps.latName, inps.lonName, inps.waterBodyRadar) + sw.toRadar(geo_file, lat_file, lon_file, rdr_file) + return rdr_file #looks.py -i watermask.msk -r 4 -a 14 -o 'waterMask.14alks_4rlks.msk' @@ -60,10 +102,11 @@ def geo2radar(inps): def main(iargs=None): - inps = cmdLineParse(iargs) - inps = download_waterMask(inps) - if inps.latName and inps.lonName: - inps = geo2radar(inps) + inps = cmdLineParse(iargs) + geo_file = download_waterMask(inps.bbox, inps.demName) + if inps.latName and inps.lonName: + geo2radar(geo_file, inps.outfile, inps.latName, inps.lonName) + return if __name__ == '__main__' : ''' From ad51158a442b045fa8beea3f50c0781363541ec8 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Tue, 31 Dec 2019 16:00:58 -0800 Subject: [PATCH 22/27] stripmapStack: add function createWaterMask to config_master_* Stack.py: add function call to create water mask file in radar coordinates in the step one "run_1_master_*" createWaterMask.py: fill pixels without DEM data with value of -1, same as water body --- contrib/stack/stripmapStack/Stack.py | 34 ++++++++++++++----- .../stack/stripmapStack/createWaterMask.py | 8 +++-- 2 files changed, 32 insertions(+), 10 deletions(-) diff --git a/contrib/stack/stripmapStack/Stack.py b/contrib/stack/stripmapStack/Stack.py index 818393e..0caa9b2 100755 --- a/contrib/stack/stripmapStack/Stack.py +++ b/contrib/stack/stripmapStack/Stack.py @@ -74,7 +74,18 @@ class config(object): else: self.f.write('useGPU : False\n') self.f.write('##########################'+'\n') - + + def createWaterMask(self, function): + + self.f.write('##########################'+'\n') + self.f.write(function+'\n') + self.f.write('createWaterMask : '+'\n') + self.f.write('dem_file : ' + self.dem +'\n') + self.f.write('lat_file : ' + self.latFile +'\n') + self.f.write('lon_file : ' + self.lonFile +'\n') + self.f.write('output : ' + self.waterMaskFile + '\n') + self.f.write('##########################'+'\n') + def geo2rdr(self, function): self.f.write('##########################'+'\n') @@ -311,8 +322,7 @@ class run(object): self.runf.write(self.text_cmd+'stripmapWrapper.py -c '+ configName+'\n') def master_focus_split_geometry(self, stackMaster, config_prefix, split=False, focus=True, native=True): - ######## - # focusing master and producing geometry files + """focusing master and producing geometry files""" configName = os.path.join(self.configDir, config_prefix + stackMaster) configObj = config(configName) configObj.configure(self) @@ -329,11 +339,19 @@ class run(object): counter += 1 if split: - configObj.slc = os.path.join(configObj.slcDir,stackMaster+self.raw_string+'.slc') - configObj.outDir = configObj.slcDir - configObj.shelve = os.path.join(configObj.slcDir, 'data') - configObj.splitRangeSpectrum('[Function-{0}]'.format(counter)) - + configObj.slc = os.path.join(configObj.slcDir,stackMaster+self.raw_string+'.slc') + configObj.outDir = configObj.slcDir + configObj.shelve = os.path.join(configObj.slcDir, 'data') + configObj.splitRangeSpectrum('[Function-{0}]'.format(counter)) + counter += 1 + + # generate water mask in radar coordinates + configObj.latFile = os.path.join(self.workDir, 'geom_master/lat.rdr') + configObj.lonFile = os.path.join(self.workDir, 'geom_master/lon.rdr') + configObj.waterMaskFile = os.path.join(self.workDir, 'geom_master/waterMask.rdr') + configObj.createWaterMask('[Function-{0}]'.format(counter)) + counter += 1 + configObj.finalize() del configObj self.runf.write(self.text_cmd+'stripmapWrapper.py -c '+ configName+'\n') diff --git a/contrib/stack/stripmapStack/createWaterMask.py b/contrib/stack/stripmapStack/createWaterMask.py index 97f9875..5946655 100755 --- a/contrib/stack/stripmapStack/createWaterMask.py +++ b/contrib/stack/stripmapStack/createWaterMask.py @@ -40,6 +40,7 @@ def createParser(): help='output filename of water mask in radar coordinates') return parser + def cmdLineParse(iargs = None): ''' Command line parser. @@ -85,11 +86,12 @@ def download_waterMask(bbox, dem_file): #inps.waterBodyGeo = sw.defaultName(inps.bbox) sw.outputFile = os.path.join(out_dir, sw.defaultName(bbox)) sw._noFilling = False - #sw._fillingValue = -1.0 - sw._fillingValue = 0.0 + sw._fillingValue = -1.0 #fill pixels without DEM data with value of -1, same as water body + #sw._fillingValue = 0.0 sw.stitch(bbox[0:2], bbox[2:]) return sw.outputFile + def geo2radar(geo_file, rdr_file, lat_file, lon_file): #inps.waterBodyRadar = inps.waterBodyGeo + '.rdr' sw = SWBDStitcher() @@ -100,6 +102,7 @@ def geo2radar(geo_file, rdr_file, lat_file, lon_file): #imageMath.py -e='a*b' --a=filt_20100911_20101027.int --b=watermask.14alks_4rlks.msk -o filt_20100911_20101027_masked.int -t cfloat -s BIL + def main(iargs=None): inps = cmdLineParse(iargs) @@ -108,6 +111,7 @@ def main(iargs=None): geo2radar(geo_file, inps.outfile, inps.latName, inps.lonName) return + if __name__ == '__main__' : ''' creates a water mask and transforms to radar Coord From 1b70f52f914d1e16d4ed0eef67ce23b99e6444f4 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Mon, 18 Mar 2019 19:32:08 -0400 Subject: [PATCH 23/27] stripmapStack: convert README.txt to README.md + convert README.txt to README.md + stackStripMap.py: group iono options in createParser() - group all iono options in createParser() for more easily navigation in --help. - fix typos in the descriptions for the script and --dem option. --- contrib/stack/stripmapStack/README.md | 85 ++++++++++++++++++++ contrib/stack/stripmapStack/README.txt | 64 --------------- contrib/stack/stripmapStack/stackStripMap.py | 65 ++++++++------- contrib/stack/stripmapStack/topo.py | 63 +++++++-------- 4 files changed, 151 insertions(+), 126 deletions(-) create mode 100644 contrib/stack/stripmapStack/README.md delete mode 100644 contrib/stack/stripmapStack/README.txt diff --git a/contrib/stack/stripmapStack/README.md b/contrib/stack/stripmapStack/README.md new file mode 100644 index 0000000..4bf9260 --- /dev/null +++ b/contrib/stack/stripmapStack/README.md @@ -0,0 +1,85 @@ +## StripMap stack processor + +The detailed algorithms and workflow for stack processing of stripmap SAR data can be found here: + ++ Fattahi, H., M. Simons, and P. Agram (2017), InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique, IEEE Transactions on Geoscience and Remote Sensing, 55(10), 5984-5996, doi:[10.1109/TGRS.2017.2718566](https://ieeexplore.ieee.org/abstract/document/7987747/). + +----------------------------------- + +To use the stripmap stack processor, make sure to add the path of your `contrib/stack/stripmapStack` folder to your `$PATH` environment varibale. + +Currently supported workflows include a coregistered stack of SLC, interferograms, ionospheric delays. + +Here are some notes to get started with processing stacks of stripmap data with ISCE. + +#### 1. Create your project folder somewhere + +``` +mkdir MauleAlosDT111 +cd MauleAlosDT111 +``` + +#### 2. Prepare DEM +a) create a folder for DEM; +b) create a DEM using dem.py with SNWE of your study area in integer; +d) Keep only ".dem.wgs84", ".dem.wgs84.vrt" and ".dem.wgs84.xml" and remove unnecessary files; +d) fix the path of the file in the xml file of the DEM by using fixImageXml.py. + +``` +mkdir DEM; cd DEM +dem.py -a stitch -b -37 -31 -72 -69 -r -s 1 -c +rm demLat*.dem demLat*.dem.xml demLat*.dem.vrt +fixImageXml.py -f -i demLat*.dem.wgs84 +cd .. +``` + +#### 3. Download data + +##### 3.1 create a folder to download SAR data (i.e. ALOS-1 data from ASF) + +``` +mkdir download +cd download +``` + +##### 3.2 Download the data that that you want to process to the "downlowd" directory + +#### 4. Prepare SAR data + +Once all data have been downloaded, we need to unzip them and move them to different folders and getting ready for unpacking and then SLC generation. This can be done by running the following command in a directory above "download": + +``` +prepRawALOS.py -i download/ -o SLC +``` + +This command generates an empty SLC folder and a run file called: "run_unPackALOS". + +You could also run prepRawSensor.py which aims to recognize the sensor data automatically followed by running the sensor specific preparation script. For now we include support for ALOS and CSK raw data, but it is trivial to expand and include other sensors as unpacking routines are already included in the distribution. + +``` +prepRawSensor.py -i download/ -o SLC +``` + +#### 5. Execute the commands in "run_unPackALOS" file + +If you have a cluster that you can submit jobs, you can submit each line of command to a processor. The commands are independent and can be run in parallel. + +After successfully running the previous step, you should see acquisition dates in the SLC folder and the ".raw" files for each acquisition. + +Note: For ALOS-1, If there is an acquisition that does not include .raw file, this is most likely due to PRF change between frames and can not be currently handled by ISCE. You have to ignore those. + +#### 6. Run "stackStripmap.py" + +This will generate many config and run files that need to be executed. Here is an example: + +``` +stackStripMap.py -s SLC/ -d DEM/demLat*.dem.wgs84 -t 250 -b 1000 -a 14 -r 4 -u snaphu +``` + +This will produce: +a) baseline folder, which contains baseline information +b) pairs.png which is a baseline-time plot of the network of interferograms +c) configs: which contains the configuration parameter to run different InSAR processing steps +d) run_files: a folder that includes several run and job files that needs to be run in order + +#### 7. Execute the commands in run files (run_1*, run_2*, etc) in the "run_files" folder diff --git a/contrib/stack/stripmapStack/README.txt b/contrib/stack/stripmapStack/README.txt deleted file mode 100644 index e6ae915..0000000 --- a/contrib/stack/stripmapStack/README.txt +++ /dev/null @@ -1,64 +0,0 @@ - -The detailed algorithms for stack processing of stripmap data can be found here: - -H. Fattahi, M. Simons, and P. Agram, "InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique", IEEE Trans. Geosci. Remote Sens., vol. 55, no. 10, 5984-5996, 2017. (https://ieeexplore.ieee.org/abstract/document/7987747/) - - ------------------------------------ - -Notes on stripmap stack processor: - -Here are some notes to get started with processing stacks of stripmap data with ISCE. - - -1- create a folder somewhere for your project - -mkdir MauleT111 -cd MauleT111 - -2- create a DEM: - -dem.py -a stitch -b -37 -31 -72 -69 -r -s 1 -c - -3- Keep only ".dem.wgs84", ".dem.wgs84.vrt" and ".dem.wgs84.xml" and remove unnecessary files - -4- fix the path of the file in the xml file of the DEM by using this command: - -fixImageXml.py -f -i demLat_S37_S31_Lon_W072_W069.dem.wgs84 - -5- create a folder to download the ALOS-1 data from ASF: - -mkdir download -cd download - -6- Download the data that that you want to process to the downlowd directory. - -7- once all data have been downloaded, we need to unzip them and move them to different folders and getting ready for unpacking and then SLC generation. -This can be done by running the following command in a directory above "download": - -prepRawALOS.py -i download/ -o SLC - -This command generates an empty SLC folder and a run file called: "run_unPackALOS". -You could also run prepRawSensor.py which aims to recognize the sensor data automatically followed by running the sensor specific preparation script. For now we include support for ALOS and CSK raw data, but it is trivial to expand and include other sensors as unpacking routines are already included in the distribution. - -prepRawSensor.py -i download/ -o SLC - -8- execute the commands inside run_unPackALOS file. If you have a cluster that you can submit jobs, you can submit each line of command to a processor. The commands are independent and can be run in parallel. - -9- After successfully running the previous step, you should see acquisition dates in the SLC folder and the ".raw" files for each acquisition - -Note: For ALOS-1, If there is an acquisition that does not include .raw file, this is most likely due to PRF change between frames and can not be currently handled by ISCE. You have to ignore those. - -10- run stackStripmap.py which will generate many config and run files that need to be executed. Here is an example: - -stackStripMap.py -s SLC/ -d demLat_S37_S31_Lon_W072_W069.dem.wgs84 -t 250 -b 1000 -a 14 -r 4 -u snaphu - -This will produce: -a) baseline folder, which contains baseline information -b) pairs.png which is a baseline-time plot of the network of interferograms -c) configs: which contains the configuration parameter to run different InSAR processing steps -d) run_files: a folder that includes several run and job files that needs to be run in order - -11- execute the commands in run files (run_1, run_2, etc) in the run_files folder - - diff --git a/contrib/stack/stripmapStack/stackStripMap.py b/contrib/stack/stripmapStack/stackStripMap.py index 07c1f68..728ffe0 100755 --- a/contrib/stack/stripmapStack/stackStripMap.py +++ b/contrib/stack/stripmapStack/stackStripMap.py @@ -5,7 +5,7 @@ import os, imp, sys, glob import argparse import configparser -import datetime +import datetime import numpy as np import shelve import isce @@ -20,7 +20,7 @@ defoMax = '2' maxNodes = 72 def createParser(): - parser = argparse.ArgumentParser( description='Preparing the directory structure and config files for stack processing of Sentinel data') + parser = argparse.ArgumentParser( description='Preparing the directory structure and config files for stack processing of StripMap data') parser.add_argument('-s', '--slc_directory', dest='slcDir', type=str, required=True, help='Directory with all stripmap SLCs') @@ -31,7 +31,7 @@ def createParser(): help='Working directory ') parser.add_argument('-d', '--dem', dest='dem', type=str, required=True, - help='Directory with the DEM (.xml and .vrt files)') + help='DEM file (with .xml and .vrt files)') parser.add_argument('-m', '--master_date', dest='masterDate', type=str, default=None, help='Directory with master acquisition') @@ -43,47 +43,54 @@ def createParser(): help='Baseline threshold (max bperp in meters)') parser.add_argument('-a', '--azimuth_looks', dest='alks', type=str, default='10', - help='Number of looks in azimuth (automaticly computed as AspectR*looks when "S" or "sensor" is defined to give approximately square multi-look pixels)') + help='Number of looks in azimuth (automaticly computed as AspectR*looks when ' + '"S" or "sensor" is defined to give approximately square multi-look pixels)') parser.add_argument('-r', '--range_looks', dest='rlks', type=str, default='10', help='Number of looks in range') parser.add_argument('-S', '--sensor', dest='sensor', type=str, required=False, help='SAR sensor used to define square multi-look pixels') - parser.add_argument('-L', '--low_band_frequency', dest='fL', type=str, default=None, - help='low band frequency') - parser.add_argument('-H', '--high_band_frequency', dest='fH', type=str, default=None, - help='high band frequency') - parser.add_argument('-B', '--subband_bandwidth ', dest='bandWidth', type=str, default=None, - help='sub-band band width') - parser.add_argument('-u', '--unw_method', dest='unwMethod', type=str, default='snaphu' - , help='unwrapping method (icu, snaphu, or snaphu2stage)') + + parser.add_argument('-u', '--unw_method', dest='unwMethod', type=str, default='snaphu', + help='unwrapping method (icu, snaphu, or snaphu2stage)') parser.add_argument('-f','--filter_strength', dest='filtStrength', type=str, default=filtStrength, help='strength of Goldstein filter applied to the wrapped phase before spatial coherence estimation.' ' Default: {}'.format(filtStrength)) - parser.add_argument('--filter_sigma_x', dest='filterSigmaX', type=str, default='100' - , help='filter sigma for gaussian filtering the dispersive and nonDispersive phase') + iono = parser.add_argument_group('Ionosphere', 'Configurationas for ionospheric correction') + iono.add_argument('-L', '--low_band_frequency', dest='fL', type=str, default=None, + help='low band frequency') + iono.add_argument('-H', '--high_band_frequency', dest='fH', type=str, default=None, + help='high band frequency') + iono.add_argument('-B', '--subband_bandwidth ', dest='bandWidth', type=str, default=None, + help='sub-band band width') - parser.add_argument('--filter_sigma_y', dest='filterSigmaY', type=str, default='100.0', - help='sigma of the gaussian filter in Y direction, default=100') + iono.add_argument('--filter_sigma_x', dest='filterSigmaX', type=str, default='100', + help='filter sigma for gaussian filtering the dispersive and nonDispersive phase') - parser.add_argument('--filter_size_x', dest='filterSizeX', type=str, default='800.0', - help='size of the gaussian kernel in X direction, default = 800') + iono.add_argument('--filter_sigma_y', dest='filterSigmaY', type=str, default='100.0', + help='sigma of the gaussian filter in Y direction, default=100') - parser.add_argument('--filter_size_y', dest='filterSizeY', type=str, default='800.0', - help='size of the gaussian kernel in Y direction, default=800') + iono.add_argument('--filter_size_x', dest='filterSizeX', type=str, default='800.0', + help='size of the gaussian kernel in X direction, default = 800') - parser.add_argument('--filter_kernel_rotation', dest='filterKernelRotation', type=str, default='0.0', - help='rotation angle of the filter kernel in degrees (default = 0.0)') + iono.add_argument('--filter_size_y', dest='filterSizeY', type=str, default='800.0', + help='size of the gaussian kernel in Y direction, default=800') - parser.add_argument('-W', '--workflow', dest='workflow', type=str, default='slc' - , help='The InSAR processing workflow : (slc, interferogram, ionosphere)') + iono.add_argument('--filter_kernel_rotation', dest='filterKernelRotation', type=str, default='0.0', + help='rotation angle of the filter kernel in degrees (default = 0.0)') - parser.add_argument('-z', '--zero', dest='zerodop', action='store_true', default=False, help='Use zero doppler geometry for processing - Default : No') - parser.add_argument('--nofocus', dest='nofocus', action='store_true', default=False, help='If input data is already focused to SLCs - Default : do focus') - parser.add_argument('-c', '--text_cmd', dest='text_cmd', type=str, default='' - , help='text command to be added to the beginning of each line of the run files. Example : source ~/.bash_profile;') - parser.add_argument('-useGPU', '--useGPU', dest='useGPU',action='store_true', default=False, help='Allow App to use GPU when available') + parser.add_argument('-W', '--workflow', dest='workflow', type=str, default='slc', + help='The InSAR processing workflow : (slc, interferogram, ionosphere)') + + parser.add_argument('-z', '--zero', dest='zerodop', action='store_true', default=False, + help='Use zero doppler geometry for processing - Default : No') + parser.add_argument('--nofocus', dest='nofocus', action='store_true', default=False, + help='If input data is already focused to SLCs - Default : do focus') + parser.add_argument('-c', '--text_cmd', dest='text_cmd', type=str, default='', + help='text command to be added to the beginning of each line of the run files. Example : source ~/.bash_profile;') + parser.add_argument('-useGPU', '--useGPU', dest='useGPU',action='store_true', default=False, + help='Allow App to use GPU when available') parser.add_argument('--summary', dest='summary', action='store_true', default=False, help='Show summary only') return parser diff --git a/contrib/stack/stripmapStack/topo.py b/contrib/stack/stripmapStack/topo.py index b730af5..2b175a6 100755 --- a/contrib/stack/stripmapStack/topo.py +++ b/contrib/stack/stripmapStack/topo.py @@ -1,17 +1,17 @@ #!/usr/bin/env python3 + +import os import argparse +import shelve +import datetime +import shutil +import numpy as np import isce import isceobj -import numpy as np -import shelve -import os -import datetime -import shutil from isceobj.Constants import SPEED_OF_LIGHT from isceobj.Util.Poly2D import Poly2D from mroipac.looks.Looks import Looks - def createParser(): ''' Command line parser. @@ -48,7 +48,7 @@ class Dummy(object): def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): - + from isceobj.Planet.Planet import Planet from zerodop.GPUtopozero.GPUtopozero import PyTopozero from isceobj import Constants as CN @@ -84,14 +84,14 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): omethod = 2 # LEGENDRE INTERPOLATION else: omethod = 0 # HERMITE INTERPOLATION - + # tracking doppler specifications if nativedop and (dop is not None): try: coeffs = dop._coeffs except: coeffs = dop - + polyDoppler = Poly2D() polyDoppler.setWidth(width) polyDoppler.setLength(length) @@ -109,10 +109,10 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): polyDoppler.createPoly2D() - # dem + # dem demImage.setCaster('read','FLOAT') demImage.createImage() - + # slant range file slantRangeImage = Poly2D() slantRangeImage.setWidth(width) @@ -130,12 +130,12 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): dataType = 'DOUBLE' latImage.initImage(latFilename,accessMode,width,dataType) latImage.createImage() - + # lon file lonImage = isceobj.createImage() lonImage.initImage(lonFilename,accessMode,width,dataType) lonImage.createImage() - + # LOS file losImage = isceobj.createImage() dataType = 'FLOAT' @@ -144,7 +144,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): losImage.initImage(losFilename,accessMode,width,dataType,bands=bands,scheme=scheme) losImage.setCaster('write','DOUBLE') losImage.createImage() - + # height file heightImage = isceobj.createImage() dataType = 'DOUBLE' @@ -158,7 +158,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): incImage.initImage(incFilename,accessMode,width,dataType,bands=bands,scheme=scheme) incImage.createImage() incImagePtr = incImage.getImagePointer() - + maskImage = isceobj.createImage() dataType = 'BYTE' bands = 1 @@ -168,7 +168,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): else: incImagePtr = 0 maskImagePtr = 0 - + # initalize planet elp = Planet(pname='Earth').ellipsoid @@ -214,14 +214,14 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): topo.set_orbitBasis(1) # Is this ever different? topo.createOrbit() # Initializes the empty orbit to the right allocated size count = 0 - + for sv in info.orbit.stateVectors.list: td = DTU.seconds_since_midnight(sv.getTime()) pos = sv.getPosition() vel = sv.getVelocity() topo.set_orbitVector(count,td,pos[0],pos[1],pos[2],vel[0],vel[1],vel[2]) count += 1 - + # run topo topo.runTopo() @@ -244,13 +244,13 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): # los file descr = '''Two channel Line-Of-Sight geometry image (all angles in degrees). Represents vector drawn from target to platform. - Channel 1: Incidence angle measured from vertical at target (always +ve). + Channel 1: Incidence angle measured from vertical at target (always +ve). Channel 2: Azimuth angle measured from North in Anti-clockwise direction.''' losImage.setImageType('bil') losImage.addDescription(descr) losImage.finalizeImage() losImage.renderHdr() - + # dem/ height file demImage.finalizeImage() @@ -259,7 +259,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): descr = '''Two channel angle file. Channel 1: Angle between ray to target and the vertical at the sensor Channel 2: Local incidence angle accounting for DEM slope at target''' - + incImage.addDescription(descr) incImage.finalizeImage() incImage.renderHdr() @@ -268,7 +268,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): maskImage.addDescription(descr) maskImage.finalizeImage() maskImage.renderHdr() - + if slantRangeImage: try: slantRangeImage.finalizeImage() @@ -276,7 +276,7 @@ def runTopoGPU(info, demImage, dop=None, nativedop=False, legendre=False): pass -def runTopoCPU(info, demImage, dop=None, +def runTopoCPU(info, demImage, dop=None, nativedop=False, legendre=False): from zerodop.topozero import createTopozero from isceobj.Planet.Planet import Planet @@ -298,7 +298,7 @@ def runTopoCPU(info, demImage, dop=None, topo.numberRangeLooks = info.numberRangeLooks topo.numberAzimuthLooks = info.numberAzimuthLooks topo.lookSide = info.lookSide - topo.sensingStart = info.sensingStart + datetime.timedelta(seconds = ((info.numberAzimuthLooks - 1) /2) / info.prf) + topo.sensingStart = info.sensingStart + datetime.timedelta(seconds = ((info.numberAzimuthLooks - 1) /2) / info.prf) topo.rangeFirstSample = info.rangeFirstSample + ((info.numberRangeLooks - 1)/2) * info.slantRangePixelSpacing topo.demInterpolationMethod='BIQUINTIC' @@ -331,9 +331,10 @@ def runTopoCPU(info, demImage, dop=None, topo.topo() return + def runSimamp(outdir, hname='z.rdr'): from iscesys.StdOEL.StdOELPy import create_writer - + #####Run simamp stdWriter = create_writer("log","",True,filename='sim.log') objShade = isceobj.createSimamplitude() @@ -461,9 +462,9 @@ def extractInfo(frame, inps): info.prf = frame.PRF info.radarWavelength = frame.radarWavelegth info.orbit = frame.getOrbit() - - info.width = frame.getNumberOfSamples() - info.length = frame.getNumberOfLines() + + info.width = frame.getNumberOfSamples() + info.length = frame.getNumberOfLines() info.sensingStop = frame.getSensingStop() info.outdir = inps.outdir @@ -472,7 +473,7 @@ def extractInfo(frame, inps): def main(iargs=None): - + inps = cmdLineParse(iargs) # see if the user compiled isce with GPU enabled @@ -502,11 +503,9 @@ def main(iargs=None): doppler = db['doppler'] except: doppler = frame._dopplerVsPixel - db.close() - ####Setup dem demImage = isceobj.createDemImage() demImage.load(inps.dem + '.xml') @@ -522,7 +521,6 @@ def main(iargs=None): info.incFilename = os.path.join(info.outdir, 'incLocal.rdr') info.maskFilename = os.path.join(info.outdir, 'shadowMask.rdr') - runTopo(info,demImage,dop=doppler,nativedop=inps.nativedop, legendre=inps.legendre) runSimamp(os.path.dirname(info.heightFilename),os.path.basename(info.heightFilename)) @@ -540,4 +538,3 @@ if __name__ == '__main__': Main driver. ''' main() - From 62871d9806e2ce5f11093771ced116eccf81ea57 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Fri, 22 Nov 2019 23:27:40 -0800 Subject: [PATCH 24/27] topsStack: convert README.txt to README.md + contrib/stack/topsStack: - convert README.txt to README.md - add AUX_CAL file download instruction - add note for DEM preparation from stripmapStack/README.md + contrib/stack: - merge References.txt and INSTALL.txt to README.md - remove the duplicated README_stripmap/topsStack.txt files. --- contrib/stack/INSTALL.txt | 17 --- contrib/stack/README.md | 34 +++++ contrib/stack/README_stripmapStack.txt | 64 --------- contrib/stack/README_topsStack.txt | 117 ---------------- contrib/stack/References.txt | 11 -- .../stack/topsStack/{README.txt => README.md} | 129 ++++++++++++------ 6 files changed, 125 insertions(+), 247 deletions(-) delete mode 100644 contrib/stack/INSTALL.txt create mode 100644 contrib/stack/README.md delete mode 100644 contrib/stack/README_stripmapStack.txt delete mode 100644 contrib/stack/README_topsStack.txt delete mode 100644 contrib/stack/References.txt rename contrib/stack/topsStack/{README.txt => README.md} (59%) diff --git a/contrib/stack/INSTALL.txt b/contrib/stack/INSTALL.txt deleted file mode 100644 index 0f94988..0000000 --- a/contrib/stack/INSTALL.txt +++ /dev/null @@ -1,17 +0,0 @@ -To use the TOPS or Stripmap stack processors you need to: - -1- Install ISCE as usual -2- Depending on which stack processor you need to try, add the path of the folder containing the python scripts to your $PATH environment variable as follows: - - to use the topsStack for processing a stack of Sentinel-1 tops data add the full path of your "contrib/stack/topsStack" to your $PATH environemnt variable - - to use the stripmapStack for processing a stack of Stripmap data, add the full path of your "contrib/stack/stripmapStack" to your $PATH environemnt variableu - -NOTE: - - The stack processors do not show up in the install directory of your isce software. They can be found in the isce source directory. - -Important Note: - - There might be conflicts between topsStack and stripmapStack scripts (due to comman names of different scripts). Therefore users MUST only have the path of one stack processor in their $PATH environment at a time, to avoid conflicts between the two stack processors. - - - diff --git a/contrib/stack/README.md b/contrib/stack/README.md new file mode 100644 index 0000000..389ddf3 --- /dev/null +++ b/contrib/stack/README.md @@ -0,0 +1,34 @@ +## Stack Processors + +Read the document for each stack processor for details. + ++ [stripmapStack](./stripmapStack/README.md) ++ [topsStack](./topsStack/README.md) + +### Installation + +To use the TOPS or Stripmap stack processors you need to: + +1. Install ISCE as usual + +2. Depending on which stack processor you need to try, add the path of the folder containing the python scripts to your `$PATH` environment variable as follows: + - add the full path of your **contrib/stack/topsStack** to `$PATH` to use the topsStack for processing a stack of Sentinel-1 TOPS data + - add the full path of your **contrib/stack/stripmapStack** to `$PATH` to use the stripmapStack for processing a stack of StripMap data + +Note: The stack processors do not show up in the install directory of your isce software. They can be found in the isce source directory. + +#### Important Note: #### + +There might be conflicts between topsStack and stripmapStack scripts (due to comman names of different scripts). Therefore users **MUST only** have the path of **one stack processor in their $PATH environment at a time**, to avoid conflicts between the two stack processors. + +### References + +Users who use the stack processors may refer to the following literatures: + +For StripMap stack processor and ionospheric phase estimation: + ++ H. Fattahi, M. Simons, and P. Agram, "InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique", IEEE Trans. Geosci. Remote Sens., vol. 55, no. 10, 5984-5996, 2017. (https://ieeexplore.ieee.org/abstract/document/7987747/) + +For TOPS stack processing: + ++ H. Fattahi, P. Agram, and M. Simons, “A network-based enhanced spectral diversity approach for TOPS time-series analysis,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 2, pp. 777–786, Feb. 2017. (https://ieeexplore.ieee.org/abstract/document/7637021/) diff --git a/contrib/stack/README_stripmapStack.txt b/contrib/stack/README_stripmapStack.txt deleted file mode 100644 index e6ae915..0000000 --- a/contrib/stack/README_stripmapStack.txt +++ /dev/null @@ -1,64 +0,0 @@ - -The detailed algorithms for stack processing of stripmap data can be found here: - -H. Fattahi, M. Simons, and P. Agram, "InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique", IEEE Trans. Geosci. Remote Sens., vol. 55, no. 10, 5984-5996, 2017. (https://ieeexplore.ieee.org/abstract/document/7987747/) - - ------------------------------------ - -Notes on stripmap stack processor: - -Here are some notes to get started with processing stacks of stripmap data with ISCE. - - -1- create a folder somewhere for your project - -mkdir MauleT111 -cd MauleT111 - -2- create a DEM: - -dem.py -a stitch -b -37 -31 -72 -69 -r -s 1 -c - -3- Keep only ".dem.wgs84", ".dem.wgs84.vrt" and ".dem.wgs84.xml" and remove unnecessary files - -4- fix the path of the file in the xml file of the DEM by using this command: - -fixImageXml.py -f -i demLat_S37_S31_Lon_W072_W069.dem.wgs84 - -5- create a folder to download the ALOS-1 data from ASF: - -mkdir download -cd download - -6- Download the data that that you want to process to the downlowd directory. - -7- once all data have been downloaded, we need to unzip them and move them to different folders and getting ready for unpacking and then SLC generation. -This can be done by running the following command in a directory above "download": - -prepRawALOS.py -i download/ -o SLC - -This command generates an empty SLC folder and a run file called: "run_unPackALOS". -You could also run prepRawSensor.py which aims to recognize the sensor data automatically followed by running the sensor specific preparation script. For now we include support for ALOS and CSK raw data, but it is trivial to expand and include other sensors as unpacking routines are already included in the distribution. - -prepRawSensor.py -i download/ -o SLC - -8- execute the commands inside run_unPackALOS file. If you have a cluster that you can submit jobs, you can submit each line of command to a processor. The commands are independent and can be run in parallel. - -9- After successfully running the previous step, you should see acquisition dates in the SLC folder and the ".raw" files for each acquisition - -Note: For ALOS-1, If there is an acquisition that does not include .raw file, this is most likely due to PRF change between frames and can not be currently handled by ISCE. You have to ignore those. - -10- run stackStripmap.py which will generate many config and run files that need to be executed. Here is an example: - -stackStripMap.py -s SLC/ -d demLat_S37_S31_Lon_W072_W069.dem.wgs84 -t 250 -b 1000 -a 14 -r 4 -u snaphu - -This will produce: -a) baseline folder, which contains baseline information -b) pairs.png which is a baseline-time plot of the network of interferograms -c) configs: which contains the configuration parameter to run different InSAR processing steps -d) run_files: a folder that includes several run and job files that needs to be run in order - -11- execute the commands in run files (run_1, run_2, etc) in the run_files folder - - diff --git a/contrib/stack/README_topsStack.txt b/contrib/stack/README_topsStack.txt deleted file mode 100644 index 262c78e..0000000 --- a/contrib/stack/README_topsStack.txt +++ /dev/null @@ -1,117 +0,0 @@ -The detailed algorithm for stack processing of TOPS data can be find here: - -H. Fattahi, P. Agram, and M. Simons, “A network-based enhanced spectral diversity approach for TOPS time-series analysis,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 2, pp. 777–786, Feb. 2017. (https://ieeexplore.ieee.org/abstract/document/7637021/) - -<<<<<< Sentinel-1 TOPS stack processor >>>>>> - - -To use the sentinel stack processor, make sure to add the path of your "contrib/stack/topsStack" folder to your $PATH environment varibale. - - -The scripts provides support for Sentinel-1 TOPS stack processing. Currently supported workflows include a coregistered stack of SLC, interferograms, offsets, and coherence. - -stackSentinel.py generates all configuration and run files required to be executed on a stack of Sentinel-1 TOPS data. When stackSentinel.py is executed for a given workflow (-W option) a “configs” and “run_files” folder is generated. No processing is performed at this stage. Within the “run_files” folder different run_#_description files are contained which are to be executed as shell scripts in the run number order. Each of these run scripts call specific configure files contained in the “configs” folder which call ISCE in a modular fashion. The configure and run files will change depending on the selected workflow. To make run_# files executable, change the file permission accordingly (e.g., chmod +x run_1_unpack_slc). - -To see workflow examples, type “stackSentinel.py -H” -To get an overview of all the configurable parameters, type “stackSentinel.py -h” - -Required parameters of stackSentinel.py include: --s SLC_DIRNAME A folder with downloaded Sentinel-1 SLC’s. --o ORBIT_DIRNAME A folder containing the Sentinel-1 orbits. -Missing orbit files will be downloaded automatically --a AUX_DIRNAME A folder containing the Sentinel-1 Auxiliary files --d DEM A DEM (Digital Elevation Model) referenced to wgs84 - - -In the following, different workflow examples are provided. Note that stackSentinel.py only generates the run and configure files. To perform the actual processing, the user will need to execute each run file in their numbered order. - -In all workflows, coregistration (-C option) can be done using only geometry (set option = geometry) or with geometry plus refined azimuth offsets through NESD (set option = NESD) approach, the latter being the default. For the NESD coregistrstion the user can control the ESD coherence threshold (-e option) and the number of overlap interferograms (-O) to be used in NESD estimation. - ------------------------------- Example 1: Coregistered stack of SLC ---------------------------- -Generate the run and configure files needed to generate a coregistered stack of SLCs. -In this example, a pre-defined bounding box is specified. Note, if the bounding box is not provided it is set by default to the common SLC area among all SLCs. We recommend that user always set the processing bounding box. Since ESA does not have a fixed frame definition, we suggest to download data for a larger bounding box compared to the actual bounding box used in stackSentinel.py. This way user can ensure to have required data to cover the region of interest. Here is an example command to create configuration files for a stack of SLCs: - -stackSentinel.py -s ../SLC/ -d ../DEM/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -a ../../AuxDir/ -o ../../Orbits -b '19 20 -99.5 -98.5' -W slc -by running the command above, the configs and run_files folders are created. User needs to execute each run file in order. The order is specified by the index number of the run file name. For the example above, the run_files folder includes the following files: -- run_1_unpack_slc_topo_master -- run_2_average_baseline -- run_3_extract_burst_overlaps -- run_4_overlap_geo2rdr_resample -- run_5_pairs_misreg -- run_6_timeseries_misreg -- run_7_geo2rdr_resample -- run_8_extract_stack_valid_region -- run_9_merge -- run_10_grid_baseline - -The generated run files are self descriptive. Below is a short explanation on what each run_file does: - -***run_1_unpack_slc_topo_master:*** -Includes commands to unpack Sentinel-1 TOPS SLCs using ISCE readers. For older SLCs which need antenna elevation pattern correction, the file is extracted and written to disk. For newer version of SLCs which don’t need the elevation antenna pattern correction, only a gdal virtual “vrt” file (and isce xml file) is generated. The “.vrt” file points to the Sentinel SLC file and reads them whenever required during the processing. If a user wants to write the “.vrt” SLC file to disk, it can be done easily using gdal_translate (e.g. gdal_translate –of ENVI File.vrt File.slc). -The “run_1_unpack_slc_topo_master” also includes a command that refers to the config file of the stack master, which includes configuration for running topo for the stack master. Note that in the pair-wise processing strategy one should run topo (mapping from range-Doppler to geo coordinate) for all pairs. However, with stackSentinel, topo needs to be run only one time for the master in the stack. - -***run_2_average_baseline: *** -Computes average baseline for the stack. These baselines are not used for processing anywhere. They are only an approximation and can be used for plotting purposes. A more precise baseline grid is estimated later in run_10. - -***run_3_extract_burst_overlaps: *** -Burst overlaps are extracted for estimating azimuth misregistration using NESD technique. If coregistration method is chosen to be “geometry”, then this run file won’t exist and the overlaps are not extracted. - -***run_4_overlap_geo2rdr_resample: *** -Running geo2rdr to estimate geometrical offsets between slave burst overlaps and the stack master burst overlaps. The slave burst overlaps are then resampled to the stack master burst overlaps. - -***run_5_pairs_misreg: *** -Using the coregistered stack burst overlaps generated from the previous step, differential overlap interferograms are generated and are used for estimating azimuth misregistration using Enhanced Spectral Diversity (ESD) technique. - -***run_6_timeseries_misreg: *** -A time-series of azimuth and range misregistration is estimated with respect to the stack master. The time-series is a least squares esatimation from the pair misregistration from the previous step. - -***run_7_geo2rdr_resample: *** -Using orbit and DEM, geometrical offsets among all slave SLCs and the stack master is computed. The goometrical offsets, together with the misregistration time-series (from previous step) are used for precise coregistration of each burst SLC. - -***run_8_extract_stack_valid_region: *** -The valid region between burst SLCs at the overlap area of the bursts slightly changes for different acquisitions. Therefore we need to keep track of these overlaps which will be used during merging bursts. Without these knowledge, lines of invalid data may appear in the merged products at the burst overlaps. - -***run_9_merge: *** -Merges all bursts for the master and coregistered SLCs. The geometry files are also merged including longitude, latitude, shadow and layer mask, line-of-sight files, etc. . - -***run_10_grid_baseline: *** -A coarse grid of baselines between each slave SLC and the stack master is generated. This is not used in any computation. - - --------- Example 2: Coregistered stack of SLC with modified parameters ----------- -In the following example, the same stack generation is requested but where the threshold of the default coregistration approach (NESD) is relaxed from its default 0.85 value 0.7. - -stackSentinel.py -s ../SLC/ -d ../DEM/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -a ../../AuxDir/ -o ../../Orbits -b '19 20 -99.5 -98.5' -W slc -e 0.7 - -When running all the run files, the final products are located in the merge folder which has subdirectories “geom_master”, “baselines” and “SLC”. The “geom_master” folder contains geometry products such as longitude, latitude, height, local incidence angle, look angle, heading, and shadowing/layover mask files. The “baselines” folder contains sparse grids of the perpendicular baseline for each acquisition, while the “SLC” folder contains the coregistered SLCs - - ------------------------------- Example 3: Stack of interferograms ------------------------------ -Generate the run and configure files needed to generate a stack of interferograms. -In this example, a stack of interferograms is requested for which up to 2 nearest neighbor connections are included. - -stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -b '19 20 -99.5 -98.5' -a ../../AuxDir/ -o ../../Orbits -c 2 - -In the following example, all possible interferograms are being generated and in which the coregistration approach is set to use geometry and not the default NESD. - -stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -b '19 20 -99.5 -98.5' -a ../../AuxDir/ -o ../../Orbits -C geometry -c all - -When executing all the run files, a coregistered stack of slcs are produced, the burst interferograms are generated and then merged. Merged interferograms are multilooked, filtered and unwrapped. Geocoding is not applied. If users need to geocode any product, they can use the geocodeGdal.py script. - - --------------------- Example 4: Correlation stack example ---------------------------- -Generate the run and configure files needed to generate a stack of coherence. -In this example, a correlation stack is requested considering all possible coherence pairs and where the coregistration approach is done using geometry only. - -stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -b '19 20 -99.5 -98.5' -a ../../AuxDir/ -o ../../Orbits -C geometry -c all -W correlation - -This workflow is basically similar to the previous one. The difference is that the interferograms are not unwrapped. - - ------------------------------------ DEM download example ----------------------------------- -Download of DEM (need to use wgs84 version) using the ISCE DEM download script. -dem.py -a stitch -b 18 20 -100 -97 -r -s 1 –c - -Updating DEM’s wgs84 xml to include full path to the DEM -fixImageXml.py -f -i demLat_N18_N20_Lon_W100_W097.dem.wgs84 - diff --git a/contrib/stack/References.txt b/contrib/stack/References.txt deleted file mode 100644 index 80e87a4..0000000 --- a/contrib/stack/References.txt +++ /dev/null @@ -1,11 +0,0 @@ -Users who use the stack processors may refer to the following literatures: - -for stripmap stack processor and ionospheric phase estimation: - -H. Fattahi, M. Simons, and P. Agram, "InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique", IEEE Trans. Geosci. Remote Sens., vol. 55, no. 10, 5984-5996, 2017. (https://ieeexplore.ieee.org/abstract/document/7987747/) - -For TOPS stack processing: - -H. Fattahi, P. Agram, and M. Simons, “A network-based enhanced spectral diversity approach for TOPS time-series analysis,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 2, pp. 777–786, Feb. 2017. (https://ieeexplore.ieee.org/abstract/document/7637021/) - - diff --git a/contrib/stack/topsStack/README.txt b/contrib/stack/topsStack/README.md similarity index 59% rename from contrib/stack/topsStack/README.txt rename to contrib/stack/topsStack/README.md index 62ad5f4..67b6ecc 100644 --- a/contrib/stack/topsStack/README.txt +++ b/contrib/stack/topsStack/README.md @@ -1,38 +1,80 @@ +## Sentinel-1 TOPS stack processor + The detailed algorithm for stack processing of TOPS data can be find here: -H. Fattahi, P. Agram, and M. Simons, “A network-based enhanced spectral diversity approach for TOPS time-series analysis,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 2, pp. 777–786, Feb. 2017. (https://ieeexplore.ieee.org/abstract/document/7637021/) ++ Fattahi, H., P. Agram, and M. Simons (2016), A Network-Based Enhanced Spectral Diversity Approach for TOPS Time-Series Analysis, IEEE Transactions on Geoscience and Remote Sensing, 55(2), 777-786, doi:[10.1109/TGRS.2016.2614925](https://ieeexplore.ieee.org/abstract/document/7637021). +----------------------------------- -<<<<<< Sentinel-1 TOPS stack processor >>>>>> - -To use the sentinel stack processor, make sure to add the path of your "contrib/stack/topsStack" folder to your $PATH environment varibale. - +To use the sentinel stack processor, make sure to add the path of your `contrib/stack/topsStack` folder to your `$PATH` environment varibale. The scripts provides support for Sentinel-1 TOPS stack processing. Currently supported workflows include a coregistered stack of SLC, interferograms, offsets, and coherence. -stackSentinel.py generates all configuration and run files required to be executed on a stack of Sentinel-1 TOPS data. When stackSentinel.py is executed for a given workflow (-W option) a “configs” and “run_files” folder is generated. No processing is performed at this stage. Within the “run_files” folder different run_#_description files are contained which are to be executed as shell scripts in the run number order. Each of these run scripts call specific configure files contained in the “configs” folder which call ISCE in a modular fashion. The configure and run files will change depending on the selected workflow. To make run_# files executable, change the file permission accordingly (e.g., chmod +x run_1_unpack_slc). +`stackSentinel.py` generates all configuration and run files required to be executed on a stack of Sentinel-1 TOPS data. When stackSentinel.py is executed for a given workflow (-W option) a **configs** and **run_files** folder is generated. No processing is performed at this stage. Within the run_files folder different run\_#\_description files are contained which are to be executed as shell scripts in the run number order. Each of these run scripts call specific configure files contained in the “configs” folder which call ISCE in a modular fashion. The configure and run files will change depending on the selected workflow. To make run_# files executable, change the file permission accordingly (e.g., `chmod +x run_1_unpack_slc`). -To see workflow examples, type “stackSentinel.py -H” -To get an overview of all the configurable parameters, type “stackSentinel.py -h” +```bash +stackSentinel.py -H #To see workflow examples, +stackSentinel.py -h #To get an overview of all the configurable parameters +``` Required parameters of stackSentinel.py include: --s SLC_DIRNAME A folder with downloaded Sentinel-1 SLC’s. --o ORBIT_DIRNAME A folder containing the Sentinel-1 orbits. -Missing orbit files will be downloaded automatically --a AUX_DIRNAME A folder containing the Sentinel-1 Auxiliary files --d DEM A DEM (Digital Elevation Model) referenced to wgs84 +```cfg +-s SLC_DIRNAME #A folder with downloaded Sentinel-1 SLC’s. +-o ORBIT_DIRNAME #A folder containing the Sentinel-1 orbits. Missing orbit files will be downloaded automatically +-a AUX_DIRNAME #A folder containing the Sentinel-1 Auxiliary files +-d DEM_FILENAME #A DEM (Digital Elevation Model) referenced to wgs84 +``` In the following, different workflow examples are provided. Note that stackSentinel.py only generates the run and configure files. To perform the actual processing, the user will need to execute each run file in their numbered order. In all workflows, coregistration (-C option) can be done using only geometry (set option = geometry) or with geometry plus refined azimuth offsets through NESD (set option = NESD) approach, the latter being the default. For the NESD coregistrstion the user can control the ESD coherence threshold (-e option) and the number of overlap interferograms (-O) to be used in NESD estimation. ------------------------------- Example 1: Coregistered stack of SLC ---------------------------- -Generate the run and configure files needed to generate a coregistered stack of SLCs. -In this example, a pre-defined bounding box is specified. Note, if the bounding box is not provided it is set by default to the common SLC area among all SLCs. We recommend that user always set the processing bounding box. Since ESA does not have a fixed frame definition, we suggest to download data for a larger bounding box compared to the actual bounding box used in stackSentinel.py. This way user can ensure to have required data to cover the region of interest. Here is an example command to create configuration files for a stack of SLCs: +#### AUX_CAL file download #### +The following calibration auxliary (AUX_CAL) file is used for **antenna pattern correction** to compensate the range phase offset of SAFE products with **IPF verison 002.36** (mainly for images acquired before March 2015). If all your SAFE products are from another IPF version, then no AUX files are needed. Check [ESA document](https://earth.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction) for details. + +Run the command below to download the AUX_CAL file once and store it somewhere (_i.e._ ~/aux/aux_cal) so that you can use it all the time, for `stackSentinel.py -a` or `auxiliary data directory` in `topsApp.py`. + +``` +wget https://qc.sentinel1.eo.esa.int/product/S1A/AUX_CAL/20140908T000000/S1A_AUX_CAL_V20140908T000000_G20190626T100201.SAFE.TGZ +tar zxvf S1A_AUX_CAL_V20140908T000000_G20190626T100201.SAFE.TGZ +rm S1A_AUX_CAL_V20140908T000000_G20190626T100201.SAFE.TGZ +``` + +#### 1. Create your project folder somewhere #### + +``` +mkdir MexicoSenAT72 +cd MexicoSenAT72 +``` + +#### 2. Prepare DEM #### + +Download of DEM (need to use wgs84 version) using the ISCE DEM download script. + +``` +mkdir DEM; cd DEM +dem.py -a stitch -b 18 20 -100 -97 -r -s 1 –c +rm demLat*.dem demLat*.dem.xml demLat*.dem.vrt +fixImageXml.py -f -i demLat*.dem.wgs84 #Updating DEM’s wgs84 xml to include full path to the DEM +cd .. +``` + +#### 3. Download Sentinel-1 data to SLC #### + + + +#### 4.1 Example workflow: Coregistered stack of SLC #### + +Generate the run and configure files needed to generate a coregistered stack of SLCs. In this example, a pre-defined bounding box is specified. Note, if the bounding box is not provided it is set by default to the common SLC area among all SLCs. We recommend that user always set the processing bounding box. Since ESA does not have a fixed frame definition, we suggest to download data for a larger bounding box compared to the actual bounding box used in stackSentinel.py. This way user can ensure to have required data to cover the region of interest. Here is an example command to create configuration files for a stack of SLCs: + +``` stackSentinel.py -s ../SLC/ -d ../DEM/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -a ../../AuxDir/ -o ../../Orbits -b '19 20 -99.5 -98.5' -W slc +``` + by running the command above, the configs and run_files folders are created. User needs to execute each run file in order. The order is specified by the index number of the run file name. For the example above, the run_files folder includes the following files: + - run_1_unpack_slc_topo_master - run_2_average_baseline - run_3_extract_burst_overlaps @@ -46,72 +88,83 @@ by running the command above, the configs and run_files folders are created. Use The generated run files are self descriptive. Below is a short explanation on what each run_file does: -***run_1_unpack_slc_topo_master:*** +**run_1_unpack_slc_topo_master:** + Includes commands to unpack Sentinel-1 TOPS SLCs using ISCE readers. For older SLCs which need antenna elevation pattern correction, the file is extracted and written to disk. For newer version of SLCs which don’t need the elevation antenna pattern correction, only a gdal virtual “vrt” file (and isce xml file) is generated. The “.vrt” file points to the Sentinel SLC file and reads them whenever required during the processing. If a user wants to write the “.vrt” SLC file to disk, it can be done easily using gdal_translate (e.g. gdal_translate –of ENVI File.vrt File.slc). The “run_1_unpack_slc_topo_master” also includes a command that refers to the config file of the stack master, which includes configuration for running topo for the stack master. Note that in the pair-wise processing strategy one should run topo (mapping from range-Doppler to geo coordinate) for all pairs. However, with stackSentinel, topo needs to be run only one time for the master in the stack. -***run_2_average_baseline: *** +**run_2_average_baseline:** + Computes average baseline for the stack. These baselines are not used for processing anywhere. They are only an approximation and can be used for plotting purposes. A more precise baseline grid is estimated later in run_10. -***run_3_extract_burst_overlaps: *** +**run_3_extract_burst_overlaps:** + Burst overlaps are extracted for estimating azimuth misregistration using NESD technique. If coregistration method is chosen to be “geometry”, then this run file won’t exist and the overlaps are not extracted. -***run_4_overlap_geo2rdr_resample: *** +**run_4_overlap_geo2rdr_resample:*** + Running geo2rdr to estimate geometrical offsets between slave burst overlaps and the stack master burst overlaps. The slave burst overlaps are then resampled to the stack master burst overlaps. -***run_5_pairs_misreg: *** +**run_5_pairs_misreg:** + Using the coregistered stack burst overlaps generated from the previous step, differential overlap interferograms are generated and are used for estimating azimuth misregistration using Enhanced Spectral Diversity (ESD) technique. -***run_6_timeseries_misreg: *** +**run_6_timeseries_misreg:** + A time-series of azimuth and range misregistration is estimated with respect to the stack master. The time-series is a least squares esatimation from the pair misregistration from the previous step. -***run_7_geo2rdr_resample: *** +**run_7_geo2rdr_resample:** + Using orbit and DEM, geometrical offsets among all slave SLCs and the stack master is computed. The goometrical offsets, together with the misregistration time-series (from previous step) are used for precise coregistration of each burst SLC. -***run_8_extract_stack_valid_region: *** +**run_8_extract_stack_valid_region:** + The valid region between burst SLCs at the overlap area of the bursts slightly changes for different acquisitions. Therefore we need to keep track of these overlaps which will be used during merging bursts. Without these knowledge, lines of invalid data may appear in the merged products at the burst overlaps. -***run_9_merge: *** +**run_9_merge:** + Merges all bursts for the master and coregistered SLCs. The geometry files are also merged including longitude, latitude, shadow and layer mask, line-of-sight files, etc. . -***run_10_grid_baseline: *** +**run_10_grid_baseline:** + A coarse grid of baselines between each slave SLC and the stack master is generated. This is not used in any computation. +#### 4.2 Example workflow: Coregistered stack of SLC with modified parameters #### --------- Example 2: Coregistered stack of SLC with modified parameters ----------- In the following example, the same stack generation is requested but where the threshold of the default coregistration approach (NESD) is relaxed from its default 0.85 value 0.7. +``` stackSentinel.py -s ../SLC/ -d ../DEM/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -a ../../AuxDir/ -o ../../Orbits -b '19 20 -99.5 -98.5' -W slc -e 0.7 +``` -When running all the run files, the final products are located in the merge folder which has subdirectories “geom_master”, “baselines” and “SLC”. The “geom_master” folder contains geometry products such as longitude, latitude, height, local incidence angle, look angle, heading, and shadowing/layover mask files. The “baselines” folder contains sparse grids of the perpendicular baseline for each acquisition, while the “SLC” folder contains the coregistered SLCs +When running all the run files, the final products are located in the merge folder which has subdirectories **geom_master**, **baselines** and **SLC**. The **geom_master** folder contains geometry products such as longitude, latitude, height, local incidence angle, look angle, heading, and shadowing/layover mask files. The **baselines** folder contains sparse grids of the perpendicular baseline for each acquisition, while the **SLC** folder contains the coregistered SLCs +#### 4.3 Example workflow: Stack of interferograms #### ------------------------------- Example 3: Stack of interferograms ------------------------------ Generate the run and configure files needed to generate a stack of interferograms. In this example, a stack of interferograms is requested for which up to 2 nearest neighbor connections are included. +``` stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -b '19 20 -99.5 -98.5' -a ../../AuxDir/ -o ../../Orbits -c 2 +``` In the following example, all possible interferograms are being generated and in which the coregistration approach is set to use geometry and not the default NESD. +``` stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -b '19 20 -99.5 -98.5' -a ../../AuxDir/ -o ../../Orbits -C geometry -c all +``` When executing all the run files, a coregistered stack of slcs are produced, the burst interferograms are generated and then merged. Merged interferograms are multilooked, filtered and unwrapped. Geocoding is not applied. If users need to geocode any product, they can use the geocodeGdal.py script. +#### 4.4 Example workflow: Stack of correlation #### --------------------- Example 4: Correlation stack example ---------------------------- Generate the run and configure files needed to generate a stack of coherence. In this example, a correlation stack is requested considering all possible coherence pairs and where the coregistration approach is done using geometry only. +``` stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem.wgs84 -b '19 20 -99.5 -98.5' -a ../../AuxDir/ -o ../../Orbits -C geometry -c all -W correlation +``` This workflow is basically similar to the previous one. The difference is that the interferograms are not unwrapped. - ------------------------------------ DEM download example ----------------------------------- -Download of DEM (need to use wgs84 version) using the ISCE DEM download script. -dem.py -a stitch -b 18 20 -100 -97 -r -s 1 –c - -Updating DEM’s wgs84 xml to include full path to the DEM -fixImageXml.py -f -i demLat_N18_N20_Lon_W100_W097.dem.wgs84 - +#### 5. Execute the commands in run files (run_1*, run_2*, etc) in the "run_files" folder #### From 35ca7db85e34639bc4faf0dd925322fe535e0fc7 Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Wed, 22 Jan 2020 10:07:18 -0800 Subject: [PATCH 25/27] Fix trackCRs script * Strip leading zero from LLH list * Cast slice bounds to ints * Use comma delimiter for reading CSVs --- contrib/stack/stripmapStack/trackCRs.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/contrib/stack/stripmapStack/trackCRs.py b/contrib/stack/stripmapStack/trackCRs.py index e841d4a..b98d1bb 100755 --- a/contrib/stack/stripmapStack/trackCRs.py +++ b/contrib/stack/stripmapStack/trackCRs.py @@ -73,8 +73,7 @@ def makeOnePlot(filename, pos): minx = np.clip(np.min(pos[:,2])-win, 0, npix-1) maxx = np.clip(np.max(pos[:,2])+win, 0, npix-1) - - box = np.power(np.abs(data[miny:maxy, minx:maxx]), 0.4) + box = np.power(np.abs(data[int(miny):int(maxy), int(minx):int(maxx)]), 0.4) plt.figure('CR analysis') @@ -104,7 +103,7 @@ def getAzRg(frame,llh): pol._normRange = frame.instrument.rangePixelSize pol.initPoly(azimuthOrder=0, rangeOrder=len(coeffs)-1, coeffs=[coeffs]) - taz, rgm = frame.orbit.geo2rdr(list(llh), side=frame.instrument.platform.pointingDirection, + taz, rgm = frame.orbit.geo2rdr(list(llh)[1:], side=frame.instrument.platform.pointingDirection, doppler=pol, wvl=frame.instrument.getRadarWavelength()) line = (taz - frame.sensingStart).total_seconds() * frame.PRF @@ -145,7 +144,7 @@ if __name__ == '__main__': # frame.startingRange = frame.startingRange + 100.0 ###Load CRS positions - llhs = np.loadtxt(inps.posfile) + llhs = np.loadtxt(inps.posfile, delimiter=',') crs = [] From 5d6a731753d2976739f92d6c2613cf22017d51f1 Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Thu, 30 Jan 2020 22:26:46 -0800 Subject: [PATCH 26/27] Remove internal usage of ISCE_HOME env variable ISCE_HOME was only used to get the location of the default logging config. Lots of scripts were using boilerplate to set up this config, so I added an `isce.logging` helper module which is the same as builtin python logging but already has the configuration defaults set up for isce. ISCE_HOME setup is retained in the toplevel `__init__.py` but can now be removed without affecting functionality. --- .circleci/config.yml | 2 +- __init__.py | 7 ++++++- applications/CalculatePegPoint.py | 6 +----- applications/calculateBaseline.py | 6 +----- applications/createGeneric.py | 6 +----- applications/extractHDROrbit.py | 6 +----- applications/focus.py | 6 +----- applications/insarApp.py | 8 +------- applications/isceApp.py | 7 +------ applications/make_raw.py | 10 +--------- applications/rtcApp.py | 9 +-------- applications/scansarApp.py | 11 +---------- applications/stripmapApp.py | 9 +-------- applications/topsApp.py | 9 +-------- applications/topsOffsetApp.py | 9 +-------- applications/viewMetadata.py | 9 +-------- components/isceobj/Location/Offset.py | 6 +----- components/isceobj/Scene/test/testTrack.py | 9 +-------- components/isceobj/Util/geo/__init__.py | 4 ---- components/iscesys/Component/Configurable.py | 5 +---- components/iscesys/DataRetriever/DataRetriever.py | 7 +------ components/stdproc/orbit/orbitLib/CalcSchHeightVel.py | 5 ----- contrib/demUtils/demstitcher/DemStitcher.py | 7 +------ contrib/demUtils/demstitcher/DemStitcherV3.py | 6 +----- contrib/demUtils/swbdstitcher/SWBDStitcher.py | 6 +----- contrib/demUtils/watermask/WaterMask.py | 7 +------ contrib/issi/applications/ISSI.py | 5 +---- 27 files changed, 30 insertions(+), 157 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 98bb5bd..25ba14b 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -72,7 +72,7 @@ jobs: set -ex pwd . /opt/conda/bin/activate root - export ISCE_HOME=/root/project/install/isce + ISCE_HOME=/root/project/install/isce export PATH="$ISCE_HOME/bin:$ISCE_HOME/applications:/opt/conda/bin:$PATH" export PYTHONPATH="/root/project/install:$PYTHONPATH" export LD_LIBRARY_PATH="/opt/conda/lib:$LD_LIBRARY_PATH" diff --git a/__init__.py b/__init__.py index a6e2c5b..f2f76a8 100755 --- a/__init__.py +++ b/__init__.py @@ -32,7 +32,12 @@ version = release_history # compatibility alias __version__ = release_version import sys, os -isce_path = os.path.split(os.path.abspath(__file__))[0] +isce_path = os.path.dirname(os.path.abspath(__file__)) + +import logging +from logging.config import fileConfig as _fc +_fc(os.path.join(isce_path, 'defaults', 'logging', 'logging.conf')) + sys.path.insert(1,isce_path) sys.path.insert(1,os.path.join(isce_path,'applications')) sys.path.insert(1,os.path.join(isce_path,'components')) diff --git a/applications/CalculatePegPoint.py b/applications/CalculatePegPoint.py index 5f52ca4..a300c0e 100755 --- a/applications/CalculatePegPoint.py +++ b/applications/CalculatePegPoint.py @@ -31,12 +31,8 @@ -import os import math -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from iscesys.Compatibility import Compatibility Compatibility.checkPythonVersion() from isceobj.Location.Peg import Peg diff --git a/applications/calculateBaseline.py b/applications/calculateBaseline.py index 375fec3..a6eeff5 100755 --- a/applications/calculateBaseline.py +++ b/applications/calculateBaseline.py @@ -30,11 +30,7 @@ -import os -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from iscesys.Compatibility import Compatibility Compatibility.checkPythonVersion() from iscesys.Component.FactoryInit import FactoryInit diff --git a/applications/createGeneric.py b/applications/createGeneric.py index aebb6ba..c7f706b 100755 --- a/applications/createGeneric.py +++ b/applications/createGeneric.py @@ -30,11 +30,7 @@ -import os -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging import isceobj from iscesys.Component.FactoryInit import FactoryInit diff --git a/applications/extractHDROrbit.py b/applications/extractHDROrbit.py index 421b60b..fd7278c 100755 --- a/applications/extractHDROrbit.py +++ b/applications/extractHDROrbit.py @@ -30,12 +30,8 @@ -import os import datetime -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from iscesys.Compatibility import Compatibility Compatibility.checkPythonVersion() from iscesys.Component.FactoryInit import FactoryInit diff --git a/applications/focus.py b/applications/focus.py index 280a9d2..53f64ed 100755 --- a/applications/focus.py +++ b/applications/focus.py @@ -30,12 +30,8 @@ -import os import math -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging import isceobj from iscesys.Component.FactoryInit import FactoryInit from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTU diff --git a/applications/insarApp.py b/applications/insarApp.py index 0ab2f0a..8136d1d 100755 --- a/applications/insarApp.py +++ b/applications/insarApp.py @@ -34,8 +34,7 @@ from __future__ import print_function import time import os import sys -import logging -import logging.config +from isce import logging import isce import isceobj @@ -46,11 +45,6 @@ from iscesys.Component.Configurable import SELF import isceobj.InsarProc as InsarProc from isceobj.Scene.Frame import FrameMixin -logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') -) - logger = logging.getLogger('isce.insar') diff --git a/applications/isceApp.py b/applications/isceApp.py index 7aec323..f67d08d 100755 --- a/applications/isceApp.py +++ b/applications/isceApp.py @@ -41,8 +41,7 @@ import datetime import os import sys import math -import logging -import logging.config +from isce import logging import isce import isceobj @@ -1439,10 +1438,6 @@ class IsceApp(Application, FrameMixin): os.chdir(self.outputDir) ##change working directory to given output directory ##read configfile only here so that log path is in output directory - logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') - ) logger = logging.getLogger('isce.isceProc') logger.info(self.intromsg) self._isce.dataDirectory = self.outputDir diff --git a/applications/make_raw.py b/applications/make_raw.py index a3b1f44..3c2af20 100755 --- a/applications/make_raw.py +++ b/applications/make_raw.py @@ -27,16 +27,8 @@ # Author: Walter Szeliga #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - - -import os -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) - import isce +from isce import logging from iscesys.Compatibility import Compatibility from iscesys.Component.Component import Component, Port from isceobj.Planet.Ellipsoid import Ellipsoid diff --git a/applications/rtcApp.py b/applications/rtcApp.py index 836151e..05d4575 100755 --- a/applications/rtcApp.py +++ b/applications/rtcApp.py @@ -30,10 +30,8 @@ import time -import os import sys -import logging -import logging.config +from isce import logging import isce import isceobj @@ -44,11 +42,6 @@ from iscesys.Component.Configurable import SELF from isceobj import RtcProc from isceobj.Util.decorators import use_api -logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') -) - logger = logging.getLogger('isce.grdsar') diff --git a/applications/scansarApp.py b/applications/scansarApp.py index 55fe196..a5338ce 100755 --- a/applications/scansarApp.py +++ b/applications/scansarApp.py @@ -27,13 +27,9 @@ # Authors: Giangi Sacco, Eric Gurrola #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - import time -import os import sys -import logging -import logging.config +from isce import logging import isce import isceobj @@ -43,11 +39,6 @@ from iscesys.Compatibility import Compatibility from iscesys.Component.Configurable import SELF from isceobj import ScansarProc -logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') -) - logger = logging.getLogger('isce.insar') diff --git a/applications/stripmapApp.py b/applications/stripmapApp.py index 6f76b66..3793ac6 100755 --- a/applications/stripmapApp.py +++ b/applications/stripmapApp.py @@ -35,10 +35,8 @@ from __future__ import print_function import time -import os import sys -import logging -import logging.config +from isce import logging import isce import isceobj @@ -50,11 +48,6 @@ import isceobj.StripmapProc as StripmapProc from isceobj.Scene.Frame import FrameMixin from isceobj.Util.decorators import use_api -logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') -) - logger = logging.getLogger('isce.insar') diff --git a/applications/topsApp.py b/applications/topsApp.py index 7bc2dc9..7f76a37 100755 --- a/applications/topsApp.py +++ b/applications/topsApp.py @@ -34,10 +34,8 @@ import time -import os import sys -import logging -import logging.config +from isce import logging import isce import isceobj @@ -47,11 +45,6 @@ from iscesys.Compatibility import Compatibility from iscesys.Component.Configurable import SELF from isceobj import TopsProc -logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') -) - logger = logging.getLogger('isce.insar') diff --git a/applications/topsOffsetApp.py b/applications/topsOffsetApp.py index f2049a8..52b1e3b 100755 --- a/applications/topsOffsetApp.py +++ b/applications/topsOffsetApp.py @@ -30,10 +30,8 @@ import time -import os import sys -import logging -import logging.config +from isce import logging import isce import isceobj @@ -42,11 +40,6 @@ from isce.applications.topsApp import TopsInSAR from iscesys.Component.Application import Application from isceobj.Util.decorators import use_api -logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', 'logging', - 'logging.conf') -) - logger = logging.getLogger('isce.insar') WINDOW_SIZE_WIDTH = Application.Parameter( diff --git a/applications/viewMetadata.py b/applications/viewMetadata.py index 3051193..8021ff0 100755 --- a/applications/viewMetadata.py +++ b/applications/viewMetadata.py @@ -27,14 +27,7 @@ # Author: Walter Szeliga #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - - -import os -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from iscesys.Compatibility import Compatibility Compatibility.checkPythonVersion() from iscesys.Component.FactoryInit import FactoryInit diff --git a/components/isceobj/Location/Offset.py b/components/isceobj/Location/Offset.py index 8bcff9b..96960f4 100755 --- a/components/isceobj/Location/Offset.py +++ b/components/isceobj/Location/Offset.py @@ -29,11 +29,7 @@ import math -import os -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from isceobj.Util.decorators import type_check, force, pickled, logged import numpy as np diff --git a/components/isceobj/Scene/test/testTrack.py b/components/isceobj/Scene/test/testTrack.py index b67cc5c..51883d1 100755 --- a/components/isceobj/Scene/test/testTrack.py +++ b/components/isceobj/Scene/test/testTrack.py @@ -27,14 +27,7 @@ # Author: Walter Szeliga #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - - -import os -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from isceobj.Sensor.ERS import ERS from isceobj.Scene.Track import Track logger = logging.getLogger("testTrack") diff --git a/components/isceobj/Util/geo/__init__.py b/components/isceobj/Util/geo/__init__.py index 85f13b9..0c7ffe3 100755 --- a/components/isceobj/Util/geo/__init__.py +++ b/components/isceobj/Util/geo/__init__.py @@ -45,10 +45,6 @@ ellipsoid oblate ellipsoid of revolution (e.g, WGS84) with all the See mainpage.txt for a complete dump of geo's philosophy-- otherwise, use the docstrings. """ -import os -isce_path = os.getenv("ISCE_HOME") ## \namespace geo Vector- and Affine-spaces, on Earth __all__ = ['euclid', 'coordinates', 'ellipsoid', 'charts', 'affine', 'motion'] - - diff --git a/components/iscesys/Component/Configurable.py b/components/iscesys/Component/Configurable.py index 5e9b9ff..3cff30c 100755 --- a/components/iscesys/Component/Configurable.py +++ b/components/iscesys/Component/Configurable.py @@ -32,10 +32,7 @@ from __future__ import print_function import os import sys import operator -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging from iscesys.DictUtils.DictUtils import DictUtils as DU from iscesys.Compatibility import Compatibility Compatibility.checkPythonVersion() diff --git a/components/iscesys/DataRetriever/DataRetriever.py b/components/iscesys/DataRetriever/DataRetriever.py index 4f1a4b2..1599ff8 100755 --- a/components/iscesys/DataRetriever/DataRetriever.py +++ b/components/iscesys/DataRetriever/DataRetriever.py @@ -37,8 +37,7 @@ import isce import zipfile import os import sys -import logging -import logging.config +from isce import logging from iscesys.Component.Component import Component import shutil from urllib import request @@ -325,8 +324,4 @@ class DataRetriever(Component): # logger not defined until baseclass is called if not self.logger: - logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf') - ) self.logger = logging.getLogger('isce.iscesys.DataRetriever') diff --git a/components/stdproc/orbit/orbitLib/CalcSchHeightVel.py b/components/stdproc/orbit/orbitLib/CalcSchHeightVel.py index 684cbeb..53da593 100755 --- a/components/stdproc/orbit/orbitLib/CalcSchHeightVel.py +++ b/components/stdproc/orbit/orbitLib/CalcSchHeightVel.py @@ -29,10 +29,8 @@ -import os import logging import math -import logging.config from iscesys.Compatibility import Compatibility @@ -40,9 +38,6 @@ from isceobj.Planet import Planet from isceobj import Constants as CN from iscesys.Component.Component import Component, Port -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) - RANGE_SAMPLING_RATE = Component.Parameter('rangeSamplingRate', public_name='range sampling rate', type=float, diff --git a/contrib/demUtils/demstitcher/DemStitcher.py b/contrib/demUtils/demstitcher/DemStitcher.py index ea3d51f..ba3e1e3 100755 --- a/contrib/demUtils/demstitcher/DemStitcher.py +++ b/contrib/demUtils/demstitcher/DemStitcher.py @@ -43,8 +43,7 @@ import os import sys import math import urllib.request, urllib.parse, urllib.error -import logging -import logging.config +from isce import logging from iscesys.Component.Component import Component import xml.etree.ElementTree as ET @@ -1013,10 +1012,6 @@ class DemStitcher(Component): # logger not defined until baseclass is called if not self.logger: - logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf') - ) self.logger = logging.getLogger('isce.contrib.demUtils.DemStitcher') url = property(getUrl,setUrl) diff --git a/contrib/demUtils/demstitcher/DemStitcherV3.py b/contrib/demUtils/demstitcher/DemStitcherV3.py index c2b486d..6c5f02c 100755 --- a/contrib/demUtils/demstitcher/DemStitcherV3.py +++ b/contrib/demUtils/demstitcher/DemStitcherV3.py @@ -39,8 +39,7 @@ from ctypes import cdll import os import sys import urllib.request, urllib.error, urllib.parse -import logging -import logging.config +from isce import logging from iscesys.Component.Component import Component from contrib.demUtils.DemStitcher import DemStitcher as DS #Parameters definitions @@ -291,7 +290,4 @@ class DemStitcher(DS): #it's /srtm/version2_1/SRTM(1,3) self._remove = ['.jpg','.xml'] if not self.logger: - logging.config.fileConfig( - os.environ['ISCE_HOME'] + '/library/applications/logging.conf' - ) self.logger = logging.getLogger('isce.contrib.demUtils.DemStitcherV3') diff --git a/contrib/demUtils/swbdstitcher/SWBDStitcher.py b/contrib/demUtils/swbdstitcher/SWBDStitcher.py index cca55fb..e299e9b 100755 --- a/contrib/demUtils/swbdstitcher/SWBDStitcher.py +++ b/contrib/demUtils/swbdstitcher/SWBDStitcher.py @@ -39,9 +39,8 @@ from ctypes import cdll import numpy as np import os import sys -import logging +from isce import logging import math -import logging.config import urllib.request, urllib.parse, urllib.error from iscesys.Component.Component import Component from contrib.demUtils.DemStitcher import DemStitcher @@ -315,9 +314,6 @@ class SWBDStitcher(DemStitcher): #it's /srtm/version2_1/SRTM(1,3) self._remove = ['.jpg','.xml'] if not self.logger: - logging.config.fileConfig( - os.environ['ISCE_HOME'] + '/library/applications/logging.conf' - ) self.logger = logging.getLogger('isce.contrib.demUtils.SWBDStitcher') self.parameter_list = self.parameter_list + super(DemStitcher,self).parameter_list diff --git a/contrib/demUtils/watermask/WaterMask.py b/contrib/demUtils/watermask/WaterMask.py index 766060b..57b5603 100755 --- a/contrib/demUtils/watermask/WaterMask.py +++ b/contrib/demUtils/watermask/WaterMask.py @@ -35,8 +35,7 @@ import sys import math from html.parser import HTMLParser import urllib.request, urllib.parse, urllib.error -import logging -import logging.config +from isce import logging from iscesys.Component.Component import Component import zipfile import os @@ -979,10 +978,6 @@ class MaskStitcher(Component): # logger not defined until baseclass is called if not self.logger: - logging.config.fileConfig( - os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf') - ) self.logger = logging.getLogger('isce.contrib.demUtils.MaskStitcher') utl = property(getUrl,setUrl) diff --git a/contrib/issi/applications/ISSI.py b/contrib/issi/applications/ISSI.py index 330905b..9146035 100755 --- a/contrib/issi/applications/ISSI.py +++ b/contrib/issi/applications/ISSI.py @@ -32,10 +32,7 @@ import os import math -import logging -import logging.config -logging.config.fileConfig(os.path.join(os.environ['ISCE_HOME'], 'defaults', - 'logging', 'logging.conf')) +from isce import logging import isce from iscesys.Component.FactoryInit import FactoryInit From dbeb7f512109025173e55deda747827797a5b27e Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Tue, 4 Feb 2020 16:03:25 -0800 Subject: [PATCH 27/27] Remove old comment --- applications/isceApp.py | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/isceApp.py b/applications/isceApp.py index f67d08d..4dd6c36 100755 --- a/applications/isceApp.py +++ b/applications/isceApp.py @@ -1437,7 +1437,6 @@ class IsceApp(Application, FrameMixin): sys.exit("Could not find the output directory: %s" % self.outputDir) os.chdir(self.outputDir) ##change working directory to given output directory - ##read configfile only here so that log path is in output directory logger = logging.getLogger('isce.isceProc') logger.info(self.intromsg) self._isce.dataDirectory = self.outputDir