updates

examples/real_data_example/src/real_data_example.cc

@@ -1,14 +1,15 @@
 #include <mdio/mdio.h>
-#include <png.h>
 
-#include <algorithm>
 #include <indicators/cursor_control.hpp>
 #include <indicators/indeterminate_progress_bar.hpp>
 #include <indicators/progress_bar.hpp>
-#include <sstream>
-#include <thread>
 
 #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>{indicators::FontStyle::bold}}};
+          std::vector<indicators::FontStyle>{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': '<f4', 'fortran_order': False, 'shape': (" << width
-         << ", " << height << "), }";
-
-  // Pad header to multiple of 16 bytes
-  int header_len = header.str().length() + 1;  // +1 for null terminator
-  int padding =
-      16 - ((header_len + 10) % 16);  // 10 is magic + version + header_len
-  for (int i = 0; i < padding; i++) {
-    header << " ";
-  }
-  header << "\n";
-
-  // Write header length (little endian)
-  uint16_t header_size = header.str().length();
-
-  // Write all the header components
-  outfile.write(magic, 6);
-  outfile.write(version, 2);
-  outfile.write(reinterpret_cast<char*>(&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<const char*>(&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<uint8_t> 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<uint8_t>((val - min_val) * scale);
-    }
-    png_write_row(png, reinterpret_cast<png_bytep>(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();
 }