function [Q,cortex,base,csf,Pt]=vdp2(ImRef,ImDeg,varargin)
%ImRef et ImDeg entre 0 et 1
%ImRef et ImDeg doivent être des images carrés de même tailles.
%Pour accélerer lorsque l'image a toujours la même taille on 
%peut utiliser la fonction Q=vdp2(ImRef,ImDeg,cortex,base,csf); 
%où cortex, base et csf ont été calculés lors d'une précédente utilisation
%Sinon l'utilisation classique est Q=vdp2(ImRef,ImDeg);
%La carte de visibilité Pt est disponible avec la 
%[Q,cortex,base,csf,Pt]=vdp2(ImRef,ImDeg);
%modification pour la version matlab 6.5.1 : 
 %ifft2 avec l'option symmétrique n'existe apparamment pas

n=size(ImRef,1);
if (nargin>=5) cortex=varargin{1}; base=varargin{2}; csf=varargin{3};
else [cortex,base]=Dalycortex(n); csf=Dalycsf(n);
end;
if (max(max(abs(ImRef-ImDeg)))<1e-14) Q=0;
else
  if (n-2*floor(n/2)~=0) sprintf('le programme ne fonctionne pas pour des images de tailles impaires'), Q=NaN; 
  else 
    Lumimag=DalyLum(ImRef);
  imagettes=FiltrageCortex(Lumimag,cortex,base,n,csf);
  bk=mean(mean(imagettes(1,:,:)));   % la luminance moyenne de l image de base
  beta=2; k1=exp(-7*log(6)/3); k2=exp(10*log(6)/3); s=0.7; b=4;
  tk1=(1+(k1*(k2*abs(imagettes)).^s).^b).^(1/b);
  Lumimag=DalyLum(ImDeg);
  imagettes2=FiltrageCortex(Lumimag,cortex,base,n,csf);
  tk2=(1+(k1*(k2*abs(imagettes2)).^s).^b).^(1/b);
  ckl=abs((imagettes- imagettes2)/bk);
  tkem=min(tk1,tk2);
  Pk=1-exp(-(ckl./tkem).^beta);
  Pt=1-prod(1-Pk,1);
  Q=sum(sum(Pt));
  end;
end;

function [cortex,base]=Dalycortex(n);
% function [cortex,base]=Dalycortex;
% cortex(k,l,x,y) k: dom l: fan  Les 30 filtres cortex
% modification : cortex(x,y,k,l)
% base: le filtre bande de base.
% n: dimension de l'image nombre de pixels
%
%
%
% voir aussi:
%           FiltrageCortex    Decompositioncortex    Dalycsf    DalyLuminance   Daly   ListeDaly   

[x,y]=meshgrid(1:n,1:n);
x1=(x-(n+1)/2); y1=(y-(n+1)/2);
theta=(angle(x1+sqrt(-1)*y1))*180/pi;
%on souhaite des fans symétriques
theta=theta-180*(theta>90)+180*(theta<-90);
f=(x1.^2+y1.^2)/(n/2)^2;


%fk(1:5)=2/3./2.^(0:4);
%tw=2/3*fk;

%Nombre de bandes
K=5;
% filtre Bande de Base: BaseBand
rho=1/2^(K); tw=2/3*rho; sigma=1/3*(rho+tw/2);
base=exp(-f.^2/2/sigma^2).*(f<rho+tw/2);

% filtres Dom: 5 filtres
rho=1/2^(K-1); tw=2/3*rho; sigma=1/3*(rho+tw/2);
mesa=zeros(size(x,1),size(y,1),5);
dom=zeros(size(x,1),size(y,1),5);
mesa(:,:,5)=(f<rho-tw/2)+(1+cos(pi*((f-(rho-tw/2))/tw)))/2.*(f>=rho-tw/2).*(f<rho+tw/2);
dom(:,:,5)=mesa(:,:,5)-base;

for i=1:4,
  rho=1/2^(K-1-i); tw=2/3*rho;
  mesa(:,:,5-i)=(f<rho-tw/2)+(1+cos(pi*((f-(rho-tw/2))/tw)))/2.*(f>=rho-tw/2).*(f<rho+tw/2);
  dom(:,:,5-i)=mesa(:,:,5-i)-mesa(:,:,6-i);
end

