13
0

VBAP backend re-work (part two): speaker positioning

* clean up source (whitespace)
* fix speaker 3x3 matrix iteration
* update math to go along with Ardour Cartesian -- fixes rounding errors
* fix division by zero in cross_prod()
* disable old debug output

(NB PBD::spherical_to_cartesian() returns
  3.7494e-33, 6.12323e-17, 1 for azimuth 90 elevation 90 distance 1
while it should return
  0.000000, 0.000000, 1 for azimuth 90 elevation 90 distance 1
IOW  cos(90.0 * 2.0 * M_PI / 360.0) != 0
Cause unknown. This is currently worked around check in vec_length()
)
This commit is contained in:
Robin Gareus 2014-01-11 23:29:23 +01:00 committed by Paul Davis
parent 112de00841
commit 1bf9c4c990
2 changed files with 130 additions and 122 deletions

View File

@ -1,4 +1,4 @@
/*
/*
This software is being provided to you, the licensee, by Ville Pulkki,
under the following license. By obtaining, using and/or copying this
software, you agree that you have read, understood, and will comply
@ -10,15 +10,15 @@
the disclaimer, and that the same appear on ALL copies of the software
and documentation, including modifications that you make for internal
use or for distribution:
Copyright 1998 by Ville Pulkki, Helsinki University of Technology. All
rights reserved.
rights reserved.
The software may be used, distributed, and included to commercial
products without any charges. When included to a commercial product,
the method "Vector Base Amplitude Panning" and its developer Ville
Pulkki must be referred to in documentation.
This software is provided "as is", and Ville Pulkki or Helsinki
University of Technology make no representations or warranties,
expressed or implied. By way of example, but not limitation, Helsinki
@ -31,11 +31,6 @@
of the software.
*/
#ifdef COMPILER_MSVC
#pragma warning ( disable : 4244 )
#endif
#include <vector>
#include <cmath>
#include <algorithm>
#include <stdlib.h>
@ -50,13 +45,6 @@ using namespace std;
const double VBAPSpeakers::MIN_VOL_P_SIDE_LGTH = 0.01;
typedef std::vector<double> DoubleVector;
typedef std::vector<float> FloatVector;
typedef std::vector<bool> BoolVector;
typedef std::vector<int> IntVector;
typedef std::vector<IntVector> IntVector2D;
typedef std::vector<DoubleVector> DoubleVector2D;
VBAPSpeakers::VBAPSpeakers (boost::shared_ptr<Speakers> s)
: _dimension (2)
, _parent (s)
@ -102,36 +90,46 @@ VBAPSpeakers::update ()
}
}
void
VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
void
VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
{
/* Selects the loudspeaker triplets, and
calculates the inversion matrices for each selected triplet.
A line (connection) is drawn between each loudspeaker. The lines
denote the sides of the triangles. The triangles should not be
intersecting. All crossing connections are searched and the
denote the sides of the triangles. The triangles should not be
intersecting. All crossing connections are searched and the
longer connection is erased. This yields non-intesecting triangles,
which can be used in panning.
*/
#if 0 // DEVEL/DEBUG
for (vector<Speaker>::iterator i = _speakers.begin(); i != _speakers.end(); ++i) {
cout << "Speaker " << (*i).id << " @ "
<< (*i).coords().x << ", " << (*i).coords().y << ", " << (*i).coords().z
<< " azimuth " << (*i).angles().azi
<< " elevation " << (*i).angles().ele
<< " distance " << (*i).angles().length
<< endl;
}
#endif
int i,j,k,l,table_size;
int n_speakers = _speakers.size ();
int connections[n_speakers][n_speakers];
float distance_table[((n_speakers * (n_speakers - 1)) / 2)];
int distance_table_i[((n_speakers * (n_speakers - 1)) / 2)];
int distance_table_j[((n_speakers * (n_speakers - 1)) / 2)];
float distance;
struct ls_triplet_chain *trip_ptr, *prev, *tmp_ptr;
if (n_speakers < 1) {
if (n_speakers == 0) {
return;
}
FloatVector distance_table(((n_speakers * (n_speakers - 1)) / 2));
IntVector distance_table_i(((n_speakers * (n_speakers - 1)) / 2));
IntVector distance_table_j(((n_speakers * (n_speakers - 1)) / 2));
IntVector2D connections(n_speakers, IntVector(n_speakers));
float distance;
struct ls_triplet_chain *trip_ptr, *prev, *tmp_ptr;
for (i = 0; i < n_speakers; i++) {
for (j = i+1; j < n_speakers; j++) {
for(k=j+1;k<n_speakers;k++) {
if (vol_p_side_lgth(i,j, k, _speakers) > MIN_VOL_P_SIDE_LGTH){
for(k = j+1; k < n_speakers; k++) {
if (vol_p_side_lgth(i, j, k, _speakers) > MIN_VOL_P_SIDE_LGTH) {
connections[i][j]=1;
connections[j][i]=1;
connections[i][k]=1;
@ -145,13 +143,13 @@ VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
}
/*calculate distancies between all speakers and sorting them*/
table_size =(((n_speakers - 1) * (n_speakers)) / 2);
table_size =(((n_speakers - 1) * (n_speakers)) / 2);
for (i = 0; i < table_size; i++) {
distance_table[i] = 100000.0;
}
for (i = 0;i < n_speakers; i++) {
for (j = i+1; j < n_speakers; j++) {
for (i = 0;i < n_speakers; i++) {
for (j = i+1; j < n_speakers; j++) {
if (connections[i][j] == 1) {
distance = fabs(vec_angle(_speakers[i].coords(),_speakers[j].coords()));
k=0;
@ -180,8 +178,8 @@ VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
if (connections[fst_ls][sec_ls] == 1) {
for (j = 0; j < n_speakers; j++) {
for (k = j+1; k < n_speakers; k++) {
if ((j!=fst_ls) && (k != sec_ls) && (k!=fst_ls) && (j != sec_ls)){
if (lines_intersect(fst_ls, sec_ls, j,k) == 1){
if ((j != fst_ls) && (k != sec_ls) && (k != fst_ls) && (j != sec_ls)) {
if (lines_intersect(fst_ls, sec_ls, j, k) == 1){
connections[j][k] = 0;
connections[k][j] = 0;
}
@ -199,8 +197,8 @@ VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
i = trip_ptr->ls_nos[0];
j = trip_ptr->ls_nos[1];
k = trip_ptr->ls_nos[2];
if (connections[i][j] == 0 ||
connections[i][k] == 0 ||
if (connections[i][j] == 0 ||
connections[i][k] == 0 ||
connections[j][k] == 0 ||
any_ls_inside_triplet(i,j,k) == 1 ){
if (prev != 0) {
@ -222,7 +220,7 @@ VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
}
}
int
int
VBAPSpeakers::any_ls_inside_triplet(int a, int b, int c)
{
/* returns 1 if there is loudspeaker(s) inside given ls triplet */
@ -240,12 +238,12 @@ VBAPSpeakers::any_ls_inside_triplet(int a, int b, int c)
lp1 = &(_speakers[a].coords());
lp2 = &(_speakers[b].coords());
lp3 = &(_speakers[c].coords());
/* matrix inversion */
invdet = 1.0 / ( lp1->x * ((lp2->y * lp3->z) - (lp2->z * lp3->y))
- lp1->y * ((lp2->x * lp3->z) - (lp2->z * lp3->x))
+ lp1->z * ((lp2->x * lp3->y) - (lp2->y * lp3->x)));
invmx[0] = ((lp2->y * lp3->z) - (lp2->z * lp3->y)) * invdet;
invmx[3] = ((lp1->y * lp3->z) - (lp1->z * lp3->y)) * -invdet;
invmx[6] = ((lp1->y * lp2->z) - (lp1->z * lp2->y)) * invdet;
@ -255,7 +253,7 @@ VBAPSpeakers::any_ls_inside_triplet(int a, int b, int c)
invmx[2] = ((lp2->x * lp3->y) - (lp2->y * lp3->x)) * invdet;
invmx[5] = ((lp1->x * lp3->y) - (lp1->y * lp3->x)) * -invdet;
invmx[8] = ((lp1->x * lp2->y) - (lp1->y * lp2->x)) * invdet;
any_ls_inside = false;
for (i = 0; i < n_speakers; i++) {
if (i != a && i!=b && i != c) {
@ -278,7 +276,7 @@ VBAPSpeakers::any_ls_inside_triplet(int a, int b, int c)
}
void
void
VBAPSpeakers::add_ldsp_triplet(int i, int j, int k, struct ls_triplet_chain **ls_triplets)
{
/* adds i,j,k triplet to triplet chain*/
@ -286,7 +284,7 @@ VBAPSpeakers::add_ldsp_triplet(int i, int j, int k, struct ls_triplet_chain **ls
struct ls_triplet_chain *trip_ptr, *prev;
trip_ptr = *ls_triplets;
prev = 0;
while (trip_ptr != 0){
prev = trip_ptr;
trip_ptr = trip_ptr->next;
@ -306,51 +304,51 @@ VBAPSpeakers::add_ldsp_triplet(int i, int j, int k, struct ls_triplet_chain **ls
trip_ptr->ls_nos[2] = k;
}
float
double
VBAPSpeakers::vec_angle(CartesianVector v1, CartesianVector v2)
{
float inner= ((v1.x*v2.x + v1.y*v2.y + v1.z*v2.z)/
double inner= ((v1.x*v2.x + v1.y*v2.y + v1.z*v2.z)/
(vec_length(v1) * vec_length(v2)));
if (inner > 1.0) {
inner= 1.0;
inner = 1.0;
}
if (inner < -1.0) {
inner = -1.0;
}
return fabsf((float) acos((double) inner));
return fabs(acos(inner));
}
float
double
VBAPSpeakers::vec_length(CartesianVector v1)
{
return (sqrt(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z));
double rv = sqrt(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z);
if (rv > 1e-14) return rv;
return 0;
}
float
double
VBAPSpeakers::vec_prod(CartesianVector v1, CartesianVector v2)
{
return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
}
float
VBAPSpeakers::vol_p_side_lgth(int i, int j,int k, const vector<Speaker>& speakers)
double
VBAPSpeakers::vol_p_side_lgth(int i, int j, int k, const vector<Speaker>& speakers)
{
/* calculate volume of the parallelepiped defined by the loudspeaker
direction vectors and divide it with total length of the triangle sides.
direction vectors and divide it with total length of the triangle sides.
This is used when removing too narrow triangles. */
float volper, lgth;
double volper, lgth;
CartesianVector xprod;
cross_prod (speakers[i].coords(), speakers[j].coords(), &xprod);
volper = fabsf (vec_prod(xprod, speakers[k].coords()));
lgth = (fabsf (vec_angle(speakers[i].coords(), speakers[j].coords()))
+ fabsf (vec_angle(speakers[i].coords(), speakers[k].coords()))
+ fabsf (vec_angle(speakers[j].coords(), speakers[k].coords())));
volper = fabs (vec_prod(xprod, speakers[k].coords()));
lgth = ( fabs (vec_angle(speakers[i].coords(), speakers[j].coords()))
+ fabs (vec_angle(speakers[i].coords(), speakers[k].coords()))
+ fabs (vec_angle(speakers[j].coords(), speakers[k].coords())));
if (lgth > 0.00001) {
return volper / lgth;
} else {
@ -358,28 +356,34 @@ VBAPSpeakers::vol_p_side_lgth(int i, int j,int k, const vector<Speaker>& speaker
}
}
void
VBAPSpeakers::cross_prod(CartesianVector v1,CartesianVector v2, CartesianVector *res)
void
VBAPSpeakers::cross_prod(CartesianVector v1,CartesianVector v2, CartesianVector *res)
{
float length;
double length;
res->x = (v1.y * v2.z) - (v1.z * v2.y);
res->y = (v1.z * v2.x) - (v1.x * v2.z);
res->z = (v1.x * v2.y) - (v1.y * v2.x);
res->x = (v1.y * v2.z ) - (v1.z * v2.y);
res->y = (v1.z * v2.x ) - (v1.x * v2.z);
res->z = (v1.x * v2.y ) - (v1.y * v2.x);
length = vec_length(*res);
res->x /= length;
res->y /= length;
res->z /= length;
if (length > 0) {
res->x /= length;
res->y /= length;
res->z /= length;
} else {
res->x = 0;
res->y = 0;
res->z = 0;
}
}
int
int
VBAPSpeakers::lines_intersect (int i, int j, int k, int l)
{
/* checks if two lines intersect on 3D sphere
/* checks if two lines intersect on 3D sphere
see theory in paper Pulkki, V. Lokki, T. "Creating Auditory Displays
with Multiple Loudspeakers Using VBAP: A Case Study with
DIVA Project" in International Conference on
DIVA Project" in International Conference on
Auditory Displays -98. E-mail Ville.Pulkki@hut.fi
if you want to have that paper.
*/
@ -389,11 +393,11 @@ VBAPSpeakers::lines_intersect (int i, int j, int k, int l)
CartesianVector v3, neg_v3;
float dist_ij,dist_kl,dist_iv3,dist_jv3,dist_inv3,dist_jnv3;
float dist_kv3,dist_lv3,dist_knv3,dist_lnv3;
cross_prod(_speakers[i].coords(),_speakers[j].coords(),&v1);
cross_prod(_speakers[k].coords(),_speakers[l].coords(),&v2);
cross_prod(v1,v2,&v3);
neg_v3.x= 0.0 - v3.x;
neg_v3.y= 0.0 - v3.y;
neg_v3.z= 0.0 - v3.z;
@ -410,15 +414,14 @@ VBAPSpeakers::lines_intersect (int i, int j, int k, int l)
dist_lnv3 = (vec_angle(neg_v3,_speakers[l].coords()));
/* if one of loudspeakers is close to crossing point, don't do anything*/
if(fabsf(dist_iv3) <= 0.01 || fabsf(dist_jv3) <= 0.01 ||
if(fabsf(dist_iv3) <= 0.01 || fabsf(dist_jv3) <= 0.01 ||
fabsf(dist_kv3) <= 0.01 || fabsf(dist_lv3) <= 0.01 ||
fabsf(dist_inv3) <= 0.01 || fabsf(dist_jnv3) <= 0.01 ||
fabsf(dist_inv3) <= 0.01 || fabsf(dist_jnv3) <= 0.01 ||
fabsf(dist_knv3) <= 0.01 || fabsf(dist_lnv3) <= 0.01 ) {
return(0);
}
/* if crossing point is on line between both loudspeakers return 1 */
if (((fabsf(dist_ij - (dist_iv3 + dist_jv3)) <= 0.01 ) &&
(fabsf(dist_kl - (dist_kv3 + dist_lv3)) <= 0.01)) ||
((fabsf(dist_ij - (dist_inv3 + dist_jnv3)) <= 0.01) &&
@ -429,9 +432,9 @@ VBAPSpeakers::lines_intersect (int i, int j, int k, int l)
}
}
void
void
VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
{
{
/* Calculates the inverse matrices for 3D */
float invdet;
const CartesianVector* lp1;
@ -443,7 +446,7 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
int triplet;
assert (tr_ptr);
/* counting triplet amount */
while (tr_ptr != 0) {
@ -451,7 +454,9 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
tr_ptr = tr_ptr->next;
}
cerr << "@@@ triplets generate " << triplet_count << " of speaker tuples\n";
#if 0 // DEVEL/DEBUG
cerr << "@@@ VBAP triplets generate " << triplet_count << " of speaker tuples\n";
#endif
triplet = 0;
@ -463,17 +468,18 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
_speaker_tuples.push_back (tmatrix());
}
tr_ptr = ls_triplets;
while (tr_ptr != 0) {
lp1 = &(_speakers[tr_ptr->ls_nos[0]].coords());
lp2 = &(_speakers[tr_ptr->ls_nos[1]].coords());
lp3 = &(_speakers[tr_ptr->ls_nos[2]].coords());
/* matrix inversion */
invmx = tr_ptr->inv_mx;
invdet = 1.0 / ( lp1->x * ((lp2->y * lp3->z) - (lp2->z * lp3->y))
- lp1->y * ((lp2->x * lp3->z) - (lp2->z * lp3->x))
+ lp1->z * ((lp2->x * lp3->y) - (lp2->y * lp3->x)));
invmx[0] = ((lp2->y * lp3->z) - (lp2->z * lp3->y)) * invdet;
invmx[3] = ((lp1->y * lp3->z) - (lp1->z * lp3->y)) * -invdet;
invmx[6] = ((lp1->y * lp2->z) - (lp1->z * lp2->y)) * invdet;
@ -483,7 +489,7 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
invmx[2] = ((lp2->x * lp3->y) - (lp2->y * lp3->x)) * invdet;
invmx[5] = ((lp1->x * lp3->y) - (lp1->y * lp3->x)) * -invdet;
invmx[8] = ((lp1->x * lp2->y) - (lp1->y * lp2->x)) * invdet;
/* copy the matrix */
_matrices[triplet][0] = invmx[0];
@ -500,10 +506,12 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
_speaker_tuples[triplet][1] = tr_ptr->ls_nos[1];
_speaker_tuples[triplet][2] = tr_ptr->ls_nos[2];
cerr << "Triplet[" << triplet << "] = "
<< tr_ptr->ls_nos[0] << " + "
<< tr_ptr->ls_nos[1] << " + "
#if 0 // DEVEL/DEBUG
cerr << "Triplet[" << triplet << "] = "
<< tr_ptr->ls_nos[0] << " + "
<< tr_ptr->ls_nos[1] << " + "
<< tr_ptr->ls_nos[2] << endl;
#endif
triplet++;
@ -511,55 +519,55 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
}
}
void
void
VBAPSpeakers::choose_speaker_pairs (){
/* selects the loudspeaker pairs, calculates the inversion
matrices and stores the data to a global array
*/
const int n_speakers = _speakers.size();
if (n_speakers < 1) {
return;
}
IntVector sorted_speakers(n_speakers);
BoolVector exists(n_speakers);
DoubleVector2D inverse_matrix(n_speakers, DoubleVector(4));
const double AZIMUTH_DELTA_THRESHOLD_DEGREES = (180.0/M_PI) * (M_PI - 0.175);
const double AZIMUTH_DELTA_THRESHOLD_DEGREES = (180.0/M_PI) * (M_PI - 0.175);
int sorted_speakers[n_speakers];
bool exists[n_speakers];
double inverse_matrix[n_speakers][4];
int expected_pairs = 0;
int pair;
int speaker;
if (n_speakers == 0) {
return;
}
for (speaker = 0; speaker < n_speakers; ++speaker) {
exists[speaker] = false;
}
/* sort loudspeakers according their aximuth angle */
sort_2D_lss (&sorted_speakers[0]);
sort_2D_lss (sorted_speakers);
/* adjacent loudspeakers are the loudspeaker pairs to be used.*/
for (speaker = 0; speaker < n_speakers-1; speaker++) {
if ((_speakers[sorted_speakers[speaker+1]].angles().azi -
if ((_speakers[sorted_speakers[speaker+1]].angles().azi -
_speakers[sorted_speakers[speaker]].angles().azi) <= AZIMUTH_DELTA_THRESHOLD_DEGREES) {
if (calc_2D_inv_tmatrix( _speakers[sorted_speakers[speaker]].angles().azi,
_speakers[sorted_speakers[speaker+1]].angles().azi,
&inverse_matrix[speaker][0]) != 0){
if (calc_2D_inv_tmatrix( _speakers[sorted_speakers[speaker]].angles().azi,
_speakers[sorted_speakers[speaker+1]].angles().azi,
inverse_matrix[speaker]) != 0){
exists[speaker] = true;
expected_pairs++;
}
}
}
if (((6.283 - _speakers[sorted_speakers[n_speakers-1]].angles().azi)
if (((6.283 - _speakers[sorted_speakers[n_speakers-1]].angles().azi)
+_speakers[sorted_speakers[0]].angles().azi) <= AZIMUTH_DELTA_THRESHOLD_DEGREES) {
if (calc_2D_inv_tmatrix(_speakers[sorted_speakers[n_speakers-1]].angles().azi,
_speakers[sorted_speakers[0]].angles().azi,
&inverse_matrix[n_speakers-1][0]) != 0) {
if (calc_2D_inv_tmatrix(_speakers[sorted_speakers[n_speakers-1]].angles().azi,
_speakers[sorted_speakers[0]].angles().azi,
inverse_matrix[n_speakers-1]) != 0) {
exists[n_speakers-1] = true;
expected_pairs++;
}
}
}
pair = 0;
@ -585,7 +593,7 @@ VBAPSpeakers::choose_speaker_pairs (){
pair++;
}
}
if (exists[n_speakers-1]) {
_matrices[pair][0] = inverse_matrix[speaker][0];
_matrices[pair][1] = inverse_matrix[speaker][1];
@ -597,7 +605,7 @@ VBAPSpeakers::choose_speaker_pairs (){
}
}
void
void
VBAPSpeakers::sort_2D_lss (int* sorted_speakers)
{
vector<Speaker> tmp = _speakers;
@ -612,7 +620,7 @@ VBAPSpeakers::sort_2D_lss (int* sorted_speakers)
}
}
int
int
VBAPSpeakers::calc_2D_inv_tmatrix (double azi1, double azi2, double* inverse_matrix)
{
double x1,x2,x3,x4;
@ -625,7 +633,7 @@ VBAPSpeakers::calc_2D_inv_tmatrix (double azi1, double azi2, double* inverse_mat
det = (x1 * x4) - ( x3 * x2 );
if (fabs(det) <= 0.001) {
inverse_matrix[0] = 0.0;
inverse_matrix[1] = 0.0;
inverse_matrix[2] = 0.0;

View File

@ -84,11 +84,11 @@ private:
struct ls_triplet_chain *next;
};
static float vec_angle(PBD::CartesianVector v1, PBD::CartesianVector v2);
static float vec_length(PBD::CartesianVector v1);
static float vec_prod(PBD::CartesianVector v1, PBD::CartesianVector v2);
static float vol_p_side_lgth(int i, int j,int k, const std::vector<Speaker>&);
static void cross_prod(PBD::CartesianVector v1,PBD::CartesianVector v2, PBD::CartesianVector *res);
static double vec_angle(PBD::CartesianVector v1, PBD::CartesianVector v2);
static double vec_length(PBD::CartesianVector v1);
static double vec_prod(PBD::CartesianVector v1, PBD::CartesianVector v2);
static double vol_p_side_lgth(int i, int j,int k, const std::vector<Speaker>&);
static void cross_prod(PBD::CartesianVector v1,PBD::CartesianVector v2, PBD::CartesianVector *res);
void update ();
int any_ls_inside_triplet (int a, int b, int c);