Fixing merge problem

include/ITAPropagationPathSim/CombinedModel/PropagationEngine.h

@@ -24,24 +24,25 @@ namespace ITAPropagationPathSim
 			CPathEngine();
 
 			void Configure(const bool bFilterNotNeighbouredEdges = false, const bool bFilterIlluminatedRegionDiffraction = false,
- 
-			                     const bool bFilterEdgeToEdgeIntersectedPaths = false,
- const bool bFilterEmitterToEdgeIntersectedPaths = false, 
-			                     const bool bFilterSensorToEdgeIntersectedPaths = false,
- const bool bFilterIntersectedPaths = false,
-			                     const int iNumberIterationApexCalculation = 5,
- const float fIntersectionTestResolution = 0.001f,
-			                     const float fAccumulatedAngleThreshold = -1.0f);
+				const bool bFilterEdgeToEdgeIntersectedPaths = false, const bool bFilterEmitterToEdgeIntersectedPaths = false, 
+				const bool bFilterSensorToEdgeIntersectedPaths = false, const bool bFilterIntersectedPaths = false,
+				const int iNumberIterationApexCalculation = 5, const float fIntersectionTestResolution = 0.001f);
+
+			void SetAbortionCriteria(const int iMaxDiffractions, const int iMaxReflections, const int iMaxCombinedOrder,
+				const float fMaxLevelReduction=-1.f, const float fReflectionPenalty=0.f, const float fDiffractionPenalty=0.f,
+				const float fAccumulatedAngleThreshold = -1.0f);
 
 			void InitializePathEnvironment(shared_ptr<const Halfedge::CMeshModelList> pMeshModelList);
 
-			void ApplyEmitter(shared_ptr<CEmitter> pEmitter, const int iMaxDiffractions, const int iMaxReflections, const int iMaxCombinedOrder);
+			void ApplyEmitter(shared_ptr<CEmitter> pEmitter);
 
 			void ApplySensor(shared_ptr<CSensor> pSensor);
 
-			void ConstructPropagationPaths(CPropagationPathList& oPaths);
 
+			void ConstructPropagationPathsWithStopWatch(CPropagationPathList& oPaths);
+			void ConstructPropagationPaths(CPropagationPathList& oPaths);
 
+			size_t GetNumberPropagationPathCandidates();
 
 		private:
 
@@ -49,13 +50,20 @@ namespace ITAPropagationPathSim
 			shared_ptr<const ITAGeo::Halfedge::CMeshModelList> m_pMeshModelList; //!< Copy of propagation model list
 			CShapesMap m_mvpShapeVisibilityMap; //!< Map of which shape can illuminate which ones
 			vector<CPropagationShapeShared> m_vpPropagationTree; //!< Vector of all propagation shapes (as well as its children) visible to the emitter
-			vector<CPropagationShapeShared> m_vpPropagationLists; //!< Vector of all propagation shape visible to the emitter with each one pointing to up to one child with the last one visible to the sensor
+			vector<CPropagationShapeShared> m_vpPropagationCandidates; //!< Vector of all propagation shape visible to the emitter with each one pointing to up to one child with the last one visible to the sensor
 			shared_ptr<CEmitter> m_pEmitter; //!< Const pointer to emitter
 			shared_ptr<CSensor> m_pSensor; //!< Const pointer to sensor
+
+			//Abortion criteria
 			unique_ptr<const int> m_pMaxReflectionOrder; //!< Const pointer to maximum reflection order
 			unique_ptr<const int> m_pMaxDiffractionOrder; //!< Const pointer to maximum diffraction order
 			unique_ptr<const int> m_pMaxCombinedOrder; //!< Const pointer to maximum order of combined diffractions and reflections
