From 3218051dcf8637ed9c1299486e66c66cba33a9f7 Mon Sep 17 00:00:00 2001 From: blasscoc Date: Wed, 13 Nov 2024 16:07:32 -0600 Subject: [PATCH] updates --- .../src/real_data_example.cc | 150 +++--------------- 1 file changed, 21 insertions(+), 129 deletions(-) diff --git a/examples/real_data_example/src/real_data_example.cc b/examples/real_data_example/src/real_data_example.cc index c2fa213..6202a3a 100644 --- a/examples/real_data_example/src/real_data_example.cc +++ b/examples/real_data_example/src/real_data_example.cc @@ -1,14 +1,15 @@ #include -#include -#include #include #include #include -#include -#include #include "interpolation.h" +#include "seismic_numpy.h" +#include "seismic_png.h" +#include "tensorstore/tensorstore.h" + +#define MDIO_RETURN_IF_ERROR(...) TENSORSTORE_RETURN_IF_ERROR(__VA_ARGS__) using Index = mdio::Index; @@ -25,8 +26,9 @@ absl::Status Run() { // Select a inline slice mdio::SliceDescriptor desc1 = {"inline", inline_index, inline_index + 1, 1}; + mdio::SliceDescriptor desc2 = {"crossline", 500, 700, 1}; - MDIO_ASSIGN_OR_RETURN(auto inline_slice, dataset.isel(desc1)) + MDIO_ASSIGN_OR_RETURN(auto inline_slice, dataset.isel({desc1, desc2})) // Add seismic data reading example MDIO_ASSIGN_OR_RETURN(auto seismic_var, @@ -42,9 +44,11 @@ absl::Status Run() { indicators::option::PostfixText{"Loading seismic data..."}, indicators::option::ForegroundColor{indicators::Color::yellow}, indicators::option::FontStyles{ - std::vector{indicators::FontStyle::bold}}}; + std::vector{indicators::FontStyle::bold}}, + }; - // Hide cursor while bar is active + // Clear any existing output and hide cursor + std::cout << "\033[2K\r" << std::flush; indicators::show_console_cursor(false); // Start the async read operation @@ -61,7 +65,8 @@ absl::Status Run() { time_str << "Loading seismic data... " << elapsed << "s"; bar.set_option(indicators::option::PostfixText{time_str.str()}); bar.tick(); - std::this_thread::sleep_for(std::chrono::milliseconds(100)); + // Reduce update frequency slightly + std::this_thread::sleep_for(std::chrono::milliseconds(200)); } // Get the result and stop the bar @@ -88,128 +93,15 @@ absl::Status Run() { auto depth_exclusive_max = seismic_var.get_store().domain()[2].interval().exclusive_max(); - // Define these variables once at the start of the function, before both PNG - // and numpy sections - const int width = xline_exclusive_max - xline_inclusive_min; - const int height = depth_exclusive_max - depth_inclusive_min; - - // - - // Dump seismic data to numpy .npy format - std::ofstream outfile("seismic_data.npy", std::ios::binary); - if (!outfile) { - return absl::InvalidArgumentError("Could not open numpy file for writing"); - } - - // Write the numpy header - // Format described at: - // https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html - const char magic[] = "\x93NUMPY"; - const char version[] = "\x01\x00"; - - // Create the header string - std::stringstream header; - header << "{'descr': '(&header_size), 2); - outfile.write(header.str().c_str(), header_size); - - // Write the actual data - for (Index i = xline_inclusive_min; i < xline_exclusive_max; ++i) { - for (Index j = depth_inclusive_min; j < depth_exclusive_max; ++j) { - float val = seismic_accessor({inline_index, i, j}); - outfile.write(reinterpret_cast(&val), sizeof(float)); - } - } - outfile.close(); - - // For PNG section, keep using width/height - const float upsample = 2.0f; - const int output_width = width * upsample; - const int output_height = height * upsample; - - // Create PNG file - FILE* fp = fopen("seismic_slice.png", "wb"); - if (!fp) return absl::InvalidArgumentError("Could not open file for writing"); - - png_structp png = - png_create_write_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr); - if (!png) return absl::InternalError("Could not create PNG write struct"); - - png_infop info = png_create_info_struct(png); - if (!info) return absl::InternalError("Could not create PNG info struct"); - - if (setjmp(png_jmpbuf(png))) return absl::InternalError("PNG error occurred"); - - png_init_io(png, fp); - - // Set image attributes with 8-bit depth instead of 16 - png_set_IHDR(png, info, output_width, output_height, 8, PNG_COLOR_TYPE_GRAY, - PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, - PNG_FILTER_TYPE_DEFAULT); - png_write_info(png, info); - - // Update the row vector to use 8-bit values - std::vector row(output_width); - - // seismic wavelet has a zero mean .... - double mean = 0.0; - - // Calculate standard deviation - double sum_squared_diff = 0.0; - double count = 0.0; - for (Index i = xline_inclusive_min; i < xline_exclusive_max; ++i) { - for (Index j = depth_inclusive_min; j < depth_exclusive_max; ++j) { - float val = seismic_accessor({inline_index, i, j}); - double diff = val - mean; - sum_squared_diff += diff * diff; - count += 1.0; - } - } - double std_dev = std::sqrt(sum_squared_diff / count); - - // Set clipping bounds to ±2 standard deviations from mean - float min_val = mean - 3 * std_dev; - float max_val = mean + 3 * std_dev; - - // Update the scale factor for 8-bit range (0-255) - float scale = 255.0f / (max_val - min_val); - - for (int j = 0; j < output_height; ++j) { - float y = depth_inclusive_min + (j / upsample); - - for (int i = 0; i < output_width; ++i) { - float x = xline_inclusive_min + (i / upsample); - - float val = bilinear_interpolate( - seismic_accessor, x, y, inline_index, xline_inclusive_min, - xline_exclusive_max, depth_inclusive_min, depth_exclusive_max); - val = std::max(min_val, std::min(max_val, val)); - row[i] = static_cast((val - min_val) * scale); - } - png_write_row(png, reinterpret_cast(row.data())); - } + // Write numpy file + MDIO_RETURN_IF_ERROR(WriteNumpy(seismic_accessor, inline_index, + xline_inclusive_min, xline_exclusive_max, + depth_inclusive_min, depth_exclusive_max)); - // Cleanup - png_write_end(png, nullptr); - png_destroy_write_struct(&png, &info); - fclose(fp); + // Write PNG file + MDIO_RETURN_IF_ERROR(WritePNG(seismic_accessor, inline_index, + xline_inclusive_min, xline_exclusive_max, + depth_inclusive_min, depth_exclusive_max)); return absl::OkStatus(); }