% Reconstrucción IP con SR, para r=2s-1. Sin filtros, con cualquier grado
% para el polinomio interpolador. Se eliminan todos los coeficicentes de
% escala.
%       preliipvpsrgenultim(xx,ff,grau,nivells)
%       Entradas:
%           xx      = nodos de la funcion inicial
%           ff      = f(xx)
%           grau    = grado del polinomio interpolador (r)
%           nivells = niveles de multirresolucion 
%       Salidas:
%           Plot de la reconstruccion(-), funcion en el nivel mas grosero(o), funcion original con mas nodos (--), diferencia original y aproximacion(:) 
%
% Autor: José Jaime Noguera. +INFO: www.uv.es/anims/noguera

function preliipvpsrgenultim(xx,ff,grau,nivells)

nivells=nivells+1;
[xn,xm]=size(xx);
n=log(xm-1)/log(2);
f=zeros(nivells,xm);
x=zeros(nivells,xm);
f(nivells,1:xm)=ff(1:xm);
x(nivells,1:xm)=xx(1:xm);


%Discretizar por valores puntuales
for k=nivells:-1:2
   for i=1:1+2^(n-nivells+k-1)
      f(k-1,i)=f(k,2*i-1);
      x(k-1,i)=x(k,2*i-1);
   end 
end 
fr=zeros(nivells,xm); %fr significa la reconstruccion
fr(1,1:2^(n-nivells+1)+1)=f(1,1:2^(n-nivells+1)+1);


%Reconstruccion interpolatoria+eno+sr
for k=2:nivells   
     ultim=1+2^(n-nivells+k-1);
        %%%%%%a%%%%%%%%%% asignacion eno para vp
        for i=2:ultim
            t=0;
            for l=i-grau+1:i
                t=t+1; %para evitar indices negativos en dif(t)
                if(l-1>=1 && l+grau-1<=ultim)
                   dd=difdiv(x(k-1,l-1:l+grau-1),fr(k-1,l-1:l+grau-1),grau+1);
                   dif(t)=abs(dd(grau+1,grau+1));
                else
                    dif(t)=1e100;
                end    
            end
        [minin,pos]=min(dif(1:grau));
        is(i)=i+pos-grau;
        end
       %%%%%%%%%%%%%%%% fin-eno
        
        for i=grau+2:ultim-grau
            if(is(i-1)+grau-1 < is(i+1)-1) %CELDA SOSPECHOSA
                   Gl=intseg(x(k-1,1:2^(n-(nivells-k+1))),fr(k-1,1:2^(n-(nivells-k+1))),is,i+1,grau,x(k-1,i-1))-intseg(x(k-1,1:2^(n-(nivells-k+1))),fr(k-1,1:2^(n-(nivells-k+1))),is,i-1,grau,x(k-1,i-1)); %G_i(x^{k-1}_{i-1})
                   Gc=intseg(x(k-1,1:2^(n-(nivells-k+1))),fr(k-1,1:2^(n-(nivells-k+1))),is,i+1,grau,x(k,2*i-2))-intseg(x(k-1,1:2^(n-(nivells-k+1))),fr(k-1,1:2^(n-(nivells-k+1))),is,i-1,grau,x(k,2*i-2));%G_i(x^k_{2i-2})
                   Gr=intseg(x(k-1,1:2^(n-(nivells-k+1))),fr(k-1,1:2^(n-(nivells-k+1))),is,i+1,grau,x(k-1,i))-intseg(x(k-1,1:2^(n-(nivells-k+1))),fr(k-1,1:2^(n-(nivells-k+1))),is,i-1,grau,x(k-1,i)); %G_i(x^{k-1}_i)
                   if(Gl*Gr<0) % CELDA SINGULAR
                        if(Gc*Gr<0) %estamos a la izquierda del salto
                           fr(k,2*i-2)=interpola(x(k-1,i-grau-1:i-1),fr(k-1,i-grau-1:i-1),x(k,2*i-2));
                        else %estamos a la derecha del salto
                           fr(k,2*i-2)=interpola(x(k-1,i:i+grau),fr(k-1,i:i+grau),x(k,2*i-2));
                        end
                   else %CELDA NO SINGULAR
                        %hagamos eno
                           fr(k,2*i-2)=interpola(x(k-1,is(i)-1:is(i)+grau-1),fr(k-1,is(i)-1:is(i)+grau-1),x(k,2*i-2));  
                   end
            else %CELDA NO SOPECHOSA
                 %hagamos eno
                    fr(k,2*i-2)=interpola(x(k-1,is(i)-1:is(i)+grau-1),fr(k-1,is(i)-1:is(i)+grau-1),x(k,2*i-2));
            end     
            fr(k,2*i-1)=fr(k-1,i);
        end
  
        %FRONTERAS
         fr(k,1)=fr(k-1,1);
          for i=2:grau+1
              fr(k,2*i-2)=interpola(x(k-1,is(i)-1:is(i)+grau-1),fr(k-1,is(i)-1:is(i)+grau-1),x(k,2*i-2));  
              fr(k,2*i-1)=fr(k-1,i);
          end
          for i=ultim-grau+1:ultim
              fr(k,2*i-2)=interpola(x(k-1,is(i)-1:is(i)+grau-1),fr(k-1,is(i)-1:is(i)+grau-1),x(k,2*i-2));  
              fr(k,2*i-1)=fr(k-1,i);
          end
        %fin-fronteras
end   


[xor,yor]=funtreballtesi(257,0); 

hold on
plot(xor, yor,'--k','LineWidth',2)                                                        %original discretizada con 257 pts 
plot(x(1,1:1+2^(n-(nivells-2+1))),f(1,1:1+2^(n-(nivells-2+1))),'ko')                      %original en el nivel mas grosero
plot(xx,(ff-fr(nivells,1:xm))-(0.1)*ones(1,xm),':k','LineWidth',3)                        %original-reconstruccion
plot(xx,fr(nivells,1:xm),'-k','LineWidth',3)                                              %reconstruccion
axis([0 1 -0.2 1.1 ]);
set(gca,'Xtick',[0 .5 1], 'FontSize',15)
set(gca,'Ytick',[0 .5 1], 'FontSize',15)
box
hold off