-			
+			unique_ptr<const float> m_pMaxLevelReduction; //!< Const pointer to maximum level reduction
+			unique_ptr<const float> m_pMaxPropagationRange; //!< Const pointer to maximum propagation range [m]
+			unique_ptr<const float> m_pReflectionPenalty; //!< Const pointer to reflection penalty 
+			unique_ptr<const float> m_pDiffractionPenalty; //!< Const pointer to angle independent diffraction penalty
+			unique_ptr <const float> m_pAccumulatedAngleThreshold; //!< Threshold for the accumulated diffraction angle
+
 			//Filter member variables
 			unique_ptr<const bool> m_pFilterNotNeighbouredEdges; //!< Filter not neighboured edges for faster filter calculation 
 			unique_ptr<const bool> m_pFilterIlluminatedRegionDiffraction; //!< Boolean pointer whether to ignore diffractions that can only diffract in illuminated regions(and have therefore a low impact to the whole IR)
@@ -64,7 +72,6 @@ namespace ITAPropagationPathSim
 			unique_ptr<const bool> m_pFilterSensorToEdgeIntersectedPaths; //!< Boolean pointer for filtering paths, where always an intersection between the sensor and an edge occurs
 			unique_ptr<const bool> m_pFilterIntersectedPaths; //!< Boolean pointer for for filtering paths if intersection occurs
 			unique_ptr<const float> m_pIntersectionTestResolution;  //!< Resolution for intersection test
-			unique_ptr <const float> m_pAccumulatedAngleThreshold; //!< Threshold for the accumulated diffraction angle
 			unique_ptr <const size_t> m_pNumberIterationApexCalc; //!< Number of iterations for the calculation of the aperture points
 
 			
@@ -80,7 +87,7 @@ namespace ITAPropagationPathSim
 
 			//Create the propagation tree
 			void CreatePropagationTree();
-			void RecursiveAddChildrenToTree(const CPropagationShapeShared & pPropagationShapeChildIn, CPropagationShapeShared& vpShapeChildCopyOut, int iReflectionOrder, int iDiffractionOrder, int iCombinedOrder);
+			void RecursiveAddChildrenToTree(const CPropagationShapeShared & pPropagationShapeChildIn, CPropagationShapeShared& vpShapeChildCopyOut, int iReflectionOrder, int iDiffractionOrder, int iCombinedOrder, float fAccumulatedPenalty, float fFirstDistance, float fSecondDistance, float fLevelDropAtFirstEdge);
 			
 			//Create the propagation list
 			void CreatePropagationLists();
@@ -89,6 +96,9 @@ namespace ITAPropagationPathSim
 			//Filter not visible, intersected, paths
 			void FilterVisiblePaths(const vector<CPropagationShapeShared> vpAllPathsIn, vector<CPropagationShapeShared>& vpVisiblePathsOut);
 
+			//Culling of unperceived paths
+			void PerceptionalCulling(const vector<CPropagationShapeShared> vpAllPathsIn, vector<CPropagationShapeShared>& vpAudiblePathsOut);
+
 
 			//Convert the propagation list
 			void ConvertShapeListsToPropagationPaths(ITAGeo::CPropagationPathList& oPathsOut);

include/ITAPropagationPathSim/CombinedModel/PropagationShapes.h

@@ -25,11 +25,12 @@ namespace ITAPropagationPathSim
 		typedef shared_ptr<CPropagationFace> CPropagationFaceShared;
 		typedef shared_ptr<CPropagationEdge> CPropagationEdgeShared;
 		typedef OpenMesh::PolyMesh_ArrayKernelT<> CITAMesh;
-		typedef map <CPropagationShapeShared, vector<CPropagationShapeShared>> CShapesMap;
+		typedef map <size_t, vector<CPropagationShapeShared>> CShapesMap;
 
 		typedef string MeshModelHandle;
 		typedef pair < MeshModelHandle, CITAMesh::FaceHandle> PropagationFaceHandle;
 
