Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • remote_modules/diapotb
  • geninv/diapotb
  • lhermitte/diapotb
3 results
Show changes
Commits on Source (275)
Showing
with 878 additions and 900 deletions
/build/*
/build*/*
*/build/*
*/build*/*
/scripts/diapOTBEnv.sh
/data/*.omd
/data/*.ovr
......@@ -9,3 +11,8 @@
*.pyc
.\#*
\#*\#
/doc_cookbook/rst_tmp/*
# python packaging
*egg-info/
dist
\ No newline at end of file
DiapOTB is an official OTB remote module. It is available from OTB distributions.
All notable changes to DiapOTB will be documented in this file.
### OTB 8.x/v2.0.0
Coming, see [v2.0.0](https://gitlab.orfeo-toolbox.org/remote_modules/diapotb/-/tree/diapotb-2.0)
### OTB 7.4.1/v1.1.0
Coming with :
#### Added
* Unit tests with ctest
* TSX/PAZ/TDX support
* Refactoring of python_src
#### Fixed
* Correction on SARDEMGrid with Ossim projection
* Update utils script to retrieve fine orbit files
#### Removed
* old python chains
### OTB 7.3/v1.0.1
#### Added
* New utils scripts in python_src:
* Add GCPs on tiff images
* Generate configuration file from user's QA
* Retrieve fine orbit files from ESA website
* Improve memory consumption
* New post-processing : Goldstein filtering
#### Fixed
* Correction on projetion for `SARCartesianMeanEstimation` application
* Imporve fine orbit selection
* Update ESD to handle if only one burst is processed
### OTB 7.2/v1.0.0
#### Added
* Processing chains with S1 (mode IW and SM) and Cosmo (mode SM adn Spotligth) support.
* Burst sélection (S1IW) between reference and secondary image.
* Use fine orbit files (if available) in processing.
......@@ -6,13 +6,12 @@ if(NOT OTB_SOURCE_DIR)
list(APPEND CMAKE_MODULE_PATH ${OTB_CMAKE_DIR})
include(${OTB_USE_FILE})
include(OTBModuleExternal)
else()
else()
otb_module_impl()
endif()
# Execute and Copy during installation
install(DIRECTORY python_src DESTINATION ${CMAKE_INSTALL_PREFIX})
install(DIRECTORY json_schemas DESTINATION ${CMAKE_INSTALL_PREFIX})
#install(DIRECTORY python_src DESTINATION ${CMAKE_INSTALL_PREFIX})
if(UNIX)
install(CODE "execute_process(COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/scripts/generate_diapOTBEnv.sh ${CMAKE_INSTALL_PREFIX})")
install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/scripts/diapOTBEnv.sh PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ DESTINATION ${CMAKE_INSTALL_PREFIX})
......
......@@ -2,58 +2,140 @@
This is a module for the Diapason processing with OTB
The module must be used with OTB Orfeo Toolbox (https://www.orfeo-toolbox.org/).
It is designed to work with OTB v6.2 (or higher).
It is designed to work with OTB v7.0 (or higher).
This module can be built as a remote module or as a standalone module
# Getting Started
DiapOTB is an OTB's official remote module. Its applications and processing chains are directly included inside OTB's package.
Thus, two main ways can be used to install DiapOTB :
* From OTB standalone package (the easy way)
* From scratch (source code)
## Installation
1. Clone the repository
### Installation from OTB standalone package
This installation provides standalone binaries for several OS. You can choose the wanted binaries at https://www.orfeo-toolbox.org/download/
On all platforms, the OTB standalone package contains everything: applications for command line and graphical user interface, python bindings, Monteverdi and also official remote modules (with **DiapOTB**). This installation is extermely simplify by following given instructions : https://www.orfeo-toolbox.org/CookBook/Installation.html
### Installation from scratch
To install DiapOTB from scratch (from source code), binaries and libraries have to be overrided inside the existing installation of OTB. This override involves installing a “personal” OTB.
Therefore, two main steps are required for a complete installation of DiapOTB form scratch :
**1. OTB installation**
Find further information at https://www.orfeo-toolbox.org/CookBook/CompilingOTBFromSource.html
`
NB: A Superbuild is available for OTB installation to handle dependencies. With this build, dependencies are directly downloading inside a given repository (-DDOWNLOAD_LOCATION) before being installed.
`
* Clone the repository and choose a branch (usually develop by default)
```
%git clone
git clone https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb.git
git checkout <develop>
```
2. Go to source directory
* Create a build directory (choose a place)
```
%cd DiapOTBModule
% mkdir build
% cd build
```
* Setup environment for all dependencies and cmake
* Configure for a normal build
Cmake command makes configuration with an example given below (<install_OTB_directory> represents the selected repository for the OTB installation and <root_OTB> the root for the source code):
```
% CC=$GCCHOME/bin/gcc CXX=$GCCHOME/bin/g++ cmake <root_OTB> -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_COMPILER:FILEPATH=$GCCHOME/bin/g++ -DCMAKE_C_COMPILER:FILEPATH=$GCCHOME/bin/gcc -DCMAKE_CXX_FLAGS="-std=c++14" -DCMAK
E_INSTALL_PREFIX=<install_directory> -DOTB_USE_QWT=ON -DOTB_USE_OPENGL=ON -DOTB_USE_GLEW=ON -DOTB_USE_GLUT=ON -DOTB_USE_QT=ON -DOTB_USE_MUPARSER=ON -DOTB_USE_MUPAR
SERX=ON -DOTB_WRAP_PYTHON=ON -DOTB_USE_LIBSVM=ON -DBUILD_TESTING=OFF -DOTB_I18N_MERGE_TS=OFF -DOTB_USE_CURL=OFF -DOTB_USE_GLFW=OFF -DOTB_USE_GSL=OFF -DOTB_USE_LIBKML=OFF -DOTB_USE_OPENCV=OFF -DOTB_USE_SHARK=OFF
-DOTB_USE_SPTW=OFF
```
`
NB: Configuration is sligthly different with superbuild. Please, find futher information at https://www.orfeo-toolbox.org/CookBook/CompilingOTBFromSource.html?highlight=superbuild
`
3. Create build directory
* Compile and install
```
% make
% make install
```
**2. DiapOTB installation**
* Clone the repository
```
% git clone
For OTB Gitlab : git clone https://gitlab.orfeo-toolbox.org/remote_modules/diapotb.git
```
* Go to source directory
```
% cd DiapOTB
```
* Create build directory
```
% mkdir build
% cd build
```
4. Setup OTB environment
* Setup environment for all dependencies and cmake
* Setup OTB environment by setting the following variables (with x.x = OTB version)
```
% export OTB_DIR=<install_OTB_directory>/lib/cmake/OTB-x.x
% export LD_LIBRARY_PATH=<install_OTB_directory>/lib:$ LD_LIBRARY_PATH
```
5. Configure, compile and install
* Configure, compile and install (<install_DiapOTB_directory> represents the selected repository for the DiapOTB installation)
```
% CC=${GCCHOME}/bin/gcc CXX=${GCCHOME}/bin/g++ cmake -DOTB_BUILD_MODULE_AS_STANDALONE=ON -DCMAKE_INSTALL_PREFIX=<install_directory> -DCMAKE_BUILD_TYPE=<Release or Debug> ..
% CC=${GCCHOME}/bin/gcc CXX=${GCCHOME}/bin/g++ cmake -DOTB_BUILD_MODULE_AS_STANDALONE=ON -DCMAKE_INSTALL_PREFIX=<install_DiapOTB_directory> -DCMAKE_BUILD_TYPE=<Release or Debug> ..
% make
% make install
```
6. Configure with tests, compile and run tests
* Prepare execution by setting:
```
% CC=${GCCHOME}/bin/gcc CXX=${GCCHOME}/bin/g++ cmake -DBUILD_TESTING=ON -DOTB_BUILD_MODULE_AS_STANDALONE=ON -DCMAKE_INSTALL_PREFIX=<install_directory> -DCMAKE_BUILD_TYPE=<Release or Debug> ..
% make
% ctest
% export OTB_APPLICATION_PATH=<install_OTB_directory>/lib/otb/applications:<install_DiapOTB_directory>/lib
% export LD_LIBRARY_PATH=<install_OTB_directory>/lib/:<install_DiapOTB_directory>/lib/:$LD_LIBRARY_PATH
% export PYTHONPATH=<install_OTB_directory>/lib:<install_OTB_directory>/lib/otb/python:$PATH
% export PATH=<install_OTB_directory>/bin:$PATH
```
Other environment variables can be setted such as OTB_GEOID_FILE or ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
## Usage to launch one application (ex : SARMultiLook)
### C++ application
1. Load environment
1. Load environment with
* otbenv for the first kind of installation (standalone package)
* by setting the following variable for the installation from scratch
```
export OTB_APPLICATION_PATH=<install_directory>/lib:$OTB_APPLICATION_PATH
export LD_LIBRARY_PATH=<install_directory>/lib/:$LD_LIBRARY_PATH
% export OTB_APPLICATION_PATH=<install_OTB_directory>/lib/otb/applications:<install_DiapOTB_directory>/lib
% export LD_LIBRARY_PATH=<install_OTB_directory>/lib/:<install_DiapOTB_directory>/lib/:$LD_LIBRARY_PATH
% export PYTHONPATH=<install_OTB_directory>/lib:<install_OTB_directory>/lib/otb/python:$PATH
% export PATH=<install_OTB_directory>/bin:$PATH
```
A file (diapOTBEnv.sh) was creating into install_directory during the installation.
This file can be used to load the environment. Other environment variables can be
setted such as OTB_GEOID_FILE or ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
Other environment variables can be setted such as OTB_GEOID_FILE or ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
2. Get help
```
......@@ -62,26 +144,30 @@ otbApplicationLauncherCommandLine SARMultiLook
3. Execute
```
otbApplicationLauncherCommandLine SARMultiLookApp -in ./Data_In/sar/Reunion/S1A_S4_SLC__without_calibration.SAFE/measurement/s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002.tiff -out ./Data_Out/out_SAR_Diap/s1a-s4-ml-vv.tif -mlran 3 -mlazi 3 -mlgain 0.1
otbApplicationLauncherCommandLine SARMultiLookApp -incomplex ./Data_In/sar/Reunion/S1A_S4_SLC__without_calibration.SAFE/measurement/s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002.tiff -out ./Data_Out/out_SAR_Diap/s1a-s4-ml-vv.tif -mlran 3 -mlazi 3 -mlgain 0.1
```
### Python
1. Load environment
1. Load environment
* otbenv for the first kind of installation (standalone package)
* by setting the following variable for the installation from scratch
```
export OTB_APPLICATION_PATH=<install_directory>/lib:$OTB_APPLICATION_PATH
export LD_LIBRARY_PATH=<install_directory>/lib/:$LD_LIBRARY_PATH
% export OTB_APPLICATION_PATH=<install_OTB_directory>/lib/otb/applications:<install_DiapOTB_directory>/lib
% export LD_LIBRARY_PATH=<install_OTB_directory>/lib/:<install_DiapOTB_directory>/lib/:$LD_LIBRARY_PATH
% export PYTHONPATH=<install_OTB_directory>/lib:<install_OTB_directory>/lib/otb/python:$PATH
% export PATH=<install_OTB_directory>/bin:$PATH
```
A file (diapOTBEnv.sh) was creating into install_directory during the installation.
This file can be used to load the environment. Other environment variables can be
setted such as OTB_GEOID_FILE or ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
Other environment variables can be setted such as OTB_GEOID_FILE or ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
2. Write a Python file
```
import otbApplication as otb
app = otb.Registry.CreateApplication("SARMultiLook")
app.SetParameterString("in", "s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002.tiff")
app.SetParameterString("incomplex", "s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002.tiff")
app.SetParameterString("out", "output.tif")
app.SetParameterInt("mlran",3)
app.SetParameterInt("mlazi",3)
......@@ -96,19 +182,35 @@ python <script_python.py>
## Usage to launch DiapOTB processing
1. Load environment
1. Configure Python (if needed)
Some packages are required to run DiapOTB's processing chains with :
* GDAL
* h5py
* jsonschema
* numpy
2. Load environment
* otbenv for the first kind of installation (standalone package)
* by setting the following variable for the installation from scratch
```
source <install_directory>/diapOTBEnv.sh
% export OTB_APPLICATION_PATH=<install_OTB_directory>/lib/otb/applications:<install_DiapOTB_directory>/lib
% export LD_LIBRARY_PATH=<install_OTB_directory>/lib/:<install_DiapOTB_directory>/lib/:$LD_LIBRARY_PATH
% export PYTHONPATH=<install_OTB_directory>/lib:<install_OTB_directory>/lib/otb/python::<install_DiapOTB_directory>/python_src/:$PYTHONPATH
% export PATH=<install_OTB_directory>/bin:$PATH
```
Other environment variables can be setted such as OTB_GEOID_FILE or
ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
2. Configure your processing with a configuration file
Other environment variables can be setted such as OTB_GEOID_FILE or ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS.
3. Configure your processing with a configuration file
```
Example into <install_directory>/python_src/ex_config.ini
Example into <install_directory>/python_src/ex_config repository
or
Launch python <install_directory>/python_src/utils/generateConfigFile.py that corresponds to a user's Q&A.
```
3. Launch the processing
4. Launch the processing with for example
```
python $DIAPOTB_INSTALL/python_src/diapOTB.py <configuration_file>
python <install_directory>/python_src/diapOTB.py <configuration_file>
```
Further information are available in python_src README to install/configure python interperter with a setup.py.
......@@ -44,11 +44,6 @@ OTB_CREATE_APPLICATION(NAME SARFineDeformationGrid
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
OTB_CREATE_APPLICATION(NAME SARFineMetadata
SOURCES otbSARFineMetadata.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
OTB_CREATE_APPLICATION(NAME SARCoRegistration
SOURCES otbSARCoRegistration.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
......@@ -84,6 +79,7 @@ OTB_CREATE_APPLICATION(NAME SARDeramp
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
OTB_CREATE_APPLICATION(NAME SARGridOffset
SOURCES otbSARGridOffset.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
......@@ -99,16 +95,12 @@ OTB_CREATE_APPLICATION(NAME SARCorrelationRough
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
OTB_CREATE_APPLICATION(NAME SARMetadataCorrection
SOURCES otbSARMetadataCorrection.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
OTB_CREATE_APPLICATION(NAME SAROrthoInterferogram
SOURCES otbSAROrthoInterferogram.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
OTB_CREATE_APPLICATION(NAME SARAltAmbig
SOURCES otbSARAltAmbig.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
......
......@@ -81,12 +81,10 @@ namespace otb
//Parameter declarations
AddParameter(ParameterType_InputImage, "ininterfamp", "Input Interferogram with correct amplitude");
SetParameterDescription("ininterfamp", "Input Interferogram with correct amplitude.");
MandatoryOff("ininterfamp");
AddParameter(ParameterType_InputImage, "incomplexamp", "Input Compensated complex image");
SetParameterDescription("incomplexamp", "Input Compensated complex image.");
MandatoryOff("incomplexamp");
AddParameter(ParameterType_InputImage, "ininterf", "Input Interferogram");
SetParameterDescription("ininterf", "Input Interferogram.");
......@@ -121,7 +119,23 @@ namespace otb
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
// Handle the dependency between incomplexamp and ininterfamp (two kinds of pipelines)
if (GetParameterByKey("incomplexamp")->HasValue())
{
DisableParameter("ininterfamp");
MandatoryOff("ininterfamp");
otbAppLogINFO(<<"Complex input");
}
else
{
EnableParameter("ininterfamp");
MandatoryOn("ininterfamp");
MandatoryOff("incomplexamp");
otbAppLogINFO(<<"Interferogram input");
}
}
void DoExecute() override
......@@ -150,17 +164,17 @@ namespace otb
{
FloatVectorImageType::Pointer inInterfAmpPtr = GetParameterFloatVectorImage("ininterfamp");
// Check ML factors (must be 1x1) thanks to the keyWordlist
otb::ImageKeywordlist KWL = inInterfAmpPtr->GetImageKeywordlist();
// Check ML factors (must be 1x1) thanks to the Metadata
otb::ImageMetadata MD = inInterfAmpPtr->GetImageMetadata();
int mlran_interf = 1;
int mlazi_interf = 1;
if (KWL.HasKey("support_data.ml_ran"))
if (MD.Has("ml_ran"))
{
mlran_interf = std::atoi(KWL.GetMetadataByKey("support_data.ml_ran").c_str());
mlran_interf = std::atoi(MD["ml_ran"].c_str());
}
if (KWL.HasKey("support_data.ml_azi"))
if (MD.Has("ml_azi"))
{
mlazi_interf = std::atoi(KWL.GetMetadataByKey("support_data.ml_azi").c_str());
mlazi_interf = std::atoi(MD["ml_azi"].c_str());
}
if (mlran_interf != 1 || mlazi_interf != 1)
......
......@@ -22,7 +22,8 @@
#include "otbWrapperApplicationFactory.h"
#include "otbImageFileReader.h"
#include "otbSarSensorModelAdapter.h"
#include "otbSarSensorModel.h"
#include "otbGeocentricTransform.h"
#include <iostream>
#include <iomanip>
......@@ -39,7 +40,7 @@ class SARAltAmbig : public Application
{
public:
typedef SARAltAmbig Self;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
itkTypeMacro(SARAltAmbig, otb::Wrapper::Application);
......@@ -67,16 +68,16 @@ private:
AddParameter(ParameterType_InputImageList, "inlist", "Input images list (for estimations)");
SetParameterDescription("inlist", "The list of images to estimate required values (local vecotrs, ambiguity ...).");
AddParameter(ParameterType_Float, "lat", "Tartget lat (degree)");
AddParameter(ParameterType_Float, "lat", "Target lat (degree)");
SetParameterDescription("lat", "Target Latitude (degree)");
SetMinimumParameterFloatValue("lat", -90.);
SetMaximumParameterFloatValue("lat", 90.);
AddParameter(ParameterType_Float, "lon", "Target long (degree)");
SetParameterDescription("lon", "Target Longitude (degree)");
SetMinimumParameterFloatValue("lon", -360.);
SetMaximumParameterFloatValue("lon", 360.);
AddParameter(ParameterType_Float, "height", "Target height (m)");
SetParameterDescription("height", "Target Height (m)");
SetMinimumParameterFloatValue("height", -100.);
......@@ -96,7 +97,7 @@ private:
SetDocExampleParameterValue("lon", "5.8");
SetDocExampleParameterValue("height", "45.9");
SetDocExampleParameterValue("outfile", "alt_ambig.log");
}
void DoUpdateParameters() override
......@@ -106,50 +107,53 @@ private:
void DoExecute() override
{
//////////// Get lat, lon and height (same for all couples) ////////////
float lat = GetParameterFloat("lat");
otbAppLogINFO(<<"Lat : "<<lat);
float lon = GetParameterFloat("lon");
otbAppLogINFO(<<"Lon : "<<lon);
if (lon < 0)
lon = 360.0 + lon;
float height = GetParameterFloat("height");
otbAppLogINFO(<<"Height : "<<height);
bool bistatic = IsParameterEnabled("bistatic");
bool bistatic = GetParameterInt("bistatic");
otbAppLogINFO(<<"Bistatic : " << bistatic);
// Estimation on target
SarSensorModelAdapter::Point3DType demGeoPoint;
SarSensorModelAdapter::Point3DType xyzCart;
SarSensorModel::Point3DType demGeoPoint;
SarSensorModel::Point3DType xyzCart;
demGeoPoint[0] = lon;
demGeoPoint[1] = lat;
demGeoPoint[2] = height;
// Get List parameters
FloatVectorImageListType::Pointer inList = this->GetParameterImageList("inlist");
std::vector<std::string> inListNames = GetParameterStringList("inlist");
// Cartesian conversion
SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
FloatVectorImageType::Pointer SARPtr = inList->GetNthElement(0);
SARPtr->UpdateOutputInformation();
xyzCart = Projection::WorldToEcef(demGeoPoint);
otbAppLogDEBUG(<<"Cartesian coords : " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2]);
//////////// Create couples and global storage ///////////////
std::vector<double *> altAmbigVector;
std::map<int, double *> groundMap;
std::map<int, std::string> namesMap;
// Get List parameters
FloatVectorImageListType::Pointer inList = this->GetParameterImageList("inlist");
std::vector<std::string> inListNames = GetParameterStringList("inlist");
//////////// For each couple ///////////////
// Nested lopp to create all couples
for (unsigned int i = 0; i < (inList->Size() - 1); i++)
for (unsigned int i = 0; i < (inList->Size() - 1); i++)
{
for (unsigned int j = i+1; j < inList->Size(); j++)
for (unsigned int j = i+1; j < inList->Size(); j++)
{
// Start the target pipeline (Read SAR image metadata)
FloatVectorImageType::Pointer SARPtr = inList->GetNthElement(i);
......@@ -159,23 +163,22 @@ private:
FloatVectorImageType::Pointer SARPtr_slave = inList->GetNthElement(j);
SARPtr_slave->UpdateOutputInformation();
// Get image keywordlists
otb::ImageKeywordlist sarKWL = SARPtr->GetImageKeywordlist();
otb::ImageKeywordlist sarKWL_slave = SARPtr_slave->GetImageKeywordlist();
// Get image metadata
auto sarKWL = SARPtr->GetImageMetadata();
auto sarKWL_slave = SARPtr_slave->GetImageMetadata();
// Create and Initilaze the SarSensorModelAdapter
SarSensorModelAdapter::Pointer m_SarSensorModelAdapter = SarSensorModelAdapter::New();
bool loadOk = m_SarSensorModelAdapter->LoadState(sarKWL);
// Create and Initilaze the SarSensorModel
SarSensorModel sarSensorModelMaster(sarKWL);
SarSensorModelAdapter::Point3DType satellitePosition[2];
SarSensorModelAdapter::Point3DType satelliteVelocity[2];
SarSensorModel::Point3DType satellitePosition[2];
SarSensorModel::Point3DType satelliteVelocity[2];
// Results allocations
double * resultsAltAmbig = new double[6];
altAmbigVector.push_back(resultsAltAmbig);
// Image names
// Image names
std::string master = getFilename(inListNames[i]);
std::string slave = getFilename(inListNames[j]);
......@@ -183,21 +186,16 @@ private:
int master_orbit = 0;
int slave_orbit = 0;
if (sarKWL.HasKey("support_data.abs_orbit") && sarKWL_slave.HasKey("support_data.abs_orbit"))
if (sarKWL.Has(MDNum::OrbitNumber) && sarKWL_slave.Has(MDNum::OrbitNumber))
{
master_orbit = std::atoi(sarKWL.GetMetadataByKey("support_data.abs_orbit").c_str());
slave_orbit = std::atoi(sarKWL_slave.GetMetadataByKey("support_data.abs_orbit").c_str());
master_orbit = sarKWL[MDNum::OrbitNumber];
slave_orbit = sarKWL_slave[MDNum::OrbitNumber];
}
else if (sarKWL.HasKey("manifest_data.abs_orbit") && sarKWL_slave.HasKey("manifest_data.abs_orbit"))
else
{
master_orbit = std::atoi(sarKWL.GetMetadataByKey("manifest_data.abs_orbit").c_str());
slave_orbit = std::atoi(sarKWL_slave.GetMetadataByKey("manifest_data.abs_orbit").c_str());
// Send an exception
itkExceptionMacro(<<"Impossible to get orbit number, please check your input images");
}
else
{
// Send an exception
itkExceptionMacro(<<"Impossible to get orbit number, please check your input images");
}
namesMap[master_orbit] = master;
namesMap[slave_orbit] = slave;
......@@ -208,46 +206,48 @@ private:
// Satellite position and velocity
// Master
otbAppLogDEBUG(<<"Inmaster Image : ");
m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[0], satelliteVelocity[0]);
sarSensorModelMaster.WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[0],
satelliteVelocity[0]);
otbAppLogDEBUG(<<"Satellite Position : " << satellitePosition[0][0] << " " << satellitePosition[0][1] << " " << satellitePosition[0][2]);
otbAppLogDEBUG(<<"Satellite Velocity : " << satelliteVelocity[0][0] << " " << satelliteVelocity[0][1] << " " << satelliteVelocity[0][2]);
// Slave
loadOk = m_SarSensorModelAdapter->LoadState(sarKWL_slave);
SarSensorModel sarSensorModelSlave(sarKWL_slave);
otbAppLogDEBUG(<<"Target Image : ");
m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[1], satelliteVelocity[1]);
sarSensorModelSlave.WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[1], satelliteVelocity[1]);
otbAppLogDEBUG(<<"Satellite Position : " << satellitePosition[1][0] << " " << satellitePosition[1][1] << " " << satellitePosition[1][2]);
otbAppLogDEBUG(<<"Satellite Velocity : " << satelliteVelocity[1][0] << " " << satelliteVelocity[1][1] << " " << satelliteVelocity[1][2]);
// Calculate lambda
double radarFreq = atof(sarKWL.GetMetadataByKey("support_data.radar_frequency").c_str());
double radarFreq = sarKWL[MDNum::RadarFrequency];
const double C = 299792458.;
float lambda = C/radarFreq;
otbAppLogDEBUG(<<"Lambda : " << lambda);
// Estimations : mean incidence, ambig_heigth, radial coord, lateral coord
otbAppLogDEBUG(<< "Compute Alt Ambig");
computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, master,
computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, master,
slave, &resultsAltAmbig[2]);
// Local vector estimation (if not already done)
if (groundMap.find(master_orbit) == groundMap.end())
{
{
double * resultsLocalVector = new double[3];
computeLocalVector(satellitePosition[0], xyzCart, resultsLocalVector);
groundMap[master_orbit] = resultsLocalVector;
}
if (groundMap.find(slave_orbit) == groundMap.end())
{
{
double * resultsLocalVector = new double[3];
computeLocalVector(satellitePosition[1], xyzCart, resultsLocalVector);
groundMap[slave_orbit] = resultsLocalVector;
}
}
} // End of nested loop
......@@ -270,8 +270,8 @@ private:
output << std::left;
for(std::map<int,std::string>::iterator it = namesMap.begin();
it != namesMap.end(); ++it)
for(std::map<int,std::string>::iterator it = namesMap.begin();
it != namesMap.end(); ++it)
{
output << std::setw(70) << it->second
<< std::setw(12) << it->first << std::endl;
......@@ -287,7 +287,7 @@ private:
<< std::setw(15) << "Radial"
<< std::setw(15) << "Lateral" << std::endl;
for(unsigned int i = 0; i < altAmbigVector.size(); i++ )
for(unsigned int i = 0; i < altAmbigVector.size(); i++ )
{
output << std::left << std::setprecision(6)
<< std::setw(15) << static_cast<int>(altAmbigVector[i][0])
......@@ -300,7 +300,7 @@ private:
}
// Write satellite vectors array
output << std::setprecision(6) <<
output << std::setprecision(6) <<
std::endl << "Ground target --> Satellite Vectors : " << std::endl;
output << std::left
<< std::setw(15) << "Image Orbit"
......@@ -308,8 +308,8 @@ private:
<< std::setw(18) << "East Projection"
<< std::setw(18) << "Vertical Projection" << std::endl;
for(std::map<int,double *>::iterator it = groundMap.begin();
it != groundMap.end(); ++it)
for(std::map<int,double *>::iterator it = groundMap.begin();
it != groundMap.end(); ++it)
{
double * localVector = it->second;
......@@ -319,21 +319,21 @@ private:
<< std::setw(18) << localVector[1]
<< std::setw(18) << localVector[2] << std::endl;
}
// Close the output file
output.close();
// Free Memory
for(unsigned int i = 0; i < altAmbigVector.size(); i++ )
// Free Memory
for(unsigned int i = 0; i < altAmbigVector.size(); i++ )
{
delete altAmbigVector[i];
altAmbigVector[i] = 0;
}
altAmbigVector.clear();
for(std::map<int,double *>::iterator it = groundMap.begin();
it != groundMap.end(); ++it)
for(std::map<int,double *>::iterator it = groundMap.begin();
it != groundMap.end(); ++it)
{
delete it->second;
it->second = 0;
......@@ -343,16 +343,16 @@ private:
namesMap.clear();
}
/////// getFilename function ////////
std::string getFilename(std::string filePath)
{
// Remove HDF5 datasets into filePath
// Remove HDF5 datasets into filePath
std::size_t HDF5Pos = filePath.rfind("HDF5:");
std::size_t DatasetPos = filePath.rfind("://");
std::string fname = filePath;
if(HDF5Pos != std::string::npos)
......@@ -376,7 +376,7 @@ private:
std::size_t dotPos = fname.rfind(".");
std::size_t sepPos = fname.rfind("/");
return fname.substr(sepPos + 1, dotPos - sepPos - 1);
}
......@@ -384,43 +384,43 @@ private:
///////// compute functions /////////
// Incidence
double computeIncidence(SarSensorModelAdapter::Point3DType satellitePosition,
SarSensorModelAdapter::Point3DType xyzCart,
SarSensorModelAdapter::Point3DType rm)
double computeIncidence(SarSensorModel::Point3DType satellitePosition,
SarSensorModel::Point3DType xyzCart,
SarSensorModel::Point3DType rm)
{
// Incidence
double normeS = sqrt(satellitePosition[0]*satellitePosition[0] +
double normeS = sqrt(satellitePosition[0]*satellitePosition[0] +
satellitePosition[1]*satellitePosition[1] +
satellitePosition[2]*satellitePosition[2]);
double normeCible = sqrt(xyzCart[0] * xyzCart[0] + xyzCart[1] * xyzCart[1] + xyzCart[2] * xyzCart[2]);
double normeRm = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]);
double incidence = acos((normeS * normeS - normeRm * normeRm - normeCible * normeCible) /
(2 * normeCible * normeRm) ) * 180. / M_PI;
return incidence;
}
// Local Vector
void computeLocalVector(SarSensorModelAdapter::Point3DType satellitePosition,
SarSensorModelAdapter::Point3DType xyzCart,
void computeLocalVector(SarSensorModel::Point3DType satellitePosition,
SarSensorModel::Point3DType xyzCart,
double * resultsLocalVector)
{
SarSensorModelAdapter::Point3DType rm;
SarSensorModel::Point3DType rm;
rm[0] = -(xyzCart[0] - satellitePosition[0]);
rm[1] = -(xyzCart[1] - satellitePosition[1]);
rm[2] = -(xyzCart[2] - satellitePosition[2]);
double normeRm = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]);
// Calculate the local vectors
SarSensorModelAdapter::Point3DType vecti;
SarSensorModelAdapter::Point3DType vectj;
SarSensorModelAdapter::Point3DType vectk;
SarSensorModel::Point3DType vecti;
SarSensorModel::Point3DType vectj;
SarSensorModel::Point3DType vectk;
// lon an lat into radians
double lon = GetParameterFloat("lon")*M_PI/180.;
double lat = GetParameterFloat("lat")*M_PI/180.;
vecti[0] = -sin(lat)*cos(lon);
vecti[1] = -sin(lat)*sin(lon);
vecti[2] = cos(lat);
......@@ -449,22 +449,22 @@ private:
}
// Alt Ambig
void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart,
SarSensorModelAdapter::Point3DType* satellitePosition,
SarSensorModelAdapter::Point3DType* satelliteVelocity,
double lambda, std::string master, std::string slave,
void computeAltAmbig(SarSensorModel::Point3DType xyzCart,
SarSensorModel::Point3DType* satellitePosition,
SarSensorModel::Point3DType* satelliteVelocity,
double lambda, std::string master, std::string slave,
double * results)
{
// Adjust the formula in case of non-monostatic mode
double factorBperp = (!IsParameterEnabled("bistatic"))? 1.0:2.0;
double factorHamb = (!IsParameterEnabled("bistatic"))? 2.0:1.0;
SarSensorModelAdapter::Point3DType rm;
double factorBperp = (!GetParameterInt("bistatic"))? 1.0:2.0;
double factorHamb = (!GetParameterInt("bistatic"))? 2.0:1.0;
SarSensorModel::Point3DType rm;
rm[0] = -(xyzCart[0] - satellitePosition[0][0]);
rm[1] = -(xyzCart[1] - satellitePosition[0][1]);
rm[2] = -(xyzCart[2] - satellitePosition[0][2]);
SarSensorModelAdapter::Point3DType re;
SarSensorModel::Point3DType re;
re[0] = -(xyzCart[0] - satellitePosition[1][0]);
re[1] = -(xyzCart[1] - satellitePosition[1][1]);
re[2] = -(xyzCart[2] - satellitePosition[1][2]);
......@@ -480,7 +480,7 @@ private:
// unitary vector with the same direction as the master velocity
// this vector is orthogonal to the master zero doppler plane
SarSensorModelAdapter::Point3DType vmUnit;
SarSensorModel::Point3DType vmUnit;
vmUnit[0] = satelliteVelocity[0][0] / vmNorm;
vmUnit[1] = satelliteVelocity[0][1] / vmNorm;
vmUnit[2] = satelliteVelocity[0][2] / vmNorm;
......
......@@ -139,7 +139,7 @@ private:
// to estimate amplitude image)
FilterType::Pointer filterAmplitudeEstimation = FilterType::New();
m_Ref.push_back(filterAmplitudeEstimation.GetPointer());
filterAmplitudeEstimation->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
filterAmplitudeEstimation->SetSARImageMetadata(SARPtr->GetImageMetadata());
filterAmplitudeEstimation->SetSARImagePtr(SARPtr);
filterAmplitudeEstimation->SetDEMImagePtr(DEMPtr);
filterAmplitudeEstimation->SetDEMInformation(gain, DEMScanDirectionC, DEMScanDirectionL);
......
......@@ -140,7 +140,7 @@ private:
// to estimate cartesian coordonates mean image)
FilterType::Pointer filterCartesianMeanEstimation = FilterType::New();
m_Ref.push_back(filterCartesianMeanEstimation.GetPointer());
filterCartesianMeanEstimation->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
filterCartesianMeanEstimation->SetSARImageMetadata(SARPtr->GetImageMetadata());
filterCartesianMeanEstimation->SetSARImagePtr(SARPtr);
filterCartesianMeanEstimation->SetDEMImagePtr(DEMPtr);
filterCartesianMeanEstimation->SetDEMInformation(0, DEMScanDirectionC, DEMScanDirectionL);
......
......@@ -43,14 +43,14 @@ namespace otb
{
public:
typedef SARCoRegistration Self;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
itkTypeMacro(SARCoRegistration, otb::Wrapper::Application);
// Typedefs
typedef ComplexFloatImageType ComplexCoRegistrationImageType;
// Function
typedef otb::Function::SARTilesCoregistrationFunctor<ComplexCoRegistrationImageType, ComplexCoRegistrationImageType, FloatVectorImageType> TilesCoRegistrationFunctorType;
......@@ -93,12 +93,12 @@ namespace otb
SetParameterDescription("gridstepazimut", "Grid Step for azimut dimension into SLC geometry.");
SetDefaultParameterInt("gridstepazimut", 150);
MandatoryOff("gridstepazimut");
AddParameter(ParameterType_Float, "doppler0", "Doppler 0 in azimut for SAR Slave Image");
SetParameterDescription("doppler0", "Doppler 0 in azimut for SAR Slave Image.");
SetDefaultParameterFloat("doppler0", 0.);
MandatoryOff("doppler0");
AddParameter(ParameterType_Int, "sizetiles", "Size Tiles for output cut");
SetParameterDescription("sizetiles", "Size Tiles for output cut.");
SetDefaultParameterInt("sizetiles", 50);
......@@ -136,7 +136,7 @@ namespace otb
}
void DoExecute() override
{
{
// Get numeric parameters : sizeTiles, margin, gridStep, nbRamps and Doppler0
double doppler0 = GetParameterFloat("doppler0");
int gridStepRange = GetParameterInt("gridsteprange");
......@@ -144,18 +144,18 @@ namespace otb
int nbRamps = GetParameterInt("nbramps");
int margin = GetParameterInt("margin");
int sizeTiles = GetParameterInt("sizetiles");
otbAppLogINFO(<<"Doppler 0 in azimut : "<<doppler0);
otbAppLogINFO(<<"Grid Step Range : "<<gridStepRange);
otbAppLogINFO(<<"Grid Step Azimut : "<<gridStepAzimut);
otbAppLogINFO(<<"Number of Ramps for filtering : "<<nbRamps);
otbAppLogINFO(<<"Margin for window extraction into slave input image: "<<margin);
otbAppLogINFO(<<"Size of tiles for ouput (slave coregistrated image) cutting : "<<sizeTiles);
// Start the first pipeline (Read SAR Master image metedata)
ComplexCoRegistrationImageType::Pointer SARMasterPtr;
SARMasterPtr = GetParameterComplexFloatImage("insarmaster");
// Retrive input grid
FloatVectorImageType::Pointer GridPtr = GetParameterFloatVectorImage("ingrid");
......@@ -169,12 +169,12 @@ namespace otb
// Retrieve min/max into intput grid for each dimension
filterGridInfo->Update();
double minRan, minAzi, maxRan, maxAzi;
double minRan, minAzi, maxRan, maxAzi;
filterGridInfo->GetMinMax(minRan, minAzi, maxRan, maxAzi);
double maxAbsRan = std::max(std::abs(minRan), std::abs(maxRan));
double maxAbsAzi = std::max(std::abs(minAzi), std::abs(maxAzi));
double maxAbsAzi = std::max(std::abs(minAzi), std::abs(maxAzi));
///////////////////////////////////// CoRegistration Filter ///////////////////////////////////////////////
///////////////////////////////////// CoRegistration Filter ///////////////////////////////////////////////
// Function Type : TilesCoRegistation
unsigned int sizeWindow = sizeTiles + 2*margin;
TilesCoRegistrationFunctorType::Pointer functor = TilesCoRegistrationFunctorType::New();
......@@ -184,10 +184,8 @@ namespace otb
functor->SetNbRampes(nbRamps);
functor->SetDopAzi(doppler0);
functor->SetGridParameters(GridPtr, gridStepRange, gridStepAzimut);
if (IsParameterEnabled("cint16"))
{
functor->SetConvToCInt16(true);
}
functor->SetConvToCInt16(GetParameterInt("cint16"));
// Instanciate the Coregistration filter
TilesFilterType::Pointer filterCoregistration = TilesFilterType::New();
......@@ -200,7 +198,7 @@ namespace otb
// Master metadata
filterCoregistration->SetMasterImagePtr(SARMasterPtr);
// Define the main pipeline (controlled with extended FileName in order to obtain better performances)
//std::string origin_FileName = GetParameterString("out"); // Doesn't work with ComplexOutputImage
Parameter* param = GetParameterByKey("out");
......@@ -209,22 +207,22 @@ namespace otb
// Check if FileName is extended (with the ? caracter)
// If not extended then override the FileName
if (origin_FileName.find("?") == std::string::npos && !origin_FileName.empty())
if (origin_FileName.find("?") == std::string::npos && !origin_FileName.empty())
{
std::string extendedFileName = origin_FileName;
// Get the ram value (in MB)
int ram = GetParameterInt("ram");
// Define with the ram value, the number of lines for the streaming
int nbColSAR = SARMasterPtr->GetLargestPossibleRegion().GetSize()[0];
int nbLinesSAR = SARMasterPtr->GetLargestPossibleRegion().GetSize()[1];
// To determine the number of lines :
// To determine the number of lines :
// nbColSAR * nbLinesStreaming * sizeof(OutputPixel) = RamValue/2 (/2 to be sure)
// OutputPixel are float => sizeof(OutputPixel) = 4 (Bytes)
long int ram_In_KBytes = (ram/2) * 1024;
long int nbLinesStreaming = (ram_In_KBytes / (nbColSAR)) * (1024/4);
// Check the value of nbLinesStreaming
int nbLinesStreamingMax = 2000;
if (nbLinesStreamingMax > nbLinesSAR)
......@@ -234,16 +232,16 @@ namespace otb
if (nbLinesStreaming <= 0 || nbLinesStreaming > nbLinesStreamingMax)
{
nbLinesStreaming = nbLinesStreamingMax;
}
}
// Construct the extendedPart
std::ostringstream os;
os << "?&streaming:type=stripped&streaming:sizevalue=" << nbLinesStreaming;
// Add the extendedPart
std::string extendedPart = os.str();
extendedFileName.append(extendedPart);
// Set the new FileName with extended options
SetParameterString("out", extendedFileName);
}
......@@ -259,7 +257,7 @@ namespace otb
// Output
SetParameterOutputImage("out", filterCoregistration->GetOutput());
}
// Vector for filters
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
};
......@@ -269,5 +267,5 @@ namespace otb
}
OTB_APPLICATION_EXPORT(otb::Wrapper::SARCoRegistration)
......@@ -213,9 +213,6 @@ void DoExecute() override
double diffMeanRange = meanRangeCor - meanRangeDEM;
double diffMeanAzimut = meanAzimutCor - meanAzimutDEM;
std::cout << "diffMeanRange = " << diffMeanRange << std::endl;
std::cout << "diffMeanAzimut = " << diffMeanAzimut << std::endl;
// Correction Filter to replace value into DEMGrid with no occurences (high gap and high threshold)
// by mean value
CorrectionGridFilterType::Pointer filterCorrectionGrid = CorrectionGridFilterType::New();
......
......@@ -24,6 +24,8 @@
#include "itkMacro.h"
#include "itkFFTNormalizedCorrelationImageFilter.h"
#include "otbMultiToMonoChannelExtractROI.h"
#include "otbSARStreamingMaximumMinimumImageFilter.h"
#include "otbSARTemporalCorrelationGridImageFilter.h"
......@@ -49,7 +51,8 @@ public:
typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
typedef otb::SARStreamingMaximumMinimumImageFilter< FloatImageType> MinMaxFilterType;
typedef otb::SARTemporalCorrelationGridImageFilter<FloatImageType, FloatVectorImageType> CorGridFilterType;
typedef otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
FloatImageType::InternalPixelType> ExtractROIFilterType;
private:
void DoInit() override
......@@ -137,40 +140,67 @@ void DoExecute() override
otbAppLogINFO(<<"Grid Step for azimut : "<<grid_step_azi);
// Get master and slave image
FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
FloatVectorImageType::Pointer MasterPtr = GetParameterFloatVectorImage("inmaster");
FloatVectorImageType::Pointer SlavePtr = GetParameterFloatVectorImage("inslave");
// Check input size (Master size) to protect memory.
// If size > 15000 raise an exception and indicate more appropriate ml factors
if (MasterPtr->GetLargestPossibleRegion().GetSize()[0] > 15000 ||
MasterPtr->GetLargestPossibleRegion().GetSize()[1] > 15000)
{
// integer divisions
int intermediate_mlran = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/15000 + 1;
int intermediate_mlazi = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/15000 + 1;
// Initialize values
float shiftML_range = 0;
float shiftML_azimut = 0;
int startx = 0;
int starty = 0;
int sizex = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[0],
SlavePtr->GetLargestPossibleRegion().GetSize()[0]);
int sizey = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[1],
SlavePtr->GetLargestPossibleRegion().GetSize()[1]);
otbAppLogCRITICAL(<<"ML Factors are not appropriate for these estimations. Please use other factors such as : "<< intermediate_mlran << " x " << intermediate_mlazi );
// Extract ROI from Master and Slave Image to estimate global shifts without a too large amount
// of RAM => size = 3000*3000 at the center
if (sizex > 3000 &&
sizey > 3000)
{
startx = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/2 - 1500;
starty = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/2 - 1500;
itkExceptionMacro(<<"ML Factors are not appropriate for these estimations");
sizex = 3000;
sizey = 3000;
}
ExtractROIFilterType::Pointer extractROIFilter_master = ExtractROIFilterType::New();
m_Ref.push_back(extractROIFilter_master.GetPointer());
extractROIFilter_master->SetInput(MasterPtr);
extractROIFilter_master->SetChannel(1);
extractROIFilter_master->SetStartX(startx);
extractROIFilter_master->SetStartY(starty);
extractROIFilter_master->SetSizeX(sizex);
extractROIFilter_master->SetSizeY(sizey);
ExtractROIFilterType::Pointer extractROIFilter_slave = ExtractROIFilterType::New();
m_Ref.push_back(extractROIFilter_slave.GetPointer());
extractROIFilter_slave->SetInput(SlavePtr);
extractROIFilter_slave->SetChannel(1);
extractROIFilter_slave->SetStartX(startx);
extractROIFilter_slave->SetStartY(starty);
extractROIFilter_slave->SetSizeX(sizex);
extractROIFilter_slave->SetSizeY(sizey);
// Correlation Filter
CorFilterType::Pointer correlationFilter = CorFilterType::New();
m_Ref.push_back(correlationFilter.GetPointer());
correlationFilter->SetFixedImage(MasterPtr);
correlationFilter->SetMovingImage(SlavePtr);
correlationFilter->SetFixedImage(extractROIFilter_master->GetOutput());
correlationFilter->SetMovingImage(extractROIFilter_slave->GetOutput());
MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
minMaxFilter->SetInput(correlationFilter->GetOutput());
// Adapt streaming with ram parameter (default 256 MB)
minMaxFilter->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
minMaxFilter->Update();
float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
minMaxFilter->GetIndexOfMax()[0]);
float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-minMaxFilter->GetIndexOfMax()[1]);
shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
minMaxFilter->GetIndexOfMax()[0]);
shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-minMaxFilter->GetIndexOfMax()[1]);
// Correlation Grid Filter
CorGridFilterType::Pointer filterCorrelationGrid = CorGridFilterType::New();
m_Ref.push_back(filterCorrelationGrid.GetPointer());
......@@ -188,8 +218,18 @@ void DoExecute() override
// Define the main pipeline
filterCorrelationGrid->SetMasterInput(MasterPtr);
filterCorrelationGrid->SetSlaveInput(SlavePtr);
ExtractROIFilterType::Pointer extractROIFilter_masterfull = ExtractROIFilterType::New();
m_Ref.push_back(extractROIFilter_masterfull.GetPointer());
extractROIFilter_masterfull->SetInput(MasterPtr);
extractROIFilter_masterfull->SetChannel(1);
ExtractROIFilterType::Pointer extractROIFilter_slavefull = ExtractROIFilterType::New();
m_Ref.push_back(extractROIFilter_slavefull.GetPointer());
extractROIFilter_slavefull->SetInput(SlavePtr);
extractROIFilter_slavefull->SetChannel(1);
filterCorrelationGrid->SetMasterInput(extractROIFilter_masterfull->GetOutput());
filterCorrelationGrid->SetSlaveInput(extractROIFilter_slavefull->GetOutput());
SetParameterOutputImage("out", filterCorrelationGrid->GetOutput());
}
// Vector for filters
......
......@@ -24,6 +24,8 @@
#include "itkMacro.h"
#include "itkFFTNormalizedCorrelationImageFilter.h"
#include "otbMultiToMonoChannelExtractROI.h"
#include "otbSARStreamingMaximumMinimumImageFilter.h"
#include <iostream>
......@@ -47,6 +49,8 @@ public:
// Filters
typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
typedef otb::SARStreamingMaximumMinimumImageFilter< FloatImageType> MinMaxFilterType;
typedef otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
FloatImageType::InternalPixelType> ExtractROIFilterType;
private:
void DoInit() override
......@@ -129,46 +133,73 @@ void DoExecute() override
otbAppLogINFO(<<"RAM : "<< GetParameterInt("ram"));
// Get master and slave image
FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
FloatVectorImageType::Pointer MasterPtr = GetParameterFloatVectorImage("inmaster");
FloatVectorImageType::Pointer SlavePtr = GetParameterFloatVectorImage("inslave");
// Initialize values
float shiftML_range = 0;
float shiftML_azimut = 0;
float shiftSLC_range = 0;
float shiftSLC_azimut = 0;
int startx = 0;
int starty = 0;
int sizex = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[0],
SlavePtr->GetLargestPossibleRegion().GetSize()[0]);
int sizey = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[1],
SlavePtr->GetLargestPossibleRegion().GetSize()[1]);
// Check input size (Master size) to protect memory.
// If size > 15000 raise an exception and indicate more appropriate ml factors
if (MasterPtr->GetLargestPossibleRegion().GetSize()[0] > 15000 ||
MasterPtr->GetLargestPossibleRegion().GetSize()[1] > 15000)
// Extract ROI from Master and Slave Image to estimate global shifts without a too large amount
// of RAM => size = 3000*3000 at the center
if (sizex > 3000 &&
sizey > 3000)
{
// integer divisions
int intermediate_mlran = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/15000 + 1;
int intermediate_mlazi = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/15000 + 1;
startx = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/2 - 1500;
starty = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/2 - 1500;
otbAppLogCRITICAL(<<"ML Factors are not appropriate for these estimations. Please use other factors such as : "<< intermediate_mlran << " x " << intermediate_mlazi );
itkExceptionMacro(<<"ML Factors are not appropriate for these estimations");
sizex = 3000;
sizey = 3000;
}
ExtractROIFilterType::Pointer extractROIFilter_master = ExtractROIFilterType::New();
m_Ref.push_back(extractROIFilter_master.GetPointer());
extractROIFilter_master->SetInput(MasterPtr);
extractROIFilter_master->SetChannel(1);
extractROIFilter_master->SetStartX(startx);
extractROIFilter_master->SetStartY(starty);
extractROIFilter_master->SetSizeX(sizex);
extractROIFilter_master->SetSizeY(sizey);
ExtractROIFilterType::Pointer extractROIFilter_slave = ExtractROIFilterType::New();
m_Ref.push_back(extractROIFilter_slave.GetPointer());
extractROIFilter_slave->SetInput(SlavePtr);
extractROIFilter_slave->SetChannel(1);
extractROIFilter_slave->SetStartX(startx);
extractROIFilter_slave->SetStartY(starty);
extractROIFilter_slave->SetSizeX(sizex);
extractROIFilter_slave->SetSizeY(sizey);
// Correlation Filter
CorFilterType::Pointer correlationFilter = CorFilterType::New();
m_Ref.push_back(correlationFilter.GetPointer());
correlationFilter->SetFixedImage(MasterPtr);
correlationFilter->SetMovingImage(SlavePtr);
correlationFilter->SetFixedImage(extractROIFilter_master->GetOutput());
correlationFilter->SetMovingImage(extractROIFilter_slave->GetOutput());
MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
minMaxFilter->SetInput(correlationFilter->GetOutput());
// Adapt streaming with ram parameter (default 256 MB)
minMaxFilter->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
minMaxFilter->Update();
shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
minMaxFilter->GetIndexOfMax()[0]) * static_cast<float>(factorML_ran);
shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-minMaxFilter->GetIndexOfMax()[1]) * static_cast<float>(factorML_azi);
float shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
minMaxFilter->GetIndexOfMax()[0]) * static_cast<float>(factorML_ran);
float shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-minMaxFilter->GetIndexOfMax()[1]) * static_cast<float>(factorML_azi);
float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
minMaxFilter->GetIndexOfMax()[0]);
float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-minMaxFilter->GetIndexOfMax()[1]);
shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
minMaxFilter->GetIndexOfMax()[0]);
shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-minMaxFilter->GetIndexOfMax()[1]);
// Assigne Outputs
SetParameterFloat("shiftranslc", shiftSLC_range);
......
......@@ -122,7 +122,24 @@ private:
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
// Handle the dependency between indemproj* and indem (two kinds of pipelines)
if (GetParameterByKey("indemprojmaster")->HasValue() && GetParameterByKey("indemprojslave")->HasValue())
{
DisableParameter("indem");
MandatoryOff("indem");
otbAppLogINFO(<<"DEM Projection already estimated");
}
else
{
EnableParameter("indem");
MandatoryOn("indem");
MandatoryOff("indemprojmaster");
MandatoryOff("indemprojslave");
otbAppLogINFO(<<"DEM Projection must be estimated");
}
}
void DoExecute() override
......@@ -176,8 +193,8 @@ void DoExecute() override
m_Ref.push_back(filterDEMProjMaster.GetPointer());
m_Ref.push_back(filterDEMProjSlave.GetPointer());
filterDEMProjMaster->SetSARImageKeyWorList(MasterSARPtr->GetImageKeywordlist());
filterDEMProjSlave->SetSARImageKeyWorList(SlaveSARPtr->GetImageKeywordlist());
filterDEMProjMaster->SetSARImageMetadata(MasterSARPtr->GetImageMetadata());
filterDEMProjSlave->SetSARImageMetadata(SlaveSARPtr->GetImageMetadata());
// Start the Pipeline
filterDEMProjMaster->SetInput(inputDEM);
......
......@@ -39,7 +39,7 @@ class SARDEMProjection : public Application
{
public:
typedef SARDEMProjection Self;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
itkTypeMacro(SARDEMProjection, otb::Wrapper::Application);
......@@ -48,7 +48,7 @@ public:
typedef otb::SARDEMProjectionImageFilter<FloatImageType, FloatVectorImageType> DEMFilterType;
typedef otb::SARStreamingDEMInformationFilter<FloatImageType> DEMInfoFilterType;
typedef otb::SARStreamingDEMCheckFilter<FloatVectorImageType> DEMCheckFilterType;
private:
void DoInit() override
......@@ -75,17 +75,13 @@ private:
AddParameter(ParameterType_InputImage, "insar", "Input SAR image");
SetParameterDescription("insar", "SAR Image to extract SAR geometry.");
AddParameter(ParameterType_InputFilename, "infilemetadata", "Input file to get precise metadata");
SetParameterDescription("infilemetadata","Input file to get precise metadata.");
MandatoryOff("infilemetadata");
// Empty parameters to allow more components into the projeted DEM
AddParameter(ParameterType_Bool, "withh", "Set H components into projection");
SetParameterDescription("withh", "If true, then H components into projection.");
AddParameter(ParameterType_Bool, "withxyz", "Set XYZ Cartesian components into projection");
SetParameterDescription("withxyz", "If true, then XYZ Cartesian components into projection.");
AddParameter(ParameterType_Bool, "withsatpos", "Set SatPos components into projection");
SetParameterDescription("withsatpos", "If true, then SatPos components into projection.");
......@@ -125,74 +121,28 @@ private:
}
void DoExecute() override
{
{
// Get No data value
int noData = GetParameterInt("nodata");
otbAppLogINFO(<<"No Data : "<<noData);
// Start the first pipeline (Read SAR image metedata)
ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("insar");
// Add into SAR Image keyword list, the precise metadata
std::string inFile = GetParameterString("infilemetadata");
if (!inFile.empty())
{
std::ifstream preciseMetadataFile(inFile);
if (preciseMetadataFile.is_open())
{
std::string line;
std::string delimiter = " : ";
std::string token;
std::string value;
otb::ImageKeywordlist sarKWL = SARPtr->GetImageKeywordlist();
while ( std::getline (preciseMetadataFile,line) )
{
token = line.substr(0, line.find(delimiter));
value = line.substr(line.find(delimiter) + delimiter.length(), line.length());
if (token == "PRECISE TIME OF FIRST LINE")
{
sarKWL.AddKey("header.first_line_time", value);
sarKWL.AddKey("support_data.first_line_time ", value);
}
else if (token == "PRECISE SLANT RANGE")
{
sarKWL.AddKey("support_data.slant_range_to_first_pixel", value);
}
else
{
otbAppLogWARNING(<<"Wrong Metadata into input file ");
}
}
preciseMetadataFile.close();
// Uodate the SAR KeywordList
otbAppLogINFO(<<"SAR KeyWordList is updated with precise metadata");
SARPtr->SetImageKeywordList(sarKWL);
}
}
// DEMProjection Filter
DEMFilterType::Pointer filterDEMProjection = DEMFilterType::New();
m_Ref.push_back(filterDEMProjection.GetPointer());
filterDEMProjection->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
filterDEMProjection->SetSARImageMetadata(SARPtr->GetImageMetadata());
// Streaming Gain Filter
DEMInfoFilterType::Pointer filterStreamingDEMInfo = DEMInfoFilterType::New();
m_Ref.push_back(filterStreamingDEMInfo.GetPointer());
filterStreamingDEMInfo->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
filterStreamingDEMInfo->SetSARImageMetadata(SARPtr->GetImageMetadata());
// Streaming Check Filter (To check if DEM cover the SAR Image)
DEMCheckFilterType::Pointer filterStreamingDEMCheck = DEMCheckFilterType::New();
m_Ref.push_back(filterStreamingDEMCheck.GetPointer());
int sizeC = SARPtr->GetLargestPossibleRegion().GetSize()[0];
int sizeL = SARPtr->GetLargestPossibleRegion().GetSize()[1];
int originC = static_cast<int>(SARPtr->GetOrigin()[0] - 0.5);
......@@ -200,38 +150,33 @@ void DoExecute() override
filterStreamingDEMCheck->setSARSizeAndOrigin(sizeC, sizeL, originC, originL);
// Define the main pipeline
// Define the main pipeline
FloatImageType::Pointer inputDEM = GetParameterFloatImage("indem");
// Set Input and parameters
filterDEMProjection->SetInput(inputDEM);
if (IsParameterEnabled("withh"))
{
filterDEMProjection->SetwithH(true);
}
if (IsParameterEnabled("withxyz"))
{
filterDEMProjection->SetwithXYZ(true);
}
if (IsParameterEnabled("withsatpos"))
{
filterDEMProjection->SetwithSatPos(true);
}
filterDEMProjection->SetwithH(GetParameterInt("withh"));
filterDEMProjection->SetwithXYZ(GetParameterInt("withxyz"));
filterDEMProjection->SetwithSatPos(GetParameterInt("withsatpos"));
filterDEMProjection->SetNoData(noData);
// Define the pipeline to estimate the gain
filterStreamingDEMInfo->SetInput(inputDEM);
// Adapt streaming with ram parameter (default 256 MB)
filterStreamingDEMInfo->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
filterStreamingDEMInfo->Update();
double gain;
int directionDEM_forlines, directionDEM_forcolunms;
int directionDEM_forlines, directionDEM_forcolunms;
filterStreamingDEMInfo->GetDEMInformation(gain, directionDEM_forcolunms, directionDEM_forlines);
// Pass to DEMProjection filter, dem information to save into metadata
filterDEMProjection->SetdirectionToScanDEMC(directionDEM_forcolunms);
filterDEMProjection->SetdirectionToScanDEML(directionDEM_forlines);
filterDEMProjection->Setgain(gain);
// Check if DEM cover SAR Image (Counter must be > 0)
filterStreamingDEMCheck->SetInput(filterDEMProjection->GetOutput());
// Adapt streaming with ram parameter (default 256 MB)
......@@ -239,7 +184,7 @@ void DoExecute() override
filterStreamingDEMCheck->Update();
long int counterOfDEMPts_intoSARGeo = 0;
filterStreamingDEMCheck->GetDEMCheck(counterOfDEMPts_intoSARGeo);
if (counterOfDEMPts_intoSARGeo <= 0)
{
std::cout << "Input DEM does not cover the SAR geometry. Please check your DEM" << std::endl;
......@@ -254,7 +199,7 @@ void DoExecute() override
SetParameterInt("directiontoscandemc", directionDEM_forcolunms);
SetParameterInt("directiontoscandeml", directionDEM_forlines);
}
// Vector for filters
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
};
......
......@@ -21,7 +21,8 @@
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbSARDerampImageFilter.h"
#include "otbSARDerampS1IWImageFilter.h"
#include "otbSARDerampTSXImageFilter.h"
#include "otbWrapperOutputImageParameter.h"
#include "otbWrapperTypes.h"
......@@ -40,24 +41,29 @@ namespace otb
{
public:
typedef SARDeramp Self;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
itkTypeMacro(SARDeramp, otb::Wrapper::Application);
// Filter
typedef otb::SARDerampImageFilter<ComplexFloatImageType> DerampFilterType;
typedef otb::SARDerampS1IWImageFilter<ComplexFloatImageType> DerampS1IWFilterType;
typedef otb::SARDerampTSXImageFilter<ComplexFloatImageType> DerampTSXFilterType;
// MetaDataInterface
typedef otb::ImageMetadataInterfaceBase ImageMetadataInterfaceType;
private:
void DoInit() override
{
SetName("SARDeramp");
SetDescription("Deramping/Reramping of S1 IW Burst.");
SetDescription("Deramping/Reramping of S1 IW Burst and TSX mono-burst.");
SetDocLongDescription("This application does the deramping or reramping of S1 Iw burst.");
SetDocLongDescription("This application does the deramping or reramping of S1 Iw burst and TSX mono-burst.");
//Optional descriptors
SetDocLimitations("Only Sentinel 1 (IW and StripMap mode) and Cosmo products are supported for now.");
SetDocLimitations("Only Sentinel 1 (IW and StripMap mode), TSX and Cosmo products are supported for now.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso(" ");
AddDocTag(Tags::SAR);
......@@ -65,18 +71,16 @@ namespace otb
//Parameter declarations
AddParameter(ParameterType_InputImage, "in", "Input burst");
SetParameterDescription("in", "Input burst (from S1 Iw product).");
SetParameterDescription("in", "Input burst (from S1 Iw or TSX product).");
AddParameter(ParameterType_InputImage, "inslave", "Input SAR slave image (for interferometry processing");
SetParameterDescription("inslave", "Slave SAR Image to extract SAR geometry (interferometry).");
MandatoryOff("inslave");
AddParameter(ParameterType_OutputImage, "out", "Output burst");
SetParameterDescription("out","Output burst (deramped or reramped).");
AddParameter(ParameterType_InputImage, "ingrid", "Input deformation grid");
SetParameterDescription("ingrid", "Input deformation grid.");
MandatoryOff("ingrid");
AddParameter(ParameterType_Int, "gridsteprange", "Grid Step for range dimension into SLC geometry");
SetParameterDescription("gridsteprange", "Grid Step for range dimension into SLC geometry.");
......@@ -103,17 +107,33 @@ namespace otb
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
// Handle inputs following the mode chose by the user : shift mode
if (GetParameterInt("shift"))
{
// If shift mode : grid and slave image mandatory
EnableParameter("ingrid");
EnableParameter("inslave");
MandatoryOn("ingrid");
MandatoryOn("inslave");
}
else
{
// Classic deramp or reramp without deformation grid of slave metadata
DisableParameter("ingrid");
DisableParameter("inslave");
MandatoryOff("ingrid");
MandatoryOff("inslave");
}
}
void DoExecute() override
{
{
// Check if reramp is enabled
bool derampMode = true;
if (IsParameterEnabled("reramp"))
{
derampMode = false;
}
bool derampMode = !GetParameterInt("reramp");
if (derampMode)
{
......@@ -126,10 +146,10 @@ namespace otb
// Check if shift is enabled
bool shiftMode = false;
if (IsParameterEnabled("shift"))
if (GetParameterInt("shift"))
{
shiftMode = true;
// Check if required argument are there
if (!GetParameterByKey("ingrid")->HasValue() || !GetParameterByKey("gridsteprange")->HasValue() ||
!GetParameterByKey("gridstepazimut")->HasValue() || !GetParameterByKey("inslave")->HasValue())
......@@ -144,28 +164,58 @@ namespace otb
otbAppLogINFO(<<"Shift Mode");
}
// Get which sensor the image comes from
ImageMetadata metadata = GetParameterComplexFloatImage("in")->GetImageMetadata();
std::string sensorID = metadata[MDStr::SensorID];
// Instanciate the Interferogram filter
DerampFilterType::Pointer filterDeramp = DerampFilterType::New();
m_Ref.push_back(filterDeramp.GetPointer());
filterDeramp->SetDerampMode(derampMode);
filterDeramp->SetShiftMode(shiftMode);
if (shiftMode)
{
if(sensorID.find("TSX-1") != std::string::npos || sensorID.find("PAZ-1") != std::string::npos ||
sensorID.find("TDX-1") != std::string::npos){
// It is a TSX product so we apply the SARDerampTSXImageFilter
otbAppLogINFO(<<"TSX product");
DerampTSXFilterType::Pointer filterDeramp = DerampTSXFilterType::New();
m_Ref.push_back(filterDeramp.GetPointer());
filterDeramp->SetDerampMode(derampMode);
filterDeramp->SetShiftMode(shiftMode);
if (shiftMode)
{
filterDeramp->SetGridInput(GetParameterFloatVectorImage("ingrid"));
filterDeramp->SetSARImageKeyWorList(GetParameterComplexFloatImage("inslave")->GetImageKeywordlist());
filterDeramp->SetSARImageKeyWorList(GetParameterComplexFloatImage("inslave")->GetImageMetadata());
// Get numerical parameters
int gridStepRan = GetParameterInt("gridsteprange");
int gridStepAzi = GetParameterInt("gridstepazimut");
filterDeramp->SetGridStep(gridStepRan, gridStepAzi);
}
int gridStepRan = GetParameterInt("gridsteprange");
int gridStepAzi = GetParameterInt("gridstepazimut");
filterDeramp->SetGridStep(gridStepRan, gridStepAzi);
}
// Pipeline
filterDeramp->SetImageInput(GetParameterComplexFloatImage("in"));
SetParameterOutputImage("out", filterDeramp->GetOutput());
}
else{
// It is S1 IW product so we apply the SARDerampS1IWImageFilter
otbAppLogINFO(<<"S1 IW product");
DerampS1IWFilterType::Pointer filterDeramp = DerampS1IWFilterType::New();
m_Ref.push_back(filterDeramp.GetPointer());
filterDeramp->SetDerampMode(derampMode);
filterDeramp->SetShiftMode(shiftMode);
if (shiftMode)
{
filterDeramp->SetGridInput(GetParameterFloatVectorImage("ingrid"));
filterDeramp->SetSARImageKeyWorList(GetParameterComplexFloatImage("inslave")->GetImageMetadata());
// Pipeline
filterDeramp->SetImageInput(GetParameterComplexFloatImage("in"));
SetParameterOutputImage("out", filterDeramp->GetOutput());
// Get numerical parameters
int gridStepRan = GetParameterInt("gridsteprange");
int gridStepAzi = GetParameterInt("gridstepazimut");
filterDeramp->SetGridStep(gridStepRan, gridStepAzi);
}
// Pipeline
filterDeramp->SetImageInput(GetParameterComplexFloatImage("in"));
SetParameterOutputImage("out", filterDeramp->GetOutput());
}
}
// Vector for filters
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
};
......@@ -175,5 +225,5 @@ namespace otb
}
OTB_APPLICATION_EXPORT(otb::Wrapper::SARDeramp)
......@@ -182,15 +182,15 @@ void DoExecute() override
}
// Check invalid Pixel Key
const bool inputWithInvalidPixels1 = InterfUpPtr->GetImageKeywordlist().HasKey("support_data.invalid_pixels")
&& InterfUpPtr->GetImageKeywordlist().GetMetadataByKey("support_data.invalid_pixels") == "yes";
const bool inputWithInvalidPixels2 = InterfLowPtr->GetImageKeywordlist().HasKey("support_data.invalid_pixels")
&& InterfLowPtr->GetImageKeywordlist().GetMetadataByKey("support_data.invalid_pixels") == "yes";
const bool inputWithInvalidPixels1 = InterfUpPtr->GetImageMetadata().Has("invalid_pixels")
&& InterfUpPtr->GetImageMetadata()["invalid_pixels"] == "yes";
const bool inputWithInvalidPixels2 = InterfLowPtr->GetImageMetadata().Has("invalid_pixels")
&& InterfLowPtr->GetImageMetadata()["invalid_pixels"] == "yes";
if (inputWithInvalidPixels1!= inputWithInvalidPixels2)
{
// Throw an execption
otbAppLogFATAL(<< "Incoherency between input images (for support_data.invalid_pixels key).");
otbAppLogFATAL(<< "Incoherency between input images (for invalid_pixels key).");
}
......@@ -199,7 +199,7 @@ void DoExecute() override
std::pair<unsigned long,unsigned long> linesLow;
std::pair<unsigned long,unsigned long> samplesUp;
std::pair<unsigned long,unsigned long> samplesLow;
meanPhase->getOverlapLinesAndSamples(GetParameterComplexFloatImage("insar")->GetImageKeywordlist(),
meanPhase->getOverlapLinesAndSamples(GetParameterComplexFloatImage("insar")->GetImageMetadata(),
linesUp, linesLow, samplesUp, samplesLow, burstindex,
inputWithInvalidPixels1);
......@@ -243,23 +243,6 @@ void DoExecute() override
unsigned long startS_Low = samplesLow.first;
unsigned long sizeS_Low = minSamples_Low - samplesLow.first + 1;
std::cout << "startL_Up = " << startL_Up << std::endl;
std::cout << "sizeL_Up = " << sizeL_Up << std::endl;
std::cout << "startS_Up = " << startS_Up << std::endl;
std::cout << "sizeS_Up = " << sizeS_Up << std::endl;
std::cout << "startL_Low = " << startL_Low << std::endl;
std::cout << "sizeL_Low = " << sizeL_Low << std::endl;
std::cout << "startS_Low = " << startS_Low << std::endl;
std::cout << "sizeS_Low = " << sizeS_Low << std::endl;
std::cout << "samplesLow.first = " << samplesLow.first << std::endl;
std::cout << "samplesUp.first = " << samplesUp.first << std::endl;
std::cout << "samplesLow.second = " << samplesLow.second << std::endl;
std::cout << "samplesUp.second = " << samplesUp.second << std::endl;
std::cout << "minSamples_Low = " << minSamples_Low << std::endl;
std::cout << "minSamples_Up = " << minSamples_Up << std::endl;
////////// Extract Phase for overlap area //////////
// Set the channel to extract : Phase = 2
extractorUp->SetInput(InterfUpPtr);
......@@ -333,7 +316,6 @@ void DoExecute() override
meanPhase->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
meanPhase->Update();
std::cout << "meanPhase->GetMean() = " << meanPhase->GetMean() << std::endl;
// Doppler Centroid Filter
dcfUp->SetInput(extractorUp->GetOutput());
......@@ -358,7 +340,8 @@ void DoExecute() override
// //////////// ESD Filter ////////////
esdFilter->SetPhaseDiffMean(meanPhase->GetMean());
double aziInterval = std::stod(InterfUpPtr->GetImageKeywordlist().GetMetadataByKey("support_data.line_time_interval"));
//double aziInterval = std::stod(InterfUpPtr->GetImageMetadata().GetMetadataByKey("line_time_interval"));
double aziInterval = boost::any_cast<const otb::SARParam&>(InterfUpPtr->GetImageMetadata()[otb::MDGeom::SAR]).azimuthTimeInterval.TotalSeconds();
esdFilter->SetAziTimeInt(aziInterval);
//esdFilter->SetInput(subPhase->GetOutput());
esdFilter->SetInput1(dcfUp->GetOutput());
......@@ -378,8 +361,6 @@ void DoExecute() override
meanAziShift->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
meanAziShift->Update();
std::cout << "meanAziShift->GetMean() = " << meanAziShift->GetMean() << std::endl;
////////// Outputs : Image of azimut shift (for overlap area) and mean of azimut shift ////////////
//SetParameterOutputImage("out", dcfLow->GetOutput());
SetParameterOutputImage("out", subPhase->GetOutput());
......
/*
* Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/
*
* 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.
*/
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbSARStreamingGridHistogramFilter.h"
#include "otbSARUpdateMetadataImageFilter.h"
#include <iostream>
#include <string>
#include <fstream>
#include <complex>
#include <cmath>
// include ossim
#if defined(__GNUC__) || defined(__clang__)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include "ossim/ossimTimeUtilities.h"
#include "ossim/base/ossimKeywordlist.h"
#include "ossim/base/ossimString.h"
#include "ossim/base/ossimDate.h"
#pragma GCC diagnostic pop
#else
#include "ossim/ossimTimeUtilities.h"
#include "ossim/base/ossimKeywordlist.h"
#include "ossim/base/ossimString.h"
#include "ossim/base/ossimDate.h"
#endif
namespace otb
{
namespace Wrapper
{
class SARFineMetadata : public Application
{
public:
typedef SARFineMetadata Self;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
itkTypeMacro(SARFineMetadata, otb::Wrapper::Application);
// Filter
typedef otb::SARStreamingGridHistogramFilter<FloatVectorImageType> GridHistogramType;
typedef otb::SARUpdateMetadataImageFilter<ComplexFloatImageType> UpdateMetadataeFilter;
// Time/Date ossim
typedef ossimplugins::time::ModifiedJulianDate TimeType;
typedef ossimplugins::time::Duration DurationType;
private:
void DoInit() override
{
SetName("SARFineMetadata");
SetDescription("SAR metadata correction.");
SetDocLongDescription("This application corrects two metadata : the time "
"of the first line and the slant near range.");
//Optional descriptors
SetDocLimitations("Only Sentinel 1 (IW and StripMap mode) and Cosmo products are supported for now.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso(" ");
AddDocTag(Tags::SAR);
AddDocTag("DiapOTB");
//Parameter declarations
AddParameter(ParameterType_InputImage, "ingrid", "Input correlation grid");
SetParameterDescription("ingrid", "Correlation grid with shifts for the two dimensions and associated correlation rate between two images (for instance : ML Master and simulated amplitude image).");
AddParameter(ParameterType_InputImage, "insar", "Input SAR image");
SetParameterDescription("insar", "SAR Image to extract initial metadata.");
AddParameter(ParameterType_OutputFilename, "outfile", "Output file to store new metadata values");
SetParameterDescription("outfile","Output file to store new metadata values.");
MandatoryOff("outfile");
AddParameter(ParameterType_Float, "stepmax", "Max step for histogram");
SetParameterDescription("stepmax", "Max step for histogram.");
SetDefaultParameterFloat("stepmax", 0.1);
MandatoryOff("stepmax");
AddParameter(ParameterType_Float, "threshold", "Threshold on correlation rate");
SetParameterDescription("threshold", "Threshold on correlation rate.");
SetDefaultParameterFloat("threshold", 0.3);
MandatoryOff("threshold");
// Output parameters
AddParameter(ParameterType_String,"timefirstline","Time of the first line (output of this application)");
SetParameterDescription("timefirstline", "Time of the first line (output of this application).");
SetParameterRole("timefirstline", Role_Output);
AddParameter(ParameterType_Float,"slantrange","Slant near range (output of this application)");
SetParameterDescription("slantrange", "Slant near range (output of this application).");
SetParameterRole("slantrange", Role_Output);
AddParameter(ParameterType_OutputImage, "out", "Output image with precise metadata");
SetParameterDescription("out","Output image with precise metadata.");
AddRAMParameter();
SetDocExampleParameterValue("ingrid","corGrid.tiff");
SetDocExampleParameterValue("outfile","fine_metadata.txt");
}
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
}
void DoExecute() override
{
// Get numeric parameters : stepMax and threshold
float stepMax = GetParameterFloat("stepmax");
//float Threshold = GetParameterFloat("threshold");
// Get inputs
ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("insar");
FloatVectorImageType::Pointer CorGridPtr = GetParameterFloatVectorImage("ingrid");
// Get Sar keyWordList
otb::ImageKeywordlist sar_kwl = SARPtr->GetImageKeywordlist();
// Define C
const double C = 299792458;
//////////////////// Slant Range (range dimension) /////////////////////////
// Initial Parameters
double stepHistRange = -1.;
double meanShiftRange = 0.;
do {
// Configure Filter
GridHistogramType::Pointer statisticsRange = GridHistogramType::New();
m_Ref.push_back(statisticsRange.GetPointer());
statisticsRange->SetInput(CorGridPtr);
statisticsRange->SetDimension(0);
statisticsRange->SetStepHist(stepHistRange);
statisticsRange->SetThreshold(0.3);
statisticsRange->SetMeanShift(meanShiftRange);
statisticsRange->SetNumberHist(256);
// Execute pipeline
statisticsRange->Update();
//double step, meanShift;
statisticsRange->GetGridHistogramInformation(stepHistRange, meanShiftRange);
} while(stepHistRange > stepMax);
// New slant range
double slantRange = atof(sar_kwl.GetMetadataByKey("support_data.slant_range_to_first_pixel").c_str());
double pixelSpacing_Range = atof(sar_kwl.GetMetadataByKey("support_data.range_spacing").c_str());
double fineSlantRange = slantRange + ((2*meanShiftRange*pixelSpacing_Range)/C);
//////////////////// Time of the first line (azimut dimension) /////////////////////////
// Initial Parameters
double stepHistAzimut = -1.;
double meanShiftAzimut = 0.;
do {
// Configure Filter
GridHistogramType::Pointer statisticsAzimut = GridHistogramType::New();
m_Ref.push_back(statisticsAzimut.GetPointer());
statisticsAzimut->SetInput(CorGridPtr);
statisticsAzimut->SetDimension(1);
statisticsAzimut->SetStepHist(stepHistAzimut);
statisticsAzimut->SetThreshold(0.3);
statisticsAzimut->SetMeanShift(meanShiftAzimut);
statisticsAzimut->SetNumberHist(256);
// Execute pipeline
statisticsAzimut->Update();
//double step, meanShift;
statisticsAzimut->GetGridHistogramInformation(stepHistAzimut, meanShiftAzimut);
} while(stepHistAzimut > stepMax);
// New time for the first line
std::string timeFirstLine = sar_kwl.GetMetadataByKey("support_data.first_line_time");
double prf = atof(sar_kwl.GetMetadataByKey("support_data.pulse_repetition_frequency").c_str());
// Change into julian day
ossimplugins::time::ModifiedJulianDate time_first_pixel_Date = ossimplugins::time::toModifiedJulianDate(timeFirstLine);
ossimplugins::time::ModifiedJulianDate fine_Time_first_pixel_Date = time_first_pixel_Date + ossimplugins::time::seconds(meanShiftAzimut/prf);
// Change into std::string
std::string fineTimeFirstLine = to_simple_string(fine_Time_first_pixel_Date);
//////////////////// Output File to store fine metadata /////////////////////////
// if outfile parameter, fill the file with new metadata values
std::string outFile = GetParameterString("outfile");
if (!outFile.empty())
{
std::ofstream MetadataFile(outFile, std::ios::app);
if (MetadataFile)
{
MetadataFile << std::setprecision(11);
MetadataFile << "PRECISE TIME OF FIRST LINE : " << fineTimeFirstLine << std::endl;
MetadataFile << "PRECISE SLANT RANGE : " << fineSlantRange << std::endl;
MetadataFile.close();
}
}
// Assigne timefirstline and slantrange
SetParameterString("timefirstline", fineTimeFirstLine);
SetParameterFloat("slantrange", fineSlantRange);
// Copie input image and update the metadata into a naw geom file
UpdateMetadataeFilter::Pointer updateMetadataFilter = UpdateMetadataeFilter::New();
m_Ref.push_back(updateMetadataFilter.GetPointer());
updateMetadataFilter->SetTimeFirstLine(fineTimeFirstLine);
updateMetadataFilter->SetSlantRange(fineSlantRange);
updateMetadataFilter->SetInput(SARPtr);
SetParameterOutputImage("out", updateMetadataFilter->GetOutput());
}
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::SARFineMetadata)
This diff is collapsed.
......@@ -73,11 +73,9 @@ private:
//Parameter declarations
AddParameter(ParameterType_InputImage, "incomplex", "Input Complex image");
SetParameterDescription("incomplex", "Complex Image to perform computation on.");
MandatoryOff("incomplex");
AddParameter(ParameterType_InputImage, "inreal", "Input image");
SetParameterDescription("inreal", "Image to perform computation on.");
MandatoryOff("inreal");
AddParameter(ParameterType_Int, "mlran", "Averaging on distance");
SetParameterDescription("mlran","Averaging on distance");
......@@ -111,7 +109,23 @@ private:
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
// Handle the dependency between incomplex and inreal (two kinds of pipelines)
if (GetParameterByKey("incomplex")->HasValue())
{
DisableParameter("inreal");
MandatoryOff("inreal");
otbAppLogINFO(<<"Complex input");
}
else
{
EnableParameter("inreal");
MandatoryOn("inreal");
MandatoryOff("incomplex");
otbAppLogINFO(<<"Real input");
}
}
void DoExecute() override
......