Added calculation of diffraction aperture points as well as the image source...

Added calculation of diffraction aperture points as well as the image source construction of secondary sources created by the edges.

include/ITAPropagationPathSim/CombinatedModel/DiffractionApertures.h → include/ITAPropagationPathSim/CombinatedModel/DiffractionLocator.h

 #ifndef INCLUDE_WATCHER_DIFFRACTION_APERTURES
 #define INCLUDE_WATCHER_DIFFRACTION_APERTURES
 
+
+// ITA includes
+#include <ITAGeo/Base.h>
+#include <ITAGeo/Halfedge/MeshModel.h>
+#include <ITAGeo/Utils.h>
+#include<ITAPropagationPathSim/CombinatedModel/PropagationShapes.h>
+
+
 namespace ITAPropagationPathSim
 {
 	namespace CombinatedModel
@@ -8,9 +16,9 @@ namespace ITAPropagationPathSim
 		using namespace ITAGeo;
 		using namespace std;
 
-		namespace  ImageConstruction
+		namespace  Diffraction
 		{
-
+			bool ConstructAperturePoints(shared_ptr<const CEmitter> pEmitter, shared_ptr<const CSensor> pSensor, vector<CPropagationShapeShared> pPropagationLists);
 		}
 
 	}

include/ITAPropagationPathSim/CombinatedModel/ImageConstructor.h

@@ -22,6 +22,10 @@ namespace ITAPropagationPathSim
 			ITA_PROPAGATION_PATH_SIM_API bool ConstructImageSources(shared_ptr<const ITAGeo::CEmitter> pEmitter, vector<CPropagationShapeShared>& vpPropagationTree);
 			ITA_PROPAGATION_PATH_SIM_API bool ConstructImageReceivers(shared_ptr<const CSensor> pSensor, vector<CPropagationShapeShared>& vpPropagationLists, CShapesMap vpVisibilityMap, const int iMaxReflectionOrder);
 
+			ITA_PROPAGATION_PATH_SIM_API bool ConstructImageApertures(vector<CPropagationShapeShared>& vpPropagationLists);
+
+
+
 			ITA_PROPAGATION_PATH_SIM_API bool ConstructImageEdges(shared_ptr<const ITAGeo::CEmitter> pEmitter, vector<CPropagationShapeShared> vpPropagationTree);
 		}
 	}

include/ITAPropagationPathSim/CombinatedModel/PropagationEngine.h

@@ -25,6 +25,7 @@ namespace ITAPropagationPathSim
 
 			void ApplySensor(shared_ptr<const CSensor> pSensor);
 
+			void ConstructPropagationPaths(CPropagationPathList& oPaths);
 
 
 		private:

include/ITAPropagationPathSim/CombinatedModel/_SourceFiles.cmake

@@ -5,12 +5,12 @@ set( RelativeDir "include/ITAPropagationPathSim/CombinatedModel" )
 set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim\\CombinatedModel" )
 
 set( DirFiles 
-	DiffractionApertures.h
+	DiffractionLocator.h
 	ImageConstructor.h
 	PropagationEngine.h
 	PropagationConstructor.h
 	PropagationShapes.h
-	ReflectionApertures.h
+	ReflectionLocator.h
 	)
 		
 set( DirFiles_SourceGroup "${RelativeSourceGroup}" )

src/ITAPropagationPathSim/CombinatedModel/DiffractionApertures.cpp

-#include <ITAPropagationPathSim\CombinatedModel\PropagationEngine.h>

src/ITAPropagationPathSim/CombinatedModel/DiffractionLocator.cpp