+		//Old implementation[TODO: remove]
 		//struct ImageReceiverMap : map<PropagationFaceHandle, shared_ptr<ImageReceiverMap>>
 	/*	{
 			shared_ptr<VistaVector3D> v3ImageReceiver;
@@ -46,7 +47,7 @@ namespace ITAPropagationPathSim
 		class CPropagationShape
 		{
 		public:
-			//Public member variables
+			//---Public member variables-----------------------------------------------------------------------
 			
 			//!< Type of shape(e.g. edge or face)
 			enum EShapeType
@@ -65,6 +66,12 @@ namespace ITAPropagationPathSim
 			//!< Point of interaction
 			unique_ptr<VistaVector3D> v3InteractionPoint;
 
+			//!< Shape handle
+			OpenMesh::BaseHandle hShape;
+
+			//!< Vertices of shape
+			vector<shared_ptr<VistaVector3D>> vv3Vertices;
+
 			//!< Parent shape(previous shape)
 			weak_ptr<CPropagationShape> pParent;
 
@@ -74,24 +81,52 @@ namespace ITAPropagationPathSim
 
 			shared_ptr<string> sMeshModelName;
 
-			//Public member functions
+			//!< Barycenter of shape (used for boundary sphere)
+			shared_ptr<VistaVector3D> v3Barycenter;
+
+			//!< Radius of boundary sphere (longest distance from barycenter to vertex)
+			float fRadius;
+
+			//!< Minimum distance to previous shape
+			float fMinimumDistance = 0.0f;
+
+			//!< Distance level drop
+			float fDistanceLevelDrop = 0.0f;
+			//---Public member functions--------------------------------------------------------------------
 
+			//!<Copy shape
 			void CopyFrom(const CPropagationShape& oPropagationShapeIn);
 
-		private:
+			void GetBoundingBoxAxesAligned(VistaVector3D& v3Min, VistaVector3D& v3Max) const;
+
+			//!< Set barycenter and radius of the boundary sphere of the shape
+			bool SetBoundarySphere();
+
+			bool SetMinimumDistance(const CPropagationShape& oPropagationShape);
+
+			//!< Set fMinimumDistance to point (e.g. emitter)
+			bool SetMinimumDistance(const VistaVector3D& v3Point);
+			
+			//!< Calculate minimum distance to point (e.g. emitter) and return value
+			float CalculateMinimumDistance(const VistaVector3D& v3Point);
+
+			//!< Calculate the level drop due to the propagation of the sound
+			float CalculateLevelDrop(float& fFirstDistance, float& fSecondDistance, float& fLevelDropAtFirstEdge);
+
+			size_t GetIdentifier() const;
+
+
 		};
 
 		class CPropagationFace : public CPropagationShape
 		{
 		public:
-			//Public member variables
+			//---Public member variables--------------------------------------------------------------------
 			
 			CITAMesh::FaceHandle hFace; //Face handle corresponding to propagation face
 
 			shared_ptr<VistaPlane> pPlane; //Plane built up by face
 
-			vector<shared_ptr<VistaVector3D>> vv3Vertices;
-
 			shared_ptr<VistaVector3D> v3ImageSource;
 
 			shared_ptr<VistaVector3D> v3ImageEdgeSourceStart;
@@ -99,15 +134,26 @@ namespace ITAPropagationPathSim
 
 			shared_ptr<bool> pHasValidImageSource;
 
-			//Public member functions
+			//---Public member functions----------------------------------------------------------------------
+
+			//!< Copy face
 			void CopyFrom(const CPropagationFace& oPropagationFaceIn);
 
+			// Set minimum distance to shape
+			bool SetMinimumDistance(const CPropagationShape& oPropagationShape);
+
+			// Set fMinimumDistance to point (e.g. emitter)
+			bool SetMinimumDistance(const VistaVector3D& v3Point);
+
+			// Calculate minimum distance to point (e.g. emitter) and return value
+			float CalculateMinimumDistance(const VistaVector3D& v3Point);
+
 		};
 
 		class CPropagationEdge : public CPropagationShape
 		{
 		public:
-			//Public member variables
+			//---Public member variables-------------------------------------------------------------------------
 			CITAMesh::HalfedgeHandle hHalfedge;
 			CITAMesh::EdgeHandle hEdge;
 			CITAMesh::FaceHandle hMainFace;
@@ -129,9 +175,22 @@ namespace ITAPropagationPathSim
 			//!< If the relative position lies between 0.0(at v3FromVertex) and 1.0(at v3ToVertex), the interaction point is on the edge 
 			unique_ptr<float> pRelativeInteractionPoint;
 
-			//Public member functions
+			//---Public member functions-----------------------------------------------------------------------
+
+			//!< Copy edge
 			void CopyFrom(const CPropagationEdge& oPropagationEdgeIn);
 
+			// Set minimum distance to shape
+			bool SetMinimumDistance(const CPropagationShape& oPropagationShape);
+
+			// Set fMinimumDistance to point (e.g. emitter)
+			bool SetMinimumDistance(const VistaVector3D& v3Point);
+
+			// Calculate minimum distance to point (e.g. emitter) and return value
+			float CalculateMinimumDistance(const VistaVector3D& v3Point);
+
+			//!< Calculate the level drop due to the propagation of the sound
+			float CalculateLevelDrop(float& fFirstDistance, float& fSecondDistance, float& fLevelDropAtFirstEdge);
 		};
 	}
 }

include/ITAPropagationPathSim/CombinedModel/RTree.h

+#ifndef INCLUDE_WATCHER_PROPAGATION_RTREE
+#define INCLUDE_WATCHER_PROPAGATION_RTREE
+
+
+// ITAGeo
+#include <ITAGeo/Definitions.h>
+#include <ITAGeo/Utils.h>
+#include <ITAGeo/Base.h>
+#include <ITAGeo/ModelBase.h>
+
+#include <ITAGeo/Halfedge/MeshModel.h>
+
+#include<ITAPropagationPathSim/CombinedModel/PropagationShapes.h>
+
+using namespace std;
+
+namespace ITAPropagationPathSim
+{
+	namespace CombinedModel
+	{
+			// Typedefs
+			typedef OpenMesh::PolyMesh_ArrayKernelT<> CITAMesh;
+
+			//! Buildings RTree class
+			/**
+			* @todo
+			*
+			*/
+			class ITA_GEO_API CShapesTree
+			{
+			public:
+				//---Constructor and Destructor---------------------------------------------
+				CShapesTree();
+				
+				~CShapesTree();
+
+				//---Public member functions------------------------------------------------
+
+				void Create(const vector<CPropagationShapeShared>& vpShapesIn);
+
+				//TODO
+				bool IsLineSegmentIntersected(const VistaVector3D& oLineStart, const VistaVector3D& oLineEnd) const;
+
+
+				//TODO
+				//void GetNonIntersectedFaces(vector<CPropagationShapeShared>& vpFacesOut, const VistaVector3D& oLineStart, const VistaVector3D& oLineStart);
+
+				/*TODO: For faster illumination tests
+				void GetIlluminableShapes(vector<CPropagationShapeShared>& vpShapesOut, const VistaPlane& oPlaneIn);
+				
+				void GetIlluminableShapes(vector<CPropagationShapeShared>& vpShapesOut, const CPropagationFaceShared& pFace);
+
+				void GetIlluminableShapes(vector<CPropagationShapeShared>& vpShapesOut, const CPropagationEdgeShared& pEdge);
+				*/
+
+
+				//---Public member variables------------------------------------------------
+
+			private:
+
+				//---Private member functions-----------------------------------------------
+
+				void SetVertices();
+
+				void CreateBranches();
+
+				//---Private member variables-----------------------------------------------
+
+				//Minimum and maximum of bounding box as well as vertices of box
+				VistaVector3D m_v3TreeMin;
+				VistaVector3D m_v3TreeMax;
+
+				//Vector of shapes within bounding box
+				vector<CPropagationShapeShared> m_vpShapes;
+
+				//Branches
+				vector<shared_ptr<CShapesTree>> m_vpBranches;
+
+
+			};
+	}
+}
+
+
+#endif // INCLUDE_WATCHER_PROPAGATION_RTREE
\ No newline at end of file

