seisimage.m

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

function seisimage(smat,t,x,clrmap,clip)

% seisimage(smat,t,x,clrmap,clip)
%
% SEISIMAGE does a quick plot of a seismic matrix in a figure window
% (made by seisimage). It plots the the seismic matrix as a color image
% using the colormap generated by seisclrs.m or the one provided as a fourth
% argument. It also provides interactive spectral analysis capabilities
% using fxtran.  FXTRAN amplitude psectra, phase spectra, and average
% amplitude spectra may be computed and displayed over any window. The
% most resent average amplitude spectrum may be accessed via the global
% variables FAVE, AAVE, AMAX which are: the frequencies, the db amplitudes,
% and the db reference amplitude.
%
%	smat ... the seismic matrix to be plotted. Traces are assumed stored in
%       the columns smat.
%	t ... time coordinates of traces.
%       ****** default 0:.002:(nrows-1)*.002 where nrows = number of rows
%       in smat ****
%	x ... x coordinates of the traces
%       ****** default 1:ncols where ncols=number of columns *****
%  NOTE: if only a simgle argument is provided, it is assumed to be a fleximat
%  from which seis,t,and x can be extracted.
%   clrmap ... the colormap to be used for seismic display
%   clip ... data clipping level expressed as a number of standard
%       deviations from the mean
%   ********* default = 4 ******
%   Specify nan for no clipping
%
% G.F. Margrave, Aug 1995
%
% 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( nargin < 1 | ~isstr(smat) )
	action='init';
else
	action = smat;
end

if(strcmp(action,'init'))
	if(nargin<1)
		% do a demo
			%Make a fake reflectivity
			t=0:.002:1.0;
			r=randn(size(t)).^5;
			%make a ricker wavelet
			tw=-.1:.002:.1;
			arg=(pi*15*tw).^2;
			w=(1-2.*arg).*exp(-arg);
			%convolve
			s=conv(r,w);
			s=s(51:length(t)+50)';
			s=s/max(s); %normalize
		
			smat=s*ones(1,20);
	end
	if(nargin<3)
		ncols=size(smat,2);
		x=1:ncols;
	end
	if(nargin<2)
		t=fmget(smat,'y');
		x=fmget(smat,'x');
		smat=fmget(smat,'mat');
	end
	if(nargin<4)
		clrmap=seisclrs(64);
	end
	if(nargin<5)
		clip=1;
	end

	if(length(x)>1)
		bnds=(max(x)-min(x))/(length(x)+1);
	else
		bnds=max(smat-min(smat))/2;
	end

	figure;
	p=get(gcf,'position');
	%set(gcf,'position',[p(1:2) 880 666]);
	%
	colormap(clrmap)

	%set(gca,'ydir','reverse');

	%scale the image
	clrmap=get(gcf,'colormap');
	[nkols,m]=size(clrmap);
	ilive=find(~isnan(smat));
	mxs=max(abs(smat(ilive)));

	%determine clipping
	if(~isnan(clip))
		smean=mean(smat(ilive(:)));
		stddev=sqrt( sum((smat(ilive(:))-smean).^2 )/length(smat(ilive(:))));
		mxs=min([smean+clip*stddev,mxs]);
	end
	mns=-mxs;

	seis = (smat -mns)/(mxs-mns)*(nkols-1)+1;
	% note the above is really:
	% seis= (smat+mxs)/(2*mxs)*(nkols-1)+1
	% so
	% smat = 2*mxs*(seis-1)/(nkols-1) - mxs
	%
	smat=[];

 [x,ix]=sort(x);

	hi=image(x,t,seis(:,ix));
	set(gcf,'userdata',hi);
 %install zooming
 selboxinit('seisimage(''zoom'')',1);
 whitefig;
 %set(gca,'xdir','reverse');
 %set(gca,'ylim',[0 2.6])
 %xlabel('Shotpoint')
 %ylabel('Seconds')


 set(gcf,'name','Seismic Image Plot, Simplezooming installed (Use MB1)')

 %make a few buttons
 sep=.005;
 ht=.05;
 wd=.15;
 hflip=uicontrol('style','pushbutton','string','Flip Polarity',...
 'units','normalized',...
	'position',[0 0 wd ht],'callback','seisimage(''flip'')',...
	'userdata',[1 hi mxs nkols]);
	 x=wd+sep;
	 wd=.1;
 hbrighten=uicontrol('style','pushbutton','string','Brighten',...
 'units','normalized',...
	'position',[x 0 wd ht],'callback','brighten(.5)');
	 x=x+wd+sep;
 hdarken=uicontrol('style','pushbutton','string','Darken',...
 'units','normalized',...
	'position',[x 0 wd ht],'callback','brighten(-.5)');
	 x=x+wd+sep;
	 wd=.2;
 hpolmsg=uicontrol('style','text','string','Polarity Normal',...
 'units','normalized',...
	'position',[x 0 wd ht]);
	 x=x+wd+sep;
	 wd=.4;
 hmsg=uicontrol('style','text','string','',...
 'units','normalized',...
	'position',[x 0 wd ht]);

	%build a menu bar
	htools=uimenu(gcf,'label','Tools');
	hfxtran=uimenu(htools,'label','fxtran',...
		'callback','seisimage(''fxtran'')');

 fnyq=1/(2*(t(2)-t(1)));
	hfxopt=uimenu(htools,'label','Set fxtran parameters',...
		'callback','seisimage(''fxopt'')','userdata',[1 20 20 fnyq 1 1]);
	hzone=uimenu(htools,'label','Set tool space/time zone',...
		'callback','seisimage(''zone'')');

	set(gcf,'userdata',[hflip,hbrighten,hdarken,hpolmsg,hfxtran,hfxopt,...
		hzone,hmsg,htools])

	%
	% userdata
	% hflip ... ipol hi mxs nkols
	% hbrighten ... handle of average spectra window
	% hdarken ... not used
	% hplomsg ... not used
	% hfxtran ... not used
	% hfxopt ... fx parms [secflag pctaper trand fmax dbplot] 
	%  sectflag = 1 if both
	% 	amp and phase, 1 for phase only, 2 for amp only
	% hzone ... if null, then the current zone is taken from the axis 
	%	limits. Otherwise, its: [xmin xmax tmin1 tmax1 tmin2 tmax2]
	% hmsg ... not used
	% htools ... not used

 %colorview(gca,hi,mns,mxs,0)
 
 return;
