diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 151f0f4f5952232c9698010e7e135452e4339a6d..82547d13a14bd3f08fc3d986dce6086765074055 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -108,3 +108,8 @@ 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}
+)
diff --git a/app/otbSARAltAmbig.cxx b/app/otbSARAltAmbig.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..49e3a90a1da2aacaa6fa578643fc4fb5c6dc005e
--- /dev/null
+++ b/app/otbSARAltAmbig.cxx
@@ -0,0 +1,342 @@
+/*
+ * 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 "otbImageFileReader.h"
+#include "otbSarSensorModelAdapter.h"
+
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <fstream>
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class SARAltAmbig : public Application
+{
+public:
+ typedef SARAltAmbig Self;
+ typedef itk::SmartPointer<Self> Pointer;
+
+ itkNewMacro(Self);
+ itkTypeMacro(SARAltAmbig, otb::Wrapper::Application);
+
+
+private:
+ void DoInit() override
+ {
+ SetName("SARAltAmbig");
+ SetDescription("Evaluation of altitudes of ambiguity of orbits");
+
+ SetDocLongDescription("Calculate the angle alpha between "
+ "ground-master and ground-slave vectors. "
+ "Calculates the local incidence.\nThe ambiguity "
+ "altitude is given by: "
+ "lambda * sin (incidence) / 2 * tg (alpha)");
+
+ //Optional descriptors
+ SetDocAuthors("OTB-Team");
+ SetDocSeeAlso(" ");
+ AddDocTag(Tags::SAR);
+ AddDocTag("DiapOTB");
+
+ //Parameter declarations
+ AddParameter(ParameterType_InputImage, "inmaster", "Input Master geosar image");
+ SetParameterDescription("inmaster", "SAR Image master");
+
+ AddParameter(ParameterType_Float, "lat", "Tartget lat");
+ SetParameterDescription("lat", "Target Latitude");
+ SetMinimumParameterFloatValue("lat", -90.);
+ SetMaximumParameterFloatValue("lat", 90.);
+
+ AddParameter(ParameterType_Float, "lon", "Target long");
+ SetParameterDescription("lon", "Target Longitude");
+ SetMinimumParameterFloatValue("lon", -360.);
+ SetMaximumParameterFloatValue("lon", 360.);
+
+ AddParameter(ParameterType_Float, "height", "Target height");
+ SetParameterDescription("height", "Target Height");
+ SetMinimumParameterFloatValue("height", -100.);
+ SetMaximumParameterFloatValue("height", 9000.);
+
+ AddParameter(ParameterType_InputImage, "inslave", "Input Slave geosar image");
+ SetParameterDescription("inslave", "SAR image slave");
+
+ AddParameter(ParameterType_OutputFilename, "outfile", "Output file to store the results");
+ SetParameterDescription("outfile","Output file to store the results");
+
+ AddParameter(ParameterType_Bool, "bistatic", "Activate bistatic Mode");
+ SetParameterDescription("bistatic", "If true, set the formula in bistatic mode");
+ MandatoryOff("bistatic");
+
+ AddRAMParameter();
+
+ SetDocExampleParameterValue("inmaster", "s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
+ SetDocExampleParameterValue("inslave", " s1a_s6_slc_0180416T195057_vv_cap_verde.tif");
+ SetDocExampleParameterValue("lat", "10.5");
+ SetDocExampleParameterValue("lon", "5.8");
+ SetDocExampleParameterValue("height", "45.9");
+ SetDocExampleParameterValue("outfile", "alt_ambig.log");
+
+ }
+
+ void DoUpdateParameters() override
+ {
+ // Nothing to do here : all parameters are independent
+ }
+
+ void DoExecute() override
+ {
+ // Start the first pipeline (Read SAR image metadata)
+ ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("inmaster");
+
+ // Start the target pipeline (Read SAR image metadata)
+ ComplexFloatImageType::Pointer SARPtr_slave = GetParameterComplexFloatImage("inslave");
+
+ // Open the output file
+ std::ofstream output(GetParameterString("outfile"), std::ios::out | std::ios::trunc);
+
+ // Get lat, lon and height
+ float lat = GetParameterFloat("lat");
+ otbAppLogINFO(<<"Lat : "<<lat);
+ output << "Latitude of the target = " << lat << std::endl;
+
+ float lon = GetParameterFloat("lon");
+ otbAppLogINFO(<<"Lon : "<<lon);
+ output << "Longitude of the target = " << lon << std::endl;
+ if (lon < 0)
+ lon = 360.0 + lon;
+
+ float height = GetParameterFloat("height");
+ otbAppLogINFO(<<"Height : "<<height);
+ output << "Altitude of the target = " << height << std::endl;
+
+ bool bistatic = IsParameterEnabled("bistatic");
+ otbAppLogINFO(<<"Bistatic : " << bistatic);
+
+ otb::ImageKeywordlist sarKWL = SARPtr->GetImageKeywordlist();
+ otb::ImageKeywordlist sarKWL_slave = SARPtr_slave->GetImageKeywordlist();
+
+ // Create and Initilaze the SarSensorModelAdapter
+ SarSensorModelAdapter::Pointer m_SarSensorModelAdapter = SarSensorModelAdapter::New();
+ bool loadOk = m_SarSensorModelAdapter->LoadState(sarKWL);
+
+ SarSensorModelAdapter::Point3DType demGeoPoint;
+ SarSensorModelAdapter::Point3DType xyzCart;
+ SarSensorModelAdapter::Point3DType satellitePosition[2];
+ SarSensorModelAdapter::Point3DType satelliteVelocity[2];
+ double incidence[2];
+
+ demGeoPoint[0] = lon;
+ demGeoPoint[1] = lat;
+ demGeoPoint[2] = height;
+
+ SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
+ otbAppLogINFO(<<"Cartesian coords : " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2]);
+ output << "GROUND COORDINATES = " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2] << std::endl;
+
+ // Master
+ otbAppLogINFO(<<"Inmaster Image : ");
+ m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[0], satelliteVelocity[0]);
+
+ otbAppLogINFO(<<"Satellite Position : " << satellitePosition[0][0] << " " << satellitePosition[0][1] << " " << satellitePosition[0][2]);
+ otbAppLogINFO(<<"Satellite Velocity : " << satelliteVelocity[0][0] << " " << satelliteVelocity[0][1] << " " << satelliteVelocity[0][2]);
+
+ // Calculate lambda
+ double radarFreq = atof(sarKWL.GetMetadataByKey("support_data.radar_frequency").c_str());
+ const double C = 299792458.;
+
+ float lambda = C/radarFreq;
+ otbAppLogINFO(<<"Lambda : " << lambda);
+
+ // Target
+ loadOk = m_SarSensorModelAdapter->LoadState(sarKWL_slave);
+ otbAppLogINFO(<<"Target Image : ");
+ m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[1], satelliteVelocity[1]);
+ otbAppLogINFO(<<"Satellite Position : " << satellitePosition[1][0] << " " << satellitePosition[1][1] << " " << satellitePosition[1][2]);
+ otbAppLogINFO(<<"Satellite Velocity : " << satelliteVelocity[1][0] << " " << satelliteVelocity[1][1] << " " << satelliteVelocity[1][2]);
+
+ otbAppLogINFO(<< "Compute Alt Ambig");
+ otbAppLogINFO(<< std::left
+ << std::setw(60) << "MASTER IMAGE"
+ << std::setw(60) << "SLAVE IMAGE"
+ << std::setw(12) << "Incidence"
+ << std::setw(12) << "ALT AMBIG"
+ << std::setw(12) << "Radial"
+ << std::setw(12) << "Lateral");
+ output << std::endl << std::left
+ << std::setw(60) << "MASTER IMAGE"
+ << std::setw(60) << "SLAVE IMAGE"
+ << std::setw(12) << "Incidence"
+ << std::setw(12) << "ALT AMBIG"
+ << std::setw(12) << "Radial"
+ << std::setw(12) << "Lateral" << std::endl;
+ computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, output);
+
+ // Close the output file
+ output.close();
+ }
+
+ std::string getFilename(std::string filePath){
+ std::size_t dotPos = filePath.rfind(".");
+ std::size_t sepPos = filePath.rfind("/");
+ return filePath.substr(sepPos + 1, dotPos - sepPos - 1);
+ }
+
+ double computeIncidence(SarSensorModelAdapter::Point3DType satellitePosition, SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType rm){
+ // Incidence
+ 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;
+ }
+
+ void searchPassageCloserToGround(SarSensorModelAdapter::Point3DType satellitePosition, SarSensorModelAdapter::Point3DType xyzCart, std::ofstream& output, SarSensorModelAdapter::Point3DType rm){
+ 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;
+ 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);
+
+ vectj[0] = -sin(lon);
+ vectj[1] = cos(lon);
+ vectj[2] = 0.;
+
+ vectk[0] = cos(lat)*cos(lon);
+ vectk[1] = cos(lat)*sin(lon);
+ vectk[2] = sin(lat);
+
+ // Rm projection
+ output << "Ground target --> Satellite Vector" << std::endl;
+ double northProjection = -(vecti[0]*(-rm[0])+vecti[1]*(-rm[1])+vecti[2]*(-rm[2]))/normeRm;
+ double eastProjection = -(vectj[0]*(-rm[0])+vectj[1]*(-rm[1])+vectj[2]*(-rm[2]))/normeRm;
+ double verticalProjection = -(vectk[0]*(-rm[0])+vectk[1]*(-rm[1])+vectk[2]*(-rm[2]))/normeRm;
+ otbAppLogINFO(<< "North Projection : " << northProjection);
+ otbAppLogINFO(<< "East Projection : " << eastProjection);
+ otbAppLogINFO(<< "Vertical Projection : " << verticalProjection);
+ output << "North Projection : " << northProjection << std::endl;
+ output << "East Projection : " << eastProjection << std::endl;
+ output << "Vertical Projection : " << verticalProjection << std::endl;
+ }
+
+ void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType* satellitePosition, SarSensorModelAdapter::Point3DType* satelliteVelocity, double lambda, std::ofstream& output){
+ // 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;
+ rm[0] = -(xyzCart[0] - satellitePosition[0][0]);
+ rm[1] = -(xyzCart[1] - satellitePosition[0][1]);
+ rm[2] = -(xyzCart[2] - satellitePosition[0][2]);
+
+ SarSensorModelAdapter::Point3DType re;
+ re[0] = -(xyzCart[0] - satellitePosition[1][0]);
+ re[1] = -(xyzCart[1] - satellitePosition[1][1]);
+ re[2] = -(xyzCart[2] - satellitePosition[1][2]);
+
+ double incidenceRm = computeIncidence(satellitePosition[0], xyzCart, rm);
+ double incidenceRe = computeIncidence(satellitePosition[1], xyzCart, re);
+
+ // project the slave orbit slant range vector
+ // onto the master orbit zero doppler plane
+ double vmNorm = sqrt(satelliteVelocity[0][0] * satelliteVelocity[0][0] +
+ satelliteVelocity[0][1] * satelliteVelocity[0][1] +
+ satelliteVelocity[0][2] * satelliteVelocity[0][2]);
+
+ // unitary vector with the same direction as the master velocity
+ // this vector is orthogonal to the master zero doppler plane
+ SarSensorModelAdapter::Point3DType vmUnit;
+ vmUnit[0] = satelliteVelocity[0][0] / vmNorm;
+ vmUnit[1] = satelliteVelocity[0][1] / vmNorm;
+ vmUnit[2] = satelliteVelocity[0][2] / vmNorm;
+
+ // compute the component along vmUnit of the slave slant-range vector
+ // scalar product Re . vmUnit
+ double reVmunit = re[0] * vmUnit[0] + re[1] * vmUnit[1] + re[2] * vmUnit[2];
+
+ // remove the along-track component
+ // this gives the projection of the slave
+ // slant-range vector onto the master zero-doppler plane
+ re[0] -= reVmunit * vmUnit[0];
+ re[1] -= reVmunit * vmUnit[1];
+ re[2] -= reVmunit * vmUnit[2];
+
+ double normeRm = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]);
+ double normeRe = sqrt(re[0] * re[0] + re[1] * re[1] + re[2] * re[2]);
+
+ // Scalar product RM.RE
+ double scalarRmRe = rm[0] * re[0] + rm[1] * re[1] + rm[2] * re[2];
+ double theta = acos(scalarRmRe / (normeRm * normeRe));
+
+ double ye = normeRm * tan(theta);
+ double xe = normeRe - normeRm / cos(theta);
+
+ if (incidenceRe > incidenceRm){
+ ye *= -1.0;
+ }
+
+ double alpha = (incidenceRe - incidenceRm) * M_PI / 180.;
+ double inc_moy = ((incidenceRe + incidenceRm) / 2) * M_PI / 180.;
+
+ double alt_amb = lambda * sin(inc_moy) / (factorHamb * tan(alpha));
+ std::string master = getFilename(GetParameterString("inmaster"));
+ std::string slave = getFilename(GetParameterString("inslave"));
+ otbAppLogINFO(<< std::left
+ << std::setw(60) << master
+ << std::setw(60) << slave
+ << std::setw(12) << inc_moy * 180. / M_PI
+ << std::setw(12) << alt_amb
+ << std::setw(12) << xe
+ << std::setw(12) << ye / factorBperp);
+ output << std::left
+ << std::setw(60) << master
+ << std::setw(60) << slave
+ << std::setw(12) << inc_moy * 180. / M_PI
+ << std::setw(12) << alt_amb
+ << std::setw(12) << xe
+ << std::setw(12) << ye / factorBperp << std::endl;
+
+ searchPassageCloserToGround(satellitePosition[0], xyzCart, output, rm);
+ }
+
+};
+
+}
+
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SARAltAmbig)