% multiplePointDistort.m
% 
% Zach Gildersleeve
% Januaru 26, 2007
%
% This script displays an image, and allows the user to select starting
% points and a landmark points.  The image is transformed to warp the
% starting points to the landmark points using a gaussian radial basis
% function. The user can select as many points as desired, and hits the
% return key to transform the image.  If an uneven number of coordinates
% are entered, the last coordinate is ignored.


% Read in the image
myImage = imread('../images/photo1-2.png');

% Convert the image to double between 0 and 1
myImage = (double(myImage)/255);

% Define the cordinate axis of the image
% For simplicity assume that the image Domain is [-1 1]
[X,Y] = meshgrid(linspace(-1,1,512),linspace(-1,1,512));

% Define sigma
sigma = 0.25;

% Create the figure and draw the original image
figure('Name', 'Single Point Warp', 'Position', [0 0 512 512]);
imagesc(linspace(-1, 1, 512), linspace(-1, 1, 512), myImage);
colormap('gray');
axis('square');

% Get multiple coordinates from the user's mouse
% in sets of two, where the first click point 1,
% the second click is the translated point 1, and so on.
[L1, L2] = ginput();

% Ignore the last coordinate if an uneven number are entered
if rem(size(L2, 1), 2) ~= 0
    L1(size(L2, 1)) = [];
    L2(size(L2, 1)) = [];
end

% Initialize U1 and U2 to zeros, just for safety
U1 = zeros(512, 512);
U2 = zeros(512, 512);

for i=1:2:size(L2, 1)
    % Compute the alpha values
    alpha1 = (L1(i+1) - L1(i));
    alpha2 = (L2(i+1) - L2(i));
    
    % Calculate the transform [U1, U2] and add
    % it to the previous transform
    U1 = U1 + alpha1 .* exp(-((X - L1(i+1)).^2 + (Y - L2(i+1)).^2) / sigma^2);
    U2 = U2 + alpha2 .* exp(-((X - L1(i+1)).^2 + (Y - L2(i+1)).^2) / sigma^2);
end

% Interpolate the transformed image
% We use X - U1, Y - U2 to transform in correct direction
distortedImage = interp2(X, Y, myImage, X - U1, Y - U2, 'linear');
imagesc(distortedImage);

colormap('gray');
axis('square');

% Draw line from landmark point location to translated point
if 1
    newL1 = L1 * 256 + 256;
    newL2 = L2 * 256 + 256;

    hold on;

    for i=1:2:size(L2, 1);
        plot(newL1(i:i+1), newL2(i:i+1), 'r--x');
    end
end