diff --git a/.circleci/config.yml b/.circleci/config.yml index 444447f..607ca5b 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -77,7 +77,7 @@ jobs: pwd mkdir config build install . /opt/conda/bin/activate root - conda install --yes cython gdal h5py libgdal pytest numpy fftw scipy scons hdf4 hdf5 netcdf4 libgcc libstdcxx-ng cmake astropy + conda install --yes cython gdal h5py libgdal pytest numpy fftw scipy scons hdf4 hdf5 netcdf4 libgcc libstdcxx-ng cmake astropy opencv yum install -y uuid-devel x11-devel motif-devel jq gcc-gfortran ln -s /opt/conda/bin/cython /opt/conda/bin/cython3 cd /opt/conda/lib diff --git a/CMakeLists.txt b/CMakeLists.txt index c526ec4..7529e8a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,6 +18,8 @@ find_package(Python 3.5 REQUIRED COMPONENTS Interpreter Development find_package(FFTW REQUIRED) find_package(Motif) find_package(OpenMP REQUIRED COMPONENTS C CXX Fortran) +find_package(OpenCV COMPONENTS core highgui imgproc) + # Find these, and create IMPORTED INTERFACE libraries for them if they exist include(TargetGDAL) include(TargetMotif) diff --git a/contrib/CMakeLists.txt b/contrib/CMakeLists.txt index e1c8031..4395db0 100644 --- a/contrib/CMakeLists.txt +++ b/contrib/CMakeLists.txt @@ -4,6 +4,7 @@ add_subdirectory(demUtils) add_subdirectory(frameUtils) #add_subdirectory(unwUtils) add_subdirectory(downsample_unwrapper) +add_subdirectory(geo_autoRIFT) add_subdirectory(PyCuAmpcor) add_subdirectory(splitSpectrum) diff --git a/contrib/SConscript b/contrib/SConscript index bfc08a5..11366b4 100644 --- a/contrib/SConscript +++ b/contrib/SConscript @@ -80,6 +80,4 @@ SConscript('splitSpectrum/SConscript') SConscript('alos2filter/SConscript') SConscript('alos2proc/SConscript') SConscript('alos2proc_f/SConscript') - -if os.path.exists('geo_autoRIFT'): - SConscript('geo_autoRIFT/SConscript') +SConscript('geo_autoRIFT/SConscript') diff --git a/contrib/geo_autoRIFT/CMakeLists.txt b/contrib/geo_autoRIFT/CMakeLists.txt new file mode 100644 index 0000000..45b35cf --- /dev/null +++ b/contrib/geo_autoRIFT/CMakeLists.txt @@ -0,0 +1,12 @@ +# Early exit if prereqs not found +if(NOT TARGET Python::NumPy +OR NOT TARGET GDAL::GDAL +OR NOT OpenCV_FOUND + ) + return() +endif() + +InstallSameDir(__init__.py) + +add_subdirectory(autoRIFT) +add_subdirectory(geogrid) diff --git a/contrib/geo_autoRIFT/README.md b/contrib/geo_autoRIFT/README.md new file mode 100644 index 0000000..4e5a84d --- /dev/null +++ b/contrib/geo_autoRIFT/README.md @@ -0,0 +1,240 @@ +# autoRIFT (autonomous Repeat Image Feature Tracking) + + + +[![Language](https://img.shields.io/badge/python-3.6%2B-blue.svg)](https://www.python.org/) +[![Latest version](https://img.shields.io/badge/latest%20version-v1.0.6-yellowgreen.svg)](https://github.com/leiyangleon/autoRIFT/releases) +[![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://github.com/leiyangleon/autoRIFT/blob/master/LICENSE) +[![Citation](https://img.shields.io/badge/DOI-10.5281/zenodo.3756192-blue)](https://doi.org/10.5281/zenodo.3756192) + + + +**A Python module of a fast and intelligent algorithm for finding the pixel displacement between two images** + +**autoRIFT can be installed as a standalone Python module (only supports Cartesian coordinates) either manually or as a conda install (https://github.com/conda-forge/autorift-feedstock). To allow support for both Cartesian and radar coordinates, autoRIFT must be installed with the InSAR Scientific Computing Environment (ISCE: https://github.com/isce-framework/isce2)** + +**Use cases include all dense feature tracking applications, including the measurement of surface displacements occurring between two repeat satellite images as a result of glacier flow, large earthquake displacements, and land slides** + +**autoRIFT can be used for dense feature tracking between two images over a grid defined in an arbitrary geographic Cartesian (northing/easting) coordinate projection when used in combination with the sister Geogrid Python module (https://github.com/leiyangleon/Geogrid). Example applications include searching radar-coordinate imagery on a polar stereographic grid and searching Universal Transverse Mercator (UTM) imagery at a specified geographic Cartesian (northing/easting) coordinate grid** + + +Copyright (C) 2019 California Institute of Technology. Government Sponsorship Acknowledged. + +Link: https://github.com/leiyangleon/autoRIFT + +Citation: https://doi.org/10.5281/zenodo.3756192 + + +## 1. Authors + +Alex Gardner (JPL/Caltech; alex.s.gardner@jpl.nasa.gov) first described the algorithm "auto-RIFT" in (Gardner et al., 2018) and developed the first version in MATLAB; + +Yang Lei (GPS/Caltech; ylei@caltech.edu) translated it to Python, further optimized and incoporated to the ISCE software; + +Piyush Agram (JPL/Caltech; piyush.agram@jpl.nasa.gov) set up the installation as a standalone module and further cleaned the code. + +Reference: Gardner, A. S., Moholdt, G., Scambos, T., Fahnstock, M., Ligtenberg, S., van den Broeke, M., & Nilsson, J. (2018). Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. The Cryosphere, 12(2), 521–547. https://doi.org/10.5194/tc-12-521-2018 + + +## 2. Acknowledgement + +This effort was funded by the NASA MEaSUREs program in contribution to the Inter-mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) project (https://its-live.jpl.nasa.gov/) and through Alex Gardner’s participation in the NASA NISAR Science Team + + +## 3. Features + +* fast algorithm that finds displacement between the two images (where source and template image patches are extracted and correlated) using sparse search and iteratively progressive chip (template) sizes +* for intelligent use, user can specify unstructured spatially-varying search centers (center of the source image patch), search offsets (center displacement of the template image patch compared to the source image patch), chip sizes (size of template image patch), search ranges (size of the source image patch) +* faster than the conventional dense feature tracking "ampcor"/"denseampcor" algorithm included in ISCE by nearly 2 orders of magnitude +* supports various preprocessing modes for specified image pair, e.g. either the raw image (both texture and topography) or the texture only (high-frequency components without the topography) can be used with various choices of the high-pass filter options +* supports image data format of unsigned integer 8 (uint8; faster) and single-precision float (float32) +* user can adjust all of the relevant parameters, as detailed below in the instructions +* a Normalized Displacement Coherence (NDC) filter is developed (Gardner et al., 2018) to filter the chip displacemnt results based on displacement difference thresholds that are scaled to the local search range +* a sparse integer displacement search is used to narrow the template search range and to eliminate areas of "low coherence" for fine decimal displacement search and following chip size iterations +* novel nested grid design allows chip size progressively increase iteratively where a smaller chip size failed +* a light interpolation is done to fill the missing (unreliable) chip displacement results using median filtering and bicubic-mode interpolation (that can remove pixel discrepancy when using other modes) and an interpolation mask is returned +* the core displacement estimator (normalized cross correlation) and image pre-processing call OpenCV's Python and/or C++ functions for efficiency +* sub-pixel displacement estimation using a fast Gaussian pyramid upsampling algorithm that is robust to sub-pixel locking; the intelligent (chip size-dependent) selection of oversampling ratio is also supported for a tradeoff between efficiency and accuracy +* the current version can be installed with the ISCE software (that supports both Cartesian and radar coordinates) or as a standalone Python module (Cartesian coordinates only) +* when used in combination with the Geogrid Python module (https://github.com/leiyangleon/Geogrid), autoRIFT can be used for feature tracking between image pair over a grid defined in an arbitrary geographic Cartesian (northing/easting) coordinate projection +* when the grid is provided in geographic Cartesian (northing/easting) coordinates, outputs are returned in geocoded GeoTIFF image file format with the same EPSG projection code as input search grid + +## 4. Possible Future Development + +* for radar (SAR) images, it is yet to include the complex correlation of the two images, i.e. the current version only uses the amplitude correlation +* the current version works for single-core CPU while the multi-core parallelization or GPU implementation would be useful to extend + + +## 5. Demo + +### 5.1 Optical image over regular grid in imaging coordinates + + + +***Output of "autoRIFT" module for a pair of Landsat-8 images (20170708-20170724; same as the Demo dataset at https://github.com/leiyangleon/Geogrid) in Greenland over a regular-spacing grid: (a) estimated horizontal pixel displacement, (b) estimated vertical pixel displacement, (c) chip size used, (d) light interpolation mask.*** + +This is done by implementing the following command line: + +Standalone: + + testautoRIFT.py -m I1 -s I2 -fo 1 + +With ISCE: + + testautoRIFT_ISCE.py -m I1 -s I2 -fo 1 + +where "I1" and "I2" are the reference and test images as defined in the instructions below. The "-fo" option indicates whether or not to read optical image data. + + +### 5.2 Optical image over user-defined geographic Cartesian coordinate grid + + + +***Output of "autoRIFT" module for a pair of Landsat-8 images (20170708-20170724; same as the Demo dataset at https://github.com/leiyangleon/Geogrid) in Greenland over user-defined geographic Cartesian (northing/easting) coordinate grid (same grid used in the Demo at https://github.com/leiyangleon/Geogrid): (a) estimated horizontal pixel displacement, (b) estimated vertical pixel displacement, (c) light interpolation mask, (d) chip size used. Notes: all maps are established exactly over the same geographic Cartesian coordinate grid from input.*** + +This is done by implementing the following command line: + +Standalone: + + testautoRIFT.py -m I1 -s I2 -g winlocname -o winoffname -sr winsrname -csmin wincsminname -csmax wincsmaxname -vx winro2vxname -vy winro2vyname -fo 1 + +With ISCE: + + testautoRIFT_ISCE.py -m I1 -s I2 -g winlocname -o winoffname -sr winsrname -csmin wincsminname -csmax wincsmaxname -vx winro2vxname -vy winro2vyname -fo 1 + +where "I1" and "I2" are the reference and test images as defined in the instructions below, and the optional inputs "winlocname", "winoffname", "winsrname", "wincsminname", "wincsmaxname", "winro2vxname", "winro2vyname" are outputs from running "testGeogrid_ISCE.py" or "testGeogridOptical.py" as defined at https://github.com/leiyangleon/Geogrid. When "winlocname" (grid location) is specified, each of the rest optional input ("win * name") can be either used or omitted. + + + + +***Final motion velocity results by combining outputs from "Geogrid" (i.e. matrix of conversion coefficients from the Demo at https://github.com/leiyangleon/Geogrid) and "autoRIFT" modules: (a) estimated motion velocity from Landsat-8 data (x-direction; in m/yr), (b) motion velocity from input data (x-direction; in m/yr), (c) estimated motion velocity from Landsat-8 data (y-direction; in m/yr), (d) motion velocity from input data (y-direction; in m/yr). Notes: all maps are established exactly over the same geographic Cartesian (northing/easting) coordinate grid from input.*** + + +### 5.3 Radar image over regular grid in imaging coordinates + + + +***Output of "autoRIFT" module for a pair of Sentinel-1A/B images (20170221-20170227; same as the Demo dataset at https://github.com/leiyangleon/Geogrid) at Jakobshavn Glacier of Greenland over a regular-spacing grid: (a) estimated range pixel displacement, (b) estimated azimuth pixel displacement, (c) chip size used, (d) light interpolation mask.*** + + +This is obtained by implementing the following command line: + +With ISCE: + + testautoRIFT_ISCE.py -m I1 -s I2 + +where "I1" and "I2" are the reference and test images as defined in the instructions below. + + +### 5.4 Radar image over user-defined geographic Cartesian coordinate grid + + + +***Output of "autoRIFT" module for a pair of Sentinel-1A/B images (20170221-20170227; same as the Demo dataset at https://github.com/leiyangleon/Geogrid) at Jakobshavn Glacier of Greenland over user-defined geographic Cartesian (northing/easting) coordinate grid (same grid used in the Demo at https://github.com/leiyangleon/Geogrid): (a) estimated range pixel displacement, (b) estimated azimuth pixel displacement, (c) light interpolation mask, (d) chip size used. Notes: all maps are established exactly over the same geographic Cartesian coordinate grid from input.*** + +This is done by implementing the following command line: + +With ISCE: + + testautoRIFT_ISCE.py -m I1 -s I2 -g winlocname -o winoffname -sr winsrname -csmin wincsminname -csmax wincsmaxname -vx winro2vxname -vy winro2vyname + +where "I1" and "I2" are the reference and test images as defined in the instructions below, and the optional inputs "winlocname", "winoffname", "winsrname", "wincsminname", "wincsmaxname", "winro2vxname", "winro2vyname" are outputs from running "testGeogrid.py" as defined at https://github.com/leiyangleon/Geogrid. When "winlocname" (grid location) is specified, each of the rest optional input ("win * name") can be either used or omitted. + +**Runtime comparison for this test (on an OS X system with 2.9GHz Intel Core i7 processor and 16GB RAM):** +* __autoRIFT (single core): 10 mins__ +* __Dense ampcor from ISCE (8 cores): 90 mins__ + + + + +***Final motion velocity results by combining outputs from "Geogrid" (i.e. matrix of conversion coefficients from the Demo at https://github.com/leiyangleon/Geogrid) and "autoRIFT" modules: (a) estimated motion velocity from Sentinel-1 data (x-direction; in m/yr), (b) motion velocity from input data (x-direction; in m/yr), (c) estimated motion velocity from Sentinel-1 data (y-direction; in m/yr), (d) motion velocity from input data (y-direction; in m/yr). Notes: all maps are established exactly over the same geographic Cartesian (northing/easting) coordinate grid from input.*** + + +## 6. Instructions + +**Note:** + +* When the grid is provided in geographic Cartesian (northing/easting) coordinates (geographic coordinates lat/lon is not supported), it is required to run the "Geogrid" module (https://github.com/leiyangleon/Geogrid) first before running "autoRIFT". In other words, the outputs from "testGeogrid_ISCE.py" or "testGeogridOptical.py" (a.k.a "winlocname", "winoffname", "winsrname", "wincsminname", "wincsmaxname", "winro2vxname", "winro2vyname") will serve as optional inputs for running "autoRIFT" at a geographic grid, which is required to generate the final motion velocity maps. +* When the outputs from running the "Geogrid" module are not provided, a regular grid in the imaging coordinates will be automatically assigned + +**For quick use:** + +* Refer to the file "testautoRIFT.py" (standalone) and "testautoRIFT_ISCE.py" (with ISCE) for the usage of the module and modify it for your own purpose +* Input files include the reference image (required), test image (required), and the outputs from running "testGeogrid_ISCE.py" or "testGeogridOptical.py" (optional; a.k.a "winlocname", "winoffname", "winsrname", "wincsminname", "wincsmaxname", "winro2vxname", "winro2vyname"). When "winlocname" (grid location) is specified, each of the rest optional input ("win * name") can be either used or omitted. +* Output files include 1) estimated horizontal displacement (equivalent to range for radar), 2) estimated vertical displacement (equivalent to minus azimuth for radar), 3) light interpolation mask, 4) iteratively progressive chip size used. + +_Note: These four output files will be stored in a file named "offset.mat" that can be viewed in Python and MATLAB. When the grid is provided in geographic Cartesian (northing/easting) coordinates, a 4-band GeoTIFF with the same EPSG code as input grid will be created as well and named "offset.tif"; a 2-band GeoTIFF of the final converted motion velocity in geographic x- (easting) and y- (northing) coordinates will be created and named "velocity.tif". Also, it is possible to save the outputs in netCDF standard format by adding the "-nc" option to the "testautoRIFT.py" (standalone) and "testautoRIFT_ISCE.py" (with ISCE) command._ + +**For modular use:** + +* In Python environment, type the following to import the "autoRIFT" module and initialize the "autoRIFT" object + +_With ISCE:_ + + import isce + from contrib.geo_autoRIFT.autoRIFT import autoRIFT_ISCE + obj = autoRIFT_ISCE() + obj.configure() + +_Standalone:_ + + from autoRIFT import autoRIFT + obj = autoRIFT() + +* The "autoRIFT" object has several inputs that have to be assigned (listed below; can also be obtained by referring to "testautoRIFT.py"): + + ------------------input------------------ + I1 reference image (extracted image patches defined as "source") + I2 test image (extracted image patches defined as "template"; displacement = motion vector of I2 relative to I1) + xGrid [units = integer image pixels] horizontal pixel index at each grid point + yGrid [units = integer image pixels] vertical pixel index at each grid point + (if xGrid and yGrid not provided, a regular grid spanning the entire image will be automatically set up, which is similar to the conventional ISCE module, "ampcor" or "denseampcor") + Dx0 [units = integer image pixels] horizontal "downstream" reach location (that specifies the horizontal pixel displacement of the template's search center relative to the source's) at each grid point + Dy0 [units = integer image pixels] vertical "downstream" reach location (that specifies the vertical pixel displacement of the template's search center relative to the source's) at each grid point + (if Dx0 and Dy0 not provided, an array with zero values will be automatically assigned and there will be no offsets of the search centers) + +* After the inputs are specified, run the module as below + + obj.preprocess_filt_XXX() or obj.preprocess_db() + obj.uniform_data_type() + obj.runAutorift() + +where "XXX" can be "wal" for the Wallis filter, "hps" for the trivial high-pass filter, "sob" for the Sobel filter, "lap" for the Laplacian filter, and also a logarithmic operator without filtering is adopted for occasions where low-frequency components (i.e. topography) are desired, i.e. "obj.preprocess_db()". + +* The "autoRIFT" object has the following four outputs: + + ------------------output------------------ + Dx [units = decimal image pixels] estimated horizontal pixel displacement + Dy [units = decimal image pixels] estimated vertical pixel displacement + InterpMask [unitless; boolean data type] light interpolation mask + ChipSizeX [units = integer image pixels] iteratively progressive chip size in horizontal direction (different chip size allowed for vertical) + +* The "autoRIFT" object has many parameters that can be flexibly tweaked by the users for their own purpose (listed below; can also be obtained by referring to "geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py"): + + ------------------parameter list: general function------------------ + ChipSizeMinX [units = integer image pixels] Minimum size (in horizontal direction) of the template (chip) to correlate (default = 32; could be scalar or array with same dimension as xGrid) + ChipSizeMaxX [units = integer image pixels] Maximum size (in horizontal direction) of the template (chip) to correlate (default = 64; could be scalar or array with same dimension as xGrid) + ChipSize0X [units = integer image pixels] Minimum acceptable size (in horizontal direction) of the template (chip) to correlate (default = 32) + ScaleChipSizeY [unitless; integer data type] Scaling factor to get the vertically-directed chip size in reference to the horizontally-directed size (default = 1) + SearchLimitX [units = integer image pixels] Range (in horizontal direction) to search for displacement in the source (default = 25; could be scalar or array with same dimension as xGrid; when provided in array, set its elements to 0 if excluded for finding displacement) + SearchLimitY [units = integer image pixels] Range (in vertical direction) to search for displacement in the source (default = 25; could be scalar or array with same dimension as xGrid; when provided in array, set its elements to 0 if excluded for finding displacement) + SkipSampleX [units = integer image pixels] Number of samples to skip between search windows in horizontal direction for automatically-created grid if not specified by the user (default = 32) + SkipSampleY [units = integer image pixels] Number of lines to skip between search windows in vertical direction for automatically-created grid if not specified by the user (default = 32) + minSearch [units = integer image pixels] Minimum search limit (default = 6) + + ------------------parameter list: about Normalized Displacement Coherence (NDC) filter ------------------ + FracValid Fraction of valid displacements (default = 8/25) to be multiplied by filter window size "FiltWidth^2" and then used for thresholding the number of chip displacements that have passed the "FracSearch" disparity filtering + FracSearch Fraction of search limit used as threshold for disparity filtering of the chip displacement difference that is normalized by the search limit (default = 0.25) + FiltWidth Disparity filter width (default = 5) + Iter Number of iterations (default = 3) + MadScalar Scalar to be multiplied by Mad used as threshold for disparity filtering of the chip displacement deviation from the median (default = 4) + + ------------------parameter list: miscellaneous------------------ + WallisFilterWidth Width of the filter to be used for the preprocessing (default = 21) + fillFiltWidth Light interpolation filling filter width (default = 3) + sparseSearchSampleRate downsampling rate for sparse search (default = 4) + BuffDistanceC Buffer coarse correlation mask by this many pixels for use as fine search mask (default = 8) + CoarseCorCutoff Coarse correlation search cutoff (default = 0.01) + OverSampleRatio Factor for pyramid upsampling for sub-pixel level offset refinement (default = 16; can be scalar or Python dictionary for intelligent use, i.e. chip size-dependent, such as {ChipSize_1:OverSampleRatio_1,..., ChipSize_N:OverSampleRatio_N}) + DataTypeInput Image data type: 0 -> uint8, 1 -> float32 (default = 0) + zeroMask Force the margin (no data) to zeros which is useful for Wallis-filter-preprocessed images (default = None; 1 for Wallis filter) diff --git a/contrib/geo_autoRIFT/SConscript b/contrib/geo_autoRIFT/SConscript new file mode 100755 index 0000000..e49f703 --- /dev/null +++ b/contrib/geo_autoRIFT/SConscript @@ -0,0 +1,36 @@ +#!/usr/bin/env python + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. +# United States Government Sponsorship acknowledged. Any commercial use must be +# negotiated with the Office of Technology Transfer at the California Institute of +# Technology. This software is subject to U.S. export control laws and regulations +# and has been classified as EAR99. By accepting this software, the user agrees to +# comply with all applicable U.S. export laws and regulations. User has the +# responsibility to obtain export licenses, or other export authority as may be +# required before exporting such information to foreign countries or providing +# access to foreign persons. +# +# Author: Yang Lei +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + + + +import os +Import('envcontrib') +package = 'geo_autoRIFT' +envgeoAutorift = envcontrib.Clone() +envgeoAutorift['PACKAGE'] = envcontrib['PACKAGE'] + '/' + package +install = envcontrib['PRJ_SCONS_INSTALL'] + '/' + envgeoAutorift['PACKAGE'] +listFiles = ['__init__.py'] +envgeoAutorift.Install(install,listFiles) +envgeoAutorift.Alias('install',install) +Export('envgeoAutorift') + +autorift='autoRIFT/SConscript' +SConscript(autorift) + +geogrid='geogrid/SConscript' +SConscript(geogrid) diff --git a/contrib/geo_autoRIFT/__init__.py b/contrib/geo_autoRIFT/__init__.py new file mode 100755 index 0000000..f363046 --- /dev/null +++ b/contrib/geo_autoRIFT/__init__.py @@ -0,0 +1,37 @@ +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. +# United States Government Sponsorship acknowledged. Any commercial use must be +# negotiated with the Office of Technology Transfer at the California Institute of +# Technology. This software is subject to U.S. export control laws and regulations +# and has been classified as EAR99. By accepting this software, the user agrees to +# comply with all applicable U.S. export laws and regulations. User has the +# responsibility to obtain export licenses, or other export authority as may be +# required before exporting such information to foreign countries or providing +# access to foreign persons. +# +# Author: Yang Lei +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + + + +def createAutorift(name=''): + from contrib.geo_autoRIFT.Autorift import Autorift + return Autorift(name=name) + +def createGeogrid(name=''): + from contrib.geo_autoRIFT.Geogrid import Geogrid + return Geogrid(name=name) + + +def getFactoriesInfo(): + return {'Autorift': + { + 'factory':'createAutorift' + }, + 'Geogrid': + { + 'factory':'createGeogrid' + } + } diff --git a/contrib/geo_autoRIFT/autoRIFT/CMakeLists.txt b/contrib/geo_autoRIFT/autoRIFT/CMakeLists.txt new file mode 100644 index 0000000..55ec5a0 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/CMakeLists.txt @@ -0,0 +1,18 @@ +Python_add_library(autoriftcore MODULE + bindings/autoriftcoremodule.cpp + ) +target_include_directories(autoriftcore PRIVATE + include + ${OpenCV_INCLUDE_DIRS} + ) +target_link_libraries(autoriftcore PRIVATE + Python::NumPy + ${OpenCV_LIBS} + ) + +InstallSameDir( + autoriftcore + __init__.py + autoRIFT_ISCE.py + autoRIFT.py + ) diff --git a/contrib/geo_autoRIFT/autoRIFT/SConscript b/contrib/geo_autoRIFT/autoRIFT/SConscript new file mode 100755 index 0000000..a9099d9 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/SConscript @@ -0,0 +1,35 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Yang Lei +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeoAutorift') +envautorift = envgeoAutorift.Clone() +package = envautorift['PACKAGE'] +project = 'autoRIFT' +envautorift['PROJECT'] = project +install = envautorift['PRJ_SCONS_INSTALL'] + '/' + package + '/' + project +initFile = '__init__.py' + + +listFiles = ['autoRIFT.py','autoRIFT_ISCE.py',initFile] +envautorift.Install(install,listFiles) +envautorift.Alias('install',install) +Export('envautorift') +bindingsScons = 'bindings/SConscript' +SConscript(bindingsScons,variant_dir = envautorift['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/bindings') +includeScons = 'include/SConscript' +SConscript(includeScons) +#srcScons = 'src/SConscript' +#SConscript(srcScons,variant_dir = envautorift['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/src') diff --git a/contrib/geo_autoRIFT/autoRIFT/__init__.py b/contrib/geo_autoRIFT/autoRIFT/__init__.py new file mode 100755 index 0000000..08a6351 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/__init__.py @@ -0,0 +1,14 @@ +#!/usr/bin/env python + +#def SplitRangeSpectrum(): +# from .splitSpectrum import PySplitRangeSpectrum +# return PySplitRangeSpectrum() + +# should always work - standalone or with ISCE +from .autoRIFT import autoRIFT + +try: + from .autoRIFT_ISCE import autoRIFT_ISCE +except ImportError: + # this means ISCE support not available. Don't raise error. Allow standalone use + pass diff --git a/contrib/geo_autoRIFT/autoRIFT/autoRIFT.py b/contrib/geo_autoRIFT/autoRIFT/autoRIFT.py new file mode 100755 index 0000000..976a29b --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/autoRIFT.py @@ -0,0 +1,1093 @@ +#!/usr/bin/env python3 + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2019 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, +# end user, or in support of a prohibited end use). By downloading this software, +# the user agrees to comply with all applicable U.S. export laws and regulations. +# The user has the responsibility to obtain export licenses, or other export +# authority as may be required before exporting this software to any 'EAR99' +# embargoed foreign country or citizen of those countries. +# +# Author: Yang Lei +# +# Note: this is based on the MATLAB code, "auto-RIFT", written by Alex Gardner, +# and has been translated to Python and further optimized. +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + +import pdb +import subprocess +import re +import string +import sys + + + +class autoRIFT: + ''' + Class for mapping regular geographic grid on radar imagery. + ''' + + + + def preprocess_filt_wal(self): + ''' + Do the pre processing using wallis filter (10 min vs 15 min in Matlab). + ''' + import cv2 + import numpy as np +# import scipy.io as sio + + + if self.zeroMask is not None: + self.zeroMask = (self.I1 == 0) + + kernel = np.ones((self.WallisFilterWidth,self.WallisFilterWidth), dtype=np.float32) + + m = cv2.filter2D(self.I1,-1,kernel,borderType=cv2.BORDER_CONSTANT)/np.sum(kernel) + + m2 = (self.I1)**2 + + m2 = cv2.filter2D(m2,-1,kernel,borderType=cv2.BORDER_CONSTANT)/np.sum(kernel) + + s = np.sqrt(m2 - m**2) * np.sqrt(np.sum(kernel)/(np.sum(kernel)-1.0)) + + self.I1 = (self.I1 - m) / s + +# pdb.set_trace() + + m = cv2.filter2D(self.I2,-1,kernel,borderType=cv2.BORDER_CONSTANT)/np.sum(kernel) + + m2 = (self.I2)**2 + + m2 = cv2.filter2D(m2,-1,kernel,borderType=cv2.BORDER_CONSTANT)/np.sum(kernel) + + s = np.sqrt(m2 - m**2) * np.sqrt(np.sum(kernel)/(np.sum(kernel)-1.0)) + + self.I2 = (self.I2 - m) / s + + + +# #### obsolete definition of "preprocess_filt_hps" +# def preprocess_filt_hps(self): +# ''' +# Do the pre processing using (orig - low-pass filter) = high-pass filter filter (3.9/5.3 min). +# ''' +# import cv2 +# import numpy as np +# +# if self.zeroMask is not None: +# self.zeroMask = (self.I1 == 0) +# +# kernel = np.ones((self.WallisFilterWidth,self.WallisFilterWidth), dtype=np.float32) +# +# lp = cv2.filter2D(self.I1,-1,kernel,borderType=cv2.BORDER_CONSTANT)/np.sum(kernel) +# +# self.I1 = (self.I1 - lp) +# +# lp = cv2.filter2D(self.I2,-1,kernel,borderType=cv2.BORDER_CONSTANT)/np.sum(kernel) +# +# self.I2 = (self.I2 - lp) + + + def preprocess_filt_hps(self): + ''' + Do the pre processing using (orig - low-pass filter) = high-pass filter filter (3.9/5.3 min). + ''' + import cv2 + import numpy as np + + + if self.zeroMask is not None: + self.zeroMask = (self.I1 == 0) + + kernel = -np.ones((self.WallisFilterWidth,self.WallisFilterWidth), dtype=np.float32) + + kernel[int((self.WallisFilterWidth-1)/2),int((self.WallisFilterWidth-1)/2)] = kernel.size - 1 + + kernel = kernel / kernel.size + +# pdb.set_trace() + + self.I1 = cv2.filter2D(self.I1,-1,kernel,borderType=cv2.BORDER_CONSTANT) + + self.I2 = cv2.filter2D(self.I2,-1,kernel,borderType=cv2.BORDER_CONSTANT) + + + + def preprocess_db(self): + ''' + Do the pre processing using db scale (4 min). + ''' + import cv2 + import numpy as np + + + if self.zeroMask is not None: + self.zeroMask = (self.I1 == 0) + +# pdb.set_trace() + + self.I1 = 20.0 * np.log10(self.I1) + + self.I2 = 20.0 * np.log10(self.I2) + + + + + def preprocess_filt_sob(self): + ''' + Do the pre processing using sobel filter (4.5/5.8 min). + ''' + import cv2 + import numpy as np + + + + if self.zeroMask is not None: + self.zeroMask = (self.I1 == 0) + + sobelx = cv2.getDerivKernels(1,0,self.WallisFilterWidth) + + kernelx = np.outer(sobelx[0],sobelx[1]) + + sobely = cv2.getDerivKernels(0,1,self.WallisFilterWidth) + + kernely = np.outer(sobely[0],sobely[1]) + + kernel = kernelx + kernely + + self.I1 = cv2.filter2D(self.I1,-1,kernel,borderType=cv2.BORDER_CONSTANT) + + self.I2 = cv2.filter2D(self.I2,-1,kernel,borderType=cv2.BORDER_CONSTANT) + + + + + + + + + + def preprocess_filt_lap(self): + ''' + Do the pre processing using Laplacian filter (2.5 min / 4 min). + ''' + import cv2 + import numpy as np + + + if self.zeroMask is not None: + self.zeroMask = (self.I1 == 0) + + self.I1 = 20.0 * np.log10(self.I1) + self.I1 = cv2.Laplacian(self.I1,-1,ksize=self.WallisFilterWidth,borderType=cv2.BORDER_CONSTANT) + + self.I2 = 20.0 * np.log10(self.I2) + self.I2 = cv2.Laplacian(self.I2,-1,ksize=self.WallisFilterWidth,borderType=cv2.BORDER_CONSTANT) + + + + + + + + + + def uniform_data_type(self): + + import numpy as np + + if self.DataType == 0: +# validData = np.logical_not(np.isnan(self.I1)) + validData = np.isfinite(self.I1) + S1 = np.std(self.I1[validData])*np.sqrt(self.I1[validData].size/(self.I1[validData].size-1.0)) + M1 = np.mean(self.I1[validData]) + self.I1 = (self.I1 - (M1 - 3*S1)) / (6*S1) * (2**8 - 0) + +# self.I1[np.logical_not(np.isfinite(self.I1))] = 0 + self.I1 = np.round(np.clip(self.I1, 0, 255)).astype(np.uint8) + +# validData = np.logical_not(np.isnan(self.I2)) + validData = np.isfinite(self.I2) + S2 = np.std(self.I2[validData])*np.sqrt(self.I2[validData].size/(self.I2[validData].size-1.0)) + M2 = np.mean(self.I2[validData]) + self.I2 = (self.I2 - (M2 - 3*S2)) / (6*S2) * (2**8 - 0) + +# self.I2[np.logical_not(np.isfinite(self.I2))] = 0 + self.I2 = np.round(np.clip(self.I2, 0, 255)).astype(np.uint8) + + if self.zeroMask is not None: + self.I1[self.zeroMask] = 0 + self.I2[self.zeroMask] = 0 + self.zeroMask = None + + elif self.DataType == 1: + self.I1[np.logical_not(np.isfinite(self.I1))] = 0 + self.I2[np.logical_not(np.isfinite(self.I2))] = 0 + self.I1 = self.I1.astype(np.float32) + self.I2 = self.I2.astype(np.float32) + + if self.zeroMask is not None: + self.I1[self.zeroMask] = 0 + self.I2[self.zeroMask] = 0 + self.zeroMask = None + + else: + sys.exit('invalid data type for the image pair which must be unsigned integer 8 or 32-bit float') + + + def autorift(self): + ''' + Do the actual processing. + ''' + import numpy as np + import cv2 + from scipy import ndimage + + + ChipSizeUniX = np.unique(np.append(np.unique(self.ChipSizeMinX), np.unique(self.ChipSizeMaxX))) + ChipSizeUniX = np.delete(ChipSizeUniX,np.where(ChipSizeUniX == 0)[0]) + + if np.any(np.mod(ChipSizeUniX,self.ChipSize0X) != 0): + sys.exit('chip sizes must be even integers of ChipSize0') + + ChipRangeX = self.ChipSize0X * np.array([1,2,4,8,16,32,64],np.float32) +# ChipRangeX = ChipRangeX[ChipRangeX < (2**8 - 1)] + if np.max(ChipSizeUniX) > np.max(ChipRangeX): + sys.exit('max each chip size is out of range') + + ChipSizeUniX = ChipRangeX[(ChipRangeX >= np.min(ChipSizeUniX)) & (ChipRangeX <= np.max(ChipSizeUniX))] + + maxScale = np.max(ChipSizeUniX) / self.ChipSize0X + + if (np.mod(self.xGrid.shape[0],maxScale) != 0)|(np.mod(self.xGrid.shape[1],maxScale) != 0): + message = 'xgrid and ygrid have an incorect size ' + str(self.xGrid.shape) + ' for nested search, they must have dimensions that an interger multiple of ' + str(maxScale) + sys.exit(message) + + self.xGrid = self.xGrid.astype(np.float32) + self.yGrid = self.yGrid.astype(np.float32) + + if np.size(self.Dx0) == 1: + self.Dx0 = np.ones(self.xGrid.shape, np.float32) * np.round(self.Dx0) + else: + self.Dx0 = self.Dx0.astype(np.float32) + if np.size(self.Dy0) == 1: + self.Dy0 = np.ones(self.xGrid.shape, np.float32) * np.round(self.Dy0) + else: + self.Dy0 = self.Dy0.astype(np.float32) + if np.size(self.SearchLimitX) == 1: + self.SearchLimitX = np.ones(self.xGrid.shape, np.float32) * np.round(self.SearchLimitX) + else: + self.SearchLimitX = self.SearchLimitX.astype(np.float32) + if np.size(self.SearchLimitY) == 1: + self.SearchLimitY = np.ones(self.xGrid.shape, np.float32) * np.round(self.SearchLimitY) + else: + self.SearchLimitY = self.SearchLimitY.astype(np.float32) + if np.size(self.ChipSizeMinX) == 1: + self.ChipSizeMinX = np.ones(self.xGrid.shape, np.float32) * np.round(self.ChipSizeMinX) + else: + self.ChipSizeMinX = self.ChipSizeMinX.astype(np.float32) + if np.size(self.ChipSizeMaxX) == 1: + self.ChipSizeMaxX = np.ones(self.xGrid.shape, np.float32) * np.round(self.ChipSizeMaxX) + else: + self.ChipSizeMaxX = self.ChipSizeMaxX.astype(np.float32) + + ChipSizeX = np.zeros(self.xGrid.shape, np.float32) + InterpMask = np.zeros(self.xGrid.shape, np.bool) + Dx = np.empty(self.xGrid.shape, dtype=np.float32) + Dx.fill(np.nan) + Dy = np.empty(self.xGrid.shape, dtype=np.float32) + Dy.fill(np.nan) + + Flag = 3 + + for i in range(ChipSizeUniX.__len__()): + + # Nested grid setup: chip size being ChipSize0X no need to resize, otherwise has to resize the arrays + if self.ChipSize0X != ChipSizeUniX[i]: + Scale = self.ChipSize0X / ChipSizeUniX[i] + dstShape = (int(self.xGrid.shape[0]*Scale),int(self.xGrid.shape[1]*Scale)) + xGrid0 = cv2.resize(self.xGrid.astype(np.float32),dstShape[::-1],interpolation=cv2.INTER_AREA) + yGrid0 = cv2.resize(self.yGrid.astype(np.float32),dstShape[::-1],interpolation=cv2.INTER_AREA) + + if np.mod(ChipSizeUniX[i],2) == 0: + xGrid0 = np.round(xGrid0+0.5)-0.5 + yGrid0 = np.round(yGrid0+0.5)-0.5 + else: + xGrid0 = np.round(xGrid0) + yGrid0 = np.round(yGrid0) + + M0 = (ChipSizeX == 0) & (self.ChipSizeMinX <= ChipSizeUniX[i]) & (self.ChipSizeMaxX >= ChipSizeUniX[i]) + M0 = colfilt(M0.copy(), (int(1/Scale*6), int(1/Scale*6)), 0) + M0 = cv2.resize(np.logical_not(M0).astype(np.uint8),dstShape[::-1],interpolation=cv2.INTER_NEAREST).astype(np.bool) + + SearchLimitX0 = colfilt(self.SearchLimitX.copy(), (int(1/Scale), int(1/Scale)), 0) + colfilt(self.Dx0.copy(), (int(1/Scale), int(1/Scale)), 4) + SearchLimitY0 = colfilt(self.SearchLimitY.copy(), (int(1/Scale), int(1/Scale)), 0) + colfilt(self.Dy0.copy(), (int(1/Scale), int(1/Scale)), 4) + Dx00 = colfilt(self.Dx0.copy(), (int(1/Scale), int(1/Scale)), 2) + Dy00 = colfilt(self.Dy0.copy(), (int(1/Scale), int(1/Scale)), 2) + + SearchLimitX0 = np.ceil(cv2.resize(SearchLimitX0,dstShape[::-1])) + SearchLimitY0 = np.ceil(cv2.resize(SearchLimitY0,dstShape[::-1])) + SearchLimitX0[M0] = 0 + SearchLimitY0[M0] = 0 + Dx00 = np.round(cv2.resize(Dx00,dstShape[::-1],interpolation=cv2.INTER_NEAREST)) + Dy00 = np.round(cv2.resize(Dy00,dstShape[::-1],interpolation=cv2.INTER_NEAREST)) +# pdb.set_trace() + else: + SearchLimitX0 = self.SearchLimitX.copy() + SearchLimitY0 = self.SearchLimitY.copy() + Dx00 = self.Dx0.copy() + Dy00 = self.Dy0.copy() + xGrid0 = self.xGrid.copy() + yGrid0 = self.yGrid.copy() +# M0 = (ChipSizeX == 0) & (self.ChipSizeMinX <= ChipSizeUniX[i]) & (self.ChipSizeMaxX >= ChipSizeUniX[i]) +# SearchLimitX0[np.logical_not(M0)] = 0 +# SearchLimitY0[np.logical_not(M0)] = 0 + + if np.logical_not(np.any(SearchLimitX0 != 0)): + continue + + idxZero = (SearchLimitX0 <= 0) | (SearchLimitY0 <= 0) + SearchLimitX0[idxZero] = 0 + SearchLimitY0[idxZero] = 0 + SearchLimitX0[(np.logical_not(idxZero)) & (SearchLimitX0 < self.minSearch)] = self.minSearch + SearchLimitY0[(np.logical_not(idxZero)) & (SearchLimitY0 < self.minSearch)] = self.minSearch + + if ((xGrid0.shape[0] - 2)/self.sparseSearchSampleRate < 5) | ((xGrid0.shape[1] - 2)/self.sparseSearchSampleRate < 5): + Flag = 2 + return Flag + + # Setup for coarse search: sparse sampling / resize + rIdxC = slice(self.sparseSearchSampleRate-1,xGrid0.shape[0],self.sparseSearchSampleRate) + cIdxC = slice(self.sparseSearchSampleRate-1,xGrid0.shape[1],self.sparseSearchSampleRate) + xGrid0C = xGrid0[rIdxC,cIdxC] + yGrid0C = yGrid0[rIdxC,cIdxC] + +# pdb.set_trace() + + if np.remainder(self.sparseSearchSampleRate,2) == 0: + filtWidth = self.sparseSearchSampleRate + 1 + else: + filtWidth = self.sparseSearchSampleRate + + SearchLimitX0C = colfilt(SearchLimitX0.copy(), (int(filtWidth), int(filtWidth)), 0) + SearchLimitY0C = colfilt(SearchLimitY0.copy(), (int(filtWidth), int(filtWidth)), 0) + SearchLimitX0C = SearchLimitX0C[rIdxC,cIdxC] + SearchLimitY0C = SearchLimitY0C[rIdxC,cIdxC] + + Dx0C = Dx00[rIdxC,cIdxC] + Dy0C = Dy00[rIdxC,cIdxC] + + # Coarse search + SubPixFlag = False + ChipSizeXC = ChipSizeUniX[i] + ChipSizeYC = np.float32(np.round(ChipSizeXC*self.ScaleChipSizeY/2)*2) + + if type(self.OverSampleRatio) is dict: + overSampleRatio = self.OverSampleRatio[ChipSizeUniX[i]] + else: + overSampleRatio = self.OverSampleRatio + +# pdb.set_trace() + + if self.I1.dtype == np.uint8: + DxC, DyC = arImgDisp_u(self.I2.copy(), self.I1.copy(), xGrid0C.copy(), yGrid0C.copy(), ChipSizeXC, ChipSizeYC, SearchLimitX0C.copy(), SearchLimitY0C.copy(), Dx0C.copy(), Dy0C.copy(), SubPixFlag, overSampleRatio) + elif self.I1.dtype == np.float32: + DxC, DyC = arImgDisp_s(self.I2.copy(), self.I1.copy(), xGrid0C.copy(), yGrid0C.copy(), ChipSizeXC, ChipSizeYC, SearchLimitX0C.copy(), SearchLimitY0C.copy(), Dx0C.copy(), Dy0C.copy(), SubPixFlag, overSampleRatio) + else: + sys.exit('invalid data type for the image pair which must be unsigned integer 8 or 32-bit float') + +# pdb.set_trace() + + # M0C is the mask for reliable estimates after coarse search, MC is the mask after disparity filtering, MC2 is the mask after area closing for fine search + M0C = np.logical_not(np.isnan(DxC)) + + DispFiltC = DISP_FILT() + if ChipSizeUniX[i] == ChipSizeUniX[0]: + DispFiltC.FracValid = self.FracValid + DispFiltC.FracSearch = self.FracSearch + DispFiltC.FiltWidth = self.FiltWidth + DispFiltC.Iter = self.Iter + DispFiltC.MadScalar = self.MadScalar + DispFiltC.Iter = DispFiltC.Iter - 1 + MC = DispFiltC.filtDisp(DxC.copy(), DyC.copy(), SearchLimitX0C.copy(), SearchLimitY0C.copy(), M0C.copy(), overSampleRatio) + + MC[np.logical_not(M0C)] = False + + ROIC = (SearchLimitX0C > 0) + CoarseCorValidFac = np.sum(MC[ROIC]) / np.sum(M0C[ROIC]) + if (CoarseCorValidFac < self.CoarseCorCutoff): + continue + + MC2 = ndimage.distance_transform_edt(np.logical_not(MC)) < self.BuffDistanceC + dstShape = (int(MC2.shape[0]*self.sparseSearchSampleRate),int(MC2.shape[1]*self.sparseSearchSampleRate)) + + MC2 = cv2.resize(MC2.astype(np.uint8),dstShape[::-1],interpolation=cv2.INTER_NEAREST).astype(np.bool) +# pdb.set_trace() + if np.logical_not(np.all(MC2.shape == SearchLimitX0.shape)): + rowAdd = SearchLimitX0.shape[0] - MC2.shape[0] + colAdd = SearchLimitX0.shape[1] - MC2.shape[1] + if rowAdd>0: + MC2 = np.append(MC2,MC2[-rowAdd:,:],axis=0) + if colAdd>0: + MC2 = np.append(MC2,MC2[:,-colAdd:],axis=1) + + SearchLimitX0[np.logical_not(MC2)] = 0 + SearchLimitY0[np.logical_not(MC2)] = 0 + + # Fine Search + SubPixFlag = True + ChipSizeXF = ChipSizeUniX[i] + ChipSizeYF = np.float32(np.round(ChipSizeXF*self.ScaleChipSizeY/2)*2) +# pdb.set_trace() + if self.I1.dtype == np.uint8: + DxF, DyF = arImgDisp_u(self.I2.copy(), self.I1.copy(), xGrid0.copy(), yGrid0.copy(), ChipSizeXF, ChipSizeYF, SearchLimitX0.copy(), SearchLimitY0.copy(), Dx00.copy(), Dy00.copy(), SubPixFlag, overSampleRatio) + elif self.I1.dtype == np.float32: + DxF, DyF = arImgDisp_s(self.I2.copy(), self.I1.copy(), xGrid0.copy(), yGrid0.copy(), ChipSizeXF, ChipSizeYF, SearchLimitX0.copy(), SearchLimitY0.copy(), Dx00.copy(), Dy00.copy(), SubPixFlag, overSampleRatio) + else: + sys.exit('invalid data type for the image pair which must be unsigned integer 8 or 32-bit float') + +# pdb.set_trace() + + DispFiltF = DISP_FILT() + if ChipSizeUniX[i] == ChipSizeUniX[0]: + DispFiltF.FracValid = self.FracValid + DispFiltF.FracSearch = self.FracSearch + DispFiltF.FiltWidth = self.FiltWidth + DispFiltF.Iter = self.Iter + DispFiltF.MadScalar = self.MadScalar + + + M0 = DispFiltF.filtDisp(DxF.copy(), DyF.copy(), SearchLimitX0.copy(), SearchLimitY0.copy(), np.logical_not(np.isnan(DxF)), overSampleRatio) +# pdb.set_trace() + DxF[np.logical_not(M0)] = np.nan + DyF[np.logical_not(M0)] = np.nan + + # Light interpolation with median filtered values: DxFM (filtered) and DxF (unfiltered) + DxFM = colfilt(DxF.copy(), (self.fillFiltWidth, self.fillFiltWidth), 3) + DyFM = colfilt(DyF.copy(), (self.fillFiltWidth, self.fillFiltWidth), 3) + + # M0 is mask for original valid estimates, MF is mask for filled ones, MM is mask where filtered ones exist for filling + MF = np.zeros(M0.shape, dtype=np.bool) + MM = np.logical_not(np.isnan(DxFM)) + + for j in range(3): + foo = MF | M0 # initial valid estimates + foo1 = (cv2.filter2D(foo.astype(np.float32),-1,np.ones((3,3)),borderType=cv2.BORDER_CONSTANT) >= 6) | foo # 1st area closing followed by the 2nd (part of the next line calling OpenCV) +# pdb.set_trace() + fillIdx = np.logical_not(bwareaopen(np.logical_not(foo1).astype(np.uint8), 5)) & np.logical_not(foo) & MM + MF[fillIdx] = True + DxF[fillIdx] = DxFM[fillIdx] + DyF[fillIdx] = DyFM[fillIdx] + + # Below is for replacing the valid estimates with the bicubic filtered values for robust and accurate estimation + if self.ChipSize0X == ChipSizeUniX[i]: + Dx = DxF + Dy = DyF + ChipSizeX[M0|MF] = ChipSizeUniX[i] + InterpMask[MF] = True +# pdb.set_trace() + else: +# pdb.set_trace() + Scale = ChipSizeUniX[i] / self.ChipSize0X + dstShape = (int(Dx.shape[0]/Scale),int(Dx.shape[1]/Scale)) + + # DxF0 (filtered) / Dx (unfiltered) is the result from earlier iterations, DxFM (filtered) / DxF (unfiltered) is that of the current iteration + # first colfilt nans within 2-by-2 area (otherwise 1 nan will contaminate all 4 points) + DxF0 = colfilt(Dx.copy(),(int(Scale+1),int(Scale+1)),2) + # then resize to half size using area (similar to averaging) to match the current iteration + DxF0 = cv2.resize(DxF0,dstShape[::-1],interpolation=cv2.INTER_AREA) + DyF0 = colfilt(Dy.copy(),(int(Scale+1),int(Scale+1)),2) + DyF0 = cv2.resize(DyF0,dstShape[::-1],interpolation=cv2.INTER_AREA) + + # Note this DxFM is almost the same as DxFM (same variable) in the light interpolation (only slightly better); however, only small portion of it will be used later at locations specified by M0 and MF that are determined in the light interpolation. So even without the following two lines, the final Dx and Dy result is still the same. + # to fill out all of the missing values in DxF + DxFM = colfilt(DxF.copy(), (5,5), 3) + DyFM = colfilt(DyF.copy(), (5,5), 3) + + # fill the current-iteration result with previously determined reliable estimates that are not searched in the current iteration + idx = np.isnan(DxF) & np.logical_not(np.isnan(DxF0)) + DxFM[idx] = DxF0[idx] + DyFM[idx] = DyF0[idx] + + # Strong interpolation: use filtered estimates wherever the unfiltered estimates do not exist + idx = np.isnan(DxF) & np.logical_not(np.isnan(DxFM)) + DxF[idx] = DxFM[idx] + DyF[idx] = DyFM[idx] + + dstShape = (Dx.shape[0],Dx.shape[1]) + DxF = cv2.resize(DxF,dstShape[::-1],interpolation=cv2.INTER_CUBIC) + DyF = cv2.resize(DyF,dstShape[::-1],interpolation=cv2.INTER_CUBIC) + MF = cv2.resize(MF.astype(np.uint8),dstShape[::-1],interpolation=cv2.INTER_NEAREST).astype(np.bool) + M0 = cv2.resize(M0.astype(np.uint8),dstShape[::-1],interpolation=cv2.INTER_NEAREST).astype(np.bool) + + idxRaw = M0 & (ChipSizeX == 0) + idxFill = MF & (ChipSizeX == 0) + ChipSizeX[idxRaw | idxFill] = ChipSizeUniX[i] + InterpMask[idxFill] = True + Dx[idxRaw | idxFill] = DxF[idxRaw | idxFill] + Dy[idxRaw | idxFill] = DyF[idxRaw | idxFill] + + Flag = 1 + ChipSizeY = np.round(ChipSizeX * self.ScaleChipSizeY /2) * 2 + self.Dx = Dx + self.Dy = Dy + self.InterpMask = InterpMask + self.Flag = Flag + self.ChipSizeX = ChipSizeX + self.ChipSizeY = ChipSizeY + + + + + + + + + def runAutorift(self): + ''' + quick processing routine which calls autorift main function (user can define their own way by mimicing the workflow here). + ''' + import numpy as np + + + # truncate the grid to fit the nested grid + if np.size(self.ChipSizeMaxX) == 1: + chopFactor = self.ChipSizeMaxX / self.ChipSize0X + else: + chopFactor = np.max(self.ChipSizeMaxX) / self.ChipSize0X + rlim = int(np.floor(self.xGrid.shape[0] / chopFactor) * chopFactor) + clim = int(np.floor(self.xGrid.shape[1] / chopFactor) * chopFactor) + self.origSize = self.xGrid.shape +# pdb.set_trace() + self.xGrid = np.round(self.xGrid[0:rlim,0:clim]) + 0.5 + self.yGrid = np.round(self.yGrid[0:rlim,0:clim]) + 0.5 + + # truncate the initial offset as well if they exist + if np.size(self.Dx0) != 1: + self.Dx0 = self.Dx0[0:rlim,0:clim] + self.Dy0 = self.Dy0[0:rlim,0:clim] + + # truncate the search limits as well if they exist + if np.size(self.SearchLimitX) != 1: + self.SearchLimitX = self.SearchLimitX[0:rlim,0:clim] + self.SearchLimitY = self.SearchLimitY[0:rlim,0:clim] + + # truncate the chip sizes as well if they exist + if np.size(self.ChipSizeMaxX) != 1: + self.ChipSizeMaxX = self.ChipSizeMaxX[0:rlim,0:clim] + self.ChipSizeMinX = self.ChipSizeMinX[0:rlim,0:clim] + + # call autoRIFT main function + self.autorift() + + + + + + def __init__(self): + + super(autoRIFT, self).__init__() + + ##Input related parameters + self.I1 = None + self.I2 = None + self.xGrid = None + self.yGrid = None + self.Dx0 = 0 + self.Dy0 = 0 + self.origSize = None + self.zeroMask = None + + ##Output file + self.Dx = None + self.Dy = None + self.InterpMask = None + self.Flag = None + self.ChipSizeX = None + self.ChipSizeY = None + + ##Parameter list + self.WallisFilterWidth = 21 + self.ChipSizeMinX = 32 + self.ChipSizeMaxX = 64 + self.ChipSize0X = 32 + self.ScaleChipSizeY = 1 + self.SearchLimitX = 25 + self.SearchLimitY = 25 + self.SkipSampleX = 32 + self.SkipSampleY = 32 + self.fillFiltWidth = 3 + self.minSearch = 6 + self.sparseSearchSampleRate = 4 + self.FracValid = 8/25 + self.FracSearch = 0.25 + self.FiltWidth = 5 + self.Iter = 3 + self.MadScalar = 4 + self.BuffDistanceC = 8 + self.CoarseCorCutoff = 0.01 + self.OverSampleRatio = 16 + self.DataType = 0 + + + + + + +class AUTO_RIFT_CORE: + def __init__(self): + ##Pointer to C + self._autoriftcore = None + + + +def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, SearchLimitY, Dx0, Dy0, SubPixFlag, oversample): + + import numpy as np + from . import autoriftcore + + core = AUTO_RIFT_CORE() + if core._autoriftcore is not None: + autoriftcore.destroyAutoRiftCore_Py(core._autoriftcore) + + core._autoriftcore = autoriftcore.createAutoRiftCore_Py() + + +# if np.size(I1) == 1: +# if np.logical_not(isinstance(I1,np.uint8) & isinstance(I2,np.uint8)): +# sys.exit('input images must be uint8') +# else: +# if np.logical_not((I1.dtype == np.uint8) & (I2.dtype == np.uint8)): +# sys.exit('input images must be uint8') + + if np.size(SearchLimitX) == 1: + if np.logical_not(isinstance(SearchLimitX,np.float32) & isinstance(SearchLimitY,np.float32)): + sys.exit('SearchLimit must be float') + else: + if np.logical_not((SearchLimitX.dtype == np.float32) & (SearchLimitY.dtype == np.float32)): + sys.exit('SearchLimit must be float') + + if np.size(Dx0) == 1: + if np.logical_not(isinstance(Dx0,np.float32) & isinstance(Dy0,np.float32)): + sys.exit('Search offsets must be float') + else: + if np.logical_not((Dx0.dtype == np.float32) & (Dy0.dtype == np.float32)): + sys.exit('Search offsets must be float') + + if np.size(ChipSizeX) == 1: + if np.logical_not(isinstance(ChipSizeX,np.float32) & isinstance(ChipSizeY,np.float32)): + sys.exit('ChipSize must be float') + else: + if np.logical_not((ChipSizeX.dtype == np.float32) & (ChipSizeY.dtype == np.float32)): + sys.exit('ChipSize must be float') + + + + if np.any(np.mod(ChipSizeX,2) != 0) | np.any(np.mod(ChipSizeY,2) != 0): +# if np.any(np.mod(xGrid-0.5,1) == 0) & np.any(np.mod(yGrid-0.5,1) == 0): +# sys.exit('for an even chip size ImgDisp returns displacements centered at pixel boundaries so xGrid and yGrid must an inter - 1/2 [example: if you want the velocity centered between pixel (1,1) and (2,2) then specify a grid center of (1.5, 1.5)]') +# else: +# xGrid = np.ceil(xGrid) +# yGrid = np.ceil(yGrid) + sys.exit('it is better to have ChipSize = even number') + + if np.any(np.mod(SearchLimitX,1) != 0) | np.any(np.mod(SearchLimitY,1) != 0): + sys.exit('SearchLimit must be an integar value') + + if np.any(SearchLimitX < 0) | np.any(SearchLimitY < 0): + sys.exit('SearchLimit cannot be negative') + + if np.any(np.mod(ChipSizeX,4) != 0) | np.any(np.mod(ChipSizeY,4) != 0): + sys.exit('ChipSize should be evenly divisible by 4') + + if np.size(Dx0) == 1: + Dx0 = np.ones(xGrid.shape, dtype=np.float32) * Dx0 + + if np.size(Dy0) == 1: + Dy0 = np.ones(xGrid.shape, dtype=np.float32) * Dy0 + + if np.size(SearchLimitX) == 1: + SearchLimitX = np.ones(xGrid.shape, dtype=np.float32) * SearchLimitX + + if np.size(SearchLimitY) == 1: + SearchLimitY = np.ones(xGrid.shape, dtype=np.float32) * SearchLimitY + + if np.size(ChipSizeX) == 1: + ChipSizeX = np.ones(xGrid.shape, dtype=np.float32) * ChipSizeX + + if np.size(ChipSizeY) == 1: + ChipSizeY = np.ones(xGrid.shape, dtype=np.float32) * ChipSizeY + + # convert from cartesian X-Y to matrix X-Y: X no change, Y from up being positive to down being positive + Dy0 = -Dy0 + + SLx_max = np.max(SearchLimitX + np.abs(Dx0)) + Px = np.int(np.max(ChipSizeX)/2 + SLx_max + 2) + SLy_max = np.max(SearchLimitY + np.abs(Dy0)) + Py = np.int(np.max(ChipSizeY)/2 + SLy_max + 2) + + I1 = np.lib.pad(I1,((Py,Py),(Px,Px)),'constant') + I2 = np.lib.pad(I2,((Py,Py),(Px,Px)),'constant') + + # adjust center location by the padarray size and 0.5 is added because we need to extract the chip centered at X+1 with -chipsize/2:chipsize/2-1, which equivalently centers at X+0.5 (X is the original grid point location). So for even chipsize, always returns offset estimates at (X+0.5). + xGrid += (Px + 0.5) + yGrid += (Py + 0.5) + + Dx = np.empty(xGrid.shape,dtype=np.float32) + Dx.fill(np.nan) + Dy = Dx.copy() + + + for jj in range(xGrid.shape[1]): + if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,jj] == 0): + continue + Dx1 = Dx[:,jj] + Dy1 = Dy[:,jj] + for ii in range(xGrid.shape[0]): + if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0): + continue + + # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference) + clx = np.floor(ChipSizeX[ii,jj]/2) + ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj])) + cly = np.floor(ChipSizeY[ii,jj]/2) + ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj])) + ChipI = I2[ChipRangeY,ChipRangeX] + + SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj])) + SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj])) + RefI = I1[SearchRangeY,SearchRangeX] + + minChipI = np.min(ChipI) + if minChipI < 0: + ChipI = ChipI - minChipI + if np.all(ChipI == ChipI[0,0]): + continue + + minRefI = np.min(RefI) + if minRefI < 0: + RefI = RefI - minRefI + if np.all(RefI == RefI[0,0]): + continue + + + if SubPixFlag: + # call C++ + Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arSubPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(),oversample)) +# # call Python +# Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI) + else: + # call C++ + Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel())) +# # call Python +# Dx1[ii], Dy1[ii] = arPixDisp(ChipI,RefI) + + + # add back 1) I1 (RefI) relative to I2 (ChipI) initial offset Dx0 and Dy0, and + # 2) RefI relative to ChipI has a left/top boundary offset of -SearchLimitX and -SearchLimitY + idx = np.logical_not(np.isnan(Dx)) + Dx[idx] += (Dx0[idx] - SearchLimitX[idx]) + Dy[idx] += (Dy0[idx] - SearchLimitY[idx]) + + # convert from matrix X-Y to cartesian X-Y: X no change, Y from down being positive to up being positive + Dy = -Dy + + autoriftcore.destroyAutoRiftCore_Py(core._autoriftcore) + core._autoriftcore = None + + return Dx, Dy + + + + + + +def arImgDisp_s(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, SearchLimitY, Dx0, Dy0, SubPixFlag, oversample): + + import numpy as np + from . import autoriftcore + + core = AUTO_RIFT_CORE() + if core._autoriftcore is not None: + autoriftcore.destroyAutoRiftCore_Py(core._autoriftcore) + + core._autoriftcore = autoriftcore.createAutoRiftCore_Py() + + +# if np.size(I1) == 1: +# if np.logical_not(isinstance(I1,np.uint8) & isinstance(I2,np.uint8)): +# sys.exit('input images must be uint8') +# else: +# if np.logical_not((I1.dtype == np.uint8) & (I2.dtype == np.uint8)): +# sys.exit('input images must be uint8') + + if np.size(SearchLimitX) == 1: + if np.logical_not(isinstance(SearchLimitX,np.float32) & isinstance(SearchLimitY,np.float32)): + sys.exit('SearchLimit must be float') + else: + if np.logical_not((SearchLimitX.dtype == np.float32) & (SearchLimitY.dtype == np.float32)): + sys.exit('SearchLimit must be float') + + if np.size(Dx0) == 1: + if np.logical_not(isinstance(Dx0,np.float32) & isinstance(Dy0,np.float32)): + sys.exit('Search offsets must be float') + else: + if np.logical_not((Dx0.dtype == np.float32) & (Dy0.dtype == np.float32)): + sys.exit('Search offsets must be float') + + if np.size(ChipSizeX) == 1: + if np.logical_not(isinstance(ChipSizeX,np.float32) & isinstance(ChipSizeY,np.float32)): + sys.exit('ChipSize must be float') + else: + if np.logical_not((ChipSizeX.dtype == np.float32) & (ChipSizeY.dtype == np.float32)): + sys.exit('ChipSize must be float') + + + + if np.any(np.mod(ChipSizeX,2) != 0) | np.any(np.mod(ChipSizeY,2) != 0): +# if np.any(np.mod(xGrid-0.5,1) == 0) & np.any(np.mod(yGrid-0.5,1) == 0): +# sys.exit('for an even chip size ImgDisp returns displacements centered at pixel boundaries so xGrid and yGrid must an inter - 1/2 [example: if you want the velocity centered between pixel (1,1) and (2,2) then specify a grid center of (1.5, 1.5)]') +# else: +# xGrid = np.ceil(xGrid) +# yGrid = np.ceil(yGrid) + sys.exit('it is better to have ChipSize = even number') + + if np.any(np.mod(SearchLimitX,1) != 0) | np.any(np.mod(SearchLimitY,1) != 0): + sys.exit('SearchLimit must be an integar value') + + if np.any(SearchLimitX < 0) | np.any(SearchLimitY < 0): + sys.exit('SearchLimit cannot be negative') + + if np.any(np.mod(ChipSizeX,4) != 0) | np.any(np.mod(ChipSizeY,4) != 0): + sys.exit('ChipSize should be evenly divisible by 4') + + if np.size(Dx0) == 1: + Dx0 = np.ones(xGrid.shape, dtype=np.float32) * Dx0 + + if np.size(Dy0) == 1: + Dy0 = np.ones(xGrid.shape, dtype=np.float32) * Dy0 + + if np.size(SearchLimitX) == 1: + SearchLimitX = np.ones(xGrid.shape, dtype=np.float32) * SearchLimitX + + if np.size(SearchLimitY) == 1: + SearchLimitY = np.ones(xGrid.shape, dtype=np.float32) * SearchLimitY + + if np.size(ChipSizeX) == 1: + ChipSizeX = np.ones(xGrid.shape, dtype=np.float32) * ChipSizeX + + if np.size(ChipSizeY) == 1: + ChipSizeY = np.ones(xGrid.shape, dtype=np.float32) * ChipSizeY + + # convert from cartesian X-Y to matrix X-Y: X no change, Y from up being positive to down being positive + Dy0 = -Dy0 + + SLx_max = np.max(SearchLimitX + np.abs(Dx0)) + Px = np.int(np.max(ChipSizeX)/2 + SLx_max + 2) + SLy_max = np.max(SearchLimitY + np.abs(Dy0)) + Py = np.int(np.max(ChipSizeY)/2 + SLy_max + 2) + + I1 = np.lib.pad(I1,((Py,Py),(Px,Px)),'constant') + I2 = np.lib.pad(I2,((Py,Py),(Px,Px)),'constant') + + # adjust center location by the padarray size and 0.5 is added because we need to extract the chip centered at X+1 with -chipsize/2:chipsize/2-1, which equivalently centers at X+0.5 (X is the original grid point location). So for even chipsize, always returns offset estimates at (X+0.5). + xGrid += (Px + 0.5) + yGrid += (Py + 0.5) + + Dx = np.empty(xGrid.shape,dtype=np.float32) + Dx.fill(np.nan) + Dy = Dx.copy() + + + for jj in range(xGrid.shape[1]): + if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,jj] == 0): + continue + Dx1 = Dx[:,jj] + Dy1 = Dy[:,jj] + for ii in range(xGrid.shape[0]): + if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0): + continue + + # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference) + clx = np.floor(ChipSizeX[ii,jj]/2) + ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj])) + cly = np.floor(ChipSizeY[ii,jj]/2) + ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj])) + ChipI = I2[ChipRangeY,ChipRangeX] + + SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj])) + SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj])) + RefI = I1[SearchRangeY,SearchRangeX] + + minChipI = np.min(ChipI) + if minChipI < 0: + ChipI = ChipI - minChipI + if np.all(ChipI == ChipI[0,0]): + continue + + minRefI = np.min(RefI) + if minRefI < 0: + RefI = RefI - minRefI + if np.all(RefI == RefI[0,0]): + continue + + + if SubPixFlag: + # call C++ + Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arSubPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(), oversample)) +# # call Python +# Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI) + else: + # call C++ + Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel())) +# # call Python +# Dx1[ii], Dy1[ii] = arPixDisp(ChipI,RefI) + + + # add back 1) I1 (RefI) relative to I2 (ChipI) initial offset Dx0 and Dy0, and + # 2) RefI relative to ChipI has a left/top boundary offset of -SearchLimitX and -SearchLimitY + idx = np.logical_not(np.isnan(Dx)) + Dx[idx] += (Dx0[idx] - SearchLimitX[idx]) + Dy[idx] += (Dy0[idx] - SearchLimitY[idx]) + + # convert from matrix X-Y to cartesian X-Y: X no change, Y from down being positive to up being positive + Dy = -Dy + + autoriftcore.destroyAutoRiftCore_Py(core._autoriftcore) + core._autoriftcore = None + + return Dx, Dy + + + + +def colfilt(A, kernelSize, option): + + from skimage.util import view_as_windows as viewW + import numpy as np + + A = np.lib.pad(A,((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(int((kernelSize[1]-1)/2),int((kernelSize[1]-1)/2))),mode='constant',constant_values=np.nan) + + B = viewW(A, kernelSize).reshape(-1,kernelSize[0]*kernelSize[1]).T[:,::1] + + output_size = (A.shape[0]-kernelSize[0]+1,A.shape[1]-kernelSize[1]+1) + C = np.zeros(output_size,dtype=A.dtype) + if option == 0:# max + C = np.nanmax(B,axis=0).reshape(output_size) + elif option == 1:# min + C = np.nanmin(B,axis=0).reshape(output_size) + elif option == 2:# mean + C = np.nanmean(B,axis=0).reshape(output_size) + elif option == 3:# median + C = np.nanmedian(B,axis=0).reshape(output_size) + elif option == 4:# range + C = np.nanmax(B,axis=0).reshape(output_size) - np.nanmin(B,axis=0).reshape(output_size) + elif option == 6:# MAD (Median Absolute Deviation) + m = B.shape[0] + D = np.abs(B - np.dot(np.ones((m,1),dtype=A.dtype), np.array([np.nanmedian(B,axis=0)]))) + C = np.nanmedian(D,axis=0).reshape(output_size) + elif option[0] == 5:# displacement distance count with option[1] being the threshold + m = B.shape[0] + c = int(np.round((m + 1) / 2)-1) +# c = 0 + D = np.abs(B - np.dot(np.ones((m,1),dtype=A.dtype), np.array([B[c,:]]))) + C = np.sum(D= dToleranceX) & (colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), (5,self.FracSearch)) >= dToleranceY) + +# if self.Iter == 3: +# pdb.set_trace() + + for i in range(np.max([self.Iter-1,1])): + Dx[np.logical_not(M)] = np.nan + Dy[np.logical_not(M)] = np.nan + + DxMad = colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), 6) + DyMad = colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), 6) + + DxM = colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), 3) + DyM = colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), 3) + + M = (np.abs(Dx - DxM) <= np.maximum(self.MadScalar * DxMad, DxMadmin)) & (np.abs(Dy - DyM) <= np.maximum(self.MadScalar * DyMad, DyMadmin)) & M + + return M + + + + +def bwareaopen(image,size1): + + import numpy as np + from skimage import measure + + # now identify the objects and remove those above a threshold + labels, N = measure.label(image,connectivity=2,return_num=True) + label_size = [(labels == label).sum() for label in range(N + 1)] + + # now remove the labels + for label,size in enumerate(label_size): + if size < size1: + image[labels == label] = 0 + + return image + + + + diff --git a/contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py b/contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py new file mode 100755 index 0000000..7345e59 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py @@ -0,0 +1,236 @@ +#!/usr/bin/env python3 + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2019 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, +# end user, or in support of a prohibited end use). By downloading this software, +# the user agrees to comply with all applicable U.S. export laws and regulations. +# The user has the responsibility to obtain export licenses, or other export +# authority as may be required before exporting this software to any 'EAR99' +# embargoed foreign country or citizen of those countries. +# +# Author: Yang Lei +# +# Note: this is based on the MATLAB code, "auto-RIFT", written by Alex Gardner, +# and has been translated to Python and further optimized. +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import isce +from iscesys.Component.Component import Component +import pdb +import subprocess +import re +import string +import sys + +WALLIS_FILTER_WIDTH = Component.Parameter('WallisFilterWidth', + public_name='WALLIS_FILTER_WIDTH', + default = 21, + type = int, + mandatory = False, + doc = 'Width of the Wallis filter to be used for the pre-processing') + +CHIP_SIZE_MIN_X = Component.Parameter('ChipSizeMinX', + public_name='CHIP_SIZE_MIN', + default = 32, + type = int, + mandatory = False, + doc = 'Minimum size (in X direction) of the reference data window to be used for correlation') + +CHIP_SIZE_MAX_X = Component.Parameter('ChipSizeMaxX', + public_name='CHIP_SIZE_MAX', + default = 64, + type = int, + mandatory = False, + doc = 'Maximum size (in X direction) of the reference data window to be used for correlation') + +CHIP_SIZE_0X = Component.Parameter('ChipSize0X', + public_name='CHIP_SIZE_0X', + default = 32, + type = int, + mandatory = False, + doc = 'Minimum acceptable size (in X direction) of the reference data window to be used for correlation') + +SCALE_CHIP_SIZE_Y = Component.Parameter('ScaleChipSizeY', + public_name='SCALE_CHIP_SIZE_Y', + default = 1, + type = float, + mandatory = False, + doc = 'Scaling factor to get the Y-directed chip size in reference to the X-directed sizes') + +SEARCH_LIMIT_X = Component.Parameter('SearchLimitX', + public_name='SEARCH_LIMIT_X', + default = 25, + type = int, + mandatory = False, + doc = 'Limit (in X direction) of the search data window to be used for correlation') + +SEARCH_LIMIT_Y = Component.Parameter('SearchLimitY', + public_name='SEARCH_LIMIT_Y', + default = 25, + type = int, + mandatory = False, + doc = 'Limit (in Y direction) of the search data window to be used for correlation') + +SKIP_SAMPLE_X = Component.Parameter('SkipSampleX', + public_name = 'SKIP_SAMPLE_X', + default = 32, + type = int, + mandatory = False, + doc = 'Number of samples to skip between windows in X (range) direction.') + +SKIP_SAMPLE_Y = Component.Parameter('SkipSampleY', + public_name = 'SKIP_SAMPLE_Y', + default = 32, + type = int, + mandatory=False, + doc = 'Number of lines to skip between windows in Y ( "-" azimuth) direction.') + +FILL_FILT_WIDTH = Component.Parameter('fillFiltWidth', + public_name = 'FILL_FILT_WIDTH', + default = 3, + type = int, + mandatory=False, + doc = 'light interpolation Fill Filter width') + +MIN_SEARCH = Component.Parameter('minSearch', + public_name = 'MIN_SEARCH', + default = 6, + type = int, + mandatory=False, + doc = 'minimum search limit') + +SPARSE_SEARCH_SAMPLE_RATE = Component.Parameter('sparseSearchSampleRate', + public_name = 'SPARSE_SEARCH_SAMPLE_RATE', + default = 4, + type = int, + mandatory=False, + doc = 'sparse search sample rate') + +FRAC_VALID = Component.Parameter('FracValid', + public_name = 'FRAC_VALID', + default = 8/25, + type = float, + mandatory=False, + doc = 'Fraction of valid displacements') + +FRAC_SEARCH = Component.Parameter('FracSearch', + public_name = 'FRAC_SEARCH', + default = 0.25, + type = float, + mandatory=False, + doc = 'Fraction of search') + +FILT_WIDTH = Component.Parameter('FiltWidth', + public_name = 'FILT_WIDTH', + default = 5, + type = int, + mandatory=False, + doc = 'Disparity Filter width') + +ITER = Component.Parameter('Iter', + public_name = 'ITER', + default = 3, + type = int, + mandatory=False, + doc = 'Number of iterations') + +MAD_SCALAR = Component.Parameter('MadScalar', + public_name = 'MAD_SCALAR', + default = 4, + type = int, + mandatory=False, + doc = 'Mad Scalar') + +BUFF_DISTANCE_C = Component.Parameter('BuffDistanceC', + public_name = 'BUFF_DISTANCE_C', + default = 8, + type = int, + mandatory=False, + doc = 'buffer coarse corr mask by this many pixels for use as fine search mask') + +COARSE_COR_CUTOFF = Component.Parameter('CoarseCorCutoff', + public_name = 'COARSE_COR_CUTOFF', + default = 0.01, + type = float, + mandatory=False, + doc = 'coarse correlation search cutoff') + +OVER_SAMPLE_RATIO = Component.Parameter('OverSampleRatio', + public_name = 'OVER_SAMPLE_RATIO', + default = 16, + type = int, + mandatory=False, + doc = 'factor for pyramid up sampling for sub-pixel level offset refinement') + +DATA_TYPE = Component.Parameter('DataType', + public_name = 'DATA_TYPE', + default = 0, + type = int, + mandatory=False, + doc = 'Input data type: 0 -> uint8, 1 -> float32') + + +try: + # Try Autorift within ISCE first + from .autoRIFT import autoRIFT +except ImportError: + # Try global Autorift + from autoRIFT import autoRIFT +except: + raise Exception('Autorift does not appear to be installed.') + + + +class autoRIFT_ISCE(autoRIFT, Component): + ''' + Class for mapping regular geographic grid on radar imagery. + ''' + + parameter_list = (WALLIS_FILTER_WIDTH, + CHIP_SIZE_MIN_X, + CHIP_SIZE_MAX_X, + CHIP_SIZE_0X, + SCALE_CHIP_SIZE_Y, + SEARCH_LIMIT_X, + SEARCH_LIMIT_Y, + SKIP_SAMPLE_X, + SKIP_SAMPLE_Y, + FILL_FILT_WIDTH, + MIN_SEARCH, + SPARSE_SEARCH_SAMPLE_RATE, + FRAC_VALID, + FRAC_SEARCH, + FILT_WIDTH, + ITER, + MAD_SCALAR, + BUFF_DISTANCE_C, + COARSE_COR_CUTOFF, + OVER_SAMPLE_RATIO, + DATA_TYPE) + + + + + def __init__(self): + + super(autoRIFT_ISCE, self).__init__() + + + diff --git a/contrib/geo_autoRIFT/autoRIFT/bindings/SConscript b/contrib/geo_autoRIFT/autoRIFT/bindings/SConscript new file mode 100755 index 0000000..ee66f6c --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/bindings/SConscript @@ -0,0 +1,30 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Yang Lei +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envautorift') +package = envautorift['PACKAGE'] +project = envautorift['PROJECT'] +install = envautorift['PRJ_SCONS_INSTALL'] + '/' + package + '/' + project +build = envautorift['PRJ_SCONS_BUILD'] + '/' + package + '/' + project +#libList = ['gomp','autorift','combinedLib','gdal'] +#libList = ['gomp','combinedLib','gdal','opencv_core','opencv_highgui','opencv_imgproc'] +libList = ['opencv_core','opencv_highgui','opencv_imgproc'] +envautorift.PrependUnique(LIBS = libList) +module = envautorift.LoadableModule(target = 'autoriftcore.so', source = 'autoriftcoremodule.cpp') +envautorift.Install(install,module) +envautorift.Alias('install',install) +envautorift.Install(build,module) +envautorift.Alias('build',build) diff --git a/contrib/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp b/contrib/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp new file mode 100755 index 0000000..94238d9 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp @@ -0,0 +1,369 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 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, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Author: Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + + + + + + +#include +#include +//#include "autoriftcore.h" +#include "autoriftcoremodule.h" + + +#include "stdio.h" +#include "iostream" +#include "numpy/arrayobject.h" +#include "opencv2/imgproc/imgproc.hpp" +#include "opencv2/highgui/highgui.hpp" +#include "opencv2/core/core.hpp" + + +using namespace cv; +using namespace std; + +struct autoRiftCore{ +// This empty structure "autoRiftCore" in C++ is assgined to "self._autoriftcore" in python, which can take a set of variables in this file (declare here or in "autoriftcore.h" and set states below). For example, +// ((autoRiftCore*)(ptr))->widC = widC; +// ((autoRiftCore*)(ptr))->arPixDisp() +// If taking all the variables here in the structure, the complicated computation can be performed in another C++ file, "autoriftcore.cpp" (that includes functions like void autoRiftCore::arPixDisp()). +}; + + +static const char * const __doc__ = "Python extension for autoriftcore"; + + +PyModuleDef moduledef = { + //header + PyModuleDef_HEAD_INIT, + //name of the module + "autoriftcore", + //module documentation string + __doc__, + //size of the per-interpreter state of the module; + -1, + autoriftcore_methods, +}; + +//Initialization function for the module +PyMODINIT_FUNC +PyInit_autoriftcore() +{ + PyObject* module = PyModule_Create(&moduledef); + if (!module) + { + return module; + } + return module; +} + +PyObject* createAutoRiftCore(PyObject* self, PyObject *args) +{ + autoRiftCore* ptr = new autoRiftCore; + return Py_BuildValue("K", (uint64_t) ptr); +} + +PyObject* destroyAutoRiftCore(PyObject *self, PyObject *args) +{ + uint64_t ptr; + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + if (((autoRiftCore*)(ptr))!=NULL) + { + delete ((autoRiftCore*)(ptr)); + } + return Py_BuildValue("i", 0); +} + + + +PyObject* arPixDisp_u(PyObject *self, PyObject *args) +{ + uint64_t ptr; + PyArrayObject *ChipI, *RefI; + int widC, lenC; + int widR, lenR; + if (!PyArg_ParseTuple(args, "KiiOiiO", &ptr, &widC, &lenC, &ChipI, &widR, &lenR, &RefI)) + { + return NULL; + } + + uint8_t my_arrC[widC * lenC]; + + for(int i=0; i<(widC*lenC); i++){ + my_arrC[i] = (*(uint8_t *)PyArray_GETPTR1(ChipI,i)); + } + + uint8_t my_arrR[widR * lenR]; + + for(int i=0; i<(widR*lenR); i++){ + my_arrR[i] = (*(uint8_t *)PyArray_GETPTR1(RefI,i)); + } + + + cv::Mat my_imgC = cv::Mat(cv::Size(widC, lenC), CV_8UC1, &my_arrC); + cv::Mat my_imgR = cv::Mat(cv::Size(widR, lenR), CV_8UC1, &my_arrR); + + int result_cols = widR - widC + 1; + int result_rows = lenR - lenC + 1; + + cv::Mat result; + result.create( result_rows, result_cols, CV_32FC1 ); + + cv::matchTemplate( my_imgR, my_imgC, result, CV_TM_CCOEFF_NORMED ); + + cv::Point maxLoc; + cv::minMaxLoc(result, NULL, NULL, NULL, &maxLoc); + + double x = maxLoc.x; + double y = maxLoc.y; + + + return Py_BuildValue("dd", x, y); +} + + + + + +PyObject* arSubPixDisp_u(PyObject *self, PyObject *args) +{ + uint64_t ptr; + PyArrayObject *ChipI, *RefI; + int widC, lenC; + int widR, lenR; + int overSampleNC; + if (!PyArg_ParseTuple(args, "KiiOiiOi", &ptr, &widC, &lenC, &ChipI, &widR, &lenR, &RefI, &overSampleNC)) + { + return NULL; + } + + uint8_t my_arrC[widC * lenC]; + + for(int i=0; i<(widC*lenC); i++){ + my_arrC[i] = (*(uint8_t *)PyArray_GETPTR1(ChipI,i)); + } + + uint8_t my_arrR[widR * lenR]; + + for(int i=0; i<(widR*lenR); i++){ + my_arrR[i] = (*(uint8_t *)PyArray_GETPTR1(RefI,i)); + } + + + cv::Mat my_imgC = cv::Mat(cv::Size(widC, lenC), CV_8UC1, &my_arrC); + cv::Mat my_imgR = cv::Mat(cv::Size(widR, lenR), CV_8UC1, &my_arrR); + + int result_cols = widR - widC + 1; + int result_rows = lenR - lenC + 1; + + cv::Mat result; + result.create( result_rows, result_cols, CV_32FC1 ); + + cv::matchTemplate( my_imgR, my_imgC, result, CV_TM_CCOEFF_NORMED ); + + cv::Point maxLoc; + cv::minMaxLoc(result, NULL, NULL, NULL, &maxLoc); + + + // refine the offset at the sub-pixel level using image upsampling (pyramid algorithm): extract 5x5 small image at the coarse offset location + int x_start, y_start, x_count, y_count; + + x_start = cv::max(maxLoc.x-2, 0); + x_start = cv::min(x_start, result_cols-5); + x_count = 5; + + y_start = cv::max(maxLoc.y-2, 0); + y_start = cv::min(y_start, result_rows-5); + y_count = 5; + + + cv::Mat result_small (result, cv::Rect(x_start, y_start, x_count, y_count)); + + int cols = result_small.cols; + int rows = result_small.rows; + int overSampleFlag = 1; + cv::Mat predecessor_small = result_small; + cv::Mat foo; + + while (overSampleFlag < overSampleNC){ + cols *= 2; + rows *= 2; + overSampleFlag *= 2; + foo.create(cols, rows, CV_32FC1); + cv::pyrUp(predecessor_small, foo, cv::Size(cols, rows)); + predecessor_small = foo; + } + + cv::Point maxLoc_small; + cv::minMaxLoc(foo, NULL, NULL, NULL, &maxLoc_small); + + + double x = ((maxLoc_small.x + 0.0)/overSampleNC + x_start); + double y = ((maxLoc_small.y + 0.0)/overSampleNC + y_start); + + + return Py_BuildValue("dd", x, y); +} + + + + +PyObject* arPixDisp_s(PyObject *self, PyObject *args) +{ + uint64_t ptr; + PyArrayObject *ChipI, *RefI; + int widC, lenC; + int widR, lenR; + if (!PyArg_ParseTuple(args, "KiiOiiO", &ptr, &widC, &lenC, &ChipI, &widR, &lenR, &RefI)) + { + return NULL; + } + + float my_arrC[widC * lenC]; + + for(int i=0; i<(widC*lenC); i++){ + my_arrC[i] = (*(float *)PyArray_GETPTR1(ChipI,i)); + } + + float my_arrR[widR * lenR]; + + for(int i=0; i<(widR*lenR); i++){ + my_arrR[i] = (*(float *)PyArray_GETPTR1(RefI,i)); + } + + + cv::Mat my_imgC = cv::Mat(cv::Size(widC, lenC), CV_32FC1, &my_arrC); + cv::Mat my_imgR = cv::Mat(cv::Size(widR, lenR), CV_32FC1, &my_arrR); + + int result_cols = widR - widC + 1; + int result_rows = lenR - lenC + 1; + + cv::Mat result; + result.create( result_rows, result_cols, CV_32FC1 ); + + cv::matchTemplate( my_imgR, my_imgC, result, CV_TM_CCORR_NORMED ); + + cv::Point maxLoc; + cv::minMaxLoc(result, NULL, NULL, NULL, &maxLoc); + + double x = maxLoc.x; + double y = maxLoc.y; + + + return Py_BuildValue("dd", x, y); +} + + + + + +PyObject* arSubPixDisp_s(PyObject *self, PyObject *args) +{ + uint64_t ptr; + PyArrayObject *ChipI, *RefI; + int widC, lenC; + int widR, lenR; + int overSampleNC; + if (!PyArg_ParseTuple(args, "KiiOiiOi", &ptr, &widC, &lenC, &ChipI, &widR, &lenR, &RefI, &overSampleNC)) + { + return NULL; + } + + float my_arrC[widC * lenC]; + + for(int i=0; i<(widC*lenC); i++){ + my_arrC[i] = (*(float *)PyArray_GETPTR1(ChipI,i)); + } + + float my_arrR[widR * lenR]; + + for(int i=0; i<(widR*lenR); i++){ + my_arrR[i] = (*(float *)PyArray_GETPTR1(RefI,i)); + } + + + cv::Mat my_imgC = cv::Mat(cv::Size(widC, lenC), CV_32FC1, &my_arrC); + cv::Mat my_imgR = cv::Mat(cv::Size(widR, lenR), CV_32FC1, &my_arrR); + + int result_cols = widR - widC + 1; + int result_rows = lenR - lenC + 1; + + cv::Mat result; + result.create( result_rows, result_cols, CV_32FC1 ); + + cv::matchTemplate( my_imgR, my_imgC, result, CV_TM_CCORR_NORMED ); + + cv::Point maxLoc; + cv::minMaxLoc(result, NULL, NULL, NULL, &maxLoc); + + + // refine the offset at the sub-pixel level using image upsampling (pyramid algorithm): extract 5x5 small image at the coarse offset location + int x_start, y_start, x_count, y_count; + + x_start = cv::max(maxLoc.x-2, 0); + x_start = cv::min(x_start, result_cols-5); + x_count = 5; + + y_start = cv::max(maxLoc.y-2, 0); + y_start = cv::min(y_start, result_rows-5); + y_count = 5; + + + cv::Mat result_small (result, cv::Rect(x_start, y_start, x_count, y_count)); + + int cols = result_small.cols; + int rows = result_small.rows; + int overSampleFlag = 1; + cv::Mat predecessor_small = result_small; + cv::Mat foo; + + while (overSampleFlag < overSampleNC){ + cols *= 2; + rows *= 2; + overSampleFlag *= 2; + foo.create(cols, rows, CV_32FC1); + cv::pyrUp(predecessor_small, foo, cv::Size(cols, rows)); + predecessor_small = foo; + } + + cv::Point maxLoc_small; + cv::minMaxLoc(foo, NULL, NULL, NULL, &maxLoc_small); + + + double x = ((maxLoc_small.x + 0.0)/overSampleNC + x_start); + double y = ((maxLoc_small.y + 0.0)/overSampleNC + y_start); + + + + return Py_BuildValue("dd", x, y); +} diff --git a/contrib/geo_autoRIFT/autoRIFT/include/SConscript b/contrib/geo_autoRIFT/autoRIFT/include/SConscript new file mode 100755 index 0000000..6b447f9 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/include/SConscript @@ -0,0 +1,24 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Yang Lei +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envautorift') +package = envautorift['PACKAGE'] +project = envautorift['PROJECT'] +build = envautorift['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/include' +envautorift.AppendUnique(CPPPATH = [build]) +listFiles = ['autoriftcoremodule.h'] +envautorift.Install(build,listFiles) +envautorift.Alias('build',build) diff --git a/contrib/geo_autoRIFT/autoRIFT/include/autoriftcoremodule.h b/contrib/geo_autoRIFT/autoRIFT/include/autoriftcoremodule.h new file mode 100755 index 0000000..5cb8fb2 --- /dev/null +++ b/contrib/geo_autoRIFT/autoRIFT/include/autoriftcoremodule.h @@ -0,0 +1,58 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 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, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Author: Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + + +#ifndef autoriftcoremodule_h +#define autoriftcoremodule_h + +#include +#include + + +extern "C" +{ + PyObject * createAutoRiftCore(PyObject*, PyObject*); + PyObject * destroyAutoRiftCore(PyObject*, PyObject*); + PyObject * arPixDisp_u(PyObject *, PyObject *); + PyObject * arSubPixDisp_u(PyObject *, PyObject *); + PyObject * arPixDisp_s(PyObject *, PyObject *); + PyObject * arSubPixDisp_s(PyObject *, PyObject *); +} + +static PyMethodDef autoriftcore_methods[] = +{ + {"createAutoRiftCore_Py", createAutoRiftCore, METH_VARARGS, " "}, + {"destroyAutoRiftCore_Py", destroyAutoRiftCore, METH_VARARGS, " "}, + {"arPixDisp_u_Py", arPixDisp_u, METH_VARARGS, " "}, + {"arSubPixDisp_u_Py", arSubPixDisp_u, METH_VARARGS, " "}, + {"arPixDisp_s_Py", arPixDisp_s, METH_VARARGS, " "}, + {"arSubPixDisp_s_Py", arSubPixDisp_s, METH_VARARGS, " "}, + {NULL, NULL, 0, NULL} +}; +#endif //autoRiftCoremodule_h + diff --git a/contrib/geo_autoRIFT/geogrid/CMakeLists.txt b/contrib/geo_autoRIFT/geogrid/CMakeLists.txt new file mode 100644 index 0000000..d155e28 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/CMakeLists.txt @@ -0,0 +1,18 @@ +Python_add_library(geogrid MODULE + bindings/geogridmodule.cpp + src/geogrid.cpp + ) +target_include_directories(geogrid PRIVATE + include + ) +target_link_libraries(geogrid PRIVATE + GDAL::GDAL + combinedLib + ) + +InstallSameDir( + geogrid + __init__.py + Geogrid.py + GeogridOptical.py + ) diff --git a/contrib/geo_autoRIFT/geogrid/Geogrid.py b/contrib/geo_autoRIFT/geogrid/Geogrid.py new file mode 100755 index 0000000..77f9c6a --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/Geogrid.py @@ -0,0 +1,345 @@ +#!/usr/bin/env python3 + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2019 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, +# end user, or in support of a prohibited end use). By downloading this software, +# the user agrees to comply with all applicable U.S. export laws and regulations. +# The user has the responsibility to obtain export licenses, or other export +# authority as may be required before exporting this software to any 'EAR99' +# embargoed foreign country or citizen of those countries. +# +# Authors: Piyush Agram, Yang Lei +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + + +import isce +from iscesys.Component.Component import Component +import pdb +import subprocess +import re +import string + +class Geogrid(Component): + ''' + Class for mapping regular geographic grid on radar imagery. + ''' + + def geogrid(self): + ''' + Do the actual processing. + ''' + import isce + from components.contrib.geo_autoRIFT.geogrid import geogrid + + ##Determine appropriate EPSG system + self.epsg = self.getProjectionSystem() + + ###Determine extent of data needed + bbox = self.determineBbox() + + ###Load approrpriate DEM from database + if self.demname is None: + self.demname, self.dhdxname, self.dhdyname, self.vxname, self.vyname, self.srxname, self.sryname, self.csminxname, self.csminyname, self.csmaxxname, self.csmaxyname, self.ssmname = self.getDEM(bbox) + + + ##Create and set parameters + self.setState() + + ##Run + geogrid.geogrid_Py(self._geogrid) + + ##Clean up + self.finalize() + + def getProjectionSystem(self): + ''' + Testing with Greenland. + ''' + if not self.demname: + raise Exception('At least the DEM parameter must be set for geogrid') + + from osgeo import gdal, osr + ds = gdal.Open(self.demname, gdal.GA_ReadOnly) + srs = osr.SpatialReference() + srs.ImportFromWkt(ds.GetProjection()) + srs.AutoIdentifyEPSG() + ds = None +# pdb.set_trace() + + if srs.IsGeographic(): + epsgstr = srs.GetAuthorityCode('GEOGCS') + elif srs.IsProjected(): + epsgstr = srs.GetAuthorityCode('PROJCS') + elif srs.IsLocal(): + raise Exception('Local coordinate system encountered') + else: + raise Exception('Non-standard coordinate system encountered') + if not epsgstr: #Empty string->use shell command gdalsrsinfo for last trial + cmd = 'gdalsrsinfo -o epsg {0}'.format(self.demname) + epsgstr = subprocess.check_output(cmd, shell=True) +# pdb.set_trace() + epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0] +# pdb.set_trace() + if not epsgstr: #Empty string + raise Exception('Could not auto-identify epsg code') +# pdb.set_trace() + epsgcode = int(epsgstr) +# pdb.set_trace() + return epsgcode + + def determineBbox(self, zrange=[-200,4000]): + ''' + Dummy. + ''' + import numpy as np + import datetime + from osgeo import osr,gdal + +# import pdb +# pdb.set_trace() + +# rng = self.startingRange + np.linspace(0, self.numberOfSamples, num=21) + rng = self.startingRange + np.linspace(0, self.numberOfSamples-1, num=21) * self.rangePixelSize + deltat = np.linspace(0, 1., num=21)[1:-1] + + lonlat = osr.SpatialReference() + lonlat.ImportFromEPSG(4326) + + coord = osr.SpatialReference() + coord.ImportFromEPSG(self.epsg) + + trans = osr.CoordinateTransformation(lonlat, coord) + + llhs = [] + xyzs = [] + + + ###First range line + for rr in rng: + for zz in zrange: + llh = self.orbit.rdr2geo(self.sensingStart, rr, side=self.lookSide, height=zz) + llhs.append(llh) + if gdal.__version__[0] == '2': + x,y,z = trans.TransformPoint(llh[1], llh[0], llh[2]) + else: + x,y,z = trans.TransformPoint(llh[0], llh[1], llh[2]) + xyzs.append([x,y,z]) + + ##Last range line + sensingStop = self.sensingStart + datetime.timedelta(seconds = (self.numberOfLines-1) / self.prf) + for rr in rng: + for zz in zrange: + llh = self.orbit.rdr2geo(sensingStop, rr, side=self.lookSide, height=zz) + llhs.append(llh) + if gdal.__version__[0] == '2': + x,y,z = trans.TransformPoint(llh[1], llh[0], llh[2]) + else: + x,y,z = trans.TransformPoint(llh[0], llh[1], llh[2]) + xyzs.append([x,y,z]) + + + ##For each line in middle, consider the edges + for frac in deltat: + sensingTime = self.sensingStart + datetime.timedelta(seconds = frac * (self.numberOfLines-1)/self.prf) +# print('sensing Time: %f %f %f'%(sensingTime.minute,sensingTime.second,sensingTime.microsecond)) + for rr in [rng[0], rng[-1]]: + for zz in zrange: + llh = self.orbit.rdr2geo(sensingTime, rr, side=self.lookSide, height=zz) + llhs.append(llh) + if gdal.__version__[0] == '2': + x,y,z = trans.TransformPoint(llh[1], llh[0], llh[2]) + else: + x,y,z = trans.TransformPoint(llh[0], llh[1], llh[2]) + xyzs.append([x,y,z]) + + + llhs = np.array(llhs) + xyzs = np.array(xyzs) + + + self._xlim = [np.min(xyzs[:,0]), np.max(xyzs[:,0])] + self._ylim = [np.min(xyzs[:,1]), np.max(xyzs[:,1])] + + + def getIncidenceAngle(self, zrange=[-200,4000]): + ''' + Dummy. + ''' + import numpy as np + import datetime + from osgeo import osr,gdal + from isceobj.Util.geo.ellipsoid import Ellipsoid + from isceobj.Planet.Planet import Planet + + planet = Planet(pname='Earth') + refElp = Ellipsoid(a=planet.ellipsoid.a, e2=planet.ellipsoid.e2, model='WGS84') + + deg2rad = np.pi/180.0 + + thetas = [] + + midrng = self.startingRange + (np.floor(self.numberOfSamples/2)-1) * self.rangePixelSize + midsensing = self.sensingStart + datetime.timedelta(seconds = (np.floor(self.numberOfLines/2)-1) / self.prf) + masterSV = self.orbit.interpolateOrbit(midsensing, method='hermite') + mxyz = np.array(masterSV.getPosition()) + + for zz in zrange: + llh = self.orbit.rdr2geo(midsensing, midrng, side=self.lookSide, height=zz) + targxyz = np.array(refElp.LLH(llh[0], llh[1], llh[2]).ecef().tolist()) + los = (mxyz-targxyz) / np.linalg.norm(mxyz-targxyz) + n_vec = np.array([np.cos(llh[0]*deg2rad)*np.cos(llh[1]*deg2rad), np.cos(llh[0]*deg2rad)*np.sin(llh[1]*deg2rad), np.sin(llh[0]*deg2rad)]) + theta = np.arccos(np.dot(los, n_vec)) + thetas.append([theta]) + + thetas = np.array(thetas) + + self.incidenceAngle = np.mean(thetas) + + + def getDEM(self, bbox): + ''' + Look up database and return values. + ''' + + return "", "", "", "", "" + + def setState(self): + ''' + Create C object and populate. + ''' + from components.contrib.geo_autoRIFT.geogrid import geogrid + from iscesys import DateTimeUtil as DTU + + if self._geogrid is not None: + geogrid.destroyGeoGrid_Py(self._geogrid) + + self._geogrid = geogrid.createGeoGrid_Py() + geogrid.setRadarImageDimensions_Py( self._geogrid, self.numberOfSamples, self.numberOfLines) + geogrid.setRangeParameters_Py( self._geogrid, self.startingRange, self.rangePixelSize) + geogrid.setAzimuthParameters_Py( self._geogrid, DTU.seconds_since_midnight(self.sensingStart), self.prf) + geogrid.setRepeatTime_Py(self._geogrid, self.repeatTime) + + geogrid.setEPSG_Py(self._geogrid, self.epsg) + geogrid.setIncidenceAngle_Py(self._geogrid, self.incidenceAngle) + geogrid.setChipSizeX0_Py(self._geogrid, self.chipSizeX0) + + geogrid.setXLimits_Py(self._geogrid, self._xlim[0], self._xlim[1]) + geogrid.setYLimits_Py(self._geogrid, self._ylim[0], self._ylim[1]) + if self.demname: + geogrid.setDEM_Py(self._geogrid, self.demname) + + if (self.dhdxname is not None) and (self.dhdyname is not None): + geogrid.setSlopes_Py(self._geogrid, self.dhdxname, self.dhdyname) + + if (self.vxname is not None) and (self.vyname is not None): + geogrid.setVelocities_Py(self._geogrid, self.vxname, self.vyname) + + if (self.srxname is not None) and (self.sryname is not None): + geogrid.setSearchRange_Py(self._geogrid, self.srxname, self.sryname) + + if (self.csminxname is not None) and (self.csminyname is not None): + geogrid.setChipSizeMin_Py(self._geogrid, self.csminxname, self.csminyname) + + if (self.csmaxxname is not None) and (self.csmaxyname is not None): + geogrid.setChipSizeMax_Py(self._geogrid, self.csmaxxname, self.csmaxyname) + + if (self.ssmname is not None): + geogrid.setStableSurfaceMask_Py(self._geogrid, self.ssmname) + + geogrid.setWindowLocationsFilename_Py( self._geogrid, self.winlocname) + geogrid.setWindowOffsetsFilename_Py( self._geogrid, self.winoffname) + geogrid.setWindowSearchRangeFilename_Py( self._geogrid, self.winsrname) + geogrid.setWindowChipSizeMinFilename_Py( self._geogrid, self.wincsminname) + geogrid.setWindowChipSizeMaxFilename_Py( self._geogrid, self.wincsmaxname) + geogrid.setWindowStableSurfaceMaskFilename_Py( self._geogrid, self.winssmname) + geogrid.setRO2VXFilename_Py( self._geogrid, self.winro2vxname) + geogrid.setRO2VYFilename_Py( self._geogrid, self.winro2vyname) + geogrid.setLookSide_Py(self._geogrid, self.lookSide) + geogrid.setNodataOut_Py(self._geogrid, self.nodata_out) + + self._orbit = self.orbit.exportToC() + geogrid.setOrbit_Py(self._geogrid, self._orbit) + + def finalize(self): + ''' + Clean up all the C pointers. + ''' + + from components.contrib.geo_autoRIFT.geogrid import geogrid + from isceobj.Util import combinedlibmodule + + combinedlibmodule.freeCOrbit(self._orbit) + self._orbit = None + + geogrid.destroyGeoGrid_Py(self._geogrid) + self._geogrid = None + + def __init__(self): + super(Geogrid, self).__init__() + + ##Radar image related parameters + self.orbit = None + self.sensingStart = None + self.startingRange = None + self.prf = None + self.rangePixelSize = None + self.numberOfSamples = None + self.numberOfLines = None + self.lookSide = None + self.repeatTime = None + self.incidenceAngle = None + self.chipSizeX0 = None + + ##Input related parameters + self.demname = None + self.dhdxname = None + self.dhdyname = None + self.vxname = None + self.vyname = None + self.srxname = None + self.sryname = None + self.csminxname = None + self.csminyname = None + self.csmaxxname = None + self.csmaxyname = None + self.ssmname = None + + ##Output related parameters + self.winlocname = None + self.winoffname = None + self.winsrname = None + self.wincsminname = None + self.wincsmaxname = None + self.winssmname = None + self.winro2vxname = None + self.winro2vyname = None + + + + ##Coordinate system + self.epsg = None + self._xlim = None + self._ylim = None + self.nodata_out = None + + ##Pointer to C + self._geogrid = None + self._orbit = None diff --git a/contrib/geo_autoRIFT/geogrid/GeogridOptical.py b/contrib/geo_autoRIFT/geogrid/GeogridOptical.py new file mode 100755 index 0000000..2832532 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/GeogridOptical.py @@ -0,0 +1,757 @@ +#!/usr/bin/env python3 + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2019 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, +# end user, or in support of a prohibited end use). By downloading this software, +# the user agrees to comply with all applicable U.S. export laws and regulations. +# The user has the responsibility to obtain export licenses, or other export +# authority as may be required before exporting this software to any 'EAR99' +# embargoed foreign country or citizen of those countries. +# +# Authors: Piyush Agram, Yang Lei +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + + + + +import pdb +import subprocess +import re +import string + +class GeogridOptical(): + ''' + Class for mapping regular geographic grid on radar imagery. + ''' + + def runGeogrid(self): + ''' + Do the actual processing. + ''' + + ##Determine appropriate EPSG system + self.epsgDem = self.getProjectionSystem(self.demname) + self.epsgDat = self.getProjectionSystem(self.dat1name) + + ###Determine extent of data needed + bbox = self.determineBbox() + + + ##Run + self.geogrid() + + + def getProjectionSystem(self, filename): + ''' + Testing with Greenland. + ''' + if not filename: + raise Exception('File {0} does not exist'.format(filename)) + + from osgeo import gdal, osr + ds = gdal.Open(filename, gdal.GA_ReadOnly) + srs = osr.SpatialReference() + srs.ImportFromWkt(ds.GetProjection()) + srs.AutoIdentifyEPSG() + ds = None +# pdb.set_trace() + + if srs.IsGeographic(): + epsgstr = srs.GetAuthorityCode('GEOGCS') + elif srs.IsProjected(): + epsgstr = srs.GetAuthorityCode('PROJCS') + elif srs.IsLocal(): + raise Exception('Local coordinate system encountered') + else: + raise Exception('Non-standard coordinate system encountered') + if not epsgstr: #Empty string->use shell command gdalsrsinfo for last trial + cmd = 'gdalsrsinfo -o epsg {0}'.format(filename) + epsgstr = subprocess.check_output(cmd, shell=True) +# pdb.set_trace() + epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0] +# pdb.set_trace() + if not epsgstr: #Empty string + raise Exception('Could not auto-identify epsg code') +# pdb.set_trace() + epsgcode = int(epsgstr) +# pdb.set_trace() + return epsgcode + + def determineBbox(self, zrange=[-200,4000]): + ''' + Dummy. + ''' + import numpy as np + import datetime + from osgeo import osr + +# import pdb +# pdb.set_trace() + + + samples = self.startingX + np.array([0, self.numberOfSamples]) * self.XSize + lines = self.startingY + np.array([0, self.numberOfLines]) * self.YSize + + coordDat = osr.SpatialReference() + if self.epsgDat: + coordDat.ImportFromEPSG(self.epsgDat) + else: + raise Exception('EPSG code does not exist for image data') + + + coordDem = osr.SpatialReference() + if self.epsgDem: + coordDem.ImportFromEPSG(self.epsgDem) + else: + raise Exception('EPSG code does not exist for DEM') + + + trans = osr.CoordinateTransformation(coordDat, coordDem) + + + + utms = [] + xyzs = [] + + + ### Four corner coordinates + for ss in samples: + for ll in lines: + for zz in zrange: + utms.append([ss,ll,zz]) + x,y,z = trans.TransformPoint(ss, ll, zz) + xyzs.append([x,y,z]) + + utms = np.array(utms) + xyzs = np.array(xyzs) + + self._xlim = [np.min(xyzs[:,0]), np.max(xyzs[:,0])] + self._ylim = [np.min(xyzs[:,1]), np.max(xyzs[:,1])] + + + + + + + def geogrid(self): + + # For now print inputs that were obtained + + print("Optical Image parameters: ") + print("X-direction coordinate: " + str(self.startingX) + " " + str(self.XSize)) + print("Y-direction coordinate: " + str(self.startingY) + " " + str(self.YSize)) + print("Dimensions: " + str(self.numberOfSamples) + " " + str(self.numberOfLines)) + + print("Map inputs: ") + print("EPSG: " + str(self.epsgDem)) + print("Smallest Allowable Chip Size in m: " + str(self.chipSizeX0)) + print("Repeat Time: " + str(self.repeatTime)) + print("XLimits: " + str(self._xlim[0]) + " " + str(self._xlim[1])) + print("YLimits: " + str(self._ylim[0]) + " " + str(self._ylim[1])) + print("Extent in km: " + str((self._xlim[1]-self._xlim[0])/1000.0) + " " + str((self._ylim[1]-self._ylim[0])/1000.0)) + print("DEM: " + str(self.demname)) + print("Velocities: " + str(self.vxname) + " " + str(self.vyname)) + print("Search Range: " + str(self.srxname) + " " + str(self.sryname)) + print("Chip Size Min: " + str(self.csminxname) + " " + str(self.csminyname)) + print("Chip Size Max: " + str(self.csmaxxname) + " " + str(self.csmaxyname)) + print("Stable Surface Mask: " + str(self.ssmname)) + print("Slopes: " + str(self.dhdxname) + " " + str(self.dhdyname)) + print("Output Nodata Value: " + str(self.nodata_out)) + + print("Outputs: ") + print("Window locations: " + str(self.winlocname)) + print("Window offsets: " + str(self.winoffname)) + print("Window search range: " + str(self.winsrname)) + print("Window chip size min: " + str(self.wincsminname)) + print("Window chip size max: " + str(self.wincsmaxname)) + print("Window stable surface mask: " + str(self.winssmname)) + print("Window rdr_off2vel_x vector: " + str(self.winro2vxname)) + print("Window rdr_off2vel_y vector: " + str(self.winro2vyname)) + + + + print("Starting processing .... ") + + + + + from osgeo import gdal, osr + import numpy as np + import struct + +# pdb.set_trace() + demDS = gdal.Open(self.demname, gdal.GA_ReadOnly) + + if (self.dhdxname != ""): + sxDS = gdal.Open(self.dhdxname, gdal.GA_ReadOnly) + syDS = gdal.Open(self.dhdyname, gdal.GA_ReadOnly) + + if (self.vxname != ""): + vxDS = gdal.Open(self.vxname, gdal.GA_ReadOnly) + vyDS = gdal.Open(self.vyname, gdal.GA_ReadOnly) + + if (self.srxname != ""): + srxDS = gdal.Open(self.srxname, gdal.GA_ReadOnly) + sryDS = gdal.Open(self.sryname, gdal.GA_ReadOnly) + + if (self.csminxname != ""): + csminxDS = gdal.Open(self.csminxname, gdal.GA_ReadOnly) + csminyDS = gdal.Open(self.csminyname, gdal.GA_ReadOnly) + + if (self.csmaxxname != ""): + csmaxxDS = gdal.Open(self.csmaxxname, gdal.GA_ReadOnly) + csmaxyDS = gdal.Open(self.csmaxyname, gdal.GA_ReadOnly) + + if (self.ssmname != ""): + ssmDS = gdal.Open(self.ssmname, gdal.GA_ReadOnly) + + if demDS is None: + raise Exception('Error opening DEM file {0}'.format(self.demname)) + + if (self.dhdxname != ""): + if (sxDS is None): + raise Exception('Error opening x-direction slope file {0}'.format(self.dhdxname)) + if (syDS is None): + raise Exception('Error opening y-direction slope file {0}'.format(self.dhdyname)) + + if (self.vxname != ""): + if (vxDS is None): + raise Exception('Error opening x-direction velocity file {0}'.format(self.vxname)) + if (vyDS is None): + raise Exception('Error opening y-direction velocity file {0}'.format(self.vyname)) + + if (self.srxname != ""): + if (srxDS is None): + raise Exception('Error opening x-direction search range file {0}'.format(self.srxname)) + if (sryDS is None): + raise Exception('Error opening y-direction search range file {0}'.format(self.sryname)) + + if (self.csminxname != ""): + if (csminxDS is None): + raise Exception('Error opening x-direction chip size min file {0}'.format(self.csminxname)) + if (csminyDS is None): + raise Exception('Error opening y-direction chip size min file {0}'.format(self.csminyname)) + + if (self.csmaxxname != ""): + if (csmaxxDS is None): + raise Exception('Error opening x-direction chip size max file {0}'.format(self.csmaxxname)) + if (csmaxyDS is None): + raise Exception('Error opening y-direction chip size max file {0}'.format(self.csmaxyname)) + + if (self.ssmname != ""): + if (ssmDS is None): + raise Exception('Error opening stable surface mask file {0}'.format(self.ssmname)) + + geoTrans = demDS.GetGeoTransform() + demXSize = demDS.RasterXSize + demYSize = demDS.RasterYSize + + + # Get offsets and size to read from DEM + lOff = int(np.max( [np.floor((self._ylim[1] - geoTrans[3])/geoTrans[5]), 0.])) +# pdb.set_trace() + lCount = int(np.min([ np.ceil((self._ylim[0] - geoTrans[3])/geoTrans[5]), demYSize-1.]) - lOff) + + pOff = int(np.max([ np.floor((self._xlim[0] - geoTrans[0])/geoTrans[1]), 0.])) + pCount = int(np.min([ np.ceil((self._xlim[1] - geoTrans[0])/geoTrans[1]), demXSize-1.]) - pOff) + + print("Xlimits : " + str(geoTrans[0] + pOff * geoTrans[1]) + " " + str(geoTrans[0] + (pOff + pCount) * geoTrans[1])) + + print("Ylimits : " + str(geoTrans[3] + (lOff + lCount) * geoTrans[5]) + " " + str(geoTrans[3] + lOff * geoTrans[5])) + + print("Dimensions of geogrid: " + str(pCount) + " x " + str(lCount)) + + projDem = osr.SpatialReference() + if self.epsgDem: + projDem.ImportFromEPSG(self.epsgDem) + else: + raise Exception('EPSG code does not exist for DEM') + + projDat = osr.SpatialReference() + if self.epsgDat: + projDat.ImportFromEPSG(self.epsgDat) + else: + raise Exception('EPSG code does not exist for image data') + + fwdTrans = osr.CoordinateTransformation(projDem, projDat) + invTrans = osr.CoordinateTransformation(projDat, projDem) + + if (self.vxname != ""): + nodata = vxDS.GetRasterBand(1).GetNoDataValue() + + nodata_out = self.nodata_out + + + pszFormat = "GTiff" + adfGeoTransform = ( self._xlim[0], (self._xlim[1]-self._xlim[0])/pCount, 0, self._ylim[1], 0, (self._ylim[0]-self._ylim[1])/lCount ) + oSRS = osr.SpatialReference() + pszSRS_WKT = projDem.ExportToWkt() + + + + poDriver = gdal.GetDriverByName(pszFormat) + if( poDriver is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilename = self.winlocname + poDstDS = poDriver.Create(pszDstFilename, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) + poDstDS.SetGeoTransform( adfGeoTransform ) + poDstDS.SetProjection( pszSRS_WKT ) + + poBand1 = poDstDS.GetRasterBand(1) + poBand2 = poDstDS.GetRasterBand(2) + poBand1.SetNoDataValue(nodata_out) + poBand2.SetNoDataValue(nodata_out) + + + + + poDriverOff = gdal.GetDriverByName(pszFormat) + if( poDriverOff is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameOff = self.winoffname + poDstDSOff = poDriverOff.Create(pszDstFilenameOff, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) + poDstDSOff.SetGeoTransform( adfGeoTransform ) + poDstDSOff.SetProjection( pszSRS_WKT ) + + poBand1Off = poDstDSOff.GetRasterBand(1) + poBand2Off = poDstDSOff.GetRasterBand(2) + poBand1Off.SetNoDataValue(nodata_out) + poBand2Off.SetNoDataValue(nodata_out) + + + + poDriverSch = gdal.GetDriverByName(pszFormat) + if( poDriverSch is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameSch = self.winsrname + poDstDSSch = poDriverSch.Create(pszDstFilenameSch, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) + poDstDSSch.SetGeoTransform( adfGeoTransform ) + poDstDSSch.SetProjection( pszSRS_WKT ) + + poBand1Sch = poDstDSSch.GetRasterBand(1) + poBand2Sch = poDstDSSch.GetRasterBand(2) + poBand1Sch.SetNoDataValue(nodata_out) + poBand2Sch.SetNoDataValue(nodata_out) + + + poDriverMin = gdal.GetDriverByName(pszFormat) + if( poDriverMin is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameMin = self.wincsminname + poDstDSMin = poDriverMin.Create(pszDstFilenameMin, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) + poDstDSMin.SetGeoTransform( adfGeoTransform ) + poDstDSMin.SetProjection( pszSRS_WKT ) + + poBand1Min = poDstDSMin.GetRasterBand(1) + poBand2Min = poDstDSMin.GetRasterBand(2) + poBand1Min.SetNoDataValue(nodata_out) + poBand2Min.SetNoDataValue(nodata_out) + + + poDriverMax = gdal.GetDriverByName(pszFormat) + if( poDriverMax is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameMax = self.wincsmaxname + poDstDSMax = poDriverMax.Create(pszDstFilenameMax, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) + poDstDSMax.SetGeoTransform( adfGeoTransform ) + poDstDSMax.SetProjection( pszSRS_WKT ) + + poBand1Max = poDstDSMax.GetRasterBand(1) + poBand2Max = poDstDSMax.GetRasterBand(2) + poBand1Max.SetNoDataValue(nodata_out) + poBand2Max.SetNoDataValue(nodata_out) + + + + poDriverMsk = gdal.GetDriverByName(pszFormat) + if( poDriverMsk is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameMsk = self.winssmname + poDstDSMsk = poDriverMsk.Create(pszDstFilenameMsk, xsize=pCount, ysize=lCount, bands=1, eType=gdal.GDT_Int32) + poDstDSMsk.SetGeoTransform( adfGeoTransform ) + poDstDSMsk.SetProjection( pszSRS_WKT ) + + poBand1Msk = poDstDSMsk.GetRasterBand(1) + poBand1Msk.SetNoDataValue(nodata_out) + + + + + + poDriverRO2VX = gdal.GetDriverByName(pszFormat) + if( poDriverRO2VX is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameRO2VX = self.winro2vxname + poDstDSRO2VX = poDriverRO2VX.Create(pszDstFilenameRO2VX, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Float64) + poDstDSRO2VX.SetGeoTransform( adfGeoTransform ) + poDstDSRO2VX.SetProjection( pszSRS_WKT ) + + poBand1RO2VX = poDstDSRO2VX.GetRasterBand(1) + poBand2RO2VX = poDstDSRO2VX.GetRasterBand(2) + poBand1RO2VX.SetNoDataValue(nodata_out) + poBand2RO2VX.SetNoDataValue(nodata_out) + + + + + poDriverRO2VY = gdal.GetDriverByName(pszFormat) + if( poDriverRO2VY is None ): + raise Exception('Cannot create gdal driver for output') + + pszDstFilenameRO2VY = self.winro2vyname + poDstDSRO2VY = poDriverRO2VY.Create(pszDstFilenameRO2VY, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Float64) + poDstDSRO2VY.SetGeoTransform( adfGeoTransform ) + poDstDSRO2VY.SetProjection( pszSRS_WKT ) + + poBand1RO2VY = poDstDSRO2VY.GetRasterBand(1) + poBand2RO2VY = poDstDSRO2VY.GetRasterBand(2) + poBand1RO2VY.SetNoDataValue(nodata_out) + poBand2RO2VY.SetNoDataValue(nodata_out) + + + + raster1 = np.zeros(pCount,dtype=np.int32) + raster2 = np.zeros(pCount,dtype=np.int32) + raster11 = np.zeros(pCount,dtype=np.int32) + raster22 = np.zeros(pCount,dtype=np.int32) + sr_raster11 = np.zeros(pCount,dtype=np.int32) + sr_raster22 = np.zeros(pCount,dtype=np.int32) + csmin_raster11 = np.zeros(pCount,dtype=np.int32) + csmin_raster22 = np.zeros(pCount,dtype=np.int32) + csmax_raster11 = np.zeros(pCount,dtype=np.int32) + csmax_raster22 = np.zeros(pCount,dtype=np.int32) + ssm_raster = np.zeros(pCount,dtype=np.int32) + raster1a = np.zeros(pCount,dtype=np.float64) + raster1b = np.zeros(pCount,dtype=np.float64) + raster2a = np.zeros(pCount,dtype=np.float64) + raster2b = np.zeros(pCount,dtype=np.float64) + + + + # X- and Y-direction pixel size + X_res = np.abs(self.XSize) + Y_res = np.abs(self.YSize) + print("X-direction pixel size: " + str(X_res)) + print("Y-direction pixel size: " + str(Y_res)) + + ChipSizeX0_PIX_X = np.ceil(self.chipSizeX0 / X_res / 4) * 4 + ChipSizeX0_PIX_Y = np.ceil(self.chipSizeX0 / Y_res / 4) * 4 + + + + + + for ii in range(lCount): + y = geoTrans[3] + (lOff+ii+0.5) * geoTrans[5] + demLine = demDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + demLine = struct.unpack('d' * pCount, demLine) + + if (self.dhdxname != ""): + sxLine = sxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + sxLine = struct.unpack('d' * pCount, sxLine) + syLine = syDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + syLine = struct.unpack('d' * pCount, syLine) + + if (self.vxname != ""): + vxLine = vxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + vxLine = struct.unpack('d' * pCount, vxLine) + vyLine = vyDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + vyLine = struct.unpack('d' * pCount, vyLine) + + if (self.srxname != ""): + srxLine = srxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + srxLine = struct.unpack('d' * pCount, srxLine) + sryLine = sryDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + sryLine = struct.unpack('d' * pCount, sryLine) + + if (self.csminxname != ""): + csminxLine = csminxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + csminxLine = struct.unpack('d' * pCount, csminxLine) + csminyLine = csminyDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + csminyLine = struct.unpack('d' * pCount, csminyLine) + + if (self.csmaxxname != ""): + csmaxxLine = csmaxxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + csmaxxLine = struct.unpack('d' * pCount, csmaxxLine) + csmaxyLine = csmaxyDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + csmaxyLine = struct.unpack('d' * pCount, csmaxyLine) + + if (self.ssmname != ""): + ssmLine = ssmDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + ssmLine = struct.unpack('d' * pCount, ssmLine) + + for jj in range(pCount): + xyzs = np.array([geoTrans[0] + (jj+pOff+0.5)*geoTrans[1], y, demLine[jj]]) + targxyz0 = xyzs.copy() + if (self.dhdxname != ""): + slp = np.array([sxLine[jj], syLine[jj], -1.0]) + if (self.vxname != ""): + vel = np.array([vxLine[jj], vyLine[jj], 0.0]) + if (self.srxname != ""): + schrng1 = np.array([srxLine[jj], sryLine[jj], 0.0]) + schrng2 = np.array([-srxLine[jj], sryLine[jj], 0.0]) + targutm0 = np.array(fwdTrans.TransformPoint(targxyz0[0],targxyz0[1],targxyz0[2])) + xind = np.round((targutm0[0] - self.startingX) / self.XSize) + 1. + yind = np.round((targutm0[1] - self.startingY) / self.YSize) + 1. + + # x-direction vector + targutm = targutm0.copy() + targutm[0] = targutm0[0] + self.XSize + targxyz = np.array(invTrans.TransformPoint(targutm[0],targutm[1],targutm[2])) + xunit = (targxyz-targxyz0) / np.linalg.norm(targxyz-targxyz0) + + # y-direction vector + targutm = targutm0.copy() + targutm[1] = targutm0[1] + self.YSize + targxyz = np.array(invTrans.TransformPoint(targutm[0],targutm[1],targutm[2])) + yunit = (targxyz-targxyz0) / np.linalg.norm(targxyz-targxyz0) + + # local normal vector + if (self.dhdxname != ""): + normal = -slp / np.linalg.norm(slp) + else: + normal = np.array([0., 0., 0.]) + + if (self.vxname != ""): + vel[2] = -(vel[0]*normal[0]+vel[1]*normal[1])/normal[2] + + if (self.srxname != ""): + schrng1[2] = -(schrng1[0]*normal[0]+schrng1[1]*normal[1])/normal[2] + schrng2[2] = -(schrng2[0]*normal[0]+schrng2[1]*normal[1])/normal[2] + + + + if ((xind > self.numberOfSamples)|(xind < 1)|(yind > self.numberOfLines)|(yind < 1)): +# pdb.set_trace() + raster1[jj] = nodata_out + raster2[jj] = nodata_out + raster11[jj] = nodata_out + raster22[jj] = nodata_out + + sr_raster11[jj] = nodata_out + sr_raster22[jj] = nodata_out + csmin_raster11[jj] = nodata_out + csmin_raster22[jj] = nodata_out + csmax_raster11[jj] = nodata_out + csmax_raster22[jj] = nodata_out + ssm_raster[jj] = nodata_out + + raster1a[jj] = nodata_out + raster1b[jj] = nodata_out + raster2a[jj] = nodata_out + raster2b[jj] = nodata_out + else: + raster1[jj] = xind; + raster2[jj] = yind; +# pdb.set_trace() +# if ((self.vxname != "")&(vel[0] != nodata)): +## pdb.set_trace() +# raster11[jj] = np.round(np.dot(vel,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1) +# raster22[jj] = np.round(np.dot(vel,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1) +# else: +# raster11[jj] = 0. +# raster22[jj] = 0. + if (self.vxname == ""): +# pdb.set_trace() + raster11[jj] = 0. + raster22[jj] = 0. + elif (vel[0] == nodata): + raster11[jj] = 0. + raster22[jj] = 0. + else: + raster11[jj] = np.round(np.dot(vel,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1) + raster22[jj] = np.round(np.dot(vel,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1) + + + if (self.srxname == ""): + sr_raster11[jj] = nodata_out + sr_raster22[jj] = nodata_out + elif (vel[0] == nodata): + sr_raster11[jj] = 0 + sr_raster22[jj] = 0 + else: + sr_raster11[jj] = np.abs(np.round(np.dot(schrng1,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)) + sr_raster22[jj] = np.abs(np.round(np.dot(schrng1,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)) + if (np.abs(np.round(np.dot(schrng2,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)) > sr_raster11[jj]): + sr_raster11[jj] = np.abs(np.round(np.dot(schrng2,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)) + if (np.abs(np.round(np.dot(schrng2,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)) > sr_raster22[jj]): + sr_raster22[jj] = np.abs(np.round(np.dot(schrng2,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)) + if (sr_raster11[jj] == 0): + sr_raster11[jj] = 1 + if (sr_raster22[jj] == 0): + sr_raster22[jj] = 1 + + if (self.csminxname != ""): + csmin_raster11[jj] = csminxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X + csmin_raster22[jj] = csminyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y + else: + csmin_raster11[jj] = nodata_out + csmin_raster22[jj] = nodata_out + + if (self.csmaxxname != ""): + csmax_raster11[jj] = csmaxxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X + csmax_raster22[jj] = csmaxyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y + else: + csmax_raster11[jj] = nodata_out + csmax_raster22[jj] = nodata_out + + + if (self.ssmname != ""): + ssm_raster[jj] = ssmLine[jj] + else: + ssm_raster[jj] = nodata_out + + + + raster1a[jj] = normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[1]-normal[1]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + raster1b[jj] = -normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[1]-normal[1]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + raster2a[jj] = -normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[0]-normal[0]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + raster2b[jj] = normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[0]-normal[0]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + + +# pdb.set_trace() + + poBand1.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand2.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand1Off.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand2Off.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand1Sch.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=sr_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand2Sch.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=sr_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand1Min.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmin_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand2Min.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmin_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand1Max.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmax_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand2Max.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmax_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand1Msk.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=ssm_raster.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) + poBand1RO2VX.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1a.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + poBand2RO2VX.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1b.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + poBand1RO2VY.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2a.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + poBand2RO2VY.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2b.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) + + + poDstDS = None + + poDstDSOff = None + + poDstDSSch = None + + poDstDSMin = None + + poDstDSMax = None + + poDstDSMsk = None + + poDstDSRO2VX = None + + poDstDSRO2VY = None + + demDS = None + + if (self.dhdxname != ""): + sxDS = None + syDS = None + + if (self.vxname != ""): + vxDS = None + vyDS = None + + if (self.srxname != ""): + srxDS = None + sryDS = None + + if (self.csminxname != ""): + csminxDS = None + csminyDS = None + + if (self.csmaxxname != ""): + csmaxxDS = None + csmaxyDS = None + + if (self.ssmname != ""): + ssmDS = None + + + + + + + + + + + + + + + + + def __init__(self): + super(GeogridOptical, self).__init__() + + ##Optical image related parameters + self.startingY = None + self.startingX = None + self.XSize = None + self.YSize = None + self.numberOfSamples = None + self.numberOfLines = None + self.repeatTime = None + self.chipSizeX0 = None + + ##Input related parameters + self.dat1name = None + self.demname = None + self.dhdxname = None + self.dhdyname = None + self.vxname = None + self.vyname = None + self.srxname = None + self.sryname = None + self.csminxname = None + self.csminyname = None + self.csmaxxname = None + self.csmaxyname = None + self.ssmname = None + + ##Output related parameters + self.winlocname = None + self.winoffname = None + self.winsrname = None + self.wincsminname = None + self.wincsmaxname = None + self.winssmname = None + self.winro2vxname = None + self.winro2vyname = None + + ##Coordinate system + self.epsgDem = None + self.epsgDat = None + self._xlim = None + self._ylim = None + self.nodata_out = None + + diff --git a/contrib/geo_autoRIFT/geogrid/SConscript b/contrib/geo_autoRIFT/geogrid/SConscript new file mode 100755 index 0000000..c7128ad --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/SConscript @@ -0,0 +1,38 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Authors: Piyush Agram, Yang Lei +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeoAutorift') +envgeogrid = envgeoAutorift.Clone() +package = envgeogrid['PACKAGE'] +project = 'geogrid' +envgeogrid['PROJECT'] = project +install = envgeogrid['PRJ_SCONS_INSTALL'] + '/' + package + '/' + project +initFile = '__init__.py' +if not os.path.exists(initFile): + fout = open(initFile,"w") + fout.write("#!/usr/bin/env python3") + fout.close() + +listFiles = ['Geogrid.py','GeogridOptical.py',initFile] +envgeogrid.Install(install,listFiles) +envgeogrid.Alias('install',install) +Export('envgeogrid') +bindingsScons = 'bindings/SConscript' +SConscript(bindingsScons,variant_dir = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/bindings') +includeScons = 'include/SConscript' +SConscript(includeScons) +srcScons = 'src/SConscript' +SConscript(srcScons,variant_dir = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/src') diff --git a/contrib/geo_autoRIFT/geogrid/__init__.py b/contrib/geo_autoRIFT/geogrid/__init__.py new file mode 100755 index 0000000..58dec40 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/__init__.py @@ -0,0 +1,14 @@ +#!/usr/bin/env python + +#def SplitRangeSpectrum(): +# from .splitSpectrum import PySplitRangeSpectrum +# return PySplitRangeSpectrum() + + +from .GeogridOptical import GeogridOptical + +try: + from .Geogrid import Geogrid +except ImportError: + # this means ISCE support not available. Don't raise error. Allow standalone use + pass diff --git a/contrib/geo_autoRIFT/geogrid/bindings/SConscript b/contrib/geo_autoRIFT/geogrid/bindings/SConscript new file mode 100755 index 0000000..3ae898c --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/bindings/SConscript @@ -0,0 +1,28 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeogrid') +package = envgeogrid['PACKAGE'] +project = envgeogrid['PROJECT'] +install = envgeogrid['PRJ_SCONS_INSTALL'] + '/' + package + '/' + project +build = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project +libList = ['gomp','geogrid','combinedLib','gdal'] +envgeogrid.PrependUnique(LIBS = libList) +module = envgeogrid.LoadableModule(target = 'geogrid.abi3.so', source = 'geogridmodule.cpp') +envgeogrid.Install(install,module) +envgeogrid.Alias('install',install) +envgeogrid.Install(build,module) +envgeogrid.Alias('build',build) diff --git a/contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp b/contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp new file mode 100755 index 0000000..f005f3e --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp @@ -0,0 +1,455 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 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, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + + + +#include +#include +#include "geogrid.h" +#include "geogridmodule.h" + +static const char * const __doc__ = "Python extension for geogrid"; + +PyModuleDef moduledef = { + //header + PyModuleDef_HEAD_INIT, + //name of the module + "geogrid", + //module documentation string + __doc__, + //size of the per-interpreter state of the module; + -1, + geogrid_methods, +}; + +//Initialization function for the module +PyMODINIT_FUNC +PyInit_geogrid() +{ + PyObject* module = PyModule_Create(&moduledef); + if (!module) + { + return module; + } + return module; +} + +PyObject* createGeoGrid(PyObject* self, PyObject *args) +{ + geoGrid* ptr = new geoGrid; + return Py_BuildValue("K", (uint64_t) ptr); +} + +PyObject* destroyGeoGrid(PyObject *self, PyObject *args) +{ + uint64_t ptr; + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + if (((geoGrid*)(ptr))!=NULL) + { + delete ((geoGrid*)(ptr)); + } + return Py_BuildValue("i", 0); +} + +PyObject* setEPSG(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int code; + if (!PyArg_ParseTuple(args, "Ki", &ptr, &code)) + { + return NULL; + } + ((geoGrid*)(ptr))->epsgcode = code; + return Py_BuildValue("i", 0); +} + +PyObject* setIncidenceAngle(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double incidenceAngle; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &incidenceAngle)) + { + return NULL; + } + ((geoGrid*)(ptr))->incidenceAngle = incidenceAngle; + return Py_BuildValue("i", 0); +} + +PyObject* setChipSizeX0(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double chipSizeX0; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &chipSizeX0)) + { + return NULL; + } + ((geoGrid*)(ptr))->chipSizeX0 = chipSizeX0; + return Py_BuildValue("i", 0); +} + +PyObject* setRepeatTime(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double repeatTime; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &repeatTime)) + { + return NULL; + } + ((geoGrid*)(ptr))->dt = repeatTime; + return Py_BuildValue("i", 0); +} + +PyObject* setRadarImageDimensions(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int wid, len; + if (!PyArg_ParseTuple(args, "Kii", &ptr, &wid, &len)) + { + return NULL; + } + + ((geoGrid*)(ptr))->nPixels = wid; + ((geoGrid*)(ptr))->nLines = len; + return Py_BuildValue("i", 0); +} + +PyObject* setRangeParameters(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double r0, rspace; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &r0, &rspace)) + { + return NULL; + } + + ((geoGrid*)(ptr))->startingRange = r0; + ((geoGrid*)(ptr))->dr = rspace; + return Py_BuildValue("i", 0); +} + +PyObject* setAzimuthParameters(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double t0, prf; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &t0, &prf)) + { + return NULL; + } + + ((geoGrid*)(ptr))->sensingStart = t0; + ((geoGrid*)(ptr))->prf = prf; + return Py_BuildValue("i", 0); +} + +PyObject* setXLimits(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double x0, x1; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &x0, &x1)) + { + return NULL; + } + + ((geoGrid*)(ptr))->xmin = x0; + ((geoGrid*)(ptr))->xmax = x1; + return Py_BuildValue("i", 0); +} + +PyObject* setYLimits(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double x0, x1; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &x0, &x1)) + { + return NULL; + } + + ((geoGrid*)(ptr))->ymin = x0; + ((geoGrid*)(ptr))->ymax = x1; + return Py_BuildValue("i", 0); +} + +PyObject* setDEM(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->demname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setVelocities(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *vx; + char *vy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &vx, &vy)) + { + return NULL; + } + + ((geoGrid*)(ptr))->vxname = std::string(vx); + ((geoGrid*)(ptr))->vyname = std::string(vy); + return Py_BuildValue("i", 0); +} + +PyObject* setSearchRange(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *srx; + char *sry; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &srx, &sry)) + { + return NULL; + } + + ((geoGrid*)(ptr))->srxname = std::string(srx); + ((geoGrid*)(ptr))->sryname = std::string(sry); + return Py_BuildValue("i", 0); +} + +PyObject* setChipSizeMin(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *csminx; + char *csminy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &csminx, &csminy)) + { + return NULL; + } + + ((geoGrid*)(ptr))->csminxname = std::string(csminx); + ((geoGrid*)(ptr))->csminyname = std::string(csminy); + return Py_BuildValue("i", 0); +} + +PyObject* setChipSizeMax(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *csmaxx; + char *csmaxy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &csmaxx, &csmaxy)) + { + return NULL; + } + + ((geoGrid*)(ptr))->csmaxxname = std::string(csmaxx); + ((geoGrid*)(ptr))->csmaxyname = std::string(csmaxy); + return Py_BuildValue("i", 0); +} + +PyObject* setStableSurfaceMask(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *ssm; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &ssm)) + { + return NULL; + } + + ((geoGrid*)(ptr))->ssmname = std::string(ssm); + return Py_BuildValue("i", 0); +} + +PyObject* setSlopes(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *sx; + char *sy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &sx, &sy)) + { + return NULL; + } + + ((geoGrid*)(ptr))->dhdxname = std::string(sx); + ((geoGrid*)(ptr))->dhdyname = std::string(sy); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowLocationsFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->pixlinename = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowOffsetsFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->offsetname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowSearchRangeFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->searchrangename = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowChipSizeMinFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->chipsizeminname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowChipSizeMaxFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->chipsizemaxname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowStableSurfaceMaskFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->stablesurfacemaskname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setRO2VXFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->ro2vx_name = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setRO2VYFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGrid*)(ptr))->ro2vy_name = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setLookSide(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int side; + if (!PyArg_ParseTuple(args, "Ki", &ptr, &side)) + { + return NULL; + } + + ((geoGrid*)(ptr))->lookSide = side; + return Py_BuildValue("i", 0); +} + +PyObject* setNodataOut(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int nodata; + if (!PyArg_ParseTuple(args, "Ki", &ptr, &nodata)) + { + return NULL; + } + + ((geoGrid*)(ptr))->nodata_out = nodata; + return Py_BuildValue("i", 0); +} + +PyObject* setOrbit(PyObject *self, PyObject *args) +{ + uint64_t ptr, orb; + if (!PyArg_ParseTuple(args, "KK", &ptr, &orb)) + { + return NULL; + } + + ((geoGrid*)(ptr))->orbit = (cOrbit*)(orb); + return Py_BuildValue("i", 0); +} + + +PyObject* geogrid(PyObject* self, PyObject* args) +{ + uint64_t ptr; + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + ((geoGrid*)(ptr))->geogrid(); + return Py_BuildValue("i", 0); +} diff --git a/contrib/geo_autoRIFT/geogrid/include/SConscript b/contrib/geo_autoRIFT/geogrid/include/SConscript new file mode 100755 index 0000000..2e45367 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/include/SConscript @@ -0,0 +1,24 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeogrid') +package = envgeogrid['PACKAGE'] +project = envgeogrid['PROJECT'] +build = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/include' +envgeogrid.AppendUnique(CPPPATH = [build]) +listFiles = ['geogrid.h', 'geogridmodule.h'] +envgeogrid.Install(build,listFiles) +envgeogrid.Alias('build',build) diff --git a/contrib/geo_autoRIFT/geogrid/include/geogrid.h b/contrib/geo_autoRIFT/geogrid/include/geogrid.h new file mode 100755 index 0000000..966df6c --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/include/geogrid.h @@ -0,0 +1,107 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 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, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + +#ifndef GEOGRID_H +#define GEOGRID_H + +#include + +extern "C" +{ +#include "orbit.h" +}; + +struct geoGrid +{ + //DEM related inputs + std::string demname; //DEM + std::string dhdxname; //Slope in X + std::string dhdyname; //Slope in Y + std::string vxname; //Velocity in X + std::string vyname; //Velocity in Y + std::string srxname; //Search range in X + std::string sryname; //Search range in Y + std::string csminxname; //Chip size minimum in x + std::string csminyname; //Chip size minimum in y + std::string csmaxxname; //Chip size maximum in x + std::string csmaxyname; //Chip size maximum in y + std::string ssmname; //Stable surface mask + int epsgcode; + double chipSizeX0; + + //Bounding box related + double xmin, xmax; + double ymin, ymax; + + //Radar image related inputs + cOrbit *orbit; + double sensingStart; + double prf; + int nLines; + double startingRange; + double dr; + double dt; + int nPixels; + int lookSide; + int nodata_out; + double incidenceAngle; + + //Output file names + std::string pixlinename; + std::string offsetname; + std::string searchrangename; + std::string chipsizeminname; + std::string chipsizemaxname; + std::string stablesurfacemaskname; + std::string ro2vx_name; + std::string ro2vy_name; + + //Functions + void computeBbox(double *); + void geogrid(); +}; + +struct geoGridPoint +{ + //Map coordinates + double pos[3]; + + //DEM slope info + double slope[2]; + + //Velocity related info + double vel[3]; + double schrng[3]; + + //Outputs + double range; + double aztime; + double los[3]; +}; + +#endif diff --git a/contrib/geo_autoRIFT/geogrid/include/geogridmodule.h b/contrib/geo_autoRIFT/geogrid/include/geogridmodule.h new file mode 100755 index 0000000..22bb6d0 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/include/geogridmodule.h @@ -0,0 +1,107 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 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, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + + +#ifndef geogridmodule_h +#define geogridmodule_h + +#include +#include + +extern "C" +{ + PyObject * createGeoGrid(PyObject*, PyObject*); + PyObject * destroyGeoGrid(PyObject*, PyObject*); + PyObject * geogrid(PyObject *, PyObject *); + PyObject * setRadarImageDimensions(PyObject *, PyObject *); + PyObject * setRangeParameters(PyObject *, PyObject *); + PyObject * setAzimuthParameters(PyObject*, PyObject *); + PyObject * setRepeatTime(PyObject *, PyObject *); + + PyObject * setDEM(PyObject *, PyObject *); + PyObject * setVelocities(PyObject*, PyObject*); + PyObject * setSearchRange(PyObject*, PyObject*); + PyObject * setChipSizeMin(PyObject*, PyObject*); + PyObject * setChipSizeMax(PyObject*, PyObject*); + PyObject * setStableSurfaceMask(PyObject*, PyObject*); + PyObject * setSlopes(PyObject*, PyObject*); + PyObject * setOrbit(PyObject *, PyObject *); + PyObject * setLookSide(PyObject *, PyObject *); + PyObject * setNodataOut(PyObject *, PyObject *); + + PyObject * setWindowLocationsFilename(PyObject *, PyObject *); + PyObject * setWindowOffsetsFilename(PyObject *, PyObject *); + PyObject * setWindowSearchRangeFilename(PyObject *, PyObject *); + PyObject * setWindowChipSizeMinFilename(PyObject *, PyObject *); + PyObject * setWindowChipSizeMaxFilename(PyObject *, PyObject *); + PyObject * setWindowStableSurfaceMaskFilename(PyObject *, PyObject *); + PyObject * setRO2VXFilename(PyObject *, PyObject *); + PyObject * setRO2VYFilename(PyObject *, PyObject *); + PyObject * setEPSG(PyObject *, PyObject *); + PyObject * setIncidenceAngle(PyObject *, PyObject *); + PyObject * setChipSizeX0(PyObject *, PyObject *); + PyObject * setXLimits(PyObject *, PyObject *); + PyObject * setYLimits(PyObject *, PyObject *); +} + +static PyMethodDef geogrid_methods[] = +{ + {"createGeoGrid_Py", createGeoGrid, METH_VARARGS, " "}, + {"destroyGeoGrid_Py", destroyGeoGrid, METH_VARARGS, " "}, + {"geogrid_Py", geogrid, METH_VARARGS, " "}, + {"setRadarImageDimensions_Py", setRadarImageDimensions, METH_VARARGS, " "}, + {"setRangeParameters_Py", setRangeParameters, METH_VARARGS, " "}, + {"setAzimuthParameters_Py", setAzimuthParameters, METH_VARARGS, " "}, + {"setRepeatTime_Py", setRepeatTime, METH_VARARGS, " "}, + {"setDEM_Py", setDEM, METH_VARARGS, " "}, + {"setEPSG_Py", setEPSG, METH_VARARGS, " "}, + {"setIncidenceAngle_Py", setIncidenceAngle, METH_VARARGS, " "}, + {"setChipSizeX0_Py", setChipSizeX0, METH_VARARGS, " "}, + {"setVelocities_Py", setVelocities, METH_VARARGS, " "}, + {"setSearchRange_Py", setSearchRange, METH_VARARGS, " "}, + {"setChipSizeMin_Py", setChipSizeMin, METH_VARARGS, " "}, + {"setChipSizeMax_Py", setChipSizeMax, METH_VARARGS, " "}, + {"setStableSurfaceMask_Py", setStableSurfaceMask, METH_VARARGS, " "}, + {"setSlopes_Py", setSlopes, METH_VARARGS, " "}, + {"setOrbit_Py", setOrbit, METH_VARARGS, " "}, + {"setLookSide_Py", setLookSide, METH_VARARGS, " "}, + {"setNodataOut_Py", setNodataOut, METH_VARARGS, " "}, + {"setXLimits_Py", setXLimits, METH_VARARGS, " "}, + {"setYLimits_Py", setYLimits, METH_VARARGS, " "}, + {"setWindowLocationsFilename_Py", setWindowLocationsFilename, METH_VARARGS, " "}, + {"setWindowOffsetsFilename_Py", setWindowOffsetsFilename, METH_VARARGS, " "}, + {"setWindowSearchRangeFilename_Py", setWindowSearchRangeFilename, METH_VARARGS, " "}, + {"setWindowChipSizeMinFilename_Py", setWindowChipSizeMinFilename, METH_VARARGS, " "}, + {"setWindowChipSizeMaxFilename_Py", setWindowChipSizeMaxFilename, METH_VARARGS, " "}, + {"setWindowStableSurfaceMaskFilename_Py", setWindowStableSurfaceMaskFilename, METH_VARARGS, " "}, + {"setRO2VXFilename_Py", setRO2VXFilename, METH_VARARGS, " "}, + {"setRO2VYFilename_Py", setRO2VYFilename, METH_VARARGS, " "}, + {NULL, NULL, 0, NULL} +}; +#endif //geoGridmodule_h + diff --git a/contrib/geo_autoRIFT/geogrid/src/SConscript b/contrib/geo_autoRIFT/geogrid/src/SConscript new file mode 100755 index 0000000..5695142 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/src/SConscript @@ -0,0 +1,24 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeogrid') +build = envgeogrid['PRJ_LIB_DIR'] +envgeogrid.AppendUnique(CXXFLAGS = '-fopenmp') +envgeogrid.AppendUnique(CXXFLAGS = '-std=c++11') +listFiles = ['geogrid.cpp'] +lib = envgeogrid.Library(target = 'geogrid', source = listFiles) +envgeogrid.Install(build,lib) +envgeogrid.Alias('build',build) diff --git a/contrib/geo_autoRIFT/geogrid/src/geogrid.cpp b/contrib/geo_autoRIFT/geogrid/src/geogrid.cpp new file mode 100755 index 0000000..578d313 --- /dev/null +++ b/contrib/geo_autoRIFT/geogrid/src/geogrid.cpp @@ -0,0 +1,1325 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 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, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + +#include "geogrid.h" +#include +#include +#include +#include +#include + + +extern "C" +{ +#include "linalg3.h" +#include "geometry.h" +} + +void geoGrid::geogrid() +{ + //Some constants + double deg2rad = M_PI/180.0; + + //For now print inputs that were obtained + + std::cout << "\nRadar parameters: \n"; + std::cout << "Range: " << startingRange << " " << dr << "\n"; + std::cout << "Azimuth: " << sensingStart << " " << prf << "\n"; + std::cout << "Dimensions: " << nPixels << " " << nLines << "\n"; + std::cout << "Incidence Angle: " << incidenceAngle/deg2rad << "\n"; + + std::cout << "\nMap inputs: \n"; + std::cout << "EPSG: " << epsgcode << "\n"; + std::cout << "Smallest Allowable Chip Size in m: " << chipSizeX0 << "\n"; + std::cout << "Repeat Time: " << dt << "\n"; + std::cout << "XLimits: " << xmin << " " << xmax << "\n"; + std::cout << "YLimits: " << ymin << " " << ymax << "\n"; + std::cout << "Extent in km: " << (xmax - xmin)/1000. << " " << (ymax - ymin)/1000. << "\n"; + std::cout << "DEM: " << demname << "\n"; + std::cout << "Velocities: " << vxname << " " << vyname << "\n"; + std::cout << "Search Range: " << srxname << " " << sryname << "\n"; + std::cout << "Chip Size Min: " << csminxname << " " << csminyname << "\n"; + std::cout << "Chip Size Max: " << csmaxxname << " " << csmaxyname << "\n"; + std::cout << "Slopes: " << dhdxname << " " << dhdyname << "\n"; + std::cout << "Stable Surface Mask: " << ssmname << "\n"; + std::cout << "Output Nodata Value: " << nodata_out << "\n"; + + std::cout << "\nOutputs: \n"; + std::cout << "Window locations: " << pixlinename << "\n"; + std::cout << "Window offsets: " << offsetname << "\n"; + std::cout << "Window search range: " << searchrangename << "\n"; + std::cout << "Window chip size min: " << chipsizeminname << "\n"; + std::cout << "Window chip size max: " << chipsizemaxname << "\n"; + std::cout << "Window rdr_off2vel_x vector: " << ro2vx_name << "\n"; + std::cout << "Window rdr_off2vel_y vector: " << ro2vy_name << "\n"; + std::cout << "Window stable surface mask: " << stablesurfacemaskname << "\n"; + + + std::cout << "\nStarting processing .... \n"; + + //Startup GDAL + GDALAllRegister(); + + //DEM related information + GDALDataset* demDS = NULL; + GDALDataset* sxDS = NULL; + GDALDataset* syDS = NULL; + GDALDataset* vxDS = NULL; + GDALDataset* vyDS = NULL; + GDALDataset* srxDS = NULL; + GDALDataset* sryDS = NULL; + GDALDataset* csminxDS = NULL; + GDALDataset* csminyDS = NULL; + GDALDataset* csmaxxDS = NULL; + GDALDataset* csmaxyDS = NULL; + GDALDataset* ssmDS = NULL; + + double geoTrans[6]; + + demDS = reinterpret_cast(GDALOpenShared(demname.c_str(), GA_ReadOnly)); + if (dhdxname != "") + { + sxDS = reinterpret_cast(GDALOpenShared(dhdxname.c_str(), GA_ReadOnly)); + syDS = reinterpret_cast(GDALOpenShared(dhdyname.c_str(), GA_ReadOnly)); + } + if (vxname != "") + { + vxDS = reinterpret_cast(GDALOpenShared(vxname.c_str(), GA_ReadOnly)); + vyDS = reinterpret_cast(GDALOpenShared(vyname.c_str(), GA_ReadOnly)); + } + if (srxname != "") + { + srxDS = reinterpret_cast(GDALOpenShared(srxname.c_str(), GA_ReadOnly)); + sryDS = reinterpret_cast(GDALOpenShared(sryname.c_str(), GA_ReadOnly)); + } + if (csminxname != "") + { + csminxDS = reinterpret_cast(GDALOpenShared(csminxname.c_str(), GA_ReadOnly)); + csminyDS = reinterpret_cast(GDALOpenShared(csminyname.c_str(), GA_ReadOnly)); + } + if (csmaxxname != "") + { + csmaxxDS = reinterpret_cast(GDALOpenShared(csmaxxname.c_str(), GA_ReadOnly)); + csmaxyDS = reinterpret_cast(GDALOpenShared(csmaxyname.c_str(), GA_ReadOnly)); + } + if (ssmname != "") + { + ssmDS = reinterpret_cast(GDALOpenShared(ssmname.c_str(), GA_ReadOnly)); + } + if (demDS == NULL) + { + std::cout << "Error opening DEM file { " << demname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (dhdxname != "") + { + if (sxDS == NULL) + { + std::cout << "Error opening x-direction slope file { " << dhdxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (syDS == NULL) + { + std::cout << "Error opening y-direction slope file { " << dhdyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (vxname != "") + { + if (vxDS == NULL) + { + std::cout << "Error opening x-direction velocity file { " << vxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (vyDS == NULL) + { + std::cout << "Error opening y-direction velocity file { " << vyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (srxname != "") + { + if (srxDS == NULL) + { + std::cout << "Error opening x-direction search range file { " << srxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (sryDS == NULL) + { + std::cout << "Error opening y-direction search range file { " << sryname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (csminxname != "") + { + if (csminxDS == NULL) + { + std::cout << "Error opening x-direction chip size min file { " << csminxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (csminyDS == NULL) + { + std::cout << "Error opening y-direction chip size min file { " << csminyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (csmaxxname != "") + { + if (csmaxxDS == NULL) + { + std::cout << "Error opening x-direction chip size max file { " << csmaxxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (csmaxyDS == NULL) + { + std::cout << "Error opening y-direction chip size max file { " << csmaxyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (ssmDS == NULL) + { + std::cout << "Error opening stable surface mask file { " << ssmname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + + demDS->GetGeoTransform(geoTrans); + int demXSize = demDS->GetRasterXSize(); + int demYSize = demDS->GetRasterYSize(); + + + //Get offsets and size to read from DEM + int lOff = std::max( std::floor((ymax - geoTrans[3])/geoTrans[5]), 0.); + int lCount = std::min( std::ceil((ymin - geoTrans[3])/geoTrans[5]), demYSize-1.) - lOff; + + int pOff = std::max( std::floor((xmin - geoTrans[0])/geoTrans[1]), 0.); + int pCount = std::min( std::ceil((xmax - geoTrans[0])/geoTrans[1]), demXSize-1.) - pOff; + + + std::cout << "Xlimits : " << geoTrans[0] + pOff * geoTrans[1] << " " + << geoTrans[0] + (pOff + pCount) * geoTrans[1] << "\n"; + + + std::cout << "Ylimits : " << geoTrans[3] + (lOff + lCount) * geoTrans[5] << " " + << geoTrans[3] + lOff * geoTrans[5] << "\n"; + + std::cout << "Dimensions of geogrid: " << pCount << " x " << lCount << "\n\n"; + + + //Create GDAL Transformers + OGRSpatialReference demSRS(nullptr); + if (demSRS.importFromEPSG(epsgcode) != 0) + { + std::cout << "Could not create OGR spatial reference for EPSG code: " << epsgcode << "\n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(102); + } + + OGRSpatialReference llhSRS(nullptr); + if (llhSRS.importFromEPSG(4326) != 0) + { + std::cout << "Could not create OGR spatil reference for EPSG code: 4326 \n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(103); + } + + OGRCoordinateTransformation *fwdTrans = OGRCreateCoordinateTransformation( &demSRS, &llhSRS); + OGRCoordinateTransformation *invTrans = OGRCreateCoordinateTransformation( &llhSRS, &demSRS); + + //WGS84 ellipsoid only + cEllipsoid wgs84; + wgs84.a = 6378137.0; + wgs84.e2 = 0.0066943799901; + + //Initial guess for solution + double tmid = sensingStart + 0.5 * nLines / prf; + double satxmid[3]; + double satvmid[3]; + + if (interpolateWGS84Orbit(orbit, tmid, satxmid, satvmid) != 0) + { + std::cout << "Error with orbit interpolation for setup. \n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(104); + } +// std::cout << "Center Satellite Velocity: " << satvmid[0] << " " << satvmid[1] << " " << satvmid[2] << "\n"; +// std::cout << satxmid[0] << " " << satxmid[1] << " " << satxmid[2] << "\n"; + + std::vector demLine(pCount); + std::vector sxLine(pCount); + std::vector syLine(pCount); + std::vector vxLine(pCount); + std::vector vyLine(pCount); + std::vector srxLine(pCount); + std::vector sryLine(pCount); + std::vector csminxLine(pCount); + std::vector csminyLine(pCount); + std::vector csmaxxLine(pCount); + std::vector csmaxyLine(pCount); + std::vector ssmLine(pCount); + + GInt32 raster1[pCount]; + GInt32 raster2[pCount]; + GInt32 raster11[pCount]; + GInt32 raster22[pCount]; + + GInt32 sr_raster11[pCount]; + GInt32 sr_raster22[pCount]; + GInt32 csmin_raster11[pCount]; + GInt32 csmin_raster22[pCount]; + GInt32 csmax_raster11[pCount]; + GInt32 csmax_raster22[pCount]; + GInt32 ssm_raster[pCount]; + + double raster1a[pCount]; + double raster1b[pCount]; +// double raster1c[pCount]; + + double raster2a[pCount]; + double raster2b[pCount]; +// double raster2c[pCount]; + + double nodata; +// double nodata_out; + if (vxname != "") + { + int* pbSuccess = NULL; + nodata = vxDS->GetRasterBand(1)->GetNoDataValue(pbSuccess); + } +// nodata_out = -2000000000; + + const char *pszFormat = "GTiff"; + char **papszOptions = NULL; + std::string str = ""; + double adfGeoTransform[6] = { xmin, (xmax-xmin)/pCount, 0, ymax, 0, (ymin-ymax)/lCount }; + OGRSpatialReference oSRS; + char *pszSRS_WKT = NULL; + demSRS.exportToWkt( &pszSRS_WKT ); + + + + GDALDriver *poDriver; + poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriver == NULL ) + exit(107); + GDALDataset *poDstDS; + + str = pixlinename; + const char * pszDstFilename = str.c_str(); + poDstDS = poDriver->Create( pszDstFilename, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + + poDstDS->SetGeoTransform( adfGeoTransform ); + poDstDS->SetProjection( pszSRS_WKT ); +// CPLFree( pszSRS_WKT ); + + + GDALRasterBand *poBand1; + GDALRasterBand *poBand2; + poBand1 = poDstDS->GetRasterBand(1); + poBand2 = poDstDS->GetRasterBand(2); + poBand1->SetNoDataValue(nodata_out); + poBand2->SetNoDataValue(nodata_out); + + + + + GDALDriver *poDriverOff; + poDriverOff = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverOff == NULL ) + exit(107); + GDALDataset *poDstDSOff; + + str = offsetname; + const char * pszDstFilenameOff = str.c_str(); + poDstDSOff = poDriverOff->Create( pszDstFilenameOff, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSOff->SetGeoTransform( adfGeoTransform ); + poDstDSOff->SetProjection( pszSRS_WKT ); +// CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1Off; + GDALRasterBand *poBand2Off; + poBand1Off = poDstDSOff->GetRasterBand(1); + poBand2Off = poDstDSOff->GetRasterBand(2); + poBand1Off->SetNoDataValue(nodata_out); + poBand2Off->SetNoDataValue(nodata_out); + + + + GDALDriver *poDriverSch; + poDriverSch = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverSch == NULL ) + exit(107); + GDALDataset *poDstDSSch; + + str = searchrangename; + const char * pszDstFilenameSch = str.c_str(); + poDstDSSch = poDriverSch->Create( pszDstFilenameSch, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSSch->SetGeoTransform( adfGeoTransform ); + poDstDSSch->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1Sch; + GDALRasterBand *poBand2Sch; + poBand1Sch = poDstDSSch->GetRasterBand(1); + poBand2Sch = poDstDSSch->GetRasterBand(2); + poBand1Sch->SetNoDataValue(nodata_out); + poBand2Sch->SetNoDataValue(nodata_out); + + + + GDALDriver *poDriverMin; + poDriverMin = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverMin == NULL ) + exit(107); + GDALDataset *poDstDSMin; + + str = chipsizeminname; + const char * pszDstFilenameMin = str.c_str(); + poDstDSMin = poDriverMin->Create( pszDstFilenameMin, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSMin->SetGeoTransform( adfGeoTransform ); + poDstDSMin->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1Min; + GDALRasterBand *poBand2Min; + poBand1Min = poDstDSMin->GetRasterBand(1); + poBand2Min = poDstDSMin->GetRasterBand(2); + poBand1Min->SetNoDataValue(nodata_out); + poBand2Min->SetNoDataValue(nodata_out); + + + + + GDALDriver *poDriverMax; + poDriverMax = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverMax == NULL ) + exit(107); + GDALDataset *poDstDSMax; + + str = chipsizemaxname; + const char * pszDstFilenameMax = str.c_str(); + poDstDSMax = poDriverMax->Create( pszDstFilenameMax, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSMax->SetGeoTransform( adfGeoTransform ); + poDstDSMax->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1Max; + GDALRasterBand *poBand2Max; + poBand1Max = poDstDSMax->GetRasterBand(1); + poBand2Max = poDstDSMax->GetRasterBand(2); + poBand1Max->SetNoDataValue(nodata_out); + poBand2Max->SetNoDataValue(nodata_out); + + + + + + GDALDriver *poDriverMsk; + poDriverMsk = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverMsk == NULL ) + exit(107); + GDALDataset *poDstDSMsk; + + str = stablesurfacemaskname; + const char * pszDstFilenameMsk = str.c_str(); + poDstDSMsk = poDriverMsk->Create( pszDstFilenameMsk, pCount, lCount, 1, GDT_Int32, + papszOptions ); + + poDstDSMsk->SetGeoTransform( adfGeoTransform ); + poDstDSMsk->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1Msk; + poBand1Msk = poDstDSMsk->GetRasterBand(1); + poBand1Msk->SetNoDataValue(nodata_out); + + + + + GDALDriver *poDriverRO2VX; + poDriverRO2VX = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverRO2VX == NULL ) + exit(107); + GDALDataset *poDstDSRO2VX; + + str = ro2vx_name; + const char * pszDstFilenameRO2VX = str.c_str(); + poDstDSRO2VX = poDriverRO2VX->Create( pszDstFilenameRO2VX, pCount, lCount, 2, GDT_Float64, + papszOptions ); + + poDstDSRO2VX->SetGeoTransform( adfGeoTransform ); + poDstDSRO2VX->SetProjection( pszSRS_WKT ); +// CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1RO2VX; + GDALRasterBand *poBand2RO2VX; +// GDALRasterBand *poBand3Los; + poBand1RO2VX = poDstDSRO2VX->GetRasterBand(1); + poBand2RO2VX = poDstDSRO2VX->GetRasterBand(2); +// poBand3Los = poDstDSLos->GetRasterBand(3); + poBand1RO2VX->SetNoDataValue(nodata_out); + poBand2RO2VX->SetNoDataValue(nodata_out); +// poBand3Los->SetNoDataValue(nodata_out); + + + + GDALDriver *poDriverRO2VY; + poDriverRO2VY = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverRO2VY == NULL ) + exit(107); + GDALDataset *poDstDSRO2VY; + + str = ro2vy_name; + const char * pszDstFilenameRO2VY = str.c_str(); + poDstDSRO2VY = poDriverRO2VY->Create( pszDstFilenameRO2VY, pCount, lCount, 2, GDT_Float64, + papszOptions ); + + poDstDSRO2VY->SetGeoTransform( adfGeoTransform ); + poDstDSRO2VY->SetProjection( pszSRS_WKT ); + CPLFree( pszSRS_WKT ); + + GDALRasterBand *poBand1RO2VY; + GDALRasterBand *poBand2RO2VY; +// GDALRasterBand *poBand3Alt; + poBand1RO2VY = poDstDSRO2VY->GetRasterBand(1); + poBand2RO2VY = poDstDSRO2VY->GetRasterBand(2); +// poBand3Alt = poDstDSAlt->GetRasterBand(3); + poBand1RO2VY->SetNoDataValue(nodata_out); + poBand2RO2VY->SetNoDataValue(nodata_out); +// poBand3Alt->SetNoDataValue(nodata_out); + + + // ground range and azimuth pixel size + double grd_res, azm_res; + +// double incang = 38.0*deg2rad; + double incang = incidenceAngle; + grd_res = dr / std::sin(incang); + azm_res = norm_C(satvmid) / prf; + std::cout << "Ground range pixel size: " << grd_res << "\n"; + std::cout << "Azimuth pixel size: " << azm_res << "\n"; +// int ChipSizeX0 = 240; + double ChipSizeX0 = chipSizeX0; + int ChipSizeX0_PIX_grd = std::ceil(ChipSizeX0 / grd_res / 4) * 4; + int ChipSizeX0_PIX_azm = std::ceil(ChipSizeX0 / azm_res / 4) * 4; + + + + for (int ii=0; iiGetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (demLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from DEM file: " << demname << "\n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(105); + } + + if (dhdxname != "") + { + status = sxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (sxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction slope file: " << dhdxname << "\n"; + GDALClose(sxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = syDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (syLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction slope file: " << dhdyname << "\n"; + GDALClose(syDS); + GDALDestroyDriverManager(); + exit(105); + } + + } + + if (vxname != "") + { + status = vxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (vxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction velocity file: " << vxname << "\n"; + GDALClose(vxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = vyDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (vyLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction velocity file: " << vyname << "\n"; + GDALClose(vyDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + if (srxname != "") + { + status = srxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (srxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction search range file: " << srxname << "\n"; + GDALClose(srxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = sryDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (sryLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction search range file: " << sryname << "\n"; + GDALClose(sryDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + + if (csminxname != "") + { + status = csminxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csminxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction chip size min file: " << csminxname << "\n"; + GDALClose(csminxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = csminyDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csminyLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction chip size min file: " << csminyname << "\n"; + GDALClose(csminyDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + + if (csmaxxname != "") + { + status = csmaxxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csmaxxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction chip size max file: " << csmaxxname << "\n"; + GDALClose(csmaxxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = csmaxyDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csmaxyLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction chip size max file: " << csmaxyname << "\n"; + GDALClose(csmaxyDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + + + if (ssmname != "") + { + status = ssmDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (ssmLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from stable surface mask file: " << ssmname << "\n"; + GDALClose(ssmDS); + GDALDestroyDriverManager(); + exit(105); + } + + } + + + + + int rgind; + int azind; + + for (int jj=0; jjTransform(1, llh, llh+1, llh+2); + + //Bringing it into ISCE + if (GDAL_VERSION_MAJOR == 2) + { + llhi[0] = deg2rad * llh[1]; + llhi[1] = deg2rad * llh[0]; + } + else + { + llhi[0] = deg2rad * llh[0]; + llhi[1] = deg2rad * llh[1]; + } + + llhi[2] = llh[2]; + + //Convert to ECEF + latlon_C(&wgs84, xyz, llhi, LLH_2_XYZ); + +// if ((ii == (lCount+1)/2)&(jj == pCount/2)){ +// std::cout << xyz[0] << " " << xyz[1] << " " << xyz[2] << "\n"; +// } + + //Start the geo2rdr algorithm + double satx[3]; + double satv[3]; + double tprev; + + double tline = tmid; + double rngpix; + double los[3]; + double alt[3]; + double normal[3]; + + double dopfact; + double height; + double vhat[3], that[3], chat[3], nhat[3], delta[3], targVec[3], targXYZ[3], diffvec[3], temp[3], satvc[3], altc[3]; + double vmag; + double major, minor; + double satDist; + double alpha, beta, gamma; + double radius, hgt, zsch; + double a, b, costheta, sintheta; + double rdiff; + + for(int kk=0; kk<3; kk++) + { + satx[kk] = satxmid[kk]; + } + + for(int kk=0; kk<3; kk++) + { + satv[kk] = satvmid[kk]; + } + + //Iterations + for (int kk=0; kk<51;kk++) + { + tprev = tline; + + for(int pp=0; pp<3; pp++) + { + drpos[pp] = xyz[pp] - satx[pp]; + } + + rngpix = norm_C(drpos); + double fn = dot_C(drpos, satv); + double fnprime = -dot_C(satv, satv); + + tline = tline - fn/fnprime; + + if (interpolateWGS84Orbit(orbit, tline, satx, satv) != 0) + { + std::cout << "Error with orbit interpolation. \n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(106); + } + + } +// if ((ii==600)&&(jj==600)) +// { +// std::cout << "\n" << lOff+ii << " " << pOff+jj << " " << demLine[jj] << "\n"; +// } + rgind = std::round((rngpix - startingRange) / dr) + 1.; + azind = std::round((tline - sensingStart) * prf) + 1.; + + + //*********************Slant-range vector + + + unitvec_C(drpos, los); + + for(int pp=0; pp<3; pp++) + { + llh[pp] = xyz[pp] + los[pp] * dr; + } + + latlon_C(&wgs84, llh, llhi, XYZ_2_LLH); + + //Bringing it from ISCE into LLH + if (GDAL_VERSION_MAJOR == 2) + { + llh[0] = llhi[1] / deg2rad; + llh[1] = llhi[0] / deg2rad; + } + else + { + llh[0] = llhi[0] / deg2rad; + llh[1] = llhi[1] / deg2rad; + } + + llh[2] = llhi[2]; + + //Convert from LLH inplace to DEM coordinates + invTrans->Transform(1, llh, llh+1, llh+2); + + for(int pp=0; pp<3; pp++) + { + drpos[pp] = llh[pp] - targllh0[pp]; + } + unitvec_C(drpos, los); + + //*********************Along-track vector + + tline = tline + 1/prf; + + if (interpolateWGS84Orbit(orbit, tline, satx, satv) != 0) + { + std::cout << "Error with orbit interpolation. \n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(106); + } + //run the topo algorithm for new tline + dopfact = 0.0; + height = demLine[jj]; + unitvec_C(satv, vhat); + vmag = norm_C(satv); + + //Convert position and velocity to local tangent plane + major = wgs84.a; + minor = major * std::sqrt(1 - wgs84.e2); + + //Setup ortho normal system right below satellite + satDist = norm_C(satx); + temp[0] = (satx[0] / major); + temp[1] = (satx[1] / major); + temp[2] = (satx[2] / minor); + alpha = 1 / norm_C(temp); + radius = alpha * satDist; + hgt = (1.0 - alpha) * satDist; + + //Setup TCN basis - Geocentric + unitvec_C(satx, nhat); + for(int pp=0; pp<3; pp++) + { + nhat[pp] = -nhat[pp]; + } + cross_C(nhat,satv,temp); + unitvec_C(temp, chat); + cross_C(chat,nhat,temp); + unitvec_C(temp, that); + + + //Solve the range doppler eqns iteratively + //Initial guess + zsch = height; + + for (int kk=0; kk<10;kk++) + { + a = satDist; + b = (radius + zsch); + + costheta = 0.5 * (a / rngpix + rngpix / a - (b / a) * (b / rngpix)); + sintheta = std::sqrt(1-costheta*costheta); + + gamma = rngpix * costheta; + alpha = dopfact - gamma * dot_C(nhat,vhat) / dot_C(vhat,that); + beta = -lookSide * std::sqrt(rngpix * rngpix * sintheta * sintheta - alpha * alpha); + for(int pp=0; pp<3; pp++) + { + delta[pp] = alpha * that[pp] + beta * chat[pp] + gamma * nhat[pp]; + } + + for(int pp=0; pp<3; pp++) + { + targVec[pp] = satx[pp] + delta[pp]; + } + + latlon_C(&wgs84, targVec, llhi, XYZ_2_LLH); + llhi[2] = height; + latlon_C(&wgs84, targXYZ, llhi, LLH_2_XYZ); + + zsch = norm_C(targXYZ) - radius; + + for(int pp=0; pp<3; pp++) + { + diffvec[pp] = satx[pp] - targXYZ[pp]; + } + rdiff = rngpix - norm_C(diffvec); + } + + //Bringing it from ISCE into LLH + + if (GDAL_VERSION_MAJOR == 2) + { + llh[0] = llhi[1] / deg2rad; + llh[1] = llhi[0] / deg2rad; + } + else + { + llh[0] = llhi[0] / deg2rad; + llh[1] = llhi[1] / deg2rad; + } + + llh[2] = llhi[2]; + + //Convert from LLH inplace to DEM coordinates + invTrans->Transform(1, llh, llh+1, llh+2); + + for(int pp=0; pp<3; pp++) + { + alt[pp] = llh[pp] - targllh0[pp]; + } + unitvec_C(alt, temp); + + + if (dhdxname != "") + { + //*********************Local normal vector + unitvec_C(slp, normal); + for(int pp=0; pp<3; pp++) + { + normal[pp] = -normal[pp]; + } + } + else + { + for(int pp=0; pp<3; pp++) + { + normal[pp] = 0.0; + } + } + + if (vxname != "") + { + vel[2] = -(vel[0]*normal[0]+vel[1]*normal[1])/normal[2]; + } + + if (srxname != "") + { + schrng1[2] = -(schrng1[0]*normal[0]+schrng1[1]*normal[1])/normal[2]; + schrng2[2] = -(schrng2[0]*normal[0]+schrng2[1]*normal[1])/normal[2]; + } + + + if ((rgind > nPixels)|(rgind < 1)|(azind > nLines)|(azind < 1)) + { + raster1[jj] = nodata_out; + raster2[jj] = nodata_out; + raster11[jj] = nodata_out; + raster22[jj] = nodata_out; + + sr_raster11[jj] = nodata_out; + sr_raster22[jj] = nodata_out; + csmin_raster11[jj] = nodata_out; + csmin_raster22[jj] = nodata_out; + csmax_raster11[jj] = nodata_out; + csmax_raster22[jj] = nodata_out; + ssm_raster[jj] = nodata_out; + + raster1a[jj] = nodata_out; + raster1b[jj] = nodata_out; +// raster1c[jj] = nodata_out; + raster2a[jj] = nodata_out; + raster2b[jj] = nodata_out; +// raster2c[jj] = nodata_out; + } + else + { + raster1[jj] = rgind; + raster2[jj] = azind; + + if (vxname == "") + { + raster11[jj] = 0.; + raster22[jj] = 0.; + } + else if (vel[0] == nodata) + { + raster11[jj] = 0.; + raster22[jj] = 0.; + } + else + { + raster11[jj] = std::round(dot_C(vel,los)*dt/dr/365.0/24.0/3600.0*1); + raster22[jj] = std::round(dot_C(vel,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1); + } + + if (srxname == "") + { + sr_raster11[jj] = nodata_out; + sr_raster22[jj] = nodata_out; + } + else if (vel[0] == nodata) + { + sr_raster11[jj] = 0; + sr_raster22[jj] = 0; + } + else + { + sr_raster11[jj] = std::abs(std::round(dot_C(schrng1,los)*dt/dr/365.0/24.0/3600.0*1)); + sr_raster22[jj] = std::abs(std::round(dot_C(schrng1,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)); + if (std::abs(std::round(dot_C(schrng2,los)*dt/dr/365.0/24.0/3600.0*1)) > sr_raster11[jj]) + { + sr_raster11[jj] = std::abs(std::round(dot_C(schrng2,los)*dt/dr/365.0/24.0/3600.0*1)); + } + if (std::abs(std::round(dot_C(schrng2,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)) > sr_raster22[jj]) + { + sr_raster22[jj] = std::abs(std::round(dot_C(schrng2,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)); + } + if (sr_raster11[jj] == 0) + { + sr_raster11[jj] = 1; + } + if (sr_raster22[jj] == 0) + { + sr_raster22[jj] = 1; + } + } + + if (csminxname != "") + { + csmin_raster11[jj] = csminxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; + csmin_raster22[jj] = csminyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; + } + else + { + csmin_raster11[jj] = nodata_out; + csmin_raster22[jj] = nodata_out; + } + if (csmaxxname != "") + { + csmax_raster11[jj] = csmaxxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; + csmax_raster22[jj] = csmaxyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; + } + else + { + csmax_raster11[jj] = nodata_out; + csmax_raster22[jj] = nodata_out; + } + if (ssmname != "") + { + ssm_raster[jj] = ssmLine[jj]; + } + else + { + ssm_raster[jj] = nodata_out; + } + raster1a[jj] = normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[1]-normal[1]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + raster1b[jj] = -normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[1]-normal[1]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + raster2a[jj] = -normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[0]-normal[0]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + raster2b[jj] = normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[0]-normal[0]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + +// raster1a[jj] = los[0]*dt/dr/365.0/24.0/3600.0; +// raster1b[jj] = los[1]*dt/dr/365.0/24.0/3600.0; +// raster1c[jj] = los[2]*dt/dr/365.0/24.0/3600.0; +// raster2a[jj] = temp[0]*dt/norm_C(alt)/365.0/24.0/3600.0; +// raster2b[jj] = temp[1]*dt/norm_C(alt)/365.0/24.0/3600.0; +// raster2c[jj] = temp[2]*dt/norm_C(alt)/365.0/24.0/3600.0; + } + + +// std::cout << ii << " " << jj << "\n"; +// std::cout << rgind << " " << azind << "\n"; +// std::cout << raster1[jj][ii] << " " << raster2[jj][ii] << "\n"; +// std::cout << raster1[ii][jj] << "\n"; + } + + poBand1->RasterIO( GF_Write, 0, ii, pCount, 1, + raster1, pCount, 1, GDT_Int32, 0, 0 ); + poBand2->RasterIO( GF_Write, 0, ii, pCount, 1, + raster2, pCount, 1, GDT_Int32, 0, 0 ); + poBand1Off->RasterIO( GF_Write, 0, ii, pCount, 1, + raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Off->RasterIO( GF_Write, 0, ii, pCount, 1, + raster22, pCount, 1, GDT_Int32, 0, 0 ); + + poBand1Sch->RasterIO( GF_Write, 0, ii, pCount, 1, + sr_raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Sch->RasterIO( GF_Write, 0, ii, pCount, 1, + sr_raster22, pCount, 1, GDT_Int32, 0, 0 ); + poBand1Min->RasterIO( GF_Write, 0, ii, pCount, 1, + csmin_raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Min->RasterIO( GF_Write, 0, ii, pCount, 1, + csmin_raster22, pCount, 1, GDT_Int32, 0, 0 ); + poBand1Max->RasterIO( GF_Write, 0, ii, pCount, 1, + csmax_raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Max->RasterIO( GF_Write, 0, ii, pCount, 1, + csmax_raster22, pCount, 1, GDT_Int32, 0, 0 ); + poBand1Msk->RasterIO( GF_Write, 0, ii, pCount, 1, + ssm_raster, pCount, 1, GDT_Int32, 0, 0 ); + + poBand1RO2VX->RasterIO( GF_Write, 0, ii, pCount, 1, + raster1a, pCount, 1, GDT_Float64, 0, 0 ); + poBand2RO2VX->RasterIO( GF_Write, 0, ii, pCount, 1, + raster1b, pCount, 1, GDT_Float64, 0, 0 ); +// poBand3Los->RasterIO( GF_Write, 0, ii, pCount, 1, +// raster1c, pCount, 1, GDT_Float64, 0, 0 ); + poBand1RO2VY->RasterIO( GF_Write, 0, ii, pCount, 1, + raster2a, pCount, 1, GDT_Float64, 0, 0 ); + poBand2RO2VY->RasterIO( GF_Write, 0, ii, pCount, 1, + raster2b, pCount, 1, GDT_Float64, 0, 0 ); +// poBand3Alt->RasterIO( GF_Write, 0, ii, pCount, 1, +// raster2c, pCount, 1, GDT_Float64, 0, 0 ); + + } + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDS ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSOff ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSSch ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMin ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMax ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMsk ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSRO2VX ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSRO2VY ); + + GDALClose(demDS); + + if (dhdxname != "") + { + GDALClose(sxDS); + GDALClose(syDS); + } + + if (vxname != "") + { + GDALClose(vxDS); + GDALClose(vyDS); + } + + if (srxname != "") + { + GDALClose(srxDS); + GDALClose(sryDS); + } + + if (csminxname != "") + { + GDALClose(csminxDS); + GDALClose(csminyDS); + } + + if (csmaxxname != "") + { + GDALClose(csmaxxDS); + GDALClose(csmaxyDS); + } + + if (ssmname != "") + { + GDALClose(ssmDS); + } + + GDALDestroyDriverManager(); + +} +void geoGrid::computeBbox(double *wesn) +{ + std::cout << "\nEstimated bounding box: \n" + << "West: " << wesn[0] << "\n" + << "East: " << wesn[1] << "\n" + << "South: " << wesn[2] << "\n" + << "North: " << wesn[3] << "\n"; +} diff --git a/docker/Dockerfile b/docker/Dockerfile index 021901f..f671d4f 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -38,6 +38,7 @@ RUN set -ex \ libgcc \ libstdcxx-ng \ cmake \ + opencv \ && yum install -y uuid-devel x11-devel motif-devel jq \ && ln -sf /opt/conda/bin/cython /opt/conda/bin/cython3 \ && mkdir -p /opt/isce2/src