clear all
close all

V = gingerbread_LD_contour();
n = size(V,1);

N=40;
res = [N,N];
nbPixels = prod(res);


minX = min(V(:,1));
maxX = max(V(:,1));
minY = min(V(:,2));
maxY = max(V(:,2));

V = (V - kron([minX, minY], ones(n,1))) ./ kron([maxX, maxY]-[minX, minY], ones(n,1));
V = V .* kron(res-1, ones(n,1));

data = poly2mask(V(:,1), V(:,2), res(1), res(2));


index_fixed = reshape(data, 1, nbPixels) .* (1:nbPixels);
index_fixed = index_fixed(find(index_fixed));

index_free = reshape(1-data, 1, nbPixels) .* (1:nbPixels);
index_free = index_free(find(index_free));



nbFixedPixels = sum(data(:));
nbFreePixels = nbPixels - nbFixedPixels;


A = sparse(nbFreePixels, nbFreePixels);
B = sparse(nbFreePixels, nbFixedPixels);

BB = zeros(nbFreePixels, 1);

for k=1:nbFreePixels
  [i,j] = ind2sub(res, index_free(k));
  
  nbA = 0;
  
  if i>1
    ii = sub2ind(res, i-1, j);
    nbA = nbA + 1;
    if data(ii) == 0
      ii = find(index_free-ii == 0)(1);
      A(k, ii) = A(k, ii) - 1;
    else
      ii = find(index_fixed-ii == 0)(1);
      B(k, ii) = B(k, ii) + 1;
      BB(k) = BB(k)+1;
    end
  end
  
  if i<res(1)
    ii = sub2ind(res, i+1, j);
    nbA = nbA + 1;
    if data(ii) == 0
      ii = find(index_free-ii == 0)(1);
      A(k, ii) = A(k, ii) - 1;
    else
      ii = find(index_fixed-ii == 0)(1);
      B(k, ii) = B(k, ii) + 1;
      BB(k) = BB(k)+1;
    end
  end
  
  if j>1
    ii = sub2ind(res, i, j-1);
    nbA = nbA + 1;
    if data(ii) == 0
      ii = find(index_free-ii == 0)(1);
      A(k, ii) = A(k, ii) - 1;
    else
      ii = find(index_fixed-ii == 0)(1);
      B(k, ii) = B(k, ii) + 1;
      BB(k) = BB(k)+1;
    end
  end
  
  if j<res(2)
    ii = sub2ind(res, i, j+1);
    nbA = nbA + 1;
    if data(ii) == 0
      ii = find(index_free-ii == 0)(1);
      A(k, ii) = A(k, ii) - 1;
    else
      ii = find(index_fixed-ii == 0)(1);
      B(k, ii) = B(k, ii) + 1;
      BB(k) = BB(k)+1;
    end
  end
  
  A(k,  k) = nbA;
end

for k=1:nbFreePixels
  B(k, :) = B(k, :) / BB(k);
end

X = rand(nbFixedPixels, 1);

%X = (1:nbFixedPixels)';

I = A \ (B*X);


%borderI = find(B*ones(nbFixedPixels,1));
%
%I = zeros(nbFreePixels,1);
%lambda = 0.05;
%for i=1:10000
%  I(borderI) = (B*X)(borderI);
%  I = I - lambda*A*I;
%end


XX = zeros(res);
XX(index_fixed) = X;

II = zeros(res);
II(index_free) = I;


figure()
subplot(1,2,1)
imagesc(XX)
axis image xy
colorbar

subplot(1,2,2)
imagesc(II)
axis image xy
colorbar


%figure()
%imagesc(XX+II)
%axis image xy
%colorbar


