Commit 728289fb authored by Lukas Vollmer's avatar Lukas Vollmer
Browse files

Added basic Triangulation class + Tests for further developments

parent 4ff2fc4d
......@@ -13,6 +13,7 @@ vista_use_package( IPP QUIET )
vista_use_package( SimpleIni QUIET )
vista_use_package( SPLINE REQUIRED )
vista_find_package( libjson OPTIONAL QUIET )
vista_use_package( Eigen REQUIRED )
if ( NOT WIN32 )
vista_use_package ( Threads REQUIRED )
......
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2021
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#include <vector>
#include <VistaBase/VistaVector3D.h>
#include <VistaMath/VistaGeometries.h>
#include <ITABaseDefinitions.h>
namespace ITABase
{
namespace Math
{
class ITA_BASE_API CTriangulation
{
public:
CTriangulation(const std::vector< VistaVector3D >& vvLoudspeakerPos, const VistaVector3D& vReproductionCenterPos);
CTriangulation(const std::vector< VistaVector3D >& vvLoudspeakerPos, const std::vector< VistaVector3D >& vvTriangulationIndices, const VistaVector3D& vReproductionCenterPos = VistaVector3D(0, 0, 0));
~CTriangulation();
VistaTriangle* GetTriangleIntersectedBySourceDirection(const VistaVector3D& vSourceDir);
bool IsDelaunayTriangulation(float epsilon = 1.0e-4);
private:
/* new Implementation @LVO */
std::vector< VistaVector3D > m_vvLoudspeakerPos;
/* Implementation from VBAP Renderer*/
struct sLoudspeaker
{
std::string sIdentifier;
int iIdentifier;
int iChannel;
VistaVector3D pos;
};
struct sSection
{
int iIdentifier;
std::vector <int> iLSIdentifier;
VistaPolygon Polygon;
};
VistaVector3D m_vCenterPosition;
std::vector< sLoudspeaker > m_vLoudspeaker;
std::vector< sSection > m_vSection;
std::vector< VistaVector3D > m_vvTriangulationIndices;
std::vector< VistaTriangle > m_vTriangulation; //!< Triangulation all loudspeaker for the VBAP setup; Entries are the pos of the speaker safed in a VistaTriangle struct
std::vector< std::vector< sLoudspeaker > > m_vvTriangulationLoudspeaker; //!< Triangulation all loudspeaker for the VBAP setup; Entries are the actual loudspeaker according to m_voLoudspeaker
void Triangulate();
};
}
}
......@@ -8,6 +8,7 @@ set( DirFiles
LinearAlgebra.h
PiecewisePolynomial.h
Spline.h
Triangulation.h
)
......
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2021
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#include <algorithm>
#include <ITAFastMath.h>
#include <ITABase/Math/Triangulation.h>
using namespace ITABase;
using namespace ITABase::Math;
/* Constructors */
CTriangulation::CTriangulation(const std::vector< VistaVector3D >& vvLoudspeakerPos, const VistaVector3D& vReproductionCenterPos)
{
}
CTriangulation::CTriangulation(const std::vector< VistaVector3D >& vvLoudspeakerPos, const std::vector< VistaVector3D >& vvTriangulationIndices, const VistaVector3D& vReproductionCenterPos)
:m_vvTriangulationIndices(vvTriangulationIndices), m_vCenterPosition(vReproductionCenterPos)
{
// Reference LS positions to vReproductionCenterPos, so that the new center is (0, 0, 0)
VistaVector3D vReferencedLoudspeakerPos;
for (int i = 0; i < vvLoudspeakerPos.size(); i++)
{
vReferencedLoudspeakerPos = vvLoudspeakerPos[i] - vReproductionCenterPos;
//vReferencedLoudspeakerPos.Normalize();
m_vvLoudspeakerPos.push_back(vReferencedLoudspeakerPos);
}
// Construct triangles from given Triangulation
VistaVector3D currentTriangle;
for (int i = 0; i < m_vvTriangulationIndices.size(); i++)
{
currentTriangle = m_vvTriangulationIndices[i];
m_vTriangulation.push_back(VistaTriangle(m_vvLoudspeakerPos[currentTriangle[0]], m_vvLoudspeakerPos[currentTriangle[1]], m_vvLoudspeakerPos[currentTriangle[2]]));
}
}
/* Destructor */
CTriangulation::~CTriangulation() {}
/* Methods */
VistaTriangle* CTriangulation::GetTriangleIntersectedBySourceDirection(const VistaVector3D& vSourceDir)
{
// Get closest speaker
VistaVector3D vRelSourceDir = vSourceDir - m_vCenterPosition;
double dCurrentDistance;
double dMinDistance;
int iClosestSpeaker = -1;
for (int i = 0; i < m_vvLoudspeakerPos.size(); i++)
{
dCurrentDistance = (m_vvLoudspeakerPos[i] - vRelSourceDir).GetLength();
if (i == 0 || dCurrentDistance < dMinDistance)
{
dMinDistance = dCurrentDistance;
iClosestSpeaker = i;
}
}
// Get triangles that contain closest speaker
VistaTriangle* vtIntersectedTriangle;
for (int i = 0; i < m_vTriangulation.size(); i++)
{
if (m_vvTriangulationIndices[i][0] == iClosestSpeaker || m_vvTriangulationIndices[i][1] == iClosestSpeaker || m_vvTriangulationIndices[i][2] == iClosestSpeaker)
{
float bary_ab, bary_ac;
if (m_vTriangulation[i].IntersectionPoint(vRelSourceDir, bary_ab, bary_ac, 1e-3))
return &m_vTriangulation[i];
}
}
return nullptr;
}
float calculateCircumcenter3D(VistaTriangle& tTriangle, VistaVector3D& vCenter)
{
// Taken from https://gamedev.stackexchange.com/questions/60630/how-do-i-find-the-circumcenter-of-a-triangle-in-3d
/*
* Let A, B, C, M be points in R and T = {A, B, C} be a Triangle. Then the circumcenter M of T is calculated as
*
* ||C-A|| * [(B-A) x (C-A)] x (B-A) + ||B-A|| * (C-A) x [(B-A) x (C-A)]
* M = A + -------------------------------------------------------------------------.
* 2 * ||(B-A) x (C-A)||
*
*
*/
VistaVector3D vA2B = tTriangle.GetB() - tTriangle.GetA();
VistaVector3D vA2C = tTriangle.GetC() - tTriangle.GetA();
VistaVector3D vABxAC = vA2B.Cross(vA2C);
VistaVector3D vA2M = (vABxAC.Cross(vA2B) * vA2C.GetLengthSquared() + vA2C.Cross(vABxAC) * vA2B.GetLengthSquared()) / (2 * vABxAC.GetLengthSquared());
vCenter = tTriangle.GetA() + vA2M;
return vA2M.GetLength();
}
float sphereEquation(VistaVector3D& vCenter, VistaVector3D& vPoint)
{
/*
* r (<)= (x - m_x) + (y - m_y) + (z - m_z)
* = ||P - M||
*/
return (vPoint - vCenter).GetLength();
}
bool CTriangulation::IsDelaunayTriangulation(float epsilon)
{
VistaVector3D vCenter;
float fRadius;
bool bBelongsToTriangle;
for (int i = 0; i < m_vTriangulation.size(); ++i)
{
fRadius = calculateCircumcenter3D(m_vTriangulation[i], vCenter);
for (int j = 0; j < m_vvLoudspeakerPos.size(); ++j)
{
bBelongsToTriangle = m_vvTriangulationIndices[i][1] == j || m_vvTriangulationIndices[i][2] == j || m_vvTriangulationIndices[i][3] == j;
if (~bBelongsToTriangle && (sphereEquation(vCenter, m_vvLoudspeakerPos[j]) - fRadius) < -epsilon)
return false;
}
}
return true;
}
......@@ -8,6 +8,7 @@ set( DirFiles
LinearAlgebra.cpp
PiecewisePolynomial.cpp
Spline.cpp
Triangulation.cpp
)
......
......@@ -94,4 +94,24 @@ vista_create_default_info_file( SplineTest )
set_property( TARGET SplineTest PROPERTY FOLDER "ITACoreLibs/Tests/ITABase" )
add_executable( NumericUtilsTest NumericUtilsTest.cpp )
target_link_libraries( NumericUtilsTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( NumericUtilsTest )
vista_install( NumericUtilsTest )
vista_create_default_info_file( NumericUtilsTest )
set_property( TARGET NumericUtilsTest PROPERTY FOLDER "ITACoreLibs/Tests/ITABase" )
add_executable( TriangulationTest TriangulationTest.cpp )
target_link_libraries( TriangulationTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( TriangulationTest )
vista_install( TriangulationTest )
vista_create_default_info_file( TriangulationTest )
set_property( TARGET TriangulationTest PROPERTY FOLDER "ITACoreLibs/Tests/ITABase" )
add_subdirectory( "VistaTests" )
#include <ITANumericUtils.h>
#include <ITAConstants.h>
#include <VistaBase/VistaVector3D.h>
#include <vector>
#include <math.h>
#include "Eigen/Dense"
using namespace std;
double calcAzimuthDirection(VistaVector3D& view, VistaVector3D& up, VistaVector3D dir)
{
dir.Normalize();
VistaVector3D right = up.Cross(view);
const double azimuth = atan2(dir.Dot(right), dir.Dot(view));
return (azimuth < 0 ? azimuth + 2*ITAConstants::PI_D : azimuth);
}
double calcZenithDirection(VistaVector3D& view, VistaVector3D& up, VistaVector3D dir)
{
dir.Normalize();
const double zenith = acos(dir.Dot(up));
return zenith;
}
// ------------------------ Tests start here ------------------------
void testSphericalHarmonics(VistaVector3D& view, VistaVector3D& up, vector<VistaVector3D>& vv3LsPos)
{
const int N = floor(sqrt(vv3LsPos.size()) - 1),
hoaChannels = (N + 1) * (N + 1);
Eigen::MatrixXd matY(hoaChannels, vv3LsPos.size());
Eigen::VectorXd vecMagnitude(hoaChannels);
vecMagnitude.setZero();
for (int i = 0; i < vv3LsPos.size(); i++)
{
VistaVector3D v3LsPos = vv3LsPos[i];
double phi = calcAzimuthDirection(view, up, v3LsPos),
theta = calcZenithDirection(view, up, v3LsPos);
vector< double > Y_i = SHRealvaluedBasefunctions(theta, phi, N);
for (int j = 0; j < Y_i.size(); j++)
{
matY(j, i) = Y_i[j];
vecMagnitude[j] += Y_i[j] * Y_i[j];
}
}
vecMagnitude /= vv3LsPos.size();
if (abs(vecMagnitude.mean() - 1.0 / (4 * ITAConstants::PI_D)) > 1e-3)
{
cout << "Error: Mean magnitude for Williams normalization deviates from 1/(4*pi) by more than 1e-3" << endl;
}
}
int main( int, char** )
{
VistaVector3D up(0, 1, 0),
view(0, 0, -1);
vector<VistaVector3D> vv3LsPos{ // 4-Design, 14 points
VistaVector3D(0, 0, -1),
VistaVector3D(-0.17765, 0.78406, -0.59472),
VistaVector3D(0.76784, -0.23818, -0.59472),
VistaVector3D(-0.59019, -0.54588, -0.59472),
VistaVector3D(-0.17765, -0.78406, 0.59472),
VistaVector3D(-0.59019, 0.54588, 0.59472),
VistaVector3D(0.76784, 0.23818, 0.59472),
VistaVector3D(-0.87947, 0.36847, -0.30125),
VistaVector3D(0.75884, 0.57741, -0.30125),
VistaVector3D(0.12063, -0.94588, -0.30125),
VistaVector3D(-0.87947, -0.36847, 0.30125),
VistaVector3D(0.12063, 0.94588, 0.30125),
VistaVector3D(0.75884, -0.57741, 0.30125),
VistaVector3D(-1.2246e-16, 0, 1)
};
cout << "--- Testing Spherical Harmonics ---" << endl;
testSphericalHarmonics(view, up, vv3LsPos);
cout << endl << endl << "--- All tests finished ---";
}
#include <ITAStringUtils.h>
#include <iostream>
using namespace std;
int main( int, char** )
{
std::vector< int > v = StringToIntVec( "1,2,3,4" );
for( int n : v )
cout << n << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment