throughput.c (18912B)
1 /* See LICENSE for license details. */ 2 /* TODO(rnp): 3 * [ ]: for finer grained evaluation of throughput latency just queue a data upload 4 * without replacing the data. 5 * [ ]: bug: we aren't inserting rf data between each frame 6 */ 7 8 #define BEAMFORMER_LIB_EXPORT function 9 #include "ogl_beamformer_lib.c" 10 11 #include <signal.h> 12 #include <stdarg.h> 13 #include <stdio.h> 14 #include <stdlib.h> 15 #include <zstd.h> 16 17 global iv3 g_output_points = {{512, 1, 1024}}; 18 global v2 g_axial_extent = {{ 10e-3f, 165e-3f}}; 19 global v2 g_lateral_extent = {{-60e-3f, 60e-3f}}; 20 global f32 g_f_number = 0.5f; 21 22 typedef struct { 23 b32 loop; 24 u32 frame_number; 25 26 char **remaining; 27 i32 remaining_count; 28 } Options; 29 30 #include "external/zemp_bp.h" 31 32 typedef struct { 33 ZBP_DataKind kind; 34 ZBP_DataCompressionKind compression_kind; 35 s8 bytes; 36 } ZBP_Data; 37 38 global b32 g_should_exit; 39 40 #define die(...) die_((char *)__func__, __VA_ARGS__) 41 function no_return void 42 die_(char *function_name, char *format, ...) 43 { 44 if (function_name) 45 fprintf(stderr, "%s: ", function_name); 46 47 va_list ap; 48 49 va_start(ap, format); 50 vfprintf(stderr, format, ap); 51 va_end(ap); 52 53 os_exit(1); 54 } 55 56 #if OS_LINUX 57 58 #include <fcntl.h> 59 #include <sys/stat.h> 60 #include <unistd.h> 61 62 function s8 63 os_read_file_simp(char *fname) 64 { 65 s8 result; 66 i32 fd = open(fname, O_RDONLY); 67 if (fd < 0) 68 die("couldn't open file: %s\n", fname); 69 70 struct stat st; 71 if (stat(fname, &st) < 0) 72 die("couldn't stat file\n"); 73 74 result.len = st.st_size; 75 result.data = malloc((uz)st.st_size); 76 if (!result.data) 77 die("couldn't alloc space for reading\n"); 78 79 iz rlen = read(fd, result.data, (u32)st.st_size); 80 close(fd); 81 82 if (rlen != st.st_size) 83 die("couldn't read file: %s\n", fname); 84 85 return result; 86 } 87 88 #elif OS_WINDOWS 89 90 function s8 91 os_read_file_simp(char *fname) 92 { 93 s8 result; 94 iptr h = CreateFileA(fname, GENERIC_READ, 0, 0, OPEN_EXISTING, 0, 0); 95 if (h == INVALID_FILE) 96 die("couldn't open file: %s\n", fname); 97 98 w32_file_info fileinfo; 99 if (!GetFileInformationByHandle(h, &fileinfo)) 100 die("couldn't get file info\n", stderr); 101 102 result.len = fileinfo.nFileSizeLow; 103 result.data = malloc(fileinfo.nFileSizeLow); 104 if (!result.data) 105 die("couldn't alloc space for reading\n"); 106 107 i32 rlen = 0; 108 if (!ReadFile(h, result.data, (i32)fileinfo.nFileSizeLow, &rlen, 0) && rlen != (i32)fileinfo.nFileSizeLow) 109 die("couldn't read file: %s\n", fname); 110 CloseHandle(h); 111 112 return result; 113 } 114 115 #else 116 #error Unsupported Platform 117 #endif 118 119 function void 120 stream_ensure_termination(Stream *s, u8 byte) 121 { 122 b32 found = 0; 123 if (!s->errors && s->widx > 0) 124 found = s->data[s->widx - 1] == byte; 125 if (!found) { 126 s->errors |= s->cap - 1 < s->widx; 127 if (!s->errors) 128 s->data[s->widx++] = byte; 129 } 130 } 131 132 function void * 133 decompress_zstd_data(s8 raw) 134 { 135 uz requested_size = ZSTD_getFrameContentSize(raw.data, (uz)raw.len); 136 void *out = malloc(requested_size); 137 if (out) { 138 uz decompressed = ZSTD_decompress(out, requested_size, raw.data, (uz)raw.len); 139 if (decompressed != requested_size) { 140 free(out); 141 out = 0; 142 } 143 } 144 return out; 145 } 146 147 function b32 148 beamformer_simple_parameters_from_zbp_file(BeamformerSimpleParameters *bp, char *path, ZBP_Data *raw_data) 149 { 150 s8 raw = os_read_file_simp(path); 151 if (raw.len < (iz)sizeof(ZBP_BaseHeader) || ((ZBP_BaseHeader *)raw.data)->magic != ZBP_HeaderMagic) 152 return 0; 153 154 switch (((ZBP_BaseHeader *)raw.data)->major) { 155 156 case 1:{ 157 ZBP_HeaderV1 *header = (ZBP_HeaderV1 *)raw.data; 158 159 bp->sample_count = header->sample_count; 160 bp->channel_count = header->channel_count; 161 bp->acquisition_count = header->receive_event_count; 162 163 bp->sampling_mode = BeamformerSamplingMode_4X; 164 bp->acquisition_kind = header->beamform_mode; 165 bp->decode_mode = header->decode_mode; 166 bp->sampling_frequency = header->sampling_frequency; 167 bp->demodulation_frequency = header->sampling_frequency / 4; 168 bp->speed_of_sound = header->speed_of_sound; 169 bp->time_offset = header->time_offset; 170 171 mem_copy(bp->channel_mapping, header->channel_mapping, sizeof(*bp->channel_mapping) * bp->channel_count); 172 mem_copy(bp->xdc_transform.E, header->transducer_transform_matrix, sizeof(bp->xdc_transform)); 173 mem_copy(bp->xdc_element_pitch.E, header->transducer_element_pitch, sizeof(bp->xdc_element_pitch)); 174 // NOTE(rnp): ignores emission count and ensemble count 175 mem_copy(bp->raw_data_dimensions.E, header->raw_data_dimension, sizeof(bp->raw_data_dimensions)); 176 177 bp->data_kind = (BeamformerDataKind)ZBP_DataKind_Int16; 178 raw_data->kind = ZBP_DataKind_Int16; 179 raw_data->compression_kind = ZBP_DataCompressionKind_ZSTD; 180 181 read_only local_persist u8 transmit_mode_to_orientation[] = { 182 [0] = (ZBP_RCAOrientation_Rows << 4) | ZBP_RCAOrientation_Rows, 183 [1] = (ZBP_RCAOrientation_Rows << 4) | ZBP_RCAOrientation_Columns, 184 [2] = (ZBP_RCAOrientation_Columns << 4) | ZBP_RCAOrientation_Rows, 185 [3] = (ZBP_RCAOrientation_Columns << 4) | ZBP_RCAOrientation_Columns, 186 }; 187 if (header->transmit_mode >= countof(transmit_mode_to_orientation)) 188 return 0; 189 190 bp->transmit_receive_orientation = transmit_mode_to_orientation[header->transmit_mode]; 191 192 ZBP_AcquisitionKind acquisition_kind = header->beamform_mode; 193 if (acquisition_kind == ZBP_AcquisitionKind_FORCES || 194 acquisition_kind == ZBP_AcquisitionKind_HERCULES || 195 acquisition_kind == ZBP_AcquisitionKind_UFORCES || 196 acquisition_kind == ZBP_AcquisitionKind_UHERCULES) 197 { 198 bp->single_focus = 1; 199 bp->single_orientation = 1; 200 bp->focal_vector.E[0] = header->steering_angles[0]; 201 bp->focal_vector.E[1] = header->focal_depths[0]; 202 } 203 204 if (acquisition_kind == ZBP_AcquisitionKind_UFORCES || 205 acquisition_kind == ZBP_AcquisitionKind_UHERCULES) 206 { 207 mem_copy(bp->sparse_elements, header->sparse_elements, sizeof(*bp->sparse_elements) * bp->acquisition_count); 208 } 209 210 if (acquisition_kind == ZBP_AcquisitionKind_RCA_TPW || 211 acquisition_kind == ZBP_AcquisitionKind_RCA_VLS) 212 { 213 mem_copy(bp->focal_depths, header->focal_depths, sizeof(*bp->focal_depths) * bp->acquisition_count); 214 mem_copy(bp->steering_angles, header->steering_angles, sizeof(*bp->steering_angles) * bp->acquisition_count); 215 for EachIndex(bp->acquisition_count, it) 216 bp->transmit_receive_orientations[it] = bp->transmit_receive_orientation; 217 } 218 219 bp->emission_parameters.kind = BeamformerEmissionKind_Sine; 220 bp->emission_parameters.sine.cycles = 2; 221 bp->emission_parameters.sine.frequency = bp->demodulation_frequency; 222 }break; 223 224 case 2:{ 225 ZBP_HeaderV2 *header = (ZBP_HeaderV2 *)raw.data; 226 227 bp->sample_count = header->sample_count; 228 bp->channel_count = header->channel_count; 229 bp->acquisition_count = header->receive_event_count; 230 231 read_only local_persist BeamformerSamplingMode zbp_sampling_mode_to_beamformer[] = { 232 [ZBP_SamplingMode_Standard] = BeamformerSamplingMode_4X, 233 [ZBP_SamplingMode_Bandpass] = BeamformerSamplingMode_2X, 234 }; 235 bp->sampling_mode = zbp_sampling_mode_to_beamformer[header->sampling_mode]; 236 237 bp->acquisition_kind = header->acquisition_mode; 238 bp->decode_mode = header->decode_mode; 239 bp->sampling_frequency = header->sampling_frequency; 240 bp->demodulation_frequency = header->demodulation_frequency; 241 bp->speed_of_sound = header->speed_of_sound; 242 bp->time_offset = header->time_offset; 243 244 bp->contrast_mode = header->contrast_mode; 245 246 if (header->channel_mapping_offset != -1) { 247 mem_copy(bp->channel_mapping, raw.data + header->channel_mapping_offset, 248 sizeof(*bp->channel_mapping) * bp->channel_count); 249 } else { 250 for EachIndex(bp->channel_count, it) 251 bp->channel_mapping[it] = it; 252 } 253 254 mem_copy(bp->xdc_transform.E, header->transducer_transform_matrix, sizeof(bp->xdc_transform)); 255 mem_copy(bp->xdc_element_pitch.E, header->transducer_element_pitch, sizeof(bp->xdc_element_pitch)); 256 // NOTE(rnp): ignores group count and ensemble count 257 mem_copy(bp->raw_data_dimensions.E, header->raw_data_dimension, sizeof(bp->raw_data_dimensions)); 258 259 bp->data_kind = header->raw_data_kind; 260 raw_data->kind = header->raw_data_kind; 261 raw_data->compression_kind = header->raw_data_compression_kind; 262 263 if (header->raw_data_offset != -1) { 264 raw_data->bytes.data = raw.data + header->raw_data_offset; 265 if (raw_data->compression_kind == ZBP_DataCompressionKind_ZSTD) { 266 // NOTE(rnp): limitation in the header format 267 raw_data->bytes.len = raw.len - header->raw_data_offset; 268 } else { 269 raw_data->bytes.len = header->raw_data_dimension[0] * header->raw_data_dimension[1] * 270 header->raw_data_dimension[2] * header->raw_data_dimension[3]; 271 raw_data->bytes.len *= beamformer_data_kind_byte_size[header->raw_data_kind]; 272 } 273 } 274 275 // NOTE(rnp): only look at the first emission descriptor, other cases aren't currently relevant 276 { 277 ZBP_EmissionDescriptor *ed = (ZBP_EmissionDescriptor *)(raw.data + header->emission_descriptors_offset); 278 switch (ed->emission_kind) { 279 280 case ZBP_EmissionKind_Sine:{ 281 ZBP_EmissionSineParameters *ep = (ZBP_EmissionSineParameters *)(raw.data + ed->parameters_offset); 282 bp->emission_parameters.kind = BeamformerEmissionKind_Sine; 283 bp->emission_parameters.sine.cycles = ep->cycles; 284 bp->emission_parameters.sine.frequency = ep->frequency; 285 }break; 286 287 case ZBP_EmissionKind_Chirp:{ 288 ZBP_EmissionChirpParameters *ep = (ZBP_EmissionChirpParameters *)(raw.data + ed->parameters_offset); 289 bp->emission_parameters.kind = BeamformerEmissionKind_Chirp; 290 bp->emission_parameters.chirp.duration = ep->duration; 291 bp->emission_parameters.chirp.min_frequency = ep->min_frequency; 292 bp->emission_parameters.chirp.max_frequency = ep->max_frequency; 293 }break; 294 295 InvalidDefaultCase; 296 static_assert(ZBP_EmissionKind_Count == (ZBP_EmissionKind_Chirp + 1), ""); 297 } 298 } 299 300 switch (header->acquisition_mode) { 301 case ZBP_AcquisitionKind_FORCES:{}break; 302 303 case ZBP_AcquisitionKind_HERCULES:{ 304 ZBP_HERCULESParameters *p = (ZBP_HERCULESParameters *)(raw.data + header->acquisition_parameters_offset); 305 bp->transmit_receive_orientation = p->transmit_focus.transmit_receive_orientation; 306 bp->focal_vector.E[0] = p->transmit_focus.steering_angle; 307 bp->focal_vector.E[1] = p->transmit_focus.focal_depth; 308 309 bp->single_focus = 1; 310 bp->single_orientation = 1; 311 }break; 312 313 case ZBP_AcquisitionKind_UFORCES:{ 314 ZBP_uFORCESParameters *p = (ZBP_uFORCESParameters *)(raw.data + header->acquisition_parameters_offset); 315 mem_copy(bp->sparse_elements, raw.data + p->sparse_elements_offset, 316 sizeof(*bp->sparse_elements) * bp->acquisition_count); 317 }break; 318 319 case ZBP_AcquisitionKind_UHERCULES:{ 320 ZBP_uHERCULESParameters *p = (ZBP_uHERCULESParameters *)(raw.data + header->acquisition_parameters_offset); 321 bp->transmit_receive_orientation = p->transmit_focus.transmit_receive_orientation; 322 bp->focal_vector.E[0] = p->transmit_focus.steering_angle; 323 bp->focal_vector.E[1] = p->transmit_focus.focal_depth; 324 325 bp->single_focus = 1; 326 bp->single_orientation = 1; 327 328 mem_copy(bp->sparse_elements, raw.data + p->sparse_elements_offset, 329 sizeof(*bp->sparse_elements) * bp->acquisition_count); 330 }break; 331 332 case ZBP_AcquisitionKind_RCA_TPW:{ 333 ZBP_TPWParameters *p = (ZBP_TPWParameters *)(raw.data + header->acquisition_parameters_offset); 334 335 mem_copy(bp->transmit_receive_orientations, raw.data + p->transmit_receive_orientations_offset, 336 sizeof(*bp->transmit_receive_orientations) * bp->acquisition_count); 337 mem_copy(bp->steering_angles, raw.data + p->tilting_angles_offset, 338 sizeof(*bp->steering_angles) * bp->acquisition_count); 339 340 for EachIndex(bp->acquisition_count, it) 341 bp->focal_depths[it] = F32_INFINITY; 342 }break; 343 344 case ZBP_AcquisitionKind_RCA_VLS:{ 345 ZBP_VLSParameters *p = (ZBP_VLSParameters *)(raw.data + header->acquisition_parameters_offset); 346 347 mem_copy(bp->transmit_receive_orientations, raw.data + p->transmit_receive_orientations_offset, 348 sizeof(*bp->transmit_receive_orientations) * bp->acquisition_count); 349 350 f32 *focal_depths = (f32 *)(raw.data + p->focal_depths_offset); 351 f32 *origin_offsets = (f32 *)(raw.data + p->origin_offsets_offset); 352 353 for EachIndex(bp->acquisition_count, it) { 354 f32 sign = Sign(focal_depths[it]); 355 f32 depth = focal_depths[it]; 356 f32 origin = origin_offsets[it]; 357 bp->steering_angles[it] = atan2_f32(origin, -depth) * 180.0f / PI; 358 bp->focal_depths[it] = sign * sqrt_f32(depth * depth + origin * origin); 359 } 360 }break; 361 362 InvalidDefaultCase; 363 } 364 365 }break; 366 367 default:{return 0;}break; 368 } 369 370 return 1; 371 } 372 373 #define shift_n(v, c, n) v += n, c -= n 374 #define shift(v, c) shift_n(v, c, 1) 375 376 function void 377 usage(char *argv0) 378 { 379 die("%s [--loop] [--frame n] parameters_file\n" 380 " --loop: reupload data forever\n" 381 " --frame n: use frame n of the data for display\n", 382 argv0); 383 } 384 385 function Options 386 parse_argv(i32 argc, char *argv[]) 387 { 388 Options result = {0}; 389 390 char *argv0 = argv[0]; 391 shift(argv, argc); 392 393 while (argc > 0) { 394 s8 arg = c_str_to_s8(*argv); 395 396 if (s8_equal(arg, s8("--loop"))) { 397 shift(argv, argc); 398 result.loop = 1; 399 } else if (s8_equal(arg, s8("--frame"))) { 400 shift(argv, argc); 401 if (argc) { 402 result.frame_number = (u32)atoi(*argv); 403 shift(argv, argc); 404 } 405 } else if (arg.len > 0 && arg.data[0] == '-') { 406 usage(argv0); 407 } else { 408 break; 409 } 410 } 411 412 result.remaining = argv; 413 result.remaining_count = argc; 414 415 return result; 416 } 417 418 function b32 419 send_frame(void *restrict data, BeamformerSimpleParameters *restrict bp) 420 { 421 u32 data_size = bp->raw_data_dimensions.E[0] * bp->raw_data_dimensions.E[1] 422 * beamformer_data_kind_byte_size[bp->data_kind]; 423 b32 result = beamformer_push_data_with_compute(data, data_size, BeamformerViewPlaneTag_XZ, 0); 424 if (!result && !g_should_exit) printf("lib error: %s\n", beamformer_get_last_error_string()); 425 426 return result; 427 } 428 429 function void 430 execute_study(Arena arena, Stream path, Options *options) 431 { 432 i32 path_work_index = path.widx; 433 stream_ensure_termination(&path, 0); 434 435 ZBP_Data raw_data = {0}; 436 BeamformerSimpleParameters bp = {0}; 437 if (!beamformer_simple_parameters_from_zbp_file(&bp, (char *)path.data, &raw_data)) 438 die("failed to load parameters file: %s\n", (char *)path.data); 439 440 v3 min_coordinate = (v3){{g_lateral_extent.x, g_axial_extent.x, 0}}; 441 v3 max_coordinate = (v3){{g_lateral_extent.y, g_axial_extent.y, 0}}; 442 bp.das_voxel_transform = das_transform(min_coordinate, max_coordinate, &g_output_points); 443 444 bp.output_points.xyz = g_output_points; 445 bp.output_points.w = 1; 446 447 bp.f_number = g_f_number; 448 bp.interpolation_mode = BeamformerInterpolationMode_Cubic; 449 450 bp.decimation_rate = 1; 451 452 if (bp.data_kind != BeamformerDataKind_Float32Complex && 453 bp.data_kind != BeamformerDataKind_Int16Complex) 454 { 455 bp.compute_stages[bp.compute_stages_count++] = BeamformerShaderKind_Demodulate; 456 } 457 bp.compute_stages[bp.compute_stages_count++] = BeamformerShaderKind_Decode; 458 bp.compute_stages[bp.compute_stages_count++] = BeamformerShaderKind_DAS; 459 460 { 461 BeamformerFilterParameters filter = {.sampling_frequency = bp.sampling_frequency / 2}; 462 463 BeamformerEmissionParameters *ep = &bp.emission_parameters; 464 switch (bp.emission_parameters.kind) { 465 466 case BeamformerEmissionKind_Sine:{ 467 filter.kind = BeamformerFilterKind_Kaiser; 468 filter.kaiser.beta = 5.65f; 469 filter.kaiser.cutoff_frequency = 0.5f * ep->sine.frequency; 470 filter.kaiser.length = 36; 471 }break; 472 473 case BeamformerEmissionKind_Chirp:{ 474 filter.kind = BeamformerFilterKind_MatchedChirp; 475 filter.matched_chirp.duration = ep->chirp.duration; 476 filter.matched_chirp.min_frequency = ep->chirp.min_frequency - bp.demodulation_frequency; 477 filter.matched_chirp.max_frequency = ep->chirp.max_frequency - bp.demodulation_frequency; 478 filter.complex = 1; 479 480 //bp.time_offset += ep->chirp.duration / 2; 481 }break; 482 483 InvalidDefaultCase; 484 } 485 486 beamformer_create_filter(&filter, 0, 0); 487 488 bp.compute_stage_parameters[0] = 0; 489 } 490 491 beamformer_push_simple_parameters(&bp); 492 493 beamformer_set_global_timeout(1000); 494 495 void *data = 0; 496 if (raw_data.bytes.len == 0) { 497 // NOTE(rnp): strip ".bp" 498 stream_reset(&path, path_work_index - 3); 499 500 stream_append_byte(&path, '_'); 501 stream_append_u64_width(&path, options->frame_number, 2); 502 stream_append_s8(&path, s8(".zst")); 503 stream_ensure_termination(&path, 0); 504 s8 compressed_data = os_read_file_simp((char *)path.data); 505 506 data = decompress_zstd_data(compressed_data); 507 if (!data) 508 die("failed to decompress data: %s\n", path.data); 509 free(compressed_data.data); 510 } else { 511 if (raw_data.compression_kind == ZBP_DataCompressionKind_ZSTD) { 512 data = decompress_zstd_data(raw_data.bytes); 513 if (!data) 514 die("failed to decompress data: %s\n", path.data); 515 } else { 516 data = raw_data.bytes.data; 517 } 518 } 519 520 if (options->loop) { 521 BeamformerLiveImagingParameters lip = {.active = 1, .save_enabled = 1}; 522 s8 short_name = s8("Throughput"); 523 mem_copy(lip.save_name_tag, short_name.data, (uz)short_name.len); 524 lip.save_name_tag_length = (i32)short_name.len; 525 beamformer_set_live_parameters(&lip); 526 527 u32 frame = 0; 528 f32 times[32] = {0}; 529 f32 data_size = (f32)(bp.raw_data_dimensions.E[0] * bp.raw_data_dimensions.E[1] 530 * beamformer_data_kind_byte_size[bp.data_kind]); 531 u64 start = os_timer_count(); 532 f64 frequency = os_timer_frequency(); 533 for (;!g_should_exit;) { 534 if (send_frame(data, &bp)) { 535 u64 now = os_timer_count(); 536 f64 delta = (now - start) / frequency; 537 start = now; 538 539 if ((frame % 16) == 0) { 540 f32 sum = 0; 541 for (u32 i = 0; i < countof(times); i++) 542 sum += times[i] / countof(times); 543 printf("Frame Time: %8.3f [ms] | 32-Frame Average: %8.3f [ms] | %8.3f GB/s\n", 544 delta * 1e3, sum * 1e3, data_size / (sum * (GB(1)))); 545 } 546 547 times[frame % countof(times)] = delta; 548 frame++; 549 } 550 i32 flag = beamformer_live_parameters_get_dirty_flag(); 551 if (flag != -1 && (1 << flag) == BeamformerLiveImagingDirtyFlags_StopImaging) 552 break; 553 } 554 555 lip.active = 0; 556 beamformer_set_live_parameters(&lip); 557 } else { 558 send_frame(data, &bp); 559 } 560 } 561 562 function void 563 sigint(i32 _signo) 564 { 565 g_should_exit = 1; 566 } 567 568 extern i32 569 main(i32 argc, char *argv[]) 570 { 571 Options options = parse_argv(argc, argv); 572 573 if (options.remaining_count != 1) 574 usage(argv[0]); 575 576 signal(SIGINT, sigint); 577 578 Arena arena = os_alloc_arena(KB(8)); 579 Stream path = stream_alloc(&arena, KB(4)); 580 stream_append_s8(&path, c_str_to_s8(options.remaining[0])); 581 582 execute_study(arena, path, &options); 583 584 return 0; 585 }