# this script generates the examplesa for Riga talk
NB = 20;
[Fp,H9] = hexnet(NB,64);
Fp /= 1.23;

lz = length(Fp);
lF = (lz-1)/6;
#gplot(H9,[real(Fp),imag(Fp)])

# find neighbours
for j= 1:lz;
  nbj= find(H9(j,:));
  Fdj = Fp(nbj) - Fp(j);
  if(length(nbj) == 4 ) nbj=[nbj,0]; Fdj=[Fdj;Fp(j)/abs(Fp(j))]; endif
  [~,ix] = sort( arg(Fdj)); 
  nb{j}= nbj(ix);
endfor

function b = baryc0(a,p)
# compute barycentric coordinateso of p in triangle a
  p = p.';
  b = [real(a(2)-p).*imag(a(3)-p)-real(a(3)-p).*imag(a(2)-p);
    real(a(3)-p).*imag(a(1)-p)-real(a(1)-p).*imag(a(3)-p);
    real(a(1)-p).*imag(a(2)-p)-real(a(2)-p).*imag(a(1)-p)] / ...
	(real(a(2)-a(1)).*imag(a(3)-a(1))-real(a(3)-a(1)).*imag(a(2)-a(1)));
endfunction

function a = Ar(dexs, Fp0, Fp)
 if(all(dexs)); F = Fp(dexs)-Fp0 ;
 else; F = ifelse((dexs>0)',Fp(max(1,dexs))-Fp0,Fp0/abs(Fp0)*ones(5,1)) ; endif;
 a = sum(imag(([F(2:end);F(1)].*conj(F)) -F.*conj([F(2:end);F(1)])));
endfunction

function a = Ar3(Fp0)
 F1 = Fp0(:,1)-Fp0(:,3);
 F2 = Fp0(:,2)-Fp0(:,3);
 a = abs(real(F1).*imag(F2)-imag(F1).*real(F2));
endfunction

Fpa = cellfun("Ar",nb,num2cell(Fp)',{Fp}) / 12;

tri= trigon(H9);

# zugehöriger Laplace-Operator
H9L = H9;
H9L = (diag(sum(H9L)) -H9L)/6;

# Entfernen der Fünfeckmittelpunkte
Dred = sparse(eye(length(Fpa)));
if(1)
tmp = find(sum(H9)==5);
for i = tmp; Dred(i,nb{i}) = 0.2; Dred(i,i) = 0; endfor
endif

# TS is the area element for teh positions
TSqi = sqrt(Fpa'/655.36);
TS = 1./Fpa'/655.36;
TS(isinf(TS)) = 0;
TSq = sqrt(TS);
H9L = diag(TSq)*H9L*diag(TSq)*655.36/1.15;
H9L = Dred' * H9L *Dred;
H9L = HLap(H9L,Fp);
Fpa = Fpa*Dred;
TSqi = sqrt(Fpa'/655.36);
TS = 1./Fpa'/655.36;
TS(isinf(TS)) = 0;
TSq = sqrt(TS);
# Hamilton-Funktion H = H9L + diag(abs(Fp*pi/10).^2)

# read the default picture, if I0 is not defined
if( ! exist( "I0") ) 
I=imread("lena.ppm") ;
I0 = double(I)/256-0.5;
endif
i512 = max(size(I0));

function xydex = nextxy ( Fp,Ip )
# nextxy finds nearest neighbour in a Cartesian grid
i512 = size(Ip)(1);
Nx = i512*abs(Ip(2)-Ip(1));
sr= roundb((real(Fp)/Nx+0.5)*i512-0.5);
si= roundb((imag(Fp)/Nx+0.5)*i512-0.5);
si(si> i512-1) = i512-1; si(si<0) =0;
sr(sr> i512-1) = i512-1; sr(sr<0) = 0;
xydex = 1+si+i512*sr;
endfunction

Nx = ceil(NB*sqrt(4/3));# Nx2 = 2*Nx+1;
Nx = 20;
rr= (real(Fp)/Nx+0.5)*i512; sr= floor(rr); rr -= sr;
ri= (imag(Fp)/Nx+0.5)*i512; si= floor(ri); ri -= si;
si(si> i512-2) = i512-2; si(si<0) =0;
sr(sr> i512-2) = i512-2; sr(sr<0) = 0;
I2D = sparse( 1:length(Fp), 1+si+i512*sr, (1-ri).*(1-rr),length(Fp),length(vec(I0)));
I2D += sparse( 1:length(Fp), 1+si+1+i512*sr, (1-ri).*rr,length(Fp),length(vec(I0)));
I2D += sparse( 1:length(Fp), 1+si+i512*(sr+1), (1-rr).*ri,length(Fp),length(vec(I0)));
I2D += sparse( 1:length(Fp), 1+si+1+i512*(sr+1), ri.*rr,length(Fp),length(vec(I0)));

Id = I2D * vec(I0);

#v24FF;

#It = mean( Id(tri),1).';
Ipmax = Nx/2;
#Ip = ones(i512,1) * ((0:i512-1)-floor(i512/2)+0.5);
Ip = (1:i512)-(i512+1)/2;
Ip += 1i* Ip';
Ip *= Ipmax/256;
I2D1 = sparse( 1:length(Fp), nextxy ( Fp,Ip ), 1,length(Fp),length(vec(I0)));
D2I= toI2(Fp,tri,Ip);
#D2IL = D2I *H9L';
#D2Ix = D2I + diag((1-sqrt(sumsq(D2I,2)))./sqrt(sumsq(D2IL,2)))*D2IL;

Fpsel = find(sum(D2I)>0);

[~,~,H4] = laplacian(size(I0));
magic7853=length(Fp);
H4 *= 0.25*655.36;
D2Ix = interpj(Fp,Ip,H9L*H9L*D2I',H9)';
nshow(reshape(vec(I0)'*D2Ix*diag(TS)*D2Ix',size(I0)))
# try to improve teh interpolation
D2Ix = smoothie(D2Ix,H9L,H4);
#D2Ix = ortho(D2Ix,TS);
