1 /*
2     vbap.c:
3 
4     Copyright (C) 2000 Ville Pulkki
5 
6     This file is part of Csound.
7 
8     The Csound Library is free software; you can redistribute it
9     and/or modify it under the terms of the GNU Lesser General Public
10     License as published by the Free Software Foundation; either
11     version 2.1 of the License, or (at your option) any later version.
12 
13     Csound is distributed in the hope that it will be useful,
14     but WITHOUT ANY WARRANTY; without even the implied warranty of
15     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16     GNU Lesser General Public License for more details.
17 
18     You should have received a copy of the GNU Lesser General Public
19     License along with Csound; if not, write to the Free Software
20     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
21     02110-1301 USA
22 */
23 
24 /* vbap.c
25 
26 assisting functions for VBAP
27 functions for loudspeaker table initialization
28 Re-written to take flexible number of outputs by JPff 2012 */
29 
30 
31 #include "csoundCore.h"
32 #include "interlocks.h"
33 #include "vbap.h"
34 #include <math.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include "interlocks.h"
38 
39 #define MATSIZE (4)
40 #define ATORAD  (TWOPI_F / FL(360.0))
41 
42 /* static void choose_ls_triplets(CSOUND *csound, ls *lss, */
43 /*                                ls_triplet_chain **ls_triplets, */
44 /*                                int32_t ls_amount, int32_t channels); */
45 static int32_t any_ls_inside_triplet(int32_t, int32_t, int32_t, ls[], int32_t);
46 static void add_ldsp_triplet(CSOUND *csound, int32_t i, int32_t j, int32_t k,
47                              ls_triplet_chain **ls_triplets,
48                              ls *lss);
49 static void calculate_3x3_matrixes(CSOUND *csound,
50                                    ls_triplet_chain *ls_triplets,
51                                    ls lss[], int32_t ls_amount, int32_t ind);
52 static void choose_ls_tuplets(CSOUND *csound, ls lss[],
53                               ls_triplet_chain **ls_triplets,
54                               int32_t ls_amount, int32_t ind);
55 static void sort_2D_lss(ls lss[], int32_t sorted_lss[],
56                         int32_t ls_amount);
57 
vec_prod(CART_VEC v1,CART_VEC v2)58 static inline MYFLT vec_prod(CART_VEC v1, CART_VEC v2)
59 {
60     return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
61 }
62 
vec_length(CART_VEC v1)63 static inline MYFLT vec_length(CART_VEC v1)
64 {
65     return SQRT(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z);
66 }
67 
create_ls_table(CSOUND * csound,size_t cnt,int32_t ind)68 static MYFLT *create_ls_table(CSOUND *csound, size_t cnt, int32_t ind)
69 {
70     char name[24];
71     snprintf(name, 24, "vbap_ls_table_%d", ind);
72     csound->DestroyGlobalVariable(csound, name);
73     if (UNLIKELY(csound->CreateGlobalVariable(csound, name,
74                                               cnt * sizeof(MYFLT)) != 0)) {
75       csound->ErrorMsg(csound, Str("vbap: error allocating loudspeaker table"));
76       return NULL;
77     }
78     return (MYFLT*) (csound->QueryGlobalVariableNoCheck(csound, name));
79 }
80 
calc_vbap_gns(int32_t ls_set_am,int dim,LS_SET * sets,MYFLT * gains,int32_t ls_amount,CART_VEC cart_dir)81 void calc_vbap_gns(int32_t ls_set_am, int dim, LS_SET *sets,
82                    MYFLT *gains, int32_t ls_amount,
83                    CART_VEC cart_dir)
84      /* Selects a vector base of a virtual source.
85         Calculates gain factors in that base. */
86 {
87     int32_t i,j,k, tmp2;
88     MYFLT vec[3], tmp;
89     /* direction of the virtual source in cartesian coordinates*/
90     vec[0] = cart_dir.x;
91     vec[1] = cart_dir.y;
92     vec[2] = cart_dir.z;
93 
94 
95     for (i=0; i< ls_set_am; i++) {
96       sets[i].set_gains[0] = FL(0.0);
97       sets[i].set_gains[1] = FL(0.0);
98       sets[i].set_gains[2] = FL(0.0);
99       sets[i].smallest_wt  = FL(1000.0);
100       sets[i].neg_g_am = 0;
101   }
102 
103     for (i=0; i< ls_set_am; i++) {
104       for (j=0; j< dim; j++) {
105         for (k=0; k< dim; k++) {
106           sets[i].set_gains[j] += vec[k] * sets[i].ls_mx[((dim * j )+ k)];
107         }
108         if (sets[i].smallest_wt > sets[i].set_gains[j])
109           sets[i].smallest_wt = sets[i].set_gains[j];
110         if (sets[i].set_gains[j] < -FL(0.05))
111           sets[i].neg_g_am++;
112       }
113     }
114 
115 
116 
117     j=0;
118     tmp = sets[0].smallest_wt;
119     tmp2=sets[0].neg_g_am;
120     for (i=1; i< ls_set_am; i++) {
121       if (sets[i].neg_g_am < tmp2) {
122         tmp = sets[i].smallest_wt;
123         tmp2=sets[i].neg_g_am;
124         j=i;
125       }
126       else if (sets[i].neg_g_am == tmp2) {
127         if (sets[i].smallest_wt > tmp) {
128           tmp = sets[i].smallest_wt;
129           tmp2=sets[i].neg_g_am;
130           j=i;
131         }
132       }
133     }
134 
135 
136     if (sets[j].set_gains[0]<=FL(0.0) &&
137         sets[j].set_gains[1]<=FL(0.0) &&
138         sets[j].set_gains[2]<=FL(0.0)) {
139       sets[j].set_gains[0] = FL(1.0);
140       sets[j].set_gains[1] = FL(1.0);
141       sets[j].set_gains[2] = FL(1.0);
142     }
143 
144     memset(gains, 0, ls_amount*sizeof(MYFLT));
145 
146     gains[sets[j].ls_nos[0]-1] = sets[j].set_gains[0];
147     gains[sets[j].ls_nos[1]-1] = sets[j].set_gains[1];
148     if (dim==3) gains[sets[j].ls_nos[2]-1] = sets[j].set_gains[2];
149 
150     for (i=0;i<ls_amount;i++) {
151       if (gains[i]<LOWEST_ACCEPTABLE_WT)
152         gains[i]=FL(0.0);
153     }
154 }
155 
scale_angles(ANG_VEC * avec)156 void scale_angles(ANG_VEC *avec)
157      /* -180 < azi < 180
158         -90 < ele < 90 */
159 {
160     while (avec->azi > FL(180.0))
161       avec->azi -= FL(360.0);
162     while (avec->azi < -FL(180.0))
163       avec->azi += FL(360.0);
164     if (avec->ele > FL(90.0))
165       avec->ele = FL(90.0);
166     if (avec->ele < -FL(90.0))
167       avec->ele = -FL(90.0);
168 }
169 
normalize_wts(OUT_WTS * wts)170 void normalize_wts(OUT_WTS *wts)
171      /* performs equal-power normalization to gain factors*/
172 {
173     double tmp;
174     MYFLT tmp1;
175     if (wts->wt1 < 0) wts->wt1 = FL(0.0);
176     if (wts->wt2 < 0) wts->wt2 = FL(0.0);
177     if (wts->wt3 < 0) wts->wt3 = FL(0.0);
178 
179     tmp  = (double)wts->wt1 * wts->wt1;
180     tmp += (double)wts->wt2 * wts->wt2;
181     tmp += (double)wts->wt3 * wts->wt3;
182 
183     tmp = sqrt(tmp);
184     tmp1 = (MYFLT)(1.0 / tmp);
185     wts->wt1 *= tmp1;
186     wts->wt2 *= tmp1;
187     wts->wt3 *= tmp1;
188 }
189 
angle_to_cart(ANG_VEC avec,CART_VEC * cvec)190 void angle_to_cart(ANG_VEC avec, CART_VEC *cvec)
191      /* conversion */
192 {
193     /* length unattended */
194     //MYFLT atorad = (TWOPI_F / FL(360.0));
195     cvec->x = (MYFLT) (cos((double) (avec.azi * ATORAD)) *
196                        cos((double) (avec.ele * ATORAD)));
197     cvec->y = (MYFLT) (sin((double) (avec.azi * ATORAD)) *
198                        cos((double) (avec.ele * ATORAD)));
199     cvec->z = (MYFLT) (sin((double) (avec.ele * ATORAD)));
200 }
201 
cart_to_angle(CART_VEC cvec,ANG_VEC * avec)202 void cart_to_angle(CART_VEC cvec, ANG_VEC *avec)
203      /* conversion */
204 {
205     MYFLT tmp, tmp2, tmp3, tmp4;
206     //MYFLT atorad = (TWOPI_F / FL(360.0));
207 
208     tmp3 = SQRT(FL(1.0) - cvec.z*cvec.z);
209     if (FABS(tmp3) > FL(0.001)) {
210       tmp4 = (cvec.x / tmp3);
211       if (tmp4 > FL(1.0)) tmp4 = FL(1.0);
212       if (tmp4 < -FL(1.0)) tmp4 = -FL(1.0);
213       tmp = ACOS(tmp4 );
214     }
215     else {
216       tmp = FL(10000.0);
217     }
218     if (FABS(cvec.y) <= FL(0.001))
219       tmp2 = FL(1.0);
220     else
221       tmp2 = cvec.y / FABS(cvec.y);
222     tmp *= tmp2;
223     if (FABS(tmp) <= PI_F) {
224       avec->azi =  tmp;
225       avec->azi /= ATORAD;
226     }
227     avec->ele = ASIN(cvec.z);
228     avec->length = SQRT(cvec.x * cvec.x + cvec.y * cvec.y + cvec.z * cvec.z);
229     avec->ele /= ATORAD;
230 }
231 
angle_to_cart_II(ANG_VEC * from,CART_VEC * to)232 void angle_to_cart_II(ANG_VEC *from, CART_VEC *to)
233      /* conversion, double*/
234 {
235     MYFLT ang2rad = TWOPI_F / FL(360.0);
236     to->x= COS(from->azi * ang2rad) * COS(from->ele * ang2rad);
237     to->y= SIN(from->azi * ang2rad) * COS(from->ele * ang2rad);
238     to->z= SIN(from->ele * ang2rad);
239 }
240 
vol_p_side_lgth(int32_t i,int32_t j,int32_t k,ls lss[])241 MYFLT vol_p_side_lgth(int32_t i, int32_t j,int32_t k, ls  lss[] )
242 {
243   /* calculate volume of the parallelepiped defined by the loudspeaker
244      direction vectors and divide it with total length of the triangle sides.
245      This is used when removing too narrow triangles. */
246 
247     MYFLT volper, lgth;
248     CART_VEC xprod;
249     cross_prod(lss[i].coords, lss[j].coords, &xprod);
250     volper = FABS(vec_prod(xprod, lss[k].coords));
251     lgth =    FABS(vec_angle(lss[i].coords,lss[j].coords))
252             + FABS(vec_angle(lss[i].coords,lss[k].coords))
253             + FABS(vec_angle(lss[j].coords,lss[k].coords));
254     if (LIKELY(lgth>FL(0.00001)))
255       return volper / lgth;
256     else
257       return FL(0.0);
258 }
259 
choose_ls_triplets(CSOUND * csound,ls * lss,struct ls_triplet_chain ** ls_triplets,int32_t ls_amount)260 static void choose_ls_triplets(CSOUND *csound, ls *lss,
261                                struct ls_triplet_chain **ls_triplets,
262                                int32_t ls_amount)
263   /* Selects the loudspeaker triplets, and
264      calculates the inversion matrices for each selected triplet.
265      A line (connection) is drawn between each loudspeaker. The lines
266      denote the sides of the triangles. The triangles should not be
267      intersecting. All crossing connections are searched and the
268      longer connection is erased. This yields non-intesecting triangles,
269      which can be used in panning.*/
270 {
271     int32_t i, j, k, l, table_size;
272     int32_t *connections;
273 /*  int32_t *i_ptr; */
274     MYFLT *distance_table;
275     int32_t *distance_table_i;
276     int32_t *distance_table_j;
277     MYFLT distance;
278     struct ls_triplet_chain *trip_ptr, *prev, *tmp_ptr;
279 
280     if (UNLIKELY(ls_amount == 0)) {
281       csound->ErrorMsg(csound, Str("Number of loudspeakers is zero\nExiting"));
282       return;
283     }
284 
285     connections = csound->Calloc(csound, ls_amount * ls_amount * sizeof(int32_t));
286     distance_table =
287       csound->Calloc(csound, ((ls_amount * (ls_amount - 1)) / 2)* sizeof(MYFLT));
288     distance_table_i =
289       csound->Calloc(csound, ((ls_amount * (ls_amount - 1)) / 2)* sizeof(int32_t));
290     distance_table_j =
291       csound->Calloc(csound, ((ls_amount * (ls_amount - 1)) / 2)* sizeof(int32_t));
292 
293 /*  i_ptr = (int32_t *) connections; */
294 /*  for (i=0;i< ((CHANNELS) * (CHANNELS )); i++) */
295 /*    *(i_ptr++) = 0; */
296 
297     for (i=0;i<ls_amount;i++)
298     for (j=i+1;j<ls_amount;j++)
299       for (k=j+1;k<ls_amount;k++) {
300         if (vol_p_side_lgth(i,j, k, lss) > MIN_VOL_P_SIDE_LGTH) {
301           connections[i+ls_amount*j]=1;
302           connections[j+ls_amount*i]=1;
303           connections[i+ls_amount*k]=1;
304           connections[k+ls_amount*i]=1;
305           connections[j+ls_amount*k]=1;
306           connections[k+ls_amount*j]=1;
307           add_ldsp_triplet(csound, i, j, k, ls_triplets, lss);
308         }
309       }
310 
311   /*calculate distancies between all lss and sorting them*/
312     table_size =(((ls_amount - 1) * (ls_amount)) / 2);
313     for (i=0;i<table_size; i++)
314       distance_table[i] = FL(100000.0);
315 
316     for (i=0;i<ls_amount;i++) {
317       for (j=(i+1);j<ls_amount; j++) {
318         if (connections[i+ls_amount*j] == 1) {
319           distance = FABS(vec_angle(lss[i].coords,lss[j].coords));
320           k=0;
321 
322           while (distance_table[k] < distance)
323             k++;
324           for (l=(table_size - 1);l > k;l--) {
325             distance_table[l] = distance_table[l-1];
326             distance_table_i[l] = distance_table_i[l-1];
327             distance_table_j[l] = distance_table_j[l-1];
328           }
329           distance_table[k] = distance;
330           distance_table_i[k] = i;
331           distance_table_j[k] = j;
332         }
333         else
334           table_size--;
335       }
336     }
337 
338     /* disconnecting connections which are crossing shorter ones,
339        starting from shortest one and removing all that cross it,
340        and proceeding to next shortest */
341     for (i=0; i<(table_size); i++) {
342       int32_t fst_ls = distance_table_i[i];
343       int32_t sec_ls = distance_table_j[i];
344       if (connections[fst_ls+ls_amount*sec_ls] == 1)
345         for (j=0; j<ls_amount; j++)
346           for (k=j+1; k<ls_amount; k++)
347             if ( (j!=fst_ls) && (k != sec_ls) && (k!=fst_ls) && (j != sec_ls))
348               if (lines_intersect(fst_ls, sec_ls, j,k,lss) == 1) {
349                 connections[j+ls_amount*k] = 0;
350                 connections[k+ls_amount*j] = 0;
351               }
352     }
353 
354     /* remove triangles which had crossing sides
355        with smaller triangles or include loudspeakers*/
356 
357     trip_ptr = *ls_triplets;
358     prev = NULL;
359     while (trip_ptr != NULL) {
360       i = trip_ptr->ls_nos[0];
361       j = trip_ptr->ls_nos[1];
362       k = trip_ptr->ls_nos[2];
363       if (connections[i+ls_amount*j] == 0 ||
364           connections[i+ls_amount*k] == 0 ||
365           connections[j+ls_amount*k] == 0 ||
366           any_ls_inside_triplet(i,j,k,lss,ls_amount) == 1) {
367         if (prev != NULL) {
368           prev->next = trip_ptr->next;
369           tmp_ptr = trip_ptr;
370           trip_ptr = trip_ptr->next;
371           csound->Free(csound, tmp_ptr);
372         }
373         else {
374           *ls_triplets = trip_ptr->next;
375           tmp_ptr = trip_ptr;
376           trip_ptr = trip_ptr->next;
377           csound->Free(csound, tmp_ptr);
378         }
379       }
380       else {
381         prev = trip_ptr;
382         trip_ptr = trip_ptr->next;
383       }
384     }
385     csound->Free(csound,connections);
386     csound->Free(csound,distance_table);
387     csound->Free(csound,distance_table_i);
388     csound->Free(csound,distance_table_j);
389 }
390 
391 /* returns 1 if there is loudspeaker(s) inside given ls triplet */
392 
any_ls_inside_triplet(int32_t a,int32_t b,int32_t c,ls lss[],int32_t ls_amount)393 static int32_t any_ls_inside_triplet(int32_t a, int32_t b, int32_t c, ls lss[],
394                                  int32_t ls_amount)
395 {
396     MYFLT invdet;
397     CART_VEC *lp1, *lp2, *lp3;
398     MYFLT invmx[9];
399     int32_t i,j;
400     MYFLT tmp;
401     int32_t any_ls_inside, this_inside;
402 
403     lp1 =  &(lss[a].coords);
404     lp2 =  &(lss[b].coords);
405     lp3 =  &(lss[c].coords);
406 
407     /* matrix inversion */
408     invdet = FL(1.0) / (  lp1->x * ((lp2->y * lp3->z) - (lp2->z * lp3->y))
409                        - lp1->y * ((lp2->x * lp3->z) - (lp2->z * lp3->x))
410                        + lp1->z * ((lp2->x * lp3->y) - (lp2->y * lp3->x)));
411 
412     invmx[0] = ((lp2->y * lp3->z) - (lp2->z * lp3->y)) * invdet;
413     invmx[3] = ((lp1->y * lp3->z) - (lp1->z * lp3->y)) * -invdet;
414     invmx[6] = ((lp1->y * lp2->z) - (lp1->z * lp2->y)) * invdet;
415     invmx[1] = ((lp2->x * lp3->z) - (lp2->z * lp3->x)) * -invdet;
416     invmx[4] = ((lp1->x * lp3->z) - (lp1->z * lp3->x)) * invdet;
417     invmx[7] = ((lp1->x * lp2->z) - (lp1->z * lp2->x)) * -invdet;
418     invmx[2] = ((lp2->x * lp3->y) - (lp2->y * lp3->x)) * invdet;
419     invmx[5] = ((lp1->x * lp3->y) - (lp1->y * lp3->x)) * -invdet;
420     invmx[8] = ((lp1->x * lp2->y) - (lp1->y * lp2->x)) * invdet;
421 
422     any_ls_inside = 0;
423     for (i=0; i< ls_amount; i++) {
424       if (i != a && i!=b && i != c) {
425         this_inside = 1;
426         for (j=0; j< 3; j++) {
427           tmp = lss[i].coords.x * invmx[0 + j*3];
428           tmp += lss[i].coords.y * invmx[1 + j*3];
429           tmp += lss[i].coords.z * invmx[2 + j*3];
430           if (tmp < -FL(0.001))
431             this_inside = 0;
432         }
433         if (this_inside == 1)
434           any_ls_inside=1;
435       }
436     }
437     return any_ls_inside;
438 }
439 
add_ldsp_triplet(CSOUND * csound,int32_t i,int32_t j,int32_t k,struct ls_triplet_chain ** ls_triplets,ls lss[])440 static void add_ldsp_triplet(CSOUND *csound, int32_t i, int32_t j, int32_t k,
441                              struct ls_triplet_chain **ls_triplets,
442                              ls lss[])
443 {
444     IGN(lss);
445     struct ls_triplet_chain *ls_ptr, *prev;
446     ls_ptr = *ls_triplets;
447     prev = NULL;
448 
449   /*printf("Adding triangle %d %d %d %x... \n",i,j,k,ls_ptr);*/
450     while (ls_ptr != NULL) {
451       /*printf("ls_ptr %x %x\n",ls_ptr,ls_ptr->next);*/
452       prev = ls_ptr;
453       ls_ptr = ls_ptr->next;
454     }
455     ls_ptr = (struct ls_triplet_chain*)
456       csound->Malloc(csound, sizeof(struct ls_triplet_chain));
457     if (prev == NULL)
458       *ls_triplets = ls_ptr;
459     else
460       prev->next = ls_ptr;
461     ls_ptr->next = NULL;
462     ls_ptr->ls_nos[0] = i;
463     ls_ptr->ls_nos[1] = j;
464     ls_ptr->ls_nos[2] = k;
465     /*printf("added.\n");*/
466 }
467 
angle_in_base(CART_VEC vb1,CART_VEC vb2,CART_VEC vec)468 MYFLT angle_in_base(CART_VEC vb1,CART_VEC vb2,CART_VEC vec)
469 {
470     MYFLT tmp1,tmp2;
471     tmp1 = vec_prod(vec,vb2);
472     if (FABS(tmp1) <= FL(0.001))
473       tmp2 = FL(1.0);
474     else
475       tmp2 = tmp1 / FABS(tmp1);
476     return (vec_angle(vb1,vec) * tmp2);
477 }
478 
vec_angle(CART_VEC v1,CART_VEC v2)479 MYFLT vec_angle(CART_VEC v1, CART_VEC v2)
480 {
481     MYFLT inner= ((v1.x*v2.x + v1.y*v2.y + v1.z*v2.z)/
482                   (vec_length(v1) * vec_length(v2)));
483     if (inner > FL(1.0))
484       inner= FL(1.0);
485     if (inner < -FL(1.0))
486       inner = -FL(1.0);
487     return ACOS(inner);
488 }
489 
vec_mean(CART_VEC v1,CART_VEC v2,CART_VEC * v3)490 void vec_mean(CART_VEC v1, CART_VEC v2, CART_VEC *v3)
491 {
492     v3->x=(v1.x+v2.x)*FL(0.5);
493     v3->y=(v1.y+v2.y)*FL(0.5);
494     v3->z=(v1.z+v2.z)*FL(0.5);
495 }
496 
cross_prod(CART_VEC v1,CART_VEC v2,CART_VEC * res)497 void cross_prod(CART_VEC v1,CART_VEC v2,
498                 CART_VEC *res)
499 {
500     MYFLT length;
501     res->x = (v1.y * v2.z ) - (v1.z * v2.y);
502     res->y = (v1.z * v2.x ) - (v1.x * v2.z);
503     res->z = (v1.x * v2.y ) - (v1.y * v2.x);
504 
505     length= vec_length(*res);
506     res->x /= length;
507     res->y /= length;
508     res->z /= length;
509 }
510 
vec_print(CSOUND * csound,CART_VEC v)511 void vec_print(CSOUND *csound, CART_VEC v)
512 {
513     csound->Message(csound, "vec_print %f %f %f\n", v.x, v.y,v.z);
514 
515 }
516 
lines_intersect(int32_t i,int32_t j,int32_t k,int32_t l,ls lss[])517 int32_t lines_intersect(int32_t i,int32_t j,int32_t k,int32_t l,ls  lss[])
518   /* checks if two lines intersect on 3D sphere
519      see theory in paper Pulkki, V. Lokki, T. "Creating Auditory Displays
520      with Multiple Loudspeakers Using VBAP: A Case Study with
521      DIVA Project" in International Conference on
522      Auditory Displays -98. E-mail Ville.Pulkki@hut.fi
523      if you want to have that paper.*/
524 {
525     CART_VEC v1;
526     CART_VEC v2;
527     CART_VEC v3, neg_v3;
528     MYFLT dist_ij,dist_kl,dist_iv3,dist_jv3,dist_inv3,dist_jnv3;
529     MYFLT dist_kv3,dist_lv3,dist_knv3,dist_lnv3;
530 
531     cross_prod(lss[i].coords,lss[j].coords,&v1);
532     cross_prod(lss[k].coords,lss[l].coords,&v2);
533     cross_prod(v1,v2,&v3);
534 
535     neg_v3.x= FL(0.0) - v3.x;
536     neg_v3.y= FL(0.0) - v3.y;
537     neg_v3.z= FL(0.0) - v3.z;
538 
539     dist_ij = (vec_angle(lss[i].coords,lss[j].coords));
540     dist_kl = (vec_angle(lss[k].coords,lss[l].coords));
541     dist_iv3 = (vec_angle(lss[i].coords,v3));
542     dist_jv3 = (vec_angle(v3,lss[j].coords));
543     dist_inv3 = (vec_angle(lss[i].coords,neg_v3));
544     dist_jnv3 = (vec_angle(neg_v3,lss[j].coords));
545     dist_kv3 = (vec_angle(lss[k].coords,v3));
546     dist_lv3 = (vec_angle(v3,lss[l].coords));
547     dist_knv3 = (vec_angle(lss[k].coords,neg_v3));
548     dist_lnv3 = (vec_angle(neg_v3,lss[l].coords));
549 
550     /* if one of loudspeakers is close to crossing point, don't do anything*/
551     if (FABS(dist_iv3) <= FL(0.01) || FABS(dist_jv3) <= FL(0.01) ||
552         FABS(dist_kv3) <= FL(0.01) || FABS(dist_lv3) <= FL(0.01) ||
553         FABS(dist_inv3) <= FL(0.01) || FABS(dist_jnv3) <= FL(0.01) ||
554         FABS(dist_knv3) <= FL(0.01) || FABS(dist_lnv3) <= FL(0.01) )
555       return(0);
556 
557     if (((FABS(dist_ij - (dist_iv3 + dist_jv3)) <= FL(0.01) ) &&
558          (FABS(dist_kl - (dist_kv3 + dist_lv3))  <= FL(0.01))) ||
559         ((FABS(dist_ij - (dist_inv3 + dist_jnv3)) <= FL(0.01))  &&
560          (FABS(dist_kl - (dist_knv3 + dist_lnv3)) <= FL(0.01) ))) {
561       return (1);
562     }
563     else {
564       return (0);
565     }
566 }
567 
vbap_ls_init_sr(CSOUND * csound,int32_t dim,int32_t count,MYFLT ** f,int32_t layout)568 static inline int32_t vbap_ls_init_sr (CSOUND *csound, int32_t dim, int32_t count,
569                             MYFLT **f, int32_t layout)
570      /* Inits the loudspeaker data. Calls choose_ls_tuplets or _triplets
571         according to current dimension. The inversion matrices are
572         stored in transposed form to ease calculation at run time.*/
573 {
574     struct ls_triplet_chain *ls_triplets = NULL;
575     ls *lss = malloc(sizeof(ls)*count);
576 
577     ANG_VEC a_vector;
578     CART_VEC c_vector;
579     int32_t i=0,j;
580 
581     //dim = (int32_t) *p->dim;
582     csound->Message(csound, "dim : %d\n",dim);
583     if (UNLIKELY(!((dim==2) || (dim == 3)))) {
584       free(lss);
585       csound->ErrorMsg(csound,
586                        Str("Error in loudspeaker dimension. %d not permitted"),
587                        dim);
588       return NOTOK;
589     }
590     //count = (int32_t) *p->ls_amount;
591     for (j=1;j<=count;j++) {
592       if (dim == 3) {
593         a_vector.azi= (MYFLT) *f[2*j-2];
594         a_vector.ele= (MYFLT) *f[2*j-1];
595       }
596       else if (dim == 2) {
597         a_vector.azi= (MYFLT) *f[j-1];
598         a_vector.ele=FL(0.0);
599       }
600       angle_to_cart_II(&a_vector,&c_vector);
601       lss[i].coords.x = c_vector.x;
602       lss[i].coords.y = c_vector.y;
603       lss[i].coords.z = c_vector.z;
604       lss[i].angles.azi = a_vector.azi;
605       lss[i].angles.ele = a_vector.ele;
606       lss[i].angles.length = FL(1.0);
607       /* printf("**** lss[%d]: (%g %g %g) %g %g\n", i, lss[i].coords.x, */
608       /*        lss[i].coords.y, lss[i].coords.z, a_vector.azi, a_vector.ele); */
609       i++;
610     }
611     //ls_amount = (int32_t)*p->ls_amount;
612     if (UNLIKELY(count < dim)) {
613       free(lss);
614       csound->ErrorMsg(csound, Str("Too few loudspeakers"));
615       return NOTOK;
616     }
617 
618     if (dim == 3) {
619       choose_ls_triplets(csound, lss, &ls_triplets, count);
620       calculate_3x3_matrixes(csound, ls_triplets, lss, count, layout);
621     }
622     else if (dim ==2) {
623       choose_ls_tuplets(csound, lss, &ls_triplets, count, layout);
624     }
625     free(lss);
626     return OK;
627 }
628 
vbap_ls_init(CSOUND * csound,VBAP_LS_INIT * p)629 int32_t vbap_ls_init (CSOUND *csound, VBAP_LS_INIT *p)
630 {
631     int32_t dim = (int32_t) *p->dim;
632     MYFLT  layout = (*p->dim-dim)*100;
633     return vbap_ls_init_sr(csound, dim, (int32_t) *p->ls_amount,
634                            p->f, round(layout));
635 }
636 
vbap_ls_inita(CSOUND * csound,VBAP_LS_INITA * p)637 int32_t vbap_ls_inita (CSOUND *csound, VBAP_LS_INITA *p)
638 {
639     int32_t dim = (int32_t) *p->dim;
640     MYFLT  layout = (*p->dim-dim)*100;
641     int32_t i, n = (int32_t)*p->ls_amount;
642     /* if (n>CHANNELS) */
643     /*   return csound->InitError(csound, Str("Too many speakers (%n)\n"), n); */
644     if (UNLIKELY(n>p->a->sizes[0]))
645       return csound->InitError(csound, Str("Too little data speakers (%d)\n"),
646                               n>p->a->sizes[0]);
647     MYFLT  **f = csound->Malloc(csound, 2*sizeof(MYFLT*)*n);
648     // Transfer values to pointers
649     for (i=0; i<2*n; i++) f[i] = &(p->a->data[i]);
650     n = vbap_ls_init_sr(csound, dim, n, f, round(layout));
651     csound->Free(csound, f);
652     return n;
653 }
654 
calculate_3x3_matrixes(CSOUND * csound,struct ls_triplet_chain * ls_triplets,ls lss[],int32_t ls_amount,int32_t ind)655 static void calculate_3x3_matrixes(CSOUND *csound,
656                                    struct ls_triplet_chain *ls_triplets,
657                                    ls lss[], int32_t ls_amount, int32_t ind)
658      /* Calculates the inverse matrices for 3D */
659 {
660     MYFLT invdet;
661     CART_VEC *lp1, *lp2, *lp3;
662     MYFLT *ls_table, *invmx;
663     MYFLT *ptr;
664     struct ls_triplet_chain *tr_ptr = ls_triplets;
665     int32_t triplet_amount = 0, i,j,k;
666 
667     if (UNLIKELY(tr_ptr == NULL)) {
668       csound->ErrorMsg(csound, Str("Not valid 3-D configuration"));
669       return;
670     }
671 
672     /* counting triplet amount */
673     while (tr_ptr != NULL) {
674       triplet_amount++;
675       tr_ptr = tr_ptr->next;
676     }
677 
678     /* calculations and data storage to a global array */
679     ls_table = create_ls_table(csound, triplet_amount * 12 + 3, ind);
680     ls_table[0] = FL(3.0);  /* dimension */
681     ls_table[1] = (MYFLT) ls_amount;
682     ls_table[2] = (MYFLT) triplet_amount;
683     tr_ptr = ls_triplets;
684     ptr = (MYFLT *) &(ls_table[3]);
685     while (tr_ptr != NULL) {
686       lp1 =  &(lss[tr_ptr->ls_nos[0]].coords);
687       lp2 =  &(lss[tr_ptr->ls_nos[1]].coords);
688       lp3 =  &(lss[tr_ptr->ls_nos[2]].coords);
689 
690       /* matrix inversion */
691       invmx = tr_ptr->inv_mx;
692       invdet = FL(1.0) / (  lp1->x * ((lp2->y * lp3->z) - (lp2->z * lp3->y))
693                          - lp1->y * ((lp2->x * lp3->z) - (lp2->z * lp3->x))
694                          + lp1->z * ((lp2->x * lp3->y) - (lp2->y * lp3->x)));
695 
696       invmx[0] = ((lp2->y * lp3->z) - (lp2->z * lp3->y)) * invdet;
697       invmx[3] = ((lp1->y * lp3->z) - (lp1->z * lp3->y)) * -invdet;
698       invmx[6] = ((lp1->y * lp2->z) - (lp1->z * lp2->y)) * invdet;
699       invmx[1] = ((lp2->x * lp3->z) - (lp2->z * lp3->x)) * -invdet;
700       invmx[4] = ((lp1->x * lp3->z) - (lp1->z * lp3->x)) * invdet;
701       invmx[7] = ((lp1->x * lp2->z) - (lp1->z * lp2->x)) * -invdet;
702       invmx[2] = ((lp2->x * lp3->y) - (lp2->y * lp3->x)) * invdet;
703       invmx[5] = ((lp1->x * lp3->y) - (lp1->y * lp3->x)) * -invdet;
704       invmx[8] = ((lp1->x * lp2->y) - (lp1->y * lp2->x)) * invdet;
705       for (i=0;i<3;i++) {
706         *(ptr++) = (MYFLT) tr_ptr->ls_nos[i]+1;
707       }
708       for (i=0;i<9;i++) {
709         *(ptr++) = (MYFLT) invmx[i];
710       }
711       tr_ptr = tr_ptr->next;
712     }
713 
714     k = 3;
715     csound->Warning(csound, Str("\nConfigured loudspeakers\n"));
716     for (i = 0; i < triplet_amount; i++) {
717       csound->Warning(csound, Str("Triplet %d Loudspeakers: "), i);
718       for (j = 0; j < 3; j++) {
719         csound->Warning(csound, "%d ", (int32_t) ls_table[k++]);
720       }
721       csound->Warning(csound, "\n");
722 
723     /* printf("\nMatrix ");  */
724     /*   for (j = 0; j < 9; j++) { */
725     /*     printf("%f ", ls_table[k]);  */
726     /*     k++; */
727     /*   } */
728     /* printf("\n\n"); */
729     }
730 }
731 
choose_ls_tuplets(CSOUND * csound,ls lss[],ls_triplet_chain ** ls_triplets,int32_t ls_amount,int32_t ind)732 static void choose_ls_tuplets(CSOUND *csound,
733                               ls lss[],
734                               ls_triplet_chain **ls_triplets,
735                               int32_t ls_amount, int32_t ind)
736      /* selects the loudspeaker pairs, calculates the inversion
737         matrices and stores the data to a global array */
738 {
739     IGN(ls_triplets);
740     int32_t i, j, k;
741     int32_t *sorted_lss = (int32_t*)malloc(sizeof(int32_t)*ls_amount);
742     int32_t *exist = (int32_t*)calloc(1,sizeof(int32_t)*ls_amount);
743     int32_t amount = 0;
744     MYFLT *inv_mat = (MYFLT*)malloc(MATSIZE*sizeof(MYFLT)*ls_amount),
745           *ls_table, *ptr;
746     //int32_t ftable_size;
747 
748     /* sort loudspeakers according their aximuth angle */
749     sort_2D_lss(lss,sorted_lss,ls_amount);
750 
751     /* adjacent loudspeakers are the loudspeaker pairs to be used.*/
752     for (i=0;i<(ls_amount-1);i++) {
753       csound->Message(csound, "***%d %d %f %f\n",sorted_lss[i],sorted_lss[i+1],
754                       lss[sorted_lss[i]].angles.azi,
755                       lss[sorted_lss[i+1]].angles.azi);
756       if (LIKELY((lss[sorted_lss[i+1]].angles.azi -
757                   lss[sorted_lss[i]].angles.azi) <= (PI - 0.0175))) {
758         if (LIKELY(calc_2D_inv_tmatrix( lss[sorted_lss[i]].angles.azi,
759                                         lss[sorted_lss[i+1]].angles.azi,
760                                         &inv_mat[MATSIZE*i]) != 0)) {
761           exist[i]=1;
762           amount++;
763         }
764       }
765       else  csound->Warning(csound, Str("Pair of speakers at %f and %f ignored\n"),
766                             lss[sorted_lss[i]].angles.azi*FL(180.0)/PI_F,
767                             lss[sorted_lss[i+1]].angles.azi*FL(180.0)/PI_F);
768     }
769 
770     if (LIKELY(((TWOPI_F - lss[sorted_lss[ls_amount-1]].angles.azi)
771                 +lss[sorted_lss[0]].angles.azi) < (PI - 0.0175))) {
772       //printf("**less than PI type 2- 0.175\n");
773       if (LIKELY(calc_2D_inv_tmatrix(lss[sorted_lss[ls_amount-1]].angles.azi,
774                                      lss[sorted_lss[0]].angles.azi,
775                                      &inv_mat[MATSIZE*(ls_amount-1)]) != 0)) {
776         exist[ls_amount-1]=1;
777         amount++;
778       }
779     }
780     else  csound->Warning(csound, Str("Pair of speakers at %f and %f ignored\n"),
781                           lss[sorted_lss[ls_amount-1]].angles.azi*FL(180.0)/PI_F,
782                           lss[sorted_lss[0]].angles.azi*FL(180.0)/PI_F);
783 
784     if (UNLIKELY(amount==0)) {
785       csound->InitError(csound, Str("insufficient valid speakers"));
786       free(sorted_lss); free(exist); free(inv_mat);
787       return;
788     }
789 
790 #if 0
791     if ( amount*6 + 6 <= 16) ftable_size = 16;
792     else if ( amount*6 + 6 <= 32) ftable_size = 32;
793     else if ( amount*6 + 6 <= 64) ftable_size = 64;
794     else if ( amount*6 + 6 <= 128) ftable_size = 128;
795     else if ( amount*6 + 6 <= 256) ftable_size = 256;
796     else if ( amount*6 + 6 <= 1024) ftable_size = 1024;
797     csound->Message(csound,
798                     "Loudspeaker matrices calculated with configuration : ");
799     for (i=0; i< ls_amount; i++)
800       csound->Message(csound, "%.1f ", lss[i].angles.azi / ATORAD);
801     csound->Message(csound, "\n");
802 #endif
803     ls_table = create_ls_table(csound, amount * 6 + 3 + 100, ind);
804     ls_table[0] = FL(2.0);  /* dimension */
805     ls_table[1] = (MYFLT) ls_amount;
806     ls_table[2] = (MYFLT) amount;
807     ptr = &(ls_table[3]);
808     for (i=0;i<ls_amount - 1;i++) {
809       if (exist[i] == 1) {
810         *(ptr++) = (MYFLT)sorted_lss[i]+1;
811         *(ptr++) = (MYFLT)sorted_lss[i+1]+1;
812         for (j=0;j<MATSIZE;j++) {
813           /*printf("iv_mat i=%d a=%d [%d] %f\n",
814             i, ls_amount, i*MATSIZE+j, inv_mat[i*ls_amount+j]); */
815           *(ptr++) = inv_mat[i*MATSIZE+j];
816         }
817       }
818     }
819     if (exist[ls_amount-1] == 1) {
820       *(ptr++) = (MYFLT)sorted_lss[ls_amount-1]+1;
821       *(ptr++) = (MYFLT)sorted_lss[0]+1;
822       for (j=0;j<MATSIZE;j++) {
823 /*         printf("iv_mat[%d] %f\n", (ls_amount-1)*MATSIZE+j, */
824 /*                inv_mat[(ls_amount-1)*MATSIZE+j]); */
825         *(ptr++) = inv_mat[(ls_amount-1)*MATSIZE+j];
826       }
827     }
828     k=3;
829     csound->Message(csound, Str("\nConfigured loudspeakers\n"));
830     for (i=0; i < amount; i++) {
831       csound->Message(csound, Str("Pair %d Loudspeakers: "), i);
832       for (j=0; j < 2; j++) {
833         csound->Message(csound, "%d ", (int32_t) ls_table[k++]);
834       }
835 
836       csound->Message(csound, "\nMatrix ");
837       for (j=0; j < MATSIZE; j++) {
838         csound->Message(csound, "%f ", ls_table[k]);
839         k++;
840       }
841       csound->Message(csound, "\n\n");
842     }
843     free(sorted_lss); free(exist); free(inv_mat);
844 }
845 
sort_2D_lss(ls lss[],int32_t sorted_lss[],int32_t ls_amount)846 static void sort_2D_lss(ls lss[], int32_t sorted_lss[],
847                         int32_t ls_amount)
848 {
849     int32_t i,j,index=-1;
850     MYFLT tmp, tmp_azi;
851 
852     /* Transforming angles between -180 and 180 */
853     for (i=0; i<ls_amount; i++) {
854       angle_to_cart_II(&lss[i].angles, &lss[i].coords);
855       lss[i].angles.azi = ACOS(lss[i].coords.x);
856       if (FABS(lss[i].coords.y) <= FL(0.001))
857         tmp = FL(1.0);
858       else
859         tmp = lss[i].coords.y / FABS(lss[i].coords.y);
860       lss[i].angles.azi *= tmp;
861       //printf("***tulos %f\n",    lss[i].angles.azi);
862     }
863     for (i=0;i<ls_amount;i++) {
864       tmp = FL(2000.0);
865       for (j=0; j<ls_amount;j++) {
866         if (lss[j].angles.azi <= tmp) {
867           tmp=lss[j].angles.azi;
868           index = j;
869         }
870       }
871       sorted_lss[i]=index;
872       tmp_azi = (lss[index].angles.azi);
873       lss[index].angles.azi = (tmp_azi + FL(4000.0));
874     }
875     for (i=0;i<ls_amount;i++) {
876       tmp_azi = (lss[i].angles.azi);
877       lss[i].angles.azi = (tmp_azi - FL(4000.0));
878     }
879 }
880 
calc_2D_inv_tmatrix(MYFLT azi1,MYFLT azi2,MYFLT inv_mat[MATSIZE])881 int32_t calc_2D_inv_tmatrix(MYFLT azi1,MYFLT azi2, MYFLT inv_mat[MATSIZE])
882 {
883     MYFLT x1,x2,x3,x4; /* x1 x3 */
884     MYFLT det;
885     x1 = COS(azi1 );
886     x2 = SIN(azi1 );
887     x3 = COS(azi2 );
888     x4 = SIN(azi2 );
889     det = (x1 * x4) - ( x3 * x2 );
890     if (FABS(det) <= FL(0.001)) {
891       //printf("unusable*** pair, det %f\n",det);
892       inv_mat[0] = FL(0.0);
893       inv_mat[1] = FL(0.0);
894       inv_mat[2] = FL(0.0);
895       inv_mat[3] = FL(0.0);
896       return 0;
897     }
898     else {
899       //printf("***inv x (%f,%f,%f,%f): det=%f\n", x4, -x3, -x2, x1, det);
900       inv_mat[0] =  (x4 / det);
901       inv_mat[1] =  (-x3 / det);
902       inv_mat[2] =  (-x2 / det);
903       inv_mat[3] =  (x1 / det);
904       return 1;
905     }
906 }
907 
new_spread_dir(CART_VEC * spreaddir,CART_VEC vscartdir,CART_VEC spread_base,MYFLT azi,MYFLT spread)908 void new_spread_dir(CART_VEC *spreaddir, CART_VEC vscartdir,
909                     CART_VEC spread_base, MYFLT azi, MYFLT spread)
910 {
911     MYFLT beta,gamma;
912     MYFLT a,b;
913     MYFLT power;
914     ANG_VEC tmp;
915     gamma = ACOS(vscartdir.x * spread_base.x +
916                  vscartdir.y * spread_base.y +
917                  vscartdir.z * spread_base.z)/PI_F*FL(180.0);
918     if (FABS(gamma) < FL(1.0)) {
919       tmp.azi=azi+FL(90.0);
920       tmp.ele=FL(0.0); tmp.length=FL(1.0);
921       angle_to_cart(tmp, &spread_base);
922       gamma = ACOS(vscartdir.x * spread_base.x +
923                    vscartdir.y * spread_base.y +
924                    vscartdir.z * spread_base.z)/PI_F*FL(180.0);
925     }
926     beta = FL(180.0) - gamma;
927     b=SIN(spread * PI_F / FL(180.0)) /
928       SIN(beta * PI_F / FL(180.0));
929     a=SIN((FL(180.0)- spread - beta) * PI_F / FL(180.0)) /
930       SIN (beta * PI_F / FL(180.0));
931     spreaddir->x = a * vscartdir.x + b * spread_base.x;
932     spreaddir->y = a * vscartdir.y + b * spread_base.y;
933     spreaddir->z = a * vscartdir.z + b * spread_base.z;
934 
935     power=SQRT(spreaddir->x*spreaddir->x +
936                spreaddir->y*spreaddir->y +
937                spreaddir->z*spreaddir->z);
938     spreaddir->x /= power;
939     spreaddir->y /= power;
940     spreaddir->z /= power;
941 }
942 
new_spread_base(CART_VEC spreaddir,CART_VEC vscartdir,MYFLT spread,CART_VEC * spread_base)943 void new_spread_base(CART_VEC spreaddir, CART_VEC vscartdir,
944                      MYFLT spread, CART_VEC *spread_base)
945 {
946     MYFLT d;
947     MYFLT power;
948 
949     d = COS(spread/FL(180.0)*PI_F);
950     spread_base->x = spreaddir.x - d * vscartdir.x;
951     spread_base->y = spreaddir.y - d * vscartdir.y;
952     spread_base->z = spreaddir.z - d * vscartdir.z;
953     power=SQRT(spread_base->x*spread_base->x +
954                spread_base->y*spread_base->y +
955                spread_base->z*spread_base->z);
956     spread_base->x /= power;
957     spread_base->y /= power;
958     spread_base->z /= power;
959 }
960 
961 #define S(x)    sizeof(x)
962 
963 /* static */
964 static OENTRY vbap_localops[] = {
965   { "vbap.a",      S(VBAP),
966     TR, 3,  "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm"
967     "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm",
968     "akOOo",
969     (SUBR) vbap_init,    (SUBR) vbap                   },
970   { "vbap.A",      S(VBAPA), TR, 3,  "a[]",    "akOOo",
971     (SUBR) vbap_init_a,    (SUBR) vbap_a               },
972   { "vbap4",      S(VBAP),
973     TR|_QQ, 3,  "aaaammmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm"
974     "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm",
975     "akOOo", (SUBR) vbap_init, (SUBR) vbap },
976   { "vbap8",      S(VBAP),
977     TR|_QQ, 3,  "aaaaaaaammmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm"
978     "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm",
979     "akOOo",
980     (SUBR) vbap_init,    (SUBR) vbap                   },
981   { "vbap16",      S(VBAP),
982     TR|_QQ, 3,  "aaaaaaaaaaaaaaaammmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm"
983     "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm",
984     "akOOo", (SUBR) vbap_init,    (SUBR) vbap                   },
985   { "vbapg.a",      S(VBAP1),             TR, 3,
986     "zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz"
987     "zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz", "kOOo",
988     (SUBR) vbap1_init,         (SUBR) vbap1                                  },
989   { "vbapg.A",      S(VBAPA1),            TR, 3,
990     "k[]",  "kOOo",
991     (SUBR) vbap1_init_a,         (SUBR) vbap1a                               },
992   { "vbapz",      S(VBAP_ZAK),     ZW|TR, 3,  "",                 "iiakOOo",
993     (SUBR) vbap_zak_init,    (SUBR) vbap_zak         },
994   { "vbaplsinit",S(VBAP_LS_INIT),TR,1, "",
995     "ii"
996     "oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo"
997     "oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo",
998     (SUBR) vbap_ls_init, (SUBR) NULL, (SUBR) NULL, (SUBR) NULL         },
999   { "vbaplsinit",S(VBAP_LS_INIT),TR,1, "", "iii[]",
1000     (SUBR) vbap_ls_inita, (SUBR) NULL, (SUBR) NULL, (SUBR) NULL         },
1001   { "vbapmove.a", S(VBAP_MOVING),
1002     TR, 3,  "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm"
1003     "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm",
1004     "aiiim",
1005     (SUBR) vbap_moving_init, (SUBR) vbap_moving },
1006   { "vbapgmove.a",  S(VBAP1_MOVING),      TR, 3,
1007     "zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz"
1008     "zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz", "iiim",
1009     (SUBR) vbap1_moving_init,    (SUBR) vbap1_moving },
1010   { "vbapmove.A", S(VBAPA_MOVING),
1011     TR, 3,  "a[]",  "aiiim",
1012     (SUBR) vbap_moving_init_a, (SUBR) vbap_moving_a },
1013   { "vbapgmove.A",  S(VBAPA1_MOVING),      TR, 3,
1014     "k[]", "iiim",
1015     (SUBR) vbap1_moving_init_a,    (SUBR) vbap1_moving_a },
1016   { "vbapzmove",  S(VBAP_ZAK_MOVING),    ZW|TR, 3,  "",  "iiaiiim",
1017     (SUBR) vbap_zak_moving_init,    (SUBR) vbap_zak_moving  },
1018   { "vbap4move", S(VBAP_MOVING),   TR|_QQ, 3,  "aaaa",
1019    "aiiim",
1020     (SUBR) vbap_moving_init, (SUBR) vbap_moving },
1021   { "vbap8move", S(VBAP_MOVING),
1022     TR|_QQ, 3,  "aaaaaaaa",
1023     "aiiim",
1024     (SUBR) vbap_moving_init, (SUBR) vbap_moving }
1025 
1026 };
1027 
1028 LINKAGE_BUILTIN(vbap_localops)
1029 
1030