include/ITAPropagationPathSim/CombinedModel/_SourceFiles.cmake

@@ -9,6 +9,7 @@ set( DirFiles
 	ImageConstructor.h
 	PropagationEngine.h
 	PropagationShapes.h
+	RTree.h
 	ReflectionLocator.h
 	)
 		

include/ITAPropagationPathSim/UrbanEngine/UrbanDiffraction.h

@@ -178,7 +178,7 @@ namespace ITAPropagationPathSim
 			map <UrbanEdgeHandle, map <UrbanEdgeHandle, double>> m_mdApertureStartDifferenceMap;
 
 			//RTree of building meshes
-			shared_ptr<ITAGeo::Urban::RTree::CBuildings> m_pBuildingsTree;
+			shared_ptr<ITAGeo::Urban::RTree::CShapesTree> m_pBuildingsTree;
 
 			//Matrix of illumination of diffraction edges
 			map<UrbanEdgeHandle, set<shared_ptr<DiffractionEdge>>> m_mDiffractionMatrix;

src/ITAPropagationPathSim/CombinedModel/DiffractionLocator.cpp

@@ -28,14 +28,14 @@ namespace ITAPropagationPathSim
 					//Because of possibly occuring reflections during the paths, E_in can contain mirrored versions of the direction edge directions.
 					//VertexDiff contains the vista vectors of the difference between the FromVertices of two edges
 					vector<shared_ptr<VistaVector3D>> vpE_in, vpE_out;
