/////////////////////////
//
//  Executer ce programme  avec la commande
//
//       exec('gibbs.sci',-1); 
//
//  puis activer le bouton 'Champ'  ou 'Texture' de la fenetre graphique créée. 
//
////////////////////////
function m=gibbs(b,n,ng)
//
// Simulation d'un champ de Markov a deux niveaux de gris avec 4-voisinage.
//
// b  (def=[1;-1]) :  Les deux coefficients,
// n  (def=64)     :  Nombre de points d'un cote de la grille (carree).
// ng (def=25)     :  Nombre d'iterations de Gibbs 
//                      (1 iteration=passage de tous les points).

[lhs,rhs]=argn(0)
if rhs<3 then; ng=25;end;
if rhs<2 then; n=64;end;
if rhs<1 then; b=[1;-1];end;

// creation du booleen damier
bidon=ones(1,n);
bidon(1:2:n)=-1;
chess=bidon'*bidon;
chess(:,1)=0;
chess(:,n)=0;
chess(1,:)=0;
chess(n,:)=0;
chess=chess==1;
nchess=~chess;
//initialisation et conditions de bord constantes
m=2*round(rand(n,n))-1; 
bord=0
m(1,:)=bord;
m($,:)=bord;
m(:,1)=bord;
m(:,$)=bord;
//iterations de Gibbs
b1=b(1); b2=b(2);
lc=2:n-1;
ld=lc-1;
lg=lc+1;
lt=1:n;
for i=1:ng
 p=b1*m(lc,ld)+b1*m(lc,lg)+b2*m(ld,lc)+b2*m(lg,lc);
 p=exp(p);
 p=p./(1+p);
 p=[lc' p lc']; p=[lt;p;lt];
 z=rand(n,n);
 m(p>z & chess)=1; 
 m(p<=z & chess)=-1;
 p=b1*m(lc,ld)+b1*m(lc,lg)+b2*m(ld,lc)+b2*m(lg,lc);
 p=exp(p);
 p=p./(1+p);
 p=[lc' p lc']; p=[lt;p;lt];
 m(p>z & nchess)=1; 
 m(p<=z & nchess)=-1;
end;
endfunction

function plotmar(m)
x=1:size(m,1);y=1:size(m,2);
xset("colormap",graycolormap(2));xbasc();
//Matplot(m+1,'022');
[u,v]=find(m<0);
plot2d([1;size(m,1)+1],[0;size(m,2)],0,'020')  
u=[u;v;ones(u);ones(v)+.1]
xrects(u,ones(v))
endfunction

function texture(b,n,ng)
// creation d'une texture par imbrication de champs.
m1=gibbs(b,n,ng);
xbasc();xtitle('Premier champ');plotmar(m1);
m2=gibbs(b,n,ng)
m3=zeros([m1 m2]); 
m3(1:2:$,1:2:$)=m1(1:$/2,:); 
m3(2:2:$,1:2:$)=m1($/2+1:$,:); 
m3(:,2:2:$)=m2';
xbasc();xtitle('Texture');plotmar(m3);
endfunction


xset('font',2,3);
addmenu(0,'Champ');
str=['Taille (nb  de points de coté)';'Nombre d''itérations de Gibbs';...
  'Coefficient horizontal';'Coefficient vertical'];
defc=['128';'25';'1.5';'1.5'];
Champ_0(1)='y=x_mdialog(''Champ simple'',str,defc);y=evstr(y);.. 
  m=gibbs(y(3:4),y(1),y(2));xtitle(''Champ de Markov'');plotmar(m);'

addmenu(0,'Texture');
deft=['128';'20';'1';'1'];
Texture_0(1)='y=x_mdialog(''Champ avec texture'',str,deft);y=evstr(y);.. 
  texture(y(3:4),y(1),y(2));'

addmenu(0,'HelpGibbs');
str1=['Champ de gibbs de caractéristiques suivantes :';..
'  Chaque point de la grille nXn vaut 0 ou 1';..
'  Interactions aux quatre plus proches voisins avec coefficients b1 et b2';..
'  Le potentiel est';..
'     Somme { x(i,j)[b1.x(i+1,j)+b2.x(i,j+1)] }';..
'b1 est le coefficient horizontal';..
'n est la taille de la grille';..
'';..
'Les textures sont faites en fabricant deux champs nXn indépendants';..
'et en les imbriquant (voir le programme), d''où un champ (2n)Xn.';..
'Le premier champ est affiché.']
HelpGibbs_0(1)='x_dialog(str1)';