theodisplay.m

基于matlab的反演程序,用于地球物理勘探中射线追踪及偏移成像程序.

function theodisplay(theo,t,rr,w,tw,phs,hhors,thors,tscale,xscale,...
		title,plotmode,nreps,tint,ievery,kolor,lasheader)

% theodisplay(theo,t,rr,wavelet,tw,phs,hhors,thors,tscale,xscale,title,
%             plotmode,nreps,tint,ievery,kolor,lasheader)
%
% THEODISPLAY is called by LOGSEC to create a display of a theogram showing
% the wavelet and reflectivity (in time) that it was compute from, the
% theogram, and an arbitrary number of phase rotated copies of the theogram.
% Also shown are any named horizons which intersect the
% theogram. THEODISPLAY then controls the figure and provides facilities for
% zooming, grid, more phase rotations, altering the title, and hardcopy.
% Once launched, the figure has no further communication with LOGSEC.
%
%	theo = a column vector containing the theogram
%	t= column vector of times for theo and r
%	rr = a 2 column matrix. First column is a vector of the reflectivity samples
%           while the secone is the times for the reflectivities
%	w = a column vector with the wavelet for theo. 
%			(theo was computed by convolution of r and w)
%	tw = column vector of times for w
%	phs = vector of phase angles (in degrees). A phase rotated copy of theo
%		will be generated for each phs
%	hhors = vector of handles of the horizons which intersect theo. At the 
%		time of the creation of this display, these horizons must be
%		drawn in an active figure.
%	thors = vector of times at which the above horizons intersect theo
%	tscale = default time scale (seconds per inch) for hardcopy
%	xscale = default horizontal scale (traces per inch) for hardcopy
%	title = title
%   plotmode = flag ... 1-> wiggle trace
%                       2-> wiggle trace variable area
%   nreps =  number of times to repeat the plot of each phase rotation
%   **** default=1 ****
%   tint = interval between timing lines. 
%          Must be one of: .5 .25 .1 .05 .025 .01
%          **** default = .1 ****
%  ievery = timing line labeling frequency. The first timing line will be 
%           labeled and every ievery'th line after that. 
%           Must be one of: 1 2 4 5 10
%           **** default =1 ****
%  kolor = rgb triplet defining color of the theogram traces
%          **** default = [0 0 1] (blue) ****
%  lasheader =  string matrix containing an LAS header appropriate for the
%               theogram. If provided then the 'writelas' menu will be
%               active and allow the theogram, impedance, rcs, and wavelet
%               to be written as an LAS log in time.
%               (Log mnemonic matrix in MAKELASHEADER should be: 
%		['THEO';'IMPD';'RCS ';'WLET']
%		NULLS value should be -999.25, and units should be 'S')
%               **** default = '' (no write capability) ****
% 
% NOTE: It is illegal for you to use this software for a purpose other
% than non-profit education or research UNLESS you are employed by a CREWES
% Project sponsor. By using this software, you are agreeing to the terms
% detailed in this software's Matlab source file.
 
% BEGIN TERMS OF USE LICENSE
%
% This SOFTWARE is maintained by the CREWES Project at the Department
% of Geology and Geophysics of the University of Calgary, Calgary,
% Alberta, Canada.  The copyright and ownership is jointly held by 
% its author (identified above) and the CREWES Project.  The CREWES 
% project may be contacted via email at:  [e-mail]
% 
% The term 'SOFTWARE' refers to the Matlab source code, translations to
% any other computer language, or object code
%
% Terms of use of this SOFTWARE
%
% 1) Use of this SOFTWARE by any for-profit commercial organization is
%    expressly forbidden unless said organization is a CREWES Project
%    Sponsor.
%
% 2) A CREWES Project sponsor may use this SOFTWARE under the terms of the 
%    CREWES Project Sponsorship agreement.
%
% 3) A student or employee of a non-profit educational institution may 
%    use this SOFTWARE subject to the following terms and conditions:
%    - this SOFTWARE is for teaching or research purposes only.
%    - this SOFTWARE may be distributed to other students or researchers 
%      provided that these license terms are included.
%    - reselling the SOFTWARE, or including it or any portion of it, in any
%      software that will be resold is expressly forbidden.
%    - transfering the SOFTWARE in any form to a commercial firm or any 
%      other for-profit organization is expressly forbidden.
%
% END TERMS OF USE LICENSE


if( ~isstr(theo) )
	action='init';
else
	action=theo;
end

