Merge branch 'alt_ambiguity' into 'master'

Alt Ambiguity port to diapotb

See merge request remote_modules/diapotb!31

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}
+)

app/otbSARAltAmbig.cxx

+/*
+ * 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)