decompose_image.m

this is a very very very nice code

%Given an input image, this decomposes it into a structure image, u, and a
%texture image, v. The parameters used by the algorithm can be given by
%passing a struct for the value of param. An example of such a struct can
%be seen below in the default value of param used. To use the default
%value, simply pass a value for param that is not a struct.
function [u,v] = decompose_image(IM, param)
if sum(class(param) == 'struct') < 6
    %Lambda is the tuning parameter, mu is the step size, and num_iter is
    %the number of iterations. Note that the algorithm may degrade after
    %too many iterations. Typically you want fewer than 20 iterations.
    param = struct('lamda', .01, 'mu', .1, 'num_iter', 13);
end
if size(IM,3) > 1
    IM = rgb2gray(IM);
end
f = double(IM);

lamda = param.lamda;
mu = param.mu;
num_iter = param.num_iter;
eps = .00000001;

%Initialize u as the input image
u = f;
[fy, fx] = comp_gradient(f);
ndf = sqrt(fx.^2 + fy.^2);
ndf = ndf(:)' * ndf(:);

%Initialize g1 and g2 as defined in the report
g1 = -fx ./ (2*lamda*ndf);
g2 = -fy ./ (2*lamda*ndf);
un = zeros(size(u,1)+2,size(u,2)+2);
g1n = zeros(size(g1,1)+2,size(g1,2)+2);
g2n = zeros(size(g2,1)+2,size(g2,2)+2);

f = repl(f, 1);
for i = 1:num_iter
%     disp(i)
    u = repl(u,1);
    g1 = repl(g1,1);
    g2 = repl(g2,1);
    [M,N] = size(u);
    Ha = 1 ./ sqrt(g1.^2 + g2.^2 + eps);
    
    %Update the values for u, g1, and g2 as defined in the paper by Stanley
    %Osher and Luminita Vese, "Modeling Textures with Total Variation
    %Minimization and Oscillating Patterns in Image Processing"
    for y = 2:M-1
        for x = 2:N-1
            H = Ha(y,x);
            c1 = 1 / sqrt(  eps + ( u(y+1,x) - u(y,x))^2    + ( ( u(y,x+1)   -  u(y,x-1) ) /2)^2   );
            c2 = 1 / sqrt(  eps + ( u(y,x)   - u(y-1,x))^2  + ( ( u(y-1,x+1) -  u(y-1,x-1) ) /2)^2 );
            c3 = 1 / sqrt(  eps + ( ( u(y+1,x)   - u(y-1,x) ) /2)^2   + ( u(y,x+1) -  u(y,x) )^2   );
            c4 = 1 / sqrt(  eps + ( ( u(y+1,x-1) - u(y-1,x-1) ) /2)^2 + ( u(y,x)   -  u(y,x-1) )^2 );
            
            un(y,x) = ( f(y,x) - (g1(y+1,x) - g1(y-1,x))/2 - (g2(y,x-1) - g2(y,x-1))/2 ...
                    + ( c1*u(y+1,x) + c2*u(y-1,x) + c3*u(y,x+1) + c4*u(y,x-1) ) / (2*lamda) ) ...
                    / (1+(c1+c2+c3+c4)/(2*lamda));
            
            g1n(y,x) = ( (u(y+1,x) - u(y-1,x))/2 - (f(y+1,x) - f(y-1,x))/2 + (g1(y+1,x) + g1(y-1,x)) ...
                     + ( 2*g2(y,x) + g2(y-1,x-1) + g2(y+1,x+1) - g2(y,x-1) - g2(y-1,x) - g2(y+1,x) - g2(y,x+1) )/2 ) ...
                     * (2*lamda / (mu*H+4*lamda));
            
            g2n(y,x) = ( (u(y,x+1) - u(y,x-1))/2 - (f(y,x+1)-f(y,x-1))/2 + (g2(y,x+1)+g2(y,x-1)) ... 
                     + ( 2*g1(y,x) + g1(y-1,x-1) + g1(y+1,x+1) - g1(y,x-1) - g1(y-1,x) - g1(y+1,x) - g1(y,x+1) )/2 ) ... 
                     * (2*lamda / (mu*H+4*lamda));
        end
    end
    
    %Update u, g1, and g2
    u = un(2:M-1, 2:N-1);
    g1 = g1n(2:M-1, 2:N-1);
    g2 = g2n(2:M-1, 2:N-1);
end

%The structure image v, is defined as dg1/dx + dg2/dy
filt = [-.5, 0, .5];
v = imfilter(g1, filt) + imfilter(g2, filt');