Pycuampcor v2

contrib/PyCuAmpcor/CMakeLists.txt

@@ -8,12 +8,16 @@ endif()
 
 set(CMAKE_CUDA_STANDARD 11)
 set(CMAKE_CUDA_STANDARD_REQUIRED TRUE)
+set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --prec-div=true --prec-sqrt=true")
+
 
 pybind11_add_module(PyCuAmpcor
     src/PyCuAmpcor.cpp
     src/GDALImage.cpp
     src/SConscript
-    src/cuAmpcorChunk.cpp
+    src/cuAmpcorProcessor.cpp
+    src/cuAmpcorProcessorROIPAC.cpp
+    src/cuAmpcorProcessorGrIMP.cpp
     src/cuAmpcorController.cpp
     src/cuAmpcorParameter.cpp
     src/cuArrays.cpp

contrib/PyCuAmpcor/CMakeLists.txt.isce2

@@ -0,0 +1,50 @@
+# Early exit if prereqs not available
+if(NOT TARGET GDAL::GDAL
+OR NOT TARGET CUDA::cufft
+OR NOT pybind11_FOUND
+   )
+    return()
+endif()
+
+set(CMAKE_CUDA_STANDARD 11)
+set(CMAKE_CUDA_STANDARD_REQUIRED TRUE)
+set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --prec-div=true --prec-sqrt=true")
+
+
+pybind11_add_module(PyCuAmpcor
+    src/PyCuAmpcor.cpp
+    src/GDALImage.cpp
+    src/SConscript
+    src/cuAmpcorProcessor.cpp
+    src/cuAmpcorProcessorROIPAC.cpp
+    src/cuAmpcorProcessorGrIMP.cpp
+    src/cuAmpcorController.cpp
+    src/cuAmpcorParameter.cpp
+    src/cuArrays.cpp
+    src/cuArraysCopy.cu
+    src/cuArraysPadding.cu
+    src/cuCorrFrequency.cu
+    src/cuCorrNormalization.cu
+    src/cuCorrNormalizationSAT.cu
+    src/cuCorrNormalizer.cpp
+    src/cuCorrTimeDomain.cu
+    src/cuDeramp.cu
+    src/cuEstimateStats.cu
+    src/cuOffset.cu
+    src/cuOverSampler.cpp
+    src/cuSincOverSampler.cu
+    src/cudaError.cpp
+    src/cudaUtil.cpp
+    )
+target_include_directories(PyCuAmpcor PRIVATE
+    src
+    )
+target_link_libraries(PyCuAmpcor PRIVATE
+    CUDA::cufft
+    GDAL::GDAL
+    )
+
+InstallSameDir(
+    __init__.py
+    PyCuAmpcor
+    )

contrib/PyCuAmpcor/CMakeLists.txt.standalone

@@ -0,0 +1,96 @@
+# -*- cmake -*-
+# Example with conda:
+#   cmake .. -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX -DCMAKE_CUDA_ARCHITECTURES=native -DCMAKE_PREFIX_PATH=${CONDA_PREFIX} -DCMAKE_BUILD_TYPE=Release
+
+# minimum cmake version to support CUDA and Python
+cmake_minimum_required(VERSION 3.18)
+
+# project name
+project(PyCuAmpcor LANGUAGES CUDA CXX)
+
+# C++ Standard
+set(CMAKE_CXX_STANDARD 11)
+
+# CUDA requirement
+find_package(CUDAToolkit COMPONENTS cudart cufft REQUIRED)
+set(CMAKE_CUDA_STANDARD 11)
+set(CMAKE_CUDA_STANDARD_REQUIRED TRUE)
+
+if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
+    set(CMAKE_CUDA_ARCHITECTURES 52 60 70 75 80 86 90)
+endif()
+
+# python and pybind11 requirement
+find_package(Python 3.6 COMPONENTS Interpreter Development REQUIRED)
+set(PYBIND11_PYTHON_VERSION ${Python_VERSION})
+set(PYBIND11_FINDPYTHON ON) # Use new FindPython if available
+find_package(pybind11 CONFIG REQUIRED)
+
+# gdal requirement
+find_package(GDAL REQUIRED)
+
+# python module
+pybind11_add_module(PyCuAmpcor
+    src/PyCuAmpcor.cpp
+    src/GDALImage.cpp
+    src/cuAmpcorController.cpp
+    src/cuAmpcorParameter.cpp
+    src/cuArrays.cpp
+    src/cuArraysCopy.cu
+    src/cuArraysPadding.cu
+    src/cuCorrFrequency.cu
+    src/cuCorrNormalization.cu
+    src/cuCorrNormalizationSAT.cu
+    src/cuCorrNormalizer.cpp
+    src/cuCorrTimeDomain.cu
+    src/cuDeramp.cu
+    src/cuEstimateStats.cu
+    src/cuOffset.cu
+    src/cuOverSampler.cpp
+    src/cuSincOverSampler.cu
+    src/cudaError.cpp
+    src/cudaUtil.cpp
+    src/cuAmpcorProcessor.cpp
+    src/cuAmpcorProcessorROIPAC.cpp
+    src/cuAmpcorProcessorGrIMP.cpp
+    )
+
+# set_property(TARGET PyCuAmpcor PROPERTY CUDA_SEPARABLE_COMPILATION ON)
+
+target_include_directories(PyCuAmpcor PRIVATE
+    src
+    )
+
+target_link_libraries(PyCuAmpcor PRIVATE
+    CUDA::cufft
+    GDAL::GDAL
+    )
+
+if(NOT DEFINED PYTHON_MODULE_DIR)
+    set(PYTHON_MODULE_DIR packages CACHE PATH
+        "Python module directory (relative to install prefix)")
+endif()
+
+
+# Install with skbuild
+if (SKBUILD)
+    # Install package
+    install(TARGETS PyCuAmpcor
+            LIBRARY DESTINATION PyCuAmpcor)
+else()
+    # Install package
+    set(PYCUAMOCOR_INSTALL_DIR ${Python_SITEARCH}/PyCuAmpcor)
+
+    # Install compiled library
+    install(TARGETS PyCuAmpcor
+            LIBRARY DESTINATION ${PYCUAMOCOR_INSTALL_DIR})
+
+    # Install __init__.py
+    install(FILES ${CMAKE_SOURCE_DIR}/__init__.py
+            DESTINATION ${PYCUAMOCOR_INSTALL_DIR})
+    message(STATUS "not using skbuild, standard install.")
+endif()
+
+
+
+

contrib/PyCuAmpcor/README.md

@@ -1,4 +1,4 @@
-# PyCuAmpcor - Amplitude Cross-Correlation with GPU
+# PyCuAmpcor - Amplitude Cross-Correlation with GPU (Version 2)
 
 ## Contents
 
@@ -24,7 +24,12 @@ In practice, we
 
 A detailed formulation can be found, e.g., by J. P. Lewis with [the frequency domain approach](http://scribblethink.org/Work/nvisionInterface/nip.html).
 
-PyCuAmpcor follows the same procedure as the FORTRAN code, ampcor.F, in ROIPAC. In order to optimize the performance on GPU, some implementations are slightly different. In the [list the procedures](#5-list-of-procedures), we show the detailed steps of PyCuAmpcor, as well as their differences.
+PyCuAmpcor follows the same procedure as the FORTRAN code, ampcor.F, in ROIPAC. ROIPAC/ampcor workflow adopts a two-pass procedure. In the first pass, the cross-correlation is performed over the entire search range, but on the raw image, to obtain an estimated offset $(r_x^0, r_y^0)$ (in units of a pixel). In the second pass, another cross-correlation is performed over a much smaller search range centered at $(r_x^0, r_y^0)$, on (anti-aliasing) oversampled windows. The obtained correlation surface is further oversampled to achieve sub-pixel resolutions.  
+
+In order to optimize the performance on GPU, some implementations are slightly different. In the [list the procedures](#5-list-of-procedures), we show the detailed steps of this workflow, as well as its difference to the Fortran code.
+
+In PyCuAmpcor Version 2, we have also added a workflow following [GrIMP/SpeckleSource](https://github.com/fastice/speckleSource). It skips the first pass on ROIPAC/ampcor, while starting from the second pass with (anti-aliasing) oversampled secondary chip over the entire search range. It increases computational cost but is found to be more accurate in cases with low signal to noise ratio (SNR).   
+
 
 ## 2. Installation
 
@@ -55,10 +60,15 @@ You may also install PyCuAmpcor as a standalone package.
 
 ```bash
     # go to PyCuAmpcor source directory
-    cd contrib/PyCuAmpcor/src
-    # edit Makefile to provide the correct gdal include path and gpu architecture to NVCCFLAGS
-    # call make to compile
-    make
+    cd contrib/PyCuAmpcor/
+    cp CMakeLists.txt.standalone CMakeLists.txt
+    # follow CMake routines
+    mkdir build && cd build
+    cmake .. -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX \
+      -DCMAKE_CUDA_ARCHITECTURES=native \
+      -DCMAKE_PREFIX_PATH=${CONDA_PREFIX} \
+      -DCMAKE_BUILD_TYPE=Release 
+    make -j && make install
  ```
 
 ## 3. User Guide
@@ -83,25 +93,30 @@ ww=64  # template window width
 wh=64  # template window height
 sw=20   # (half) search range along width
 sh=20   # (half) search range along height
-kw=300   # skip between windows along width
-kh=100   # skip between windows along height
+kw=32   # skip between windows along width
+kh=32   # skip between windows along height
 mm=0   # margin to be neglected 
 gross=0  # whether to use a varying gross offset
 azshift=0 # constant gross offset along height/azimuth 
 rgshift=0 # constant gross offset along width/range
 deramp=0 # 0 for mag (TOPS), 1 for complex linear ramp, 2 for complex no deramping  
-oo=32  # correlation surface oversampling factor 
+oo=16  # correlation surface oversampling factor 
 outprefix=./merged/20151120_20151214/offset  # output prefix
 outsuffix=_ww64_wh64   # output suffix
 gpuid=0   # GPU device ID
 nstreams=2 # number of CUDA streams
 usemmap=1 # whether to use memory-map i/o
-mmapsize=8 # buffer size in GB for memory map
-nwac=32 # number of windows in a batch along width
+mmapsize=4 # buffer size in GB for memory map
+nwac=8 # number of windows in a batch along width
 nwdc=1  # number of windows in a batch along height
+workflow=0 # 0=ROIPAC/ampcor 1=GrIMP/SpeckleSource
 
 rm $outprefix$outsuffix*
-cuDenseOffsets.py --reference $reference --secondary $secondary --ww $ww --wh $wh --sw $sw --sh $sh --mm $mm --kw $kw --kh $kh --gross $gross --rr $rgshift --aa $azshift --oo $oo --deramp $deramp --outprefix $outprefix --outsuffix $outsuffix --gpuid $gpuid  --usemmap $usemmap --mmapsize $mmapsize --nwac $nwac --nwdc $nwdc 
+cuDenseOffsets.py --reference $reference --secondary $secondary --ww $ww --wh $wh --sw $sw --sh $sh \
+    --mm $mm --kw $kw --kh $kh --gross $gross --rr $rgshift --aa $azshift --oo $oo \
+    --deramp $deramp --outprefix $outprefix --outsuffix $outsuffix \
+    --gpuid $gpuid  --usemmap $usemmap --mmapsize $mmapsize \
+    --nwac $nwac --nwdc $nwdc --workflow $workflow
  ```
 
 Note that in PyCuAmpcor, the following names for directions are equivalent:
@@ -116,6 +131,24 @@ If you are interested in a particular region instead of the whole image, you may
 --startpixelac $startPixelAcross --startpixeldw $startPixelDown --nwa $numberOfWindowsAcross --nwd $numberOfWindowsDown
 ```
 
+If you prefer to provide the location of the last pixel instead of the number of windows, you may use
+
+```
+--startpixelac $startPixelAcross --startpixeldw $startPixelDown --endpixelac $endPixelAcross --endpixeldw $endPixelDown
+```
+
+The program will calculate the number of windows based on the skip size and the window size. 
+
+If you prefer to provide the locations for the center of the starting/ending window, rather than those of the first/last pixels, you may add *--use-center 1* together with the above options, e.g., 
+
+```
+--startpixelac $startPixelAcross --startpixeldw $startPixelDown --use-center 1 --nwa $numberOfWindowsAcross --nwd $numberOfWindowsDown
+```
+
+Here, *($startPixelDown, $startPixelAcross)* is the coordinate of the first window center, and will be fixed if you vary the window size or search range. 
+
+But please note, the provided starting/ending pixels may go beyond the original reference/secondary image range. The program will still run but the offset results on the edge may be incorrect.
+
 PyCuAmpcor supports two types of gross offset fields,
 * static (--gross=0), i.e., a constant shift between reference and secondary images. The static gross offsets can be passed by *--rr $rgshift --aa $azshift*. Note that the margin as well as the starting pixel may be adjusted.
 * dynamic (--gross=1), i.e., shifts between reference windows and secondary windows are varying in different locations. This is helpful to reduce the search range if you have a prior knowledge of the estimated offset fields, e.g., the velocity model of glaciers. You may prepare a BIP input file of the varying gross offsets (same format as the output offset fields), and use the option *--gross-file $grossOffsetFilename*. If you need the coordinates of reference windows, you may run *cuDenseOffsets.py* at first to find out the location of the starting pixel and the total number of windows. The coordinate for the starting pixel of the (iDown, iAcross) window will be (startPixelDown+iDown\*skipDown, startPixelAcross+iAcross\*skipAcross).
@@ -163,7 +196,7 @@ objOffset.checkPixelInImageRange()
 objOffset.runAmpcor()
 ```
 
-## 4. List of Parameters
+## 4. List of Parameters (ROIPAC/ampcor workflow)
 
 **Image Parameters**