function Image_Adapt(imgPath1,imgPath2)
%illustration of image adaptation
%this program calculates the mean color and contrasts in two input images
%two output images are then created by adapting the gamut of each image to
%the other such that the average responses of color mechanisms for one image
%equals the average responses for the other image

close all;
clearvars -except imgPath*;
if nargin ~= 2
    [imgName,imgDir] = uigetfile('*.*');
    img1 = imread(strcat(imgDir,imgName));
    [imgName,imgDir] = uigetfile('*.*');
    img2 = imread(strcat(imgDir,imgName));
else
    img1 = imread(imgPath1);
    img2 = imread(imgPath2);
end

%--------------------------------------------------------------------------
%select and define matrices for the two original images
%image 1
subplot(2,2,1), imshow(img1);
title('original image 1');
RGB1 = double(img1);
CIE1 = RGB1;
LMS1 = RGB1;
DKL1 = RGB1;
d = size(DKL1);
MECH1 = zeros(d(1),d(2),26);
DKLadap1 = zeros(d(1),d(2),3);

%image 2
subplot(2,2,2), imshow(img2);
title('original image 2');
RGB2 = double(img2);
CIE2 = RGB2;
LMS2 = RGB2;
DKL2 = RGB2;
d = size(DKL2);
MECH2 = zeros(d(1),d(2),26);
DKLadap2 = zeros(d(1),d(2),3);

%--------------------------------------------------------------------------
%color conversions
%RGB primaries to LMS cones (note this is observer/device dependent)
RGBtoLMS = [0.142336  0.435583   0.061297; ...
            0.030371  0.274839   0.055574; ...
            0.000269  0.001759   0.017736];

LMStoRGB = inv(RGBtoLMS);

%default scaling factors and mean levels for converting MacLeod-Boynton
%coordinates to cone-opponent (DKL) contrasts
LMStoDKL = [2000     5200     1.28; ... %contrast scaling
    .6392    .01976   128];             %monitor gray

%--------------------------------------------------------------------------
%post-receptoral mechanisms defined in terms of their preferred col/lum
%angles in cone-opponent space
mechsen = zeros(26,3);
imech = 1;
for lumangle = -45:45:45
    for colangle = 0:45:315
        mechsen(imech,1) = cosd(lumangle)*cosd(colangle);
        mechsen(imech,2) = cosd(lumangle)*sind(colangle);
        mechsen(imech,3) = sind(lumangle);
        imech = imech+1;
    end
end
%last 2 are pure luminance increment or decrement
mechsen(25,1) = 0;
mechsen(25,2) = 0;
mechsen(25,3) = 1;
mechsen(26,1) = 0;
mechsen(26,2) = 0;
mechsen(26,3) = -1;

%--------------------------------------------------------------------------
%forward RGB to mechanism contrasts: image 1
LMS1(:,:,1) = RGB1(:,:,1).*RGBtoLMS(1,1) + RGB1(:,:,2).*RGBtoLMS(1,2) + RGB1(:,:,3).*RGBtoLMS(1,3);
LMS1(:,:,2) = RGB1(:,:,1).*RGBtoLMS(2,1) + RGB1(:,:,2).*RGBtoLMS(2,2) + RGB1(:,:,3).*RGBtoLMS(2,3);
LMS1(:,:,3) = RGB1(:,:,1).*RGBtoLMS(3,1) + RGB1(:,:,2).*RGBtoLMS(3,2) + CIE1(:,:,3).*RGBtoLMS(3,3);

LMS1mean(1) = mean2(LMS1(:,:,1));
LMS1mean(2) = mean2(LMS1(:,:,2));
LMS1mean(3) = mean2(LMS1(:,:,3));

%cone adaptation: sets mean and contrast thresholds relative to the image
%mean, consistent with scaling (adapting) the cone sensitivities
LMmean1 = LMS1mean(1) / (LMS1mean(1)+LMS1mean(2));
Smean1 = LMS1mean(3) / (LMS1mean(1)+LMS1mean(2));
LUMmean1 = LMS1mean(1)+LMS1mean(2);

%scale contrast sensitivity for the S and LUM mean
%~independent of the LM mean
LMscale1 = LMStoDKL(1,1);
Sscale1 = LMStoDKL(1,2) * LMStoDKL(2,2) / Smean1;
LUMscale1 = LMStoDKL(2,3) / LUMmean1;

