function zp=griddata_grid_to_mesh_completer(x_m,y_m,u,index,x,y,zp);
%GRIDDATA_GRID_TO_MESH_COMPLETER - interpolates data from a grid to a mesh   
%
%   Usage:
%      zp=griddata_grid_to_mesh_completer(x_m,y_m,u,index,x,y,zp)

pos=find(isnan(u));
for n=1:100,
   pos=find(isnan(u));
   u(pos)=u(pos+1);
end

imagesc(u)
break

zp=zeros(length(x),1);
for n=1:length(x),
   if mod(n,100)==0,
      disp(['Count is ' num2str(n/length(x)*100)]);
   end
   
   pos1=find((x_m-x(n))<0);
   pos2=find((x_m-x(n))>0);
   pos3=find((y_m-y(n))<0);
   pos4=find((y_m-y(n))>0);
   
   if (isempty(pos1) | isempty(pos2) | isempty(pos3) | isempty(pos4)),
      out_of_bounds=1;
   else
      out_of_bounds=0;
   end
   
   if out_of_bounds==0,   
   %first find where (x(n),y(n)) lies.
   pos1=find((x_m-x(n))<0);
   pos2=find((x_m-x(n))>0);
   x_bound=[pos1(length(pos1)) pos2(1)];   
   
   pos1=find((y_m-y(n))<0);
   pos2=find((y_m-y(n))>0);
   y_bound=[pos1(length(pos1)) pos2(1)];   
   
   
   %we look for a linear interpolation, therefore have to decide 
   % in which triangle our point lies.
   
   A=[x_m(x_bound(1)) y_m(y_bound(1))]';
   u_A=u(x_bound(1),y_bound(1));
   B=[x_m(x_bound(2)) y_m(y_bound(1))]';
   u_B=u(x_bound(2),y_bound(1));
   C=[x_m(x_bound(2)) y_m(y_bound(2))]';
   u_C=u(x_bound(2),y_bound(2));
   D=[x_m(x_bound(1)) y_m(y_bound(2))]';
   u_D=u(x_bound(1),y_bound(2));
   
   
   
   res=(D(2)-B(2))*x(n)+(B(1)-D(1))*y(n)+D(1)*B(2)-D(2)*B(1);
   if res>0,
      xx=[B(1) C(1) D(1)];
      yy=[B(2) C(2) D(2)];
      uu=[u_B u_C u_D];
      ind=[1 2 3];
   else
      xx=[A(1) B(1) D(1)];
      yy=[A(2) B(2) D(2)];
      uu=[u_A u_B u_D];
      ind=[1 2 3];
   end
   
   
	alpha=zeros(1,3);
	beta=zeros(1,3);
	gamma=zeros(1,3);
	X=inv([xx(ind(1,:))' yy(ind(1,:))' ones(3,1)]);
	alpha(1,:)=X(1,:);
	beta(1,:)=X(2,:);
	gamma(1,:)=X(3,:);
   
   if length(find(isnan(uu)))==0,
      zp(n)=sum(uu.*(alpha.*x(n)+beta.*y(n)+gamma));
   else
      if length(find(~isnan(uu)))==0,
         zp(n)=NaN;
         else
         pos=find(~isnan(uu));
         zp(n)=uu(pos(1));
      end
   end
   
   else
   zp(n)=NaN;
end

   
end
%disp('Filling final gaps');
%zp=void_fill(zp,x,y);
