function data = col2im_overlap(rs, size1, size2, overlapPixel, maps, method)
% method represents different combination method from pieces
% 0: original, replaced. 1: Graph-cut based 3: replacement
% 2: Pixels averaging in overlapped region
% 4: only center subpatch will be extracted and pasted.
op = overlapPixel;

s1 = size1;
data = single(zeros(size2));
rs = single(rs);
ct = 1;
if (~exist('maps', 'var')), maps=[]; end
if (~exist('method', 'var')), method = 0; end
% h = fspecial('gaussian', [op s1(1)], s1(1)/2);
% h = h / max(h(:));
for i=1:s1(1)-op:size2(1)-s1(1)+1
    for j=1:s1(2)-op:size2(2)-s1(2)+1
        rrow = i:i+s1(1)-1;
        rcol = j:j+s1(2)-1;
        if ~isempty(maps)
            varify = mean2(maps(rrow, rcol));
        end
        if (~isempty(maps) && varify>0.2) || isempty(maps)
            v = rs(:, ct);
            v = reshape(v, s1);
            bk=data(rrow,rcol);
            if (method==0)
                idx1=bk~=0 & v~=0;
                bk(~idx1) = bk(~idx1) + v(~idx1);
                if (sum(idx1(:))>0)
                    % find the overlapping region
                    
                    bk(idx1) = (v(idx1) + bk(idx1))*0.5;
                    
                end
            elseif (method == 3)
                idx1 = bk==0;
                bk(idx1) = v(idx1);
            elseif (method == 1) % using grapg-cut method
                
                idx1= bk == 0;
                bk(~idx1) = bk(~idx1) + v(~idx1);
                if (op<=1)
                    bk(idx1) = (v(idx1) + bk(idx1))*0.5;
                else
                    if (sum(idx1(:))>0)
                    r1 = 1:op;
                    c1 = 1:op;
                    r2 = op+1:s1(1);
                    c2 = op+1:s1(2);
                    
                    if (i==1)
                        diff1 = abs(v(r1,:) - bk(r1,:));
                        for k=1:s1(2)
                            [~,mv] = max(diff1(:, k));
                            bk(mv:s1(1), :) = v(1:length(mv:s1(1)), :);
                        end
                    else
                        diff1 = abs(v(:,c1) - bk(:,c1));
                        for k=1:s1(1)
                            [~,mv] = max(diff1(k,:));
                            bk(:, mv:s1(2)) = v(:, 1:length(mv:s1(2)));
                        end
                        diff1 = abs(v(r1,:) - bk(r1,:));
                        for k=1:s1(2)
                            [~,mv] = max(diff1(:, k));
                            bk(mv:s1(1), :) = v(1:length(mv:s1(1)), :);
                        end
                    end
                    data(rrow,rcol) = bk;
                    end
                end
            elseif (method==2)
                bk = bk + v;
            end
            data(rrow,rcol) = bk;
            
            ct = ct + 1;
        end
    end
end
if (method==2)
    h = fspecial('gaussian', op, op*2);
    h = h./max(h(:)) * 0.5;
    start1 = 1;
    for i=1:s1(1)-op:size2(1)-s1(1)+1
        for j=1:s1(2)-op:size2(2)-s1(2)+1
            rrow = i:i+s1(1)-1;
            rcol = j:j+s1(2)-1;
            if ~isempty(maps)
                varify = mean2(maps(rrow, rcol));
            end
            if (~isempty(maps) && varify>0.2) || isempty(maps)
                bk=data(rrow,rcol);
                idx1=bk~=0 & v~=0;
                bk(~idx1) = bk(~idx1) + v(~idx1);
                if (sum(idx1(:))>0)
                    r1 = 1:op;
                    c1 = 1:op;
                    r2 = op+1:s1(1);
                    c2 = op+1:s1(2);
                    
                    if (i==1)
                        bk(r1,c2) = bk(r1,c2) .* h;
                    else
                        bk(r1,c2) = bk(r1,c2) .* h;
                        
                        bk(r2,c1) = bk(r2,c1) .* h;
                    
                    end
                    data(rrow,rcol) = bk;
                end
                
            end

        end
    end
end
ct = 1;
s2 = s1 ./ 4;

if (method==4) % Using the center block to fill the image
    for i=1:s1(1)-op:size2(1)-s1(1)+1
        for j=1:s1(2)-op:size2(2)-s1(2)+1
            rrow = i:i+s1(1)-1;
            rcol = j:j+s1(2)-1;
            if ~isempty(maps)
                varify = mean2(maps(rrow, rcol));
            end
            if (~isempty(maps) && varify>0.2) || isempty(maps)
                v = rs(:, ct);
                v = reshape(v, s1);
                % only center will be extracted
                v = v(1+s2(1):end-s2(1), 1+s2(2):end-s2(2));
                data(i+s2(1):i+s1(1)-s2(1)-1, j+s2(2):j+s1(2)-s2(2)-1) = v;
                ct = ct + 1;
                
            end

        end
    end
end

data = single(data);