RayGrid

-now sorts angle vectors on initiation
-made certain functions/members private
-can now represent the full 360° of azimuth angle
-GetNeighboringRays() now returns all phi angles if given ray points towards pole

CEquiangularRayDistribution
-added class (derived from RayGrid) for equiangular ray distribution

RayResolutionAdapter
-changes according to changes in RayGrid
Signed-off-by: Philipp Schäfer <philipp.schaefer@akustik.rwth-aachen.de>

include/ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.h

@@ -42,24 +42,31 @@ namespace ITAPropagationPathSim
 			typedef std::vector< RayPtr > RayVector;
 			typedef std::vector<RayVector> RayMatrix;
 
-			RayVector vpRays;
-			RayMatrix vvpRayMatrix;
-
-			VistaVector3D v3SourcePos;
-			std::vector<double> vdThetaDeg;
-			std::vector<double> vdPhiDeg;
+			VistaVector3D v3SourcePos; //< Origin of all rays in this ray grid.
+			std::vector<double> vdThetaDeg; //< Sorted vector of elevation angles in degrees.
+			std::vector<double> vdPhiDeg; //< Sorted vector of azimuth angles in degrees.
+			RayMatrix vvpRayMatrix; //< Matrix of pointers to ray. First index refers to elevation (theta), second to azimuth (phi) angle.
+		private:
+			RayVector vpRays; //< Vector of all rays contained by vppRayMatrix. Call InitRayVectorFromRayMatrix() if vppRayMatrix changes.
+			bool bCircularPhi = false; //< Indicates whether the azimuth vector is considered to be circular (phi covers full 360 degrees) or not.
 		
-		protected:
-			CRayGrid() {}
-			CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg);
+		private:
+			CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool& circularPhi = false);
 		public:
-			//TODO: Make sure that theta and phi are sorted
-			CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg);
+			CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool& circularPhi = false);
 		protected:
+			CRayGrid() {}
 			void InitRayVectorFromRayMatrix();
+		private:
 			int IndexToThetaIndex(const int& idx) const;
 			int IndexToPhiIndex(const int& idx) const;
 			int GetIndex(const std::shared_ptr<CRay>& pRay) const;
+			void GetNeighboringAngleIndices(const std::shared_ptr<CRay>& pRay, std::vector<int>& thetaIdxVec, std::vector<int>& phiIdxVec) const;
+			bool HasPoleDirection(const std::shared_ptr<CRay>& pRay) const;
+
+			//---GET Functions---
+		protected:
+			const std::vector<std::shared_ptr<CRay>>& ConstRayVectorReference() const { return vpRays; }
 		public:
 			std::vector<double> ThetaDeg() const { return vdThetaDeg; }
 			std::vector<double> PhiDeg() const { return vdPhiDeg; }
@@ -70,8 +77,23 @@ namespace ITAPropagationPathSim
 			bool IsEmpty() const { return NTheta() == 0 || NPhi() == 0; }
 			bool Contains(std::shared_ptr<CRay> pRay) const { return std::find(vpRays.cbegin(), vpRays.cend(), pRay) == vpRays.cend(); }
 
-			CRayGrid GetNeighboringRays(std::shared_ptr<CRay> pRay) const;
-			void SetToNeighboringRays(std::shared_ptr<CRay> pRay);
+			//! Returns a ray grid containing the rays surrounding the given ray (including this ray).
+			/**
+			  * Returns an empty ray grid if the given ray is not part of this ray grid.
+			  */
+			CRayGrid GetNeighboringRays(const std::shared_ptr<CRay>& pRay) const;
+			//! Sets the boundaries of this ray grid to the rays surrounding the given ray (including this ray).
+			/**
+			  * If the ray is not part of this ray grid, this will be empty.
+			  */
+			void SetToNeighboringRays(const std::shared_ptr<CRay>& pRay);
+		};
+
+		class ITA_PROPAGATION_PATH_SIM_API CEquiangularRayDistribution : public CRayGrid
+		{
+		public:
+			CEquiangularRayDistribution(const VistaVector3D& sourcePos, const int& nAngles) : CEquiangularRayDistribution(sourcePos, nAngles, nAngles) {}
+			CEquiangularRayDistribution(const VistaVector3D& sourcePos, const int& nTheta, const int& nPhi);
 		};
 	}
 }

src/ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.cpp

@@ -2,7 +2,7 @@
 #include <ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.h>
 
 // ITA includes
