picksle.m

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

    dx=length(xnum);
elseif((xend-xstr)<0)
    % click going from lower to higher
    dx=x1(1)-x2(end);
else
    % click going from higher to lower
    dx=x2(end)-x1(1);
end
if((yend-ystr)<0)
    dy=y1(1)-y2(end);
else
    dy=y2(end)-y1(1);
end
rto=dy/dx;
% tolerance, which is a percent of total length of depth (time or
% otherwise)
pertol=round(.05*size(seisnew,1)/2);

picksdat=[];
ys=[wvltmid];
xdist=xstr-xend;
ydist=ystr-yend;
y2=wvltmid;
dy=t(2)-t(1);
if(xend==xstr)
    allcrosscor={[ ]};
else
    allcrosscor={[100 wvltmid xfirst]};
end
yhold=[wvltmid];
xhold=[xfirst];
perhold=[100];
holdcrosscor=[];
hcc=1;
for ii=2:size(seisnew,2)
    % ii will have to be flipped if user has made line go from smallest to
    % largest
    if(swt==2)
        trc=seisnew(:,size(seisnew,2)-ii+1);
        dx=(size(seisnew,2)-ii+1)-(size(seisnew,2)-ii+2);
        jj=(size(seisnew,2)-ii+1);
        adon=1;
    else
        trc=seisnew(:,ii);
        dx=(ii)-(ii-1);
        jj=ii;
        adon=-1;
    end
    % geting point along initial line
    yup=dx*rto;
    y2=round(y2+yup);
    if(y2<=0)
        y2=1;
    elseif(y2>=length(t))
        y2=size(t,1);
    end
    % ajusting the trace for tolerance and angle 
    wndwtop=y2-round(pertol);
    wndwbot=y2+round(pertol);
    if(wndwtop<=0)
        wndwtop=1;
    end
    if(wndwbot>length(trc))
        wndwbot=length(trc);
    end
    wndw=[wndwtop:1:wndwbot];
    wn=size(wndw);
%     wndwln=line([xnew(jj) xnew(jj)],[t(wndw(1)) t(wndw(end))],'color','g',...
%         'erasemode','xor');
    % creating scansection, it is larger then the wndw because the wavelet
    % will scan through the window, top to bottom, from its center point.
    scansectop=wndwtop-(wvltmid-wvltwndw(1));
    scansecbot=wndwbot+(wvltwndw(end)-wvltmid);
    if(scansectop<=0)
        scansectop=1;
    end
    if(scansecbot>length(trc))
        scansecbot=length(trc);
    end
    scansec=[scansectop:1:scansecbot];
    if(~isempty(showprogress))
        % this will show user progress of picking
        scansecln=line([xnew(jj) xnew(jj)],[t(scansec(1)) t(scansec(end))],'color','g',...
            'erasemode','xor');
        lns=[lns scansecln];
    end
    % scaning 
    scantrc=trc(scansec);
    crosscor=[0 0 0];   % [%crosscorrelation position future]
    if(length(scansec)<=length(wvlt))
        % possibly have to shrink the wvlt size followed by scanning only
        % crosscorrelation only once.
        dsize=length(wvlt)-length(scansec);
        if(round(dsize/2)==dsize/2)
            % even number 
            wvlt2=wvlt((1+dsize/2):1:(length(wvlt)-dsize/2));
        else
            % odd number
            if((-.5+dsize/2)<=0)
                strp=1;
                adp=0;
            else
                strp=(-.5+dsize/2);
                adp=1;
            end
            wvlt2=wvlt(strp:1:(length(wvlt)-round(dsize/2)-adp));
        end
        % setting warning off, because sometimes zeros really screw with
        % junk and I don't feel like worrying about zeros right now
        warning off
        crossc=(sum(wvlt2.*scantrc))/((sum(wvlt2.^2)*sum(scantrc.^2))^0.5);
        warning on
        crosscor=[crossc scansec(round(length(scansec)/2)) xnum(jj)];
    else
        trcwndw=trc(scansec);
        for kk=1:(length(scansec)-length(wvlt))
            warning off
            % seeing if normalization helps
            mat1=wvlt/max(wvlt);
            mat2=scantrc(kk:length(wvlt)+kk-1)/max(scantrc(kk:length(wvlt)+kk-1));
            % setting warning off, because sometimes zeros really screw with
            % junk and I don't feel like worrying about zeros right now
            % crossc=(sum(wvlt.*scantrc(kk:length(wvlt)+kk-1)))/((sum(wvlt.^2)*sum(scantrc(kk:length(wvlt)+kk-1).^2))^0.5);
            crossc=(sum(mat1.*mat2))/((sum(mat1.^2)*sum(mat2.^2))^0.5);
            % crossc=(sum(scantrc(kk:length(wvlt)+kk-1).*wvlt))/((sum(scantrc(kk:length(wvlt+kk-1)).^2)*sum(wvlt.^2))^0.5);
            warning on
            % set(wvltln2,'xdata',[xnew(jj) xnew(jj)],'ydata',[scansec(kk) (scansec(length(wvlt))+kk-1)],...
            %     'erasemode','xor');
            if(crossc>crosscor(1))
                if(crossc>=5)
                else
                    crosscor=[crossc scansec(kk+round(length(wvlt)/2)) xnum(jj)];
                end
            end
        end
    end
    holdcrosscor(hcc)=crossc;
    hcc=hcc+1;
    xdt=[xnew(jj) xnew(jj)];
    ydt=[t((ys)) t((y2))];
    % ln=line(xdt,ydt,'erasemode','xor','color','k');
    ys=y2;
    %     if(hcc==10)
    %         % plotting cross correlation lines
    %         if(sum(holdcrosscor)/5>0.40)
    %             % only using 65% crosscorrelation
    %             allcrosscor{length(allcrosscor)+1}=crosscor;
    %             dat1=allcrosscor{length(allcrosscor)-1};
    %             crosslnl=line([x(dat1(3)) x(crosscor(3))],[t(dat1(2)) t(crosscor(2))],'color','b',...
    %                 'erasemode','xor','linewidth',2);
    %         end
    %         hcc=1;
    %     end
    
    % percenttolerance is a user or default value which determines the
    % percentage crosscorelation that will be accepted as a node
    if(crosscor(1)>percenttolerance)
        % checking the largest 
        allcrosscor{length(allcrosscor)+1}=crosscor;
        dat1=allcrosscor{length(allcrosscor)-1};
        if(isempty(dat1))
        else
            yhold=[yhold crosscor(2)];
            xhold=[xhold crosscor(3)];
            perhold=[perhold crosscor(1)];
            if(~isempty(showprogress))
                crosslnl=line([x(dat1(3)) x(crosscor(3))],[t(dat1(2)) t(crosscor(2))],'color','b',...
                    'erasemode','xor','linewidth',2);
                lns=[lns crosslnl];
            end
        end
    else
    end