DKL1(:,:,1) = LMscale1 * (LMS1(:,:,1) ./ (LMS1(:,:,1)+LMS1(:,:,2)) - LMmean1);
DKL1(:,:,2) = Sscale1 * (LMS1(:,:,3) ./ (LMS1(:,:,1)+LMS1(:,:,2)) - Smean1);
DKL1(:,:,3) = 1/LUMscale1 * ((LMS1(:,:,1)+LMS1(:,:,2)) - LUMmean1);

%average response of each mechanism
mechavgresp1 = zeros(26,1);
for imech = 1:26
    %response of ith mech to each pixel
    MECH1(:,:,imech) = DKL1(:,:,1).*mechsen(imech,1) + DKL1(:,:,2).*mechsen(imech,2) + ...
        DKL1(:,:,3).*mechsen(imech,3);
    %sensitivity is half-rectified cosine so truncate negative responses
    MECH1(:,:,imech) = max(0,MECH1(:,:,imech));
    %average response of each mechanism
    mechavgresp1(imech) = mean2(MECH1);
end

%--------------------------------------------------------------------------
%forward RGB to mechanism contrasts: image 2
LMS2(:,:,1) = RGB2(:,:,1).*RGBtoLMS(1,1) + RGB2(:,:,2).*RGBtoLMS(1,2) + RGB2(:,:,3).*RGBtoLMS(1,3);
LMS2(:,:,2) = RGB2(:,:,1).*RGBtoLMS(2,1) + RGB2(:,:,2).*RGBtoLMS(2,2) + RGB2(:,:,3).*RGBtoLMS(2,3);
LMS2(:,:,3) = RGB2(:,:,1).*RGBtoLMS(3,1) + RGB2(:,:,2).*RGBtoLMS(3,2) + CIE2(:,:,3).*RGBtoLMS(3,3);

LMS2mean(1) = mean2(LMS2(:,:,1));
LMS2mean(2) = mean2(LMS2(:,:,2));
LMS2mean(3) = mean2(LMS2(:,:,3));

%cone adaptation: sets mean and contrast thresholds relative to the image
%mean, consistent with scaling (adapting) the cone sensitivities
LMmean2 = LMS2mean(1) / (LMS2mean(1)+LMS2mean(2));
Smean2 = LMS2mean(3) / (LMS2mean(1)+LMS2mean(2));
LUMmean2 = LMS2mean(1)+LMS2mean(2);

%scale contrast sensitivity for the S and LUM mean
%~independent of the LM mean
LMscale2 = LMStoDKL(1,1);
Sscale2 = LMStoDKL(1,2) * LMStoDKL(2,2) / Smean2;
LUMscale2 = LMStoDKL(2,3) / LUMmean2;

DKL2(:,:,1) = LMscale2 * (LMS2(:,:,1) ./ (LMS2(:,:,1)+LMS2(:,:,2)) - LMmean2);
DKL2(:,:,2) = Sscale2 * (LMS2(:,:,3) ./ (LMS2(:,:,1)+LMS2(:,:,2)) - Smean2);
DKL2(:,:,3) = 1/LUMscale2 * ((LMS2(:,:,1)+LMS2(:,:,2)) - LUMmean2);

%average response of each mechanism
mechavgresp2 = zeros(26,1);
for imech = 1:26
    %response of ith mech to each pixel
    MECH2(:,:,imech) = DKL2(:,:,1).*mechsen(imech,1) + DKL2(:,:,2).*mechsen(imech,2) + ...
        DKL2(:,:,3).*mechsen(imech,3);
    %sensitivity is half-rectified cosine so truncate negative responses
    MECH2(:,:,imech) = max(0,MECH2(:,:,imech));
    %average response of each mechanism
    mechavgresp2(imech) = mean2(MECH2);
end

%--------------------------------------------------------------------------
%contrast adaptation: rescale the mechanism sensitivities so average response 
%of each mechanism equals its average response to other image
mechscale1 = ones(26,1);
mechscale2 = ones(26,1);
mechscale1(:) = mechavgresp2(:) ./ mechavgresp1(:);
mechscale2(:) = mechavgresp1(:) ./ mechavgresp2(:);

