% To apply affine diffusion equation u_t=div(A*grad u) on images
% where A is a positive definite matrix

% Accept Image from the user
disp('Enter Image')
im=uigetfile('*.jpg','Select the Image file');

% input image file and store in u1,u2,u3
u1=imread(im);
u1=double(u1);

% Accept the matrix from the user
a=input('Enter a: ');
b=input('Enter b: ');
c=input('Enter c: ');
T1=input('Enter Stop Time: ');

% Display Original Image
subplot(1,2,1), imshow(uint8(u1))
title('Original image');

% store size in m and n
[m,n,~]=size(u1);

% step size along t
dt=0.25;

for t = 0:dt:T1

% finite difference approximation for u_xx
u_xx = u1(:,[2:n n],:) - 2*u1 + u1(:,[1 1:n-1],:);

% finite difference approximation for u_xy
u_xy = ( u1([2:m,m],[2:n,n],:) + u1([1,1:m-1],[1,1:n-1],:) - ...
             u1([1,1:m-1],[2:n,n],:) - u1([2:m,m],[1,1:n-1],:) ) / 4;

% finite difference approximation for u_yy
u_yy = u1([2:m m],:,:) - 2*u1 + u1([1 1:m-1],:,:);

% Calculate RHS
div_mat_u=a*u_xx+2*b*u_xy+c*u_yy;

u1= u1 + dt*(div_mat_u);
temp=u1;

end
subplot(1,2,2), imshow(uint8(temp))
str1=sprintf('T=%d',T1);
title(str1);
clc;