-					vector<unique_ptr<VistaVector3D>> vpVertexDiff;
+					vector<unique_ptr<VistaVector3D>> vpStartVertexDiff;
 
 					//RelativeApexPosition contains the aperture position relative to the edge from- and to-vertices. 
 					//vfE_square contains the squared length of the edge
 					vector<float> vfRelativeApexPosition, vfE_square;
 
 					//The propagation paths always starts at the position of the emitter 
-					shared_ptr<VistaVector3D> v3LastFromVertex = make_shared<VistaVector3D>(pEmitter->v3InteractionPoint);
+					shared_ptr<VistaVector3D> v3LastStartVertex = make_shared<VistaVector3D>(pEmitter->v3InteractionPoint);
 
 					//The propagation paths always ends at the location of the sensor
 					//Because of possible faces after the last edge, a local copy of the location is needed and will eventually mirrored in a later step
@@ -58,6 +58,18 @@ namespace ITAPropagationPathSim
 						}
 						else //pShape is an edge
 						{
+							//If pEdge is the first edge and the parent face contains no valid image source, the whole path is not able to be valid
+							if (pEdge == nullptr)
+							{
+								if (!pShape->pParent.expired() && pShape->pParent.lock()->iShapeType == CPropagationShape::FACE)
+								{
+									auto pFace = static_pointer_cast<CPropagationFace>(pShape->pParent.lock());
+									if (pFace->pHasValidImageSource == false)
+										continue;
+								}
+							}
+
+
 							pEdge = static_pointer_cast<CPropagationEdge>(pShape);
 
 							//The start value for the aperture position is located on the mid of the edge
@@ -74,19 +86,19 @@ namespace ITAPropagationPathSim
 							{
 								vpE_in.push_back(make_shared<VistaVector3D>(*pEdge->v3ImageEdgeReceiverEnd - *pEdge->v3ImageEdgeReceiverStart));
 
-								vpVertexDiff.push_back(make_unique<VistaVector3D>(*v3LastFromVertex - *pEdge->v3ImageEdgeReceiverStart));
+								vpStartVertexDiff.push_back(make_unique<VistaVector3D>(*v3LastStartVertex - *pEdge->v3ImageEdgeReceiverStart));
 							}
 							else
 							{
 								//The input edge direction does not have to be mirrored and is therefore the same direction as the output edge direction 
 								vpE_in.push_back(vpE_out.back());
 
-								vpVertexDiff.push_back(make_unique<VistaVector3D>(*v3LastFromVertex - *pEdge->v3StartVertex));
+								vpStartVertexDiff.push_back(make_unique<VistaVector3D>(*v3LastStartVertex - *pEdge->v3StartVertex));
 
 							}
 
 							//Set the current from vertex as the next "last" one
-							v3LastFromVertex = pEdge->v3StartVertex;
+							v3LastStartVertex = pEdge->v3StartVertex;
 
 							//Set the next shape to the child of the current shape
 							pShape = pShape->pChild;
@@ -94,30 +106,32 @@ namespace ITAPropagationPathSim
 						}
 					}
 
