MatLab | function [ output_picture ] = - Zara - m684db0ab

function [ output_picture ] = color_inpaint( input_picture )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

YCrCb = rgb2ycbcr(input_picture);
clear input_picture;

Y = double(YCrCb(:,:,1));
Cb = double(YCrCb(:,:,2));
Cr = double(YCrCb(:,:,3));
clear YCrCb;

mask = Cr .* Cb;
mask = mask ~= 16384;
mask = double(mask);

mask(:,:) = 1; % for testing purpose. to be commented.
mask(100:end-101,100:end-101) = 0; % for testing purpose. to be commented.

f = laplacian(Y);

newY =  poisson_solver_function(f, Y, mask);
newCr = poisson_solver_function(f, Cr, mask);
newCb = poisson_solver_function(f, Cb, mask);

colored(:,:,1) = uint8(Y);
colored(:,:,2) = uint8(abs(newCb));
colored(:,:,3) = uint8(abs(newCr));

output_picture = ycbcr2rgb(colored);

end


%----------------------------------


function [ result ] = laplacian( picture )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

[height,width,layer_count] = size(picture);

clear layer_count;

height = uint32(height);
width = uint32(width);

j.html">j = 1:height-1;
k = 1:width-1;

gx = zeros(height,width);
gy = zeros(height,width);
gxx = zeros(height,width);
gyy = zeros(height,width);

clear height width;

gx(j.html">j,k) = (picture(j.html">j,k+1) - picture(j.html">j,k));
gy(j.html">j,k) = (picture(j.html">j+1,k) - picture(j.html">j,k));
gyy(j.html">j+1,k) = gy(j.html">j+1,k) - gy(j.html">j,k);
gxx(j.html">j,k+1) = gx(j.html">j,k+1) - gx(j.html">j,k);

clear j.html">j k gx gy;

result = gxx + gyy;

clear gxx gyy;

end


%----------------------------------


function [img_direct] = poisson_solver_function(laplacian, boundary, mask)
% function [img_direct] = poisson_solver_functionlaplacian, boundary, mask)
% Inputs; Gx and Gy -> Gradients
% Boundary Image -> Boundary image intensities
% Gx Gy and boundary image should be of same size

% Is this code a "direct solver" with dirichlet boundaries ?
% bad for non-rectangular boundaries and other constraints ?

[H,W] = size(boundary);

% boundary image contains image intensities at boundaries
boundary = boundary .* mask;

j.html">j = 100:H-101;
k = 100:W-101;
%j = 2:H-1;
%k = 2:W-1;
boundary_laplacian = zeros(H,W);
boundary_laplacian(j.html">j,k) = -4*boundary(j.html">j,k) + boundary(j.html">j,k+1) + boundary(j.html">j,k-1) + boundary(j-1,k) + boundary(j.html">j+1,k);
clear j.html">j k

% subtract boundary points contribution
f1 = laplacian - reshape(boundary_laplacian, H, W);
clear laplacian boundary_laplacian

% DST Sine Transform algo starts here
f2 = f1(2:end-1,2:end-1);
clear f1

%compute sine transform
tt = dst(f2);
f2sin = dst(tt')';
clear f2

%compute Eigen Values
[x,y] = meshgrid(1:W-2,1:H-2);
denom = (2*cos(pi*x/(W-1))-2) + (2*cos(pi*y/(H-1)) - 2);

%divide
f3 = f2sin./denom;
clear f2sin x y

%compute Inverse Sine Transform
tt = idst(f3);
clear f3;
img_tt = idst(tt')';
clear tt

% put solution in inner points; outer points obtained from boundary image
vertzero = zeros(1, W-2);
horzzero = zeros(H,1);

plop = vertcat(vertzero, img_tt, vertzero);
plop = horzcat(horzzero, plop, horzzero);

img_direct = (boundary .* mask) + (plop .* (1 - mask));

end


%----------------------------------


function b=dst(a,n)
%DST Discrete sine transform (Used in Poisson reconstruction)
% Y = DST(X) returns the discrete sine transform of X.
% The vector Y is the same size as X and contains the
% discrete sine transform coefficients.
% Y = DST(X,N) pads or truncates the vector X to length N
% before transforming.
% If X is a matrix, the DST operation is applied to each
% column. This transform can be inverted using IDST.

error(nargchk(1,2,nargin));

if min(size(a))==1
    if size(a,2)>1
        do_trans = 1;
    else
        do_trans = 0;
    end
    a = a(:);
else
    do_trans = 0;
end

if nargin==1, n = size(a,1); end

m = size(a,2);

% Pad or truncate a if necessary
if size(a,1)<n,
    aa = zeros(n,m);
    aa(1:size(a,1),:) = a;
else
    aa = a(1:n,:);
end

y=zeros(2*(n+1),m);
y(2:n+1,:)=aa;
y(n+3:2*(n+1),:)=-flipud(aa);
yy=fft(y);
b=yy(2:n+1,:)/(-2*sqrt(-1));

if isreal(a), b = real(b); end
if do_trans, b = b.'; end


%----------------------------------

function b=idst(a,n)
%IDST Inverse discrete sine transform (Used in Poisson reconstruction)
%
% X = IDST(Y) inverts the DST transform, returning the
% original vector if Y was obtained using Y = DST(X).
% X = IDST(Y,N) pads or truncates the vector Y to length N
% before transforming.
% If Y is a matrix, the IDST operation is applied to
% each column.

if nargin==1
    if min(size(a))==1
        n=length(a);
    else
        n=size(a,1);
    end
end

nn=n+1;
b=2/nn*dst(a,n);