Image-Processing/HeatEquationNonLinRGB.m at master · lernerjn/Image-Processing · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function [ I ] = HeatEquationNonLinRGB( I, tend, dt, strIn )
%Author: Jeremy Lerner, Stony Brook University
%This program runs the non linear heat equation on a color image,
%then saves the output image as a jpg.
%inputs:
%   I:      n x m array of the grayscale values for an image
%   tend:   end time
%   dt:     time step
%   str:    string with the name of the image, to be used in the title and
%           saved file.
%
%outputs:
%   I:      The inputted image, after being operated on by heat equation


dx = 1;
dy = 1;

[m,n,~] = size(I);

%perform all operations in double precision
I = double(I)/255;
% subplot(121);
image(I);

% Operate on the image with heat equation, using finite differences
for time=dt:dt:tend
    Iold = I;


    % main loop, using a centered finite difference to approximate the
    % second derivatives in heat equation, % dU/dt = d^2U/dx^2 + d^2U/dy^2.
    % Note, this loop is here for readability, the real work is done in
    % Matlab notation below
%     for i=2:m-1
%         for j=2:n-1
%             I(i,j) = Iold(i,j) + dt*((Iold(i+1,j) - 2*Iold(i,j) + Iold(i-1,j))*(Iold(i,j+1)-Iold(i,j-1)).^2 / (2*dy*dx^2) ...
%                 + 2*( (Iold(i+1,j+1) - Iold(i-1,j+1) + Iold(i-1,j-1) - Iold(i+1,j-1))*(Iold(i+1,j) - Iold(i-1,j))*(Iold(i,j+1)-Iold(i,j-1)))/(16*(dx^2)*(dy^2)) ...
%                 + (Iold(i,j+1) - 2*Iold(i,j) + Iold(i,j-1))*(Iold(i+1,j)-Iold(i-1,j).^2)/(2*dx*dy^2))/( ((Iold(i+1,j) - Iold(i-1,j))/(2*dx)).^2 + ((Iold(i,j+1) - Iold(i,j-1))/(2*dy)).^2);
%         end
%     end

%     for i=2:m-1
%         for j=2:n-1
%             Uxx = ( Iold(i-1,j) - 2*Iold(i,j) + Iold(i+1,j))/(dx^2);
%             Uyy = ( Iold(i,j-1) - 2*Iold(i,j) + Iold(i,j+1))/(dy^2);
%             Uxy = ( Iold(i+1,j+1) - Iold(i-1,j+1) + Iold(i-1,j-1) - Iold(i+1,j-1)) / (4*dx*dy);
%             Ux  = ( Iold(i+1,j) - Iold(i-1,j)) / (2*dx);
%             Uy  = ( Iold(i,j+1) - Iold(i,j-1)) / (2*dy);
%
%             I(i,j) = Iold(i,j) + dt*( Uxx * Uy^2 - 2*Uxy*Ux*Uy + Uyy*Ux^2)/(Ux^2+Uy^2+1e-14);
%
%         end
%     end

    %The most efficient way of implementing the for loops, using double
    %colon notation


    Uxx = ( Iold(3:m,2:n-1,:) - 2*Iold(2:m-1,2:n-1,:) + Iold(1:m-2,2:n-1,:))./(dx^2);
    Uyy = ( Iold(2:m-1,3:n,:) - 2*Iold(2:m-1,2:n-1,:) + Iold(2:m-1,1:n-2,:))./(dy^2);
    Uxy = ( Iold(3:m,3:n,:) - Iold(1:m-2,3:n,:) + Iold(1:m-2,1:n-2,:) - Iold(3:m,1:n-2,:))./(4*dx*dy);
    Uy = ( Iold(2:m-1,3:n,:) - Iold(2:m-1,1:n-2,:))./(2*dy);
    Ux = ( Iold(3:m,2:n-1,:) - Iold(1:m-2,2:n-1,:))./(2*dx);


    I(2:m-1,2:n-1,:) = Iold(2:m-1,2:n-1,:) + dt.*(Uxx .* Uy.^2 - 2*Uxy.*Ux.*Uy + Uyy.*Ux.^2)./((Ux.^2 + Uy.^2 + eps));

    I( I < 0 ) = 0;
    I(I > 1 ) = 1;


end

    %display the miage
%     figure(2);
%     image(I);
%     str = sprintf('%s at t=%f with dt=%f', strIn, tend,dt );
%     title(str);
%     h = figure(2);
%     str2 = sprintf('%s%g', strIn, tend );
%
    %save the image
    %saveas(h,str2,'jpg');

end