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