%Ejemplo de cómo generar figuras periódicas, evitando el aliasing

%Generación de redes sinusoidales acromáticas, para medida de la CSF

%Cargamos el sistema de referencia y la información sobre el monitor
[T_l,Yw,Msx]=loadsys;
[tm,a,g]=loadmon(Msx);

%PROCEDIMIENTO ERRONEO: NO SE TIENE EN CUENTA LA POSIBILIDAD DE ALIASING
%La imagen será una red sinusoidal en una ventana cuadrado de 256x256 píxels
%Definimos el dominio espacial del test, en grados subtendidos desde el observador
tamano=6.7;
dist=50;
ang=2*atan(tamano/dist/2)*180/pi;
x=linspace(-ang/2,ang/2,256);
[x,y]=meshgrid(x);



%Elegimos la frecuencia espacial en ciclos/grado
fx=32;
fy=0;

%Definimos el estímulo en el espacio ATD de Boynton
%La variación sinusoidal será alrededor de un acromático de 50 cd/m2
%Coordenadas cromáticas del blanco del monitor
Tw=val2tri([1 1 1],Yw,tm,a,g);
tYw=tri2coor(Tw,Yw);
tYw=[tYw(:,1:2) 50];
Tw=coor2tri(tYw,Yw);
%Escribimos la variación en el espacio ATD de Boynton
ATDw=xyz2atd(Tw,6);
%Definimos el contraste de la red
%C=(Amax-Amin)/(Amax+Amin)
C=1;
A=ATDw(1)*(1+C*sin(2*pi*(fx*x+fy*y)));
T=ATDw(2)*ones(256,256);
D=ATDw(3)*ones(256,256);

%Definimos una imagen de color verdadero con estos valores
ATD(:,:,1)=A;
ATD(:,:,2)=T;
ATD(:,:,3)=D;

%Esta matriz contiene 256*256 colores que tendríamos que transformar
% a valores triestímulo y luego a niveles digitales. Mucho tiempo de cálculo
%Reducimos el número de colores a 256 (máximo representable simultaneamente
%con nuestras tarjetas)
[im,palATD]=true2pal(ATD,256);

%Transformamos la paleta resultante a T y a n
palT=atd2xyz(palATD,1,6);
[n,sat,Tn]=tri2val(palT,Yw,tm,a,g,8);

%Representamos el resultado
figure;imshow(im,n)

%Aparece aliasing! Lo hubieramos podido detectar antes, representando A vs. x
figure;plot(x(1,:),A(1,:));xlabel('Posicion (g)');ylabel('A');

%¿Cómo evitamos este problema? ¡Calculando la frecuencia de Niquist!
Npixels=256;
fmuestreo=Npixels/ang;
fNiquist=fmuestreo/2;

%Para imágenes no cuadradas habría que trabajar con la fNiquist en x y en y
%fgenerable=fNiquist/(2*k+1), k entero
%Ejemplo, k=4
k=4;
fx=fNiquist/(2*k+1);
fy=0;

%Generamos la nueva red
C=1;
A=ATDw(1)*(1+C*sin(2*pi*(fx*x+fy*y)));
T=ATDw(2)*ones(256,256);
D=ATDw(3)*ones(256,256);

%Definimos una imagen de color verdadero con estos valores
ATD(:,:,1)=A;
ATD(:,:,2)=T;
ATD(:,:,3)=D;

%Esta matriz contiene 256*256 colores que tendríamos que transformar
% a valores triestímulo y luego a niveles digitales. Mucho tiempo de cálculo
%Reducimos el número de colores a 256 (máximo representable simultaneamente
%con nuestras tarjetas)
[im,palATD]=true2pal(ATD,256);

%Transformamos la paleta resultante a T y a n
palT=atd2xyz(palATD,1,6);
[n,sat,Tn]=tri2val(palT,Yw,tm,a,g,8);

%Representamos el resultado
figure;imshow(im,n)

%Aparece aliasing! Lo hubieramos podido detectar antes, representando A vs. x
figure;plot(x(1,:),A(1,:));xlabel('Posicion (g)');ylabel('A');

%Si queremos generar una frecuencia mayor a la fNiquist que hemos calculado,
%o aumentamos el número de pixels por grado de la imagen (difícil, depende del
%monitor) o disminuímos el ángulo subtendido por esta (alejando al observador)


%Una mirada más detallada al aliasing
%Representemos un conjunto de frecuencias con aliasing

k=4;
fy=0;

%Generamos la nueva red
C=1;
for i=1:2:10
 fx=fNiquist*(i-1)+fNiquist/(2*k+1)
 A=ATDw(1)*(1+C*sin(2*pi*(fx*x+fy*y)));
 T=ATDw(2)*ones(256,256);
 D=ATDw(3)*ones(256,256);
 ATD(:,:,1)=A;
 ATD(:,:,2)=T;
 ATD(:,:,3)=D;
 [im,palATD]=true2pal(ATD,256);
 palT=atd2xyz(palATD,1,6);
 [n,sat,Tn]=tri2val(palT,Yw,tm,a,g,8);

 %Representamos el resultado
 figure;imshow(im,n);title(['f=' num2str(fx)]);

 %Aparece aliasing! Lo hubieramos podido detectar antes, representando A vs. x
 figure;plot(x(1,:),A(1,:));xlabel('Posicion (g)');ylabel('A');title(['f=' num2str(fx)]);
end