% filtres FAN
delta=30;
thetac1=-90;
thetac2=90;
fan(:,:,1)=(1+cos(pi*(abs(theta-thetac1)/delta)))/2.*(abs(theta-thetac1)<=delta)+(1+cos(pi*(abs(theta-thetac2)/delta)))/2.*(abs(theta-thetac2)<=delta);
for l=2:6,
  thetac=((l-1)*delta)-90;
  fan(:,:,l)=(1+cos(pi*(abs(theta-thetac)/delta)))/2.*(abs(theta-thetac)<=delta);
 % [ind1,ind2]=find(abs(teta-tetac)>delta); 
 % if (length(ind1)>0) fan(ind1,ind2,l)=0; end;
 % [ind1,ind2]=find(abs(teta-tetac)<=delta); 
 %[length(ind1) l], 
 % if (length(ind1)>0) 
 %   fan(ind1,ind2,l)=(1+cos(pi*(abs(teta(ind1,ind2)-tetac)/delta)))/2;
 % end;
end 

 % Filtres Cortex
for k=1:5,
  for l=1:6,
    cortex(:,:,k,l)=fan(:,:,l).*dom(:,:,k);
  end
end

    


function [Lumimag]=DalyLum(im);
% function [Lumimag]=Dalyluminance(nom);
% nom: fichier image format BMP.
% r/rmax=l/(l+(12.6*l)^0.63)
% l=0.3081*ng+0.1646
%
%
%
% voir aussi:
%           FiltrageCortex    Dalycortex    Dalycsf    Decompositioncortex   Daly   ListeDaly   
im=256*im;
%n=size(im);
%dimx=n(1);
%dimy=n(2);
%rmax=0;
l=0.3081*im+0.1646;
r=l./(l+(12.6*l).^0.63);
Lumimag=r/max(r(:));
%for x=1:dimx,
%    for y=1:dimy,
%        l=0.3081*double(im(x,y))+0.1646;
%        r(x,y)=l/(l+exp(0.63*log(12.6*l)));
%        if r(x,y)>rmax rmax=r(x,y); end
%    end
%end
%Lumimag=r/rmax;

function [csf]=Dalycsf(n);
% function [csf]=Dalycsf(n);
% n: nombre de pixels
%
%
% voir aussi:
%           FiltrageCortex    Dalycortex    Decompositioncortex    DalyLuminance   Daly   ListeDaly   
p=250;
wm=0.076;   % la taille de l image en métres: 256x256 ---> 7.6 cm.
%wm=0.076/256*n;     %modification pour s'adapter à une variation de n
  %modification annulée
d=6*wm;     % la distance de visualisation égale six fois la taille de l image.
w=360*atan(1/12)/pi;   % la taille de l image en degré visuel.
la=100;      % la luminance d adaptation (cd/m^2)
ec=0;
al=0.801*(1+(0.7/la))^(-0.2);
bl=0.3*(1+(100/la))^0.15;
ra=0.856*d^0.14;
re=1/(1+0.24*ec);
%al=0.801*exp(-0.2*log(1+(0.7/la)));
%bl=0.3*exp(0.15*log(1+(100/la)));
%ra=0.856*exp(0.14*log(d));
%re=1/(1+0.24*ec);
[x,y]=meshgrid(1:n,1:n);
x1=x-((n+1)/2)+1e-2;
y1=y-((n+1)/2)+1e-2;
f=256/w*sqrt(x1.^2+y1.^2)/(n-1)/sqrt(2); 
  %modification pour que l'image reste inchangée quelle que soit la valeur de n
rteta=0.11*cos(4*angle(x1+sqrt(-1)*y1))+0.89;
  %modification de cs à partir de la documentation sur vdp
  %c'est ^(-1/5) et non ^(1/5), la différence est le comportement en 0.
cs=(1+(3.23*(f*la).^(-0.6)).^5).^(-1/5); 
s1=cs.*f*al*0.9.*exp(-0.9*bl*f).*sqrt(1+0.06*exp(-0.9*bl*f));
f2=f./(ra*re*rteta);
cs=(1+(3.23*(f2*la).^(-0.6)).^5).^(-1/5); 
s2=cs.*f2*al*0.9.*exp(-0.9*bl*f2).*sqrt(1+0.06*exp(-0.9*bl*f2));
%csf=p*s1+p*(s2-s1).*(s2<s1); %cette ligne a été changée en la suivante, 
csf=p*min(s1,s2);             %simple réécriture
%sur les graphiques la csf ne doit pas descendre en dessous de 0.4 du maximum
%à basse fréquence.


