plotimage_2.m

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

function plotimage(smat,t,x,clip,clrmap)
%
% plotimage(smat,t,x,clip,clrmap)
%
% PLOTIMAGE does a quick plot of a seismic matrix in a figure window
% (made by plotimage). It plots the the seismic matrix as a color image
% using the current color map.
%
%	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.
%  clip ... data clip level in standard deviations from the mean
%   ***** default =4 ******* (enter nan for no clipping)
%  clrmap ... the colormap to use
%    ***** defaults to seisclrs(64) ******
% 
%
% NOTE: Scaling options are controlled by defining a global variable in your
% base workspace and setting its value as needed. For example, at the matlab
% prompt, type:
% >> global scaleopt
% This defines the globale variable called "scaleopt". Plotimage recognizes the
%	following values for scaleopt:
%   scaleopt=1 ... The mean and standard deviation of all samples are computed.
%		samples>= mean+clip*stddev are clipped (set equal to mean+clip*stddev).
%		(The same is done for negative samples which are set to mean -clip*stddev).
%		These two extremes are mapped to the ends of the color map and all intermed
%		values are displayed linearly.
%   scaleopt=2 ... The maximum absolute value of the data is computed (say mxs).
%		The extremes of the colormap are then assigned to +/- mxs and all intermed
%		values are displayed linearly. There is no clipping but the display may
%		be dominated by any large values.
% It is suggested that mean scaling (scaleopt=1) be used for real data (this is the
%	default) and max scaling (scaleopt=2) be used for synthetic data. 
%
% G.F. Margrave, CREWES Project, U of Calgary, 1996

%
% Gary F. Margrave
% The CREWES Project
% Department of Geology and Geophysics
% The University of Calgary
% [e-mail]
%

global scaleopt

if ( isempty(scaleopt) ) scaleopt=1; end %default to mean scaling

if( nargin < 1 | ~ischar(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<5)
		clrmap=seisclrs(64);
	end
	if(nargin<4)
		clip=4;
	end
	
	clips=[30 25 20 15 10 9 8 7 6 5 4 3 2 1 .5 .25 .1];
	iclip=near(clips,clip);
	clip=clips(iclip); %make sure we have a sanctioned value

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

	figure('menubar','none');
	%

	colormap(clrmap);

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

	%scale the image
	clrmap=get(gcf,'colormap');
	[nkols,m]=size(clrmap);
	
	mxs=full(max(max(abs(smat))));
	%determine clipping
	smean=full(mean(mean(smat)));
	stddev=full(sqrt( sum(sum((smat-smean).^2 ) )...
				/prod(size(smat))));
							
	disp(['data maximum: ' num2str(full(mxs))])
	disp(['data mean: ' num2str(full(smean))])
	disp(['data stddev: ' num2str(full(stddev))])
	
	if( scaleopt ==1) %mean scaling
		if(~isnan(clip))
			mxsprime=min([smean+clip*stddev,mxs]);
		end
		mns=-mxsprime;
		disp(['mean scaling ']);
		disp(['data clipped outside or mean +/- '...
			num2str(clip) ' standard deviations ']);
		disp(['sigma = ' num2str(full(stddev))])
		seis = (smat -mns)/(mxsprime-mns)*(nkols-1)+1;
		clear smat
	elseif( scaleopt==2 )
		mns=-mxs;
		disp(['max scaling ']);
		seis = (smat -mns)/(mxs-mns)*(nkols-1)+1;
		clear smat
	else
		error('invalid scaling option');
	end

	smat=[];
	
   		ix=1:length(x);
		dx=x(2)-x(1);
		dt=t(2)-t(1);
		if(dx~=(x(3)-x(2)))
			dx=0;
		end
		if(dt~=(t(3)-t(2)))
			dt=0;
		end
		
		hi=image(x,t,seis(:,ix));
		

 %install zooming
 selboxinit('plotimage(''zoom'')',1);

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

 %make a few buttons
 sep=.005;
 ht=.05;
 wd=.25;
 hflip=uicontrol('style','popupmenu',...
		'string','Normal Polarity|Reverse Polarity',...
 	'units','normalized',...
	'position',[sep 0 wd ht],'callback','plotimage(''flip'')',...
	'userdata',1);
 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=.15;
 hcmap=uicontrol('style','pushbutton','string','Colormap',...
 	'units','normalized',...
	'position',[x 0 wd ht],'callback','plotimage(''colormap'')');
 x=x+wd+sep;
 wd=.17;
 hmsg=uicontrol('style','text','string','Polarity Normal',...
 	'units','normalized',...
	'position',[x 0 wd ht],'visible','off');
 %x=x+wd+sep;
 x=x+sep;
 wd=.2;
 hscale=uicontrol('style','popupmenu','string',str2mat('Mean scaling',...
 	'Max scaling'),'units','normalized','position',[x,0,wd,ht],...
 	'callback','plotimage(''rescale'')','value',scaleopt);
 vis='on'; if(scaleopt==2) vis='off'; end
 x=x+wd+sep;
 wd=.15;
 nclips=length(clips);
 clipmsg=num2strmat(clips);
 clipmsg=[ones(nclips,1)*'Cliplevel: ' num2strmat(clips)];  
 hclip=uicontrol('style','popupmenu','string',clipmsg,...
 	'units','normalized','position',[x,0,wd,ht],...
 	'callback','plotimage(''rescale'')','value',iclip,...
 	'visible',vis);

	set(gcf,'userdata',[hflip,hbrighten,hdarken,hmsg,hi,hscale,hclip,hcmap]);
	
	set(hscale,'userdata',[scaleopt mxs mns smean stddev]);
	set(hclip,'userdata',iclip);
	set(hmsg,'userdata',clips);

 %colorview(gca,hi,mns,mxs,0)

 global NAME_
 global NOSIG
 if(~NOSIG)
 	signature(NAME_);
 end

 return;
end

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

 h=get(gcf,'userdata');
 hi=h(5);
 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);
 hmsg=h(4);
 pol=get(hflip,'userdata');
 pol=-1*pol;
	clr=get(gcf,'colormap');
	colormap(flipud(clr));
	set(hflip,'userdata',pol);
	if(pol==1)
		set(hmsg,'string','Polarity Normal');
	elseif(pol==-1)
		set(hmsg,'string','Polarity Reversed');
	end
	return;
