1 #include <errno.h>
2 #include <string.h>
3
4 #include <segyio/segy.h>
5 #include "segyutil.h"
6
7 #include "matrix.h"
8 #include "mex.h"
9
mexFunction(int nlhs,mxArray * plhs[],int nrhs,const mxArray * prhs[])10 void mexFunction(int nlhs, mxArray *plhs[],
11 int nrhs, const mxArray *prhs[]) {
12
13 char* msg1;
14 char* msg2;
15 int err;
16
17 segy_file* fp = segyfopen( prhs[ 0 ], "rb" );
18 int first_trace = mxGetScalar( prhs[ 1 ] );
19 int last_trace = mxGetScalar( prhs[ 2 ] );
20 int notype = mxGetScalar( prhs[ 3 ] );
21
22 struct segy_file_format fmt = filefmt( fp );
23
24 char binary[ SEGY_BINARY_HEADER_SIZE ];
25 err = segy_binheader( fp, binary );
26 if( err != 0 ) {
27 msg1 = "segy:get_traces:binary";
28 msg2 = strerror( errno );
29 goto cleanup;
30 }
31
32 // if last_trace was defaulted we assign it to the last trace in the file
33 if( last_trace == -1 )
34 last_trace = fmt.traces - 1;
35
36 int traces = 1 + (last_trace - first_trace);
37 long long bufsize = (long long)fmt.samples * traces;
38
39 plhs[0] = mxCreateNumericMatrix( fmt.samples, traces, mxSINGLE_CLASS, mxREAL );
40 float* out = mxGetData( plhs[ 0 ] );
41
42 if( first_trace > last_trace ) {
43 msg1 = "segy:get_traces:bounds";
44 msg2 = "first trace must be smaller than last trace";
45 goto cleanup;
46 }
47
48 for( int i = first_trace; i <= last_trace; ++i ) {
49 err = segy_readtrace( fp, i, out, fmt.trace0, fmt.trace_bsize );
50 out += fmt.samples;
51
52 if( err != 0 ) {
53 msg1 = "segy:get_traces:segy_readtrace";
54 msg2 = strerror( errno );
55 goto cleanup;
56 }
57 }
58
59 segy_close( fp );
60
61 if( notype != -1 )
62 fmt.format = notype;
63
64 segy_to_native( fmt.format, bufsize, mxGetData( plhs[ 0 ] ) );
65
66 int interval;
67 segy_get_bfield( binary, SEGY_BIN_INTERVAL, &interval );
68 plhs[ 1 ] = mxCreateDoubleScalar( interval );
69 plhs[ 2 ] = mxCreateDoubleScalar( fmt.format );
70
71 return;
72
73 cleanup:
74 segy_close( fp );
75
76 mexErrMsgIdAndTxt( msg1, msg2 );
77 }
78