+#include <ITAPropagationPathSim\CombinatedModel\DiffractionApertures.h>
+
+using namespace std;
+using namespace ITAGeo;
+
+bool ITAPropagationPathSim::CombinatedModel::Diffraction::ConstructAperturePoints(shared_ptr<const CEmitter> pEmitter, shared_ptr<const CSensor> pSensor, vector<CPropagationShapeShared> pPropagationLists)
+{
+	for (auto& pPropagationStartPoint : pPropagationLists)
+	{
+		if (pPropagationStartPoint == nullptr)
+			continue;
+
+		//E_in and E_out are direction vectors of the edges. 
+		//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, vpVertexDiff;
+		
+		//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 location of the emitter 
+		shared_ptr<VistaVector3D> v3LastFromVertex = 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
+		VistaVector3D v3SensorPosition = pSensor->v3InteractionPoint;
+
+		
+		CPropagationShapeShared pShape = pPropagationStartPoint;
+		CPropagationEdgeShared pEdge;
+		while (pShape != nullptr)
+		{
+			//Diffraction can only occur on edges
+			if (pShape->iShapeType != CPropagationShape::EDGE)
+			{
+				//If the shape is a face and is the last shape of the current list, break the while loop
+				//The face is needed for the calculation of the image receiver
+				if (pShape->pChild == nullptr)
+					break;
+
+				pShape = pShape->pChild;
+			}
+			else //pShape is an edge
+			{
+				pEdge = static_pointer_cast<CPropagationEdge>(pShape);
+
+				//The start value for the aperture position is on the mid of the edge
+				vfRelativeApexPosition.push_back(0.5);
+
+				//The edge direction for the outgoing part is the difference of the end points of the edge
+				vpE_out.push_back(make_shared<VistaVector3D>(*pEdge->v3ToVertex - *pEdge->v3FromVertex));
+
+				//Get the squared length
+				vfE_square.push_back(vpE_out.back()->GetLengthSquared());
+
+				//The first shape has no parent. If the parent is a face, the edge must be mirrored(already done in a previous step)
+				if (pEdge->pParent != nullptr && pEdge->pParent->iShapeType == CPropagationShape::FACE)
+				{
+					vpE_in.push_back(make_shared<VistaVector3D>(*pEdge->v3ToVertexMirrored - *pEdge->v3FromVertexMirrored));
+
+					vpVertexDiff.push_back(make_shared<VistaVector3D>(*v3LastFromVertex - *pEdge->v3FromVertexMirrored));
+				}
+				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_shared<VistaVector3D>(*v3LastFromVertex - *pEdge->v3FromVertex));
+
+				}
+
+				//Set the current from vertex as the next "last" one
+				v3LastFromVertex = pEdge->v3FromVertex;
+
+				//Set the next shape to the child of the current shape
+				pShape = pShape->pChild;
+
+			}
+		}
+
+		//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;
+			}
+			else //If the (imaged) sensor is not visible, the whole path is invalid
+			{
+				pPropagationStartPoint = nullptr;
+				break;
+			}
+		}
+
+		//In the last step, the path candidate could be declared as invalid. If it is invalid, go to the next path candidate
+		if (pPropagationStartPoint == nullptr)
+			continue;
+
+		//Add the (maybe imaged) receiver to the differences vector
+		vpVertexDiff.push_back(make_shared<VistaVector3D>(*v3LastFromVertex - v3SensorPosition));
+
+		//Iterative a given number of times to get the exact position of the diffraction apertures
+		const size_t iNumIterations = 5;
+
+		for(size_t iCurrentIteration = 0; iCurrentIteration< iNumIterations; iCurrentIteration++)
+		{
+			//Length of direction vector between two adjacent points in the shortest path.
+			vector<float> vfDirectionVectorLength;
+
+			//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());
+
+			//Get the length of the shortest paths between two adjecent aperture points 
+			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());
+			}
+
+			//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());
+
+			// The terms can be written in short form to:
+			//  b[i-1] *  vfRelativeApexPosition[i-1] + a[i] *  vfRelativeApexPosition[i] + b[i] *  vfRelativeApexPosition[i+1] = c[i]
+			// with b[-1] *  vfRelativeApexPosition[-1] = 0 and b[n-1] *  vfRelativeApexPosition[n] =0 and with n, the number of diffractions
+			//┌                                                               ┐ ┌                            ┐ ┌       ┐
+			//| a[0]   b[0]   0      0      0  ... 0      0      0      0     | | vfRelativeApexPosition[0]  | | c[0]  |
+			//| b[0]   a[1]   b[1]   0      0  ... 0      0      0      0     | | vfRelativeApexPosition[1]  | | c[1]  |
+			//| 0      b[1]   a[2]   b[2]   0  ... 0      0      0      0     | | vfRelativeApexPosition[2]  | | c[2]  |
+			//|	...    ...    ...    ...    ...    ...    ...    ...    ...   |*| ...                        |=| ...   |
+			//| 0      0      0      0      0  ... 0      b[n-3] a[n-2] b[n-2]| | vfRelativeApexPosition[n-2]| | c[n-2]|
+			//| 0      0      0      0      0  ... 0      0      b[n-2] a[n-1]| | vfRelativeApexPosition[n-1]| | c[n-1]|
+			//└                                                               ┘ └                            ┘ └       ┘
+			vector<float> a, b, c;
+			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]);
+			}
+			for (int i = 0; i < vfE_square.size() - 1; i++)
+			{
+				b.push_back(vpE_out[i]->Dot(*vpE_in[i + 1] / vfDirectionVectorLength[i + 1]));
+			}
+
+			//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]  |
+			//|	...    ...    ...    ...    ...    ...    ...   |*| ...                        |=| ...   |
+			//| 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]|
+			//└                                                 ┘ └                            ┘ └       ┘
+			vector<float> d, e;
+			
+			//The first row remains the same
+			d.push_back(a.front());
+			e.push_back(c.front());
+
+			//
+			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]);
+			}
+
+			//Solve the equation system beginning at the last row
+			vfRelativeApexPosition.back() = e.back() / d.back();
+			for (int i = vfRelativeApexPosition.size() - 2; i >= 0; i--)
+			{
+				vfRelativeApexPosition[i] = (e[i] - b[i] * vfRelativeApexPosition[i + 1]) / d[i];
+			}
+		}
+
+		//After the last iteration check vfRelativeApexPosition for valid values [0;1]
+		for (auto& fRelApexPos : vfRelativeApexPosition)
+		{
+			//If one apex point does not lie on the edge, the whole path is invalid
+			if (fRelApexPos > 1 || fRelApexPos < 0)
+			{
+				pPropagationStartPoint = nullptr;
+				break;
+			}
+		}
+
+		//Set the aperture points of the propagation candidate starting with the last edge
+		while (pEdge != nullptr)
+		{
+			//Set the interaction point to the aperture point location
+			//v3InteractionPoint = from_vertex + vfRelativeApexPosition * (to_vertex - from_vertex)
+			pEdge->v3InteractionPoint = make_unique <VistaVector3D>(*pEdge->v3FromVertex + vfRelativeApexPosition.back()*(*pEdge->v3ToVertex - *pEdge->v3FromVertex));
+
+			//Pop the last element
+			vfRelativeApexPosition.pop_back();
+
+			//Get the previous shape
+			pShape = pEdge->pParent;
+
+			//Loop until the next edge is found or no further shape are before the current one
+			while (pShape != nullptr && pShape->iShapeType != CPropagationShape::EDGE)
+			{
+				pShape == pEdge->pParent;
+			}
+
+			//Break if no further edge are found
+			if (pShape == nullptr)
+				break;
+
+			//If an edge is found set pEdge
+			pEdge = static_pointer_cast<CPropagationEdge>(pShape);
+
+		}
+	}
+
+}