end;

if(strcmp(action,'colormap'))
	h=get(gcf,'userdata');
	hmsg=h(4);
	hi=h(5);
	hscale=h(6);
	hclip=h(7);
	
	%get the data
	seis=get(hi,'cdata');
	
	%determine old scaling
	dat=get(hscale,'userdata');
	oldscaleopt=dat(1);
	mxs=dat(2); mns=dat(3); smean=dat(4); stddev=dat(5);
	
	colorview(gca,hi,mns,mxs,0);
	
	return;
end

if(strcmp(action,'rescale'))
	h=get(gcf,'userdata');
	hmsg=h(4);
	hi=h(5);
	hscale=h(6);
	hclip=h(7);
	
	%get the data
	seis=get(hi,'cdata');
	
	%determine old scaling
	dat=get(hscale,'userdata');
	oldscaleopt=dat(1);
	mxs=dat(2); mns=dat(3); smean=dat(4); stddev=dat(5);
	%new opt
	newscaleopt=get(hscale,'value');
	dat(1)=newscaleopt;
	set(hscale,'userdata',dat);
	
	%get clip value
	inewclip=get(hclip,'value');
	ioldclip=get(hclip,'userdata');
	set(hclip,'userdata',inewclip);
	clips=get(hmsg,'userdata');
	clipold=clips(ioldclip);
	clipnew=clips(inewclip);
	
	%get number of columns in colormap
	clrmap=get(gcf,'colormap');
	nkols=size(clrmap,1);
	
	flag=computer;
	flag='shit';
	
	if(strcmp(flag,'MAC2'))
		ntr=size(seis,2);
		for k=1:ntr
			%undo the old scaling
			if( oldscaleopt == 1 ) %undo mean scaling
				if(k==1) mxsprime = min([smean+clipold*stddev,mxs]);
					mns=-mxsprime;
				end
				tmp = (seis(:,k)-1)*(mxsprime-mns)/(nkols-1) + mns;
			elseif( oldscaleopt == 2) %undo max scaling
				if(k==1)	mns=-mxs; end
				tmp = (seis(:,k)-1)*(mxs-mns)/(nkols-1) + mns;
			end
			
			%apply new scaling
			if(newscaleopt==1)
				if(k==1) disp(['mean scaling ']);
					mxsprime = min([smean+clipnew*stddev,mxs]);
					mns=-mxsprime;
					set(hclip,'visible','on');
				end
				seis(:,k) = (tmp -mns)/(mxsprime-mns)*(nkols-1)+1;
			elseif(newscaleopt==2)
				if(k==1) disp(['max scaling ']);
					mns=-mxs;
					set(hclip,'visible','off');
				end
				seis(:,k) = (tmp -mns)/(mxs-mns)*(nkols-1)+1;
				
			end
		end
	else
		%undo the old scaling
		if( oldscaleopt == 1 ) %undo mean scaling
			mxsprime = min([smean+clipold*stddev,mxs]);
			mns=-mxsprime;
			smat = (seis-1)*(mxsprime-mns)/(nkols-1) + mns;
		elseif( oldscaleopt == 2) %undo max scaling
			mns=-mxs;
			smat = (seis-1)*(mxs-mns)/(nkols-1) + mns;
		end
		
		%apply new scaling
		if(newscaleopt==1)
			mxsprime = min([smean+clipnew*stddev,mxs]);
			mns=-mxsprime;
			disp(['mean scaling ']);
			seis = (smat -mns)/(mxsprime-mns)*(nkols-1)+1;
			set(hclip,'visible','on');
		elseif(newscaleopt==2)
			disp(['max scaling ']);
			mns=-mxs;
			seis = (smat -mns)/(mxs-mns)*(nkols-1)+1;
			set(hclip,'visible','off');
		end
	end

	
	set(hi,'cdata',seis);
	return;
end