function [ res ] = embed( dest, patch, subx, suby, blendwidth )
%embed Blends the patch into the destination image

res = dest;

if (subx < 1)
    subx = 1;
end

if (suby < 1)
    suby = 1;
end

[rows, cols, channels] = size(dest);
[prows, pcols, pchannels] = size(patch);

%shortestDist = round(min([prows pcols]) / 2);

% let's blend better!
mask = fspecial('gaussian', [64 64], 15);
%mask = mask / mask(32,round(blendwidth * 31) + 1);
mask = mask / max(mask(:));
if (blendwidth == 0)
    mask = ones(64, 64);
else
    mask = mask / blendwidth;
end
mask(find(mask > 1)) = 1;
mask = imresize(mask, [prows pcols], 'bilinear');
alphaMask(:,:,1) = mask;
alphaMask(:,:,2) = mask;
alphaMask(:,:,3) = mask;

clear mask;

x = round(subx - (pcols/2) + 1);
y = round(suby - (prows/2) + 1);

if (x < 1)
    x = 1;
elseif (x > cols + 1 - pcols)
    x = cols + 1 - pcols;
end

if (y < 1)
    y = 1;
elseif (y > rows + 1 - prows)
    y = rows + 1 - prows;
end

%background = res(y:y+prows-1, x:x+pcols-1, :);

%background = im2double(background);
%patch = im2double(patch);

%blended = alphaMask .* patch + (1 - alphaMask) .* background;

%res(y:y+prows-1, x:x+pcols-1, :) = blended;


res(y:y+prows-1, x:x+pcols-1, :) = alphaMask .* patch + (1 - alphaMask) .* res(y:y+prows-1, x:x+pcols-1, :);
%res(y:y+prows-1, x:x+pcols-1, :) = patch;

%{
for x = 1 : 1 : pcols - 1
    for y = 1 : 1 : prows - 1
        
        p = min([x, y, pcols - x, prows - y]);
        if (p > blendwidth)
            alpha = 0;
        else
            alpha = (blendwidth - (p - 1)) / blendwidth;
        end
        
        
        %{
        xDist = x - round(pcols / 2);
        yDist = y - round(prows / 2);
        distToCenter = sqrt(xDist * xDist + yDist * yDist);
        distFromEdge = shortestDist - distToCenter        
        %}

        tx = round(subx - (pcols/2) + x);
        ty = round(suby - (prows/2) + y);
        %tx = subx + x;
        %ty = suby + y;
        %{
        res(ty, tx, 1) = alpha * patch(y, x, 1) + (1-alpha) * dest(ty, tx, 1);
        res(ty, tx, 2) = alpha * patch(y, x, 2) + (1-alpha) * dest(ty, tx, 2);
        res(ty, tx, 3) = alpha * patch(y, x, 3) + (1-alpha) * dest(ty, tx, 3);
        %}
        if (tx > 0 && ty > 0 && tx < cols && ty < rows)
            res(ty, tx, :) = (1-alpha) * patch(y, x, :) + alpha * dest(ty, tx, :);
        end
        
    end
end
%}


end