[segyio] 157/376: mex: get/put_ps_line, parse_ps_segycube
Jørgen Kvalsvik
jokva-guest at moszumanska.debian.org
Wed Sep 20 08:04:25 UTC 2017
This is an automated email from the git hooks/post-receive script.
jokva-guest pushed a commit to branch debian
in repository segyio.
commit 28336eca22c135d6ddcd2f3ee50e1d5a868f7d43
Author: Jørgen Kvalsvik <jokva at statoil.com>
Date: Thu Jan 12 08:59:28 2017 +0100
mex: get/put_ps_line, parse_ps_segycube
Matlab support for prestack files.
---
mex/CMakeLists.txt | 4 ++
mex/Segy.m | 105 ++++++++++++++++++++++++++++++--
mex/SegySpec.m | 4 ++
mex/segy_get_offsets_mex.c | 63 +++++++++++++++++++
mex/segy_get_ps_line_mex.c | 62 +++++++++++++++++++
mex/segy_read_write_ps_line_mex.c | 125 ++++++++++++++++++++++++++++++++++++++
tests/test_segy_mex.m | 50 ++++++++++++++-
7 files changed, 408 insertions(+), 5 deletions(-)
diff --git a/mex/CMakeLists.txt b/mex/CMakeLists.txt
index 8c406bc..44fe589 100644
--- a/mex/CMakeLists.txt
+++ b/mex/CMakeLists.txt
@@ -21,6 +21,7 @@ configure_file(TraceField.m TraceField.m)
mexo(segyutil)
mex(segyspec_mex)
mex(segy_read_write_line_mex segyutil)
+mex(segy_read_write_ps_line_mex segyutil)
mex(segy_get_header_mex segyutil)
mex(segy_get_traces_mex segyutil)
mex(segy_put_traces_mex segyutil)
@@ -30,10 +31,12 @@ mex(segy_get_bfield_mex segyutil)
mex(segy_get_trace_header_mex segyutil)
mex(segy_get_field_mex segyutil)
mex(segy_put_headers_mex segyutil)
+mex(segy_get_offsets_mex segyutil)
install(FILES
${CMAKE_CURRENT_BINARY_DIR}/segyspec_mex.mexa64
${CMAKE_CURRENT_BINARY_DIR}/segy_read_write_line_mex.mexa64
+ ${CMAKE_CURRENT_BINARY_DIR}/segy_read_write_ps_line_mex.mexa64
${CMAKE_CURRENT_BINARY_DIR}/segy_get_header_mex.mexa64
${CMAKE_CURRENT_BINARY_DIR}/segy_get_traces_mex.mexa64
${CMAKE_CURRENT_BINARY_DIR}/segy_put_traces_mex.mexa64
@@ -43,6 +46,7 @@ install(FILES
${CMAKE_CURRENT_BINARY_DIR}/segy_get_trace_header_mex.mexa64
${CMAKE_CURRENT_BINARY_DIR}/segy_get_field_mex.mexa64
${CMAKE_CURRENT_BINARY_DIR}/segy_put_headers_mex.mexa64
+ ${CMAKE_CURRENT_BINARY_DIR}/segy_get_offsets_mex.mexa64
SegySpec.m
Segy.m
SegySampleFormat.m
diff --git a/mex/Segy.m b/mex/Segy.m
index 0829565..d032b99 100644
--- a/mex/Segy.m
+++ b/mex/Segy.m
@@ -182,16 +182,113 @@ classdef Segy
segycube = SegySpec(filename, il_word, xl_word, t0);
end
+ % Goal:
+ % Interpret segy cube as a 3D cube and save information needed to access
+ % the segy file as a cube in terms of inline and crossline numbers.
%
+ % inputs:
+ % name meaning
+ % filename filename of segy file
+ % offset bytenumber or header word for offset number
+ % il_word bytenumber or header word for inline number
+ % xl_word bytenumber or header word for crossline number
+ % t0 Time (ms) / depth (m) of first sample. Optional (default = 0)
+ %
+ % output:
+ % name meaning
+ % segycube struct needed to access segy file as a cube. Used by
+ % get_ps_line.m, put_ps_line.m and possibly friends
+ function segycube = parse_ps_segycube(filename, offset, il_word, xl_word, t0)
+ if ~exist('offset', 'var') || isempty(offset)
+ offset = TraceField.offset;
+ end
+ if ~exist('il_word', 'var') || isempty(il_word)
+ il_word = TraceField.Inline3D;
+ end
+ if ~exist('xl_word', 'var') || isempty(xl_word)
+ xl_word = TraceField.Crossline3D;
+ end
+ if ~exist('t0', 'var') || isempty(t0)
+ t0 = 0;
+ end
+
+ if ischar(il_word)
+ il_word = TraceField.(il_word);
+ end
+
+ if ischar(xl_word)
+ xl_word = TraceField.(xl_word);
+ end
+
+ if ischar(offset)
+ offset = TraceField.(offset);
+ end
+
+ offset = int32(offset);
+ il_word = int32(il_word);
+ xl_word = int32(xl_word);
+ t0 = double(t0);
+
+ segycube = Segy.interpret_segycube(filename, il_word, xl_word, t0);
+ offsets = segy_get_offsets_mex(segycube, offset, il_word, xl_word);
+ segycube.offset = offsets;
+ end
+
% Goal:
- % Read a full cube from segycube struct.
+ % Read an inline / crosline from a cube.
%
% Inputs:
- % sc Data as an interpreted segy cube from
+ % cube Data as an interpreted segy cube from
% 'interpret_segycube.m'
- % Output:
- % data Extracted cube with same sorting as in segy.
+ % dir Direction of desired line (iline / xline) as a string
+ % n Inline / crossline number
%
+ % Output:
+ % data Extracted line
+ function data = get_ps_line(cube, dir, n)
+ if nargin < 3
+ error('Too few arguments. Usage: Segy.get_ps_line(cube, dir, n)');
+ end
+
+ if strcmpi(dir, 'iline')
+ len = max(size(cube.crossline_indexes));
+ ix = cube.inline_indexes;
+ st = cube.il_stride;
+ elseif strcmpi(dir, 'xline')
+ len = max(size(cube.inline_indexes));
+ ix = cube.crossline_indexes;
+ st = cube.xl_stride;
+ else
+ error('Only iline and xline are valid directions.');
+ end
+
+ tmp = segy_read_write_ps_line_mex(cube, n, len, ix, st);
+ tmp = reshape(tmp, [], cube.offset_count);
+ nt = length(cube.t);
+ data = permute(reshape(tmp, nt, size(tmp, 1)/nt, []), [1 3 2]);
+ end
+
+ function data = put_ps_line(cube, data, dir, n)
+ if nargin < 4
+ error('Too few arguments. Usage: Segy.put_ps_line(cube, data, dir, n)');
+ end
+
+ if strcmpi(dir, 'iline')
+ len = max(size(cube.crossline_indexes));
+ ix = cube.inline_indexes;
+ st = cube.il_stride;
+ elseif strcmpi(dir, 'xline')
+ len = max(size(cube.inline_indexes));
+ ix = cube.crossline_indexes;
+ st = cube.xl_stride;
+ else
+ error('Only iline and xline are valid directions.');
+ end
+
+ tmp = permute( data, [1 3 2] );
+ segy_read_write_ps_line_mex( cube, n, len, ix, st, tmp );
+ end
+
function data = get_cube(sc)
data = Segy.get_traces(sc.filename);
diff --git a/mex/SegySpec.m b/mex/SegySpec.m
index 35d49a6..a30eec2 100644
--- a/mex/SegySpec.m
+++ b/mex/SegySpec.m
@@ -8,6 +8,7 @@ classdef SegySpec
crossline_indexes
inline_indexes
offset_count
+ offset
first_trace_pos
il_stride
xl_stride
@@ -15,6 +16,7 @@ classdef SegySpec
t
iline
xline
+ sorting
end
methods
@@ -42,6 +44,8 @@ classdef SegySpec
obj.t = obj.sample_indexes;
obj.iline = obj.inline_indexes;
obj.xline = obj.crossline_indexes;
+ [~, s] = enumeration(obj.trace_sorting_format);
+ obj.sorting = s{uint32(obj.trace_sorting_format) + 1};
end
end
diff --git a/mex/segy_get_offsets_mex.c b/mex/segy_get_offsets_mex.c
new file mode 100644
index 0000000..00713ad
--- /dev/null
+++ b/mex/segy_get_offsets_mex.c
@@ -0,0 +1,63 @@
+#include <string.h>
+#include "mex.h"
+#include "matrix.h"
+
+#include "segyutil.h"
+
+#include <segyio/segy.h>
+
+void mexFunction(int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[]) {
+
+ int errc = SEGY_OK;
+ if (nrhs != 4) {
+ goto ERROR;
+ }
+
+ const mxArray* mx_spec = prhs[0];
+ const mxArray* mx_offset = prhs[1];
+ const mxArray* mx_il_word = prhs[2];
+ const mxArray* mx_xl_word = prhs[3];
+
+ SegySpec spec;
+ recreateSpec(&spec, mx_spec);
+
+ int offset = (int)mxGetScalar(mx_offset);
+ int il = (int)mxGetScalar(mx_il_word);
+ int xl = (int)mxGetScalar(mx_xl_word);
+
+ segy_file* fp = segy_open( spec.filename, "rb" );
+
+ if( !fp ) {
+ errc = SEGY_FOPEN_ERROR;
+ goto ERROR;
+ };
+
+ plhs[0] = mxCreateNumericMatrix(1, spec.offset_count, mxINT32_CLASS, mxREAL);
+ int* int_offsets = mxMalloc(sizeof( int ) * spec.offset_count);
+
+ errc = segy_offset_indices(fp, offset, spec.offset_count,
+ int_offsets,
+ spec.first_trace_pos, spec.trace_bsize);
+
+ if( errc != SEGY_OK ) goto CLEANUP;
+
+ int32_t* plhs0 = (int32_t*)mxGetData(plhs[0]);
+ for( unsigned int i = 0; i < spec.offset_count; ++i )
+ plhs0[i] = int_offsets[i];
+
+ mxFree( int_offsets );
+ segy_close(fp);
+ return;
+
+ CLEANUP:
+ segy_close(fp);
+ ERROR:
+ {
+ int nfields = 1;
+ const char *fnames[nfields];
+ fnames[0] = "error";
+ plhs[0] = mxCreateStructMatrix(0,0, nfields, fnames);
+ mxSetFieldByNumber(plhs[0], 0, 0, mxCreateDoubleScalar(errc));
+ }
+}
diff --git a/mex/segy_get_ps_line_mex.c b/mex/segy_get_ps_line_mex.c
new file mode 100644
index 0000000..cafdf43
--- /dev/null
+++ b/mex/segy_get_ps_line_mex.c
@@ -0,0 +1,62 @@
+#include <string.h>
+#include "mex.h"
+#include "matrix.h"
+
+#include "segyutil.h"
+
+#include <segyio/segy.h>
+
+void mexFunction(int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[]) {
+
+ int errc = 0;
+ if (nrhs != 4) {
+ goto ERROR;
+ }
+
+ const mxArray* mx_spec = prhs[0];
+ const mxArray* mx_offset = prhs[1];
+ const mxArray* mx_il_word = prhs[2];
+ const mxArray* mx_xl_word = prhs[3];
+
+ SegySpec spec;
+ recreateSpec(&spec, mx_spec);
+
+ int offset = (int)mxGetScalar(mx_offset);
+ int il = (int)mxGetScalar(mx_il_word);
+ int xl = (int)mxGetScalar(mx_xl_word);
+
+ segy_file* fp = segyio_open( spec.filename, "rb" );
+
+ if( !fp ) {
+ errc = SEGY_FOPEN_ERROR;
+ goto ERROR;
+ };
+
+ plhs[0] = mxCreateNumericMatrix(spec.sample_count, spec.offset_count, mxINT32_CLASS, mxREAL);
+ int* int_offsets = mxMalloc(sizeof( int ) * spec.offset_count);
+
+ errc = segy_offsets(fp, offset, spec.offset_count,
+ int_offsets,
+ spec.first_trace_pos, spec.trace_bsize);
+
+ if( err != SEGY_OK ) goto CLEANUP;
+
+ int32_t* plhs0 = (int32_t*)mxGetScalar(plhs[0]);
+ for( unsigned int i = 0; i < spec.sample_count; ++i )
+ plhs0[i] = int_offsets[i];
+
+ segy_close(fp);
+ return;
+
+ CLEANUP:
+ segy_close(fp);
+ ERROR:
+ {
+ int nfields = 1;
+ const char *fnames[nfields];
+ fnames[0] = "error";
+ plhs[0] = mxCreateStructMatrix(0,0, nfields, fnames);
+ mxSetFieldByNumber(plhs[0], 0, 0, mxCreateDoubleScalar(errc));
+ }
+}
diff --git a/mex/segy_read_write_ps_line_mex.c b/mex/segy_read_write_ps_line_mex.c
new file mode 100644
index 0000000..4a084df
--- /dev/null
+++ b/mex/segy_read_write_ps_line_mex.c
@@ -0,0 +1,125 @@
+#include <string.h>
+#include "mex.h"
+#include "matrix.h"
+
+#include "segyutil.h"
+
+#include <segyio/segy.h>
+
+/* The gateway function */
+void mexFunction(int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[]) {
+
+ if( nrhs < 5 || nrhs > 6 ) goto ERROR;
+ bool read = nrhs == 5;
+
+ const mxArray* mx_spec = prhs[0];
+ const mxArray* mx_index = prhs[1];
+ const mxArray* mx_line_length = prhs[2];
+ const mxArray* mx_line_indexes = prhs[3];
+ const mxArray* mx_stride = prhs[4];
+
+ SegySpec spec;
+ recreateSpec(&spec, mx_spec);
+
+ size_t index = (size_t)mxGetScalar(mx_index);
+
+ uint32_t line_length = (uint32_t)mxGetScalar(mx_line_length);
+
+ uint32_t* line_indexes = (uint32_t*)mxGetData(mx_line_indexes);
+ int n = mxGetN(mx_line_indexes);
+ int m = mxGetM(mx_line_indexes);
+ uint32_t line_count = (n>m)? n:m;
+
+ uint32_t stride = (uint32_t)mxGetScalar(mx_stride);
+ int offsets = spec.offset_count;
+
+ segy_file* fp;
+
+ unsigned int line_trace0;
+ int errc = segy_line_trace0( index, line_length, stride, offsets, line_indexes, line_count, &line_trace0 );
+ if( errc != SEGY_OK ) {
+ mexErrMsgIdAndTxt( "segy:get_ps_line:wrong_line_number",
+ "Specified line number is not in cube." );
+ return;
+ }
+
+ const char* mode = read ? "rb" : "r+b";
+ fp = segy_open( spec.filename, mode );
+ if( !fp ) {
+ mexErrMsgIdAndTxt( "segy:get:ps_line:file",
+ "unable to open file" );
+ return;
+ }
+
+ mwSize dims[] = { spec.sample_count, spec.offset_count, line_length };
+
+ const size_t tr_size = SEGY_TRACE_HEADER_SIZE + spec.trace_bsize;
+
+ if( read ) {
+ plhs[0] = mxCreateNumericArray(3, dims, mxSINGLE_CLASS, mxREAL );
+ float* buf = (float*) mxGetData(plhs[0]);
+
+ for( int i = 0; i < offsets; ++i ) {
+ errc = segy_read_line( fp,
+ line_trace0,
+ line_length,
+ stride,
+ offsets,
+ buf + (spec.sample_count * line_length * i),
+ spec.first_trace_pos + (i * tr_size),
+ spec.trace_bsize );
+
+ if( errc != 0 ) goto CLEANUP;
+ }
+
+ errc = segy_to_native( spec.sample_format,
+ offsets * line_length * spec.sample_count,
+ buf );
+
+ if( errc != 0 ) goto CLEANUP;
+ }
+ else {
+ const mxArray* mx_data = prhs[5];
+ float* buf = (float*) mxGetData(mx_data);
+
+ errc = segy_from_native( spec.sample_format,
+ offsets * line_length * spec.sample_count,
+ buf );
+
+ if( errc != 0 ) goto CLEANUP;
+
+ for( int i = 0; i < offsets; ++i ) {
+ errc = segy_write_line( fp,
+ line_trace0,
+ line_length,
+ stride,
+ offsets,
+ buf + (spec.sample_count * line_length * i),
+ spec.first_trace_pos + (i * tr_size),
+ spec.trace_bsize );
+ if( errc != 0 ) goto CLEANUP;
+ }
+
+ errc = segy_to_native( spec.sample_format,
+ offsets * line_length * spec.sample_count,
+ buf );
+ if( errc != 0 ) goto CLEANUP;
+ }
+
+ segy_close(fp);
+ return;
+
+ CLEANUP:
+ segy_close(fp);
+ ERROR:
+ {
+ int nfields = 1;
+ const char *fnames[nfields];
+ fnames[0] = "error";
+ plhs[0] = mxCreateStructMatrix(0,0, nfields, fnames);
+ mxSetFieldByNumber(plhs[0], 0, 0, mxCreateDoubleScalar(errc));
+ }
+
+
+}
diff --git a/tests/test_segy_mex.m b/tests/test_segy_mex.m
index affff32..021d6f1 100644
--- a/tests/test_segy_mex.m
+++ b/tests/test_segy_mex.m
@@ -54,7 +54,7 @@ spec = Segy.interpret_segycube(filename, 'Inline3D', 'Crossline3D', t0);
data = Segy.get_line(spec, 'iline', 4);
sample_count = length(spec.sample_indexes);
-eps = 1e-4;
+eps = 1e-6;
% first trace along xline
% first sample
assert(abs(data(1, 1) - 4.2)<eps);
@@ -228,3 +228,51 @@ Segy.get_header( filename, 'cdp' );
increasing = linspace( 1, notraces, notraces );
Segy.put_headers( filename_copy, increasing, 'offset' );
assert( all(increasing == Segy.get_header( filename_copy, 'offset' )) );
+
+%% test_traces_offset
+prestack_filename = 'test-data/small-ps.sgy';
+cube_with_offsets = Segy.parse_ps_segycube(prestack_filename);
+assert(isequal(cube_with_offsets.offset, [1 2]));
+
+try
+ ps_line1 = Segy.get_ps_line(cube_with_offsets, 'iline', 22);
+ % reading a line that doesn't exist should throw
+ assert(false);
+end
+ps_line1 = Segy.get_ps_line(cube_with_offsets, 'iline', 1);
+ps_line2 = Segy.get_ps_line(cube_with_offsets, 'xline', 1);
+
+% offset 1
+assert(abs(ps_line1(1,1,1) - 101.01) < eps);
+assert(abs(ps_line1(1,1,2) - 101.02) < eps);
+assert(abs(ps_line1(1,1,3) - 101.03) < eps);
+assert(abs(ps_line1(2,1,1) - 101.01001) < eps);
+assert(abs(ps_line1(3,1,1) - 101.01002) < eps);
+assert(abs(ps_line1(2,1,2) - 101.02001) < eps);
+
+% offset 2
+assert(abs(ps_line1(1,2,1) - 201.01) < eps);
+assert(abs(ps_line1(1,2,2) - 201.02) < eps);
+assert(abs(ps_line1(1,2,3) - 201.03) < eps);
+assert(abs(ps_line1(2,2,1) - 201.01001) < eps);
+assert(abs(ps_line1(3,2,1) - 201.01002) < eps);
+assert(abs(ps_line1(2,2,2) - 201.02001) < eps);
+
+% test writing of prestack lines
+Segy.put_ps_line( cube_with_offsets, ps_line1 - 100.00, 'iline', 1 );
+% offset 1
+wr_line1 = Segy.get_ps_line( cube_with_offsets, 'iline', 1 );
+assert(abs(wr_line1(1,1,1) - 001.01) < eps);
+assert(abs(wr_line1(1,1,2) - 001.02) < eps);
+assert(abs(wr_line1(1,1,3) - 001.03) < eps);
+assert(abs(wr_line1(2,1,1) - 001.01001) < eps);
+assert(abs(wr_line1(3,1,1) - 001.01002) < eps);
+assert(abs(wr_line1(2,1,2) - 001.02001) < eps);
+
+% offset 2
+assert(abs(wr_line1(1,2,1) - 101.01) < eps);
+assert(abs(wr_line1(1,2,2) - 101.02) < eps);
+assert(abs(wr_line1(1,2,3) - 101.03) < eps);
+assert(abs(wr_line1(2,2,1) - 101.01001) < eps);
+assert(abs(wr_line1(3,2,1) - 101.01002) < eps);
+assert(abs(wr_line1(2,2,2) - 101.02001) < eps);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/segyio.git
More information about the debian-science-commits
mailing list