end

if(strcmp(action,'zoom'))
 box=selboxfini;

 h=get(gcf,'userdata');
 hflip=h(1);
 dat=get(hflip,'userdata');
 hi=dat(2);
 if(length(box)==5)
		delete(box(5));
	end

	xmin=min([box(1) box(3)]);
	xmax=max([box(1) box(3)]);
	ymin=min([box(2) box(4)]);
	ymax=max([box(2) box(4)]);
	%get the current axis settings
	xlim=get(gca,'xlim');
	ylim=get(gca,'ylim');
	test1=xmin-xlim(1)+xmax-xlim(2)+ymin-ylim(1)+ymax-ylim(2);
	test2=(xmin-xmax)*(ymin-ymax);
	if(abs(test1)<10*eps | abs(test2)< 10*eps)
		x=get(hi,'xdata');
		y=get(hi,'ydata');
		set(gca,'xlim',[min(x) max(x)],'ylim',[min(y) max(y)]);
	else
		set(gca,'xlim',[xmin xmax],'ylim',[ymin ymax]);
	end
	return;
end

if(strcmp(action,'flip'))
 h=get(gcf,'userdata');
 hflip=h(1);
 hpolmsg=h(4);
 dat=get(hflip,'userdata');
 pol=dat(1);
 pol=-1*pol;
	clr=get(gcf,'colormap');
	colormap(flipud(clr));
	set(hflip,'userdata',[pol dat(2:length(dat))]);
	if(pol==1)
		set(hpolmsg,'string','Polarity Normal');
	elseif(pol==-1)
		set(hpolmsg,'string','Polarity Reversed');
	end
	return;
end;

