ogl_beamforming

Ultrasound Beamforming Implemented with OpenGL
git clone anongit@rnpnr.xyz:ogl_beamforming.git
Log | Files | Refs | Feed | Submodules | README | LICENSE

math.c (20102B)


      1 #include "external/cephes.c"
      2 
      3 function void
      4 fill_kronecker_sub_matrix(f32 *out, i32 out_stride, f32 scale, f32 *b, iv2 b_dim)
      5 {
      6 	f32x4 vscale = dup_f32x4((f32)scale);
      7 	for (i32 i = 0; i < b_dim.y; i++) {
      8 		for (i32 j = 0; j < b_dim.x; j += 4, b += 4) {
      9 			f32x4 vb = load_f32x4(b);
     10 			store_f32x4(out + j, mul_f32x4(vscale, vb));
     11 		}
     12 		out += out_stride;
     13 	}
     14 }
     15 
     16 /* NOTE: this won't check for valid space/etc and assumes row major order */
     17 function void
     18 kronecker_product(f32 *out, f32 *a, iv2 a_dim, f32 *b, iv2 b_dim)
     19 {
     20 	iv2 out_dim = {{a_dim.x * b_dim.x, a_dim.y * b_dim.y}};
     21 	assert(out_dim.y % 4 == 0);
     22 	for (i32 i = 0; i < a_dim.y; i++) {
     23 		f32 *vout = out;
     24 		for (i32 j = 0; j < a_dim.x; j++, a++) {
     25 			fill_kronecker_sub_matrix(vout, out_dim.y, *a, b, b_dim);
     26 			vout += b_dim.y;
     27 		}
     28 		out += out_dim.y * b_dim.x;
     29 	}
     30 }
     31 
     32 /* NOTE/TODO: to support even more hadamard sizes use the Paley construction */
     33 function f32 *
     34 make_hadamard_transpose(Arena *a, i32 dim)
     35 {
     36 	read_only local_persist	f32 hadamard_12_12_transpose[] = {
     37 		1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
     38 		1, -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1,
     39 		1,  1, -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,
     40 		1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1,  1,
     41 		1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1,
     42 		1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,
     43 		1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,
     44 		1, -1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1,
     45 		1, -1, -1,  1,  1,  1, -1,  1, -1, -1,  1, -1,
     46 		1, -1, -1, -1,  1,  1,  1, -1,  1, -1, -1,  1,
     47 		1,  1, -1, -1, -1,  1,  1,  1, -1,  1, -1, -1,
     48 		1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1, -1,
     49 	};
     50 
     51 	read_only local_persist f32 hadamard_20_20_transpose[] = {
     52 		1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
     53 		1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1,
     54 		1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1,
     55 		1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,
     56 		1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,
     57 		1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1,
     58 		1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1,
     59 		1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1,
     60 		1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1,
     61 		1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,
     62 		1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1,
     63 		1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,
     64 		1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1,
     65 		1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,
     66 		1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,
     67 		1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,
     68 		1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,
     69 		1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1,
     70 		1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1,
     71 		1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,
     72 	};
     73 
     74 	f32 *result = 0;
     75 
     76 	b32 power_of_2     = ISPOWEROF2(dim);
     77 	b32 multiple_of_12 = dim % 12 == 0;
     78 	b32 multiple_of_20 = dim % 20 == 0;
     79 	iz elements        = dim * dim;
     80 
     81 	i32 base_dim = 0;
     82 	if (power_of_2) {
     83 		base_dim  = dim;
     84 	} else if (multiple_of_20 && ISPOWEROF2(dim / 20)) {
     85 		base_dim  = 20;
     86 		dim      /= 20;
     87 	} else if (multiple_of_12 && ISPOWEROF2(dim / 12)) {
     88 		base_dim  = 12;
     89 		dim      /= 12;
     90 	}
     91 
     92 	if (ISPOWEROF2(dim) && base_dim && arena_capacity(a, f32) >= elements * (1 + (dim != base_dim))) {
     93 		result = push_array(a, f32, elements);
     94 
     95 		Arena tmp = *a;
     96 		f32 *m = dim == base_dim ? result : push_array(&tmp, f32, elements);
     97 
     98 		#define IND(i, j) ((i) * dim + (j))
     99 		m[0] = 1;
    100 		for (i32 k = 1; k < dim; k *= 2) {
    101 			for (i32 i = 0; i < k; i++) {
    102 				for (i32 j = 0; j < k; j++) {
    103 					f32 val = m[IND(i, j)];
    104 					m[IND(i + k, j)]     =  val;
    105 					m[IND(i, j + k)]     =  val;
    106 					m[IND(i + k, j + k)] = -val;
    107 				}
    108 			}
    109 		}
    110 		#undef IND
    111 
    112 		f32 *m2 = 0;
    113 		iv2 m2_dim;
    114 		switch (base_dim) {
    115 		case 12:{ m2 = hadamard_12_12_transpose; m2_dim = (iv2){{12, 12}}; }break;
    116 		case 20:{ m2 = hadamard_20_20_transpose; m2_dim = (iv2){{20, 20}}; }break;
    117 		}
    118 		if (m2) kronecker_product(result, m, (iv2){{dim, dim}}, m2, m2_dim);
    119 	}
    120 
    121 	return result;
    122 }
    123 
    124 function b32
    125 u128_equal(u128 a, u128 b)
    126 {
    127 	b32 result = a.U64[0] == b.U64[0] && a.U64[1] == b.U64[1];
    128 	return result;
    129 }
    130 
    131 function b32
    132 iv2_equal(iv2 a, iv2 b)
    133 {
    134 	b32 result = a.x == b.x && a.y == b.y;
    135 	return result;
    136 }
    137 
    138 function b32
    139 iv3_equal(iv3 a, iv3 b)
    140 {
    141 	b32 result = a.x == b.x && a.y == b.y && a.z == b.z;
    142 	return result;
    143 }
    144 
    145 function v2
    146 clamp_v2_rect(v2 v, Rect r)
    147 {
    148 	v2 result = v;
    149 	result.x = CLAMP(v.x, r.pos.x, r.pos.x + r.size.x);
    150 	result.y = CLAMP(v.y, r.pos.y, r.pos.y + r.size.y);
    151 	return result;
    152 }
    153 
    154 function v2
    155 v2_scale(v2 a, f32 scale)
    156 {
    157 	v2 result;
    158 	result.x = a.x * scale;
    159 	result.y = a.y * scale;
    160 	return result;
    161 }
    162 
    163 function v2
    164 v2_add(v2 a, v2 b)
    165 {
    166 	v2 result;
    167 	result.x = a.x + b.x;
    168 	result.y = a.y + b.y;
    169 	return result;
    170 }
    171 
    172 function v2
    173 v2_sub(v2 a, v2 b)
    174 {
    175 	v2 result = v2_add(a, v2_scale(b, -1.0f));
    176 	return result;
    177 }
    178 
    179 function v2
    180 v2_mul(v2 a, v2 b)
    181 {
    182 	v2 result;
    183 	result.x = a.x * b.x;
    184 	result.y = a.y * b.y;
    185 	return result;
    186 }
    187 
    188 function v2
    189 v2_div(v2 a, v2 b)
    190 {
    191 	v2 result;
    192 	result.x = a.x / b.x;
    193 	result.y = a.y / b.y;
    194 	return result;
    195 }
    196 
    197 function v2
    198 v2_floor(v2 a)
    199 {
    200 	v2 result;
    201 	result.x = (f32)((i32)a.x);
    202 	result.y = (f32)((i32)a.y);
    203 	return result;
    204 }
    205 
    206 function f32
    207 v2_magnitude_squared(v2 a)
    208 {
    209 	f32 result = a.x * a.x + a.y * a.y;
    210 	return result;
    211 }
    212 
    213 function f32
    214 v2_magnitude(v2 a)
    215 {
    216 	f32 result = sqrt_f32(a.x * a.x + a.y * a.y);
    217 	return result;
    218 }
    219 
    220 function v3
    221 cross(v3 a, v3 b)
    222 {
    223 	v3 result;
    224 	result.x = a.y * b.z - a.z * b.y;
    225 	result.y = a.z * b.x - a.x * b.z;
    226 	result.z = a.x * b.y - a.y * b.x;
    227 	return result;
    228 }
    229 
    230 function v3
    231 v3_from_iv3(iv3 v)
    232 {
    233 	v3 result;
    234 	result.E[0] = (f32)v.E[0];
    235 	result.E[1] = (f32)v.E[1];
    236 	result.E[2] = (f32)v.E[2];
    237 	return result;
    238 }
    239 
    240 function v3
    241 v3_from_f32_array(f32 v[3])
    242 {
    243 	v3 result;
    244 	result.E[0] = v[0];
    245 	result.E[1] = v[1];
    246 	result.E[2] = v[2];
    247 	return result;
    248 }
    249 
    250 function v3
    251 v3_abs(v3 a)
    252 {
    253 	v3 result;
    254 	result.x = ABS(a.x);
    255 	result.y = ABS(a.y);
    256 	result.z = ABS(a.z);
    257 	return result;
    258 }
    259 
    260 function v3
    261 v3_scale(v3 a, f32 scale)
    262 {
    263 	v3 result;
    264 	result.x = scale * a.x;
    265 	result.y = scale * a.y;
    266 	result.z = scale * a.z;
    267 	return result;
    268 }
    269 
    270 function v3
    271 v3_add(v3 a, v3 b)
    272 {
    273 	v3 result;
    274 	result.x = a.x + b.x;
    275 	result.y = a.y + b.y;
    276 	result.z = a.z + b.z;
    277 	return result;
    278 }
    279 
    280 function v3
    281 v3_sub(v3 a, v3 b)
    282 {
    283 	v3 result = v3_add(a, v3_scale(b, -1.0f));
    284 	return result;
    285 }
    286 
    287 function v3
    288 v3_div(v3 a, v3 b)
    289 {
    290 	v3 result;
    291 	result.x = a.x / b.x;
    292 	result.y = a.y / b.y;
    293 	result.z = a.z / b.z;
    294 	return result;
    295 }
    296 
    297 function f32
    298 v3_dot(v3 a, v3 b)
    299 {
    300 	f32 result = a.x * b.x + a.y * b.y + a.z * b.z;
    301 	return result;
    302 }
    303 
    304 function f32
    305 v3_magnitude_squared(v3 a)
    306 {
    307 	f32 result = v3_dot(a, a);
    308 	return result;
    309 }
    310 
    311 function f32
    312 v3_magnitude(v3 a)
    313 {
    314 	f32 result = sqrt_f32(v3_dot(a, a));
    315 	return result;
    316 }
    317 
    318 function v3
    319 v3_normalize(v3 a)
    320 {
    321 	v3 result = v3_scale(a, 1.0f / v3_magnitude(a));
    322 	return result;
    323 }
    324 
    325 function v4
    326 v4_scale(v4 a, f32 scale)
    327 {
    328 	v4 result;
    329 	result.x = scale * a.x;
    330 	result.y = scale * a.y;
    331 	result.z = scale * a.z;
    332 	result.w = scale * a.w;
    333 	return result;
    334 }
    335 
    336 function v4
    337 v4_add(v4 a, v4 b)
    338 {
    339 	v4 result;
    340 	result.x = a.x + b.x;
    341 	result.y = a.y + b.y;
    342 	result.z = a.z + b.z;
    343 	result.w = a.w + b.w;
    344 	return result;
    345 }
    346 
    347 function v4
    348 v4_sub(v4 a, v4 b)
    349 {
    350 	v4 result = v4_add(a, v4_scale(b, -1));
    351 	return result;
    352 }
    353 
    354 function f32
    355 v4_dot(v4 a, v4 b)
    356 {
    357 	f32 result = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
    358 	return result;
    359 }
    360 
    361 function v4
    362 v4_lerp(v4 a, v4 b, f32 t)
    363 {
    364 	v4 result = v4_add(a, v4_scale(v4_sub(b, a), t));
    365 	return result;
    366 }
    367 
    368 function m4
    369 m4_identity(void)
    370 {
    371 	m4 result;
    372 	result.c[0] = (v4){{1, 0, 0, 0}};
    373 	result.c[1] = (v4){{0, 1, 0, 0}};
    374 	result.c[2] = (v4){{0, 0, 1, 0}};
    375 	result.c[3] = (v4){{0, 0, 0, 1}};
    376 	return result;
    377 }
    378 
    379 function v4
    380 m4_row(m4 a, u32 row)
    381 {
    382 	v4 result;
    383 	result.E[0] = a.c[0].E[row];
    384 	result.E[1] = a.c[1].E[row];
    385 	result.E[2] = a.c[2].E[row];
    386 	result.E[3] = a.c[3].E[row];
    387 	return result;
    388 }
    389 
    390 function m4
    391 m4_mul(m4 a, m4 b)
    392 {
    393 	m4 result;
    394 	for (u32 i = 0; i < 4; i++) {
    395 		for (u32 j = 0; j < 4; j++) {
    396 			result.c[i].E[j] = v4_dot(m4_row(a, j), b.c[i]);
    397 		}
    398 	}
    399 	return result;
    400 }
    401 
    402 /* NOTE(rnp): based on:
    403  * https://web.archive.org/web/20131215123403/ftp://download.intel.com/design/PentiumIII/sml/24504301.pdf
    404  * TODO(rnp): redo with SIMD as given in the link (but need to rewrite for column-major)
    405  */
    406 function m4
    407 m4_inverse(m4 m)
    408 {
    409 	m4 result;
    410 	result.E[ 0] =  m.E[5] * m.E[10] * m.E[15] - m.E[5] * m.E[11] * m.E[14] - m.E[9] * m.E[6] * m.E[15] + m.E[9] * m.E[7] * m.E[14] + m.E[13] * m.E[6] * m.E[11] - m.E[13] * m.E[7] * m.E[10];
    411 	result.E[ 4] = -m.E[4] * m.E[10] * m.E[15] + m.E[4] * m.E[11] * m.E[14] + m.E[8] * m.E[6] * m.E[15] - m.E[8] * m.E[7] * m.E[14] - m.E[12] * m.E[6] * m.E[11] + m.E[12] * m.E[7] * m.E[10];
    412 	result.E[ 8] =  m.E[4] * m.E[ 9] * m.E[15] - m.E[4] * m.E[11] * m.E[13] - m.E[8] * m.E[5] * m.E[15] + m.E[8] * m.E[7] * m.E[13] + m.E[12] * m.E[5] * m.E[11] - m.E[12] * m.E[7] * m.E[ 9];
    413 	result.E[12] = -m.E[4] * m.E[ 9] * m.E[14] + m.E[4] * m.E[10] * m.E[13] + m.E[8] * m.E[5] * m.E[14] - m.E[8] * m.E[6] * m.E[13] - m.E[12] * m.E[5] * m.E[10] + m.E[12] * m.E[6] * m.E[ 9];
    414 	result.E[ 1] = -m.E[1] * m.E[10] * m.E[15] + m.E[1] * m.E[11] * m.E[14] + m.E[9] * m.E[2] * m.E[15] - m.E[9] * m.E[3] * m.E[14] - m.E[13] * m.E[2] * m.E[11] + m.E[13] * m.E[3] * m.E[10];
    415 	result.E[ 5] =  m.E[0] * m.E[10] * m.E[15] - m.E[0] * m.E[11] * m.E[14] - m.E[8] * m.E[2] * m.E[15] + m.E[8] * m.E[3] * m.E[14] + m.E[12] * m.E[2] * m.E[11] - m.E[12] * m.E[3] * m.E[10];
    416 	result.E[ 9] = -m.E[0] * m.E[ 9] * m.E[15] + m.E[0] * m.E[11] * m.E[13] + m.E[8] * m.E[1] * m.E[15] - m.E[8] * m.E[3] * m.E[13] - m.E[12] * m.E[1] * m.E[11] + m.E[12] * m.E[3] * m.E[ 9];
    417 	result.E[13] =  m.E[0] * m.E[ 9] * m.E[14] - m.E[0] * m.E[10] * m.E[13] - m.E[8] * m.E[1] * m.E[14] + m.E[8] * m.E[2] * m.E[13] + m.E[12] * m.E[1] * m.E[10] - m.E[12] * m.E[2] * m.E[ 9];
    418 	result.E[ 2] =  m.E[1] * m.E[ 6] * m.E[15] - m.E[1] * m.E[ 7] * m.E[14] - m.E[5] * m.E[2] * m.E[15] + m.E[5] * m.E[3] * m.E[14] + m.E[13] * m.E[2] * m.E[ 7] - m.E[13] * m.E[3] * m.E[ 6];
    419 	result.E[ 6] = -m.E[0] * m.E[ 6] * m.E[15] + m.E[0] * m.E[ 7] * m.E[14] + m.E[4] * m.E[2] * m.E[15] - m.E[4] * m.E[3] * m.E[14] - m.E[12] * m.E[2] * m.E[ 7] + m.E[12] * m.E[3] * m.E[ 6];
    420 	result.E[10] =  m.E[0] * m.E[ 5] * m.E[15] - m.E[0] * m.E[ 7] * m.E[13] - m.E[4] * m.E[1] * m.E[15] + m.E[4] * m.E[3] * m.E[13] + m.E[12] * m.E[1] * m.E[ 7] - m.E[12] * m.E[3] * m.E[ 5];
    421 	result.E[14] = -m.E[0] * m.E[ 5] * m.E[14] + m.E[0] * m.E[ 6] * m.E[13] + m.E[4] * m.E[1] * m.E[14] - m.E[4] * m.E[2] * m.E[13] - m.E[12] * m.E[1] * m.E[ 6] + m.E[12] * m.E[2] * m.E[ 5];
    422 	result.E[ 3] = -m.E[1] * m.E[ 6] * m.E[11] + m.E[1] * m.E[ 7] * m.E[10] + m.E[5] * m.E[2] * m.E[11] - m.E[5] * m.E[3] * m.E[10] - m.E[ 9] * m.E[2] * m.E[ 7] + m.E[ 9] * m.E[3] * m.E[ 6];
    423 	result.E[ 7] =  m.E[0] * m.E[ 6] * m.E[11] - m.E[0] * m.E[ 7] * m.E[10] - m.E[4] * m.E[2] * m.E[11] + m.E[4] * m.E[3] * m.E[10] + m.E[ 8] * m.E[2] * m.E[ 7] - m.E[ 8] * m.E[3] * m.E[ 6];
    424 	result.E[11] = -m.E[0] * m.E[ 5] * m.E[11] + m.E[0] * m.E[ 7] * m.E[ 9] + m.E[4] * m.E[1] * m.E[11] - m.E[4] * m.E[3] * m.E[ 9] - m.E[ 8] * m.E[1] * m.E[ 7] + m.E[ 8] * m.E[3] * m.E[ 5];
    425 	result.E[15] =  m.E[0] * m.E[ 5] * m.E[10] - m.E[0] * m.E[ 6] * m.E[ 9] - m.E[4] * m.E[1] * m.E[10] + m.E[4] * m.E[2] * m.E[ 9] + m.E[ 8] * m.E[1] * m.E[ 6] - m.E[ 8] * m.E[2] * m.E[ 5];
    426 
    427 	f32 determinant = m.E[0] * result.E[0] + m.E[1] * result.E[4] + m.E[2] * result.E[8] + m.E[3] * result.E[12];
    428 	determinant = 1.0f / determinant;
    429 	for(i32 i = 0; i < 16; i++)
    430 		result.E[i] *= determinant;
    431 	return result;
    432 }
    433 
    434 function m4
    435 m4_translation(v3 delta)
    436 {
    437 	m4 result;
    438 	result.c[0] = (v4){{1, 0, 0, 0}};
    439 	result.c[1] = (v4){{0, 1, 0, 0}};
    440 	result.c[2] = (v4){{0, 0, 1, 0}};
    441 	result.c[3] = (v4){{delta.x, delta.y, delta.z, 1}};
    442 	return result;
    443 }
    444 
    445 function m4
    446 m4_scale(v3 scale)
    447 {
    448 	m4 result;
    449 	result.c[0] = (v4){{scale.x, 0,       0,       0}};
    450 	result.c[1] = (v4){{0,       scale.y, 0,       0}};
    451 	result.c[2] = (v4){{0,       0,       scale.z, 0}};
    452 	result.c[3] = (v4){{0,       0,       0,       1}};
    453 	return result;
    454 }
    455 
    456 function m4
    457 m4_rotation_about_z(f32 turns)
    458 {
    459 	f32 sa = sin_f32(turns * 2 * PI);
    460 	f32 ca = cos_f32(turns * 2 * PI);
    461 	m4 result;
    462 	result.c[0] = (v4){{ca, -sa, 0, 0}};
    463 	result.c[1] = (v4){{sa,  ca, 0, 0}};
    464 	result.c[2] = (v4){{0,    0, 1, 0}};
    465 	result.c[3] = (v4){{0,    0, 0, 1}};
    466 	return result;
    467 }
    468 
    469 function m4
    470 m4_rotation_about_y(f32 turns)
    471 {
    472 	f32 sa = sin_f32(turns * 2 * PI);
    473 	f32 ca = cos_f32(turns * 2 * PI);
    474 	m4 result;
    475 	result.c[0] = (v4){{ca, 0, -sa, 0}};
    476 	result.c[1] = (v4){{0,  1,  0,  0}};
    477 	result.c[2] = (v4){{sa, 0,  ca, 0}};
    478 	result.c[3] = (v4){{0,  0,  0,  1}};
    479 	return result;
    480 }
    481 
    482 function m4
    483 y_aligned_volume_transform(v3 extent, v3 translation, f32 rotation_turns)
    484 {
    485 	m4 T = m4_translation(translation);
    486 	m4 R = m4_rotation_about_y(rotation_turns);
    487 	m4 S = m4_scale(extent);
    488 	m4 result = m4_mul(T, m4_mul(R, S));
    489 	return result;
    490 }
    491 
    492 function v4
    493 m4_mul_v4(m4 a, v4 v)
    494 {
    495 	v4 result;
    496 	result.x = v4_dot(m4_row(a, 0), v);
    497 	result.y = v4_dot(m4_row(a, 1), v);
    498 	result.z = v4_dot(m4_row(a, 2), v);
    499 	result.w = v4_dot(m4_row(a, 3), v);
    500 	return result;
    501 }
    502 
    503 function m4
    504 orthographic_projection(f32 n, f32 f, f32 t, f32 r)
    505 {
    506 	m4 result;
    507 	f32 a = -2 / (f - n);
    508 	f32 b = - (f + n) / (f - n);
    509 	result.c[0] = (v4){{1 / r, 0,     0,  0}};
    510 	result.c[1] = (v4){{0,     1 / t, 0,  0}};
    511 	result.c[2] = (v4){{0,     0,     a,  0}};
    512 	result.c[3] = (v4){{0,     0,     b,  1}};
    513 	return result;
    514 }
    515 
    516 function m4
    517 perspective_projection(f32 n, f32 f, f32 fov, f32 aspect)
    518 {
    519 	m4 result;
    520 	f32 t = tan_f32(fov / 2.0f);
    521 	f32 r = t * aspect;
    522 	f32 a = -(f + n) / (f - n);
    523 	f32 b = -2 * f * n / (f - n);
    524 	result.c[0] = (v4){{1 / r, 0,     0,  0}};
    525 	result.c[1] = (v4){{0,     1 / t, 0,  0}};
    526 	result.c[2] = (v4){{0,     0,     a, -1}};
    527 	result.c[3] = (v4){{0,     0,     b,  0}};
    528 	return result;
    529 }
    530 
    531 function m4
    532 camera_look_at(v3 camera, v3 point)
    533 {
    534 	v3 orthogonal = {{0, 1.0f, 0}};
    535 	v3 normal     = v3_normalize(v3_sub(camera, point));
    536 	v3 right      = cross(orthogonal, normal);
    537 	v3 up         = cross(normal,     right);
    538 
    539 	v3 translate;
    540 	camera      = v3_sub((v3){0}, camera);
    541 	translate.x = v3_dot(camera, right);
    542 	translate.y = v3_dot(camera, up);
    543 	translate.z = v3_dot(camera, normal);
    544 
    545 	m4 result;
    546 	result.c[0] = (v4){{right.x,     up.x,        normal.x,    0}};
    547 	result.c[1] = (v4){{right.y,     up.y,        normal.y,    0}};
    548 	result.c[2] = (v4){{right.z,     up.z,        normal.z,    0}};
    549 	result.c[3] = (v4){{translate.x, translate.y, translate.z, 1}};
    550 	return result;
    551 }
    552 
    553 /* NOTE(rnp): adapted from "Essential Mathematics for Games and Interactive Applications" (Verth, Bishop) */
    554 function f32
    555 obb_raycast(m4 obb_orientation, v3 obb_size, v3 obb_center, ray r)
    556 {
    557 	v3 p = v3_sub(obb_center, r.origin);
    558 	v3 X = obb_orientation.c[0].xyz;
    559 	v3 Y = obb_orientation.c[1].xyz;
    560 	v3 Z = obb_orientation.c[2].xyz;
    561 
    562 	/* NOTE(rnp): projects direction vector onto OBB axis */
    563 	v3 f;
    564 	f.x = v3_dot(X, r.direction);
    565 	f.y = v3_dot(Y, r.direction);
    566 	f.z = v3_dot(Z, r.direction);
    567 
    568 	/* NOTE(rnp): projects relative vector onto OBB axis */
    569 	v3 e;
    570 	e.x = v3_dot(X, p);
    571 	e.y = v3_dot(Y, p);
    572 	e.z = v3_dot(Z, p);
    573 
    574 	f32 result = 0;
    575 	f32 t[6] = {0};
    576 	for (i32 i = 0; i < 3; i++) {
    577 		if (f32_cmp(f.E[i], 0)) {
    578 			if (-e.E[i] - obb_size.E[i] > 0 || -e.E[i] + obb_size.E[i] < 0)
    579 				result = -1.0f;
    580 			f.E[i] = F32_EPSILON;
    581 		}
    582 		t[i * 2 + 0] = (e.E[i] + obb_size.E[i]) / f.E[i];
    583 		t[i * 2 + 1] = (e.E[i] - obb_size.E[i]) / f.E[i];
    584 	}
    585 
    586 	if (result != -1) {
    587 		f32 tmin = MAX(MAX(MIN(t[0], t[1]), MIN(t[2], t[3])), MIN(t[4], t[5]));
    588 		f32 tmax = MIN(MIN(MAX(t[0], t[1]), MAX(t[2], t[3])), MAX(t[4], t[5]));
    589 		if (tmax >= 0 && tmin <= tmax) {
    590 			result = tmin > 0 ? tmin : tmax;
    591 		} else {
    592 			result = -1;
    593 		}
    594 	}
    595 
    596 	return result;
    597 }
    598 
    599 function f32
    600 complex_filter_first_moment(v2 *filter, i32 length, f32 sampling_frequency)
    601 {
    602 	f32 n = 0, d = 0;
    603 	for (i32 i = 0; i < length; i++) {
    604 		f32 t = v2_magnitude_squared(filter[i]);
    605 		n += (f32)i * t;
    606 		d += t;
    607 	}
    608 	f32 result = n / d / sampling_frequency;
    609 	return result;
    610 }
    611 
    612 function f32
    613 real_filter_first_moment(f32 *filter, i32 length, f32 sampling_frequency)
    614 {
    615 	f32 n = 0, d = 0;
    616 	for (i32 i = 0; i < length; i++) {
    617 		f32 t = filter[i] * filter[i];
    618 		n += (f32)i * t;
    619 		d += t;
    620 	}
    621 	f32 result = n / d / sampling_frequency;
    622 	return result;
    623 }
    624 
    625 function f32
    626 tukey_window(f32 t, f32 tapering)
    627 {
    628 	f32 r = tapering;
    629 	f32 result = 1;
    630 	if (t < r / 2)      result = 0.5f * (1 + cos_f32(2 * PI * (t - r / 2)     / r));
    631 	if (t >= 1 - r / 2) result = 0.5f * (1 + cos_f32(2 * PI * (t - 1 + r / 2) / r));
    632 	return result;
    633 }
    634 
    635 /* NOTE(rnp): adapted from "Discrete Time Signal Processing" (Oppenheim) */
    636 function f32 *
    637 kaiser_low_pass_filter(Arena *arena, f32 cutoff_frequency, f32 sampling_frequency, f32 beta, i32 length)
    638 {
    639 	f32 *result = push_array(arena, f32, length);
    640 	f32 wc      = 2 * PI * cutoff_frequency / sampling_frequency;
    641 	f32 a       = (f32)length / 2.0f;
    642 	f32 pi_i0_b = PI * (f32)cephes_i0(beta);
    643 
    644 	for (i32 n = 0; n < length; n++) {
    645 		f32 t       = (f32)n - a;
    646 		f32 impulse = !f32_cmp(t, 0) ? sin_f32(wc * t) / t : wc;
    647 		t           = t / a;
    648 		f32 window  = (f32)cephes_i0(beta * sqrt_f32(1 - t * t)) / pi_i0_b;
    649 		result[n]   = impulse * window;
    650 	}
    651 
    652 	return result;
    653 }
    654 
    655 function f32 *
    656 rf_chirp(Arena *arena, f32 min_frequency, f32 max_frequency, f32 sampling_frequency,
    657          i32 length, b32 reverse)
    658 {
    659 	f32 *result = push_array(arena, f32, length);
    660 	for (i32 i = 0; i < length; i++) {
    661 		i32 index = reverse? length - 1 - i : i;
    662 		f32 fc    = min_frequency + (f32)i * (max_frequency - min_frequency) / (2 * (f32)length);
    663 		f32 arg   = 2 * PI * fc * (f32)i / sampling_frequency;
    664 		result[index] = sin_f32(arg) * tukey_window((f32)i / (f32)length, 0.2f);
    665 	}
    666 	return result;
    667 }
    668 
    669 function v2 *
    670 baseband_chirp(Arena *arena, f32 min_frequency, f32 max_frequency, f32 sampling_frequency,
    671                i32 length, b32 reverse, f32 scale)
    672 {
    673 	v2 *result    = push_array(arena, v2, length);
    674 	f32 conjugate = reverse ? -1 : 1;
    675 	for (i32 i = 0; i < length; i++) {
    676 		i32 index = reverse? length - 1 - i : i;
    677 		f32 fc    = min_frequency + (f32)i * (max_frequency - min_frequency) / (2 * (f32)length);
    678 		f32 arg   = 2 * PI * fc * (f32)i / sampling_frequency;
    679 		v2 sample = {{scale * cos_f32(arg), conjugate * scale * sin_f32(arg)}};
    680 		result[index] = v2_scale(sample, tukey_window((f32)i / (f32)length, 0.2f));
    681 	}
    682 	return result;
    683 }
    684 
    685 function v4
    686 hsv_to_rgb(v4 hsv)
    687 {
    688 	/* f(k(n))   = V - V*S*max(0, min(k, min(4 - k, 1)))
    689 	 * k(n)      = fmod((n + H * 6), 6)
    690 	 * (R, G, B) = (f(n = 5), f(n = 3), f(n = 1))
    691 	 */
    692 	alignas(16) f32 nval[4] = {5.0f, 3.0f, 1.0f, 0.0f};
    693 	f32x4 n   = load_f32x4(nval);
    694 	f32x4 H   = dup_f32x4(hsv.x);
    695 	f32x4 S   = dup_f32x4(hsv.y);
    696 	f32x4 V   = dup_f32x4(hsv.z);
    697 	f32x4 six = dup_f32x4(6);
    698 
    699 	f32x4 t   = add_f32x4(n, mul_f32x4(six, H));
    700 	f32x4 rem = floor_f32x4(div_f32x4(t, six));
    701 	f32x4 k   = sub_f32x4(t, mul_f32x4(rem, six));
    702 
    703 	t = min_f32x4(sub_f32x4(dup_f32x4(4), k), dup_f32x4(1));
    704 	t = max_f32x4(dup_f32x4(0), min_f32x4(k, t));
    705 	t = mul_f32x4(t, mul_f32x4(S, V));
    706 
    707 	v4 rgba;
    708 	store_f32x4(rgba.E, sub_f32x4(V, t));
    709 	rgba.a = hsv.a;
    710 	return rgba;
    711 }
    712 
    713 function f32
    714 ease_in_out_cubic(f32 t)
    715 {
    716 	f32 result;
    717 	if (t < 0.5f) {
    718 		result = 4.0f * t * t * t;
    719 	} else {
    720 		t      = -2.0f * t + 2.0f;
    721 		result =  1.0f - t * t * t / 2.0f;
    722 	}
    723 	return result;
    724 }
    725 
    726 function f32
    727 ease_in_out_quartic(f32 t)
    728 {
    729 	f32 result;
    730 	if (t < 0.5f) {
    731 		result = 8.0f * t * t * t * t;
    732 	} else {
    733 		t      = -2.0f * t + 2.0f;
    734 		result =  1.0f - t * t * t * t / 2.0f;
    735 	}
    736 	return result;
    737 }