+					//Old implementation(image edge source used instead of image receiver):
 					//If the last shape is not an edge, the image receiver is needed for the calculation of the aperture points
-					while (pShape != nullptr && pShape->iShapeType == CPropagationShape::FACE)
-					{
-						auto pPlane = static_pointer_cast<CPropagationFace>(pShape)->pPlane;
-
-						// Check for a possible visibility
-						if (ITAGeoUtils::IsPointInFrontOfPlane(*pPlane, v3SensorPosition))
-						{
-							//Mirror the (imaged) sensor until an edge is reached
-							ITAGeoUtils::MirrorPointOverPlane(v3SensorPosition, *static_pointer_cast<CPropagationFace>(pShape)->pPlane, v3SensorPosition);
-
-							pShape = pShape->pParent.lock();
-						}
-						else //If the (imaged) sensor is not visible, the whole path is invalid
-						{
-							bValidAperturePoints = false;
-							break;
-						}
-					}
+					//while (pShape != nullptr && pShape->iShapeType == CPropagationShape::FACE)
+					//{
+					//	auto pPlane = static_pointer_cast<CPropagationFace>(pShape)->pPlane;
+
+					//	// Check for a possible visibility
+					//	if (ITAGeoUtils::IsPointInFrontOfPlane(*pPlane, v3SensorPosition))
+					//	{
+					//		//Mirror the (imaged) sensor until an edge is reached
+					//		ITAGeoUtils::MirrorPointOverPlane(v3SensorPosition, *static_pointer_cast<CPropagationFace>(pShape)->pPlane, v3SensorPosition);
+
+					//		pShape = pShape->pParent.lock();
+					//	}
+					//	else //If the (imaged) sensor is not visible, the whole path is invalid
+					//	{
+					//		bValidAperturePoints = false;
+					//		break;
+					//	}
+					//}
+					//
+					////In the last step, the path candidate could be declared as invalid. 
+					////If it is invalid, go to the next path candidate
+					//if (bValidAperturePoints == false)
+					//	continue;
 
-					//In the last step, the path candidate could be declared as invalid. 
-					//If it is invalid, go to the next path candidate
-					if (bValidAperturePoints == false)
-						continue;
 
 					//Add path candidate to list and go to the next path candidate, if no edge was found
 					if (pEdge == nullptr)
@@ -126,8 +140,22 @@ namespace ITAPropagationPathSim
 						continue;
 					}
 
-					//Add the (maybe imaged) receiver to the differences vector
-					vpVertexDiff.push_back(make_unique<VistaVector3D>(*v3LastFromVertex - v3SensorPosition));
+					//If the last shape is not an edge, the image edge source of the last edge is needed for the calculation of the aperture points
+					if (pShape != nullptr && pShape->iShapeType == CPropagationShape::FACE)
+					{
+						auto pFace = static_pointer_cast<CPropagationFace>(pShape);
+
+						//If the image edge source is empty, there is no valid image edge and thus, the path is invalid
+						if (pFace->v3ImageEdgeSourceStart == nullptr)
+							continue;
+
+						v3LastStartVertex = pFace->v3ImageEdgeSourceStart;
+						vpE_out.back() = make_shared<VistaVector3D>(*pFace->v3ImageEdgeSourceEnd - *pFace->v3ImageEdgeSourceStart);
+
+					}
+
+					//Add the receiver to the differences vector
+					vpStartVertexDiff.push_back(make_unique<VistaVector3D>(*v3LastStartVertex - v3SensorPosition));
 
 
 					for (size_t iCurrentIteration = 0; iCurrentIteration < iNumberIterations; iCurrentIteration++)
@@ -139,15 +167,15 @@ namespace ITAPropagationPathSim
 						//The length is determined by the relative aperture position of the last iteration
 
 						//The first direction vector(between the source and the first (imaged) edge) does not contain a correction vector for the source
