Merge branch 'ma_2018/erraji' into develop

include/ITAGeo/Base.h

@@ -311,6 +311,7 @@ namespace ITAGeo
 		inline CSpecularReflection() 
 		{
 			iAnchorType = CPropagationAnchor::SPECULAR_REFLECTION;
+			v3FaceNormal.SetToZeroVector();
 		};
 		inline CSpecularReflection(const VistaVector3D& v3SpecularReflectionPoint)
 			: CPropagationAnchor(v3SpecularReflectionPoint)
@@ -327,6 +328,8 @@ namespace ITAGeo
 		//! Load properties from JSON formatted string
 		void LoadFromJSONString(const std::string& sContent, std::shared_ptr< CResourceManager > pResourceManager = nullptr);
 
+		VistaVector3D v3FaceNormal; //!< Face normal of reflection plane
+
 		std::shared_ptr< Material::IMaterial > pMaterial; //!< Acoustic material connection, i.e. absorption coefficient
 	};
 
@@ -408,6 +411,39 @@ namespace ITAGeo
 			return (int) this->size();
 		};
 
+		//! Get number of reflections
+		inline int GetNumReflections() const
+		{
+
+			int iNumReflections = 0;
+			for (auto& pAnchor : *this)
+			{
+				if (pAnchor->iAnchorType == ITAGeo::CPropagationAnchor::SPECULAR_REFLECTION)
+					iNumReflections++;
+			}
+
+			return iNumReflections;
+		};
+
+
+		//! Get number of diffractions
+		inline int GetNumDiffractions() const
+		{
+
+			int iNumDiffractions = 0;
+			for (auto& pAnchor : *this)
+			{
+
+				if (pAnchor->iAnchorType == ITAGeo::CPropagationAnchor::DIFFRACTION_OUTER_APEX
+					|| pAnchor->iAnchorType == ITAGeo::CPropagationAnchor::DIFFRACTION_INNER_APEX)
+				{
+					iNumDiffractions++;
+				}
+			}
+
+			return iNumDiffractions;
+		};
+
 		//! Calculates accumulated angles over diffracting anchors in radiants
 		/**
 		  * @return Accumulated angle in radiants

include/ITAGeo/SketchUp/Model.h

@@ -176,7 +176,7 @@ namespace ITAGeo
 
 			std::map< std::string, std::map<std::string, ITAGeo::CPropagationAnchor > > m_mGeoEmittersVisualization;	//!< Propagation anchors for visualization of acoustic emitters
 			std::map< std::string, std::map<std::string, ITAGeo::CPropagationAnchor > > m_mGeoSensorsVisualization;	//!< Propagation anchors for visualization of acoustic sensors
-			std::map< std::string, std::vector< ITAGeo::CPropagationPath > > m_mGeoPropPathsVisualization;	//!< Propagation paths for visualization of acoustic effects (never populated in Load() method)
+			std::map< std::string, ITAGeo::CPropagationPathList > m_mGeoPropPathsVisualization;	//!< Propagation paths for visualization of acoustic effects (never populated in Load() method)
 			std::map< std::string, std::vector< const ITAGeo::Halfedge::CMeshModel* > > m_mAcousticVisualizationMeshes;		//!< Mesh list with extra visualization data
 
 			std::vector< std::string > m_vsLayers;	//!< Layer names (also unused layers)

include/ITAGeo/Utils.h

@@ -134,6 +134,9 @@ namespace ITAGeoUtils
 	  */
 	ITA_GEO_API bool IsDiffractionAperturePointInRange( const VistaVector3D& v3StartVertex, const VistaVector3D& v3EndVertex, const VistaVector3D& v3PointOnEdge );
 
