function D2I = toI2( Fp, tri, Ip)

# Usage: D2I = toI2( Fp, tri, Ip)
# computes a linear interpolation from Fp to Ip based 
# on a triangulation tri
# this is equivalent to using griddata

[si,sr] = size(Ip);
scal = 1/abs(Ip(1)-Ip(2));
Ip1 = Ip(1) *scal;
# optimisation: determine rectangle of points possibly inside
maxr = floor(max(real(Fp(tri')))'*scal-real(Ip1));
maxi = floor(max(imag(Fp(tri')))'*scal-imag(Ip1));
minr = ceil(min(real(Fp(tri')))'*scal-real(Ip1));
mini = ceil(min(imag(Fp(tri')))'*scal-imag(Ip1));
maxr = min(maxr,sr-1);
maxi = min(maxi,si-1);
minr = max(minr,0);
mini = max(mini,0);

tmpr = arrayfun(@(a,b)  a:b,minr,maxr,"UniformOutput",false);
tmpi = arrayfun(@(a,b)  a:b,mini,maxi,"UniformOutput",false);
tmp = cellfun(@(a,b) vec(si*a(ones(size(b)),:)+b(ones(size(a)),:)')+1,tmpr,tmpi,"UniformOutput",false);

tmploc = cellfun(@(a,b) vec(a(ones(size(b)),:)+1i*b(ones(size(a)),:)'),tmpr,tmpi,"UniformOutput",false);
tmptri =mat2cell(Fp(tri)*scal-Ip1+0.000011+0.00001i,ones(1,size(tri)(1)),size(tri)(2));

tmpb =  cellfun(@baryc0,tmptri,tmploc,"UniformOutput",false);
beta = cell2mat(cellfun(@vec,tmpb,"UniformOutput", false));
beta=reshape(beta,3,[]);
sel = all(beta>=0,1);
beta0 = cellfun("size",tmpb,2);
tmpl = cell2mat(tmploc);
tmpl = real(tmpl)*si + imag(tmpl)+1;

D2I = sparse(tmpl(sel),repelems(tri(:,1),[1:length(tmptri);beta0'])(sel),beta(1,sel),length(vec(Ip)),length(Fp));
D2I += sparse(tmpl(sel),repelems(tri(:,2),[1:length(tmptri);beta0'])(sel),beta(2,sel),length(vec(Ip)),length(Fp));
D2I += sparse(tmpl(sel),repelems(tri(:,3),[1:length(tmptri);beta0'])(sel),beta(3,sel),length(vec(Ip)),length(Fp));

#D2I= sparse(kron(tmpl(sel),[1 1 1]), repelems(tri),[1:length(tmptri);beta0']) , beta',i512^2,length(Fp));
endfunction
