function im1a = colortwist(im1, im2)
%function colortwist(im1, im2, gamma)
% finds least-squares linear 3x3 color transform 'twist'
% to convert im1's colors into im2

     [h, w, nb] = size(im1);
     n = w*h;

     c1 = zeros(n, 3);
     c2 = zeros(n, 3);
     for b = 1:3,
	       c1(:,b) = double(reshape(im1(:,:,b), n, 1))/255;
     c2(:,b) = double(reshape(im2(:,:,b), n, 1))/255;
     end

%     c1 = c1 .^ gamma;
im1a = im1;

twist = zeros(3, 3);
c1a = zeros(n, 3);

for b = 1:3,
	  twist(b, :) = (c1 \ c2(:, b))';
    c1a(:, b) = c1 * twist(b, :)';
im1a(:, :, b) = reshape(uint8(255*c1a(:, b)), h, w);
end

twist

%subplot(1,4,1); imshow(im1);
%subplot(1,4,2); imshow(im2);
%subplot(1,4,3); imshow(uint8(128+10*(double(im1a)-double(im2))));
%subplot(1,4,4); imshow(im1a);
%resid = reshape(mean(mean(abs(double(im1a) - double(im2)))), 1, 3)
     