%for x=1:n,
%    for y=1:n,
        % calcul de la fréquence (cycles/pixel)
%        x1=x-((n+1)/2);
%        y1=y-((n+1)/2);
%        f=sqrt(((x1)*(x1))+((y1)*(y1)))/((n-1)*sqrt(2)); 
%         % calcul de la fréquence (cycles/degre)
%        f=(n*f)/w;
%        % calcul de l  orientation en radian
%        if (x1==0)
%            if y1>0 
%                angle=pi/2;
%            else
%                angle=-pi/2;
%            end
%        else
%        angle=atan(y1/x1);
%        end
%       rteta= 0.11*cos(4*angle)+0.89;     
%   % calcul de S1(f,la,w)     
%  f2=f;
%  cs=f*f*la*la;
%  if cs==0  cs=0;  else cs=exp(-0.3*log(cs)); end
%  cs=3.23*cs;
%  if cs==0  cs=0;  else cs=exp(5*log(cs)); end
%  cs=cs+1;
%  if cs==0  cs=0;  else cs=exp(log(cs)/5); end
%  s1=cs*f2*al*0.9*exp(-0.9*bl*f2)*sqrt(1+0.06*exp(-0.9*bl*f2));
%  
%   % calcul de S1(f/r,la,w)  
%  f2=f/(ra*re*rteta);
%  cs=f2*f2*la*la;
%  if cs==0  cs=0;  else cs=exp(-0.3*log(cs)); end
%  cs=3.23*cs;
%  if cs==0  cs=0;  else cs=exp(5*log(cs)); end
%  cs=cs+1;
%  if cs==0  cs=0;  else cs=exp(log(cs)/5); end
%  s2=cs*f2*al*0.9*exp(-0.9*bl*f2)*sqrt(1+0.06*exp(-0.9*bl*f2));
%  csf(x,y)=p*s1;
%  if s2<s1 csf(x,y)=p*s2; end
%             
%    end
%end

function [image1]=FiltrageCortex(imageorig,cortex,base,n,csf)
% function [image]=FiltrageCortex(imageorig,cortex,base,n)
% imageorig: image ou matrice origine
% cortex: les 30 filtres cortex
% base: filtre bande de base.
% n: dimension de l image
% image: 31 images 1:base  2:(1,1)  3:(1,2)   n:(k,l) avec n: k=(n-1/6)+1 l=reste(n-1/6);
% Notes:les  valeurs de "image" sont réelles format double.
%
%
%
% voir aussi:
%           DecompositionCortex    Dalycortex    Dalycsf    DalyLuminance   Daly   ListeDaly   
%tentative infructueuse pour résoudre le problème pour n impaire.
%enlarge=inline('[X zeros(size(X,1),1);zeros(1,size(X,2)+1)]','X');
%reduce=inline('X(1:end-1,1:end-1)','X');
filtre=ones(n);
F=fftshift(fft2(imageorig));
F=F.*csf;
% filtrage bande de base
F3=F.*base;
%if (n-2*floor(n/2)~=0) 
%  F2=reduce(ifft2(fftshift(enlarge(F3)),'symmetric'));
%else
  try
    F2=ifft2(fftshift(F3),'symmetric');
  catch
    F2=real(ifft2(fftshift(F3)));
  end;
%end;
image1(1,:,:)=F2;

%for i=1:5,
%    for j=1:6,
%for x=1:n,
%    for y=1:n,
%  filtre(x,y)=cortex(i,j,x,y);
%    end
%end
%F3=F.*filtre;
%F2=abs(ifft2(fftshift(F3)));
%image1((i-1)*6+j+1,:,:)=F2;
%end
%end


for i=1:5
  for j=1:6
    filtre=cortex(:,:,i,j);
    F3=F.*filtre;
    try 
      F2=ifft2(fftshift(F3),'symmetric');
    catch 
      F2=real(ifft2(fftshift(F3)));
    end;
    image1((i-1)*6+j+1,:,:)=F2;
  end;
end;