+	//! Calculates the diffraction angle
+	ITA_GEO_API bool CalculateDiffractionAngle(const VistaVector3D& v3StartPoint, const VistaVector3D& v3AperturePoint, const VistaVector3D& v3EndPoint, float& dAngle);
+
 	//! Calculates the aperture point on an finite wedge and returns true if within bounds
 	/**
 	  * @param[out] v3AperturePoint

matlab/propagation_path_eval_Snells_law_angle_error.m

+%% Load
+
+json_file_path = 'CombinedModelTest_Test.skp.json';
+txt = fileread( json_file_path );
+json_content = jsondecode( txt );
+
+N = numel( json_content.propagation_paths );
+fprintf( 'Found %i propagation paths\n', N )
+
+
+%% Eval
+
+error_angle_refl_rad = [];
+error_angle_diffr_rad = [];
+
+for n = 1:N
+    prop_path = json_content.propagation_paths( n );
+    assert( strcmpi( prop_path.class, 'propagation_path' ) )
+    
+    for a = 2:( numel( prop_path.propagation_anchors ) - 1 )
+        
+        if iscell( prop_path.propagation_anchors( a - 1 ) )
+            prop_anchor_prev = prop_path.propagation_anchors{ a - 1 };
+        else
+            prop_anchor_prev = prop_path.propagation_anchors( a - 1 );
+        end
+        
+        if iscell( prop_path.propagation_anchors( a ) )
+            prop_anchor_curr = prop_path.propagation_anchors{ a };
+        else
+            prop_anchor_curr = prop_path.propagation_anchors( a );
+        end
+        
+        if iscell( prop_path.propagation_anchors( a + 1 ) )
+            prop_anchor_next = prop_path.propagation_anchors{ a + 1 };
+        else
+            prop_anchor_next = prop_path.propagation_anchors( a + 1 );
+        end
+        
+        prop_segment_incoming_vec = prop_anchor_curr.interaction_point - prop_anchor_prev.interaction_point;
+        prop_segment_outgoing_vec = prop_anchor_next.interaction_point - prop_anchor_curr.interaction_point;
+        
+        if strcmpi( prop_anchor_curr.anchor_type, 'specular_reflection' ) && isfield( prop_anchor_curr, 'face_normal' )
+            a1 = dot( prop_segment_incoming_vec, prop_anchor_curr.face_normal ) / norm( prop_segment_incoming_vec ) / norm( prop_anchor_curr.face_normal );
+            a2 = ( -1 ) * dot( prop_segment_outgoing_vec, prop_anchor_curr.face_normal ) / norm( prop_segment_outgoing_vec ) / norm( prop_anchor_curr.face_normal );
+            error_angle_refl_rad( end + 1 ) = abs( acos( a2 ) - acos( a1 ) );
+        elseif strcmpi( prop_anchor_curr.anchor_type, 'outer_edge_diffraction' ) || strcmpi( prop_anchor_curr.anchor_type, 'inner_edge_diffraction' )
+            % @todo get edge normal, validate angles
+            vertex_dir = prop_anchor_curr.vertex_end -  prop_anchor_curr.vertex_start;
+            a1 = dot( prop_segment_incoming_vec, vertex_dir ) / norm( prop_segment_incoming_vec ) / norm( vertex_dir );
+            a2 = dot( prop_segment_outgoing_vec , vertex_dir ) / norm( prop_segment_outgoing_vec ) / norm( vertex_dir );
+            error_angle_diffr_rad( end + 1 ) = abs( acos( a2 ) - acos( a1 ) );
+        end
+    end
+end
+
+
+%% Statistics
+
+error_angle_rad = [ error_angle_refl_rad error_angle_diffr_rad ];
+fprintf( 'Propagation path angular error based on Snells law: %.1frad (%.1frad)\n', mean( error_angle_rad ), std2( error_angle_rad ) )
+fprintf( '\tspecular reflections only: %.3frad (%.3frad)\n', mean( error_angle_refl_rad ), std2( error_angle_refl_rad ) )
+fprintf( '\tspecular diffractions only: %.3frad (%.3frad)\n', mean( error_angle_diffr_rad ), std2( error_angle_diffr_rad ) )
+
+plot( error_angle_diffr_rad )
\ No newline at end of file

src/ITAGeo/Base.cpp

@@ -522,6 +522,18 @@ std::stringstream ITAGeo::CSpecularReflection::ToJSON( std::shared_ptr< const CR
 		JSONNode jnClassType("class", GetTypeStr());
 		jnRoot.push_back(jnClassType);
 	}
+	
+	if (!v3FaceNormal.GetIsZeroVector())
+	{
+		//Set interaction point array
+		JSONNode jnFaceNormalArray(JSON_ARRAY);
+		jnFaceNormalArray.set_name("face_normal");
+		jnFaceNormalArray.push_back(JSONNode("x", v3FaceNormal[Vista::X]));
+		jnFaceNormalArray.push_back(JSONNode("y", v3FaceNormal[Vista::Y]));
+		jnFaceNormalArray.push_back(JSONNode("z", v3FaceNormal[Vista::Z]));
+		jnFaceNormalArray.push_back(JSONNode("w", v3FaceNormal[Vista::W]));
+		jnRoot.push_back(jnFaceNormalArray);
+	}
 
 	//Get the acoustic material 
 	if( pMaterial && pResourceManager )

src/ITAGeo/Material/Utils.cpp

@@ -2,6 +2,8 @@
 #include <ITAException.h>
 #include <cassert>
 
+#include <cassert>
+
 float ITAGeo::Material::Utils::EstimateFilterLengthSamples( std::shared_ptr<const ITAGeo::Material::IMaterial> pMaterial, const float fSampleRate, const float fSpeedOfSound )
 {
 	// Allow for the lowest frequency a phase shift of one cycle / return samples of one wave length

src/ITAGeo/SketchUp/Model.cpp

@@ -286,6 +286,11 @@ bool CModel::Store( const std::string& sSKPFilePath, bool bOverwite/*=true*/ ) c
 	SU_EXC( ITAMeshToSUEntities( pMesh, &rEntities ) );
 
 
+	//Set unit 
+	SUModelUnits rUnits;
+	SU_EXC(SUModelGetUnits(rModel, &rUnits));
+	rUnits = SUModelUnits::SUModelUnits_Meters;
+
 	// Add Layers
 
 	std::vector< SULayerRef > vrLayer;
