From 1fdfd58e20d7f77d2936507771ecdba33f3365fb Mon Sep 17 00:00:00 2001
From: Valentin Genin <valentin.genin@cnes.fr>
Date: Mon, 27 Apr 2020 13:59:19 +0000
Subject: [PATCH] WIP : Compute Alt Ambig
---
app/otbSARAltAmbig.cxx | 108 +++++++++++++++++++++++++++++++++++++----
1 file changed, 98 insertions(+), 10 deletions(-)
diff --git a/app/otbSARAltAmbig.cxx b/app/otbSARAltAmbig.cxx
index 1ae050a..97e0b37 100644
--- a/app/otbSARAltAmbig.cxx
+++ b/app/otbSARAltAmbig.cxx
@@ -97,6 +97,9 @@ private:
// Start the first pipeline (Read SAR image metadata)
ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("inmaster");
+ // Start the target pipeline (Read SAR image metadata)
+ ComplexFloatImageType::Pointer SARPtr_target = GetParameterComplexFloatImage("target");
+
// Get lat, lon and height
float lat = GetParameterFloat("lat");
otbAppLogINFO(<<"Lat : "<<lat);
@@ -108,6 +111,7 @@ private:
otbAppLogINFO(<<"Height : "<<height);
otb::ImageKeywordlist sarKWL = SARPtr->GetImageKeywordlist();
+ otb::ImageKeywordlist sarKWL_target = SARPtr_target->GetImageKeywordlist();
// Create and Initilaze the SarSensorModelAdapter
SarSensorModelAdapter::Pointer m_SarSensorModelAdapter = SarSensorModelAdapter::New();
@@ -115,9 +119,9 @@ private:
SarSensorModelAdapter::Point3DType demGeoPoint;
SarSensorModelAdapter::Point3DType xyzCart;
- SarSensorModelAdapter::Point3DType satellitePosition;
- SarSensorModelAdapter::Point3DType satelliteVelocity;
- SarSensorModelAdapter::Point3DType R;
+ SarSensorModelAdapter::Point3DType satellitePosition[2];
+ SarSensorModelAdapter::Point3DType satelliteVelocity[2];
+ double incidence[2];
demGeoPoint[0] = lon;
demGeoPoint[1] = lat;
@@ -125,10 +129,15 @@ private:
SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
otbAppLogINFO(<<"Cartesian coords : " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2]);
-
- m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition, satelliteVelocity);
- otbAppLogINFO(<<"Satellite Position : " << satellitePosition[0] << " " << satellitePosition[1] << " " << satellitePosition[2]);
- otbAppLogINFO(<<"Satellite Velocity : " << satelliteVelocity[0] << " " << satelliteVelocity[1] << " " << satelliteVelocity[2]);
+
+ 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]);
+
+ searchPassageCloserToGround(satellitePosition[0], satelliteVelocity[0], incidence[0], xyzCart, demGeoPoint);
+ otbAppLogINFO(<<"Incidence : " << incidence[0]);
+
// Calculate lambda
double radarFreq = atof(sarKWL.GetMetadataByKey("support_data.radar_frequency").c_str());
const double C = 299792458.;
@@ -136,7 +145,22 @@ private:
float lambda = C/radarFreq;
otbAppLogINFO(<<"Lambda : " << lambda);
- // Incidence
+ loadOk = m_SarSensorModelAdapter->LoadState(sarKWL_target);
+ 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]);
+
+ searchPassageCloserToGround(satellitePosition[1], satelliteVelocity[1], incidence[1], xyzCart, demGeoPoint);
+ otbAppLogINFO(<<"Incidence : " << incidence[1]);
+
+ computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, incidence);
+
+ }
+
+ void searchPassageCloserToGround(SarSensorModelAdapter::Point3DType& satellitePosition, SarSensorModelAdapter::Point3DType& satelliteVelocity, double& incidence, SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType demGeoPoint){
+ // Incidence
+ SarSensorModelAdapter::Point3DType R;
R[0] = xyzCart[0] - satellitePosition[0];
R[1] = xyzCart[1] - satellitePosition[1];
R[2] = xyzCart[2] - satellitePosition[2];
@@ -147,9 +171,8 @@ private:
double normeCible = sqrt(xyzCart[0] * xyzCart[0] + xyzCart[1] * xyzCart[1] + xyzCart[2] * xyzCart[2]);
double normeR = sqrt(R[0] * R[0] + R[1] * R[1] + R[2] * R[2]);
- double incidence = acos((normeS * normeS - normeR * normeR - normeCible * normeCible) /
+ incidence = acos((normeS * normeS - normeR * normeR - normeCible * normeCible) /
(2 * normeCible * normeR) ) * 180. / M_PI;
- otbAppLogINFO(<<"Incidence : " << incidence);
// Calculate the local vectors
double rm[3];
@@ -161,6 +184,9 @@ private:
double vecti[3];
double vectj[3];
double vectk[3];
+ double lon = demGeoPoint[0];
+ double lat = demGeoPoint[1];
+
vecti[0] = -sin(lat)*cos(lon);
vecti[1] = -sin(lat)*sin(lon);
@@ -178,6 +204,68 @@ private:
otbAppLogINFO(<< -(vecti[0]*rm[0]+vecti[1]*rm[1]+vecti[2]*rm[2])/normeRm);
otbAppLogINFO(<< -(vectj[0]*rm[0]+vectj[1]*rm[1]+vectj[2]*rm[2])/normeRm);
otbAppLogINFO(<< -(vectk[0]*rm[0]+vectk[1]*rm[1]+vectk[2]*rm[2])/normeRm);
+ }
+
+ void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType* satellitePosition, SarSensorModelAdapter::Point3DType* satelliteVelocity, double lambda, double* incidence){
+ otbAppLogINFO(<< "Compute Alt Ambig");
+ // TODO : Get the type of phase (Bistatic by default)
+ double factorBperp = 1.0;
+ double factorHamb = 2.0;
+
+ double rm[3];
+ rm[0] = -(xyzCart[0] - satellitePosition[0][0]);
+ rm[1] = -(xyzCart[1] - satellitePosition[0][1]);
+ rm[0] = -(xyzCart[2] - satellitePosition[0][2]);
+
+ double re[3];
+ re[0] = -(xyzCart[0] - satellitePosition[1][0]);
+ re[1] = -(xyzCart[1] - satellitePosition[1][1]);
+ re[0] = -(xyzCart[2] - satellitePosition[1][2]);
+
+ // 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
+ double vmUnit[3];
+ 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 (incidence[1] > incidence[0]){
+ ye *= -1.0;
+ }
+
+ double alpha = (incidence[1] - incidence[0]) * M_PI / 180.;
+ double inc_moy = ((incidence[1] + incidence[0]) / 2) * M_PI / 180.;
+
+ double alt_amb = lambda * sin(inc_moy) / (factorHamb * tan(alpha));
+ otbAppLogINFO(<< "Incidence\tALT AMBIG\tRadial\tLateral");
+ otbAppLogINFO(<< inc_moy * 180. / M_PI << "\t\t" << alt_amb << "\t\t" << xe << "\t" << ye / factorBperp);
}
--
GitLab