diff --git a/app/otbSARAltAmbig.cxx b/app/otbSARAltAmbig.cxx
index 64a40ca656731e34a69926fecdacc2e24913c488..49e3a90a1da2aacaa6fa578643fc4fb5c6dc005e 100644
--- a/app/otbSARAltAmbig.cxx
+++ b/app/otbSARAltAmbig.cxx
@@ -25,6 +25,7 @@
#include "otbSarSensorModelAdapter.h"
#include <iostream>
+#include <iomanip>
#include <string>
#include <fstream>
@@ -67,15 +68,21 @@ private:
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, "target", "Target geosar image");
- SetParameterDescription("target", "Target geosar image");
+ 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");
@@ -87,7 +94,7 @@ private:
AddRAMParameter();
SetDocExampleParameterValue("inmaster", "s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
- SetDocExampleParameterValue("target", " s1a_s6_slc_0180416T195057_vv_cap_verde.tif");
+ SetDocExampleParameterValue("inslave", " s1a_s6_slc_0180416T195057_vv_cap_verde.tif");
SetDocExampleParameterValue("lat", "10.5");
SetDocExampleParameterValue("lon", "5.8");
SetDocExampleParameterValue("height", "45.9");
@@ -106,7 +113,7 @@ private:
ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("inmaster");
// Start the target pipeline (Read SAR image metadata)
- ComplexFloatImageType::Pointer SARPtr_target = GetParameterComplexFloatImage("target");
+ ComplexFloatImageType::Pointer SARPtr_slave = GetParameterComplexFloatImage("inslave");
// Open the output file
std::ofstream output(GetParameterString("outfile"), std::ios::out | std::ios::trunc);
@@ -119,6 +126,8 @@ private:
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);
@@ -128,14 +137,13 @@ private:
otbAppLogINFO(<<"Bistatic : " << bistatic);
otb::ImageKeywordlist sarKWL = SARPtr->GetImageKeywordlist();
- otb::ImageKeywordlist sarKWL_target = SARPtr_target->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 demGeoPointRadian;
SarSensorModelAdapter::Point3DType xyzCart;
SarSensorModelAdapter::Point3DType satellitePosition[2];
SarSensorModelAdapter::Point3DType satelliteVelocity[2];
@@ -145,26 +153,16 @@ private:
demGeoPoint[1] = lat;
demGeoPoint[2] = height;
- demGeoPointRadian[0] = lon*M_PI/180.;
- demGeoPointRadian[1] = lat*M_PI/180.;
- demGeoPointRadian[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 : ");
- output << std::endl << "Inmaster Image :" << std::endl;
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]);
- output << "SATELLITE POSITIONS = " << satellitePosition[0][0] << " " << satellitePosition[0][1] << " " << satellitePosition[0][2] << std::endl;
- output << "SATELLITE VELOCITY = " << satelliteVelocity[0][0] << " " << satelliteVelocity[0][1] << " " << satelliteVelocity[0][2] << std::endl;
-
- searchPassageCloserToGround(satellitePosition[0], incidence[0], xyzCart, demGeoPointRadian, output);
- otbAppLogINFO(<<"Incidence : " << incidence[0]);
// Calculate lambda
double radarFreq = atof(sarKWL.GetMetadataByKey("support_data.radar_frequency").c_str());
@@ -174,22 +172,28 @@ private:
otbAppLogINFO(<<"Lambda : " << lambda);
// Target
- loadOk = m_SarSensorModelAdapter->LoadState(sarKWL_target);
+ loadOk = m_SarSensorModelAdapter->LoadState(sarKWL_slave);
otbAppLogINFO(<<"Target Image : ");
- output << std::endl << "Target Image : " << std::endl;
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]);
- output << "SATELLITE POSITIONS = " << satellitePosition[1][0] << " " << satellitePosition[1][1] << " " << satellitePosition[2][2] << std::endl;
- output << "SATELLITE VELOCITY = " << satelliteVelocity[1][0] << " " << satelliteVelocity[1][1] << " " << satelliteVelocity[1][2] << std::endl;
-
- searchPassageCloserToGround(satellitePosition[1], incidence[1], xyzCart, demGeoPointRadian, output);
- otbAppLogINFO(<<"Incidence : " << incidence[1]);
otbAppLogINFO(<< "Compute Alt Ambig");
- otbAppLogINFO(<< "Master ORBIT\tSLAVE ORBIT\tIncidence\tALT AMBIG\tRadial \tLateral");
- output << std::endl << "Master ORBIT\tSLAVE ORBIT\tIncidence\tALT AMBIG\tRadial \tLateral" << std::endl;
- computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, incidence, output);
+ 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();
@@ -201,29 +205,28 @@ private:
return filePath.substr(sepPos + 1, dotPos - sepPos - 1);
}
- void searchPassageCloserToGround(SarSensorModelAdapter::Point3DType satellitePosition, double& incidence, SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType demGeoPoint, std::ofstream& output){
- SarSensorModelAdapter::Point3DType rm;
- rm[0] = xyzCart[0] - satellitePosition[0];
- rm[1] = xyzCart[1] - satellitePosition[1];
- rm[2] = xyzCart[2] - satellitePosition[2];
-
+ 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 normeR = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]);
+ double normeRm = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]);
- incidence = acos((normeS * normeS - normeR * normeR - normeCible * normeCible) /
- (2 * normeCible * normeR) ) * 180. / M_PI;
+ 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 normeRm = sqrt( rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]);
- double lon = demGeoPoint[0];
- double lat = demGeoPoint[1];
+ 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);
@@ -239,9 +242,9 @@ private:
// 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;
+ 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);
@@ -250,7 +253,7 @@ private:
output << "Vertical Projection : " << verticalProjection << std::endl;
}
- void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType* satellitePosition, SarSensorModelAdapter::Point3DType* satelliteVelocity, double lambda, double* incidence, std::ofstream& output){
+ 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;
@@ -265,6 +268,9 @@ private:
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] +
@@ -299,18 +305,32 @@ private:
double ye = normeRm * tan(theta);
double xe = normeRe - normeRm / cos(theta);
- if (incidence[1] > incidence[0]){
+ if (incidenceRe > incidenceRm){
ye *= -1.0;
}
- double alpha = (incidence[1] - incidence[0]) * M_PI / 180.;
- double inc_moy = ((incidence[1] + incidence[0]) / 2) * M_PI / 180.;
+ 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 target = getFilename(GetParameterString("target"));
- otbAppLogINFO(<< master << "\t" << target << "\t" << inc_moy * 180. / M_PI << " \t" << alt_amb << "\t" << xe << "\t" << ye / factorBperp);
- output << master << "\t" << target << "\t" << inc_moy * 180. / M_PI << " \t" << alt_amb << "\t" << xe << "\t" << ye / factorBperp << std::endl;
+ 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);
}
};