@@ -389,7 +394,7 @@ bool CModel::Store( const std::string& sSKPFilePath, bool bOverwite/*=true*/ ) c
 	
 	for( auto key : m_mGeoPropPathsVisualization )
 	{
-		std::vector< ITAGeo::CPropagationPath > voGeoPropPathsVisualization( key.second );
+		ITAGeo::CPropagationPathList voGeoPropPathsVisualization( key.second );
 
 		// Propagation path visualization
 		SULayerRef rPropLayer = SU_INVALID;
@@ -409,10 +414,41 @@ bool CModel::Store( const std::string& sSKPFilePath, bool bOverwite/*=true*/ ) c
 			SUMaterialRef rMaterial = SU_INVALID;
 			SU_EXC( SUMaterialCreate( &rMaterial ) );
 			SUColor oSUColor;
-			oSUColor.red = SUByte( 255.0f ); // direktschall
-			oSUColor.green = SUByte( 0.0f ); // reflex, yellow (red+green) = diffr
-			oSUColor.blue = SUByte( 0.0f ); // late reverb
-			oSUColor.alpha = SUByte( 255.0f * ( 1.0f - min( (float) oPath.GetLength(), 1000.0f ) / 1000.0f ) ); // @todo get max first, then map from 0 to max length (10% alpha)
+			
+			
+			//Color and alpha values
+			float fRed, fGreen, fBlue, fAlpha;
+			//Set alpha value according to length of path with 10% as minimum value and 255 as maximum value
+			float fMaxLength = voGeoPropPathsVisualization.GetMaxLength();
+			fAlpha = max(0.1f, 1.0f - (float)oPath.GetLength() / fMaxLength);
+			fAlpha = min(255.0f, fAlpha*255.0f);
+
+
+			//The blue value is always zero(reserved for late reverberations)
+			fBlue = 0.0f;
+
+			//Direct path (red)
+			if (oPath.GetNumAnchors() == 2)
+			{
+				fRed = 220.0f, fGreen = 0.0f;
+			}
+			else
+			{
+				int iNumDiffractions = oPath.GetNumDiffractions();
+				int iNumReflections = oPath.GetNumReflections();
+				float fRelNumDiffractions = (float)iNumDiffractions / (float)(iNumDiffractions + iNumReflections);
+
+
+				fRed = fRelNumDiffractions*255.0f;
+				fGreen = 155.0f + fRelNumDiffractions * 100.0f;
+
+			}
+
+			oSUColor.alpha = SUByte(fAlpha ); 
+			oSUColor.red = SUByte(fRed); // direktschall
+			oSUColor.green = SUByte(fGreen); // reflex, yellow (red+green) = diffr
+			oSUColor.blue = SUByte(fBlue); // late reverb
+			
 			SU_EXC( SUMaterialSetColor( rMaterial, &oSUColor ) );
 			SUModelAddMaterials( rModel, 1, &rMaterial );
 

src/ITAGeo/Utils.cpp

@@ -601,6 +601,16 @@ bool ITAGeoUtils::IsDiffractionAperturePointInRange( const VistaVector3D& v3Star
 		return false;
 }
 
+bool ITAGeoUtils::CalculateDiffractionAngle(const VistaVector3D & v3StartPoint, const VistaVector3D & v3AperturePoint, const VistaVector3D & v3EndPoint, float & dAngle)
+{
+	const VistaVector3D v3StartToApexDir = (v3AperturePoint - v3StartPoint).GetNormalized();
+	const VistaVector3D v3ApexToEndDir = (v3EndPoint - v3AperturePoint).GetNormalized();
+	
+	dAngle = acos(v3StartToApexDir * v3ApexToEndDir);
+	
+	return true;
+}
+
 bool ITAGeoUtils::RojectPoint( const VistaVector3D& v3VertexStart, const VistaVector3D& v3VertexEnd, const VistaVector3D& v3TargetPlaneNormal, const VistaVector3D& v3AnyPoint, VistaVector3D& v3RojectedPoint )
 {
 	// Algorithm A: use closest point on planes, check and correct amibuity and scale for constant Euklidean distance
@@ -656,6 +666,7 @@ bool ITAGeoUtils::RojectPoint( const VistaVector3D& v3VertexStart, const VistaVe
 	*/
 }
 
+
 unsigned long int ITAGeoUtils::CalculateNumberOfImages( const int iFaceCount, const int iOrder )
 {
 	if( iOrder <= 0 )

tests/PropagationPathsTest.cpp

@@ -12,6 +12,7 @@ int main( int, char** )
 
 	auto pReflection = make_shared< CSpecularReflection >();
 	pReflection->v3InteractionPoint.SetValues(-0.5f, -1.0f, -2.0f);
+	pReflection->v3FaceNormal.SetValues(1.f, 0.f, 0.f, 0.f);
 
 	auto pReceiver = make_shared< CSensor >();
 	pReceiver->v3InteractionPoint.SetValues(-1.0f, -2.0f, -3.0f);