-						vfDirectionVectorLength.push_back((*vpVertexDiff[0] - vfRelativeApexPosition[0] * *vpE_in[0]).GetLength());
+						vfDirectionVectorLength.push_back((*vpStartVertexDiff[0] - vfRelativeApexPosition[0] * *vpE_in[0]).GetLength());
 
 						for (int i = 1; i < vfRelativeApexPosition.size(); i++)
 						{
-							vfDirectionVectorLength.push_back((*vpVertexDiff[i] - (vfRelativeApexPosition[i] * *vpE_in[i]) + (vfRelativeApexPosition[i - 1] * *vpE_out[i - 1])).GetLength());
+							vfDirectionVectorLength.push_back((*vpStartVertexDiff[i] - (vfRelativeApexPosition[i] * *vpE_in[i]) + (vfRelativeApexPosition[i - 1] * *vpE_out[i - 1])).GetLength());
 						}
 
 						//The last direction vector(between the last edge and the (imaged) receiver) does not contain a correction vector for the receiver
-						vfDirectionVectorLength.push_back((*vpVertexDiff.back() - vfRelativeApexPosition.back() * *vpE_out.back()).GetLength());
+						vfDirectionVectorLength.push_back((*vpStartVertexDiff.back() + vfRelativeApexPosition.back() * *vpE_out.back()).GetLength());
 
 						// The terms can be written in tridiagonal matrix form as following:
 						//  b[i-1] *  vfRelativeApexPosition[i-1] + a[i] *  vfRelativeApexPosition[i] + b[i] *  vfRelativeApexPosition[i+1] = c[i]
@@ -164,7 +192,7 @@ namespace ITAPropagationPathSim
 						for (int i = 0; i < vfE_square.size(); i++)
 						{
 							a.push_back(vfE_square[i] * (1 / vfDirectionVectorLength[i] + 1 / vfDirectionVectorLength[i + 1]));
-							c.push_back(vpVertexDiff[i]->Dot(*vpE_in[i]) / vfDirectionVectorLength[i] - vpVertexDiff[i + 1]->Dot(*vpE_out[i]) / vfDirectionVectorLength[i + 1]);
+							c.push_back(vpStartVertexDiff[i]->Dot(*vpE_in[i]) / vfDirectionVectorLength[i] - vpStartVertexDiff[i + 1]->Dot(*vpE_out[i]) / vfDirectionVectorLength[i + 1]);
 						}
 						for (int i = 0; i < vfE_square.size() - 1; i++)
 						{
@@ -173,34 +201,35 @@ namespace ITAPropagationPathSim
 
 						//To get vfRelativeApexPosition, the Gaussian-elimination is used to delete 
 						//the term  b[i-1] * vfRelativeApexPosition[i-1] in each row, so that
-						//in the last row, only one variable is left over.
-						//Thereafter the searched vfRelativeApexPosition can be calculated.
-						//d[i] replaces a[i] and e[i] replaces c[i]
 						//┌                                                 ┐ ┌                            ┐ ┌       ┐
-						//| d[0]   b[0]   0      0  ... 0      0      0     | | vfRelativeApexPosition[0]  | | e[0]  |
-						//| 0      d[1]   b[1]   0  ... 0      0      0     | | vfRelativeApexPosition[1]  | | e[1]  |
+						//| 1      e[0]   0      0  ... 0      0      0     | | vfRelativeApexPosition[0]  | | d[0]  |
+						//| 0      1      e[1]   0  ... 0      0      0     | | vfRelativeApexPosition[1]  | | d[1]  |
 						//|	...    ...    ...    ...    ...    ...    ...   |*| ...                        |=| ...   |
-						//| 0      0      0      0  ... 0      d[n-2] b[n-2]| | vfRelativeApexPosition[n-2]| | e[n-2]|
-						//| 0      0      0      0  ... 0      0      d[n-1]| | vfRelativeApexPosition[n-1]| | e[n-1]|
+						//| 0      0      0      0  ... 0      1      e[n-2]| | vfRelativeApexPosition[n-2]| | d[n-2]|
+						//| 0      0      0      0  ... 0      0      1     | | vfRelativeApexPosition[n-1]| | d[n-1]|
 						//└                                                 ┘ └                            ┘ └       ┘
 						vector<float> d, e;
 
 						//The first row remains the same
-						d.push_back(a.front());
-						e.push_back(c.front());
+						d.push_back(c.front()/a.front());
+						
+						if(a.size()>1)
+							e.push_back(b.front()/a.front());
 
-						//
+						for (int i = 1; i < a.size()-1; i++)
+						{
+							e.push_back(b[i] / (a[i] - b[i - 1] * e[i - 1]));
+						}
 						for (int i = 1; i < a.size(); i++)
 						{
-							d.push_back(a[i] - b[i - 1] * b[i - 1] / d[i - 1]);
-							e.push_back(c[i] - b[i - 1] * c[i - 1] / d[i - 1]);
+							d.push_back((c[i] - b[i - 1] * d[i - 1]) / (a[i] - b[i-1]*e[i - 1]));
 						}
 
 						//Solve the equation system beginning at the last row
-						vfRelativeApexPosition.back() = e.back() / d.back();
+						vfRelativeApexPosition.back() = d.back();
 						for (int i = vfRelativeApexPosition.size() - 2; i >= 0; i--)
 						{
-							vfRelativeApexPosition[i] = (e[i] - b[i] * vfRelativeApexPosition[i + 1]) / d[i];
+							vfRelativeApexPosition[i] = (d[i] - e[i] * vfRelativeApexPosition[i + 1]);
 						}
 					}
 
