Skip to content
Snippets Groups Projects
Commit 1fdfd58e authored by Valentin Genin's avatar Valentin Genin
Browse files

WIP : Compute Alt Ambig

parent 9642c56b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 98 additions and 10 deletions
app/otbSARAltAmbig.cxx
+ 98
− 10
View file @ 1fdfd58e
......@@ -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);
}
......
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