open_segy_file4reading.m

实现地震勘探中

function [fid,seismic,param,parameters,text_header,binary_header] = ...
                               open_segy_file4reading(filename,varargin)
% Function opens an SEG-Y file, reads the textual header and the binary file 
% header and outputs the file ID, a preliminary seismic data structure 
% with empty fields "header" and "traces".
% In addition, it outputs parameters that allow function "read_segy_file_traces"
% to read consecutive traces.
%
% Written by: E. R.: January 6, 2007
% Last updated: February 27, 2008: bug fix in error message
%
%             [fid,seismic,param,parameters,text_header,binary_header] = ...
%                               open_segy_file4reading(filename,varargin)
% INPUT       (all input arguments are optional)
% filename    name of the  file to read
%             if the name is omitted, empty, or a file with this name is not found
%             a file selector box will pop up to allow interactive file selection
% varargin    Variable number of arguments. Each argument is a cell array whose first
%             element is a keyword and whose other elements can be strings, numeric
%             values, or other cell arrays
%             Possible keywords are:
%     format  floating point format to use; this cell array has the form
%             {'format',fp-format} where fp_format is one of the three strings:
%             'ibm'      IBM floating point format (standard and default)
%             'ieee'     IEEE format, big endian  (Sun, SGI, etc.)
%             'header'   The binary header bytes 25-26 are used to determine the format
%             Default: {'format','header'}; 
%     headers  header values to be read from binary trace header; this cell array has
%             the form {'headers',{mnem1,first,bytes,units,dimension,description}, ...
%             {mnem2,first,bytes,units,dimension,description},...} where
%             "mnem1", "mnem2", ... denote header mnemonics (such as CDP, OFFSET), 
%             "first" denotes the first byte in the binary header, 
%             "bytes" denotes the number of bytes occupied by the mnemonic (2 or 4), 
%             "units" denotes the units of measurement for the header value, and 
%             "description" is a description of the header value. 
%             Example:
%             {'headers',{'ILINE_NO',181,4,'n/a','CDP number'}, ...
%                        {'OFFSET',37,4,'m','Source-receiver distance'}}
%             See below for a list of headers retrieved by default.
%             Default: {'headers',{}}
%     times   times to output; this cell array has the form 
%             {'times',first,last}, {'times',[first,last]}, or {'times',[]} 
%             In the first two forma all samples with times between (and 
%             including) first and last (in ms) are output. In the last case
%             all samples are output.
%             Default: {'times',[]}
%     traces  select traces to output; this cell array has the form
%             {'traces',expression}; the variable "expression" is used to 
%             determine which traces to output. "expression" can be an index
%             vector specifying the traces to output
%             Examples: {'traces',1:2:100}
%                       {'traces',[2,5,7:10,22]}
%             Alternatively, it can be a string with a logical expression involving 
%             trace headers such as '10 <= cdp && 100 >= cdp'
%             Examples: {'traces','cdp == 100 && offset > 100'}
%                       {'traces','14000 < cdp  && (14660 >= cdp || 14680 <= cdp)'};
%             The variables in the logical relationships must be headers of the 
%             data set; Use of functions "fix", "mod", and "round" are permitted; all
%             other function names will be interpreted as headers and likely cause an
%             error; the case of the headers in an expression does not matter.
%             Default: {'traces',[]}
%    ignoreshift  By default this function reads byte locations 109-110 (see 
%             header "lag" below) and applies the shift to the seismic data;
%             This behavior can be overwritten by setting this parameter to true;
%             Default: {'ignoreshift',false}
%    max_mem   maximum amount of contiguous memory in megabytes (MB) bytes available 
%             to store seismic traces;
%             Default: {'max_mem',[]}
%             This means the maximum size is determined internally. 
%
% Headers retrieved by default are (any one of these headers is removed if it turns out to 
% be identically zero):
%        ds_seqno      Trace sequence number within line (1-4)
%        ffid          Original Field record number (9-12)
%        o_trace_no    Trace sequence number within original field record (13-16)
%        source        Energy source point number (17-20)
%        cdp           CDP ensemble number (21-24)
%        seq_cdp       Trace sequence number within CDP ensemble (25-28)
%        trc_type      Trace ID (1=live,2=dead,3=dummy,4=time break,...) (29-30)
%        offset        Distance from source point to receiver group (37-40)
%        rec_elev      Receiver elevation (41-44);
%        sou_elev      Surface elevation at source (45-48)
%        depth         Source depth below surface (49-52)
%        sou_h2od      Water depth at source (61-64)
%        rec_h2od      Water depth at receiver group (65-68)
%        sou_x         X coordinate of source (73-76)
%        sou_y         Y coordinate of source (77-80)
%        rec_x         X coordinate of receiver (81-84)
%        rec_y         Y coordinate of receiver (85-88)
%        lag           Lag time between shot and recording start in ms (109-110)
%                      (the value of lag is added to the start time of the 
%                      seismic; hence it can be used to simulate non-zero start
%                      time of the data)
%                      see also parameter "ignoreshift", above.
%        cdp_x         X coordinate of CDP (181-184)
%        cdp_y         Y coordinate of CDP (185-189)
%        iline_no      In-line number (189-192)
%        xline_no      Cross-line number (193-196)
%               The numbers in parentheses at the end of the line denote the location 
%               of the corresponding bytes in the SEG-Y trace header
%    
% OUTPUT
% seismic       Seismic structure
%       seismic.type               'seismic' (type of structure)
%       seismic.name               file name without exyension
%       seismic.from               Full name of the SEG-Y file
%       seismic.traces             Array of seismic traces
%       seismic.first              Start time of seismic (in ms)
%       seismic.last               End time of seismic (in ms)
%       seismic.step               Sample interval of seismic (in ms)
%       seismic.units              Time units used (ms)
%       seismic headers            Matrix with header mnemonics (one row 
%                                  per header)
%       seismic.header_info        Three-column cell array with header info 
%                                  (one row per header)
%       seismic.null               []
%       seismic.line_number        Line number (5-8)
%       seismic.reel_number        Reel number (9-12)
%       seismic.cdp_fold           CDP fold
%       seismic.traces_per_record  Data traces per record (13-14)
%       seismic.aux_per_record     Auxiliary traces per record (15-16)
%       seismic.history            A four element cell array. The first element
%                                  is the start date/time of the program that 
%                                  invoked this function; the second element is
%                                  the start date/time this function was executed;
%                                  and the last cell contains the name if the file 
%                                  that was read
%            Example:
%               seismic.offset     contains the offset for each trace
%               seismic.header_info.cdp two-element cell array {'m','Offset'}
%                                  the first element represents the units of measurement, the 
%                                  second is a description of the header
%                 
% ebcdic_header         EBCDIC reel header converted to ASCII
% binary_header         Binary reel header