@@ -223,6 +252,73 @@ namespace ITAPropagationPathSim
 						//v3InteractionPoint = from_vertex + vfRelativeApexPosition * (to_vertex - from_vertex)
 						VistaVector3D v3InteractionPoint = (1 - vfRelativeApexPosition.back()) * (*pEdge->v3StartVertex) + vfRelativeApexPosition.back()*(*pEdge->v3EndVertex);
 
+						//Interaction point must lie in valid range of parent and child of edge
+						if (!pEdge->pParent.expired())
+						{
+							if (pEdge->pParent.lock()->iShapeType == CPropagationShape::FACE)
+							{
+								auto pParentFace =  static_pointer_cast<CPropagationFace>(pEdge->pParent.lock());
+
+								bValidAperturePoints  = ITAGeoUtils::IsPointInFrontOfPlane(*pParentFace->pPlane, v3InteractionPoint);
+
+							}
+							else
+							{
+								auto pParentEdge = static_pointer_cast<CPropagationEdge>(pEdge->pParent.lock());
+
+
+								//Illumination test at from vertex position
+								VistaPlane oMainPlane, oOppositePlane;
+								oMainPlane.SetOrigin(*pParentEdge->v3StartVertex);
+								oOppositePlane.SetOrigin(*pParentEdge->v3StartVertex);
+								oMainPlane.SetNormVector(*pParentEdge->v3MainFaceNormal);
+								oOppositePlane.SetNormVector(*pParentEdge->v3OppositeFaceNormal);
+
+								bool bMainFaceIlluminated = ITAGeoUtils::IsPointInFrontOfPlane(oMainPlane, v3InteractionPoint);
+								bool bOppositeFaceIlluminated = ITAGeoUtils::IsPointInFrontOfPlane(oOppositePlane, v3InteractionPoint);
+
+								//Check whether wedge angle is >pi or <pi; TODO: export in own function
+								bool bIsOuterAngle = pParentEdge->v3MainFaceNormal->Cross(*pParentEdge->v3OppositeFaceNormal).Dot(*pParentEdge->v3EndVertex - *pParentEdge->v3StartVertex) > 0;
+
+								if (bIsOuterAngle)
+									bValidAperturePoints = bMainFaceIlluminated || bOppositeFaceIlluminated;
+								else
+									bValidAperturePoints = bMainFaceIlluminated && bOppositeFaceIlluminated;
+							}
+						}
+						if (pEdge->pChild != nullptr)
+						{