From 1bf9c4c990ccc7f308a38a2eba91fbdb4bb93292 Mon Sep 17 00:00:00 2001 From: Robin Gareus Date: Sat, 11 Jan 2014 23:29:23 +0100 Subject: [PATCH] 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() ) --- libs/panners/vbap/vbap_speakers.cc | 242 +++++++++++++++-------------- libs/panners/vbap/vbap_speakers.h | 10 +- 2 files changed, 130 insertions(+), 122 deletions(-) diff --git a/libs/panners/vbap/vbap_speakers.cc b/libs/panners/vbap/vbap_speakers.cc index f2533222bd..313fe7a5cd 100644 --- a/libs/panners/vbap/vbap_speakers.cc +++ b/libs/panners/vbap/vbap_speakers.cc @@ -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 #include #include #include @@ -50,13 +45,6 @@ using namespace std; const double VBAPSpeakers::MIN_VOL_P_SIDE_LGTH = 0.01; -typedef std::vector DoubleVector; -typedef std::vector FloatVector; -typedef std::vector BoolVector; -typedef std::vector IntVector; -typedef std::vector IntVector2D; -typedef std::vector DoubleVector2D; - VBAPSpeakers::VBAPSpeakers (boost::shared_ptr 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::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 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& speakers) +double +VBAPSpeakers::vol_p_side_lgth(int i, int j, int k, const vector& 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 } } -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 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; diff --git a/libs/panners/vbap/vbap_speakers.h b/libs/panners/vbap/vbap_speakers.h index b2f8b3c9dd..c3e90ce646 100644 --- a/libs/panners/vbap/vbap_speakers.h +++ b/libs/panners/vbap/vbap_speakers.h @@ -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&); - 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&); + 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);