s_cplot.m

实现地震勘探中

function aux=s_cplot(seismic,varargin)
% Plot seismic data in color-coded form.
%
% Written by: E. R.: May 16, 2000
% Last updated: April 22, 2009: bug fix if number of traces is 1
%
%             aux=s_cplot(seismic,varargin)
% INPUT
% seismic     seismic structure
% varargin    one or more cell arrays; the first element of each cell array
%             is a keyword, the other elements are parameters.
%             Presently, keywords are:
%     'annotation'  header mnemonic to use for the horizontal axis
%             Default: {'annotation','trace_no'}
%     'clim'  two numbers representing fractions of the difference between smallest and 
%             largest sample  value (see also keyword 'limits') to compute lower and 
%             upper limits limit of values to which colors will be assigned.
%             limits(1)=min(seismic.traces) + (max(seismic.traces)-min(seismic.traces)*clim(1)
%             limits(2)=max(seismic.traces) - (max(seismic.traces)-min(seismic.traces)*clim(2)
%             Default: {'clim',0.3333,0.3333}
%     'colorbar' plot colorbar; possible values are 'yes' and 'no'.
%             Default: {'colorbar','yes'}
%     'colormap'  colormap to use to map seismic amplitudes to color. Can be 
%             any of those predefined in MATLAB (such as 'copper','hot'); 
%             type "help graph3d" in MATLAB for a list of color maps. 
%             Also {'colormap','gray'} creates a gray-scale color map 
%             (the smaller the value the darker the color)
%             Default: {colormap,'seismic'}; this creates a blue-to-red colormap 
%                                            (function "mycolormap('seismic')")
%     'direction'  plot direction. Possible values are: left-to-right, 'l2r', 
%             and right-to-left, 'r2l'.
%             Default: {'direction','l2r') 
%     'figure'   Specifies if new figure should be created or if the seismic 
%             traces should be plotted to an existing figure. 
%             Possible values are 'new' and any other string. 
%             Default: {'figure','new'} 
%     'flip_colormap' keyword to flip the colormap (e.g. the default "blue-to-red" 
%             colormap is converted to "red-to-blue"); possible values are 'yes'
%             and any other string.
%             Default: {'flip_colormap','no'}
%     'imagemenu'  Specifies if figure should have a menu button to allow 
%             interactive change of color parameters. 
%             Possible values are: 'yes' and 'no'
%             Default: {'imagemenu','yes'}
%     'interpol'    2-element cell array {'interpol','no'|'cubic'}
%             Default: {'interpol','cubic'}
%     'limits' lower and upper limit of values to which colors will be assigned; 
%             if given it overrides parameters specified via keyword 'clim'. 
%             Default: {'limits',[],[]}
%                      specifies how data should be interpolated in time
%     'orient' Plot orientation. Possible values are: 'portrait' and 'landscape'
%             Default: {'orient','landscape'}
%     'polarity'    2-element cell array. Possible values are 1 and -1;
%             Default: {'polarity',1}
%     'scale'       2-element cell array which specifies if individual traces
%             should be scaled relative to one another. There are two scaling options: 
%             'median'  makes the median of the absolute values of each trace the same. 
%             'max'    makes the maximum of the absolute values of each trace the same. 
%             The third alternative is 'no', meaning data are displayed "as is".
%             The previously used option 'yes' is the same as 'median' and, 
%             while deprecated, still works.
%             Default: {'scale','no'}
%     'shading'  color shading (see Matlab function "shading")
%             possible values are '', 'flat', 'faceted', 'interp'
%             if not empty function "pcolor" is used to create the plot
%             with the type of shading specified. Shading 'interp" tends to 
%             give a more continuous display.
%             Otherwise, function imagesc" is used. This is faster.
%             Default: {'shading',''}
%     'times'       2-element or 3-element cell array 
%             {'times',vector of first and last time to plot} or ('times',first,last}. 
%             Default: {'times',seismic.first,seismic.last} which is
%                      equivalent to {'times',[seismic.first,seismic.last]}
%     'time_lines'  Two-element or three-element cell array. the second and
%             third element indicate time intervals at multiples of which 
%             time lines will be plotted. The time lines of the third
%             element (if given), are thicker than those of the first.
%             No time lines are plotted if {'time_lines',[]}.
%             Default: {'time_lines',500,1000}   
%     'title'       2-element cell array. The second element is a plot title.
%             Default: {'title',seismic.name}
%     'traces'      2-element or 3-element cell array. The second element can 
%             be an array of trace numbers or it can be a string. If it is a  
%             string it can be a header mnemonic or it can contain a logical
%             expression involving header values to include.
%             The "pseudo-header" 'trace_no' can also be used.
%             If the second element is a string containing a header mnemonic 
%             there must be a third element containing a vector of values 
%             (see "s_select")
%             Default:  {'traces',[]} which is equivalent to 
%                       {'traces',1:ntr} where ntr denotes the number of traces  
%                           in the input data set (ntr = size(seismic.traces,2))
%     'tracking' track cursor position; possible values are 'yes', 'no', and ''.
%             In the latter case a tracking button is created if the the
%             seismic is plotted in a new figure. Otherwise it is not.
%             Default: {'tracking',''}
% OUTPUT
% aux         structure with information about the plot
%     'figure_handle'    handle of the figure with this plot
%
% EXAMPLE
%          seismic=s_data;
%          s_cplot(seismic,{'limits',-1,1},{'colormap','copper'})
%          s_cplot(seismic,{'limits',-1,1},{'colormap','copper'},{'shading','faceted'})

% UPDATE HISTORY
%          July 25, 2007: Use a function for subset creation


if ~istype(seismic,'seismic')
   if isnumeric(seismic)
      seismic=s_convert(seismic,1,1);
      seismic.units='Samples';
   else
      error('First input argument must be a seismic dataset or a matrix.')
   end
end

run_presets_if_needed

nsamp=size(seismic.traces,1);
if nsamp == 1
   disp('Only one sample per trace; data set not plotted')
   return
end

%     Set default valuesz
param.annotation='trace_no';
param.clim={0.3333,0.3333};
param.colorbar='yes';
param.colormap='seismic';
param.direction='l2r';
param.figure='new';
param.figure_only='no';
param.flip_colormap='no';
param.flip_colors=[];   % Deprecated
param.imagemenu='yes';
param.interpol='v5cubic';
param.limits=cell(1,2);
param.npixels=1000;
param.orient='landscape';
param.polarity=1;
param.scale='no';
param.shading='';
param.time_lines={500,1000};
param.times=[];
param.title=seismic.name;
param.traces=[];
param.tracking='';

%       Decode input arguments
param=assign_input(param,varargin,'s_cplot');

if ~isempty(param.flip_colors)
   alert(' Keyword "flip_colors" has been replaced by "flip_colormap".')
   param.flip_colormap=param.flip_colors;
end

%	Select subset of seismic data
seismic=select_subset_of_seismic_traces(seismic,param.times,param.traces);

%{
if ~isempty(param.traces)  || ~isempty(param.times)
   history=S4M.history;       % Preserve value of global variable S4M.history
   S4M.history=false;

   if iscell(param.traces)
      param.traces=cell2num(param.traces);
   end

   if isempty(param.times)
      seismic=s_select(seismic,{'traces',param.traces});

   elseif iscell(param.times)
      seismic=s_select(seismic,{'traces',param.traces},{'times',param.times{1},param.times{2}});

   else
      seismic=s_select(seismic,{'traces',param.traces},{'times',param.times(1),param.times(2)});

   end

   S4M.history=history;		% Restore history setting
end
%}


[nsamp,ntr]=size(seismic.traces);
if nsamp == 1
   disp('Only one sample per trace; data set not plotted')
   return
end

%     Change polarity if necessary
if param.polarity < 0
   seismic.traces=-seismic.traces;
end

%     Interpolate data in time if necessary
if strcmpi(param.interpol,'cubic')  &&  nsamp < param.npixels
   npix=round(param.npixels/(nsamp-1))*(nsamp-1);
   dti=(seismic.last-seismic.first)/npix;
   times=(seismic.first:dti:seismic.last)';

   if isnull(seismic)   % Turn off warnings caused by NaN's in seismic traces
      warning('off','MATLAB:interp1:NaNinY')
      yi=interp1(seismic.first:seismic.step:seismic.last,seismic.traces,times,param.interpol);
      warning('on','MATLAB:interp1:NaNinY')
   else
      yi=interp1(seismic.first:seismic.step:seismic.last,seismic.traces,times,param.interpol);
   end

else
   dti=seismic.step;
   yi=seismic.traces;
end


%     Compute horizontal trace locations
xi=s_gh(seismic,param.annotation);
if min(xi) == max(xi)  &&  ntr  >  1
   error([' Header requested for annotation (',param.annotation,') is constant'])
end

%       Check if header values change uniformly from one trace to the next.
if ntr == 1 || isconstant(diff(xi),0.001)
   text2append='';
else
   text2append=' - approximate';
   alert(['Trace annotation with header "',param.annotation,'" is only an approximation.'])
end

if strcmpi(param.figure,'new')
   if strcmpi(param.orient,'portrait')
      figure_handle=pfigure;
   else
      figure_handle=lfigure;
   end
else
   figure_handle=gcf;
end
bgGray		% Gray background

if nargout > 0
   aux.figure_handle=figure_handle;
end
set(figure_handle,'DoubleBuffer','on')
if ntr > 1
   dxi=min(diff(xi))*0.5;
else
   dxi=0.5;
end

if isyes(param.colorbar)
   colorbar
end

axis([min(xi)-dxi,max(xi)+dxi,seismic.first,seismic.last])
ha=get(figure_handle,'CurrentAxes');
set(ha,'TickDir','out','YDir','reverse','XAxisLocation','top')
hold on

%    Handle reversal of plot direction
if strcmpi(param.direction,'r2l')
   set(ha,'xdir','reverse')
%   yi=fliplr(yi);
elseif ~strcmpi(param.direction,'l2r')
   error(['Keyword for plot direction is wrong (',param.direction,')'])
end

if strcmpi(param.figure_only,'yes')
   return
end

%     Scale traces relative to one another
switch param.scale

   case {'yes','median'}
      trace_max=zeros(1,ntr);
      for ii=1:ntr
         temp=abs(yi(:,ii));
         trace_max(ii)=median(temp(temp>0  &  ~isnan(temp)));
      end
      trace_max(isnan(trace_max))=1;
      yi=mrt(yi,mean(trace_max)./(trace_max+eps));

   case 'max'
      trace_max=max(abs(yi));
      yi=mrt(yi,mean(trace_max)./(trace_max+eps));
      
   otherwise
      % Do nothing
end

% 	Compute limits for color display