function K = mycolortransfer(I, J)

iLAB = RGB2Lab(J);
jLAB = RGB2Lab(I);

% Calculate the means
iLmean = mean2(iLAB(:,:,1));
iAmean = mean2(iLAB(:,:,2));
iBmean = mean2(iLAB(:,:,3));

jLmean = mean2(jLAB(:,:,1));
jAmean = mean2(jLAB(:,:,2));
jBmean = mean2(jLAB(:,:,3));

% New mean values arrays for target
%iLnew = iLAB(:,:,1) - iLmean;
%iAnew = iLAB(:,:,2) - iAmean;
%iBnew = iLAB(:,:,3) - iBmean;

jLnew = jLAB(:,:,1) - jLmean;
jAnew = jLAB(:,:,2) - jAmean;
jBnew = jLAB(:,:,3) - jBmean;

% Calculate the standard deviations
iLstd = std2(iLAB(:,:,1));
iAstd = std2(iLAB(:,:,2));
iBstd = std2(iLAB(:,:,3));

jLstd = std2(jLAB(:,:,1));
jAstd = std2(jLAB(:,:,2));
jBstd = std2(jLAB(:,:,3));

% Calculate the constants
sL = iLstd/jLstd;
sA = iAstd/jAstd;
sB = iBstd/jBstd;

% Calculate LAB prime
jLnew = jLnew * sL;
jAnew = jAnew * sA;
jBnew = jBnew * sB;

% Add in the input image mean
jLnew = jLnew + iLmean;
jAnew = jAnew + iAmean;
jBnew = jBnew + iBmean;

jLnew = double(jLnew);
jAnew = double(jAnew);
jBnew = double(jBnew);
%K = cat(3, jLnew, jAnew, jBnew);

K = Lab2RGB(jLnew, jAnew, jBnew);

end