if( strcmp(action,'init') )
 if(nargin<17)
		lasheader='';
	end
	if(nargin<16)
		kolor=[0 0 1];
	end

	if(nargin < 15)
		ievery=1;
	end
	
	if(nargin<14)
		tint=.1;
	end

	if(nargin<13)
		nreps=1;
	end
	
	it=find( [1 2 4 5 10]==ievery );
	if(isempty(it))
		error('Invalid value for ''ievery'' ');
	end
	
	it=find( [.5 .25 .1 .05 .025 .01]==tint);
	if(isempty(it))
		error('Invalid value for ''tint'' ');
	end

% make a new figure
	hfig=figure('visible','off','menubar','none','numbertitle','off','name','Theogram display');
   scrsize=get(0,'screensize');
	pos=get(hfig,'position');
   swd=scrsize(3);sht=scrsize(4);
	set(hfig,'units','pixels','position',[.1*swd .1*sht .4*swd .7*sht]);

% make sure everything is a column vector
	theo=theo(:);
	t=t(:);
	w=w(:);
	tw=tw(:);
	hhors=hhors(:);
	thors=thors(:);
	
	r=rr(:,1);
	tr=rr(:,2);

% make a menu
	hoptions=uimenu(gcf,'label','Options');
	
	% turn grid on and off... default to on
	hgrid=uimenu(hoptions,'label','Timing Lines');
	hgridon=uimenu(hgrid,'label','Show Timing Lines','checked','on','callback',...
			'theodisplay(''grid'')');
	htimeint=uimenu(hgrid,'label','Line Interval');
		if(tint==.5) chk='on'; mult=-1; else chk='off'; mult=1; end
		htime500=uimenu(htimeint,'label','500 ms','callback','theodisplay(''timeint'')',...
			'userdata',.5,'checked',chk);
		htime500=htime500*mult;

		if(tint==.25) chk='on'; mult=-1; else chk='off'; mult=1; end
		htime250=uimenu(htimeint,'label','250 ms','callback','theodisplay(''timeint'')',...
			'userdata',.250,'checked',chk);
		htime250=htime250*mult;

		if(tint==.1) chk='on'; mult=-1; else chk='off'; mult=1; end
		htime100=uimenu(htimeint,'label','100 ms','callback','theodisplay(''timeint'')',...
			'userdata',.1,'checked',chk);
		htime100=htime100*mult;

		if(tint==.05) chk='on'; mult=-1; else chk='off'; mult=1; end
		htime50=uimenu(htimeint,'label','50 ms','callback','theodisplay(''timeint'')',...
			'userdata',.050,'checked',chk);
		htime50=htime50*mult;

		if(tint==.025) chk='on'; mult=-1; else chk='off'; mult=1; end
		htime25=uimenu(htimeint,'label','25 ms','callback','theodisplay(''timeint'')',...
			'userdata',.025,'checked',chk);
		htime25=htime25*mult;

		if(tint==.01) chk='on';  mult=-1;else chk='off'; mult=1; end
		htime10=uimenu(htimeint,'label','10 ms','callback','theodisplay(''timeint'')',...
			'userdata',.010,'checked',chk);
		htime10=htime10*mult;
	set(htimeint,'userdata',[htime500 htime250 htime100 htime50 htime25 htime10]);

	htimelabel=uimenu(hgrid,'label','Line Labels');
		if(ievery==1) chk='on';mult=-1; else chk='off';mult=1; end
		hlabel1=uimenu(htimelabel,'label','Every Line','callback',...
			'theodisplay(''label'')','userdata',1,'checked',chk);
		hlabel1=hlabel1*mult;
		if(ievery==2) chk='on';mult=-1; else chk='off';mult=1; end
		hlabel2=uimenu(htimelabel,'label','Every Other Line','callback',...
			'theodisplay(''label'')','userdata',2,'checked',chk);
		hlabel2=hlabel2*mult;
		if(ievery==4) chk='on';mult=-1; else chk='off';mult=1; end
		hlabel4=uimenu(htimelabel,'label','Every Fourth Line','callback',...
			'theodisplay(''label'')','userdata',4,'checked',chk);
		hlabel4=hlabel4*mult;
		if(ievery==5) chk='on';mult=-1; else chk='off';mult=1; end
		hlabel5=uimenu(htimelabel,'label','Every Fifth Line','callback',...
			'theodisplay(''label'')','userdata',5,'checked',chk);
		hlabel5=hlabel5*mult;
		if(ievery==10) chk='on';mult=-1; else chk='off';mult=1; end
		hlabel10=uimenu(htimelabel,'label','Every Tenth Line','callback',...
			'theodisplay(''label'')','userdata',10,'checked',chk);
		hlabel10=hlabel10*mult;
	set(htimelabel,'userdata',[hlabel1 hlabel2 hlabel4 hlabel5 hlabel10]);
	
	% a new phase rotation
	hphase=uimenu(hoptions,'label','New Phase','callback',...
		'theodisplay(''newphs'')');

	% a locate action
	hlocate=uimenu(hoptions,'label','Locate','callback','theodisplay(''locate'')');
		
	% hardcopy
	hhard=uimenu(hoptions,'label','Hardcopy','callback',...
		'theodisplay(''hardcopy'')');

	%writelas
	vis='on';
	if(strcmp(lasheader,'')) vis='off'; end
	hwritelas=uimenu(hoptions,'label','Write LAS','callback',...
		'theodisplay(''writelas'')','enable',vis);
	% quit
	hquit=uimenu(hoptions,'label','Close','callback','theodisplay(''quit'')');
	
	% a message panel
	
	hmsg=uicontrol('style','text','string','MB3 drag -> Zoom ... MB3 click -> unZoom',...
		'units','normalized','position',[0 0 1 .03]);
		
	% compute the phase rotated copies of theo
	theophs=zeros(length(theo),length(phs));
	for k=1:length(phs)
		theophs(:,k)=phsrot(theo,phs(k));
	end
	
	%make a storage bucket
	hstore1=uicontrol('style','text','visible','off');
		
	% store things
	
	set(gcf,'userdata',[hmsg,hoptions,hgrid,hgridon,htimeint,hphase,hhard,hquit,...
		htimelabel,hlocate,hwritelas,hstore1]);

	%integrate the reflectivity to get impedance
	imp=integrate(r);

	% make the impedance zero mean
	imp=imp-mean(imp);
	
	%userdata assignments:
	% h=get(gcf,'userdata') =>[hmsg,hoptions,hgrid,hgridon,htimeint,hphase,hhard,hquit,
	%				htimelabel,hlocate,hwritelas]
	% hmsg = h(1) ... vectors of theogram, time, impedance, rcs:  [theo t imp r]
	% hoptions = h(2) ... [w tw] wavelet and wavelet times
	% hgrid = h(3) ... [hhors thors] vector of horizon handles in the master figure,
	%				and vector of horizon times
	% hgridon = h(4) ... [plotmode nreps abs(title)] plotmode flag, number of phase reps
	%				and title
	% htimeint = h(5) ... vector of submenu handles
	% hphase = h(6) ... vector of theogram phases
	% hhard = h(7) ... x and time scale for hardcopy
	% hquit = h(8) ... matrix of phase rotated theograms
	% htimelabel = h(9) ... vector of submenu handles
	% hlocate = h(10) ... color to plot the theograms
	% hwritelas = h(11) ... string matrix with las header
	% hstore1 = h(12) ... handles of the lines and labels of the tops:[ntops hlines hlbls]
	%
	% 
	hoot{1}=theo;
	hoot{2}=t;
	hoot{3}=imp;
	hoot{4}=r;
	hoot{5}=tr;
	set(hmsg,'userdata',hoot);
	set(hoptions,'userdata',[w tw]);
	set(hgrid,'userdata',[hhors thors]);
	set(hgridon,'userdata',[plotmode nreps abs(title)]);
	set(hphase,'userdata',phs);
	set(hhard,'userdata',[xscale tscale]);
	set(hquit,'userdata',theophs);
	set(hlocate,'userdata',kolor);
	set(hwritelas,'userdata',lasheader);
	
	%plot things
	
	theodisplay('plot');
	set(gca,'box','off');
	set(gca,'position',[.13 .11 .7 .8]);
	
	%turn on zooming with button 3
	simplezoom(3,'theodisplay(''repositops'')');

	set(hfig,'visible','on');
	
	return;
	
