function  ret = HLap(H9L,Fp)

# Usage:  ret = HLap(H9L,Fp)
# correction for raw foveate Laplacian
# make H9L zero for linear, and const for quadratic form on Fp

[i,j,v]=find(H9L);
Dxx = sparse(i,1:length(j),Fp(j));
Dxx = [real(Dxx);imag(Dxx);spones(Dxx);abs(Dxx).^2;real(Dxx).^2-imag(Dxx).^2];

#v0 = H9L*vec(Fp);
v0 = accumarray(i,v.*Fp(j));
v0=[real(v0);imag(v0);sum(H9L,2);H9L*vec(abs(Fp).^2)+1;H9L*(real(Fp).^2-imag(Fp).^2)];
v -= Dxx \ v0;
H9L = sparse(i,j,v);
v0 = accumarray(i,v.*Fp(j));
v0=[real(v0);imag(v0);sum(H9L,2);H9L*vec(abs(Fp).^2)+1;H9L*(real(Fp).^2-imag(Fp).^2)];
v -= Dxx \ v0;
ret = sparse(j,i,v);
endfunction
