Matlab is of course the perfect way to show this as the process stays linear and is a matrix fest, you could say matlab was written to do this

Edit: Tried to make smileys go away but didn't work.
%
% 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,


datamax(i) = max([D0R(i,


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


p1R{i} = hist(D1R(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,


end
D0 = relaxation * R \ (D0R_ - D0R) + D0;
end
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