From c8bb196743a0e1bdc51604d186fcd3717301794d Mon Sep 17 00:00:00 2001 From: David Paul Graham <43794491+dpgraham4401@users.noreply.github.com> Date: Sun, 5 May 2024 00:07:30 -0400 Subject: [PATCH] Velocity models (#6) * add velocity model object and allocate_velocity_model utility function * .cc file to .cpp * new velocity testing file * set up separate files for velocity and libseis * simple velocity model function for generating quick velocity model with just a starting velocity and gradient --- .gitignore | 2 + CMakeLists.txt | 9 +++- src/libseis.c | 108 +++++++++++++++++++++-------------------- src/libseis.h | 43 ++++++++-------- src/velocity.c | 28 +++++++++++ src/velocity.h | 19 ++++++++ test/libseis_test.cc | 79 ------------------------------ test/libseis_test.cpp | 44 +++++++++++++++++ test/velocity_test.cpp | 57 ++++++++++++++++++++++ 9 files changed, 236 insertions(+), 153 deletions(-) create mode 100644 src/velocity.c create mode 100644 src/velocity.h delete mode 100644 test/libseis_test.cc create mode 100644 test/libseis_test.cpp create mode 100644 test/velocity_test.cpp diff --git a/.gitignore b/.gitignore index 3f00bc7..509633e 100644 --- a/.gitignore +++ b/.gitignore @@ -34,4 +34,6 @@ include/* lib/* bin/* test/test_runner +Testing +test/Temporary *.bin diff --git a/CMakeLists.txt b/CMakeLists.txt index 044efc5..8ee3856 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,12 +15,14 @@ set(SRC_DIR "src") ### create libseis library include_directories(src) -add_library(${LIB_NAME} SHARED ${SRC_DIR}/libseis.c) +add_library(${LIB_NAME} SHARED ${SRC_DIR}/libseis.c ${SRC_DIR}/velocity.c) + set_target_properties( ${LIB_NAME} PROPERTIES VERSION "v${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}.${PROJECT_VERSION_PATCH}" PUBLIC_HEADER ${SRC_DIR}/libseis.h + PUBLIC_HEADER ${SRC_DIR}/velocity.h ) install( @@ -41,7 +43,10 @@ set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) FetchContent_MakeAvailable(googletest) enable_testing() -add_executable(libseis_test test/libseis_test.cc) +add_executable(libseis_test test/libseis_test.cpp) +add_executable(velocity_test test/velocity_test.cpp) target_link_libraries(libseis_test ${LIB_NAME} GTest::gtest_main) +target_link_libraries(velocity_test ${LIB_NAME} GTest::gtest_main) include(GoogleTest) gtest_discover_tests(libseis_test) +gtest_discover_tests(velocity_test) diff --git a/src/libseis.c b/src/libseis.c index 6a6f078..9745872 100644 --- a/src/libseis.c +++ b/src/libseis.c @@ -1,75 +1,77 @@ #include "libseis.h" +#include #include #include #include -#include double *read_double(char path[], int nt, int nx) { - FILE *fptr; - double *data; - data = (double *) malloc(nt * nx * sizeof(double)); - fptr = fopen(path, "rb"); - if (fptr == NULL) { - printf("unable to open file at %s\n", path); - exit(EXIT_FAILURE); - } - for (int i = 0; i < nt * nx; i++) { - fread(&data[i], sizeof(double), 1, fptr); - } - fclose(fptr); - return data; + FILE *fptr; + double *data; + data = (double *)malloc(nt * nx * sizeof(double)); + fptr = fopen(path, "rb"); + if (fptr == NULL) { + printf("unable to open file at %s\n", path); + exit(EXIT_FAILURE); + } + for (int i = 0; i < nt * nx; i++) { + fread(&data[i], sizeof(double), 1, fptr); + } + fclose(fptr); + return data; } void write_double(char path[], double *data, int nt, int nx) { - FILE *fptr; - fptr = fopen(path, "wb"); - if (fptr == NULL) { - printf("unable to open file at %s\n", path); - exit(EXIT_FAILURE); - } - for (int i = 0; i < nx * nt; i++) { - fwrite(&data[i], sizeof(double), 1, fptr); - } - fclose(fptr); + FILE *fptr; + fptr = fopen(path, "wb"); + if (fptr == NULL) { + printf("unable to open file at %s\n", path); + exit(EXIT_FAILURE); + } + for (int i = 0; i < nx * nt; i++) { + fwrite(&data[i], sizeof(double), 1, fptr); + } + fclose(fptr); } /** - * Creates a new gather with the same data as the input gather, but with each sample - * multiplied by a power of its corresponding time value. + * Creates a new gather with the same data as the input gather, but with each + * sample multiplied by a power of its corresponding time value. * * @param gather The input gather to be multiplied. * @param power The power to raise the time values to. - * @return A new gather with the same dimensions as the input gather, but with modified data. + * @return A new gather with the same dimensions as the input gather, but with + * modified data. */ Gather *gain_gather(Gather *gather, double power) { - Gather *new_gather = malloc(sizeof(Gather)); - if (new_gather == NULL) { - fprintf(stderr, "Error: input gather is NULL\n"); - } - memcpy(new_gather, gather, sizeof(Gather)); - new_gather->data = malloc(gather->nt * gather->nx * sizeof(double)); - if (new_gather->data == NULL) { - fprintf(stderr, "Error: input gather data is NULL\n"); - } - memcpy(new_gather->data, gather->data, gather->nt * gather->nx * sizeof(double)); - for (int trace_index = 0; trace_index < gather->nx; trace_index++) { - for (int sample_index = 0; sample_index < gather->nt; sample_index++) { - double time = (double) sample_index * new_gather->dt; - double t_pow = pow(time, power); - int data_index = (trace_index * new_gather->nt) + sample_index; - new_gather->data[data_index] = gather->data[data_index] * t_pow; - } + Gather *new_gather = malloc(sizeof(Gather)); + if (new_gather == NULL) { + fprintf(stderr, "Error: input gather is NULL\n"); + } + memcpy(new_gather, gather, sizeof(Gather)); + new_gather->data = malloc(gather->nt * gather->nx * sizeof(double)); + if (new_gather->data == NULL) { + fprintf(stderr, "Error: input gather data is NULL\n"); + } + memcpy(new_gather->data, gather->data, + gather->nt * gather->nx * sizeof(double)); + for (int trace_index = 0; trace_index < gather->nx; trace_index++) { + for (int sample_index = 0; sample_index < gather->nt; sample_index++) { + double time = (double)sample_index * new_gather->dt; + double t_pow = pow(time, power); + int data_index = (trace_index * new_gather->nt) + sample_index; + new_gather->data[data_index] = gather->data[data_index] * t_pow; } - return new_gather; + } + return new_gather; } void display_gather(Gather *gather) { - printf("--------------------\n"); - printf("Gather: id %d\n", gather->id); - printf("Sample Rate: %f seconds\n", gather->dt); - printf("Sampling Frequency: %f Hz\n", 1 / gather->dt); - printf("Nyquist Frequency: %f Hz\n", 1 / (gather->dt * 2)); - printf("Samples per trace: %d\n", gather->nt); - printf("Number of traces: %d\n", gather->nx); - printf("--------------------\n"); + printf("--------------------\n"); + printf("Gather: id %d\n", gather->id); + printf("Sample Rate: %f seconds\n", gather->dt); + printf("Sampling Frequency: %f Hz\n", 1 / gather->dt); + printf("Nyquist Frequency: %f Hz\n", 1 / (gather->dt * 2)); + printf("Samples per trace: %d\n", gather->nt); + printf("Number of traces: %d\n", gather->nx); + printf("--------------------\n"); } diff --git a/src/libseis.h b/src/libseis.h index b27dd70..d5c3298 100644 --- a/src/libseis.h +++ b/src/libseis.h @@ -3,21 +3,23 @@ // David Graham - 2023-15-01 typedef enum { - SINGLE, - DOUBLE, + SINGLE, + DOUBLE, } PRECISION; /** * @brief A collection of seismic traces and attributes * - * @details A gather holds encapsulates all the information about seismic gathers, - * including the number of traces, the time sampling, number of time samples, and the data itself. - * It's intended to be generic enough to be used by various seismic gather types, such as CMP, or shot. + * @details A gather holds encapsulates all the information about seismic + * gathers, including the number of traces, the time sampling, number of time + * samples, and the data itself. It's intended to be generic enough to be used + * by various seismic gather types, such as CMP, or shot. * * @param * id: (int) The gather's unique identifier * @param - * nt: (int) number times samples in each trace, assumed homogeneous number of samples in each trace + * nt: (int) number times samples in each trace, assumed homogeneous number of + * samples in each trace * @param * nx: (int) number traces * @param @@ -25,15 +27,16 @@ typedef enum { * @param * data: (pointer) pointer to the array of samples * @param - * precision: (enum PRECISION) the precision of the samples (e.g., double, double) + * precision: (enum PRECISION) the precision of the samples (e.g., double, + * double) */ typedef struct Gather { - int id; - int nt; - int nx; - double dt; - double *data; - PRECISION precision; + int id; + int nt; + int nx; + double dt; + double *data; + PRECISION precision; } Gather; double *read_double(char path[], int nt, int nx); @@ -44,19 +47,21 @@ void write_double(char path[], double *data, int nt, int nx); * gain_gather * * @details - * Applies a time-variant scaling to a series of seismic traces in the same gather. - * It can be used to compensate for exponential decay (constant Q) in amplitude with time - * See "Tpow: an estimator of seismic amplitude decay" by Fowler & Claerbout - * nt' = nt * (dt * t^pow) + * Applies a time-variant scaling to a series of seismic traces in the same + * gather. It can be used to compensate for exponential decay (constant Q) in + * amplitude with time See "Tpow: an estimator of seismic amplitude decay" by + * Fowler & Claerbout nt' = nt * (dt * t^pow) * * @param - * gather: A pointer to a shot, CMP, or other seismic gather with all necessary fields. + * gather: A pointer to a shot, CMP, or other seismic gather with all necessary + * fields. * @param * power: The power (double) applied to the seismic trace as a function of time. * * @returns * A pointer to a new gather with the supplied Tpow applied. - * All gather attributes, besides the seismic data, are copied from the original gather + * All gather attributes, besides the seismic data, are copied from the original + * gather * */ Gather *gain_gather(Gather *gather, double power); diff --git a/src/velocity.c b/src/velocity.c new file mode 100644 index 0000000..b90f83d --- /dev/null +++ b/src/velocity.c @@ -0,0 +1,28 @@ +#include "velocity.h" + +#include + +VelocityModel *allocate_velocity_model(const int rows, const int cols) { + VelocityModel *model = malloc(sizeof(VelocityModel)); + model->rows = rows; + model->cols = cols; + + model->data = (double **)malloc(rows * sizeof(double *)); + for (int i = 0; i < rows; i++) { + model->data[i] = (double *)malloc(cols * sizeof(double)); + } + + return model; +} + +VelocityModel *simple_velocity_model(VelocityModel *model, + const double gradient, + const double startingVelocity) { + + for (int i = 0; i < model->rows; i++) { + for (int j = 0; j < model->cols; j++) { + model->data[i][j] = startingVelocity + gradient * i; + } + }; + return model; +} diff --git a/src/velocity.h b/src/velocity.h new file mode 100644 index 0000000..5b56f6c --- /dev/null +++ b/src/velocity.h @@ -0,0 +1,19 @@ +// +// Created by dg on 5/4/24. +// + +#ifndef VELOCITY_H +#define VELOCITY_H + +typedef struct { + int rows; + int cols; + double **data; +} VelocityModel; + +VelocityModel *allocate_velocity_model(int rows, int cols); + +VelocityModel *simple_velocity_model(VelocityModel *model, double gradient, + double startingVelocity); + +#endif // VELOCITY_H diff --git a/test/libseis_test.cc b/test/libseis_test.cc deleted file mode 100644 index fc091e1..0000000 --- a/test/libseis_test.cc +++ /dev/null @@ -1,79 +0,0 @@ -#include "gtest/gtest.h" - -extern "C" { -#include -#include "libseis.h" -} - -int const nt_test = 10; -int const nx_test = 1; - -double double_trace1[nt_test] = { - 0.000, - 2.000, - 4.000, - 8.000, - 1.000, - 2.000, - 3.000, - 4.000, - 5.000, - 6.000, -}; - -double cmp_array[nt_test * 2] = { - 0.000, - 2.000, - 4.000, - 8.000, - 1.000, - 2.000, - 3.000, - 4.000, - 5.000, - 6.000, - 0.000, - 0.000, - 1.000, - 0.020, - 0.301, - -0.101, - 0.035, - 0.011, - -0.006, - 0.009, -}; - -struct Gather cmp = { - 1, - nt_test, - 2, - 1, - cmp_array, - DOUBLE -}; - -TEST(LibseisIO, ReadWriteFloat) { - char tmp_file[] = "/tmp/data.bin"; - write_double(tmp_file, double_trace1, nt_test, nx_test); - double *data = read_double(tmp_file, nt_test, nx_test); - for (int i = 0; i < nt_test; i++) { - EXPECT_EQ(data[i], double_trace1[i]); - } -} - -TEST(GainCmp, GainsCmpGather) { - double t_pow = 2; - Gather *foo = gain_gather(&cmp, t_pow); - - printf("Check that the t^power function operates correctly"); - for (int i = 0; i < foo->nt; i++) { - printf("%f * (%f^%f) = %f\n", cmp.data[i], foo->dt * i, t_pow, foo->data[i]); - double t_powered = cmp.data[i] * pow(foo->dt * i, t_pow); - EXPECT_EQ(foo->data[i], t_powered); - } -} - -TEST(Gather, GatherDisplay) { - display_gather(&cmp); -} diff --git a/test/libseis_test.cpp b/test/libseis_test.cpp new file mode 100644 index 0000000..cabc828 --- /dev/null +++ b/test/libseis_test.cpp @@ -0,0 +1,44 @@ +#include "gtest/gtest.h" + +extern "C" { +#include "libseis.h" +#include +} + +int const nt_test = 10; +int const nx_test = 1; + +double double_trace1[nt_test] = { + 0.000, 2.000, 4.000, 8.000, 1.000, 2.000, 3.000, 4.000, 5.000, 6.000, +}; + +double cmp_array[nt_test * 2] = { + 0.000, 2.000, 4.000, 8.000, 1.000, 2.000, 3.000, 4.000, 5.000, 6.000, + 0.000, 0.000, 1.000, 0.020, 0.301, -0.101, 0.035, 0.011, -0.006, 0.009, +}; + +struct Gather cmp = {1, nt_test, 2, 1, cmp_array, DOUBLE}; + +TEST(LibseisIO, ReadWriteFloat) { + char tmp_file[] = "/tmp/data.bin"; + write_double(tmp_file, double_trace1, nt_test, nx_test); + double *data = read_double(tmp_file, nt_test, nx_test); + for (int i = 0; i < nt_test; i++) { + EXPECT_EQ(data[i], double_trace1[i]); + } +} + +TEST(GainCmp, GainsCmpGather) { + double t_pow = 2; + Gather *foo = gain_gather(&cmp, t_pow); + + printf("Check that the t^power function operates correctly"); + for (int i = 0; i < foo->nt; i++) { + printf("%f * (%f^%f) = %f\n", cmp.data[i], foo->dt * i, t_pow, + foo->data[i]); + double t_powered = cmp.data[i] * pow(foo->dt * i, t_pow); + EXPECT_EQ(foo->data[i], t_powered); + } +} + +TEST(Gather, GatherDisplay) { display_gather(&cmp); } \ No newline at end of file diff --git a/test/velocity_test.cpp b/test/velocity_test.cpp new file mode 100644 index 0000000..cffe957 --- /dev/null +++ b/test/velocity_test.cpp @@ -0,0 +1,57 @@ +#include "gtest/gtest.h" +extern "C" { +#include "velocity.h" +} + +TEST(VelocityModel, AllocateVelocityModel) { + constexpr int rows = 5; + constexpr int cols = 5; + + // Call the function to test + VelocityModel *model = allocate_velocity_model(rows, cols); + + // Check if the function returned a non-null pointer + ASSERT_NE(model, nullptr); + + // Check if the rows and cols are correctly set + EXPECT_EQ(model->rows, rows); + EXPECT_EQ(model->cols, cols); + + // Check if the data array is correctly allocated + ASSERT_NE(model->data, nullptr); + + // Check each row of the data array + for (int i = 0; i < rows; i++) { + ASSERT_NE(model->data[i], nullptr); + } +} + +TEST(SimpleVelocityModel, ShouldApplyCorrectGradient) { + int rows = 5; + int cols = 10; + double gradient = 2.0; + double startingVelocity = 1.0; + VelocityModel *model = allocate_velocity_model(rows, cols); + model = simple_velocity_model(model, gradient, startingVelocity); + for (int i = 0; i < model->rows; i++) { + for (int j = 0; j < model->cols; j++) { + ASSERT_EQ(model->data[i][j], startingVelocity + gradient * i); + } + } + free(model); +} + +TEST(SimpleVelocityModel, ShouldHandleZeroGradient) { + int rows = 5; + int cols = 10; + double gradient = 0.0; + double startingVelocity = 1.0; + VelocityModel *model = allocate_velocity_model(rows, cols); + model = simple_velocity_model(model, gradient, startingVelocity); + for (int i = 0; i < model->rows; i++) { + for (int j = 0; j < model->cols; j++) { + ASSERT_EQ(model->data[i][j], startingVelocity); + } + } + free(model); +} \ No newline at end of file