Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : #include <stdint.h>
34 : #include "options.h"
35 : #include <assert.h>
36 : #include <math.h>
37 : #include "prot.h"
38 : #include "ivas_prot.h"
39 : #include "ivas_prot_rend.h"
40 : #include "ivas_stat_dec.h"
41 : #ifdef DEBUGGING
42 : #include "debug.h"
43 : #endif
44 : #include "wmc_auto.h"
45 :
46 : /*-----------------------------------------------------------------------*
47 : * Local constants
48 : *-----------------------------------------------------------------------*/
49 :
50 : /* 128 is maximum num_speaker_nodes number. This relates to memory optimization and maximum of triplets:
51 : - triplet indices are unsigned_char (see below structs) --> max triplets is 256
52 : - num_speaker_nodes_internal = num_speaker_nodes + 2 (potential virtual node channels, bottom and top)
53 : - max_num_triplets = 256 = (max_num_ls_internal - 2) * 2 = (max_num_ls) * 2
54 : --> max_num_ls = 128
55 : */
56 : #define VBAP_MAX_NUM_SPEAKER_NODES 128
57 : #define VBAP_MAX_NUM_TRIPLETS 256
58 : #define VBAP_EPSILON 0.001f /* The fairly large epsilon is for detecting planes etc and accounts for rounding issues */
59 : #define VBAP_MAX_PLANES 50
60 : #define VBAP_MAX_HORIZONTAL_GAP_FOR_PLANE_DETECTION 140.0f
61 : /* If a speaker node is found
62 : - above VBAP_NO_VIRTUAL_SPEAKER_NODE_ELE_LIMIT, no virtual node is used
63 : - above VBAP_DISTRIBUTE_VIRTUAL_SPEAKER_NODE_ELE_LIMIT, energy-spreading virtual node is used
64 : - not above VBAP_DISTRIBUTE_VIRTUAL_SPEAKER_NODE_ELE_LIMIT, energy-omitting virtual node is used
65 : Same applies for both elevations and inclinations, i.e., the two half-spheres. */
66 : #define VBAP_NO_VIRTUAL_SPEAKER_NODE_ELE_LIMIT 45.0f
67 : #define VBAP_DISTRIBUTE_VIRTUAL_SPEAKER_NODE_ELE_LIMIT 20.0f
68 : #define VBAP_VIRTUAL_BACK_ELE_LIMIT 45.0f
69 : #define VBAP_NOT_VALID_CONNECTION ( -1 )
70 :
71 : /* Maximum azimuth gap between speaker nodes for detecting zero elevation horizontal plane */
72 : #define VBAP_MAX_HORIZONTAL_GAP 170u
73 :
74 : #define VBAP_SEARCH_SECTOR_SIZE ( 360.0f / ( VBAP_NUM_SEARCH_SECTORS ) )
75 :
76 :
77 : enum VirtualSpeakerNodeType
78 : {
79 : NO_VIRTUAL_SPEAKER_NODE,
80 : VIRTUAL_SPEAKER_NODE_DISCARD_ENERGY,
81 : VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY
82 : };
83 :
84 : enum ConnectionClass
85 : {
86 : REGULAR_CONNECTION,
87 : ELEVATED_PLANE_THIN_TRIANGLE_CONNECTION,
88 : CONNECTION_WITH_SPEAKER_NODE_BEHIND
89 : };
90 :
91 : typedef struct connection_option
92 : {
93 : int16_t chA;
94 : int16_t chB;
95 : float arc;
96 : float arc_weighted;
97 : } ConnectionOption;
98 :
99 : enum SpeakerNodeGroup
100 : {
101 : SPEAKER_NODE_BOTTOM_HALF,
102 : SPEAKER_NODE_HORIZONTAL,
103 : SPEAKER_NODE_TOP_HALF,
104 : SPEAKER_NODE_BACK,
105 : SPEAKER_NODE_ALL
106 : };
107 :
108 : /* Defines a single speaker node */
109 : typedef struct vbap_speaker_node_structure
110 : {
111 : float azi_deg;
112 : float ele_deg;
113 : float unit_vec[3];
114 : enum SpeakerNodeGroup group;
115 :
116 : } VBAP_SPEAKER_NODE;
117 :
118 :
119 : /*-----------------------------------------------------------------------*
120 : * Local function prototypes
121 : *-----------------------------------------------------------------------*/
122 :
123 : static uint8_t vector_matrix_multiply_3x3( const float *src_vector, float matrix[3][3], float *result );
124 :
125 : static void init_speaker_node_direction_data( VBAP_SPEAKER_NODE *speaker_node_data, const float *speaker_node_azi_deg, const float *speaker_node_ele_deg, const int16_t num_speaker_nodes );
126 :
127 : static int16_t determine_virtual_surface_triplets( const int16_t num_speaker_nodes, const VBAP_SPEAKER_NODE *speaker_node_data, int16_t connections[][2], const int16_t max_num_connections, VBAP_VS_TRIPLET *triplets, int16_t initial_search_indices[VBAP_NUM_SEARCH_SECTORS], enum SpeakerNodeGroup allowed_group );
128 :
129 : static void determine_initial_search_indices( const int16_t num_triplets, const float triplet_azidegs[VBAP_MAX_NUM_TRIPLETS], int16_t initial_search_indices[VBAP_NUM_SEARCH_SECTORS] );
130 :
131 : static ivas_error determine_connections( const int16_t num_speaker_nodes, const VBAP_SPEAKER_NODE *speaker_node_data, int16_t connections[][2], const int16_t max_num_connections, int16_t *group1_count, int16_t *group2_start, int16_t *group2_count );
132 :
133 : static void formulate_horizontal_connections( const VBAP_SPEAKER_NODE *speaker_node_data, const int16_t num_speaker_nodes, int16_t connections[][2], int16_t *connection_write_index );
134 :
135 : static ivas_error get_half_sphere_connection_options( const VBAP_SPEAKER_NODE *speaker_node_data, const enum SpeakerNodeGroup group, const int16_t num_speaker_nodes, const int16_t num_non_crossing_planes, const float *non_crossing_plane_elevation_deg, ConnectionOption **connection_options_pr, int16_t *num_connection_options );
136 :
137 : static ivas_error formulate_half_sphere_connections( const VBAP_SPEAKER_NODE *speaker_node_data, const int16_t num_speaker_nodes, const enum SpeakerNodeGroup group, int16_t connections[][2], int16_t *connection_write_index, const int16_t max_num_connections, const int16_t num_non_crossing_planes, const float *non_crossing_plane_elevation_deg );
138 :
139 : static int16_t determine_non_crossing_planes( const int16_t num_speaker_nodes, const VBAP_SPEAKER_NODE *node_data, float *non_crossing_plane_elevation_deg );
140 :
141 : static enum VirtualSpeakerNodeType check_need_of_virtual_speaker_node( VBAP_HANDLE hVBAPdata, const float *speaker_node_azi_deg, const float *speaker_node_ele_deg, enum SpeakerNodeGroup group );
142 :
143 : static int16_t determine_best_triplet_and_gains( VBAP_SEARCH_STRUCT *search_struct, const float panning_unit_vec[3], const int16_t azi_deg, float gains[3] );
144 :
145 : static void determine_virtual_speaker_node_division_gains( const int16_t virtual_speaker_node_index, float *virtual_node_division_gains, int16_t connections[][2], const enum VirtualSpeakerNodeType type, const int16_t max_num_connections, const int16_t num_speaker_nodes, const int16_t use_object_mode );
146 :
147 : static void reorder_triplets( VBAP_VS_TRIPLET *triplets, const int16_t *target_order, const int16_t num_triplets );
148 :
149 :
150 : /*-------------------------------------------------------------------------*
151 : * vbap_init_data()
152 : *
153 : * Initialize VBAP data structure for the speaker node set
154 : *-------------------------------------------------------------------------*/
155 :
156 31140 : ivas_error vbap_init_data(
157 : VBAP_HANDLE *hVBAPdata, /* i/o: handle for VBAP data structure that will be initialized */
158 : const float *speaker_node_azi_deg, /* i : vector of speaker node azimuths (positive left) */
159 : const float *speaker_node_ele_deg, /* i : vector of speaker node elevations (positive up) */
160 : const int16_t num_speaker_nodes, /* i : number of speaker nodes in the set */
161 : const IVAS_FORMAT ivas_format /* i : IVAS format */
162 : )
163 : {
164 : /* Variables */
165 : int16_t connections[VBAP_MAX_NUM_SPEAKER_NODES][2];
166 : int16_t max_num_connections;
167 : int16_t is_success;
168 : int16_t connection_group1_count;
169 : int16_t connection_group2_start;
170 : int16_t connection_group2_count;
171 : enum VirtualSpeakerNodeType virtual_top_type;
172 : enum VirtualSpeakerNodeType virtual_bottom_type;
173 : enum VirtualSpeakerNodeType virtual_back_type;
174 : float speaker_node_azi_deg_internal[VBAP_MAX_NUM_SPEAKER_NODES];
175 : float speaker_node_ele_deg_internal[VBAP_MAX_NUM_SPEAKER_NODES];
176 : VBAP_SPEAKER_NODE speaker_node_data[VBAP_MAX_NUM_SPEAKER_NODES];
177 : VBAP_DATA *vbap;
178 : ivas_error error;
179 :
180 31140 : push_wmops( "vbap_init" );
181 :
182 : /* Basic init checks */
183 : /* If the requested layout is invalid, hVBAPdata is set to NULL and the signal will
184 : * be distributed with an equal gain into all output channels.
185 : * The surrounding code needs to handle the NULL pointer properly. */
186 31140 : if ( num_speaker_nodes > VBAP_MAX_NUM_SPEAKER_NODES || num_speaker_nodes < 3 )
187 : {
188 0 : hVBAPdata = NULL;
189 0 : pop_wmops();
190 0 : return IVAS_ERR_OK;
191 : }
192 :
193 31140 : if ( !speaker_node_azi_deg || !speaker_node_ele_deg )
194 : {
195 0 : hVBAPdata = NULL;
196 0 : return IVAS_ERR_OK;
197 : }
198 :
199 : /* Allocate VBAP structure */
200 31140 : if ( ( vbap = (VBAP_HANDLE) malloc( sizeof( VBAP_DATA ) ) ) == NULL )
201 : {
202 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
203 : }
204 :
205 31140 : is_success = 1;
206 31140 : vbap->bottom_virtual_speaker_node_index = -1;
207 31140 : vbap->top_virtual_speaker_node_index = -1;
208 31140 : vbap->back_virtual_speaker_node_index = -1;
209 31140 : vbap->bottom_virtual_speaker_node_division_gains = NULL;
210 31140 : vbap->top_virtual_speaker_node_division_gains = NULL;
211 31140 : vbap->back_virtual_speaker_node_division_gains = NULL;
212 31140 : vbap->object_mode_bottom_virtual_speaker_node_division_gains = NULL;
213 31140 : vbap->object_mode_top_virtual_speaker_node_division_gains = NULL;
214 31140 : vbap->object_mode_back_virtual_speaker_node_division_gains = NULL;
215 31140 : vbap->num_speaker_nodes = num_speaker_nodes;
216 31140 : vbap->num_speaker_nodes_internal = num_speaker_nodes;
217 :
218 : /* Check if the speaker node setup needs a virtual top or bottom node
219 : (function also increments vbap->num_speaker_nodes_internal when necessary) */
220 31140 : virtual_bottom_type = check_need_of_virtual_speaker_node( vbap, speaker_node_azi_deg, speaker_node_ele_deg, SPEAKER_NODE_BOTTOM_HALF );
221 31140 : virtual_top_type = check_need_of_virtual_speaker_node( vbap, speaker_node_azi_deg, speaker_node_ele_deg, SPEAKER_NODE_TOP_HALF );
222 31140 : virtual_back_type = check_need_of_virtual_speaker_node( vbap, speaker_node_azi_deg, speaker_node_ele_deg, SPEAKER_NODE_BACK );
223 :
224 : /* Init internal speaker node configuration, which is the original configuration
225 : potentially appended with virtual top and/or bottom loudspeakers */
226 31140 : mvr2r( speaker_node_azi_deg, speaker_node_azi_deg_internal, num_speaker_nodes );
227 31140 : mvr2r( speaker_node_ele_deg, speaker_node_ele_deg_internal, num_speaker_nodes );
228 :
229 31140 : if ( is_success && virtual_bottom_type != NO_VIRTUAL_SPEAKER_NODE )
230 : {
231 31136 : if ( ( vbap->bottom_virtual_speaker_node_division_gains = (float *) malloc( num_speaker_nodes * sizeof( float ) ) ) == NULL )
232 : {
233 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
234 : }
235 31136 : set_zero( vbap->bottom_virtual_speaker_node_division_gains, num_speaker_nodes );
236 31136 : is_success &= vbap->bottom_virtual_speaker_node_division_gains != NULL;
237 :
238 31136 : if ( ivas_format == MASA_ISM_FORMAT )
239 : {
240 11727 : if ( ( vbap->object_mode_bottom_virtual_speaker_node_division_gains = (float *) malloc( num_speaker_nodes * sizeof( float ) ) ) == NULL )
241 : {
242 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
243 : }
244 11727 : set_zero( vbap->object_mode_bottom_virtual_speaker_node_division_gains, num_speaker_nodes );
245 11727 : is_success &= vbap->object_mode_bottom_virtual_speaker_node_division_gains != NULL;
246 : }
247 :
248 31136 : speaker_node_azi_deg_internal[vbap->bottom_virtual_speaker_node_index] = 0.0f;
249 31136 : speaker_node_ele_deg_internal[vbap->bottom_virtual_speaker_node_index] = -90.0f;
250 : }
251 :
252 31140 : if ( is_success && virtual_top_type != NO_VIRTUAL_SPEAKER_NODE )
253 : {
254 31136 : if ( ( vbap->top_virtual_speaker_node_division_gains = (float *) malloc( num_speaker_nodes * sizeof( float ) ) ) == NULL )
255 : {
256 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
257 : }
258 31136 : set_zero( vbap->top_virtual_speaker_node_division_gains, num_speaker_nodes );
259 31136 : is_success &= vbap->top_virtual_speaker_node_division_gains != NULL;
260 :
261 31136 : if ( ivas_format == MASA_ISM_FORMAT )
262 : {
263 11727 : if ( ( vbap->object_mode_top_virtual_speaker_node_division_gains = (float *) malloc( num_speaker_nodes * sizeof( float ) ) ) == NULL )
264 : {
265 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
266 : }
267 11727 : set_zero( vbap->object_mode_top_virtual_speaker_node_division_gains, num_speaker_nodes );
268 11727 : is_success &= vbap->object_mode_top_virtual_speaker_node_division_gains != NULL;
269 : }
270 :
271 31136 : speaker_node_azi_deg_internal[vbap->top_virtual_speaker_node_index] = 0.0f;
272 31136 : speaker_node_ele_deg_internal[vbap->top_virtual_speaker_node_index] = 90.0f;
273 : }
274 :
275 31140 : if ( is_success && virtual_back_type != NO_VIRTUAL_SPEAKER_NODE )
276 : {
277 0 : if ( ( vbap->back_virtual_speaker_node_division_gains = (float *) malloc( num_speaker_nodes * sizeof( float ) ) ) == NULL )
278 : {
279 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
280 : }
281 0 : set_zero( vbap->back_virtual_speaker_node_division_gains, num_speaker_nodes );
282 0 : is_success &= vbap->back_virtual_speaker_node_division_gains != NULL;
283 :
284 0 : if ( ivas_format == MASA_ISM_FORMAT )
285 : {
286 0 : if ( ( vbap->object_mode_back_virtual_speaker_node_division_gains = (float *) malloc( num_speaker_nodes * sizeof( float ) ) ) == NULL )
287 : {
288 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
289 : }
290 0 : set_zero( vbap->object_mode_back_virtual_speaker_node_division_gains, num_speaker_nodes );
291 0 : is_success &= vbap->object_mode_back_virtual_speaker_node_division_gains != NULL;
292 : }
293 0 : speaker_node_azi_deg_internal[vbap->back_virtual_speaker_node_index] = 180.0f;
294 0 : speaker_node_ele_deg_internal[vbap->back_virtual_speaker_node_index] = 0.0f;
295 : }
296 :
297 31140 : init_speaker_node_direction_data( speaker_node_data, speaker_node_azi_deg_internal, speaker_node_ele_deg_internal, vbap->num_speaker_nodes_internal );
298 :
299 : /* Allocate and determine node-node connections */
300 31140 : max_num_connections = ( vbap->num_speaker_nodes_internal - 2 ) * 3; /* Theoretical maximum */
301 :
302 31140 : if ( ( error = determine_connections( vbap->num_speaker_nodes_internal, speaker_node_data, connections, max_num_connections, &connection_group1_count, &connection_group2_start, &connection_group2_count ) ) != IVAS_ERR_OK )
303 : {
304 0 : return error;
305 : }
306 :
307 : /* Allocate and determine virtual surface speaker node triplets */
308 31140 : if ( is_success )
309 : {
310 : int16_t ch;
311 31140 : int16_t speaker_nodes_group1_internal = 0;
312 31140 : int16_t speaker_nodes_group2_internal = 0;
313 31140 : int16_t speaker_nodes_horiz_internal = 0;
314 31140 : uint8_t loop_done = 0;
315 :
316 : /* Count nodes in different groups to reserve correct memory */
317 348576 : for ( ch = 0; ch < vbap->num_speaker_nodes_internal && !loop_done; ch++ )
318 : {
319 317436 : switch ( speaker_node_data[ch].group )
320 : {
321 0 : case SPEAKER_NODE_ALL:
322 : /* If there is even one speaker belonging to "all" group, then all speakers belong to the "all" group.
323 : * We can skip further counts here. */
324 0 : speaker_nodes_group1_internal = vbap->num_speaker_nodes_internal;
325 0 : loop_done = 1;
326 0 : break;
327 31184 : case SPEAKER_NODE_BOTTOM_HALF:
328 31184 : speaker_nodes_group1_internal++;
329 31184 : break;
330 103554 : case SPEAKER_NODE_TOP_HALF:
331 103554 : speaker_nodes_group2_internal++;
332 103554 : break;
333 182698 : case SPEAKER_NODE_HORIZONTAL:
334 : case SPEAKER_NODE_BACK:
335 182698 : speaker_nodes_group1_internal++;
336 182698 : speaker_nodes_group2_internal++;
337 182698 : speaker_nodes_horiz_internal++;
338 182698 : break;
339 : }
340 317436 : }
341 :
342 31140 : if ( ( vbap->search_struct[0].triplets = (VBAP_VS_TRIPLET *) malloc( ( ( speaker_nodes_group1_internal - 2 ) * 2 - ( max( 0, ( speaker_nodes_horiz_internal - 2 ) ) ) ) * sizeof( VBAP_VS_TRIPLET ) ) ) == NULL )
343 : {
344 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
345 : }
346 31140 : is_success &= vbap->search_struct[0].triplets != NULL;
347 :
348 31140 : if ( speaker_nodes_group2_internal > 0 )
349 : {
350 31140 : vbap->num_search_structs = 2;
351 31140 : if ( ( vbap->search_struct[1].triplets = (VBAP_VS_TRIPLET *) malloc( ( ( speaker_nodes_group2_internal - 2 ) * 2 - ( max( 0, ( speaker_nodes_horiz_internal - 2 ) ) ) ) * sizeof( VBAP_VS_TRIPLET ) ) ) == NULL )
352 : {
353 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
354 : }
355 31140 : is_success &= vbap->search_struct[1].triplets != NULL;
356 : }
357 : else
358 : {
359 0 : vbap->num_search_structs = 1;
360 0 : vbap->search_struct[1].triplets = NULL;
361 : }
362 : }
363 :
364 31140 : if ( is_success )
365 : {
366 31140 : if ( vbap->num_search_structs == 1 )
367 : {
368 : /* If all speaker nodes belong to ALL set, then we only create one triplet set and search structure */
369 0 : vbap->search_struct[0].num_triplets = determine_virtual_surface_triplets( vbap->num_speaker_nodes_internal, speaker_node_data, connections, max_num_connections, vbap->search_struct[0].triplets, vbap->search_struct[0].initial_search_indices, SPEAKER_NODE_ALL );
370 : }
371 : else
372 : {
373 : /* Otherwise, we have two sets and can handle them separately for more opmitized processing. */
374 31140 : vbap->search_struct[0].num_triplets = determine_virtual_surface_triplets( vbap->num_speaker_nodes_internal, speaker_node_data, connections, connection_group1_count, vbap->search_struct[0].triplets, vbap->search_struct[0].initial_search_indices, SPEAKER_NODE_BOTTOM_HALF );
375 31140 : vbap->search_struct[1].num_triplets = determine_virtual_surface_triplets( vbap->num_speaker_nodes_internal, speaker_node_data, connections + connection_group2_start, connection_group2_count, vbap->search_struct[1].triplets, vbap->search_struct[1].initial_search_indices, SPEAKER_NODE_TOP_HALF );
376 : }
377 : }
378 :
379 : /* Determine how the virtual node gains should be distributed to real nodes, if necessary (checked within function). */
380 31140 : if ( is_success )
381 : {
382 31140 : determine_virtual_speaker_node_division_gains( vbap->top_virtual_speaker_node_index, vbap->top_virtual_speaker_node_division_gains, connections, virtual_top_type, max_num_connections, num_speaker_nodes, 0 );
383 31140 : determine_virtual_speaker_node_division_gains( vbap->bottom_virtual_speaker_node_index, vbap->bottom_virtual_speaker_node_division_gains, connections, virtual_bottom_type, max_num_connections, num_speaker_nodes, 0 );
384 31140 : determine_virtual_speaker_node_division_gains( vbap->back_virtual_speaker_node_index, vbap->back_virtual_speaker_node_division_gains, connections, virtual_back_type, max_num_connections, num_speaker_nodes, 0 );
385 31140 : if ( ivas_format == MASA_ISM_FORMAT )
386 : {
387 11727 : determine_virtual_speaker_node_division_gains( vbap->top_virtual_speaker_node_index, vbap->object_mode_top_virtual_speaker_node_division_gains, connections, virtual_top_type == NO_VIRTUAL_SPEAKER_NODE ? NO_VIRTUAL_SPEAKER_NODE : VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY, max_num_connections, num_speaker_nodes, 1 );
388 11727 : determine_virtual_speaker_node_division_gains( vbap->bottom_virtual_speaker_node_index, vbap->object_mode_bottom_virtual_speaker_node_division_gains, connections, virtual_bottom_type == NO_VIRTUAL_SPEAKER_NODE ? NO_VIRTUAL_SPEAKER_NODE : VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY, max_num_connections, num_speaker_nodes, 1 );
389 11727 : determine_virtual_speaker_node_division_gains( vbap->back_virtual_speaker_node_index, vbap->object_mode_back_virtual_speaker_node_division_gains, connections, virtual_back_type == NO_VIRTUAL_SPEAKER_NODE ? NO_VIRTUAL_SPEAKER_NODE : VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY, max_num_connections, num_speaker_nodes, 1 );
390 : }
391 : }
392 :
393 31140 : pop_wmops();
394 :
395 31140 : if ( is_success )
396 : {
397 31140 : *hVBAPdata = vbap;
398 : }
399 : else
400 : {
401 0 : vbap_free_data( &vbap );
402 : }
403 :
404 31140 : return IVAS_ERR_OK;
405 : }
406 :
407 :
408 : /*-------------------------------------------------------------------------*
409 : * vbap_free_data()
410 : *
411 : * Free VBAP data structure
412 : *-------------------------------------------------------------------------*/
413 :
414 120079 : void vbap_free_data(
415 : VBAP_HANDLE *hVBAPdata /* i/o: VBAP handle to be freed */
416 : )
417 : {
418 120079 : if ( hVBAPdata == NULL || *hVBAPdata == NULL )
419 : {
420 88939 : return;
421 : }
422 :
423 31140 : if ( ( *hVBAPdata )->bottom_virtual_speaker_node_division_gains != NULL )
424 : {
425 31136 : free( ( *hVBAPdata )->bottom_virtual_speaker_node_division_gains );
426 : }
427 31140 : if ( ( *hVBAPdata )->top_virtual_speaker_node_division_gains != NULL )
428 : {
429 31136 : free( ( *hVBAPdata )->top_virtual_speaker_node_division_gains );
430 : }
431 31140 : if ( ( *hVBAPdata )->back_virtual_speaker_node_division_gains != NULL )
432 : {
433 0 : free( ( *hVBAPdata )->back_virtual_speaker_node_division_gains );
434 : }
435 31140 : if ( ( *hVBAPdata )->object_mode_bottom_virtual_speaker_node_division_gains != NULL )
436 : {
437 11727 : free( ( *hVBAPdata )->object_mode_bottom_virtual_speaker_node_division_gains );
438 : }
439 31140 : if ( ( *hVBAPdata )->object_mode_top_virtual_speaker_node_division_gains != NULL )
440 : {
441 11727 : free( ( *hVBAPdata )->object_mode_top_virtual_speaker_node_division_gains );
442 : }
443 31140 : if ( ( *hVBAPdata )->object_mode_back_virtual_speaker_node_division_gains != NULL )
444 : {
445 0 : free( ( *hVBAPdata )->object_mode_back_virtual_speaker_node_division_gains );
446 : }
447 31140 : if ( ( *hVBAPdata )->search_struct[0].triplets != NULL )
448 : {
449 31140 : free( ( *hVBAPdata )->search_struct[0].triplets );
450 : }
451 31140 : if ( ( *hVBAPdata )->num_search_structs == 2 && ( *hVBAPdata )->search_struct[1].triplets != NULL )
452 : {
453 31140 : free( ( *hVBAPdata )->search_struct[1].triplets );
454 : }
455 :
456 31140 : free( *hVBAPdata );
457 31140 : *hVBAPdata = NULL;
458 :
459 31140 : return;
460 : }
461 :
462 :
463 : /*-------------------------------------------------------------------------*
464 : * vbap_determine_gains()
465 : *
466 : * Obtain panning gains for all speaker nodes based on the given direction
467 : *-------------------------------------------------------------------------*/
468 :
469 138198182 : void vbap_determine_gains(
470 : const VBAP_HANDLE hVBAPdata, /* i : prepared VBAP structure */
471 : float *gains, /* o : gain vector for loudspeakers for given direction */
472 : const int16_t azi_deg, /* i : azimuth in degrees for panning direction (positive left) */
473 : const int16_t ele_deg, /* i : elevation in degrees for panning direction (positive up) */
474 : const int16_t use_object_mode /* i : select between object mode panning and spatial mode panning */
475 : )
476 : {
477 : /* This function formulates gains for the given angle. The triplet-selection has been pre-formulated. */
478 : int16_t ch, ch2;
479 : int16_t triplet_ch;
480 : int16_t triplet_index;
481 : float panning_unit_vec[3];
482 : float gain_triplet[3];
483 : float norm_value;
484 : float gain_ene;
485 : float azi_rad;
486 : float ele_rad;
487 : float azi_temp;
488 : float ele_temp;
489 : int16_t num_speaker_nodes;
490 : int16_t bottom_virtual_speaker_node_index;
491 : int16_t top_virtual_speaker_node_index;
492 : int16_t back_virtual_speaker_node_index;
493 : VBAP_VS_TRIPLET *selected_triplet;
494 : float *bottom_virtual_speaker_node_division_gains;
495 : float *top_virtual_speaker_node_division_gains;
496 : float *back_virtual_speaker_node_division_gains;
497 :
498 : #ifdef DEBUGGING
499 : assert( hVBAPdata != NULL && "VBAP gain determination requires initialized structure." );
500 : assert( gains != NULL && "VBAP gain determination requires reserved memory for gain output." );
501 : #endif
502 :
503 138198182 : push_wmops( "vbap_gains" );
504 138198182 : num_speaker_nodes = hVBAPdata->num_speaker_nodes;
505 138198182 : bottom_virtual_speaker_node_index = hVBAPdata->bottom_virtual_speaker_node_index;
506 138198182 : top_virtual_speaker_node_index = hVBAPdata->top_virtual_speaker_node_index;
507 138198182 : back_virtual_speaker_node_index = hVBAPdata->back_virtual_speaker_node_index;
508 138198182 : if ( use_object_mode )
509 : {
510 3355509 : bottom_virtual_speaker_node_division_gains = hVBAPdata->object_mode_bottom_virtual_speaker_node_division_gains;
511 3355509 : top_virtual_speaker_node_division_gains = hVBAPdata->object_mode_top_virtual_speaker_node_division_gains;
512 3355509 : back_virtual_speaker_node_division_gains = hVBAPdata->object_mode_back_virtual_speaker_node_division_gains;
513 : }
514 : else
515 : {
516 134842673 : bottom_virtual_speaker_node_division_gains = hVBAPdata->bottom_virtual_speaker_node_division_gains;
517 134842673 : top_virtual_speaker_node_division_gains = hVBAPdata->top_virtual_speaker_node_division_gains;
518 134842673 : back_virtual_speaker_node_division_gains = hVBAPdata->back_virtual_speaker_node_division_gains;
519 : }
520 :
521 138198182 : panning_wrap_angles( (float) azi_deg, (float) ele_deg, &azi_temp, &ele_temp );
522 138198182 : azi_rad = azi_temp * PI_OVER_180;
523 138198182 : ele_rad = ele_temp * PI_OVER_180;
524 :
525 138198182 : panning_unit_vec[0] = cosf( azi_rad ) * cosf( ele_rad );
526 138198182 : panning_unit_vec[1] = sinf( azi_rad ) * cosf( ele_rad );
527 138198182 : panning_unit_vec[2] = sinf( ele_rad );
528 :
529 : /* Find the best VS triplet and speaker node gains for the panning direction using the prepared search structures. */
530 138198182 : if ( hVBAPdata->num_search_structs == 2 && ele_deg > 0 )
531 : {
532 45891116 : triplet_index = determine_best_triplet_and_gains( &( hVBAPdata->search_struct[1] ), panning_unit_vec, azi_deg, gain_triplet );
533 45891116 : selected_triplet = &hVBAPdata->search_struct[1].triplets[triplet_index];
534 : }
535 : else
536 : {
537 92307066 : triplet_index = determine_best_triplet_and_gains( &( hVBAPdata->search_struct[0] ), panning_unit_vec, azi_deg, gain_triplet );
538 92307066 : selected_triplet = &hVBAPdata->search_struct[0].triplets[triplet_index];
539 : }
540 :
541 : /* Normalize to unit energy */
542 138198182 : gain_ene = 1e-12f; /* Add small value to avoid divide by zero. */
543 552792728 : for ( ch = 0; ch < 3; ch++ )
544 : {
545 414594546 : gain_ene += gain_triplet[ch] * gain_triplet[ch];
546 : }
547 :
548 138198182 : norm_value = inv_sqrt( gain_ene );
549 :
550 552792728 : for ( ch = 0; ch < 3; ch++ )
551 : {
552 414594546 : gain_triplet[ch] *= norm_value;
553 :
554 : /* Sanity check for rounding issues */
555 414594546 : if ( gain_triplet[ch] < 0.0f )
556 : {
557 23404 : gain_triplet[ch] = 0.0f;
558 : }
559 : }
560 :
561 : /* Flush gain target */
562 138198182 : set_zero( gains, num_speaker_nodes );
563 :
564 : /* Map gain triplet (internal speaker node configuration) to speaker node output (actual speaker node configuration) */
565 552792728 : for ( ch = 0; ch < 3; ch++ )
566 : {
567 414594546 : triplet_ch = selected_triplet->speaker_node[ch];
568 :
569 414594546 : if ( triplet_ch == bottom_virtual_speaker_node_index )
570 : {
571 761113738 : for ( ch2 = 0; ch2 < num_speaker_nodes; ch2++ )
572 : {
573 673190387 : gains[ch2] += bottom_virtual_speaker_node_division_gains[ch2] * gain_triplet[ch];
574 : }
575 : }
576 326671195 : else if ( triplet_ch == top_virtual_speaker_node_index )
577 : {
578 229460876 : for ( ch2 = 0; ch2 < num_speaker_nodes; ch2++ )
579 : {
580 200739994 : gains[ch2] += top_virtual_speaker_node_division_gains[ch2] * gain_triplet[ch];
581 : }
582 : }
583 297950313 : else if ( triplet_ch == back_virtual_speaker_node_index )
584 : {
585 0 : for ( ch2 = 0; ch2 < num_speaker_nodes; ch2++ )
586 : {
587 0 : gains[ch2] += back_virtual_speaker_node_division_gains[ch2] * gain_triplet[ch];
588 : }
589 : }
590 : else
591 : {
592 297950313 : gains[triplet_ch] += gain_triplet[ch];
593 : }
594 : }
595 :
596 138198182 : pop_wmops();
597 :
598 138198182 : return;
599 : }
600 :
601 :
602 : /*-----------------------------------------------------------------------*
603 : * Local functions
604 : *-----------------------------------------------------------------------*/
605 :
606 : /*-------------------------------------------------------------------------*
607 : * vbap_crossp()
608 : *
609 : * 3-by-3 vector cross product
610 : *-------------------------------------------------------------------------*/
611 :
612 4092870 : static void vbap_crossp(
613 : const float *vec1, /* i : input vector 1 */
614 : const float *vec2, /* i : input vector 2 */
615 : float *crossProduct /* o : cross product output */
616 : )
617 : {
618 4092870 : crossProduct[0] = ( vec1[1] * vec2[2] ) - ( vec1[2] * vec2[1] );
619 4092870 : crossProduct[1] = ( vec1[2] * vec2[0] ) - ( vec1[0] * vec2[2] );
620 4092870 : crossProduct[2] = ( vec1[0] * vec2[1] ) - ( vec1[1] * vec2[0] );
621 :
622 4092870 : return;
623 : }
624 :
625 :
626 : /*-------------------------------------------------------------------------*
627 : * vector_matrix_multiply_3x3()
628 : *
629 : * 3-by-3 vector multiply with matrix
630 : *-------------------------------------------------------------------------*/
631 :
632 : /*! r: Status result if triplet is usable for panning. Allows early exit. */
633 391100711 : static uint8_t vector_matrix_multiply_3x3(
634 : const float *src_vector, /* i : input vector */
635 : float matrix[3][3], /* i : input matrix */
636 : float *result /* o : output vector */
637 : )
638 : {
639 391100711 : result[0] = src_vector[0] * matrix[0][0];
640 391100711 : result[0] += src_vector[1] * matrix[1][0];
641 391100711 : result[0] += src_vector[2] * matrix[2][0];
642 :
643 391100711 : if ( result[0] < -0.01f )
644 : {
645 26041365 : return 0;
646 : }
647 :
648 365059346 : result[1] = src_vector[0] * matrix[0][1];
649 365059346 : result[1] += src_vector[1] * matrix[1][1];
650 365059346 : result[1] += src_vector[2] * matrix[2][1];
651 :
652 365059346 : if ( result[1] < -0.01f )
653 : {
654 162780348 : return 0;
655 : }
656 :
657 202278998 : result[2] = src_vector[0] * matrix[0][2];
658 202278998 : result[2] += src_vector[1] * matrix[1][2];
659 202278998 : result[2] += src_vector[2] * matrix[2][2];
660 :
661 202278998 : if ( result[2] < -0.01f )
662 : {
663 62064010 : return 0;
664 : }
665 :
666 140214988 : return 1;
667 : }
668 :
669 :
670 : /*----------------------------------------------------------------------------------------------*
671 : * determine_best_triplet_and_gains()
672 : *
673 : * Determine the best speaker node triplet and associated gains for panning to defined direction
674 : *----------------------------------------------------------------------------------------------*/
675 :
676 : /*! r: triplet id */
677 138198182 : static int16_t determine_best_triplet_and_gains(
678 : VBAP_SEARCH_STRUCT *search_struct, /* i : VBAP search struct */
679 : const float panning_unit_vec[3], /* i : panning unit vector */
680 : const int16_t azi_deg, /* i : panning azimuth */
681 : float gains[3] /* o : panning gains */
682 : )
683 : {
684 : int16_t i, tr, k;
685 : uint8_t triplet_ok;
686 : int16_t best_triplet;
687 : float best_min_gain;
688 : float min_gain_this;
689 : float unnormalized_gains[3];
690 : int16_t sector;
691 : int16_t first_triplet;
692 : int16_t jump;
693 : int16_t num_triplets;
694 :
695 138198182 : num_triplets = search_struct->num_triplets;
696 138198182 : best_min_gain = -999.9f;
697 138198182 : best_triplet = 0;
698 138198182 : set_zero( gains, 3 );
699 :
700 : /* Determine the correct search sector for that target panning direction using an optimized algorithm for
701 : * the chosen four sectors. */
702 138198182 : if ( abs( azi_deg ) > 90 )
703 : {
704 56474678 : sector = azi_deg < 0 ? 2 : 1;
705 : }
706 : else
707 : {
708 81723504 : sector = azi_deg < 0 ? 3 : 0;
709 : }
710 138198182 : first_triplet = search_struct->initial_search_indices[sector];
711 :
712 138198182 : tr = first_triplet;
713 138198182 : jump = 1;
714 387238655 : for ( i = 0; i < num_triplets; i++ )
715 : {
716 387215251 : triplet_ok = vector_matrix_multiply_3x3( panning_unit_vec, search_struct->triplets[tr].inverse_matrix, unnormalized_gains );
717 387215251 : if ( triplet_ok )
718 : {
719 140214988 : min_gain_this = min( ( min( unnormalized_gains[0], unnormalized_gains[1] ) ), unnormalized_gains[2] );
720 :
721 140214988 : if ( min_gain_this > best_min_gain )
722 : {
723 140202256 : best_min_gain = min_gain_this;
724 140202256 : best_triplet = tr;
725 560809024 : for ( k = 0; k < 3; k++ )
726 : {
727 420606768 : gains[k] = unnormalized_gains[k];
728 : }
729 140202256 : if ( !( best_min_gain < 0.00f ) )
730 : {
731 138174778 : return best_triplet;
732 : }
733 : }
734 : }
735 249040473 : tr = first_triplet + jump;
736 249040473 : if ( tr < 0 )
737 : {
738 2326008 : tr += num_triplets;
739 : }
740 246714465 : else if ( tr >= num_triplets )
741 : {
742 63506535 : tr -= num_triplets;
743 : }
744 :
745 249040473 : jump *= -1;
746 249040473 : if ( jump > 0 )
747 : {
748 99967817 : jump += 1;
749 : }
750 : }
751 :
752 23404 : return best_triplet;
753 : }
754 :
755 : /*-------------------------------------------------------------------------*
756 : * determine_virtual_speaker_node_division_gains()
757 : *
758 : * Determines how the virtual node gains are distributed to real nodes
759 : *-------------------------------------------------------------------------*/
760 :
761 128601 : static void determine_virtual_speaker_node_division_gains(
762 : const int16_t virtual_speaker_node_index, /* i : virtual speaker node index */
763 : float *virtual_node_division_gains, /* o : virtual speaker node division gains */
764 : int16_t connections[][2], /* i : vector of all connections */
765 : const enum VirtualSpeakerNodeType type, /* i : virtual speaker node typel */
766 : const int16_t max_num_connections, /* i : max number of connections */
767 : const int16_t num_speaker_nodes, /* i : max number of speaker nodes */
768 : const int16_t use_object_mode /* i : use VBAP in object panning mode vs. spatial panning mode */
769 : )
770 : {
771 : /* When node type is VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY, the gains of the virtual node
772 : are distributed to all neighboring real speaker nodes. An amplitude-division
773 : instead of energy division is utilized just in case to avoid excessive emphasis
774 : on the coherent distributed sound. */
775 : int16_t c, ch;
776 : float sum_val;
777 :
778 128601 : if ( type == VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY )
779 : {
780 1167516 : for ( c = 0; c < max_num_connections; c++ )
781 : {
782 1122711 : if ( connections[c][0] != VBAP_NOT_VALID_CONNECTION )
783 : {
784 1122711 : int16_t connection_node = -1;
785 1122711 : if ( connections[c][0] == virtual_speaker_node_index )
786 : {
787 0 : connection_node = connections[c][1];
788 : }
789 1122711 : else if ( connections[c][1] == virtual_speaker_node_index )
790 : {
791 208918 : connection_node = connections[c][0];
792 : }
793 :
794 : /* The second condition allows division gains only to actual loudspeakers */
795 1122711 : if ( connection_node >= 0 && ( connection_node < num_speaker_nodes ) )
796 : {
797 208918 : virtual_node_division_gains[connection_node] = 1.0f;
798 : }
799 : }
800 : }
801 :
802 44805 : sum_val = 0.0f;
803 419042 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
804 : {
805 374237 : sum_val += virtual_node_division_gains[ch];
806 : }
807 :
808 419042 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
809 : {
810 374237 : virtual_node_division_gains[ch] /= sum_val;
811 374237 : if ( use_object_mode )
812 : {
813 181606 : virtual_node_division_gains[ch] = powf( virtual_node_division_gains[ch], 0.8f );
814 : }
815 : }
816 : }
817 :
818 128601 : return;
819 : }
820 :
821 :
822 : /*-------------------------------------------------------------------------*
823 : * check_need_of_virtual_speaker_node()
824 : *
825 : * Check if virtual speaker node is required
826 : *-------------------------------------------------------------------------*/
827 :
828 : /*! r: virtual speaker node type */
829 93420 : static enum VirtualSpeakerNodeType check_need_of_virtual_speaker_node(
830 : VBAP_HANDLE hVBAPdata, /* i/o: VBAP structure */
831 : const float *speaker_node_azi_deg, /* i : vector of speaker node azimuths */
832 : const float *speaker_node_ele_deg, /* i : vector of speaker node elevations */
833 : enum SpeakerNodeGroup group /* i : group of speaker nodes where this belongs */
834 : )
835 : {
836 : int16_t ch;
837 93420 : float max_elevation = 0.0f;
838 :
839 : /* The following considers if SPEAKER_NODE_BACK virtual speaker is needed */
840 93420 : if ( group == SPEAKER_NODE_BACK )
841 : {
842 31140 : int16_t virtual_back_needed = 1;
843 31140 : const float virtual_back_epsilon = -0.0175f; /* Corresponds to approximately 91 degrees, see code below */
844 122730 : for ( ch = 0; ch < hVBAPdata->num_speaker_nodes; ch++ )
845 : {
846 122730 : if ( fabsf( speaker_node_ele_deg[ch] ) < VBAP_VIRTUAL_BACK_ELE_LIMIT )
847 : {
848 122726 : if ( cosf( speaker_node_azi_deg[ch] * PI_OVER_180 ) < virtual_back_epsilon )
849 : {
850 31140 : virtual_back_needed = 0;
851 31140 : break;
852 : }
853 : }
854 : }
855 :
856 31140 : if ( virtual_back_needed )
857 : {
858 0 : hVBAPdata->back_virtual_speaker_node_index = hVBAPdata->num_speaker_nodes_internal;
859 0 : hVBAPdata->num_speaker_nodes_internal++;
860 0 : return VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY;
861 : }
862 :
863 31140 : return NO_VIRTUAL_SPEAKER_NODE; /* No virtual back needed */
864 : }
865 :
866 : /* The following considers if TOP or BOTTOM virtual speaker is needed */
867 572608 : for ( ch = 0; ch < hVBAPdata->num_speaker_nodes; ch++ )
868 : {
869 510328 : if ( group == SPEAKER_NODE_TOP_HALF )
870 : {
871 255164 : if ( speaker_node_ele_deg[ch] > max_elevation )
872 : {
873 21347 : max_elevation = speaker_node_ele_deg[ch];
874 : }
875 : }
876 : else
877 : {
878 255164 : if ( ( -speaker_node_ele_deg[ch] ) > max_elevation )
879 : {
880 16 : max_elevation = -speaker_node_ele_deg[ch];
881 : }
882 : }
883 : }
884 :
885 62280 : if ( max_elevation > VBAP_NO_VIRTUAL_SPEAKER_NODE_ELE_LIMIT - VBAP_EPSILON )
886 : {
887 8 : return NO_VIRTUAL_SPEAKER_NODE;
888 : }
889 :
890 : /* Use virtual node */
891 62272 : if ( group == SPEAKER_NODE_BOTTOM_HALF )
892 : {
893 31136 : hVBAPdata->bottom_virtual_speaker_node_index = hVBAPdata->num_speaker_nodes_internal;
894 : }
895 : else
896 : {
897 31136 : hVBAPdata->top_virtual_speaker_node_index = hVBAPdata->num_speaker_nodes_internal;
898 : }
899 :
900 62272 : hVBAPdata->num_speaker_nodes_internal++;
901 62272 : if ( max_elevation > VBAP_DISTRIBUTE_VIRTUAL_SPEAKER_NODE_ELE_LIMIT - VBAP_EPSILON )
902 : {
903 21351 : return VIRTUAL_SPEAKER_NODE_DISTRIBUTE_ENERGY;
904 : }
905 :
906 40921 : return VIRTUAL_SPEAKER_NODE_DISCARD_ENERGY;
907 : }
908 :
909 :
910 : /*-------------------------------------------------------------------------*
911 : * init_speaker_node_direction_data()
912 : *
913 : * Initialize speaker node data
914 : *-------------------------------------------------------------------------*/
915 :
916 31140 : static void init_speaker_node_direction_data(
917 : VBAP_SPEAKER_NODE *speaker_node_data, /* o : storage for speaker node data */
918 : const float *speaker_node_azi_deg, /* i : vector of speaker node azimuths */
919 : const float *speaker_node_ele_deg, /* i : vector of speaker node elevations */
920 : const int16_t num_speaker_nodes /* i : number of speaker nodes */
921 : )
922 : {
923 : int16_t ch;
924 : float azi_rad;
925 : float ele_rad;
926 31140 : int16_t num_horiz = 0;
927 31140 : uint8_t in_all_mode = TRUE;
928 :
929 348576 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
930 : {
931 317436 : speaker_node_data[ch].azi_deg = speaker_node_azi_deg[ch];
932 317436 : azi_rad = speaker_node_azi_deg[ch] * PI_OVER_180;
933 317436 : if ( ( speaker_node_ele_deg[ch] >= -5 ) && ( speaker_node_ele_deg[ch] <= 5 ) )
934 : {
935 182698 : speaker_node_data[ch].ele_deg = 0.0f;
936 182698 : ele_rad = 0.0f;
937 182698 : speaker_node_data[ch].group = SPEAKER_NODE_HORIZONTAL;
938 182698 : num_horiz++;
939 : }
940 : else
941 : {
942 134738 : speaker_node_data[ch].ele_deg = speaker_node_ele_deg[ch];
943 134738 : ele_rad = speaker_node_ele_deg[ch] * PI_OVER_180;
944 134738 : if ( ele_rad < 0 )
945 : {
946 31184 : speaker_node_data[ch].group = SPEAKER_NODE_BOTTOM_HALF;
947 : }
948 : else
949 : {
950 103554 : speaker_node_data[ch].group = SPEAKER_NODE_TOP_HALF;
951 : }
952 : }
953 :
954 317436 : speaker_node_data[ch].unit_vec[0] = cosf( azi_rad ) * cosf( ele_rad );
955 317436 : speaker_node_data[ch].unit_vec[1] = sinf( azi_rad ) * cosf( ele_rad );
956 317436 : speaker_node_data[ch].unit_vec[2] = sinf( ele_rad );
957 : }
958 :
959 : /* Check for largest horizontal gap if there are at least 3 horizontal speaker nodes */
960 31140 : if ( num_horiz >= 3 )
961 : {
962 : int16_t i;
963 : uint16_t horiz_azi[VBAP_MAX_NUM_SPEAKER_NODES];
964 : uint16_t largest_gap;
965 : uint16_t temp;
966 :
967 31140 : i = 0;
968 213866 : for ( ch = 0; ch < num_speaker_nodes && i < num_horiz; ch++ )
969 : {
970 182726 : if ( speaker_node_data[ch].group == SPEAKER_NODE_HORIZONTAL )
971 : {
972 182698 : horiz_azi[i] = speaker_node_azi_deg[ch] < 0.0f ? (uint16_t) floorf( speaker_node_azi_deg[ch] + 360.0f ) : (uint16_t) floorf( speaker_node_azi_deg[ch] );
973 182698 : i++;
974 : }
975 : }
976 :
977 : /* Reorder horizontal azi to increasing order */
978 31140 : sort( horiz_azi, num_horiz );
979 :
980 : /* Find largest gap. Initialize with the wrap over gap. */
981 31140 : largest_gap = horiz_azi[0] - horiz_azi[num_horiz - 1] + 360;
982 182698 : for ( ch = 0; ch < num_horiz - 1; ch++ )
983 : {
984 151558 : temp = horiz_azi[ch + 1] - horiz_azi[ch];
985 151558 : if ( temp > largest_gap )
986 : {
987 60786 : largest_gap = temp;
988 : }
989 : }
990 :
991 : /* If largest gap is small enough, we have definitive zero elevation plane.
992 : * Otherwise, we should assign all speaker nodes to one group. */
993 31140 : if ( largest_gap <= VBAP_MAX_HORIZONTAL_GAP )
994 : {
995 31140 : in_all_mode = FALSE;
996 : }
997 : }
998 :
999 : /* Designate all speaker nodes to same group if there was no definitive zero
1000 : * elevation plane. */
1001 31140 : if ( in_all_mode )
1002 : {
1003 0 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
1004 : {
1005 0 : speaker_node_data[ch].group = SPEAKER_NODE_ALL;
1006 : }
1007 : }
1008 :
1009 31140 : return;
1010 : }
1011 :
1012 : /*-------------------------------------------------------------------------*
1013 : * matrix_inverse_3x3()
1014 : *
1015 : * 3-by-3 matrix inverse
1016 : *-------------------------------------------------------------------------*/
1017 :
1018 510312 : static void matrix_inverse_3x3(
1019 : const float **input_matrix, /* i : input matrix */
1020 : float inverse_matrix[3][3] /* o : output matrix */
1021 : )
1022 : {
1023 : int16_t k;
1024 : float determinant;
1025 : float cross_vec[3];
1026 :
1027 510312 : vbap_crossp( input_matrix[1], input_matrix[2], cross_vec );
1028 :
1029 510312 : determinant = dotp( input_matrix[0], cross_vec, 3 );
1030 :
1031 2041248 : for ( k = 0; k < 3; k++ )
1032 : {
1033 1530936 : inverse_matrix[k][0] = cross_vec[k] / determinant;
1034 : }
1035 :
1036 510312 : vbap_crossp( input_matrix[2], input_matrix[0], cross_vec );
1037 :
1038 2041248 : for ( k = 0; k < 3; k++ )
1039 : {
1040 1530936 : inverse_matrix[k][1] = cross_vec[k] / determinant;
1041 : }
1042 :
1043 510312 : vbap_crossp( input_matrix[0], input_matrix[1], cross_vec );
1044 :
1045 2041248 : for ( k = 0; k < 3; k++ )
1046 : {
1047 1530936 : inverse_matrix[k][2] = cross_vec[k] / determinant;
1048 : }
1049 :
1050 510312 : return;
1051 : }
1052 :
1053 : /*-------------------------------------------------------------------------*
1054 : * check_and_store_triplet()
1055 : *
1056 : * Check if the given loudspeaker triplet is a valid one and store data when
1057 : * valid triplet is found.
1058 : *-------------------------------------------------------------------------*/
1059 :
1060 510312 : static int16_t check_and_store_triplet(
1061 : const int16_t chA, /* i : first channel index that forms the loudspeaker triplet */
1062 : const int16_t chB, /* i : second channel index that forms the loudspeaker triplet */
1063 : const int16_t chC, /* i : third channel index that forms the loudspeaker triplet */
1064 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1065 : const VBAP_SPEAKER_NODE *speaker_node_data, /* i : speaker node data structure */
1066 : VBAP_VS_TRIPLET *triplets, /* o : vector of virtual surface triplets */
1067 : int16_t *triplet_index, /* i/o: index for the next free triplet slot */
1068 : float *triplet_azidegs, /* o : center azimuths of the found triplets */
1069 : int16_t *triplet_order /* o : initial order of triplet indices */
1070 : )
1071 : {
1072 : int16_t ch_check;
1073 : int16_t k;
1074 : int16_t speaker_node_found_inside_triplet;
1075 : uint8_t triplet_ok;
1076 : float inverse_matrix[3][3], unnormalized_gains[3];
1077 : const float *speaker_node_triplet_unit_vec_matrix[3];
1078 :
1079 : /* Triplet found, determine inverse matrix for VBAP formulation */
1080 510312 : speaker_node_triplet_unit_vec_matrix[0] = speaker_node_data[chA].unit_vec;
1081 510312 : speaker_node_triplet_unit_vec_matrix[1] = speaker_node_data[chB].unit_vec;
1082 510312 : speaker_node_triplet_unit_vec_matrix[2] = speaker_node_data[chC].unit_vec;
1083 510312 : matrix_inverse_3x3( speaker_node_triplet_unit_vec_matrix, inverse_matrix );
1084 :
1085 : /* Check through all speaker nodes that none of them are within the triplet.
1086 : * Node within the triplet is identified by that all three panning gains are positive.
1087 : * Epsilon-condition is for some small rounding issues.*/
1088 510312 : speaker_node_found_inside_triplet = 0;
1089 5926708 : for ( ch_check = 0; ch_check < num_speaker_nodes; ch_check++ )
1090 : {
1091 5416396 : if ( ( ch_check != chA ) && ( ch_check != chB ) && ( ch_check != chC ) )
1092 : {
1093 3885460 : triplet_ok = vector_matrix_multiply_3x3( speaker_node_data[ch_check].unit_vec, inverse_matrix, unnormalized_gains );
1094 3885460 : if ( triplet_ok && unnormalized_gains[0] > VBAP_EPSILON && unnormalized_gains[1] > VBAP_EPSILON && unnormalized_gains[2] > VBAP_EPSILON )
1095 : {
1096 0 : speaker_node_found_inside_triplet = 1;
1097 0 : break;
1098 : }
1099 : }
1100 : }
1101 :
1102 : /* No speaker node inside triplet -> appropriate triplet found, save data. */
1103 510312 : if ( speaker_node_found_inside_triplet == 0 )
1104 : {
1105 510312 : triplets[*triplet_index].speaker_node[0] = (uint8_t) chA;
1106 510312 : triplets[*triplet_index].speaker_node[1] = (uint8_t) chB;
1107 510312 : triplets[*triplet_index].speaker_node[2] = (uint8_t) chC;
1108 2041248 : for ( k = 0; k < 3; k++ )
1109 : {
1110 1530936 : mvr2r( inverse_matrix[k], triplets[*triplet_index].inverse_matrix[k], 3 );
1111 : }
1112 :
1113 : /* Get center azimuth for fast search use */
1114 510312 : triplet_azidegs[*triplet_index] = atan2f( speaker_node_data[chA].unit_vec[1] + speaker_node_data[chB].unit_vec[1] + speaker_node_data[chC].unit_vec[1],
1115 510312 : speaker_node_data[chA].unit_vec[0] + speaker_node_data[chB].unit_vec[0] + speaker_node_data[chC].unit_vec[0] ) *
1116 : _180_OVER_PI;
1117 :
1118 : /* Store increasing order indices for the later sorting step. */
1119 510312 : triplet_order[*triplet_index] = *triplet_index;
1120 :
1121 510312 : ( *triplet_index )++;
1122 :
1123 510312 : return 1;
1124 : }
1125 :
1126 : /* Triplet was not good */
1127 0 : return 0;
1128 : }
1129 :
1130 :
1131 : /*-------------------------------------------------------------------------*
1132 : * determine_virtual_surface_triplets()
1133 : *
1134 : * Determine virtual surface triples that are used for panning. This
1135 : * function is optimized for the use in cases where speaker nodes are in
1136 : * one group or divided into two separate groups divided by a horizontal
1137 : * layer.
1138 : *-------------------------------------------------------------------------*/
1139 :
1140 : /*! r: number of virtual surface triplets */
1141 62280 : static int16_t determine_virtual_surface_triplets(
1142 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1143 : const VBAP_SPEAKER_NODE *speaker_node_data, /* i : speaker node data structure */
1144 : int16_t connections[][2], /* i : vector of all connections */
1145 : const int16_t max_num_connections, /* i : max number of connections */
1146 : VBAP_VS_TRIPLET *triplets, /* o : vector of virtual surface triplets */
1147 : int16_t initial_search_indices[VBAP_NUM_SEARCH_SECTORS], /* o : initial search indices for this set of triplets corresponding to the search struct */
1148 : enum SpeakerNodeGroup allowed_group /* i : group of allowed speaker nodes for forming the triplets in this call */
1149 : )
1150 : {
1151 : int16_t chA, chB, chC, k, l, m;
1152 62280 : int16_t num_triplets = 0;
1153 : int16_t num_connected_to_chA;
1154 : int16_t connected_to_chA[VBAP_MAX_NUM_SPEAKER_NODES];
1155 : int16_t connection_uses_left[VBAP_MAX_NUM_SPEAKER_NODES];
1156 : float triplet_azidegs[VBAP_MAX_NUM_TRIPLETS];
1157 : int16_t triplet_order[VBAP_MAX_NUM_TRIPLETS];
1158 :
1159 : /* Each connection can be used exactly by two different virtual surface triplets. */
1160 62280 : set_s( connection_uses_left, 2, VBAP_MAX_NUM_SPEAKER_NODES );
1161 :
1162 697152 : for ( chA = 0; chA < num_speaker_nodes; chA++ )
1163 : {
1164 : /* Early skip if not in correct group. */
1165 634872 : if ( speaker_node_data[chA].group != allowed_group )
1166 : {
1167 500134 : continue;
1168 : }
1169 :
1170 : /* Get all connections connected to current chA that have not been used by
1171 : * two triplets yet. */
1172 134738 : num_connected_to_chA = 0;
1173 2706184 : for ( k = 0; k < max_num_connections; k++ )
1174 : {
1175 2571446 : if ( ( connections[k][0] == chA || connections[k][1] == chA ) && connection_uses_left[k] > 0 )
1176 : {
1177 560604 : connected_to_chA[num_connected_to_chA] = k;
1178 560604 : num_connected_to_chA++;
1179 : }
1180 : }
1181 :
1182 : /* Check that we have enough connections to use. We need at least two available connections to form a triplet.
1183 : * This can fail in later stages when all connections are already used. */
1184 134738 : if ( num_connected_to_chA < 2 )
1185 : {
1186 22166 : continue;
1187 : }
1188 :
1189 : /* Try to form triplets from each valid connection. */
1190 673176 : for ( k = 0; k < num_connected_to_chA; k++ )
1191 : {
1192 560604 : int16_t connect_index_k = connected_to_chA[k];
1193 560604 : chB = connections[connect_index_k][0] == chA ? connections[connect_index_k][1] : connections[connect_index_k][0];
1194 980876 : for ( l = k + 1; l < num_connected_to_chA; l++ )
1195 : {
1196 868304 : int16_t connect_index_l = connected_to_chA[l];
1197 868304 : chC = connections[connect_index_l][0] == chA ? connections[connect_index_l][1] : connections[connect_index_l][0];
1198 :
1199 : /* With chA, chB, and chC selected, we still need to find connection between chB and chC and verify that the triplet is valid */
1200 10613992 : for ( m = 0; m < max_num_connections; m++ )
1201 : {
1202 10256000 : if ( ( connections[m][0] == chB && connections[m][1] == chC ) || ( connections[m][1] == chB && connections[m][0] == chC ) )
1203 : {
1204 510312 : if ( check_and_store_triplet( chA, chB, chC, num_speaker_nodes, speaker_node_data, triplets, &num_triplets, triplet_azidegs, triplet_order ) )
1205 : {
1206 510312 : connection_uses_left[connect_index_k]--;
1207 510312 : connection_uses_left[connect_index_l]--;
1208 510312 : connection_uses_left[m]--;
1209 510312 : break;
1210 : }
1211 : }
1212 : }
1213 :
1214 : /* Check if chA-chB connection has been used in two triplets already. If yes, then break out of loop
1215 : * as this connection cannot be used for more triplets and we need to continue with another chA-chB
1216 : * connection. */
1217 868304 : if ( connection_uses_left[connect_index_k] < 1 )
1218 : {
1219 448032 : break;
1220 : }
1221 : }
1222 : }
1223 : }
1224 :
1225 : /* All triplets should be stored now. Sort them for search use and then determine the initial search indices for
1226 : * each search sector for this search struct. */
1227 62280 : v_sort_ind( triplet_azidegs, triplet_order, num_triplets );
1228 62280 : reorder_triplets( triplets, triplet_order, num_triplets );
1229 62280 : determine_initial_search_indices( num_triplets, triplet_azidegs, initial_search_indices );
1230 :
1231 62280 : return num_triplets;
1232 : }
1233 :
1234 : /*-------------------------------------------------------------------------*
1235 : * determine_initial_search_indices()
1236 : *
1237 : * Determine initial search indices used for fast search of correct triangle
1238 : *-------------------------------------------------------------------------*/
1239 :
1240 62280 : static void determine_initial_search_indices(
1241 : const int16_t num_triplets, /* i : number of triplets */
1242 : const float triplet_azidegs[VBAP_MAX_NUM_TRIPLETS], /* i : azimuths of triplets (in degrees) */
1243 : int16_t initial_search_indices[VBAP_NUM_SEARCH_SECTORS] /* o : initial search indices */
1244 : )
1245 : {
1246 : int16_t i, j;
1247 : float sector_reference_azideg;
1248 : float sector_border_start_azideg;
1249 : float sector_border_end_azideg;
1250 : int16_t best_index;
1251 : float min_azideg_diff;
1252 : float azideg_diff;
1253 :
1254 311400 : for ( i = 0; i < VBAP_NUM_SEARCH_SECTORS; i++ )
1255 : {
1256 249120 : sector_border_start_azideg = i * VBAP_SEARCH_SECTOR_SIZE;
1257 249120 : sector_border_end_azideg = ( i + 1 ) * VBAP_SEARCH_SECTOR_SIZE;
1258 249120 : sector_reference_azideg = ( sector_border_start_azideg + sector_border_end_azideg ) / 2.0f;
1259 249120 : best_index = 0;
1260 249120 : min_azideg_diff = 9999.9f;
1261 :
1262 2290368 : for ( j = 0; j < num_triplets; j++ )
1263 : {
1264 2041248 : azideg_diff = sector_reference_azideg - triplet_azidegs[j];
1265 2041248 : if ( azideg_diff > 180.0f )
1266 : {
1267 888832 : azideg_diff -= 360.0f;
1268 : }
1269 1152416 : else if ( azideg_diff < -180.0f )
1270 : {
1271 0 : azideg_diff += 360.0f;
1272 : }
1273 2041248 : azideg_diff = fabsf( azideg_diff );
1274 :
1275 2041248 : if ( azideg_diff < min_azideg_diff )
1276 : {
1277 846548 : min_azideg_diff = azideg_diff;
1278 846548 : best_index = j;
1279 : }
1280 : }
1281 :
1282 249120 : initial_search_indices[i] = best_index;
1283 : }
1284 :
1285 62280 : return;
1286 : }
1287 :
1288 :
1289 : /*-------------------------------------------------------------------------*
1290 : * determine_connections()
1291 : *
1292 : * Determine all valid connections between all speaker nodes
1293 : *-------------------------------------------------------------------------*/
1294 :
1295 31140 : static ivas_error determine_connections(
1296 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1297 : const VBAP_SPEAKER_NODE *speaker_node_data, /* i : speaker node data */
1298 : int16_t connections[][2], /* o : vector of connections */
1299 : const int16_t max_num_connections, /* i : max number of connections */
1300 : int16_t *group1_count, /* o : number of connections in first group */
1301 : int16_t *group2_start, /* o : start of second group of connections */
1302 : int16_t *group2_count /* o : number of connections in second group */
1303 : )
1304 : {
1305 : int16_t num_non_crossing_planes;
1306 : int16_t c;
1307 31140 : int16_t connection_write_index = 0;
1308 : float non_crossing_plane_elevation_deg[VBAP_MAX_PLANES];
1309 : ivas_error error;
1310 :
1311 31140 : set_f( non_crossing_plane_elevation_deg, 0.0f, VBAP_MAX_PLANES );
1312 :
1313 : /* Reset connection data */
1314 796608 : for ( c = 0; c < max_num_connections; c++ )
1315 : {
1316 765468 : connections[c][0] = VBAP_NOT_VALID_CONNECTION;
1317 : }
1318 :
1319 : /* This function determines some prominent elevated planes, that are favoured in making node-node connections. */
1320 31140 : num_non_crossing_planes = determine_non_crossing_planes( num_speaker_nodes, speaker_node_data, non_crossing_plane_elevation_deg );
1321 :
1322 : /* Process in different mode based on the grouping. It is enough to check for first node. */
1323 31140 : if ( speaker_node_data[0].group == SPEAKER_NODE_ALL )
1324 : {
1325 0 : if ( ( error = formulate_half_sphere_connections( speaker_node_data, num_speaker_nodes, SPEAKER_NODE_ALL, connections, &connection_write_index, max_num_connections, num_non_crossing_planes, non_crossing_plane_elevation_deg ) ) != IVAS_ERR_OK )
1326 : {
1327 0 : return error;
1328 : }
1329 : }
1330 : else
1331 : {
1332 : /* The node-node connections are determined in three stages: bottom, horizontal, and top. */
1333 31140 : if ( ( error = formulate_half_sphere_connections( speaker_node_data, num_speaker_nodes, SPEAKER_NODE_BOTTOM_HALF, connections, &connection_write_index, max_num_connections, num_non_crossing_planes, non_crossing_plane_elevation_deg ) ) != IVAS_ERR_OK )
1334 : {
1335 0 : return error;
1336 : }
1337 31140 : *group2_start = connection_write_index;
1338 :
1339 31140 : formulate_horizontal_connections( speaker_node_data, num_speaker_nodes, connections, &connection_write_index );
1340 31140 : *group1_count = connection_write_index;
1341 :
1342 31140 : if ( ( error = formulate_half_sphere_connections( speaker_node_data, num_speaker_nodes, SPEAKER_NODE_TOP_HALF, connections, &connection_write_index, max_num_connections, num_non_crossing_planes, non_crossing_plane_elevation_deg ) ) != IVAS_ERR_OK )
1343 : {
1344 0 : return error;
1345 : }
1346 31140 : *group2_count = connection_write_index - *group2_start;
1347 : }
1348 :
1349 31140 : return IVAS_ERR_OK;
1350 : }
1351 :
1352 :
1353 : /*-------------------------------------------------------------------------*
1354 : * determine_connection_class()
1355 : *
1356 : * Determine the type of connection
1357 : *-------------------------------------------------------------------------*/
1358 :
1359 : /*! r: type of connection */
1360 949397 : static enum ConnectionClass determine_connection_class(
1361 : const VBAP_SPEAKER_NODE *node_data, /* i : speaker node data */
1362 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1363 : const enum SpeakerNodeGroup group, /* i : speaker node group */
1364 : const int16_t chA, /* i : speaker node counter 1 */
1365 : const int16_t chB /* i : speaker node counter 2 */
1366 : )
1367 : {
1368 : int16_t ch, k;
1369 : const float *p1, *v2;
1370 : float v1v1, v1v2, v2v2, v1p1, v2p1;
1371 : float determinant;
1372 : float norm_distance_on_v1;
1373 : float vec_diff[3];
1374 : float v1[3];
1375 : float vTarget[3];
1376 : float energy_sum;
1377 : float eq_value;
1378 : float uvecdot;
1379 :
1380 : /* Check if connection passes through origin. This is not desired.
1381 : * When this happens, unit vectors point in opposite directions. */
1382 949397 : uvecdot = dotp( node_data[chA].unit_vec, node_data[chB].unit_vec, 3 ) + 1.0f;
1383 949397 : if ( uvecdot < VBAP_EPSILON && uvecdot > -VBAP_EPSILON )
1384 : {
1385 0 : return CONNECTION_WITH_SPEAKER_NODE_BEHIND;
1386 : }
1387 :
1388 : /* This loop checks the chA-chB connection with respect to all loudspeakers:
1389 : - in case there is a node behind or nearly behind the connection line. These
1390 : connections need to be discarded.
1391 : - in case there is a node that is closer to the connection line than 1/5 of the
1392 : connection length AND at the same horizontal plane. These connections need to be
1393 : weighted with a penalty (a special case, for example avoiding elevated L,R,C triplet)
1394 : */
1395 11032297 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
1396 : {
1397 : /* Select speaker_node only within TOP or BOTTOM sphere half, not being part of chA-chB pair */
1398 10185494 : if ( ( group == node_data[ch].group ) && ( ch != chA ) && ( ch != chB ) )
1399 : {
1400 : /* The following lines formulate the point on the chA-chB-connection that is
1401 : nearest to the origo-ch-line */
1402 2259028 : p1 = node_data[chA].unit_vec;
1403 9036112 : for ( k = 0; k < 3; k++ )
1404 : {
1405 6777084 : v1[k] = node_data[chB].unit_vec[k] - node_data[chA].unit_vec[k];
1406 : }
1407 2259028 : v2 = node_data[ch].unit_vec;
1408 2259028 : v1v1 = dotp( v1, v1, 3 );
1409 2259028 : v1v2 = dotp( v1, v2, 3 );
1410 2259028 : v2v2 = 1.0f; /* dotp(v2, v2) is always 1. */
1411 2259028 : v1p1 = dotp( v1, p1, 3 );
1412 2259028 : v2p1 = dotp( v2, p1, 3 );
1413 2259028 : determinant = ( v1v1 ) * ( -v2v2 ) + ( v1v2 * v1v2 );
1414 : /* Norm distance = distance parameter on line chA-chB, determines point that is
1415 : nearest to origo-ch line. Distance 0 means chA and distance 1 means chB. Can be
1416 : outside this region as well.*/
1417 2259028 : norm_distance_on_v1 = ( 1.0f / determinant ) * ( ( v2v2 ) * ( v1p1 ) + ( -v1v2 ) * ( v2p1 ) );
1418 :
1419 : /* Continue only if the nearest point is between chA and chB */
1420 2259028 : if ( norm_distance_on_v1 > 0.0f && norm_distance_on_v1 < 1.0f )
1421 : {
1422 : /* Formulate vTarget, that is an unit vector that goes through the determined point on chA-chB connection */
1423 1855706 : energy_sum = 0.0f;
1424 7422824 : for ( k = 0; k < 3; k++ )
1425 : {
1426 5567118 : vTarget[k] = p1[k] + norm_distance_on_v1 * v1[k];
1427 5567118 : energy_sum += vTarget[k] * vTarget[k];
1428 5567118 : vec_diff[k] = vTarget[k] - v2[k];
1429 : }
1430 1855706 : eq_value = sqrtf( 1.0f / energy_sum );
1431 7422824 : for ( k = 0; k < 3; k++ )
1432 : {
1433 5567118 : vTarget[k] *= eq_value;
1434 : }
1435 :
1436 : /* A check if the angle between vTarget and node_data[ch].unit_vec is less than 1 degree.
1437 : Essentially reveals if there is a speaker node too close "behind" the connection. Such
1438 : connections should be rejected.*/
1439 1855706 : if ( dotp( vTarget, v2, 3 ) > 0.9998f )
1440 : {
1441 102594 : return CONNECTION_WITH_SPEAKER_NODE_BEHIND;
1442 : }
1443 :
1444 : /* A special case, mainly accounting for ELEVATED L,R,C speaker nodes.
1445 : A triplet between these nodes is not desired if there is a top node,
1446 : a penalty is implemented to take care of this. */
1447 1753112 : if ( sqrtf( v1v1 ) > 5.0f * sqrtf( dotp( vec_diff, vec_diff, 3 ) ) )
1448 : {
1449 0 : if ( fabsf( node_data[chB].unit_vec[2] - node_data[chA].unit_vec[2] ) < VBAP_EPSILON )
1450 : {
1451 0 : return ELEVATED_PLANE_THIN_TRIANGLE_CONNECTION;
1452 : }
1453 : }
1454 : }
1455 : }
1456 : }
1457 :
1458 846803 : return REGULAR_CONNECTION;
1459 : }
1460 :
1461 :
1462 : /*-------------------------------------------------------------------------*
1463 : * formulate_horizontal_connections()
1464 : *
1465 : * Formulate connections in the horizontal plane
1466 : *-------------------------------------------------------------------------*/
1467 :
1468 31140 : static void formulate_horizontal_connections(
1469 : const VBAP_SPEAKER_NODE *speaker_node_data, /* i : speaker node data */
1470 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1471 : int16_t connections[][2], /* o : vector of all connections */
1472 : int16_t *connection_write_index )
1473 : {
1474 : int16_t ch;
1475 : int16_t chCheck;
1476 : int16_t next_index;
1477 : float min_arc_diff;
1478 : float arc_diff;
1479 :
1480 348576 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
1481 : {
1482 : /* Find next horizontal speaker node */
1483 317436 : if ( speaker_node_data[ch].group == SPEAKER_NODE_HORIZONTAL )
1484 : {
1485 182698 : next_index = -1;
1486 182698 : min_arc_diff = 9999.0f;
1487 2065946 : for ( chCheck = 0; chCheck < num_speaker_nodes; chCheck++ )
1488 : {
1489 1883248 : if ( ( ch != chCheck ) && ( speaker_node_data[chCheck].group == SPEAKER_NODE_HORIZONTAL ) )
1490 : {
1491 919316 : arc_diff = speaker_node_data[chCheck].azi_deg - speaker_node_data[ch].azi_deg;
1492 1378974 : while ( arc_diff < 0.0f )
1493 : {
1494 459658 : arc_diff += 360.0f;
1495 : }
1496 919316 : if ( arc_diff < min_arc_diff )
1497 : {
1498 415514 : min_arc_diff = arc_diff;
1499 415514 : next_index = chCheck;
1500 : }
1501 : }
1502 : }
1503 182698 : connections[*connection_write_index][0] = ch;
1504 182698 : connections[*connection_write_index][1] = next_index;
1505 182698 : ( *connection_write_index )++;
1506 : }
1507 : }
1508 :
1509 31140 : return;
1510 : }
1511 :
1512 : /*-------------------------------------------------------------------------*
1513 : * check_plane_crossing()
1514 : *
1515 : * Check crossing of non-allowed planes
1516 : *-------------------------------------------------------------------------*/
1517 :
1518 : /*! r: truth value for crossing */
1519 846803 : static int16_t check_plane_crossing(
1520 : const float ele1_deg, /* i : speaker node 1 elevation */
1521 : const float ele2_deg, /* i : speaker node 2 elevation */
1522 : const int16_t num_non_crossing_planes, /* i : number of non-crossing planes */
1523 : const float *non_crossing_plane_elevation_deg /* i : vector non-crossing plane elevations*/
1524 : )
1525 : {
1526 : /* Find if the connection crosses a non-crossing plane, with 1-degree threshold. */
1527 : int16_t plane;
1528 :
1529 1404537 : for ( plane = 0; plane < num_non_crossing_planes; plane++ )
1530 : {
1531 586072 : if ( ( ele1_deg > non_crossing_plane_elevation_deg[plane] + 1.0f ) && ( ele2_deg < non_crossing_plane_elevation_deg[plane] - 1.0f ) )
1532 : {
1533 64 : return 1;
1534 : }
1535 :
1536 586008 : if ( ( ele2_deg > non_crossing_plane_elevation_deg[plane] + 1.0f ) && ( ele1_deg < non_crossing_plane_elevation_deg[plane] - 1.0f ) )
1537 : {
1538 28274 : return 1;
1539 : }
1540 : }
1541 :
1542 818465 : return 0;
1543 : }
1544 :
1545 :
1546 : /*-------------------------------------------------------------------------*
1547 : * get_half_sphere_connection_options()
1548 : *
1549 : * Get list of all potential connections at the half-sphere
1550 : *-------------------------------------------------------------------------*/
1551 :
1552 62280 : static ivas_error get_half_sphere_connection_options(
1553 : const VBAP_SPEAKER_NODE *speaker_node_data, /* i : speaker node data */
1554 : const enum SpeakerNodeGroup group, /* i : speaker node group */
1555 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1556 : const int16_t num_non_crossing_planes, /* i : number of non-crossing planes */
1557 : const float *non_crossing_plane_elevation_deg, /* i : vector of non-crossing plane elevations */
1558 : ConnectionOption **connection_options_pr, /* o : list of connection options */
1559 : int16_t *num_connection_options /* o : number of connection options */
1560 : )
1561 : {
1562 62280 : int16_t max_num_connection_options = 0;
1563 62280 : int16_t index = 0;
1564 : int16_t node, chA, chB, c, c_cmp;
1565 : ConnectionOption *c_options, *c_options_reorder;
1566 :
1567 : /* Count max num connection options at the half sphere */
1568 697152 : for ( node = 0; node < num_speaker_nodes; node++ )
1569 : {
1570 634872 : if ( speaker_node_data[node].group == group || speaker_node_data[node].group == SPEAKER_NODE_HORIZONTAL )
1571 : {
1572 500134 : max_num_connection_options += index;
1573 500134 : index++;
1574 : }
1575 : }
1576 :
1577 : /* Init memory for connection options */
1578 62280 : if ( ( c_options = (ConnectionOption *) malloc( sizeof( ConnectionOption ) * max_num_connection_options ) ) == NULL )
1579 : {
1580 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
1581 : }
1582 :
1583 1930993 : for ( c = 0; c < max_num_connection_options; c++ )
1584 : {
1585 1868713 : c_options[c].chA = -1;
1586 1868713 : c_options[c].chB = -1;
1587 1868713 : c_options[c].arc = 1e20f;
1588 1868713 : c_options[c].arc_weighted = 1e20f;
1589 : }
1590 :
1591 : /* Determine connection options for the half-sphere */
1592 62280 : index = 0;
1593 697152 : for ( chA = 0; chA < num_speaker_nodes; chA++ )
1594 : {
1595 : /* First loudspeaker at the connection is at the half sphere */
1596 634872 : if ( speaker_node_data[chA].group == group || speaker_node_data[chA].group == SPEAKER_NODE_HORIZONTAL )
1597 : {
1598 3253779 : for ( chB = chA + 1; chB < num_speaker_nodes; chB++ )
1599 : {
1600 : /* Second loudspeaker at the connection is at the half sphere, but so that first and second are not both horizontal. */
1601 2753645 : if ( ( speaker_node_data[chB].group == group ) || ( speaker_node_data[chB].group == SPEAKER_NODE_HORIZONTAL && speaker_node_data[chA].group == group ) )
1602 : {
1603 949397 : int16_t ConnectionClass = determine_connection_class( speaker_node_data, num_speaker_nodes, group, chA, chB );
1604 : /* Connection is considered only if there is no speaker node behind it */
1605 949397 : if ( ConnectionClass != CONNECTION_WITH_SPEAKER_NODE_BEHIND )
1606 : {
1607 : /* Store connection information */
1608 846803 : c_options[index].chA = chA;
1609 846803 : c_options[index].chB = chB;
1610 846803 : c_options[index].arc = acosf( dotp( speaker_node_data[chA].unit_vec, speaker_node_data[chB].unit_vec, 3 ) );
1611 846803 : c_options[index].arc_weighted = c_options[index].arc;
1612 :
1613 : /* A special case, mainly accounting for ELEVATED L,R,C speaker nodes.
1614 : A triplet between these nodes is not desired if there is a top node,
1615 : a penalty is implemented to take care of this. */
1616 846803 : if ( ConnectionClass == ELEVATED_PLANE_THIN_TRIANGLE_CONNECTION )
1617 : {
1618 0 : c_options[index].arc_weighted *= 2.0f;
1619 : }
1620 :
1621 : /* If the connection passes a pre-determined plane of speaker nodes, then add further penalty */
1622 846803 : if ( check_plane_crossing( speaker_node_data[chA].ele_deg, speaker_node_data[chB].ele_deg, num_non_crossing_planes, non_crossing_plane_elevation_deg ) )
1623 : {
1624 28338 : c_options[index].arc_weighted *= 2.0f;
1625 : }
1626 846803 : index++;
1627 : }
1628 : }
1629 : }
1630 : }
1631 : }
1632 : /* Number of found connection options at the half sphere */
1633 62280 : *num_connection_options = index;
1634 :
1635 : /* Init memory for reordered connection options and order by arc_weighted,
1636 : * which informs of the preference order of the connections in case they cross */
1637 62280 : if ( ( c_options_reorder = (ConnectionOption *) malloc( sizeof( ConnectionOption ) * ( *num_connection_options ) ) ) == NULL )
1638 : {
1639 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for VBAP data\n" ) );
1640 : }
1641 :
1642 909083 : for ( c = 0; c < *num_connection_options; c++ )
1643 : {
1644 : float min_arc_weighted;
1645 : int16_t min_arc_index;
1646 846803 : min_arc_weighted = 1e20f;
1647 846803 : min_arc_index = -1;
1648 :
1649 20824382 : for ( c_cmp = 0; c_cmp < *num_connection_options; c_cmp++ )
1650 : {
1651 19977579 : if ( c_options[c_cmp].arc_weighted <= min_arc_weighted )
1652 : {
1653 4413923 : min_arc_weighted = c_options[c_cmp].arc_weighted;
1654 4413923 : min_arc_index = c_cmp;
1655 : }
1656 : }
1657 846803 : c_options_reorder[c].chA = c_options[min_arc_index].chA;
1658 846803 : c_options_reorder[c].chB = c_options[min_arc_index].chB;
1659 846803 : c_options_reorder[c].arc = c_options[min_arc_index].arc;
1660 846803 : c_options_reorder[c].arc_weighted = c_options[min_arc_index].arc_weighted;
1661 846803 : c_options[min_arc_index].arc_weighted = 1e20f;
1662 : }
1663 :
1664 : /* Set reordered connections as output and free temporary data */
1665 62280 : *connection_options_pr = c_options_reorder;
1666 62280 : free( c_options );
1667 :
1668 62280 : return IVAS_ERR_OK;
1669 : }
1670 :
1671 :
1672 : /*-------------------------------------------------------------------------*
1673 : * formulate_half_sphere_connections()
1674 : *
1675 : * Formulate half-sphere connections
1676 : *-------------------------------------------------------------------------*/
1677 :
1678 62280 : static ivas_error formulate_half_sphere_connections(
1679 : const VBAP_SPEAKER_NODE *speaker_node_data, /* i : speaker node data */
1680 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1681 : const enum SpeakerNodeGroup group, /* i : speaker node group */
1682 : int16_t connections[][2], /* o : vector of connections */
1683 : int16_t *connection_write_index,
1684 : const int16_t max_num_connections, /* i : max number of connections */
1685 : const int16_t num_non_crossing_planes, /* i : number of non-crossing planes */
1686 : const float *non_crossing_plane_elevation_deg /* i : vector of non-crossing plane elevations */
1687 : )
1688 : {
1689 : /* Variable initializations */
1690 : int16_t c, chA, chB, cmp_chA, cmp_chB, k, c_opt;
1691 : int16_t new_connection_is_valid;
1692 : int16_t within_first_arc;
1693 : float planeCrossingVec[3];
1694 : float new_arc;
1695 : float new_cross[3];
1696 : float tmpFloat;
1697 : float cmp_arc;
1698 : float normVal;
1699 : float angleCmp;
1700 : float connection_arc[( VBAP_MAX_NUM_SPEAKER_NODES - 2 ) * 3];
1701 : float connection_cross[( VBAP_MAX_NUM_SPEAKER_NODES - 2 ) * 3][3];
1702 : ConnectionOption *connection_options;
1703 : int16_t num_connection_options;
1704 : int16_t half_sphere_first_connection;
1705 : ivas_error error;
1706 :
1707 62280 : half_sphere_first_connection = *connection_write_index;
1708 :
1709 : /* Obtain all connection options (i.e., channel pairs) at the half sphere. The function orders them
1710 : * in terms of which connection to keep if two connections would cross each other. */
1711 62280 : if ( ( error = get_half_sphere_connection_options(
1712 : speaker_node_data,
1713 : group,
1714 : num_speaker_nodes,
1715 : num_non_crossing_planes,
1716 : non_crossing_plane_elevation_deg,
1717 : &connection_options,
1718 : &num_connection_options ) ) != IVAS_ERR_OK )
1719 : {
1720 0 : return error;
1721 : }
1722 :
1723 62280 : set_f( connection_arc, 0.0f, max_num_connections );
1724 1593216 : for ( c = 0; c < max_num_connections; c++ )
1725 : {
1726 1530936 : set_f( connection_cross[c], 0.0f, 3 );
1727 : }
1728 :
1729 : /* The following loop goes through all reasonable chA - chB pairs for the half-sphere */
1730 62280 : c_opt = 0;
1731 729538 : while ( c_opt < num_connection_options && *connection_write_index < max_num_connections )
1732 : {
1733 667258 : chA = connection_options[c_opt].chA;
1734 667258 : chB = connection_options[c_opt].chB;
1735 667258 : new_arc = connection_options[c_opt].arc;
1736 :
1737 : /* Cross-product is needed for later stages */
1738 667258 : vbap_crossp( speaker_node_data[chA].unit_vec, speaker_node_data[chB].unit_vec, new_cross );
1739 :
1740 : /* Determine if new connection between chA and chB is valid */
1741 667258 : new_connection_is_valid = 1;
1742 667258 : c = half_sphere_first_connection;
1743 4542128 : while ( ( c < *connection_write_index ) && new_connection_is_valid )
1744 : {
1745 3874870 : cmp_chA = connections[c][0];
1746 3874870 : cmp_chB = connections[c][1];
1747 : /* The connections are compared only if they don't involve same speaker nodes */
1748 3874870 : if ( ( cmp_chA != chA ) && ( cmp_chA != chB ) && ( cmp_chB != chA ) && ( cmp_chB != chB ) )
1749 : {
1750 : /* The following lines determine if the connection chA-chB crosses with the connection cmp_chA-cmp_chB.*/
1751 : /* The connections, i.e., node-pairs determine a plane. The crossing can be determined by
1752 : * studying the intersection of these planes. */
1753 1894676 : vbap_crossp( connection_cross[c], new_cross, planeCrossingVec );
1754 1894676 : tmpFloat = 1e-12f;
1755 1894676 : cmp_arc = connection_arc[c];
1756 7578704 : for ( k = 0; k < 3; k++ )
1757 : {
1758 5684028 : tmpFloat += planeCrossingVec[k] * planeCrossingVec[k];
1759 : }
1760 1894676 : normVal = sqrtf( 1.0f / tmpFloat );
1761 :
1762 7578704 : for ( k = 0; k < 3; k++ )
1763 : {
1764 5684028 : planeCrossingVec[k] *= normVal;
1765 : }
1766 :
1767 : /* If the plane intersection is between both connections, then the two connections cross. */
1768 : /* Study first if the crossing is between arc chA-chB */
1769 1894676 : angleCmp = acosf( dotp( planeCrossingVec, speaker_node_data[chA].unit_vec, 3 ) );
1770 1894676 : angleCmp += acosf( dotp( planeCrossingVec, speaker_node_data[chB].unit_vec, 3 ) );
1771 :
1772 1894676 : within_first_arc = 0;
1773 1894676 : if ( fabsf( new_arc - angleCmp ) < 0.01f )
1774 : {
1775 408527 : within_first_arc = 1;
1776 : }
1777 1486149 : else if ( fabsf( new_arc - ( 2.0f * EVS_PI - angleCmp ) ) < 0.01f )
1778 : {
1779 428893 : within_first_arc = 1;
1780 : /* In this case, the plane crossing vector is inverted. The inverse is another
1781 : * plane-crossing vector, and detected to be between chA-chB connection.*/
1782 1715572 : for ( k = 0; k < 3; k++ )
1783 : {
1784 1286679 : planeCrossingVec[k] *= -1.0f;
1785 : }
1786 : }
1787 :
1788 : /* Study if the crossing is also between arc cmp_chA-cmp_chB */
1789 1894676 : if ( within_first_arc > 0 )
1790 : {
1791 837420 : angleCmp = acosf( dotp( planeCrossingVec, speaker_node_data[cmp_chA].unit_vec, 3 ) );
1792 837420 : angleCmp += acosf( dotp( planeCrossingVec, speaker_node_data[cmp_chB].unit_vec, 3 ) );
1793 837420 : if ( fabsf( cmp_arc - angleCmp ) < 0.01f )
1794 : {
1795 : /* A crossing is detected. The new connection is not valid, because
1796 : * the connections were ordered in order of preference (arc_weighted) */
1797 84488 : new_connection_is_valid = 0;
1798 : }
1799 : }
1800 : }
1801 3874870 : c++;
1802 : }
1803 :
1804 : /* Store the new connection which has been confirmed valid */
1805 667258 : if ( new_connection_is_valid > 0 )
1806 : {
1807 582770 : connections[*connection_write_index][0] = chA;
1808 582770 : connections[*connection_write_index][1] = chB;
1809 582770 : connection_arc[*connection_write_index] = new_arc;
1810 582770 : connection_cross[*connection_write_index][0] = new_cross[0];
1811 582770 : connection_cross[*connection_write_index][1] = new_cross[1];
1812 582770 : connection_cross[*connection_write_index][2] = new_cross[2];
1813 582770 : ( *connection_write_index )++;
1814 : }
1815 667258 : c_opt++;
1816 : }
1817 :
1818 62280 : free( connection_options );
1819 :
1820 62280 : return IVAS_ERR_OK;
1821 : }
1822 :
1823 :
1824 : /*-------------------------------------------------------------------------*
1825 : * determine_non_crossing_planes()
1826 : *
1827 : * Determine non-crossing planes
1828 : *-------------------------------------------------------------------------*/
1829 :
1830 : /*! r: number of non-crossing planes */
1831 31140 : static int16_t determine_non_crossing_planes(
1832 : const int16_t num_speaker_nodes, /* i : number of speaker nodes */
1833 : const VBAP_SPEAKER_NODE *node_data, /* i : speaker node data */
1834 : float *non_crossing_plane_elevation_deg /* o : vector of non-crossing plane elevations */
1835 : )
1836 : {
1837 : float next_ele_check;
1838 : float ele_check;
1839 : float max_gap;
1840 : float gap_to_next_ls;
1841 : int16_t ch, ch_cmp;
1842 : int16_t num_planes;
1843 :
1844 31140 : ele_check = -999.9f;
1845 31140 : num_planes = 0;
1846 :
1847 : /* For each plane, check if a non-crossing plane should be determined */
1848 114779 : while ( ele_check < 90.0f )
1849 : {
1850 83643 : next_ele_check = 999.9f;
1851 : /* Find next node elevation that is not in horizontal plane */
1852 953992 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
1853 : {
1854 870349 : if ( ( node_data[ch].group != SPEAKER_NODE_HORIZONTAL ) && ( node_data[ch].ele_deg > ele_check + VBAP_EPSILON ) && ( node_data[ch].ele_deg < next_ele_check - VBAP_EPSILON ) )
1855 : {
1856 105006 : next_ele_check = node_data[ch].ele_deg;
1857 : }
1858 : }
1859 83643 : ele_check = next_ele_check;
1860 83643 : if ( ele_check > 90.0f )
1861 : {
1862 : /* When no next node elevation found, break loop */
1863 4 : break;
1864 : }
1865 :
1866 83639 : max_gap = -9999.9f;
1867 953940 : for ( ch = 0; ch < num_speaker_nodes; ch++ )
1868 : {
1869 : /* Find gap to the next speaker node at the same plane */
1870 870301 : if ( fabsf( node_data[ch].ele_deg - ele_check ) < VBAP_EPSILON )
1871 : {
1872 134738 : gap_to_next_ls = 99999.9f;
1873 1594560 : for ( ch_cmp = 0; ch_cmp < num_speaker_nodes; ch_cmp++ )
1874 : {
1875 1459822 : if ( ch_cmp != ch && fabsf( node_data[ch_cmp].ele_deg - ele_check ) < VBAP_EPSILON )
1876 : {
1877 191394 : float gap = node_data[ch_cmp].azi_deg - node_data[ch].azi_deg;
1878 287091 : while ( gap < 0 )
1879 : {
1880 95697 : gap += 360.0f;
1881 : }
1882 191394 : if ( gap < gap_to_next_ls )
1883 : {
1884 131930 : gap_to_next_ls = gap;
1885 : }
1886 : }
1887 : }
1888 : /* Find maximum gap on that plane */
1889 134738 : if ( gap_to_next_ls > max_gap )
1890 : {
1891 90967 : max_gap = gap_to_next_ls;
1892 : }
1893 : }
1894 : }
1895 :
1896 : /* If maximum gap is small enough, then a non-crossing plane is detected */
1897 83639 : if ( max_gap < ( VBAP_MAX_HORIZONTAL_GAP_FOR_PLANE_DETECTION + VBAP_EPSILON ) && max_gap > 0.0f )
1898 : {
1899 14866 : non_crossing_plane_elevation_deg[num_planes] = ele_check;
1900 14866 : num_planes++;
1901 14866 : if ( num_planes == VBAP_MAX_PLANES )
1902 : {
1903 : /* Memory init limit. Does not happen with any real speaker node configuration.
1904 : Triangulation succeeds even if number of non_crossing_planes are limited. */
1905 0 : break;
1906 : }
1907 : }
1908 : }
1909 :
1910 31140 : return num_planes;
1911 : }
1912 :
1913 :
1914 : /*-------------------------------------------------------------------------*
1915 : * reorder_triplets()
1916 : *
1917 : * Reorder virtual surface triplets into provided target order.
1918 : *-------------------------------------------------------------------------*/
1919 :
1920 62280 : static void reorder_triplets(
1921 : VBAP_VS_TRIPLET *triplets, /* i/o: VS triplets to be reordered */
1922 : const int16_t *target_order, /* i : Target order for VS triplets */
1923 : const int16_t num_triplets /* i : Number of VS triplets */
1924 : )
1925 : {
1926 : VBAP_VS_TRIPLET tempTriplets[VBAP_MAX_NUM_TRIPLETS];
1927 : int16_t c;
1928 :
1929 : /* First copy to temp array */
1930 572592 : for ( c = 0; c < num_triplets; c++ )
1931 : {
1932 510312 : tempTriplets[c] = triplets[c];
1933 : }
1934 :
1935 : /* Then move back in sorted order */
1936 572592 : for ( c = 0; c < num_triplets; c++ )
1937 : {
1938 510312 : triplets[c] = tempTriplets[target_order[c]];
1939 : }
1940 :
1941 62280 : return;
1942 : }
|