src/ITAPropagationPathSim/CombinatedModel/ImageConstructor.cpp

@@ -197,8 +197,70 @@ namespace ITAPropagationPathSim
 
 				return true;
 			}
+
+			bool ConstructImageApertures(vector<CPropagationShapeShared>& vpPropagationLists)
+			{
+				for (auto& pStartShape : vpPropagationLists)
+				{
+					CPropagationShapeShared pShape = pStartShape;//Start point
+					while (pShape != nullptr)
+					{
+						while (pShape != nullptr && pShape->iShapeType != CPropagationShape::EDGE)
+						{
+							pShape = pShape->pChild;
+						}
+
+						//No edge found, go to next start shape
+						if (pShape == nullptr)
+							continue;
+
+						//Set secondary source to the aperture point of the edge
+						VistaVector3D v3SecondarySource = *pShape->v3InteractionPoint;
+
+						//Set pShape to the child of the edge
+						pShape = pShape->pChild;
+
+						//Vector of faces between two edges
+						vector<CPropagationFaceShared> vpFaces;
+						while (pShape != nullptr && pShape->iShapeType == CPropagationShape::FACE)
+						{
+							vpFaces.push_back(static_pointer_cast<CPropagationFace>(pShape));
+							pShape = pShape->pChild;
+						}
+
+						//Construct image sources between two edges
+						if (!ConstructImageSources(v3SecondarySource, vpFaces))
+						{
+							//If no valid image source could be created, ignore the propagation candidate and set its pointer to zero
+							pStartShape = nullptr;
+							break; //break the while loop
+						}
+					}
+				}
+			}
 		
 
+			bool ConstructImageSources(const VistaVector3D& v3SourceIn, vector<CPropagationFaceShared>& pFacesOut)
+			{
+				VistaVector3D v3ImageSource = v3SourceIn;
+				for (auto& pFace : pFacesOut)
+				{
+					if (ITAGeoUtils::IsPointInFrontOfPlane(*pFace->pPlane, v3ImageSource))
+					{
+						ITAGeoUtils::MirrorPointOverPlane(v3ImageSource, *pFace->pPlane, v3ImageSource);
+						pFace->v3ImageSource = make_shared <VistaVector3D>(v3ImageSource);
+
+					}
+					else //Image source is not illuminable and thus, the whole path is invalid
+					{
+						return false;
+					}
+				}
+
+				//No invalid image source found
+				return true;
+			}
+
 
 			bool RecursiveConstructImageEdge(shared_ptr<const CPropagationFace> pParentFace, shared_ptr<CPropagationEdge> pPropagationEdge)
 			{

src/ITAPropagationPathSim/CombinatedModel/PropagationEngine.cpp

 #include <ITAPropagationPathSim\CombinatedModel\PropagationEngine.h>
 #include <ITAPropagationPathSim\CombinatedModel\ImageConstructor.h>
+#include <ITAPropagationPathSim\CombinatedModel\DiffractionApertures.h>
 
 //Typedefs
 typedef OpenMesh::PolyMesh_ArrayKernelT<> CITAMesh;
@@ -63,6 +64,20 @@ void ITAPropagationPathSim::CombinatedModel::CPathEngine::ApplySensor(shared_ptr
 	ImageConstruction::ConstructImageReceivers(m_pSensor, m_vpPropagationLists, m_mvpShapeVisibilityMap, *m_pMaxReflectionOrder);
 }
 
+void ITAPropagationPathSim::CombinatedModel::CPathEngine::ConstructPropagationPaths(CPropagationPathList& oPaths)
+{
+	//First, construct the aperture points
+	Diffraction::ConstructAperturePoints(m_pEmitter, m_pSensor, m_vpPropagationLists);
+
+	//Thereafter, with the constructed aperture points, construct the image apertures as secondary image sources 
+	//for the faces after each edge
+	ImageConstruction::ConstructImageApertures(m_vpPropagationLists);
+
+	//Construct the intersection points of the reflections
+
+}
+
+
 void ITAPropagationPathSim::CombinatedModel::CPathEngine::CreatePropagationTree()
 {
 	//Set initial values for number of diffractions/reflections

src/ITAPropagationPathSim/CombinatedModel/_SourceFiles.cmake

@@ -6,11 +6,11 @@ set( RelativeSourceGroup "Source Files\\ITAPropagationPathSim\\CombinatedModel"
 
 set( DirFiles 
 	DiffractionApertures.cpp
-	ImageConstructor.cpp
+	ImageLocator.cpp
 	PropagationEngine.cpp
 	PropagationConstructor.cpp
 	PropagationShapes.cpp
-	ReflectionApertures.cpp
+	ReflectionLocator.cpp
 	)
 		
 set( DirFiles_SourceGroup "${RelativeSourceGroup}" )