Poisson Image Editing 泊松融合 matlab代码完整

做相关的实验，发现竟然没找到的能直接用的matlab代码，于是 google一下。

原理我觉得这个就讲的挺好的：https://blog.csdn.net/hjimce/article/details/45716603

原文paper： https://www.cs.virginia.edu/~connelly/class/2014/comp_photo/proj2/poisson.pdf

matlab代码出处：https://www.mathworks.com/matlabcentral/fileexchange/62287-poisson-image-editing

 

Demo.m  > PIE_Gui > PIE

%Demo
% T3 Demo code
%
source='source6.jpg';
target='target6.jpg';
result='result6.jpg';
I1 = imread(source);            % SOURCE IMAGE
I2 = imread(target);        % DESTINATION IMAGE
PIE_Gui(I1,I2,result,1,0);

 

%%PIE_Gui
% Author: Michael S. Brown, York University
%
% Simple Matlab GUI for selection a Region of Interest from Image1, and
% pasting this into Image2.
%
% It calls a function called: PossionImageIntegration()
%       Currently this function is just a dummy function adds the ROI and
%       Dest together.  This should generate the PIE image.
%
% Main functions: PIE_Gui
%   - Call back functions: myButtonPressDown, myButtonPressUp,
%   myMouseMotion, myKeypress
%
function PIE_Gui(I1, I2,FileName,Method,Color)
%
% -------------------  SELECT SOURCE ROI ------------------
%
disp('USAGE: select a polygon region by using left mouse clicks to draw the vertices');
disp('       right nouse click to finish - you can then drag the selected region around.');
disp('       When you have it where you want it, double click left to cut it and paste to the next image.');
h = figure('MenuBar', 'none', 'Toolbar', 'none');  % open window
[BW, xi, yi] = roipoly(I1);                         % this returns a binary image with white (1) in the mask
% extract mask (crop image)
[r,c] = find(BW == 1);                      % find the max values
maxH = max(r) - min(r);                     % extract the height
maxW = max(c) - min(c);                     % extract the width
SRC = imcrop(I1,[min(c) min(r) maxW maxH]);  % crop the image in the RIO
% crop mask - make the mask RGB (3 layers)
Mc = zeros(size(SRC));                       % make a copy of Ic
Mc(:,:,1) = imcrop(BW,[min(c) min(r) maxW maxH]);
Mc(:,:,2) = imcrop(BW,[min(c) min(r) maxW maxH]);
Mc(:,:,3) = imcrop(BW,[min(c) min(r) maxW maxH]);
%
% NOW SELECT PLACE TO PASTE
%
imshow(I2);
title('Click and drag region to desired location. Press any key to integrate, press q to quit');
lh = line(xi, yi, 'Marker','.','LineStyle','-', 'Color', 'r', 'LineWidth',2);
% Set up units and callback functions
set(h, 'Units', 'pixels');
set(h,'WindowButtonDownFcn',@myButtonPressDown);
set(h,'WindowButtonUpFcn',@myButtonPressUp);
set(h, 'WindowButtonMotionFcn', @myMouseMotion);
set(h, 'KeyPressFcn', @myKeyPress);
myData.xi = xi-min(xi);
myData.yi = yi-min(yi);
myData.SRC = SRC;
myData.DEST = I2;
myData.Mc=Mc;
myData.pressDown = 0;
myData.line = lh;
myData.curX = -1;
myData.curY = -1;
myData.Method=Method;
myData.Color=Color;
myData.FileName=FileName;
set(h, 'UserData', myData);
return
%%
% When button is pressed, call this function
%
function myButtonPressDown(obj,event_obj)
myData = get(obj, 'UserData');      % get the user data (variable name does not have to be "myData"
myData.pressDown = 1;               % set mouse press = true
p = get(gca,'CurrentPoint');        % get current position of mouse on the image
curX = p(1,1);                      % extract the X position (it's a floating point value)
curY = p(1,2);                      % extract the Y positions
myData.curX = curX;
myData.curY = curY;
set(myData.line,'XData', myData.xi+curX, 'YData', myData.yi+curY);
% Save the myData variable back to the object
set(obj, 'UserData', myData);
return
%%
% When button is released, call this function
%
function myButtonPressUp(obj,event_obj)
myData = get(obj, 'UserData');  % get the user data
myData.pressDown = 0;           % set mouse press to be false
set(obj, 'UserData', myData);   % set the uer data (i.e. record mouse is not longer being pressed)
return
%%
% Called anytime the mouse is moved
%
function myMouseMotion(obj,event_obj)
myData = get(obj, 'UserData');  % get the user data
if (myData.pressDown == 1)              % we are only interested if the mouse is down
    p = get(gca,'CurrentPoint');        % get the current point from the image
    curX = p(1,1);                      % extract the point from the strange matlab datastructure return by previous line of code
    curY = p(1,2);
    set(myData.line,'XData', myData.xi+curX, 'YData', myData.yi+curY);
    myData.curX = curX;
    myData.curY = curY;
    set(obj, 'UserData', myData);
end
return
%%
% Call when key any pressed any key
%
function myKeyPress(obj, event_obj)
if (event_obj.Key == 'q')
    close(obj);
    return;
end
% Update the userdata in the object
myData = get(obj, 'UserData');
if (myData.pressDown == 0)          % if mouse is not pressed
    
    if (myData.curX == -1)
        disp('Select a location');
        return;
    end
    
    %
    % Get the source and destination image
    % Compute a new image (SImage) where the source is translated to
    % the correct position based on the last mouse position.
    %
    %
    DEST = myData.DEST;
    SRC = myData.SRC;
    tx = round(myData.curX);
    ty = round(myData.curY);
    
    [hh ww depth] = size(SRC);
    
    %use only the ROI
    TRG(:,:,1)=DEST( ty:(ty+hh-1), tx:(tx+ww-1), 1 );
    TRG(:,:,2)=DEST( ty:(ty+hh-1), tx:(tx+ww-1), 2 );
    TRG(:,:,3)=DEST( ty:(ty+hh-1), tx:(tx+ww-1), 3 );
    
    %         SImage( ty:(ty+hh-1), tx:(tx+ww-1), 1 ) =  SRC(:,:,1);
    %         SImage( ty:(ty+hh-1), tx:(tx+ww-1), 2 ) =  SRC(:,:,2);
    %         SImage( ty:(ty+hh-1), tx:(tx+ww-1), 3 ) =  SRC(:,:,3);
    
    
    Mc = rgb2gray(myData.Mc);
    Mc(1,:)=0;
    Mc(end,:)=0;
    Mc(:,1)=0;
    Mc(:,end)=0;
    se = strel('disk',5);
    Mc = imerode( Mc,se);
    % Call the PIE function.  It will returned the integrated image
    newI = PIE( TRG,SRC,Mc,myData.Method,myData.Color);
    %reconstruct
    if size(newI,3)==1
        DEST=rgb2gray(DEST);
        DEST( ty:(ty+hh-1), tx:(tx+ww-1)) =  newI(:,:);
    else
        DEST( ty:(ty+hh-1), tx:(tx+ww-1), 1 ) =  newI(:,:,1);
        DEST( ty:(ty+hh-1), tx:(tx+ww-1), 2 ) =  newI(:,:,2);
        DEST( ty:(ty+hh-1), tx:(tx+ww-1), 3 ) =  newI(:,:,3);
    end
    %PossionImageIntegration(SImage, DEST, tx, ty, ww, hh);
    imshow(DEST);
    imwrite(DEST,myData.FileName);
end
return

 

 

%PIE
function im_out = PIE( im_target,im_source,im_mask,m,c )
%PIE function: blends the source image with the target one based on the
%boundary given as a BW mask using Poisson Image Editing (PIE)
%  -Usage-
%	im_out = PIE(targetImage,sourceImage,mask,0,0); %for seamless cloning
%	(true color)
%   im_out = PIE(targetImage,sourceImage,mask,1,0); %for mixing gradients
%   (true color)
%	im_out = PIE(targetImage,sourceImage,mask,0,1); %for seamless cloning
%	(grayscale)
%   im_out = PIE(targetImage,sourceImage,mask,1,1); %for mixing gradients
%   (grayscale)
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Citation:
%   Pérez, Patrick, Michel Gangnet, and Andrew Blake.
%   "Poisson image editing." ACM Transactions on Graphics (TOG). Vol. 22.
%   No. 3. ACM, 2003.
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  -Inputs-
%	 im_target: target image
%    im_source: source image
%    im_mask: mask image
%    m: 0 for seamless cloning (default), and 1 for mixing gradients.
%    c: 0 for true color source and target images (default), and 1 for
%    grayscale source and target images.
%  -Outputs-
%    im_out: output image after blending the source image with the source
%    image based on the given mask (uint8).
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialization
if nargin<3
    error('Please, use  PIE(targetImage,sourceImage,mask)');
elseif nargin<4
    m=0; %use the default method (seamless cloning)
    c=0; %use the default image type (true color)
elseif nargin<5
    c=0; %use the default image type (true color)
end
if size(im_target,1)~= size(im_source,1) || ...
        size(im_target,1)~= size(im_mask,1) || ...
        size(im_target,2)~= size(im_source,2) || ...
        size(im_target,2)~= size(im_mask,2) || ...
        size(im_target,3)~=size(im_source,3)
    error('Please, use images with the same size');
end
%if the mask is not grayscale, convert it
if size(im_mask,3)>1
    im_mask=rgb2gray(im_mask);
end
%convert source and target images to double for more precise computations
im_target=double(im_target);
im_source=double(im_source);
%if we are working with true color, let m=3 otherwise, m=1
if c==0 %true color images
    c=3; %for the next for loop
else
    c=1;
    if size(im_source,3)>1
        im_source=rgb2gray(im_source);
    end
    if size(im_target,3)>1
        im_target=rgb2gray(im_target);
    end
    
end
%initially, output image = target image
im_out=im_target;
%create the laplacian mask for the second derivative of the source image
laplacian_mask=[0 1 0; 1 -4 1; 0 1 0];
%normalize the mask image to assure that unknown pixels = 1
im_mask=mat2gray(im_mask);
%convert it to logical to ignore any fractions (soft masks)
im_mask=im_mask==1;
%find the number of unknown pixels based on the mask
n=size(find(im_mask==1),1);
%create look up table
map=zeros(size(im_mask));
%loop through the mask image to initialize the look up table for mapping
counter=0;
for x=1:size(map,1)
    for y=1:size(map,2)
        if im_mask(x,y)==1 %is it unknow pixel?
            counter=counter+1;
            map(x,y)=counter;  %map from (x,y) to the corresponding pixel
            %in the 1D vector
        end
    end
end
for i=1:c %for each color channel
    
    % loop through the mask image again to:
    %1- initialize the coefficient matrix
    %2- initialize the B vector
    
    %if the method is seamless cloning; so, intially put B= (-) laplacian of
    %im_source,
    %otherwise (mixing gradients), B= (-) max(laplacian of im_source, laplacian
    %of im_target)
    
    %create the coefficient matrix A
    
    %At most, there are 5 coefficients per row according to eq (3)
    %in the report
    coeff_num=5;
    
    %create the sparse matrix to save memory
    A=spalloc(n,n,n*coeff_num);
    
    %create the right hand side of the linear system of equations (AX=B)
    B=zeros(n,1);
    
    if m==1  % mixing gradients
        
        %create the gradient mask for the first derivative
        grad_mask_x=[-1 1];
        grad_mask_y=[-1;1]; 
        
        %get the first derivative of the target image
        g_x_target=conv2(im_target(:,:,i),grad_mask_x, 'same');
        g_y_target=conv2(im_target(:,:,i),grad_mask_y, 'same');
        g_mag_target=sqrt(g_x_target.^2+g_y_target.^2);
        
        %get the first derivative of the source image
        g_x_source=conv2(im_source(:,:,i),grad_mask_x, 'same');
        g_y_source=conv2(im_source(:,:,i),grad_mask_y, 'same');
        g_mag_source=sqrt(g_x_source.^2+g_y_source.^2);
        
        %work with 1-D
        g_mag_target=g_mag_target(:);
        g_mag_source=g_mag_source(:);
        
        %initialize the final gradient with the source gradient
        g_x_final=g_x_source(:);
        g_y_final=g_y_source(:);
        
        %if the gradient of the target image is larger than the gradient of
        %the source image, use the target's gradient instead
        g_x_final(abs(g_mag_target)>abs(g_mag_source))=...
            g_x_target(g_mag_target>g_mag_source);
        g_y_final(abs(g_mag_target)>abs(g_mag_source))=...
            g_y_target(g_mag_target>g_mag_source);
        
        %map to 2-D
        g_x_final=reshape(g_x_final,size(im_source,1),size(im_source,2));
        g_y_final=reshape(g_y_final,size(im_source,1),size(im_source,2));
        
        %get the final laplacian of the combination between the source and
        %target images lap=second deriv of x + second deriv of y
        lap=conv2(g_x_final,grad_mask_x, 'same');
        lap=lap+conv2(g_y_final,grad_mask_y, 'same');
        
    else
        %create the laplacian of the source image
        lap=conv2(im_source(:,:,i),laplacian_mask, 'same');
    end
    counter=0;
    for x=1:size(map,1)
        for y=1:size(map,2)
            if im_mask(x,y)==1
                counter=counter+1;
                A(counter,counter)=4; %the diagonal represent the current pixel
                
                %check the boundary
                if im_mask(x-1,y)==0 %known left pixel
                    B(counter)=im_target(x-1,y,i); %add it to B
                else %unknown boundary
                    A(counter,map(x-1,y))=-1; %set its coefficient to -1
                end
                if im_mask(x+1,y)==0 %known right pixel
                    B(counter)=B(counter)+im_target(x+1,y,i); %add it to B
                else %unknown boundary
                    A(counter,map(x+1,y))=-1; %set its coefficient to -1
                end
                if im_mask(x,y-1)==0 %known bottom pixel
                    B(counter)=B(counter)+im_target(x,y-1,i); %add it to B
                else %unknown boundary
                    A(counter,map(x,y-1))=-1; %set its coefficient to -1
                end
                if im_mask(x,y+1)==0 %known top pixel
                    B(counter)=B(counter)+im_target(x,y+1,i); %add it to B
                else %unknown boundary
                    A(counter,map(x,y+1))=-1; %set its coefficient to -1
                end
                
                %update the B vector with the laplacian value
                
                B(counter)=B(counter)-lap(x,y);
                
            end
        end
    end
    
    %solve the linear system of equation
    X=A\B;
    
    
    %reshape X to restore the output image
    
    %     counter=0;
    %     for x=1:size(map,1)
    %         for y=1:size(map,2)
    %             if im_mask(x,y)==1
    %                 counter=counter+1;
    %                 im_out(x,y,i)=X(counter);
    %             end
    %         end
    %     end
    
    
    for counter=1:length(X)
        [index_x,index_y]=find(map==counter);
        im_out(index_x,index_y,i)=X(counter);
        
    end
    
    
    %release all
    clear A B X lap_source lap_target g_mag_source g_mag_target
end
im_out=uint8(im_out);
end

不谢~