%
%   simple implementation of N-Dimensional PDF Transfer 
%
%   [DR] = pdf_transferND(D0, D1, rotations);
%
%     D0, D1 = NxM matrix containing N-dimensional features
%     rotations = { {R_1}, ... , {R_n} } with R_i PxN 
%
%     note that we can use more than N projection axes. In this case P > N
%     and the inverse transformation is done by least mean square. 
%     Using more than N axes leads to a more stable (but also slower) 
%     convergence.
%
function [DR] = pdf_transferND(D0, D1, Rotations)

nb_dims = size(D0,1);

relaxation = 1;
DR = D0;

for it=1:length(Rotations)

    R = Rotations{it};    
    nb_projs = size(R,1);
    
    %% apply rotation
    
    D0R = R * D0;
    D1R = R * D1;

    %% find data range
    for i=1:nb_projs
        datamin(i) = min([D0R(i,:) D1R(i,:)]);
        datamax(i) = max([D0R(i,:) D1R(i,:)]);
    end

    %% get the marginals
    for i=1:nb_projs
        step =  (datamax(i) - datamin(i))/300;
        p0R{i} = hist(D0R(i,:), datamin(i):step:datamax(i));
        p1R{i} = hist(D1R(i,:), datamin(i):step:datamax(i));
    end

    %% match the marginals
    for i=1:nb_projs
        f{i} = pdf_transfer1D(p0R{i}, p1R{i});
        scale = (length(f{i})-1)/(datamax(i)-datamin(i));
        D0R_(i,:) = interp1(0:length(f{i})-1, f{i}', (D0R(i,:) - datamin(i))*scale)/scale + datamin(i);
    end

    D0 = relaxation * R \ (D0R_ - D0R) + D0;
end

DR = D0;

end
%% 
%% 1D - PDF Transfer
%%
function f = pdf_transfer1D(pX,pY)
    nbins = max(size(pX));

    PX = cumsum(pX);
    PX = PX/PX(end);

    PY = cumsum(pY);
    PY = PY/PY(end);

    % inversion 
    PX = [0 PX nbins] + (0:nbins+1)*1e-10;
    PY = [0 PY nbins] + (0:nbins+1)*1e-10;

    f = interp1(PY, [0 ((0:nbins-1)+1e-16) (nbins+1e-10)], PX,'linear');
    f = f(2:end-1);
end