%% minimisation d'une fonction a une variable 
x=zeros(1,100);
rho=0.9

% F(x)=(x-1)^2.  grad F (x)= 2(x-1)

for k=2:100
    x(k)=x(k-1)-rho*2*(x(k-1)-1);
end
% affichage 
plot(log(abs(x-1)));

%% voir le gradient

im=imread('peppers.png');
im=double(sum(im,3)); 
[gx,gy]=grad_im(im);
imshow(abs(gx),[]);
n=gx.^2+gy.^2;
n=n.^0.5;
figure; imshow(n,[])

%% Zoom d'une image 
im=imread('peppers.png');
im=sum(double(im),3)/3; %now a gray image
[ty,tx]=size(im);
t2x=2*tx-1; %taille de l'image zoomee
t2y=2*ty-1;
imzoom=zeros(t2y,t2x);
imzoom(1:2:end,1:2:end)=im;
[XX,YY]=meshgrid(0:t2x-1,0:t2y-1);
mask=(mod(XX,2)+mod(YY,2))>0; % on ne modifie que les pixels dont 
                              %au moins une coordonnee est impaire
rho=0.1;
for k=1:10
    grad=zeros(size(imzoom));
    % calcul de la contribution au gradient de la derivee en x
    grad(:,2:end)=0+2*(imzoom(:,2:end)-imzoom(:,1:end-1));
    grad(:,1:end-1)=grad(:,1:end-1)+2*(-imzoom(:,2:end)+imzoom(:,1:end-1));
    % calcul de la contribution au gradient de la derivee en 
    grad(2:end,:)=grad(2:end,:)+2*(imzoom(2:end,:)-imzoom(1:end-1,:));
    grad(1:end-1,:)=grad(1:end-1,:)+2*(-imzoom(2:end,:)+imzoom(1:end-1,:));
    
    imzoom(mask)=imzoom(mask)-rho*grad(mask);
end
imshow(imzoom,[])
%% Reconstruction a partir d'un gradient tres degrade
im=imread('peppers.png');
im=sum(double(im),3)/3; %now a gray image
imshow(im,[])

decim=2; 
sg=2;

[ty,tx]=size(im);
[gx,gy]=grad_im(im);
[XX,YY]=meshgrid(0:tx-1,0:ty-1);
mask=(mod(XX,decim)+mod(YY,decim))>0;
gx(mask)=0;
gy(mask)=0;


n=gx.^2+gy.^2;
n=n.^(0.5);
mask=n<sg;

a=atan2(gy,gx)/(2*pi);
a=round(a*8)/8*(2*pi); %arrondi a pi/4 pres
gx=cos(a); 
gy=sin(a);

gx(mask)=0; % effacer les petits gradients
gy(mask)=0;

%gardons huit niveaux de gradient

% Mn=max(n(:));
% n=(n./Mn).^0.5;
% n=(round(n*7)/7).^2;
% gx=gx.*n;
% gy=gy.*n;

%% affichage du champ
[XX,YY]=meshgrid(1:decim:tx,1:decim:ty);
close all
imshow(im,[]);
hold on;
quiver(XX,YY,gx(1:decim:end,1:decim:end),gy(1:decim:end,1:decim:end));
hold off
%%
% solution avec Fourier
gx=gx(1:decim:end,1:decim:end);
gy=gy(1:decim:end,1:decim:end);
[ty,tx]=size(gx);
[fx,fy]=meshgrid([0:tx/2 -tx/2+1:-1], [0:ty/2 -ty/2+1:-1]);
fx=fx/tx;
fy=fy/ty;
Kx=zeros(ty,tx); Kx(1,1)=1; Kx(1,2)=-1;
fKx=fft2(Kx);
Ky=zeros(ty,tx); Ky(1,1)=1; Ky(2,1)=-1;
fKy=fft2(Ky);
fgxo=fft2(gx);
fgyo=fft2(gy);

denom=abs(fKx).^2+abs(fKy).^2;
denom(1,1)=1; %eviter la division par zero

numer=fgxo.*conj(fKx)+fgyo.*conj(fKy);
fimappr=numer./denom;
imappr=real(ifft2(fimappr));
figure; imshow(imappr,[]);
