function Ima2 = Poisson_solve(Ima, M)
% Ima2 = Poisson_solve(Ima, M)
% F. Sur 2015-2017
% Ima: image originale à restaurer dans Omega
% M masque: 0 dans Omega, 1 ailleurs
% Ima2: image restaurée par Poisson

Indices=find(M==0);
N=numel(Indices);
b=sparse(1,N);
A=sparse(N,N);

% méthode lente mais lisible
% for i=1:N,
%     if (mod(i-1,1000)==0), disp([num2str(i), ' / ',num2str(N)]); end;
%     AA=sparse(1,N);
%     AA(i)=-4;
%     [x,y]=ind2sub(size(Ima),Indices(i));
%     bb=0;
%     if M(x-1,y) bb=bb-Ima(x-1,y); else AA(Indices==sub2ind(size(Ima),x-1,y))=1; end;
%     if M(x+1,y) bb=bb-Ima(x+1,y); else AA(Indices==sub2ind(size(Ima),x+1,y))=1; end;
%     if M(x,y-1) bb=bb-Ima(x,y-1); else AA(Indices==sub2ind(size(Ima),x,y-1))=1; end;
%     if M(x,y+1) bb=bb-Ima(x,y+1); else AA(Indices==sub2ind(size(Ima),x,y+1))=1; end;
%     A(i,:)=AA;
%     b(i)=bb;
% end

% méthode "rapide" (à améliorer...)
MM=M.';
MtransYplus1=ones(size(M));
MtransYplus1(1:end-1)=M(2:end);
IndicesYplus1=find(MtransYplus1==0);
MtransYmoins1=ones(size(M));
MtransYmoins1(2:end)=M(1:end-1);
IndicesYmoins1=find(MtransYmoins1==0);
MtransXplus1=ones(size(MM));
MtransXplus1(1:end-1)=MM(2:end);
MtransXplus1=MtransXplus1.';
IndicesXplus1=find(MtransXplus1==0);
MtransXmoins1=ones(size(MM));
MtransXmoins1(2:end)=MM(1:end-1);
MtransXmoins1=MtransXmoins1.';
IndicesXmoins1=find(MtransXmoins1==0);
Imat=Ima.';
ImatransYplus1=ones(size(Ima));
ImatransYplus1(1:end-1)=Ima(2:end);
ImatransYmoins1=ones(size(Ima));
ImatransYmoins1(2:end)=Ima(1:end-1);
ImatransXplus1=ones(size(Imat));
ImatransXplus1(1:end-1)=Imat(2:end);
ImatransXplus1=ImatransXplus1.';
ImatransXmoins1=ones(size(Imat));
ImatransXmoins1(2:end)=Imat(1:end-1);
ImatransXmoins1=ImatransXmoins1.';
bIma=-(MtransYplus1.*ImatransYplus1+MtransYmoins1.*ImatransYmoins1+MtransXplus1.*ImatransXplus1+MtransXmoins1.*ImatransXmoins1);
for i=1:N,
    if (mod(i-1,1000)==0), disp([num2str(i), ' / ',num2str(N)]); end;
    A(i,Indices==IndicesXplus1(i))=1;
    A(i,Indices==IndicesXmoins1(i))=1;
    A(i,Indices==IndicesYplus1(i))=1;
    A(i,Indices==IndicesYmoins1(i))=1;
end;
A=A-4*speye(N);
b=transpose(bIma(Indices));


Ima2=Ima;
Ima2(Indices)=b/A;
end