end

if(strcmp(action,'plot'))
	h=get(gcf,'userdata');
	hmsg=h(1);
	hoptions=h(2);
	hgrid=h(3);
	hgridon=h(4);
	htimeint=h(5);
	hphase=h(6);
	hhard=h(7);
	hquit=h(8);
	htimelabel=h(9);
	hlocate=h(10);
	hstore1=h(12);

	dat=get(hgridon,'userdata');
	plotmode=dat(1);
	nreps=dat(2);
	title=setstr(dat(3:length(dat)));
	
	labelangle=45;
	
	% get the trace etc
	trcs=get(hmsg,'userdata');
	theo=trcs{1};
	t=trcs{2};
	imp=trcs{3};
	rc=trcs{4};
	tr=trcs{5};
	
	trcs=get(hoptions,'userdata');
	w=trcs(:,1);
	tw=trcs(:,2);
	
	trcs=get(hquit,'userdata');
	
	horstuff=get(hgrid,'userdata');
	if(~isempty(horstuff))
		hors=horstuff(:,1);
		thors=horstuff(:,2);
	else
		hors=[];
		thors=[];
	end
	
	phs=get(hphase,'userdata');

	%clear the axes
	cla;
	
	% plot the wavelet
	%first truncate the wavelet if it is longer than the theogram.
	if(length(w)>length(theo))
		ind=near(tw,0);
		i1=round(max([1 ind-length(theo)/2]));
		w=w(i1:i1+length(theo)-1);
		tw=tw(i1:i1+length(theo)-1);
    else
        w=pad(w,theo,1);
        tw=t;
	end
	x=1;
	wm=max(abs(w));
	%make an appropriate time axis for w
	tw2=xcoord(t(1),tw(2)-tw(1),w);
	if(plotmode==1)
		plot(x+w/wm,tw2,'r');
	else
		wtva(x+w/wm,tw2,'r',x,1,1);
	end

	set(gca,'ydir','reverse');
	ytext=t(1)-.01*(t(length(t))-t(1));
	text('string','Wavelet','units','data','rotation',labelangle,'position',...
			[x ytext]);
	
	%plot the impedance
	bigger=2;
	x=x+2;
	im=max(abs(imp));
	line(x+bigger*imp/im,tr,'color','r');
	text('string','Impedance','units','data','rotation',labelangle,'position',...
			[x ytext]);
	
	
	%plot the rcs
	x=x+2;
	ilive=find(~isnan(rc));
	rm=max(abs(rc(ilive)));
	%if( plotmode==1 )
		line(x+bigger*rc/rm,tr,'color','k');
	%else
		%wtva(x+rc(ilive)/rm,tr(ilive),'r',x,1,1);
	%end
	text('string','RCs','units','data','rotation',labelangle,'position',...
			[x ytext]);

	kolor=get(hlocate,'userdata');
	
	%plot the theogram
	thmax=max(abs(theo));
	x=x+1;
   ilabel=floor((nreps+1)/2);
	for k=1:nreps
		x=x+1;
		if( plotmode== 1)
			line(x+theo/thmax,t,'color',kolor);
		else
			wtva(x+theo/thmax,t,kolor,x,1,1);
		end
      if(k==ilabel)
		text('string','Theogram','units','data','rotation',labelangle,'position',...
				[x ytext]);
      end
	end
	
	%plot the phase rotations
	for k=1:length(phs)
		x=x+1;
		for kk=1:nreps
			x=x+1;
			if( plotmode== 1)
				line(x+trcs(:,k)/thmax,t,'color',kolor);
			else
				wtva(x+trcs(:,k)/thmax,t,kolor,x,1,1);
			end
         if(kk==ilabel)
			text('string',[num2str(phs(k)) ' Phs'],'units','data','rotation',labelangle,...
				'position',[x ytext]);
         end
		end
	end
	
	%plot the horizons
	nhors=length(hors);
	hlines=nan*zeros(1,nhors);
	hlbls=nan*zeros(1,nhors);
	xlim=get(gca,'xlim');
	x=max(xlim);
	for k=1:nhors
		if(~isnan(thors(k)))
			%get the color
			%kol=get(hors(k),'color');
			kol='k';
			%get the name
			dat=get(hors(k),'userdata');
			it=findstr(dat,':');
			it2=findstr(dat,'depth');
			if(isempty(it2))
				it2=length(dat);
			else
				it2=it2-2;
			end
			name=dat(it+2:it2);
			if( ~strcmp(name(1:2),'__') )
			
				%plot
				hlines(k)=line(xlim,[thors(k) thors(k)],'color',kol,'linestyle','-.');
				hlbls(k)=text('string',name,'units','data',...
					'position',[x thors(k)],'color',kol,'fontsize',9);
			end
		end
	end
	ind=find(isnan(hlines));
	if(~isempty(ind))
		hlines(ind)=[];
		hlbls(ind)=[];
		nhors=nhors-length(ind);
	end
	set(hstore1,'userdata',[nhors hlines hlbls]);

	%surpress x ticks
	set(gca,'xtick',[]);

	%determine the yticks
	%determine the tick interval
	dat=get(htimeint,'userdata');