%     This code implements the recoloring method described in the paper:
%     Musialski, P., Cui, M., Ye, J., Razdan, A., & Wonka, P. (2012). 
%     A Framework for Interactive Image Color Editing. 
%     The Visual Computer, 29(11), 1173–1186. doi:10.1007/s00371-012-0761-5

%     Copyright (C) 2011-2013, Przemyslaw Musialski, pm@cg.tuwien.ac.at
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% -------------------------------------------------------------------------


function handles    = InitLandmarks(handles)
% inits the subsampling and the LLE-weights


% uncomment to fix the random seed
%RandStream.setGlobalStream(RandStream('mt19937ar','Seed',0));

ratio_target        = handles.ratio_target;
ratio_nontarget     = handles.ratio_nontarget;

mask                = handles.mask;
targetindex         = find(mask==1);
nontargetindex      = find(mask==0);

nontargetsample     = randsample(length(nontargetindex),round(length(nontargetindex)*ratio_nontarget));
targetsample        = randsample(length(targetindex),round(length(targetindex)*ratio_target));

nontargetsample     = nontargetindex(nontargetsample);
targetsample        = targetindex(targetsample);
allsample           = [targetsample;nontargetsample]; %less points will make it harder to drag to rgb

% samples are allways unique and that the convex hull is included
srcLab              = handles.srcLab;
landmarks           = srcLab(allsample,:);

% boundary of the CIALAB cube (lightness limited to 100)
pointsBoundary = [
    0, -128,-128;
    0, -128,+128;
    0, +128,-128;
    0, +128,+128;
    100, -128,-128;
    100, -128,+128;
    100, +128,-128;
    100, +128,+128;
    ];

% compute DT
fprintf(1,'Computing Delunay-Triangulation \n'); 
t1=tic;

% get samples landmarks
[~, landim, ~]      = unique(landmarks,'rows','first');
allU                = allsample((landim));
tarSample           = intersect(allU, targetindex);
nonSample           = intersect(allU, nontargetindex);
allsample           = [tarSample; nonSample];

    % compute DT
    landmarks       = srcLab(allsample,:);    
    pointsDt        = [landmarks; pointsBoundary];
    %tic
    % for matlab 2013 and newer
    %dTri            = delaunayTriangulation(pointsDt);
    % for older matlab
    dTri            = DelaunayTri(pointsDt);
    %toc
    %tic
    [SI, BC]        = pointLocation(dTri,srcLab);
    %toc
    
    handles.dTri    = dTri;
    handles.BC      = BC;
    handles.BA      = pointsBoundary;
    % for MATLAB 2013 and newer
    %handles.BI      = dTri.ConnectivityList(SI,:);
    % for older MATLAB
    handles.BI      = dTri.Triangulation(SI,:);
    X               = landmarks;

targetsample        = tarSample;
nontargetsample     = nonSample;

fprintf(1,['...took ' num2str(toc(t1)) 's.\n']);

handles.W           = SolveLLE_Weights(X',handles.K,handles.tol);
handles.targetsample= targetsample;
handles.nontargetsample = nontargetsample;
handles.allsample   = allsample;
handles.landmarks   = landmarks;


% visualize the landmarks
if(1)
    samplesimg      = zeros(size(mask));    
    %samplesimg(mask==0) = 0.25;
    samplesimg(nontargetsample) = 1.0;
    samplesimg(targetsample) = 1.0;    
 
    k = ones(3,3);
    samplesimg = conv2(samplesimg,k,'same')>0;
    
    colorimg = handles.rgbimg;
    colorimg = colorimg.*repmat(samplesimg,[1,1,3]);            
    samplesimg = im2double((1-(~samplesimg))==0);
    colorimg = colorimg+repmat(samplesimg,[1,1,3]);
    handles.samplesimg   = colorimg;   
end;



function [result] = SolveLLE_Weights(X,K,tol)
% solves for optimal LLE-weights
% based on code from: http://www.cs.nyu.edu/~roweis/lle/

[D,N] = size(X);
fprintf(1,'LLE running on %d points in %d dimensions\n',N,D);
t1 = tic;
% ----------------------------------------------------------------------
% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS 
fprintf(1,'-->Finding %d nearest neighbours.\n',K);
nnidx           = knnsearch(X',X','k',K+1)';
neighborhood    = nnidx(2:(1+K),:);

% ----------------------------------------------------------------------
% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.\n');
W = sparse(N,N);
for ii=1:N
   z = (X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K)); % shift ith pt to origin
   c = z'*z;                                        % local covariance
   c = c + eye(K)*tol;                              % regularlization (K>D)
   w = c\ones(K,1);   
   W(ii,neighborhood(:,ii)) = w/sum(w);
end;
% ----------------------------------------------------------------------
% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W)
% fprintf(1,'-->Computing embedding.\n');
% wo do not compute the eigenvalues but take the whole matrix instead
I = sparse(1:N,1:N,ones(1,N),N,N,4*K*N); 
M = I-W;
result = M'*M;
fprintf(1,['...took ' num2str(toc(t1)) 's.\n'] );

clear X
clear neighborhood
clear nndist
clear nnid
clear W
clear I
clear B