end
delete(lns);
% nodetolernace is a user or default
nodetol=nodetolerance;
if(length(yhold)<=1)
    % no matches have been found
    return
else
    % if the first and last x values are not at least 6% for the size away
    % from eachother, lines will be created
    if(~abs(xhold(end)-xhold(1))>=abs(x(end)-x(1))*0.06)
        return
    end
end
ii=ceil(length(yhold)/nodetol);
if(xstr==xend)
    yholdnew=[];
    xholdnew=[];
    yhold2=[];
    xhold2=[];
else
    % I am commenting out the following becuase, well just in case the the
    % user just wants a specfic area, not end of line
    yholdnew=[t(yhold(1))];
    xholdnew=[x(xhold(1))];
    yhold2=[yhold(1)];
    xhold2=[xhold(1)];
    %     yholdnew=[];
    %     xholdnew=[];
    %     yhold2=[];
    %     xhold2=[];
end


    for kk=1:2
        for jj=1:ii
            % this section is breaking down the data set to every nodetol
            if((nodetol*jj)>=length(yhold))
                ynd=round(length(yhold));
            else
                ynd=round(nodetol*jj);
            end
            if(jj==1)
                yst=1;
            else
                yst=round(nodetol*(jj-1));
            end
            yalign=yhold(yst:ynd);
            xalign=xhold(yst:ynd);
            if(jj~=ii)
                % this section is skipping forward to the next section to see if
                % there is enough numbers to make up its own, or if it has to be
                % incorporated into present section
                if((nodetol*jj)>=length(yhold))
                    yst=round(length(yhold));
                else
                    yst=round(nodetol*jj);
                end
                if(length(yhold(yst:end))<5)
                    % less then five, will take into present checking and halt
                    % after
                    yalign=yhold(yst:end);
                    xalign=xhold(yst:end);
                    ii=jj;
                end
            end
            ydatnew=t(yalign);
            xdatnew=x(xalign);
            if(size(xdatnew,1)~=size(ydatnew,1))
                ydatnew=ydatnew';
            end
            % doing a first order polynomial fit, good enough for a single line
            [pval,s]=polyfit(xdatnew,ydatnew,0);
            ydatnew=polyval(pval,xdatnew,s);
            yholdnew=[yholdnew ydatnew(end)];
            xholdnew=[xholdnew xdatnew(end)];
            yhold2=[yhold2 yalign(end)];
            xhold2=[xhold2 xalign(end)];
        end
        if(kk==1)
            % the above loop is going to go through again, so changeing all the
            % properties
            nodetol=nodetol/2;
            ii=ceil(length(yholdnew)/(nodetol));
            yhold=yhold2;
            xhold=xhold2;
            yholdnew=[t(yhold(1))];
            xholdnew=[x(xhold(1))];
        end
    end
    
pref.holddat=[yholdnew' xholdnew'];
set(gca,'userdata',pref);
% if(~isempty(yholdnew))
%     % ploting picks
%     picksnew=[];
%     crosslnl=line(xholdnew,...
%         yholdnew,'color','r',...
%         'erasemode','normal','linewidth',2,...
%         'buttondown',[]);
%     picksnew=[picksnew crosslnl];
% end
% 
% 
% The following is taking place to get by with not using "eval" in matlab.
% The line is being created because when it is deleted it will execute the
% the transfer fucntion.
delete(line([0 0],[0 0],'deletefcn',transfer));

%--------------
%- Data Check -
%--------------
%
% Just doing a quick data check everytime to make sure that data has some
% how not been lost from the GCA.  This section should be updated if more
% preferenaces are added.
function checkdata=PLE_checkdata(dat)
checkdata=[];
if(isempty(dat))
elseif(~ishandle(dat))
else
    pref.seisstruct=seisstruct;
    pref.axes=masteraxes;
    pref.transfer=trans;
    pref.holddat=[];
    pref.checkfields=checkfield;
    fnames=fieldnames(dat);
    for ii=1:size(checkfield)
        checked=[];
        for jj=1:size(fnames)
            if(strcmp(checkfield{ii},lower(fnames{jj})))
                checked='Yup';
                checkfield{ii}=getfield(preferances,fnames{jj});
                break
            end
        end
        if(isempty(checked))
            stringinfo={'Your preferances are not set up properly.',...
                    'Please see PICKSLE help for details in proper use',...
                    'or reset PICKSLE on your program'};
            helpdlg(stringinfo,'Check PICKSLE help file');
            return
        end
    end
end
checkdata=['Data is correct'];