% UPDARE HISTORY
%          July 15, 2007: bug fix in trace constraints

global ABORTED

%	Set default output arguments
run_presets_if_needed

%	Set default for input parameters
param.format='header';
param.headers={};
param.ignoreshift=false;
param.debug=false;
%param.precision='single';
param.times=[];
param.traces=[];
param.header_precision='single';
param.max_mem=[];   

if nargin == 0
   filename='';

else
%	Replace defaults by actual input arguments
   param=assign_input(param,varargin,'open_segy_file4reading');
   if ~ismember(lower(param.format),{'ieee','ibm','header'})
      disp([' Unknown trace format: ',param.format])
      disp(' Allowed formats are ''header'', ''ieee'', and ''ibm''.')
      drawnow
      error('Abnormal termination.')
   end
   if ~isempty(param.times)
      param=check_time_range_no0(param);  % Check if the time range is specified correctly
   end
end


%	Open SEG-Y file and get file ID
fid=open_segy_file_no1(filename);

if fid < 0
    seismic=[];
    text_header=[];
    binary_header=[];
    parameters=[];
   return
end
    
%	Read textual file header
text_header=read_textual_file_header_no2(fid);

%	Read binary file header
binary_header=read_binary_file_header_no3(fid);

%	Check FP format and save it and units of measurement for distance in 
%       structure "param"
param=check_file_headers_no4(binary_header,param);