-//#include <ITAException.h>
+#include <ITAException.h>
 
 // STD
 #include <cmath>
@@ -12,8 +12,10 @@ using namespace ITAPropagationPathSim::AtmosphericRayTracing;
 
 // --- CONSTRUCTORS ---
 // --------------------
+#pragma region Constructors
 
-CRayGrid::CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg): vvpRayMatrix(rayMatrix), vdThetaDeg(thetaDeg), vdPhiDeg(phiDeg)
+CRayGrid::CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool& circularPhi /*= false*/)
+	: vvpRayMatrix(rayMatrix), vdThetaDeg(thetaDeg), vdPhiDeg(phiDeg), bCircularPhi(circularPhi)
 {
 	InitRayVectorFromRayMatrix();
 
@@ -21,31 +23,105 @@ CRayGrid::CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaD
 		v3SourcePos = vpRays.front()->SourcePoint();
 }
 
-CRayGrid::CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg) : v3SourcePos(sourcePos), vdThetaDeg(thetaDeg), vdPhiDeg(phiDeg)
+CRayGrid::CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool& circularPhi /*= false*/)
+	: v3SourcePos(sourcePos), vdThetaDeg(thetaDeg), vdPhiDeg(phiDeg), bCircularPhi(circularPhi)
 {
+	std::sort(vdPhiDeg.begin(), vdPhiDeg.end());
+	std::sort(vdThetaDeg.begin(), vdThetaDeg.end());
+
+	if (vdThetaDeg.front() < 0.0 || vdThetaDeg.back() > 180.0)
+		ITA_EXCEPT_INVALID_PARAMETER("Values for elevation angle theta are out of bounds [0, 180]");
+	if (vdPhiDeg.front() <= -360.0 || vdPhiDeg.back() >= 360.0)
+		ITA_EXCEPT_INVALID_PARAMETER("Values for azimuth angle phi are out of bounds (-360, 360)");
+	if ( (vdPhiDeg.back() - vdPhiDeg.front()) >= 360.0)
+		ITA_EXCEPT_INVALID_PARAMETER("Delta for azimuth angle phi is >= 360");
+
 	for (int idxTheta = 0; idxTheta < thetaDeg.size(); idxTheta++)
 	{
 		vvpRayMatrix.push_back(std::vector< std::shared_ptr<CRay>>());
+		bool isPoleDirection = (abs(thetaDeg[idxTheta]) < DBL_EPSILON || abs(thetaDeg[idxTheta] - 180.0) < DBL_EPSILON);
 		for (int idxPhi = 0; idxPhi < phiDeg.size(); idxPhi++)
 		{
-			auto ray = std::make_shared<CRay>(sourcePos, thetaDeg[idxTheta], phiDeg[idxPhi]);
-			vvpRayMatrix[idxTheta].push_back(ray);
+			if (isPoleDirection && idxPhi > 0) //At poles: Insert same ray multiple times
+				vvpRayMatrix[idxTheta].push_back( vvpRayMatrix[idxTheta].front() );
+			else
+				vvpRayMatrix[idxTheta].push_back( std::make_shared<CRay>(sourcePos, thetaDeg[idxTheta], phiDeg[idxPhi]) );
 		}
 	}
 	InitRayVectorFromRayMatrix();
 }
 
-// --- PRIVATE HELPER ---
+CEquiangularRayDistribution::CEquiangularRayDistribution(const VistaVector3D& sourcePos, const int& nTheta, const int& nPhi)
+{
+	const double deltaTheta = 180.0 / (nTheta - 1);
+	for (int idx = 0; idx < nTheta - 1; idx++)
+		vdThetaDeg.push_back(idx * deltaTheta);
+	vdThetaDeg.push_back(180.0);
+
+	const double deltaPhi = 360.0 / nPhi;
+	for (int idx = 0; idx < nPhi - 1; idx++)
+		vdPhiDeg.push_back(idx * deltaPhi);
+
+	CRayGrid(sourcePos, vdThetaDeg, vdPhiDeg, true);
+}
+#pragma endregion
+
+
+// --- INITIALIZATION ---
 // ----------------------
+#pragma region Initialization
 
-void ITAPropagationPathSim::AtmosphericRayTracing::CRayGrid::InitRayVectorFromRayMatrix()
+void CRayGrid::InitRayVectorFromRayMatrix()
 {
 	vpRays.clear();
 	for each (RayVector rayVector in vvpRayMatrix)
 		for each (RayPtr pRay in rayVector)
 			vpRays.push_back(pRay);
 }
