function [h ] = lbp( imga , P , R )
% ENTREE
%[h, LBP]
% imga: l'image en entrée
% P: nombre de voisins
% R: rayon déterminant la résolution spatiale
% SORTIE
% h: histo normalisé
% LBP: l'image LBP
% PROGRAMME
n=500;% taille de l'image normalisé
m=500;
img=imresize(imga,[n,m],'bicubic'); % Normalisation
[nl nc]=size(img); % taille de l'image
% Entourer l'image d'entrée de zéros
img=[ ones(R,nc+2*R)*255;ones(nl,R)*255 img ones(nl,R)*255;ones(R,nc+2*R)*255];
u=1;
for i=(R+1):nl+R
v=1;
for j=(R+1):nc+R
LBP(u,v)=0;
s=zeros(1,P);
gc=img(i,j); % pixel central
for p=0:P-1 % Je parcoure le voisinage
x=i+R*sin((2*pi*p)/P);
y=j-R*cos((2*pi*p)/P);
% Je vérifie si une interpolation est nécessaire
if (x-floor(x)~=0 && y-floor(y)~=0)
d1x=x-floor(x);
d1y=y-floor(y);
d1=sqrt(d1x^2+d1y^2);
g1=img(floor(x),floor(y));
d2x=x-floor(x);
d2y=y-floor(y)+1;
d2=sqrt(d2x^2+d2y^2);
g2=img(floor(x),floor(y)+1);
d3x=x-floor(x)+1;
d3y=y-floor(y)+1;
d3=sqrt(d3x^2+d3y^2);
g3=img(floor(x)+1,floor(y)+1);
d4x=x-floor(x)+1;
d4y=y-floor(y);
d4=sqrt(d4x^2+d4y^2);
g4=img(floor(x)+1,floor(y));
% Valeur du pixel voisin interpolé:
gp=((g1*d1)+(g2*d2)+(g3*d3)+(g4*d4))/(d1+d2+d3+d4);
else % sinon, pas besoin d'interpoller
x=floor(x);
y=floor(y);
gp=img(x,y);
end
if(gp>=gc)
s(1,p+1)=1;
else
s(1,p+1)=0;
end
end
for p=0:P-1
LBP(u,v)=LBP(u,v)+s(1,p+1)*2^p;
end
v=v+1;
end
u=u+1;
end
%% Calcul de l'histogramme normalisé
for i=1:2^P
h(i)=0;
end
for i=1:nl
for j=1:nc
a= LBP(i,j);
a=double(a);
h(a+1)=h(a+1)+1;
end
end
h=h/(nl*nc); % normalisation
end