if(strcmp(action,'fxopt'))
	h=get(gcf,'userdata');
	% parameters to define
	% sections desired, pctaper, trand
	% time zone comes from tool zone
	hfxopt=h(6);
	hmsg=h(8);
	fxparms=get(hfxopt,'userdata');
	q=str2mat('FX sections desired:','Percentage taper:',...
		'Random fluct (pct):','Max frequency (Hz):',...
		'FX Amplitude Plot:','Window Type');
	a=str2mat('Amp and Phase|Phase Only|Amp Only|Average Amp Only',...
		int2str(fxparms(2)),...
		num2str(fxparms(3)),num2str(fxparms(4)),'Linear|Decibels',...
		'Mwindow|Hanning|Hamming|Bartlett|Gaussian');
	flags=ones(1,6);
	flags(1)=fxparms(1);
	flags(5)=fxparms(5);
	flags(6)=fxparms(6);

	set(hmsg,'string','Please wait for dialog to appear');

	askthingsinit('seisimage(''fxopt2'')',q,a,flags);
	 
	return;
end

if(strcmp(action,'fxopt2'))
	h=get(gcf,'userdata');
	hflip=h(1);
	hfxopt=h(6);
	hmsg=h(8);

	reply=askthingsfini;
	%check for cancel
	if(reply==-1)
		set(hmsg,'string','FX parms cancelled');
		return;
	end
	if(strcmp(reply(1,1:5),'Amp a'))
		fxflag=1;
	elseif(strcmp(reply(1,1:5),'Amp O'))
		fxflag=3;
	elseif(strcmp(reply(1,1:5),'Phase'))
		fxflag=2;
	else
		fxflag=4;
	end
	if(strcmp(reply(5,1:2),'Li'))
		dbflag=1;
	else
		dbflag=2;
	end
	if( strcmp(reply(6,1:5),'Mwind') )
		wintype=1;
	elseif( strcmp(reply(6,1:5),'Hanni'))
		wintype=2;
	elseif( strcmp(reply(6,1:5),'Hammi'))
		wintype=3;
	elseif( strcmp(reply(6,1:5),'Bartl'))
		wintype=4;
	elseif( strcmp(reply(6,1:5),'Gauss'))
		wintype=5;
	end

	pct=round(str2num(reply(2,:)));
	pctrand=str2num(reply(3,:));
	fmax=str2num(reply(4,:));

	dat=get(hflip,'userdata');
	hi=dat(2);
	t=get(hi,'ydata');
	fnyq=1/(2*(t(2)-t(1)));

	fxparms=[1 20 20 fnyq 1 1];

	if(pct<0 | pct>100 | pctrand <0 | pctrand > 100 | fmax<0 | ...
		fmax>fnyq+1000000*eps )
		set(hmsg,'string','Invalid parms, defaults restored');
	else
		fxparms(1)=fxflag;
		fxparms(2)=pct;
		fxparms(3)=pctrand;
		fxparms(4)=fmax;
		fxparms(5)=dbflag;
		fxparms(6)=wintype;
		set(hmsg,'string','Fx parms changed');
	end

	set(hfxopt,'userdata',fxparms);

	return;
end

