function [G,MG] = make_graphe3D(xmin,xmax,nx,ymin,ymax,ny, obstacles, diag)

% on travaille sur le rectangle [xmin,xmax]x[ymin,ymax]
% en x : discretisation avec nx points (nx-1 intervalles)
% en y : discretisation avec m points (m-1 intervalles)
% obstacles : tableau de handles de fonctions representant
%             les obstacles. Chaque obstacle est défini par
%             une fonction F_k(x,y) qui renvoie vrai si le
%             point (x,y) est à l'extérieur de l'obstacle k
%             Exemples:
%              - pour un cercle de centre (cx,cy) et de rayon r
%                  F(x,y) = (x-cx)^2 + (y-cy)^2 > r
%              - pour un rectangle [a,b]x[c,d] :
%                  F(x,y) = ~(a <= x & x <= b & c <= y & y <= d)
%             De manière à être efficace on supposera que ces
%             fonctions sont codées pour pouvoir être évaluées
%             sur des vecteurs x, y
% diag : si vrai on ajoute les diagonales des sous-rectangles comme arêtes

% G: description du usuelle du graphe structure G.n G.succ,G.pred,G.c (et G.xy)
% MG: description par matrice d'adjacence du graphe (on a pas les coordonnées des sommets) 

x = linspace(xmin,xmax,nx);
y = linspace(ymin,ymax,ny);
[X,Y] = ndgrid(x,y);

% Definition de la surface
Z = cos(2*pi*X).*sin(pi*Y);
%Z = peaks(X,Y);

P = logical(ones(nx,ny));
for k = 1:numel(obstacles)
  P = P & obstacles{k}(X,Y);
end

% combien de points dans le graphe ?
n = sum(sum(P));

% init structure graphe
G.n = n;
G.succ = cell(n,1);
G.pred = cell(n,1);
G.c = cell(n,1);
G.xy = [X(P) Y(P) Z(P)];
% Num contient le numero des sommets du graphe (0 si le point n'est
% pas un sommet).
Num = zeros(nx,ny);
Num(P) = 1:n; 

% 2 petits tableaux de travail
succ = zeros(1,8); % 8 succ max (si on utilise les diagonales)
dist = zeros(1,8); % pour stocker les distances d'un point à ces successeurs

k = 1;
for j = 1:ny
  for i = 1:nx
    if P(i,j) % ce point existe : on l'ajoute dans le graphe
      if k ~= Num(i,j) 
	     error('Un probleme !!!!!!')
      end
      xy = G.xy(k,:); % coordonnees xy du sommet k
      ksucc = 0;
      if i > 1 && P(i-1,j)
    	% on rajoute le noeud en position (i-1,j) comme successeur
        ksucc = ksucc + 1; s = Num(i-1,j);
        succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
      end
      if i < nx && P(i+1,j)
        % on rajoute le noeud en position (i+1,j) comme successeur
        ksucc = ksucc + 1; s = Num(i+1,j);
        succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
      end
      if j > 1 && P(i,j-1)
        % on rajoute le noeud en position (i,j-1) comme successeur
        ksucc = ksucc + 1; s = Num(i,j-1);
        succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
      end
      if j < ny && P(i,j+1)
        % on rajoute le noeud en position (i,j+1) comme successeur
        ksucc = ksucc + 1; s = Num(i,j+1);
        succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
      end
      if diag  % on rajoute les arêtes diagonales
        if i>1 && j>1 && P(i-1,j-1)
            % on ajoute le noeud en position (i-1,j-1) comme successeur
            ksucc = ksucc + 1; s = Num(i-1,j-1);
            succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
        end
        if i>1 && j<ny && P(i-1,j+1)
            % on ajoute le noeud en position (i-1,j+1) comme successeur
            ksucc = ksucc + 1; s = Num(i-1,j+1);
            succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
        end
        if i<nx && j>1 && P(i+1,j-1)
            % on ajoute le noeud en position (i+1,j-1) comme successeur
            ksucc = ksucc + 1; s = Num(i+1,j-1);
            succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
        end
        if i<nx && j<ny && P(i+1,j+1)
            % on ajoute le noeud en position (i+1,j+1) comme successeur
            ksucc = ksucc + 1; s = Num(i+1,j+1);
            succ(ksucc) = s; dist(ksucc) = norm(G.xy(s,:) - xy);
        end
      end
      G.succ{k} = succ(1:ksucc);
      G.pred{k} = succ(1:ksucc); % les prédecesseurs sont les mêmes sommets !
      G.c{k} = dist(1:ksucc); % longueurs des arêtes
      k = k + 1;
    end
  end
end

% calcul de la matrice d'adjacence

%tableaux de travail (surdimensionnés...)
I = zeros(8*n,1); J = zeros(8*n,1); Val = zeros(8*n,1);

i = 0;
for k = 1:n
  succ = G.succ{k};
  ns = numel(succ);
  I(i+1:i+ns) = k*ones(ns,1);
  J(i+1:i+ns) = succ(:);
  Val(i+1:i+ns) = G.c{k}(:);
  i = i+ns;
end
MG = sparse(I(1:i),J(1:i),Val(1:i), n, n);