+#pragma endregion
+
+
+
+// --- PUBLIC ---
+// --------------
+#pragma region Public
+CRayGrid CRayGrid::GetNeighboringRays(const std::shared_ptr<CRay>& pRay) const
+{
+	if (!Contains(pRay))
+		return CRayGrid();
+
+	std::vector<int> thetaIdxVec;
+	std::vector<int> phiIdxVec;
+	GetNeighboringAngleIndices(pRay, thetaIdxVec, phiIdxVec);
+
+	std::vector<double> newThetaDeg;
+	std::vector<double> newPhiDeg;
+	RayMatrix newRayMatrix;
+	for each (int idxTheta in thetaIdxVec)
+	{
+		newThetaDeg.push_back(vdThetaDeg[idxTheta]);
+		newRayMatrix.push_back(RayVector());
+		for each (int idxPhi in phiIdxVec)
+		{
+			newPhiDeg.push_back(vdPhiDeg[idxPhi]);
+			newRayMatrix[idxTheta].push_back(vvpRayMatrix[idxTheta][idxPhi]);
+		}
+	}
 
+	const bool newGridHasCircularPhi = bCircularPhi && HasPoleDirection(pRay);
+	return CRayGrid(newRayMatrix, newThetaDeg, newPhiDeg, newGridHasCircularPhi);
+}
+
+void CRayGrid::SetToNeighboringRays(const std::shared_ptr<CRay>& pRay)
+{
+	*this = GetNeighboringRays(pRay);
+}
+#pragma endregion
+
+
+// --- PROTECTED HELPERS ---
+// -------------------------
+#pragma region INDICING
 int CRayGrid::IndexToThetaIndex(const int& idx) const
 {
 	if (idx < 0 || idx >= vpRays.size())
@@ -69,40 +145,43 @@ int CRayGrid::GetIndex(const std::shared_ptr<CRay>& pRay) const
 	return std::distance(vpRays.cbegin(), it);
 }
 
-
-
-// --- PUBLIC ---
-// --------------
-
-CRayGrid CRayGrid::GetNeighboringRays(std::shared_ptr<CRay> pRay) const
+void CRayGrid::GetNeighboringAngleIndices(const std::shared_ptr<CRay>& pRay, std::vector<int>& thetaIdxVec, std::vector<int>& phiIdxVec) const
 {
-	if (!Contains(pRay))
-		return CRayGrid();
-
-	const int idx = GetIndex(pRay);
-	const int idxTheta = IndexToThetaIndex(idx);
-	const int idxPhi = IndexToPhiIndex(idx);
-
-	const int idxThetaMin = std::max(0, idxTheta - 1);
-	const int idxThetaMax = std::min(NTheta(), idxTheta + 1);
-	const int idxPhiMin = std::max(0, idxPhi - 1);
-	const int idxPhiMax = std::min(NTheta(), idxPhi + 1);
-
-	std::vector<double> newThetaDeg = std::vector<double>(vdThetaDeg.cbegin() + idxThetaMin, vdThetaDeg.cend() + idxThetaMax);
-	std::vector<double> newPhiDeg = std::vector<double>(vdPhiDeg.cbegin() + idxPhiMin, vdPhiDeg.cend() + idxPhiMax);
-
-	RayMatrix newRayMatrix;
-	for (int idxTheta = idxThetaMin; idxTheta < idxThetaMax; idxTheta++)
+	thetaIdxVec.clear();
+	phiIdxVec.clear();
+
+	const int idxRay = GetIndex(pRay);
+	if (idxRay < 0)
+		return;
+	const int idxTheta = IndexToThetaIndex(idxRay);
+	const int idxPhi = IndexToPhiIndex(idxRay);
+
+	
+	thetaIdxVec = { idxTheta - 1, idxTheta, idxTheta + 1 };
+	if (thetaIdxVec.front() < 0)
+		thetaIdxVec.erase(thetaIdxVec.begin());
+	if (thetaIdxVec.back() >= NTheta())
+		thetaIdxVec.erase(thetaIdxVec.end());
+	
+
+	if (HasPoleDirection(pRay))
+		for (int idx; idx < NPhi(); idx++)
+			phiIdxVec.push_back(idx);
+	else if (bCircularPhi)
+		phiIdxVec = { (idxPhi - 1 + NPhi()) % NPhi(), idxPhi, (idxPhi + 1) % NPhi() }; //Circular indexing;
+	else
 	{
-		newRayMatrix.push_back(RayVector());
-		for (int idxPhi = idxPhiMin; idxPhi < idxPhiMax; idxPhi++)
-			newRayMatrix[idxTheta].push_back( vvpRayMatrix[idxTheta][idxPhi] );
+		phiIdxVec = { idxPhi - 1, idxPhi, idxPhi + 1 };
+		if (phiIdxVec.front() < 0)
+			phiIdxVec.erase(phiIdxVec.begin());
+		if (phiIdxVec.back() >= NPhi())
+			phiIdxVec.erase(phiIdxVec.end());
 	}
-
-	return CRayGrid(newRayMatrix, newThetaDeg, newPhiDeg);
 }