if(strcmp(action,'fxtran'))
 h=get(gcf,'userdata');
 hflip = h(1);
 hbrighten=h(2);
 hfxopt=h(6);
 hzone = h(7);
 hmsg=h(8);
 
 global FAVE AAVE AMAX;

 set(hmsg,'string','Computing ...');
 set(gcf,'pointer','watch');
	%get some plot params from the mother figure
	xdir=get(gca,'xdir');
	pos=get(gcf,'position');
 hmasterfig=gcf;

 %get the data
 dat = get(hflip,'userdata');
 hi=dat(2);
 mxs=dat(3);
 nkols=dat(4);

 % smat = 2*mxs*(seis-1)/(nkols-1) - mxs
 smat = 2*mxs*(get(hi,'cdata')-1)/(nkols-1)-mxs;

 % get x and t
 x = get(hi,'xdata');
 t=get(hi,'ydata');
 dt=t(2)-t(1);

 %get the zone
 zone = get(hzone,'userdata');
 if(isempty(zone)) 
		xlims= sort(get(gca,'xlim'));
		tlims=sort(get(gca,'ylim'));
		if(tlims(1)<t(1)) tlims(1)=t(1); end
		if(tlims(2)>t(length(t))) tlims(2)=t(length(t)); end
		if(xlims(1)<min(x)) xlims(1)=min(x); end
		if(xlims(2)>max(x)) xlims(2)=max(x); end
		tmin=[tlims(1) tlims(1)];
		tmax=[tlims(2) tlims(2)];
	else
		xlims = zone(1:2);
		tmin = [zone(3) zone(5)];
		tmax = [zone(4) zone(6)];
	end

	tit=get(get(gca,'title'),'string');

	%get the fxparms
	fxparms=get(hfxopt,'userdata');
	isec=fxparms(1);
	pctaper=fxparms(2);
	pctrand=fxparms(3);
	fmax=fxparms(4);
	dbflag=fxparms(5);
	wintype=fxparms(6);

	% determine the colums to transform
	icols = near(x,xlims(1),xlims(2));

	fmseis = fmset([], x(icols), t, smat(:,icols) );
	[fmphs,fmamp]=fxtran(fmseis,tmin,tmax,pctaper,pctrand,...
		.5/dt,0,min(x(icols)),max(x(icols)),wintype);

	%compute average amplitude spectrum
	f=fmget(fmamp,'y');
	ampave=sum(fmget(fmamp,'mat').').';%average
	ampave=ampave/(size(fmamp,2));%normalize for number of traces
	ampave=ampave/(size(fmamp,1));%normalize for number of samples
	%get the display window handle
	ampdat=get(hbrighten,'userdata');
	if(~isempty(ampdat))
		ind=find(figs==ampdat(1));
		if(~isempty(ind))
			hampfig=ampdat(1);
			ampmax=ampdat(2);
			nline=ampdat(3)+1;
		else
			ampdat=[];
		end
	end
	if(isempty(ampdat))
		hampfig=figure;
		%simplezoom;
		ampmax=max(ampave);
		nline=1;
	end
	%%%%%% Hardwired for individual normalizations %%%%%
	ampmax=max(ampave);
	%%%%%%%%
	figure(hampfig);
	hax=get(hampfig,'currentaxes');
	kols=get(hax,'colororder');
	ampdb=todb(ampave,ampmax);
	nkol=rem(nline,size(kols,1));
	if(nkol==0) nkol=size(kols,1); end
	hline=line(f,ampdb,'color',kols(nkol,:));
	ind=find(ampdb==max(ampdb));
	xtxt=f(ind);ytxt=ampdb(ind);
	if( abs(diff(tmin)) | abs(diff(tmax)) )
		text('position',[xtxt ytxt],'color',kols(nkol,:),'string',...
			[' tlims: ' num2str(tmin(1)) ' ' num2str(tmax(1)) ...
			' ' num2str(tmin(2)) ' ' num2str(tmax(2)) ...
			' xlims: ' num2str(min(xlims)) ' ' num2str(max(xlims)) ]);
	else
		text('position',[xtxt ytxt],'color',kols(nkol,:),'string',...
			[' tlims: ' num2str(tmin(1)) ' ' num2str(tmax(1)) ...
			' xlims: ' num2str(min(xlims)) ' ' num2str(max(xlims)) ]);
	end
		
	%assign to globals
	FAVE=f;
	AAVE=ampdb;
	AMAX=ampmax;
	
	set(hbrighten,'userdata',[hampfig ampmax nline]);


	titstr = ['Time Window: ' num2str(tmin(1)) ' to ' ...
		num2str(tmax(1)) ' (sec)'];

	if( isec == 3)
		fmamp=fmfmax(fmamp,fmax);
		if(dbflag==2)
			%convert to db
			x=fmget(fmamp,'x');
			f=fmget(fmamp,'y');
			amp=fmget(fmamp,'mat');
			fmamp=[];
			amax=max(max(amp));
			ind=find(amp~=0.0);
			amp(ind)=20*log(amp(ind)./amax);
			ind=find(amp==0.0);
			if(~isempty(ind))
				amp(ind)=nan*ones(size(ind));
			end
			plotimage(amp,f,x,gray,10);
		else
			plotimage(fmamp);
		end
		%p=get(gcf,'position');
		%set(gcf,'position',[p(1:2) 880 666])
		if(isempty(tit))
			title(['Amplitude Spectrum ' titstr]);
		else
			title(tit);
			xlabel(['Amplitude Spectrum ' titstr]);
		end
		ylabel('Frequency (Hz.)')
		set(gca,'xdir',xdir);
		pa=get(gcf,'position');
		pa(3:4)=pos(3:4);
		set(gcf,'position',pa);
		whitefig;
	elseif( isec ==2 )
		fmphs=fmfmax(fmphs,fmax);
		plotimage(fmphs);
		%p=get(gcf,'position');
		%set(gcf,'position',[p(1:2) 880 666])
		if(isempty(tit))
			title(['Phase Spectrum ' titstr]);
		else
			title(tit);
			xlabel(['Phase Spectrum ' titstr]);
		end
		ylabel('Frequency (Hz.)')
		set(gca,'xdir',xdir);
		pa=get(gcf,'position');
		pa(3:4)=pos(3:4);
		set(gcf,'position',pa);
		whitefig;
	elseif(isec == 1)
		fmphs=fmfmax(fmphs,fmax);
		plotimage(fmphs);
		%p=get(gcf,'position');
		%set(gcf,'position',[p(1:2) 880 666])
		if(isempty(tit))
			title(['Phase Spectrum ' titstr]);
		else
			title(tit);
			xlabel(['Phase Spectrum ' titstr]);
		end
		ylabel('Frequency (Hz.)')
		set(gca,'xdir',xdir);
		pa=get(gcf,'position');
		pa(3:4)=pos(3:4);
		set(gcf,'position',pa);
		whitefig;
		fmamp=fmfmax(fmamp,fmax);
		if(dbflag==2)
			%convert to db
			x=fmget(fmamp,'x');
			f=fmget(fmamp,'y');
			amp=fmget(fmamp,'mat');
			fmamp=[];
			amax=max(max(amp));
			ind=find(amp~=0.0);
			amp(ind)=20*log(amp(ind)./amax);
			dbdmax=-150;
			ind=find(amp<dbdmax);
			if(~isempty(ind))
				amp(ind)=dbdmax*ones(size(ind));
				disp('max dbdown set to dbdmax')
			end
			ind=find(amp==0.0);
			if(~isempty(ind))
				amp(ind)=nan*ones(size(ind));
			end
			plotimage(amp,f,x,gray);
		else
			plotimage(fmamp);
		end
		%p=get(gcf,'position');
		%set(gcf,'position',[p(1:2) 880 666])
		if(isempty(tit))
			title(['Amplitude Spectrum ' titstr]);
		else
			title(tit);
			xlabel(['Amplitude Spectrum ' titstr]);
		end
		ylabel('Frequency (Hz.)')
		pa=get(gcf,'position');
		pa(3:4)=pos(3:4);
		set(gcf,'position',pa);
		whitefig;
		set(gca,'xdir',xdir);
	end

	set(hmsg,'string','FX spectra finished');
	set(hmasterfig,'pointer','arrow');

	return;
end

if(strcmp(action,'zone'))
 h=get(gcf,'userdata');
 hflip = h(1);
 hfxopt=h(6);
 hzone = h(7);
 hmsg=h(8);

 zone=get(hzone,'userdata');
 if(isempty(zone))
		zonestr = 'auto';
	else
		zonestr=[ num2str(zone(1)) ' ' num2str(zone(2)) ' ' ...
			num2str(zone(3)) ' ' num2str(zone(4))...
			' ' num2str(zone(5)) ' ' num2str(zone(6))];
 end
	
	q=' specify x1 x2 tmin1 tmax1 (tmin2 tmax2) :';
	a=zonestr;

	set(hmsg,'string','Please wait for dialog');

	askthingsinit('seisimage(''zone2'')',q,a,1,...
		'Define space and time window');

	return;

end

if(strcmp(action,'zone2'))
 h=get(gcf,'userdata');
 hflip = h(1);
 hfxopt=h(6);
 hzone = h(7);
 hmsg=h(8);

 a=askthingsfini;
 	%check for cancel
	if(a==-1)
		set(hmsg,'string','Space/time zone set cancelled');
		return;
	end

 if(strcmp(strunpad(a),'auto'))
		zone=[];
	else
		z=sscanf(a,'%g');
		if(length(z)<4)
			zone=[];
		else
			x1=min(z(1:2));
			x2=max(z(1:2));
			tmin1=min(z(3:4));
			tmax1=max(z(3:4));
			if(length(z)>4)
				tmin2=min(z(5:6));
				tmax2=max(z(5:6));
			else
				tmin2=tmin1;
				tmax2=tmax1;
			end
			zone=[x1 x2 tmin1 tmax1 tmin2 tmax2];
		end
	end

	set(hzone,'userdata',zone);
	if(isempty(zone))
		set(hmsg,'string','Auto zone on');
	else
		set(hmsg,'string','Manual zone on');
	end

	return;
end