diff --git a/Modules/Core/Metadata/src/otbDimapMetadataHelper.cxx b/Modules/Core/Metadata/src/otbDimapMetadataHelper.cxx
index 90573912f4cf2ec5b107c6bafc085def912aa2ff..4cd68828143eedbe1ac488c50254e66508d6b50c 100644
--- a/Modules/Core/Metadata/src/otbDimapMetadataHelper.cxx
+++ b/Modules/Core/Metadata/src/otbDimapMetadataHelper.cxx
@@ -338,7 +338,6 @@ void DimapMetadataHelper::ParseDimapV2(const MetadataSupplierInterface & mds, co
void DimapMetadataHelper::ParseSpot5Model(const MetadataSupplierInterface & mds, Spot5Param& spot5Param, const std::string & prefix){
using Point3DType = itk::Point<double, 3>;
- using Point2DType = itk::Point<double, 2>;
std::vector<double> yaw_vector;
std::vector<double> pitch_vector;
@@ -390,7 +389,7 @@ void DimapMetadataHelper::ParseSpot5Model(const MetadataSupplierInterface & mds,
ParseVector(mds, expr_attitude,"ROLL", roll_vector);
ParseVector(mds, expr_attitude,"TIME", time_vector);
- for (int i=0; i < yaw_vector.size(); i++){
+ for (size_t i=0; i < yaw_vector.size(); i++){
Point3DType point3d;
point3d[0] = pitch_vector[i];
point3d[1] = roll_vector[i];
@@ -410,15 +409,15 @@ void DimapMetadataHelper::ParseSpot5Model(const MetadataSupplierInterface & mds,
// Use look angles from Green band
// /!\ Warning check condition with SWIR band not clear in OSSIM!
bool bandFound = false;
- int i = 1;
std::string expr;
hasValue = false;
if (m_Data.BandIDs.size()>1){
- while (i < m_Data.BandIDs.size() && !bandFound ){
+ size_t i = 1;
+ while (i < m_Data.BandIDs.size() && !bandFound) {
expr = prefix + "Data_Strip.Sensor_Configuration.Instrument_Look_Angles_List.Instrument_Look_Angles_"+std::to_string(i)+".BAND_INDEX";
- mds.GetMetadataValue(expr, hasValue) == "2" ? bandFound=true:i++;
+ mds.GetMetadataValue(expr, hasValue) == "2" ? bandFound=true : i++;
}
expr = "Dimap_Document.Data_Strip.Sensor_Configuration.Instrument_Look_Angles_List.Instrument_Look_Angles_"+std::to_string(i)+".Look_Angles_List.Look_Angles";
}
@@ -441,7 +440,7 @@ void DimapMetadataHelper::ParseSpot5Model(const MetadataSupplierInterface & mds,
ParseVector(mds, expr,"Velocity.Z", vel_z);
ParseVector(mds, expr,"TIME", time_vector);
- for (int i=0; i < pos_x.size(); i++){
+ for (size_t i = 0; i < pos_x.size(); i++) {
Point3DType position;
Point3DType velocity;
position[0] = pos_x[i];
diff --git a/Modules/Core/Transform/include/otbBilinearProjection.h b/Modules/Core/Transform/include/otbBilinearProjection.h
index b20e040a4c65347e26fd83913462804fbb2aef81..6e15273c842fe04ef5f0b6e5d8b075f962f2d9ac 100644
--- a/Modules/Core/Transform/include/otbBilinearProjection.h
+++ b/Modules/Core/Transform/include/otbBilinearProjection.h
@@ -54,49 +54,52 @@ public:
/** Run-time type information (and related methods). */
itkTypeMacro(BilinearProjection, Object);
- void worldToLineSample(const Point3DType& worldPoint,
- Point2DType& lineSampPt) const;
- /*!
- * METHOD: lineSampleToWorld()
- * Performs the inverse projection from line, sample to ground (world):
- */
- void lineSampleToWorld(const Point2DType& lineSampPt,
- Point3DType& worldPt) const;
+ Point2DType worldToLineSample(const Point3DType& worldPoint) const;
+ /*!
+ * METHOD: lineSampleToWorld()
+ * Performs the inverse projection from line, sample to ground (world):
+ * @return line sample to 3D pos at height 0.0
+ */
+ Point3DType lineSampleToWorld(Point2DType lineSampPt) const;
- /*!
- * METHOD: lineSampleHeightToWorld
- * This projects the image point to the given
- * elevation above ellipsoid, thereby bypassing reference to a DEM. Useful
- * for projections that are sensitive to elevation (such as sensor models).
- */
- void lineSampleHeightToWorld(const Point2DType& lineSampPt,
- const double& heightAboveEllipsoid,
- Point3DType& worldPt) const;
+ /*!
+ * METHOD: lineSampleHeightToWorld
+ * This projects the image point to the given
+ * elevation above ellipsoid, thereby bypassing reference to a DEM. Useful
+ * for projections that are sensitive to elevation (such as sensor models).
+ *
+ * @return image point to world point at given height
+ */
+ Point3DType lineSampleHeightToWorld(Point2DType lineSampPt,
+ const double& heightAboveEllipsoid) const;
- void setTiePoints(const std::vector<Point2DType>& lsPt,
- const std::vector<Point3DType>& geoPt);
+ const std::vector<Point2DType>& getLineSamplePoints() const;
- void getTiePoints(std::vector<Point2DType>& lsPt, std::vector<Point3DType>& geoPt) const;
+ void setLineSamplePoints(const std::vector<Point2DType>& lsPt);
- /** Return true if there is any nan in the points coordinates */
- bool imgPointsHaveNan();
+ const std::vector<Point3DType>& getWorldPoints() const;
- /** Return true if there is any nan in the points coordinates */
- bool worldPointsHaveNan();
+ void setWorldPoints(const std::vector<Point3DType>& wPt);
- /** Resolve the bilinear system */
- void computeLS();
+ /** Return true if there is any nan in the points coordinates */
+ bool imgPointsHaveNan();
+
+ /** Return true if there is any nan in the points coordinates */
+ bool worldPointsHaveNan();
+
+ /** Resolve the bilinear system */
+ void computeLS();
protected:
BilinearProjection();
BilinearProjection(const Point2DType& ul,
- const Point2DType& ur,
- const Point2DType& lr,
- const Point2DType& ll,
- const Point3DType& ulg,
- const Point3DType& urg,
- const Point3DType& lrg,
- const Point3DType& llg);
+ const Point2DType& ur,
+ const Point2DType& lr,
+ const Point2DType& ll,
+ const Point3DType& ulg,
+ const Point3DType& urg,
+ const Point3DType& lrg,
+ const Point3DType& llg);
virtual ~BilinearProjection() = default;
private:
diff --git a/Modules/Core/Transform/include/otbSpot5ForwardTransform.hxx b/Modules/Core/Transform/include/otbSpot5ForwardTransform.hxx
index 19d5ef4953d3e5ab8e725823cd27eb041f3c6327..cfbd7f8007e73319700d2e47e2c6f0e591868f49 100644
--- a/Modules/Core/Transform/include/otbSpot5ForwardTransform.hxx
+++ b/Modules/Core/Transform/include/otbSpot5ForwardTransform.hxx
@@ -39,9 +39,9 @@ Spot5ForwardTransform<TScalarType, NInputDimensions, NOutputDimensions>::Transfo
sensorPoint[0] = static_cast<double>(point[0]) - 0.5;
sensorPoint[1] = static_cast<double>(point[1]) - 0.5;
if (NInputDimensions > 2)
- this->m_Transformer->LineSampleHeightToWorld(sensorPoint, static_cast<double>(point[2]), worldPoint);
+ worldPoint = this->m_Transformer->LineSampleHeightToWorld(sensorPoint, static_cast<double>(point[2]));
else
- this->m_Transformer->LineSampleToWorld(sensorPoint, worldPoint);
+ worldPoint = this->m_Transformer->LineSampleToWorld(sensorPoint);
OutputPointType pOut;
pOut[0] = static_cast<TScalarType>(worldPoint[0]);
diff --git a/Modules/Core/Transform/include/otbSpot5InverseTransform.hxx b/Modules/Core/Transform/include/otbSpot5InverseTransform.hxx
index eb9e3a1f38bbb00d3b77fb1f2203dc2fef13ddd5..ab291e1be25d49b55adcb640df09adafcf6c902e 100644
--- a/Modules/Core/Transform/include/otbSpot5InverseTransform.hxx
+++ b/Modules/Core/Transform/include/otbSpot5InverseTransform.hxx
@@ -34,7 +34,6 @@ template <class TScalarType, unsigned int NInputDimensions, unsigned int NOutput
typename Spot5InverseTransform<TScalarType, NInputDimensions, NOutputDimensions>::OutputPointType
Spot5InverseTransform<TScalarType, NInputDimensions, NOutputDimensions>::TransformPoint(const InputPointType& point) const
{
- Spot5SensorModel::Point2DType sensorPoint;
Spot5SensorModel::Point3DType worldPoint;
worldPoint[0] = static_cast<double>(point[0]);
worldPoint[1] = static_cast<double>(point[1]);
@@ -43,7 +42,7 @@ Spot5InverseTransform<TScalarType, NInputDimensions, NOutputDimensions>::Transfo
else
worldPoint[2] = otb::DEMHandler::GetInstance().GetHeightAboveEllipsoid(point[0], point[1]);
- this->m_Transformer->WorldToLineSample(worldPoint, sensorPoint);
+ Spot5SensorModel::Point2DType sensorPoint(this->m_Transformer->WorldToLineSample(worldPoint));
OutputPointType pOut;
diff --git a/Modules/Core/Transform/include/otbSpot5SensorModel.h b/Modules/Core/Transform/include/otbSpot5SensorModel.h
index 37aa3068c01267412696da48eac0284af4316020..7f239506329e2d6d9d2de90dfdf0330fbbf2aade 100644
--- a/Modules/Core/Transform/include/otbSpot5SensorModel.h
+++ b/Modules/Core/Transform/include/otbSpot5SensorModel.h
@@ -68,54 +68,52 @@ public:
* @brief Transform world point (lat,lon,hgt) to image point (col,row).
*
* @param[in] inGeoPoint ground point in WGS84 (lat, lon, hgt) coordinates
- * @param[out] outLineSample image point (col,row) coordinates
+ * @return ground point converted to image point (col,row) coordinates
*/
- void WorldToLineSample(const Point3DType& inGeoPoint,
- Point2DType& outLineSample) const;
+ Point2DType WorldToLineSample(const Point3DType& inGeoPoint) const;
/**
* @brief Transform image point (col, row) with a height to world point (lat,lon,hgt).
*
* @param[in] imPt image point (col, row)
* @param[in] heightAboveEllipsoid height
- * @param[out] worldPt world point (lat,lon,hgt)
+ * @return image point and height converted to world point (lat,lon,hgt)
*/
- void LineSampleHeightToWorld(const Point2DType& imPt,
- double heightAboveEllipsoid,
- Point3DType& worldPt) const;
+ Point3DType LineSampleHeightToWorld(Point2DType imPt,
+ double heightAboveEllipsoid) const;
/**
* @brief Transform image point (col, row) to world point (lat,lon,hgt).
*
* @param[in] imPt image point (col, row)
- * @param[out] worldPt world point (lat,lon,hgt)
+ * @return image point converted to world point (lat,lon,hgt)
*/
- void LineSampleToWorld(const Point2DType& imPt, Point3DType& worldPt) const;
+ Point3DType LineSampleToWorld(Point2DType imPt) const;
/**
* @brief Compute a ray from sensor position and image point.
*
* @param[in] imPt image point (col row)
- * @param[out] epPt ephemeris point (3D space point with position and velocity)
+ * @return ephemeris point (3D space point with position and velocity)
*/
- void ImagingRay(const Point2DType& imPt, Ephemeris& epPt) const;
+ Ephemeris ImagingRay(Point2DType imPt) const;
/**
* @brief Compute the nearest intersection (world point) between the ray and the ellipsoid.
*
* @param[in] imRay ephemeris point (3D space point with position and velocity)
* @param[in] offset double offset
- * @param[out] worldPt world point (lat,lon,hgt)
+ * @return world point (lat,lon,hgt)
*/
- void NearestIntersection(const Ephemeris& imRay, const double& offset, Point3DType& worldPt) const;
+ Point3DType NearestIntersection(const Ephemeris& imRay, const double& offset) const;
/**
* @brief Compute world point intersected by image ray.
*
* @param[in] imRay ephemeris point (3D space point with position and velocity)
- * @param[out] worldPt world point (lat,lon,hgt)
+ * @return world point (lat,lon,hgt)
*/
- void IntersectRay(const Ephemeris& imRay, Point3DType& worldPt) const;
+ Point3DType IntersectRay(const Ephemeris& imRay) const;
/**
* @brief Compute the bilinear interpolation from a 3D point vector with a double time vector.
@@ -123,12 +121,11 @@ public:
* @param[in] time time to interpolate
* @param[in] V 3D point vector
* @param[in] T Index time double vector
- * @param[out] li 3D point interpolated
+ * @return 3D point interpolated
*/
- void GetBilinearInterpolation(const double& time,
- const std::vector<Point3DType>& V,
- const std::vector<double>& T,
- Point3DType& li) const;
+ Point3DType GetBilinearInterpolation(const double& time,
+ const std::vector<Point3DType>& V,
+ const std::vector<double>& T) const;
/**
* @brief Compute the lagrange interpolation from a 3D point vector with a double time vector.
@@ -136,28 +133,27 @@ public:
* @param[in] time time to interpolate
* @param[in] V 3D point vector
* @param[in] T Index time double vector
- * @param[out] li 3D point interpolated
+ * @return 3D point interpolated
*/
- void GetLagrangeInterpolation(const double& time,
- const std::vector<Point3DType>& V,
- const std::vector<double>& T,
- Point3DType& li) const;
+ Point3DType GetLagrangeInterpolation(const double& time,
+ const std::vector<Point3DType>& V,
+ const std::vector<double>& T) const;
/**
* @brief Get the 3D position of the sensor at time (interpolation of position samples vector from metadata).
*
* @param[in] time input time
- * @param[out] ecef 3D position
+ * @return 3D position
*/
- void GetPositionEcf(const double& time, Point3DType& ecef) const;
+ Point3DType GetPositionEcf(const double& time) const;
/**
* @brief Get the 3D velocity of the sensor at time (interpolation of velocity samples vector from metadata).
*
* @param[in] time input time
- * @param[out] ecef 3D position
+ * @return 3D position
*/
- void GetVelocityEcf(const double& time, Point3DType& ecef) const;
+ Point3DType GetVelocityEcf(const double& time) const;
/**
* @brief Get look angles on X and Y axis of the sensor at line (interpolation of angles samples vector from metadata).
@@ -172,26 +168,26 @@ public:
* @brief Get the Attitude of the sensor at time (interpolation of attitude samples vector from metadata).
*
* @param[in] time input time
- * @param[out] at 3D attitude
+ * @return 3D attitude
*/
- void GetAttitude(const double& time, Point3DType& at) const;
+ Point3DType GetAttitude(const double& time) const;
/**
* @brief Compute SatToOrb matrix with an input time.
*
- * @param[out] result 3X3 matrix computed
* @param[in] t input time
+ * @return 3X3 matrix computed
*/
- void ComputeSatToOrbRotation(MatrixType& result, double t) const;
+ MatrixType ComputeSatToOrbRotation(double t) const;
+
-
/**
* @brief Update a ImageMetadata object with the stored Spot5Param and GCPs, possibly modified from the
* original metadata by the Spot5SensorModel
*
* @param[in] imd ImageMetadata to be updated
*/
- void UpdateImageMetadata(ImageMetadata & imd);
+ void UpdateImageMetadata(ImageMetadata & imd);
protected:
diff --git a/Modules/Core/Transform/src/otbBilinearProjection.cxx b/Modules/Core/Transform/src/otbBilinearProjection.cxx
index 1a156395db82900f6d6da2aba8efeb901053ec11..1cf8c9d621bc3c13c0ab646a107042616142d2f9 100644
--- a/Modules/Core/Transform/src/otbBilinearProjection.cxx
+++ b/Modules/Core/Transform/src/otbBilinearProjection.cxx
@@ -31,13 +31,13 @@ BilinearProjection::BilinearProjection()
}
BilinearProjection::BilinearProjection(const Point2DType& ul,
- const Point2DType& ur,
- const Point2DType& lr,
- const Point2DType& ll,
- const Point3DType& ulg,
- const Point3DType& urg,
- const Point3DType& lrg,
- const Point3DType& llg)
+ const Point2DType& ur,
+ const Point2DType& lr,
+ const Point2DType& ll,
+ const Point3DType& ulg,
+ const Point3DType& urg,
+ const Point3DType& lrg,
+ const Point3DType& llg)
{
m_LatFit = LSQREstimatorType::New();
m_LonFit = LSQREstimatorType::New();
@@ -56,12 +56,13 @@ BilinearProjection::BilinearProjection(const Point2DType& ul,
m_worldPoints[2] = lrg;
m_worldPoints[3] = llg;
- computeLS();
+ computeLS();
}
-void BilinearProjection::worldToLineSample(const Point3DType& worldPoint,
- Point2DType& lineSampPt) const
+itk::Point<double, 2> BilinearProjection::worldToLineSample(const Point3DType& worldPoint) const
{
+ // convert Point3D to vector double to use
+ // LeastSquareBilinearTransformEstimator::lsFitValue
itk::Vector<double,2> worldptMatrix;
worldptMatrix[0] = worldPoint[0];
worldptMatrix[1] = worldPoint[1];
@@ -70,22 +71,21 @@ void BilinearProjection::worldToLineSample(const Point3DType& worldPoint,
m_XFit->lsFitValue(worldptMatrix, x);
m_YFit->lsFitValue(worldptMatrix, y);
+
+ Point2DType lineSampPt;
lineSampPt[0] = x;
lineSampPt[1] = y;
-
+ return lineSampPt;
}
-void BilinearProjection::lineSampleToWorld(const Point2DType& lineSampPt,
- Point3DType& worldPt) const
+BilinearProjection::Point3DType BilinearProjection::lineSampleToWorld(Point2DType lineSampPt) const
{
- lineSampleHeightToWorld(lineSampPt,
- 0.0,
- worldPt);
+ return lineSampleHeightToWorld(lineSampPt, 0.0);
}
-void BilinearProjection::lineSampleHeightToWorld(const Point2DType& lineSampPt,
- const double& heightAboveEllipsoid,
- Point3DType& worldPt) const
+BilinearProjection::Point3DType BilinearProjection::lineSampleHeightToWorld(
+ Point2DType lineSampPt,
+ const double& heightAboveEllipsoid) const
{
itk::Vector<double,2> lineSampMatrix;
lineSampMatrix[0] = lineSampPt[0];
@@ -93,22 +93,33 @@ void BilinearProjection::lineSampleHeightToWorld(const Point2DType& lineSampPt,
double lon, lat;
m_LonFit->lsFitValue(lineSampMatrix, lon);
m_LatFit->lsFitValue(lineSampMatrix, lat);
+
+ Point3DType worldPt;
worldPt[0] = lon;
worldPt[1] = lat;
worldPt[2] = heightAboveEllipsoid;
+
+ return worldPt;
+}
+
+const std::vector<itk::Point<double, 2>>& BilinearProjection::getLineSamplePoints() const
+{
+ return m_LineSamplePoints;
+}
+
+void BilinearProjection::setLineSamplePoints(const std::vector<Point2DType>& lsPts)
+{
+ m_LineSamplePoints = lsPts;
}
-void BilinearProjection::setTiePoints(const std::vector<Point2DType>& lsPt,
- const std::vector<Point3DType>& geoPt)
+const std::vector<BilinearProjection::Point3DType>& BilinearProjection::getWorldPoints() const
{
- m_LineSamplePoints = lsPt;
- m_worldPoints = geoPt;
+ return m_worldPoints;
}
-void BilinearProjection::getTiePoints(std::vector<Point2DType>& lsPt, std::vector<Point3DType>& geoPt) const
+void BilinearProjection::setWorldPoints(const std::vector<Point3DType>& wPts)
{
- lsPt = m_LineSamplePoints;
- geoPt = m_worldPoints;
+ m_worldPoints = wPts;
}
void BilinearProjection::computeLS()
diff --git a/Modules/Core/Transform/src/otbSpot5SensorModel.cxx b/Modules/Core/Transform/src/otbSpot5SensorModel.cxx
index 15e12272729dbc5c788dedfd62cadd5fb08df268..064c91efd148506fb1ddd9ff27d8ace13d81e5ef 100644
--- a/Modules/Core/Transform/src/otbSpot5SensorModel.cxx
+++ b/Modules/Core/Transform/src/otbSpot5SensorModel.cxx
@@ -49,13 +49,11 @@ Spot5SensorModel::Spot5SensorModel(const ImageMetadata & imd)
void Spot5SensorModel::InitBilinearTransform(){
Point2DType ul, ur, lr, ll;
- Point3DType ulg, urg, lrg, llg;
ul[0] = 0; ul[1] = 0;
ur[0] = m_Spot5Param.ImageSize[0] - 1; ur[1] = 0;
ll[0] = 0; ll[1] = m_Spot5Param.ImageSize[1] - 1;
- lr[0] = m_Spot5Param.ImageSize[0] - 1;
- lr[1] = m_Spot5Param.ImageSize[1] - 1;
+ lr[0] = m_Spot5Param.ImageSize[0] - 1; lr[1] = m_Spot5Param.ImageSize[1] - 1;
m_ImageRect = PolygonType::New();
m_ImageRect->AddVertex(Point2DToIndex(ul));
@@ -63,10 +61,10 @@ void Spot5SensorModel::InitBilinearTransform(){
m_ImageRect->AddVertex(Point2DToIndex(lr));
m_ImageRect->AddVertex(Point2DToIndex(ll));
- LineSampleToWorld(ul, ulg);
- LineSampleToWorld(ur, urg);
- LineSampleToWorld(lr, lrg);
- LineSampleToWorld(ll, llg);
+ const Point3DType ulg(LineSampleToWorld(ul));
+ const Point3DType urg(LineSampleToWorld(ur));
+ const Point3DType lrg(LineSampleToWorld(lr));
+ const Point3DType llg(LineSampleToWorld(ll));
m_GroundRect = PolygonType::New();
m_GroundRect->AddVertex(Point3DToIndex(ulg));
@@ -77,9 +75,9 @@ void Spot5SensorModel::InitBilinearTransform(){
std::vector<Point2DType> impts {ul,ur,lr,ll};
std::vector<Point3DType> wrldpts {ulg,urg,lrg,llg};
m_BilinearProj = BilinearProjection::New();
- m_BilinearProj->setTiePoints(impts,wrldpts);
+ m_BilinearProj->setLineSamplePoints(impts);
+ m_BilinearProj->setWorldPoints(wrldpts);
m_BilinearProj->computeLS();
-
}
Spot5SensorModel::ContinuousIndexType
@@ -98,7 +96,7 @@ Spot5SensorModel::Point3DToIndex(const itk::Point<double, 3> point) const{
return indexPoint;
}
-void Spot5SensorModel::WorldToLineSample(const Point3DType& inGeoPoint, Point2DType& outLineSample) const
+itk::Point<double, 2> Spot5SensorModel::WorldToLineSample(const Point3DType& inGeoPoint) const
{
static const double PIXEL_THRESHOLD = .1; // acceptable pixel error
@@ -111,9 +109,10 @@ void Spot5SensorModel::WorldToLineSample(const Point3DType& inGeoPoint, Point2DT
int iters = 0;
Point3DType wdp (inGeoPoint);
- if (!(m_GroundRect->IsInside(Point3DToIndex(wdp))||m_GroundRect->IsOnEdge(Point3DToIndex(wdp)))){
- m_BilinearProj->worldToLineSample(wdp, outLineSample);
- return;
+ if (!(m_GroundRect->IsInside(Point3DToIndex(wdp)) ||
+ m_GroundRect->IsOnEdge(Point3DToIndex(wdp))))
+ {
+ return m_BilinearProj->worldToLineSample(wdp);
}
@@ -121,8 +120,7 @@ void Spot5SensorModel::WorldToLineSample(const Point3DType& inGeoPoint, Point2DT
// Initialize with binilear interpolation
- Point2DType ip;
- m_BilinearProj->worldToLineSample(wdp, ip);
+ Point2DType ip(m_BilinearProj->worldToLineSample(wdp));
Point2DType ip_du;
Point2DType ip_dv;
@@ -148,9 +146,9 @@ void Spot5SensorModel::WorldToLineSample(const Point3DType& inGeoPoint, Point2DT
//***
// Compute numerical partials at current guessed point:
//***
- LineSampleHeightToWorld(ip, height, gp);
- LineSampleHeightToWorld(ip_du, height, gp_du);
- LineSampleHeightToWorld(ip_dv, height, gp_dv);
+ gp = LineSampleHeightToWorld(ip, height);
+ gp_du = LineSampleHeightToWorld(ip_du, height);
+ gp_dv = LineSampleHeightToWorld(ip_dv, height);
dlat_du = gp_du[1] - gp[1]; //e
dlon_du = gp_du[0] - gp[0]; //g
@@ -187,87 +185,71 @@ void Spot5SensorModel::WorldToLineSample(const Point3DType& inGeoPoint, Point2DT
} while (iters < MAX_NUM_ITERATIONS); /* && ((std::abs(delta_u) < PIXEL_THRESHOLD)
&& (std::abs(delta_v) < PIXEL_THRESHOLD)) */
- outLineSample = ip - m_Spot5Param.SubImageOffset;
+ return ip - m_Spot5Param.SubImageOffset;
}
-void Spot5SensorModel::LineSampleHeightToWorld(const Point2DType& imPt,
- double heightAboveEllipsoid,
- Point3DType& worldPt) const
+itk::Point<double, 3> Spot5SensorModel::LineSampleHeightToWorld(Point2DType imPt,
+ double heightAboveEllipsoid) const
{
// Outside the image
if (!insideImage(imPt, 1.0 - FLT_EPSILON))
{
- m_BilinearProj->lineSampleToWorld(imPt, worldPt);
- return;
+ return m_BilinearProj->lineSampleToWorld(imPt);
}
- Ephemeris ray;
- ImagingRay(imPt, ray);
-
// Ecef coordinate point intersection
- Point3DType intersectPt;
- NearestIntersection(ray, heightAboveEllipsoid, intersectPt);
-
// Conversion to lat long coordinates
- worldPt = EcefToWorld(intersectPt);
-
+ return EcefToWorld(NearestIntersection(ImagingRay(imPt), heightAboveEllipsoid));
}
-void Spot5SensorModel::LineSampleToWorld(const Point2DType& imPt,
- Point3DType& worldPt) const
+itk::Point<double, 3> Spot5SensorModel::LineSampleToWorld(Point2DType imPt) const
{
if (!insideImage(imPt, 2.0))
{
- m_BilinearProj->lineSampleToWorld(imPt, worldPt);
- return;
+ return m_BilinearProj->lineSampleToWorld(imPt);
}
// Determine ray and intersect with terrain model
- Ephemeris ray;
- ImagingRay(imPt, ray);
- IntersectRay(ray, worldPt);
-
+ return IntersectRay(ImagingRay(imPt));
}
-void Spot5SensorModel::GetBilinearInterpolation(
- const double& time,
- const std::vector<Point3DType>& V,
- const std::vector<double>& T,
- Point3DType& li) const
+itk::Point<double, 3> Spot5SensorModel::GetBilinearInterpolation(
+ const double& time,
+ const std::vector<Point3DType>& V,
+ const std::vector<double>& T) const
{
// get iterator to first element >= time
const auto first_sup_value = std::lower_bound(T.begin(), T.end(), time);
if(first_sup_value == T.begin())
{
- li = V[0];
+ return V[0];
}
else if(first_sup_value == T.end())
{
- li = V[1];
+ return V[1];
}
else
{
const int32_t samp0 = std::distance(T.begin(), first_sup_value);
double t = (*(first_sup_value-1)-time)/
(*(first_sup_value-1) - (*first_sup_value));
- const auto it_v = V.begin() + samp0;
+ const auto& it_v = V.begin() + samp0;
- li = *(it_v-1) + ((*it_v)-*(it_v-1))*t;
+ return *(it_v-1) + ((*it_v)-*(it_v-1))*t;
}
}
-void Spot5SensorModel::GetLagrangeInterpolation(
- const double& time,
- const std::vector<Point3DType>& V,
- const std::vector<double>& T,
- Point3DType& li) const
+itk::Point<double, 3> Spot5SensorModel::GetLagrangeInterpolation(
+ const double& time,
+ const std::vector<Point3DType>& V,
+ const std::vector<double>& T) const
{
//
// Verify that t is within allowable range:
@@ -348,37 +330,34 @@ void Spot5SensorModel::GetLagrangeInterpolation(
S[2] += p[2];
}
- li = S;
-
+ return S;
}
-void Spot5SensorModel::GetPositionEcf(const double& time,
- Point3DType& ecef) const
+itk::Point<double, 3> Spot5SensorModel::GetPositionEcf(const double& time) const
{
if((m_Spot5Param.EcefPosSamples.size() < 8)||
(m_Spot5Param.EcefTimeSamples.size() < 8))
{
- GetBilinearInterpolation(time, m_Spot5Param.EcefPosSamples, m_Spot5Param.EcefTimeSamples, ecef);
+ return GetBilinearInterpolation(time, m_Spot5Param.EcefPosSamples, m_Spot5Param.EcefTimeSamples);
}
else
{
- GetLagrangeInterpolation(time, m_Spot5Param.EcefPosSamples, m_Spot5Param.EcefTimeSamples, ecef);
+ return GetLagrangeInterpolation(time, m_Spot5Param.EcefPosSamples, m_Spot5Param.EcefTimeSamples);
}
}
-void Spot5SensorModel::GetVelocityEcf(const double& time,
- Point3DType& ecef) const
+itk::Point<double, 3> Spot5SensorModel::GetVelocityEcf(const double& time) const
{
if((m_Spot5Param.EcefVelSamples.size() < 8)||
(m_Spot5Param.EcefTimeSamples.size() < 8))
{
- GetBilinearInterpolation(time, m_Spot5Param.EcefVelSamples, m_Spot5Param.EcefTimeSamples, ecef);
+ return GetBilinearInterpolation(time, m_Spot5Param.EcefVelSamples, m_Spot5Param.EcefTimeSamples);
}
else
{
- GetLagrangeInterpolation(time, m_Spot5Param.EcefVelSamples, m_Spot5Param.EcefTimeSamples, ecef);
+ return GetLagrangeInterpolation(time, m_Spot5Param.EcefVelSamples, m_Spot5Param.EcefTimeSamples);
}
}
@@ -398,12 +377,12 @@ void Spot5SensorModel::GetPixelLookAngleXY(unsigned int line,
}
-void Spot5SensorModel::GetAttitude(const double& time, Point3DType& at) const
+itk::Point<double, 3> Spot5SensorModel::GetAttitude(const double& time) const
{
- GetBilinearInterpolation(time, m_Spot5Param.AttitudesSamples, m_Spot5Param.AttitudesSamplesTimes, at);
+ return GetBilinearInterpolation(time, m_Spot5Param.AttitudesSamples, m_Spot5Param.AttitudesSamplesTimes);
}
-void Spot5SensorModel::NearestIntersection(const Ephemeris& imRay, const double& offset, Point3DType& worldPt) const
+itk::Point<double, 3> Spot5SensorModel::NearestIntersection(const Ephemeris& imRay, const double& offset) const
{
// WGS 84 parameters conversion
double wgsA = 6378137.000;
@@ -453,11 +432,11 @@ void Spot5SensorModel::NearestIntersection(const Ephemeris& imRay, const double&
t = t2;
}
- worldPt = imRay.position + imRay.velocity * t;
+ return imRay.position + imRay.velocity * t;
}
-void Spot5SensorModel::IntersectRay(const Ephemeris& imRay, Point3DType& worldPt) const
+itk::Point<double, 3> Spot5SensorModel::IntersectRay(const Ephemeris& imRay) const
{
static const double CONVERGENCE_THRESHOLD = 0.001; // meters
@@ -469,17 +448,17 @@ void Spot5SensorModel::IntersectRay(const Ephemeris& imRay, Point3DType& worldPt
double distance;
int iteration_count = 0;
- worldPt = EcefToWorld(prev_intersect_pt);
+ Point3DType worldPt(EcefToWorld(prev_intersect_pt));
// Loop to iterate on ray intersection with variable elevation surface:
do
{
-
+
// Intersect ray with ellipsoid inflated by h_ellips:
h_ellips = DEMHandler::GetInstance().GetHeightAboveEllipsoid(worldPt[0], worldPt[1]);
- NearestIntersection(imRay, h_ellips, new_intersect_pt);
+ new_intersect_pt = NearestIntersection(imRay, h_ellips);
// Assign the ground point to the latest iteration's intersection point:
worldPt = EcefToWorld(new_intersect_pt);
@@ -502,14 +481,14 @@ void Spot5SensorModel::IntersectRay(const Ephemeris& imRay, Point3DType& worldPt
<< "point. Result is probably inaccurate." << std::endl);
}
+ return worldPt;
}
-void Spot5SensorModel::ComputeSatToOrbRotation(MatrixType& result, double t) const
+itk::Matrix<double, 3, 3> Spot5SensorModel::ComputeSatToOrbRotation(double t) const
{
-
+ itk::Matrix<double, 3, 3> result;
/* Linearly interpolate attitudes angles: */
- Point3DType att;
- GetAttitude(t, att);
+ Point3DType att(GetAttitude(t));
/* Compute trig functions to populate rotation matrices: ANGLES IN RADIANS */
double cp = cos(att[0]);
@@ -529,12 +508,12 @@ void Spot5SensorModel::ComputeSatToOrbRotation(MatrixType& result, double t) con
result(2, 0) = (-sp*sy+cp*sr*cy);
result(2, 1) = (-sp*cy-cp*sr*sy);
result(2, 2) = cp*cr;
+
+ return result;
}
-void Spot5SensorModel::ImagingRay(const Point2DType& imPt, Ephemeris& imRay) const
+Ephemeris Spot5SensorModel::ImagingRay(Point2DType imPt) const
{
- Point3DType sensorPos;
- Point3DType sensorVel;
Point2DType imPtTmp;
imPtTmp[0] = imPt[0] + m_Spot5Param.SubImageOffset[0];
@@ -544,11 +523,8 @@ void Spot5SensorModel::ImagingRay(const Point2DType& imPt, Ephemeris& imRay) co
double timeLine = m_Spot5Param.RefLineTime + m_Spot5Param.LineSamplingPeriod*(imPtTmp[1] - m_Spot5Param.RefLineTimeLine);
/* 2. Interpolate ephemeris position and velocity (in ECF): */
- Point3DType ecfPos;
- Point3DType ecfVel;
-
- GetPositionEcf(timeLine, ecfPos);
- GetVelocityEcf(timeLine, ecfVel);
+ Point3DType ecfPos(GetPositionEcf(timeLine));
+ Point3DType ecfVel(GetVelocityEcf(timeLine));
/* 3. Establish the look direction in Vehicle LSR space (S_sat). ANGLES IN RADIANS*/
double psiX, psiY;
@@ -560,10 +536,7 @@ void Spot5SensorModel::ImagingRay(const Point2DType& imPt, Ephemeris& imRay) co
u_sat[2] = -1.0;
/* 4. Transform vehicle LSR space look direction vector to orbital LSR space (S_orb): */
- MatrixType satToOrbit;
- ComputeSatToOrbRotation(satToOrbit, timeLine);
-
- Vector3DType uOrb = satToOrbit*u_sat;
+ Vector3DType uOrb = ComputeSatToOrbRotation(timeLine) * u_sat;
uOrb.Normalize();
// 5. Transform orbital LSR space look direction vector to ECF.
@@ -601,8 +574,8 @@ void Spot5SensorModel::ImagingRay(const Point2DType& imPt, Ephemeris& imRay) co
Vector3DType uEcf = (orbToEcfRotation*uOrb);
/* Establish the imaging ray given direction and origin: */
- imRay.position = ecfPos;
- imRay.velocity = uEcf;
+ // Ephemeris.time , Ephemeris.position, Ephemeris.velocity
+ return {0, ecfPos, uEcf};
}