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()
)

2 files changed, 113 insertions, 92 deletions

diff --git a/libs/panners/vbap/vbap_speakers.cc b/libs/panners/vbap/vbap_speakers.cc
index b84698bbd9..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
@@ -90,18 +90,29 @@ 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];
@@ -117,8 +128,8 @@ VBAPSpeakers::choose_speaker_triplets(struct ls_triplet_chain **ls_triplets)
 
 	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;
@@ -132,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;
@@ -167,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;
 						}
@@ -186,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) {
@@ -209,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 */
@@ -227,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;
@@ -242,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) {
@@ -265,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*/
@@ -273,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;
@@ -293,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;
-	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 {
@@ -345,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.
 	*/
@@ -376,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;
@@ -397,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)  &&
@@ -416,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;
@@ -430,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) {
@@ -438,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;
 
@@ -450,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;
@@ -470,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];
@@ -487,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++;
 
@@ -498,7 +519,7 @@ VBAPSpeakers::calculate_3x3_matrixes(struct ls_triplet_chain *ls_triplets)
 	}
 }
 
-void 
+void
 VBAPSpeakers::choose_speaker_pairs (){
 
 	/* selects the loudspeaker pairs, calculates the inversion
@@ -508,7 +529,7 @@ VBAPSpeakers::choose_speaker_pairs (){
 	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]; 
+	double inverse_matrix[n_speakers][4];
 	int expected_pairs = 0;
 	int pair;
 	int speaker;
@@ -524,29 +545,29 @@ VBAPSpeakers::choose_speaker_pairs (){
 
 	/* sort loudspeakers according their aximuth angle */
 	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, 
+			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) { 
+		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;
@@ -572,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];
@@ -584,7 +605,7 @@ VBAPSpeakers::choose_speaker_pairs (){
 	}
 }
 
-void 
+void
 VBAPSpeakers::sort_2D_lss (int* sorted_speakers)
 {
 	vector<Speaker> tmp = _speakers;
@@ -599,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;
@@ -612,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<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);