if ABORTED
   seismic=[];
   return
end

%	Create seismic structure and index of the time samples to retain
[seismic,parameters.idx4times]=create_seismic_structure_no5(binary_header,param);

%	Collect info about the headers to read from binary trace header block
[parameters.header_info,parameters.indices,parameters.true4four, ...
       parameters.constraint_info,param]=select_trace_headers2read_no6(param);
if ABORTED
   seismic=[];
   return
end

param.nsamp=length(parameters.idx4times);
param.nheaders=length(parameters.indices);
param.no_samples=binary_header(8);

%	Get number of traces
param.ntraces=get_no_of_traces_no7(seismic.from,binary_header(8));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function param=check_time_range_no0(param)
% Check the compatibility of the time range specified

if ~isempty(param.times)
   if iscell(param.times)
      param.times=cell2num(param.times);
   end
   if length(param.times) == 1
      param.times=[param.times,param.times];
   elseif length(param.times) == 2
      if param.times(2) < param.times(1)
         error('Start time of time range is greater than end time.')
      end
      if param.times(1) < 0
         error('Start time of time range is less than 0.')
      end
   else
      error('Range of times to read must be specified by a start time and an end time.')
   end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fid=open_segy_file_no1(filename)
%	Open SEG=Y file and return file ID
% INPUT
% filenane filename; can be empty
% mformat  machine format (see Matlab function "fopen"
% OUPTUT
% fid      file ID

mformat='ieee-be';

%       Open the file
if ~isempty(filename)
   fid=fopen(filename,'r',mformat);
   if fid == -1 
      if ~isdeployed
         disp(['... unable to find file ',filename])
      end
   else
      filename2S4M(filename)
   end
else
   fid=-1;
end

if fid==-1     % Open file selector window
   [filename,ierr]=get_filename4r('sgy');
   if ierr
      return
   end
   fid=fopen(filename,'r',mformat);
   if fid < 0
      error(['File "',filename,'" could not be opened.'])
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function text_header=read_textual_file_header_no2(fid)
% Read EBCDIC header and output it as ASCII
% INPUT
% fid          File ID
% OUTPUT
% text_header  EBCDIC header as ASCII

global ABORTED S4M

%	Read EBCDIC reel header 1
text_header=fread(fid,3200,'uchar');
if isempty(text_header)
   if S4M.deployed
      errordlg('EBCDIC header is empty; requested file is either empty or not an SEG-Y file.', ...
      S4M.name)
      ABORTED=true;
      return
   else
      error('EBCDIC header is empty; requested file is either empty or not an SEG-Y file.')
   end
end

ABORTED=false;
text_header=char(ebcdic2ascii(reshape(text_header,80,40)'));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function binary_header=read_binary_file_header_no3(fid)
% Read binary header and initiate the seismic dataset
% INPUT
% fid   file ID
% OUTPUT
% binary_header   binary header (two-byte variables and four-byte variables
%                 are stored together in one vector

%	Read binary file header
bh=fread(fid,400,'uchar');

two_bytes=bh(1:2:399)*256+bh(2:2:400);
four_bytes=((bh(1:4:9)*256+bh(2:4:10))*256+bh(3:4:11))*256+bh(4:4:12);

binary_header=[four_bytes(1:3);two_bytes(7:200)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function param=check_file_headers_no4(binary_header,param)
%	Check binary-format code and output format and units of measurement in
%       structure param
% INPUT
% binary_header  binary header
% param      structure with input parameters for "read_segy_file"
% OUTPUT
% param      input structure "param" with theo new fields: "format" and "units"   

global ABORTED

%	Check format
if binary_header(10) == 1 
   if strcmpi(param.format,'header')  ||  strcmpi(param.format,'ibm')
%      param.fmt='ibm';
      param.format='ibm';
   else
      disp('IEEE format requested')
      mywarning(['Data apparently stored as 32-bit IBM floating point numbers; ' ...
         'data sample format code = ',num2str(binary_header(10))]);