new version

include/ITAPropagationPathSim/RayTracing/RTEngine.h

@@ -73,6 +73,16 @@ namespace ITAPropagationPathSim
 			  */
 			unsigned long int GetNumberOfFaces() const;
 
+			//! Changes direction of incident ray according to the principle of specular reflection
+			/**
+			  * @param[out] rRay: Incident ray which direction is to be changed
+			  * @param[in] pMesh: Mesh of the scene
+			  * @param[in] hFace: Face handle for identifying the corresponding reflection face
+			  * @param[in] v3ReflectionPoint: Point of specular reflection on the face
+			  */
+			//void ReflectRayOnFace(VistaRay& rRay, CITAMesh* pMesh, CITAMesh::FaceHandle hFace, VistaVector3D& v3ReflectionPoint);
+
+
 			//! Audibility test for images up to given maximum order with or without culling
 			/**
 			  * @param[[in] pSource Source anchor (sensor, receiver, listener, etc.)

src/ITAPropagationPathSim/RayTracing/RTEngine.cpp

@@ -47,12 +47,30 @@ unsigned long int CEngine::GetNumberOfFaces() const
 	return ( int ) pMesh->n_faces();
 }
 
-bool RayDetectionSphereIntersectionTest(VistaRay rRay, std::shared_ptr< ITAGeo::CPropagationAnchor > pDestination, float fDetectionSphereRadius)
+bool RayDetectionSphereIntersectionTest(VistaRay rRay, std::shared_ptr< ITAGeo::CPropagationAnchor > pDestination, float fDetectionSphereRadius, VistaVector3D& v3DetectionPoint)
 {
-	VistaVector3D v3DirReflPt2RcvPt = pDestination->v3InteractionPoint - rRay.GetOrigin();
-	float tDeviationAngle = acos(rRay.GetDir().Dot(v3DirReflPt2RcvPt.GetNormalized()));
+	VistaVector3D v3RayOrigin2RcvPos = pDestination->v3InteractionPoint - rRay.GetOrigin();
 	
-	return (abs(v3DirReflPt2RcvPt.GetLength() * sin(tDeviationAngle)) < fDetectionSphereRadius);
+	// Nearest Point of ray to detection sphere origin
+	v3DetectionPoint = rRay.GetDir().Dot(v3RayOrigin2RcvPos) * rRay.GetDir() + rRay.GetOrigin();
+
+	return ((v3DetectionPoint - pDestination->v3InteractionPoint).GetLength() <= fDetectionSphereRadius);
+
+	//float tDeviationAngle = acos(rRay.GetDir().Dot(v3RayOrigin2RcvPos.GetNormalized()));
+	//return (abs(v3RayOrigin2RcvPos.GetLength() * sin(tDeviationAngle)) < fDetectionSphereRadius);
+}
+
+void ReflectRayOnFace(VistaRay& rRay, CITAMesh* pMesh, CITAMesh::FaceHandle hFace, VistaVector3D& v3ReflectionPoint)
+{
+	assert(pMesh->has_face_normals());
+	CITAMesh::Normal normal(pMesh->normal(hFace));
+	const VistaVector3D v3FaceNormal(normal.data());
+
+	VistaVector3D v3OrthogonalComponent = rRay.GetDir().Dot(v3FaceNormal) * v3FaceNormal; // direction component of incident ray orthogonal to reflection face
+
+	//! Change direction of ray orthogonal to reflection wall into the opposite direction
+	rRay.SetDir(rRay.GetDir - 2 * v3OrthogonalComponent); 
+	rRay.SetOrigin(v3ReflectionPoint);
 }
 
 
@@ -63,123 +81,86 @@ void CEngine::Run( std::shared_ptr< ITAGeo::CPropagationAnchor > pSource, std::s
 
 	for (unsigned long int n = 0; n < GetNumberOfRays(); n++)
 	{
+		// Create and initialize ray with random direction
 		VistaRay rRay;
 		rRay.SetOrigin(pSource->v3InteractionPoint);
 
-		// using spherical coordinates for direction of ray
 		float fAzimuthRad = VistaRandomNumberGenerator::GetStandardRNG()->GenerateFloat(0.0f, 2 * ITAConstants::PI_F);
 		float fElevationRad = VistaRandomNumberGenerator::GetStandardRNG()->GenerateFloat(0.0f, ITAConstants::PI_F);
 		VistaVector3D v3RandomDirection(cos(fAzimuthRad), sin(fAzimuthRad), cos(fElevationRad));
 
 		rRay.SetDir(v3RandomDirection);
 
-		const int iNumberOfFreqBands = 10;
-		float tAbsorbtionInOctaveBands[iNumberOfFreqBands] = { 0.03f, 0.03f, 0.05f, 0.05f, 0.08f, 0.09f, 0.10f, 0.12f, 0.21f, 0.24f }; // examplary hard coded absorbtion coefficients in octave bands
+		// Initialize variables
+		float tAbsorbtionOfFace = 0.9f; //@todo: dynamic addaption of absorbtion according to material
+		float tEnergy;
+		float tAbsorbtionLossFactor = 1.0f;
+		float tPropagationLossFactor = 1.0f; // considering spherical waves (1-by-r law)
+		float tEnergyThreshold = powf(10, -6);
 
-		float tEnergy[iNumberOfFreqBands] = { 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f };	// @todo: consider directivity
+		bool bDetectionSphereHitByRay = false;
+		bool bFaceIntersectedByRay = false;
+		bool bEnergyBelowThreshold = false;
 
-		bool bTotalEnergyOfBandsBelowThreshold = false; // only true if energy of every frequency band is below threshold
-		//bool bEnergySingleBandBelowThreshold[iNumberOfFreqBands]; // check of energy threshold for each band
+		VistaVector3D v3DetectionPoint;
+		VistaVector3D v3ReflectionPoint;
 		
+		// Store source position as start point of the propagation path
 		ITAGeo::CPropagationPath oPropagationPath;
+		oPropagationPath.push_back(pSource);
 
-		while( !bTotalEnergyOfBandsBelowThreshold )
+		while( !bEnergyBelowThreshold )
 		{
-			// breaks loop if neither reveicer nor any faces were hit by ray
-			bool bFlagValidConstellation = false;
-
-			// Recever hit by ray? (-> direct incident)
-			bool bReceiverHit = RayDetectionSphereIntersectionTest(rRay, pDestination, m_fDetectionSphereRadius);
-			if (bReceiverHit)
-			{
-				bFlagValidConstellation = true;
-				lPropagationPaths.push_back(oPropagationPath);
-				// @todo: store energy of ray in histogram, terminate ray?
-			}
-
 			// Check ray propagation for each face of the mesh
 			auto it = pMesh->faces_begin();
-			while( it != pMesh->faces_end() )
+			while( it != pMesh->faces_end() && !bFaceIntersectedByRay )
 			{
 				// ->siehe MirrorImage->CEngine::AudibilityTest
-				CITAMesh::FaceHandle hFace( *it++ ); // --> mach ein polygon
-
-				//const ITAGeo::ECulling eCulling;
-				VistaVector3D v3IntersectionPoint;
-				bool bIntersects = ITAGeoUtils::RayFaceIntersectionTest(rRay, pMesh, hFace, v3IntersectionPoint);
+				CITAMesh::FaceHandle hFace( *it++ );
 
-				/*
-				if (bIntersects)
+				// Check for intersection of the ray with current face
+				bFaceIntersectedByRay = ITAGeoUtils::RayFaceIntersectionTest(rRay, pMesh, hFace, v3ReflectionPoint);
+				if (bFaceIntersectedByRay)
 				{
-					//Set reflection anchor at intersection point
-					std::shared_ptr< ITAGeo::CSpecularReflection > pReflectionAnchor = std::make_shared < ITAGeo::CSpecularReflection >();
-					pReflectionAnchor->v3InteractionPoint = v3IntersectionPoint;
+					// Reflect ray on face
+					ReflectRayOnFace(rRay, pMesh, hFace, v3ReflectionPoint); //@todo: should this function be a method of class CEngine?
 				
-					//Add reflection to path
-					lPropagationPaths[pNode->pImage->iOrder] = pReflectionAnchor;
+					// Update path
+					auto pReflection = std::make_shared<ITAGeo::CSpecularReflection>(v3ReflectionPoint);
+					pReflection->v3InteractionPoint = v3ReflectionPoint;
+					oPropagationPath.push_back(pReflection);
 
-					//Add further propagation anchors to path
-					if (pNode->HasParent())
-					{
-						ConstructReflectionsPath(pNode->pParent, pReflectionAnchor, lPropagationPaths, eCulling);
-					}
+					// Update absorbtion losses
+					tAbsorbtionLossFactor *= (1.0f - tAbsorbtionOfFace); // @todo pReflection->Material = get material from face
+
+					// Update propagation losses
+					tPropagationLossFactor = 1 / powf(oPropagationPath.GetLength(), 2);
 				}				
-				else
-				{
-					//No valid path
-					lPropagationPaths.clear();
 			}
 
-				*/
-
-				if (bIntersects)
+			// Check for recever hit by ray
+			bDetectionSphereHitByRay = RayDetectionSphereIntersectionTest(rRay, pDestination, m_fDetectionSphereRadius, v3DetectionPoint);
+			if (bDetectionSphereHitByRay)
 			{
-					// create instance reflection
-					auto pReflection = std::make_shared<ITAGeo::CSpecularReflection>(v3IntersectionPoint);
-					// @todo pReflection->Material = get material from face
-
-					bFlagValidConstellation = true;
+				// Store path of ray in path list for histogram evaluation
+				lPropagationPaths.push_back(oPropagationPath);
 
-					// consider absorbtion for each frequency band
-					for (int i = 0; i < iNumberOfFreqBands; i++)
-					{
-						tEnergy[i] *= 1.0f - tAbsorbtionInOctaveBands[i]; // @todo: get from face -> absorption material
+				// Update path
+				auto pDetection = std::make_shared<ITAGeo::CSensor>(v3DetectionPoint);
+				oPropagationPath.push_back(pDetection);
 
-						// update abortion criteria according energy threshold
-						bTotalEnergyOfBandsBelowThreshold = !bTotalEnergyOfBandsBelowThreshold && (tEnergy[i] < m_tEnergyThreshold);
+				if (bFaceIntersectedByRay)
+					// Update propagation losses
+					tPropagationLossFactor = 1 / powf(oPropagationPath.GetLength(), 2);
 			}
 
-					// Proceed with next ray if energy of every freq band is below threshold
-					if (bTotalEnergyOfBandsBelowThreshold)
-						break;
-
-					// reflect ray on face
-					assert(pMesh->has_face_normals());
-					CITAMesh::Normal normal(pMesh->normal(hFace));
-					const VistaVector3D v3FaceNormal(normal.data());
-
-					VistaVector3D vPerpendicularDirComponent = rRay.GetDir().Dot(v3FaceNormal) * v3FaceNormal; // direction component of incident ray perpendicular to reflection face
+			// Update energy of ray
+			tEnergy = tAbsorbtionLossFactor * tPropagationLossFactor;
+			bEnergyBelowThreshold = tEnergy < tEnergyThreshold;
 
-					rRay.SetDir(rRay.GetDir - 2 * vPerpendicularDirComponent);
-					rRay.SetOrigin(v3IntersectionPoint);
-
-					// Update path
-					pReflection->v3InteractionPoint = v3IntersectionPoint;
-					oPropagationPath.push_back(pReflection);
-					
-					// Recever hit by ray?
-					bool bReceiverHit = RayDetectionSphereIntersectionTest(rRay, pDestination, m_fDetectionSphereRadius);
-
-					if (bReceiverHit)
-					{
-						// Store propagation path in path list and mark contribution of ray in energy histogram
-						lPropagationPaths.push_back(oPropagationPath);
-						// @todo: store energy of ray in histogram
-						// @todo: absorbtion of ray by rcv?
-					}
-				}				
-			}
-			assert(bFlagValidConstellation);
+			// Reset flags
+			bFaceIntersectedByRay = false;
+			bDetectionSphereHitByRay = false;
 		}
 	}
 }