function Irec = muestreodimdos(func,d);


close all;
    I=imread(func);
    I=double(I)./255;
    s=size(I);
    figure;
    imshow(I);

    m=s(1);
    n=s(2);
    n=min([n,m]);
    I=I(1:1:n,1:1:n);
    imshow(I);
    str(1) = {'Imagen original'};
    text(30,30,str(1))
    pause;

    % Definimos el ancho de banda
    r=floor(n/(2*d))-1;
    b=r;

    % Calculamos la mejor aproximación 
    % en el espacio de las señales del ancho 
    % de banda dado
    k=zeros(n,n);
    k(1:r,1:r)=1;
    k((n-r):n, (n-r):n)=1;
    k((n-r):n, 1:r)=1;
    k(1:r,(n-r):n)=1;
    
    P=fft2(I);
    %figure;
    %imshow(P);pause;
    
    P=P.*k;
    %figure;
    %imshow(P);pause;

    Ibanda=ifft2(P);
    figure;
    imshow(Ibanda);
    str1(1) = {'Imagen, tras ser reducido su ancho de banda'};
    text(30,30,str1(1));pause;
    figure;imshow(1-(I-Ibanda));
    str1(2) = {'Diferencia con la Imagen Original'};
    text(30,30,str1(2))
    pause;

    %Tomamos muestras:
    Y=ones(n,n);
    Y(1:d:n, 1:d:n) = Ibanda(1:d:n, 1:d:n);
    
    figure; 
    imshow(Y);
    str1(3) = {'Muestras'};
    text(30,30,str1(3))
    pause;


    % Por fín aplicamos el algoritmo a 
    % la imagen muestreada. 
    r=(n/d)-1;
    AA=ones(r+1,r+1);
    AA = Y(d*(0:r)+1, d*(0:r)+1);

 sincvector=discretsincvector(b,n);
 Irec=Y;
 
 
  V1 = ones(r+1,r+1);
 V2 = ones(r+1,r+1);
 v = ones(1, r+1);
 aux = n - d*(0:r);
 for s=0:(n-1);
     v = sincvector(aux+s);
     V1(s+1, :) = v;
 end
 V2= V1';
 Irec = d.^2 * V1 * AA * V2;

       
figure;imshow(Irec);
str1(4) = {'Imagen recuperada'};
    text(30,30,str1(4));
%figure;imshow(fft2(Irec));


function SM = discretsincvector(M,N)
% function SM = discretsincvector(M,N)
% Calcula el analogo discreto de la funcion seno cardinal 
% para trabajar con señales de banda limitada (ancho de banda M) y 
% tamaño N dados. Calcula dicha funcion entre -(N-1) y N-1, que es donde es
% util
%
%  Parametros:   
%   M: ancho de banda 
%   N: tamaño de la señal
%   SM: vector salida
%
%  Autores: José M. Almira (2009).
%

  % Vector de abcisas
  t = (-N+1):(N-1);
  
  % Reservamos memoria: la señal SM (salida) será un vector 1xN (todo a
  % cero, inicialmente)
  SM = repmat(double(0), size(t));
  
  % Calculamos la función muestreo para valores distintos de cero
  % Ojo al direccionamiento con vectores
  
  SM( t>0 ) =  sin(pi*(2*M+1)*t(t>0)/N)./(N*sin(pi*t(t>0)/N));
  SM( t<0 ) =  sin(pi*(2*M+1)*(t(t<0)+N)/N)./(N*sin(pi*(t(t<0)+N)/N));
     
  % Asignamos el valor (en el límite) de la funcion para t = 0 
  
  SM( t==0 ) = (2*M+1)/N;