diff --git a/m/1.jpg b/m/1.jpg
new file mode 100644
index 0000000..a5002ed
Binary files /dev/null and b/m/1.jpg differ
diff --git a/m/2.jpg b/m/2.jpg
new file mode 100644
index 0000000..df59b7c
Binary files /dev/null and b/m/2.jpg differ
diff --git a/m/MatchColumns.m b/m/MatchColumns.m
new file mode 100644
index 0000000..9f37c72
--- /dev/null
+++ b/m/MatchColumns.m
@@ -0,0 +1,71 @@
+function [U_t,A_t] = MatchColumns(us,ut,at)
+% Ensures consistency between the target & source image rotation matrices.
+
+% This routine matches columns in the target image rotation matrix to those
+% in the source image rotation matrix.
+
+% Each rotation matrix is derived by undertaking a singular value
+% decomposition of the respective cross covariance matrices. The
+% outcome of the decomposition is ordered within the rotation matrix in 
+% the order of the descending singular values.  This often leads to 
+% compatible rotation matrices but that is not guaranteed.  Sometimes the 
+% colour ordering of one rotation matrix may be different from the other.
+%
+% Additionally even when the matrices do correspond in orientation, they
+% need not correspond in direction. The direction axis rotation in one
+% matrix may be the negative of the rotation in the other.
+
+% Rotations of the individual colour axes are given by the columns of the
+% rotation matrices.  In the following processing, all rearrangements are 
+% considered of the columns in the target rotation matrix 'ut', to find 
+% the arrangement that best matches the source rotation matrix 'us'. 
+% Matching is measured by taking the vector dot products of the 
+% corresponding matrix columns and finding the arrangement with the 
+% largest sum of absolute dot product values.  Absolute values are used 
+% to accommodate axis pairs that have similar orientations but different 
+% directions.
+%
+% Once the best match has been found this is taken as the correct target
+% image rotation matrix.  Columns are negated where the vector cross
+% product has a negative value, to ensure direction compatibility.  The 
+% singular value matrix for the target image is reordered to match the 
+% reordering of the rotation matrix.
+%
+% This processing method and routine copyright Dr T E Johnson 2020.
+% terence.johnson@gmail.com
+
+% All 6 permutations of the three columns need to be addressed.
+perm = perms([1 2 3]);
+     
+max=0.0;
+
+for i=1:6
+    % Compute the sum of the absolute dot products for the matrix columns. 
+    sum=abs(ut(:,perm(i,1)).'*us(:,1)) + ...
+        abs(ut(:,perm(i,2)).'*us(:,2)) + ...
+        abs(ut(:,perm(i,3)).'*us(:,3)) ;
+
+    if(sum>max)
+        % Best aligned axes give the biggest sum.
+        max=sum;
+        bestperm=i;
+    end
+end
+    % Now recover the best permutation and implement it.
+    b1=perm(bestperm,1);
+    b2=perm(bestperm,2);
+    b3=perm(bestperm,3);
+ 
+    % Set the rotation matrix columns to the new order, changing the 
+    % direction of rotations if necessary by negating column values.
+    U_t(:,1)=ut(:,b1)*sign(ut(:,b1).'*us(:,1));
+    U_t(:,2)=ut(:,b2)*sign(ut(:,b2).'*us(:,2));
+    U_t(:,3)=ut(:,b3)*sign(ut(:,b3).'*us(:,3));    
+    
+    % Similarly reorder the singular values.
+    A_t(1,1)=at(b1,b1);
+    A_t(2,2)=at(b2,b2);
+    A_t(3,3)=at(b3,b3);     
+end
+
+
diff --git a/m_ruggedisation_update/1.jpg b/m_ruggedisation_update/1.jpg
new file mode 100644
index 0000000..a5002ed
Binary files /dev/null and b/m_ruggedisation_update/1.jpg differ
diff --git a/m_ruggedisation_update/2.jpg b/m_ruggedisation_update/2.jpg
new file mode 100644
index 0000000..df59b7c
Binary files /dev/null and b/m_ruggedisation_update/2.jpg differ
diff --git a/m_ruggedisation_update/3.jpg b/m_ruggedisation_update/3.jpg
new file mode 100644
index 0000000..1464ccd
Binary files /dev/null and b/m_ruggedisation_update/3.jpg differ
diff --git a/m_ruggedisation_update/4.jpg b/m_ruggedisation_update/4.jpg
new file mode 100644
index 0000000..7ec6a9c
Binary files /dev/null and b/m_ruggedisation_update/4.jpg differ
diff --git a/m_ruggedisation_update/MatchColumns.m b/m_ruggedisation_update/MatchColumns.m
new file mode 100644
index 0000000..9f37c72
--- /dev/null
+++ b/m_ruggedisation_update/MatchColumns.m
@@ -0,0 +1,71 @@
+function [U_t,A_t] = MatchColumns(us,ut,at)
+% Ensures consistency between the target & source image rotation matrices.
+
+% This routine matches columns in the target image rotation matrix to those
+% in the source image rotation matrix.
+
+% Each rotation matrix is derived by undertaking a singular value
+% decomposition of the respective cross covariance matrices. The
+% outcome of the decomposition is ordered within the rotation matrix in 
+% the order of the descending singular values.  This often leads to 
+% compatible rotation matrices but that is not guaranteed.  Sometimes the 
+% colour ordering of one rotation matrix may be different from the other.
+%
+% Additionally even when the matrices do correspond in orientation, they
+% need not correspond in direction. The direction axis rotation in one
+% matrix may be the negative of the rotation in the other.
+
+% Rotations of the individual colour axes are given by the columns of the
+% rotation matrices.  In the following processing, all rearrangements are 
+% considered of the columns in the target rotation matrix 'ut', to find 
+% the arrangement that best matches the source rotation matrix 'us'. 
+% Matching is measured by taking the vector dot products of the 
+% corresponding matrix columns and finding the arrangement with the 
+% largest sum of absolute dot product values.  Absolute values are used 
+% to accommodate axis pairs that have similar orientations but different 
+% directions.
+%
+% Once the best match has been found this is taken as the correct target
+% image rotation matrix.  Columns are negated where the vector cross
+% product has a negative value, to ensure direction compatibility.  The 
+% singular value matrix for the target image is reordered to match the 
+% reordering of the rotation matrix.
+%
+% This processing method and routine copyright Dr T E Johnson 2020.
+% terence.johnson@gmail.com
+
+% All 6 permutations of the three columns need to be addressed.
+perm = perms([1 2 3]);
+     
+max=0.0;
+
+for i=1:6
+    % Compute the sum of the absolute dot products for the matrix columns. 
+    sum=abs(ut(:,perm(i,1)).'*us(:,1)) + ...
+        abs(ut(:,perm(i,2)).'*us(:,2)) + ...
+        abs(ut(:,perm(i,3)).'*us(:,3)) ;
+
+    if(sum>max)
+        % Best aligned axes give the biggest sum.
+        max=sum;
+        bestperm=i;
+    end
+end
+    % Now recover the best permutation and implement it.
+    b1=perm(bestperm,1);
+    b2=perm(bestperm,2);
+    b3=perm(bestperm,3);
+ 
+    % Set the rotation matrix columns to the new order, changing the 
+    % direction of rotations if necessary by negating column values.
+    U_t(:,1)=ut(:,b1)*sign(ut(:,b1).'*us(:,1));
+    U_t(:,2)=ut(:,b2)*sign(ut(:,b2).'*us(:,2));
+    U_t(:,3)=ut(:,b3)*sign(ut(:,b3).'*us(:,3));    
+    
+    % Similarly reorder the singular values.
+    A_t(1,1)=at(b1,b1);
+    A_t(2,2)=at(b2,b2);
+    A_t(3,3)=at(b3,b3);     
+end
+
+
diff --git a/m_ruggedisation_update/cf_Xiao06.m b/m_ruggedisation_update/cf_Xiao06.m
new file mode 100644
index 0000000..7663889
--- /dev/null
+++ b/m_ruggedisation_update/cf_Xiao06.m
@@ -0,0 +1,70 @@
+function est_im = cf_Xiao06(source,target)
+%CF_XIAO computes Reinhard's image colour transfer
+%
+%   CF_XIAO(SOURCE,TARGET) returns the colour transfered source
+%   image SOURCE according to the target image TARGET.
+%
+
+%   References:
+%   Xiao, Xuezhong, and Lizhuang Ma. "Color transfer in correlated color
+%   space." % In Proceedings of the 2006 ACM international conference on
+%   Virtual reality continuum and its applications, pp. 305-309. ACM, 2006.
+
+%   TODO: Swatch-based transfer is not implemented (but I think it is not
+%   important) Anyone interested is welcome to contribute. :-)
+%                                                           -- Han Gong
+
+%   Copyright 2015 Han Gong <gong@fedoraproject.org>, University of East
+%   Anglia.
+
+rgb_s = reshape(im2double(source),[],3)';
+rgb_t = reshape(im2double(target),[],3)';
+
+% compute mean
+mean_s = mean(rgb_s,2);
+mean_t = mean(rgb_t,2);
+
+% compute covariance
+cov_s = cov(rgb_s');
+cov_t = cov(rgb_t');
+
+% decompose covariances
+[U_s,A_s,~] = svd(cov_s);
+[U_t,A_t,~] = svd(cov_t);
+
+%**********************************************************
+% Processing modification by T E Johnson.
+% Set true (default) to implement ruggedisation processing.
+% Set false for standard (original) processing.
+ruggedisation=true;
+
+if (ruggedisation)
+    [U_t,A_t]=MatchColumns(U_s,U_t,A_t);
+end 
+%**********************************************************
+
+rgbh_s = [rgb_s;ones(1,size(rgb_s,2))];
+
+% compute transforms
+% translations
+T_t = eye(4); T_t(1:3,4) = mean_t;
+T_s = eye(4); T_s(1:3,4) = -mean_s;
+% rotations
+R_t = blkdiag(U_t,1); R_s = blkdiag(inv(U_s),1);
+% scalings
+% NOTE!
+% I believe the author has a typo in his original paper
+% The following is the original way to compute S_t, but
+% the result does not seem correct
+% S_t = diag([diag(A_t);1]); 
+% I added a 0.5 power to correct it.
+S_t = diag([diag(A_t).^(0.5);1]);
+S_s = diag([diag(A_s).^(-0.5);1]);
+
+rgbh_e = T_t*R_t*S_t*S_s*R_s*T_s*rgbh_s; % estimated RGBs
+rgbh_e = bsxfun(@rdivide, rgbh_e, rgbh_e(4,:));
+rgb_e = rgbh_e(1:3,:);
+
+est_im = reshape(rgb_e',size(source));
+
+end
diff --git a/m_ruggedisation_update/demo.m b/m_ruggedisation_update/demo.m
new file mode 100644
index 0000000..7af5106
--- /dev/null
+++ b/m_ruggedisation_update/demo.m
@@ -0,0 +1,22 @@
+% Demonstration of Xiao's Image Colour Transfer
+
+% Copyright 2015 Han Gong <gong@fedoraproject.org>, University of East
+% Anglia.
+
+% Now ruggedised by T E Johnson (April 2020) to ensure consistent 
+% processing.
+
+% References:
+% Xiao, Xuezhong, and Lizhuang Ma. "Color transfer in correlated color
+% space." % In Proceedings of the 2006 ACM international conference on
+% Virtual reality continuum and its applications, pp. 305-309. ACM, 2006.
+
+I0 = im2double(imread('4.jpg'));
+I1 = im2double(imread('3.jpg'));
+
+IR = cf_Xiao06(I0,I1);
+
+figure; 
+subplot(1,3,1); imshow(I0); title('Original Image'); axis off
+subplot(1,3,2); imshow(I1); title('Target Palette'); axis off
+subplot(1,3,3); imshow(IR); title('Result After Colour Transfer'); axis off