griddata_grid_to_mesh_nan

PURPOSE ^

GRIDDATA_GRID_TO_MESH_NAN - interpolates data from a grid to a mesh

SYNOPSIS ^

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

DESCRIPTION ^

GRIDDATA_GRID_TO_MESH_NAN - interpolates data from a grid to a mesh   

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

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