-/*
+/*
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
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
using namespace PBD;
using namespace std;
+const double VBAPSpeakers::MIN_VOL_P_SIDE_LGTH = 0.01;
+
VBAPSpeakers::VBAPSpeakers (boost::shared_ptr<Speakers> s)
: _dimension (2)
, _parent (s)
}
}
-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)];
+
+ if (n_speakers < 3) {
+ fprintf(stderr, "VBAP: at least 3 speakers need to be defined.");
+ return;
+ }
+
+ /* variable length arrays arrived in C99, became optional in C11, and
+ are only planned for C++14. Use alloca which is functionally
+ identical (but uglier to read).
+ */
+ int* connections = (int*) alloca (sizeof (int) * n_speakers * n_speakers);
+ float* distance_table = (float *) alloca (sizeof (float) * ((n_speakers * (n_speakers - 1)) / 2));
+ int* distance_table_i = (int *) alloca (sizeof (int) * ((n_speakers * (n_speakers - 1)) / 2));
+ int* distance_table_j = (int *) alloca (sizeof (int) * ((n_speakers * (n_speakers - 1)) / 2));
float distance;
struct ls_triplet_chain *trip_ptr, *prev, *tmp_ptr;
- if (n_speakers == 0) {
- return;
+ for (i = 0; i < n_speakers * n_speakers; i++) {
+ connections[i] = 0;
}
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){
- connections[i][j]=1;
- connections[j][i]=1;
- connections[i][k]=1;
- connections[k][i]=1;
- connections[j][k]=1;
- connections[k][j]=1;
+ for(k = j+1; k < n_speakers; k++) {
+ if (vol_p_side_lgth(i, j, k, _speakers) > MIN_VOL_P_SIDE_LGTH) {
+ connections[(i*n_speakers)+j]=1;
+ connections[(j*n_speakers)+i]=1;
+ connections[(i*n_speakers)+k]=1;
+ connections[(k*n_speakers)+i]=1;
+ connections[(j*n_speakers)+k]=1;
+ connections[(k*n_speakers)+j]=1;
add_ldsp_triplet(i,j,k,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++) {
- if (connections[i][j] == 1) {
+ for (i = 0;i < n_speakers; i++) {
+ for (j = i+1; j < n_speakers; j++) {
+ if (connections[(i*n_speakers)+j] == 1) {
distance = fabs(vec_angle(_speakers[i].coords(),_speakers[j].coords()));
k=0;
while(distance_table[k] < distance) {
for (i = 0; i < table_size; i++) {
int fst_ls = distance_table_i[i];
int sec_ls = distance_table_j[i];
- if (connections[fst_ls][sec_ls] == 1) {
+ if (connections[(fst_ls*n_speakers)+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){
- connections[j][k] = 0;
- connections[k][j] = 0;
+ 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*n_speakers)+k] = 0;
+ connections[(k*n_speakers)+j] = 0;
}
}
}
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 ||
- connections[j][k] == 0 ||
+ if (connections[(i*n_speakers)+j] == 0 ||
+ connections[(i*n_speakers)+k] == 0 ||
+ connections[(j*n_speakers)+k] == 0 ||
any_ls_inside_triplet(i,j,k) == 1 ){
if (prev != 0) {
prev->next = trip_ptr->next;
}
}
-int
+int
VBAPSpeakers::any_ls_inside_triplet(int a, int b, int c)
{
/* returns 1 if there is loudspeaker(s) inside given ls triplet */
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;
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) {
}
-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*/
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;
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;
- CartesianVector xprod;
+ 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 {
}
}
-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.
*/
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;
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) &&
}
}
-void
+void
VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
-{
+{
/* Calculates the inverse matrices for 3D */
float invdet;
const CartesianVector* lp1;
int triplet;
assert (tr_ptr);
-
+
/* counting triplet amount */
while (tr_ptr != 0) {
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;
_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;
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];
_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++;
}
}
-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 < 2) {
+ fprintf(stderr, "VBAP: at least 2 speakers need to be defined.");
+ return;
+ }
+
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];
+ /* variable length arrays arrived in C99, became optional in C11, and
+ are only planned for C++14. Use alloca which is functionally
+ identical (but uglier to read).
+ */
+ int* sorted_speakers = (int*) alloca (sizeof (int) * n_speakers);
+ bool* exists = (bool*) alloca (sizeof(bool) * n_speakers);
+ double* inverse_matrix = (double*) alloca (sizeof (double) * 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 */
+#ifdef __clang_analyzer__
+ // sort_2D_lss() assigns values to all of sorted_speakers
+ // "uninitialized value"
+ memset(sorted_speakers, 0, sizeof(*sorted_speakers));
+#endif
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){
+ if (calc_2D_inv_tmatrix( _speakers[sorted_speakers[speaker]].angles().azi,
+ _speakers[sorted_speakers[speaker+1]].angles().azi,
+ &inverse_matrix[4 * 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) {
+ if (calc_2D_inv_tmatrix(_speakers[sorted_speakers[n_speakers-1]].angles().azi,
+ _speakers[sorted_speakers[0]].angles().azi,
+ &inverse_matrix[4*(n_speakers-1)]) != 0) {
exists[n_speakers-1] = true;
expected_pairs++;
- }
+ }
}
pair = 0;
for (speaker = 0; speaker < n_speakers - 1; speaker++) {
if (exists[speaker]) {
- _matrices[pair][0] = inverse_matrix[speaker][0];
- _matrices[pair][1] = inverse_matrix[speaker][1];
- _matrices[pair][2] = inverse_matrix[speaker][2];
- _matrices[pair][3] = inverse_matrix[speaker][3];
+ _matrices[pair][0] = inverse_matrix[(speaker*4)+0];
+ _matrices[pair][1] = inverse_matrix[(speaker*4)+1];
+ _matrices[pair][2] = inverse_matrix[(speaker*4)+2];
+ _matrices[pair][3] = inverse_matrix[(speaker*4)+3];
_speaker_tuples[pair][0] = sorted_speakers[speaker];
_speaker_tuples[pair][1] = sorted_speakers[speaker+1];
pair++;
}
}
-
+
if (exists[n_speakers-1]) {
- _matrices[pair][0] = inverse_matrix[speaker][0];
- _matrices[pair][1] = inverse_matrix[speaker][1];
- _matrices[pair][2] = inverse_matrix[speaker][2];
- _matrices[pair][3] = inverse_matrix[speaker][3];
+ _matrices[pair][0] = inverse_matrix[(speaker*4)+0];
+ _matrices[pair][1] = inverse_matrix[(speaker*4)+1];
+ _matrices[pair][2] = inverse_matrix[(speaker*4)+2];
+ _matrices[pair][3] = inverse_matrix[(speaker*4)+3];
_speaker_tuples[pair][0] = sorted_speakers[n_speakers-1];
_speaker_tuples[pair][1] = sorted_speakers[0];
}
}
-void
+void
VBAPSpeakers::sort_2D_lss (int* sorted_speakers)
{
vector<Speaker> tmp = _speakers;
vector<Speaker>::iterator s;
azimuth_sorter sorter;
- int n;
+ unsigned int n;
sort (tmp.begin(), tmp.end(), sorter);
for (n = 0, s = tmp.begin(); s != tmp.end(); ++s, ++n) {
sorted_speakers[n] = (*s).id;
}
+ assert(n == _speakers.size ());
}
-int
+int
VBAPSpeakers::calc_2D_inv_tmatrix (double azi1, double azi2, double* inverse_matrix)
{
double x1,x2,x3,x4;
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;
projects / ardour.git / blobdiff