Skip to content
Snippets Groups Projects
Commit 9e23ee8c authored by Gaëlle USSEGLIO's avatar Gaëlle USSEGLIO
Browse files

ENH : Small chgts into SARAmbigAlt application

parent 87d5231d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 120 additions and 63 deletions
app/otbSARAltAmbig.cxx
+ 120
− 63
View file @ 9e23ee8c
......@@ -28,6 +28,7 @@
#include <iomanip>
#include <string>
#include <fstream>
#include <vector>
namespace otb
{
......@@ -48,12 +49,12 @@ private:
void DoInit() override
{
SetName("SARAltAmbig");
SetDescription("Evaluation of altitudes of ambiguity of orbits");
SetDescription("Evaluation of ambiguity height");
SetDocLongDescription("Calculate the angle alpha between "
"ground-master and ground-slave vectors. "
"Calculates the local incidence.\nThe ambiguity "
"altitude is given by: "
"Calculates the local incidence.\n The ambiguity "
"heigth is given by: "
"lambda * sin (incidence) / 2 * tg (alpha)");
//Optional descriptors
......@@ -63,7 +64,7 @@ private:
AddDocTag("DiapOTB");
//Parameter declarations
AddParameter(ParameterType_InputImage, "inmaster", "Input Master geosar image");
AddParameter(ParameterType_InputImage, "inmaster", "Input Master image");
SetParameterDescription("inmaster", "SAR Image master");
AddParameter(ParameterType_Float, "lat", "Tartget lat");
......@@ -81,7 +82,7 @@ private:
SetMinimumParameterFloatValue("height", -100.);
SetMaximumParameterFloatValue("height", 9000.);
AddParameter(ParameterType_InputImage, "inslave", "Input Slave geosar image");
AddParameter(ParameterType_InputImage, "inslave", "Input Slave image");
SetParameterDescription("inslave", "SAR image slave");
AddParameter(ParameterType_OutputFilename, "outfile", "Output file to store the results");
......@@ -115,27 +116,24 @@ private:
// 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);
// Get image keywordlists
otb::ImageKeywordlist sarKWL = SARPtr->GetImageKeywordlist();
otb::ImageKeywordlist sarKWL_slave = SARPtr_slave->GetImageKeywordlist();
......@@ -147,56 +145,94 @@ private:
SarSensorModelAdapter::Point3DType xyzCart;
SarSensorModelAdapter::Point3DType satellitePosition[2];
SarSensorModelAdapter::Point3DType satelliteVelocity[2];
double incidence[2];
demGeoPoint[0] = lon;
demGeoPoint[1] = lat;
demGeoPoint[2] = height;
// Results allocations
double * resultsLocalVector = new double[3];
double * resultsAltAmbig = new double[4];
// Cartesian conversion
SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
otbAppLogINFO(<<"Cartesian coords : " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2]);
output << "GROUND COORDINATES = " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2] << std::endl;
otbAppLogDEBUG(<<"Cartesian coords : " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2]);
// Image names
std::string master = getFilename(GetParameterString("inmaster"));
std::string slave = getFilename(GetParameterString("inslave"));
std::vector<std::string> imgNames;
imgNames.push_back(master);
imgNames.push_back(slave);
// Satellite position and velocity
// Master
otbAppLogINFO(<<"Inmaster Image : ");
otbAppLogDEBUG(<<"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]);
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);
otbAppLogDEBUG(<<"Target Image : ");
m_SarSensorModelAdapter->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());
const double C = 299792458.;
float lambda = C/radarFreq;
otbAppLogINFO(<<"Lambda : " << lambda);
otbAppLogDEBUG(<<"Lambda : " << lambda);
// Estimations : mean incidence, ambig_heigth, radial coord, lateral coord and local projections for master
otbAppLogDEBUG(<< "Compute Alt Ambig");
computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, imgNames, resultsLocalVector, resultsAltAmbig);
// Write into output file
// Open the output file
std::ofstream output(GetParameterString("outfile"), std::ios::out | std::ios::trunc);
output << "Latitude of the target = " << lat << std::endl;
output << "Longitude of the target = " << lon << std::endl;
output << "Altitude of the target = " << height << std::endl;
output << "GROUND COORDINATES = " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2] << std::endl;
// 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(70) << "MASTER IMAGE"
<< std::setw(70) << "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);
output << std::left
<< std::setw(70) << master
<< std::setw(70) << slave
<< std::setw(12) << resultsAltAmbig[0]
<< std::setw(12) << resultsAltAmbig[1]
<< std::setw(12) <<resultsAltAmbig[2]
<< std::setw(12) << resultsAltAmbig[3] << std::endl;
output << "Ground target --> Satellite Vector" << std::endl;
output << "North Projection : " << resultsLocalVector[0] << std::endl;
output << "East Projection : " << resultsLocalVector[1] << std::endl;
output << "Vertical Projection : " << resultsLocalVector[2] << std::endl;
// Close the output file
output.close();
// Free Memory
delete resultsLocalVector;
delete resultsAltAmbig;
resultsLocalVector = 0;
resultsAltAmbig = 0;
}
std::string getFilename(std::string filePath){
......@@ -205,7 +241,13 @@ private:
return filePath.substr(sepPos + 1, dotPos - sepPos - 1);
}
double computeIncidence(SarSensorModelAdapter::Point3DType satellitePosition, SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType rm){
///////// compute functions /////////
// Incidence
double computeIncidence(SarSensorModelAdapter::Point3DType satellitePosition,
SarSensorModelAdapter::Point3DType xyzCart,
SarSensorModelAdapter::Point3DType rm)
{
// Incidence
double normeS = sqrt(satellitePosition[0]*satellitePosition[0] +
satellitePosition[1]*satellitePosition[1] +
......@@ -218,13 +260,20 @@ private:
return incidence;
}
void searchPassageCloserToGround(SarSensorModelAdapter::Point3DType satellitePosition, SarSensorModelAdapter::Point3DType xyzCart, std::ofstream& output, SarSensorModelAdapter::Point3DType rm){
// Local Vector
void computeLocalVector(SarSensorModelAdapter::Point3DType satellitePosition,
SarSensorModelAdapter::Point3DType xyzCart,
SarSensorModelAdapter::Point3DType rm,
double * resultsLocalVector)
{
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;
// lon an lat into radians
double lon = GetParameterFloat("lon")*M_PI/180.;
double lat = GetParameterFloat("lat")*M_PI/180.;
......@@ -241,19 +290,27 @@ private:
vectk[2] = sin(lat);
// Rm projection
output << "Ground target --> Satellite Vector" << std::endl;
// Use -Rm because the given rm, here, is egal to : -(xyzCart - satellitePosition)
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;
otbAppLogDEBUG(<< "North Projection : " << northProjection);
otbAppLogDEBUG(<< "East Projection : " << eastProjection);
otbAppLogDEBUG(<< "Vertical Projection : " << verticalProjection);
// Assign output
resultsLocalVector[0] = northProjection;
resultsLocalVector[1] = eastProjection;
resultsLocalVector[2] = verticalProjection;
}
void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType* satellitePosition, SarSensorModelAdapter::Point3DType* satelliteVelocity, double lambda, std::ofstream& output){
// Alt Ambig
void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart,
SarSensorModelAdapter::Point3DType* satellitePosition,
SarSensorModelAdapter::Point3DType* satelliteVelocity,
double lambda, std::vector<std::string> imgNames,
double * resultsLocalVector, 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;
......@@ -313,24 +370,24 @@ private:
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::string master = imgNames[0];
std::string slave = imgNames[1];
otbAppLogDEBUG(<< std::left
<< std::setw(70) << master
<< std::setw(70) << 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);
// Assign output
results[0] = inc_moy * 180. / M_PI; // Mean incidence (in degrees)
results[1] = alt_amb; // Ambiguity height
results[2] = xe; // Radial coordinate
results[3] = ye / factorBperp; // Lateral coordinate
// Local vector estimation (for master only)
computeLocalVector(satellitePosition[0], xyzCart, rm, resultsLocalVector);
}
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment