f_noisy = double(imread('noisy_input.png')); [nx, ny, c] = size(f_noisy); f = f_noisy(:); n = nx*ny*c; %% construct the color gradient operator e = ones(ny, 1); Dx_tilde = spdiags([e -e], -1:0, ny, ny); Dx_tilde(:, end) = 0; e = ones(nx, 1); Dy_tilde = spdiags([-e e], 0:1, nx, nx); Dy_tilde(end, end) = 0; Dx = kron(Dx_tilde', speye(nx)); Dy = kron(speye(ny), Dy_tilde); D = cat(1, kron(speye(3), Dx), kron(speye(3), Dy)); norm_D = normest(D); % smoothing of TV epsilon = 0.1; % weight of the dataterm lambda = 0.05; % Lipschitz constant of the gradient L = lambda + norm_D^2/epsilon; % compute step size tau = 2/(lambda+L); % number of iterations max_iter = 1500; % allocate memory for current iterate u_k and energy u_k = zeros(n, 1); energy = zeros(max_iter, 1); %% do gradient descent iterations for k = 1:max_iter ['Current iteration: ', num2str(k)] % compute the gradient of u_k at the k-th iteration Du = D*u_k; % compute the TV TV = sqrt(Du.^2 + epsilon^2); % compute the gradient of the whole energy Grad = lambda*(u_k-f) + D'*(Du./TV); % update u_k = u_k - tau*Grad; % compute the energy E(u_{k+1}) of the k+1-th iterate u_k energy(k) = 0.5 * lambda * norm(u_k - f)^2 + sum(TV); end %% plot result figure; plot(1:max_iter, energy); figure; imshow(uint8(reshape(u_k, nx, ny, 3)));