%--------------------------------------------------------------------------
%rebuild the cardinal axis responses from the adapted mechanisms: image 1
%division by 4 or 5 adjusts for mult mechs sensitive to each axis
for imech = 1:26
    DKLadap1(:,:,1) = DKLadap1(:,:,1) + (MECH1(:,:,imech).* mechscale1(imech).* ...
        mechsen(imech,1) / 4);
    DKLadap1(:,:,2) = DKLadap1(:,:,2) + (MECH1(:,:,imech) .* mechscale1(imech).* ...
        mechsen(imech,2) / 4);
    DKLadap1(:,:,3) = DKLadap1(:,:,3) + (MECH1(:,:,imech) .* mechscale1(imech).* ...
        mechsen(imech,3) / 5);
end

%backward DKL to RGB
LMS1(:,:,1) = (DKLadap1(:,:,1)./LMscale2+ LMmean2).*(DKLadap1(:,:,3)*LUMscale2+LUMmean2);
LMS1(:,:,2) = (1-(DKLadap1(:,:,1)./LMscale2+ LMmean2)).*(DKLadap1(:,:,3)*LUMscale2+LUMmean2);
LMS1(:,:,3) = (DKLadap1(:,:,2)./Sscale2+Smean2).*(DKLadap1(:,:,3)*LUMscale2+LUMmean2);

%truncate negative cone excitations (for luminance decrements < 0)
LMS1(:,:,:) = max(0,LMS1(:,:,:));

RGB1(:,:,1) = LMS1(:,:,1).*LMStoRGB(1,1) + LMS1(:,:,2).*LMStoRGB(1,2) + LMS1(:,:,3).*LMStoRGB(1,3);
RGB1(:,:,2) = LMS1(:,:,1).*LMStoRGB(2,1) + LMS1(:,:,2).*LMStoRGB(2,2) + LMS1(:,:,3).*LMStoRGB(2,3);
RGB1(:,:,3) = LMS1(:,:,1).*LMStoRGB(3,1) + LMS1(:,:,2).*LMStoRGB(3,2) + LMS1(:,:,3).*LMStoRGB(3,3);

imgadap1 = uint8(RGB1);

subplot(2,2,3),imshow(imgadap1);
title('image 1 adapted to image 2');
imwrite(imgadap1,'adapted1.jpg');

%--------------------------------------------------------------------------
%rebuild the cardinal axis responses from the adapted mechanisms: image 2
%division by 4 or 5 adjusts for mult mechs sensitive to each axis
for imech = 1:26
    DKLadap2(:,:,1) = DKLadap2(:,:,1) + (MECH2(:,:,imech).* mechscale2(imech).* ...
        mechsen(imech,1) / 4);
    DKLadap2(:,:,2) = DKLadap2(:,:,2) + (MECH2(:,:,imech) .* mechscale2(imech).* ...
        mechsen(imech,2) / 4);
    DKLadap2(:,:,3) = DKLadap2(:,:,3) + (MECH2(:,:,imech) .* mechscale2(imech).* ...
        mechsen(imech,3) / 5);
end

%backward DKL to RGB
LMS2(:,:,1) = (DKLadap2(:,:,1)./LMscale1+ LMmean1).*(DKLadap2(:,:,3)*LUMscale1+LUMmean1);
LMS2(:,:,2) = (1-(DKLadap2(:,:,1)./LMscale1+ LMmean1)).*(DKLadap2(:,:,3)*LUMscale1+LUMmean1);
LMS2(:,:,3) = (DKLadap2(:,:,2)./Sscale1+Smean1).*(DKLadap2(:,:,3)*LUMscale1+LUMmean1);

%truncate negative cone excitations(for luminance decrements < 0)
LMS2(:,:,:) = max(0,LMS2(:,:,:));

RGB2(:,:,1) = LMS2(:,:,1).*LMStoRGB(1,1) + LMS2(:,:,2).*LMStoRGB(1,2) + LMS2(:,:,3).*LMStoRGB(1,3);
RGB2(:,:,2) = LMS2(:,:,1).*LMStoRGB(2,1) + LMS2(:,:,2).*LMStoRGB(2,2) + LMS2(:,:,3).*LMStoRGB(2,3);
RGB2(:,:,3) = LMS2(:,:,1).*LMStoRGB(3,1) + LMS2(:,:,2).*LMStoRGB(3,2) + LMS2(:,:,3).*LMStoRGB(3,3);

imgadap2 = uint8(RGB2);
subplot(2,2,4),imshow(imgadap2);
title('image 2 adapted to image 1');
imwrite(imgadap2,'adapted2.jpg');
end