clear;
close all;

%% visualization of rotational covariance

[x,y] = meshgrid(-10:.1:10);
[xc,yc] = meshgrid(-10:10);

map = [[linspace(0,1,128)';ones(128,1)],[linspace(0,1,128)';linspace(1,0,128)'],[ones(128,1);linspace(1,0,128)']];

%% one-dimensional function

%f = @(x,y) .1*(x-3).^2.*y-10*y+y.^2;
%gx = @(x,y) .2*(x-3).*y;
%gy = @(x,y) .1*(x-3).^2-10+2*y;

% with the following the rotation is easier to see, since all gradients are
% poiting in the same directions
f = @(x,y) x.^2;
gx = @(x,y) 2*x;
gy = @(x,y) zeros(size(x));


subplot(121)
imagesc(-10:10,-10:10,f(x,y))
hold on
quiver(xc,yc,gx(xc,yc),gy(xc,yc),'k')
hold off
axis equal, axis([-10 10 -10 10])
%title('f(x,y) = 0.1(x-3)^2y-10y+y^2')
title('f(u) = f(x,y) = x^2')

colormap(map)

phi = -20; % angle in degrees
c = cosd(phi); s = sind(phi);

subplot(122)
imagesc(-10:10,-10:10,f(c*x-s*y,s*x+c*y))
hold on
quiver(xc,yc,c*gx(c*xc-s*yc,s*xc+c*yc)+s*gy(c*xc-s*yc,s*xc+c*yc),-s*gx(c*xc-s*yc,s*xc+c*yc)+c*gy(c*xc-s*yc,s*xc+c*yc),'w')
hold off
axis equal, axis([-10 10 -10 10])
title(['f rotated by ',num2str(phi),' degrees: f(Ru)'])

pause

% overlay gradient of original f
imagesc(-10:10,-10:10,f(c*x-s*y,s*x+c*y))
hold on
quiver(xc,yc,c*gx(c*xc-s*yc,s*xc+c*yc)+s*gy(c*xc-s*yc,s*xc+c*yc),-s*gx(c*xc-s*yc,s*xc+c*yc)+c*gy(c*xc-s*yc,s*xc+c*yc),'w')
quiver(xc,yc,gx(xc,yc),gy(xc,yc),'k')
hold off
axis equal, axis([-10 10 -10 10])
title(['gradient of original f overlayed: (grad f)(u)'])

pause

% rotate domain for gradient of original f
imagesc(-10:10,-10:10,f(c*x-s*y,s*x+c*y))
hold on
quiver(xc,yc,c*gx(c*xc-s*yc,s*xc+c*yc)+s*gy(c*xc-s*yc,s*xc+c*yc),-s*gx(c*xc-s*yc,s*xc+c*yc)+c*gy(c*xc-s*yc,s*xc+c*yc),'w')
quiver(xc,yc,gx(c*xc-s*yc,s*xc+c*yc),gy(c*xc-s*yc,s*xc+c*yc),'k')
hold off
axis equal, axis([-10 10 -10 10])
title(['gradient of original f with rotated domain: (grad f)(Ru)'])

pause

% rotate gradient vector field of rotated domain with R'
imagesc(-10:10,-10:10,f(c*x-s*y,s*x+c*y))
hold on
quiver(xc,yc,c*gx(c*xc-s*yc,s*xc+c*yc)+s*gy(c*xc-s*yc,s*xc+c*yc),-s*gx(c*xc-s*yc,s*xc+c*yc)+c*gy(c*xc-s*yc,s*xc+c*yc),'w')
quiver(xc,yc,c*gx(c*xc-s*yc,s*xc+c*yc)+s*gy(c*xc-s*yc,s*xc+c*yc),-s*gx(c*xc-s*yc,s*xc+c*yc)+c*gy(c*xc-s*yc,s*xc+c*yc),'k')
hold off
axis equal, axis([-10 10 -10 10])
title(["rotated gradient of original f with domain rotated: R' (grad f)(Ru)"])

pause

%%

clear f gx gy c s
close all

%% white box on black background

phi = 20; % angle in degrees
c = cosd(phi); s = sind(phi);

L_filt = [0 1 0;1 -4 1;0 1 0];

f = @(x,y) x>-2 & x<4 & y>0 & y<5;

% smoothen f
f_s = conv2(f(x,y),fspecial('gaussian',[5 5],2),'same');

subplot(231)
imagesc(-10:10,-10:10,f_s,[-1 1])
axis image
title('f: box with dimensions 6x5')

colormap(map)

pause

Lf = conv2(f_s,L_filt,'same');

subplot(232)
imagesc(-10:10,-10:10,Lf,[-.35 .35])
axis image
title('\Delta f')

pause

% smoothen rotated f
f_s_rot = conv2(f(c*x-s*y,s*x+c*y),fspecial('gaussian',[5 5],2),'same');

subplot(234)
imagesc(-10:10,-10:10,f_s_rot,[-1 1])
axis image
title(['f \cdot R, \phi=',num2str(phi),'deg'])

Lf_rot = conv2(f_s_rot,L_filt,'same');

pause

subplot(235)
imagesc(-10:10,-10:10,Lf_rot,[-.35,.35])
axis image
title('\Delta (f \cdot R)')

pause

subplot(133)
imagesc(-10:10,-10:10,interp2(x,y,Lf,c*x-s*y,s*x+c*y,'spline',0),[-.35,.35])
axis image
title('(\Delta f) \cdot R')

pause

%%

clear f f_s f_s_rot c s L_filt Lf Lf_rot