griddata_grid_to_mesh

PURPOSE ^

GRIDDATA_GRID_TO_MESH - interpolates data from a grid to a mesh

SYNOPSIS ^

function zp=griddata_grid_to_mesh(x_m,y_m,u,index,x,y);

DESCRIPTION ^

GRIDDATA_GRID_TO_MESH - interpolates data from a grid to a mesh   

   Usage:
      zp=griddata_grid_to_mesh(x_m,y_m,u,index,x,y)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function zp=griddata_grid_to_mesh(x_m,y_m,u,index,x,y);
0002 %GRIDDATA_GRID_TO_MESH - interpolates data from a grid to a mesh
0003 %
0004 %   Usage:
0005 %      zp=griddata_grid_to_mesh(x_m,y_m,u,index,x,y)
0006 
0007 pos=find(isnan(u));
0008 for n=1:10,
0009    pos=find(isnan(u));
0010    u(pos)=u(pos+1);
0011 end
0012 
0013 zp=zeros(length(x),1);
0014 for n=1:length(x),
0015    if mod(n,100)==0,
0016       disp(['Count is ' num2str(n/length(x)*100)]);
0017    end
0018    
0019    pos1=find((x_m-x(n))<0);
0020    pos2=find((x_m-x(n))>0);
0021    pos3=find((y_m-y(n))<0);
0022    pos4=find((y_m-y(n))>0);
0023    
0024    if (isempty(pos1) | isempty(pos2) | isempty(pos3) | isempty(pos4)),
0025       out_of_bounds=1;
0026    else
0027       out_of_bounds=0;
0028    end
0029    
0030    if out_of_bounds==0,   
0031    %first find where (x(n),y(n)) lies.
0032    pos1=find((x_m-x(n))<0);
0033    pos2=find((x_m-x(n))>0);
0034    x_bound=[pos1(length(pos1)) pos2(1)];   
0035    
0036    pos1=find((y_m-y(n))<0);
0037    pos2=find((y_m-y(n))>0);
0038    y_bound=[pos1(length(pos1)) pos2(1)];
0039   
0040    
0041    
0042    %we look for a linear interpolation, therefore have to decide
0043    % in which triangle our point lies.
0044    A=[x_m(x_bound(1)) y_m(y_bound(1))]';
0045    u_A=u(x_bound(1),y_bound(1));
0046    B=[x_m(x_bound(2)) y_m(y_bound(1))]';
0047    u_B=u(x_bound(2),y_bound(1));
0048    C=[x_m(x_bound(2)) y_m(y_bound(2))]';
0049    u_C=u(x_bound(2),y_bound(2));
0050    D=[x_m(x_bound(1)) y_m(y_bound(2))]';
0051    u_D=u(x_bound(1),y_bound(2));
0052    
0053    
0054    
0055    res=(D(2)-B(2))*x(n)+(B(1)-D(1))*y(n)+D(1)*B(2)-D(2)*B(1);
0056    if res>0,
0057       xx=[B(1) C(1) D(1)];
0058       yy=[B(2) C(2) D(2)];
0059       uu=[u_B u_C u_D];
0060       ind=[1 2 3];
0061    else
0062       xx=[A(1) B(1) D(1)];
0063       yy=[A(2) B(2) D(2)];
0064       uu=[u_A u_B u_D];
0065       ind=[1 2 3];
0066    end
0067    
0068    
0069     alpha=zeros(1,3);
0070     beta=zeros(1,3);
0071     gamma=zeros(1,3);
0072     X=inv([xx(ind(1,:))' yy(ind(1,:))' ones(3,1)]);
0073     alpha(1,:)=X(1,:);
0074     beta(1,:)=X(2,:);
0075     gamma(1,:)=X(3,:);
0076    
0077    if length(find(isnan(uu)))==0,
0078       zp(n)=sum(uu.*(alpha.*x(n)+beta.*y(n)+gamma));
0079    else
0080       if length(find(~isnan(uu)))==0,
0081          zp(n)=NaN;
0082          else
0083          pos=find(~isnan(uu));
0084          zp(n)=uu(pos(1));
0085       end
0086    end
0087    
0088    else
0089    zp(n)=NaN;
0090 end
0091 
0092    
0093 end
0094 %disp('Filling final gaps');
0095 %zp=void_fill(zp,x,y);

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003