+#pragma endregion
 
-void CRayGrid::SetToNeighboringRays(std::shared_ptr<CRay> pRay)
+bool CRayGrid::HasPoleDirection(const std::shared_ptr<CRay>& pRay) const
 {
-	*this = GetNeighboringRays(pRay);
-}
+	const VistaVector3D n0 = pRay->InitialDirection();
+	return abs(n0[Vista::X]) < DBL_EPSILON && abs(n0[Vista::Y]) < DBL_EPSILON; //x- and y- component are zero
+}
\ No newline at end of file

src/ITAPropagationPathSim/AtmosphericRayTracing/Simulation/RayResolutionAdapter.cpp

@@ -10,22 +10,30 @@ using namespace ITAPropagationPathSim::AtmosphericRayTracing;
 
 bool Simulation::CRayResolutionAdapter::ZoomIntoRay(std::shared_ptr<CRay> pRay)
 {
+	vpInsertedRays.clear();
+
 	if(!Contains(pRay))
 		return false;
+	if (ConstRayVectorReference().size() <= 1)
+		return false;
 
 	SetToNeighboringRays(pRay);
+
 	DoubleRayResolution();
 	return true;
 }
 
 void Simulation::CRayResolutionAdapter::DoubleRayResolution()
 {
-	vdThetaDeg = DoubleAngleResolution(vdThetaDeg);
-	vdPhiDeg = DoubleAngleResolution(vdPhiDeg);
+	if (ConstRayVectorReference().size() <= 1)
+		return;
+
+	vdThetaDeg = DoubleAngularResolution(vdThetaDeg);
+	vdPhiDeg = DoubleAngularResolution(vdPhiDeg);
 
 	RayMatrix newRayMatrix;
 	vpInsertedRays.clear();
-	RayVector::const_iterator iteratorOldRays = vpRays.cbegin();
+	RayVector::const_iterator iteratorOldRays = ConstRayVectorReference().cbegin();
 	for (int idxTheta = 0; idxTheta < vdThetaDeg.size(); idxTheta++)
 	{
 		const bool isNewTheta = (idxTheta % 2) == 1;
@@ -44,13 +52,12 @@ void Simulation::CRayResolutionAdapter::DoubleRayResolution()
 
 			newRayMatrix[idxTheta].push_back(pRay);
 		}
-
-		vvpRayMatrix = newRayMatrix;
-		InitRayVectorFromRayMatrix();
 	}
+	vvpRayMatrix = newRayMatrix;
+	InitRayVectorFromRayMatrix();
 }
 
-std::vector<double> Simulation::CRayResolutionAdapter::DoubleAngleResolution(const std::vector<double>& angleVector)
+std::vector<double> Simulation::CRayResolutionAdapter::DoubleAngularResolution(const std::vector<double>& angleVector) const
 {
 	if (angleVector.size() < 2)
 		return angleVector;

src/ITAPropagationPathSim/AtmosphericRayTracing/Simulation/RayResolutionAdapter.h

@@ -21,9 +21,6 @@
 
 #include <ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.h>
 
-// ITA includes
-//#include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
-
 
 namespace ITAPropagationPathSim
 {
@@ -43,8 +40,7 @@ namespace ITAPropagationPathSim
 
 			private:
 				void DoubleRayResolution();
-				std::vector<double> DoubleAngleResolution(const std::vector<double>& angleVector);
-
+				std::vector<double> DoubleAngularResolution(const std::vector<double>& angleVector) const;
 			};
 		}
 	}