function BitsPlusImagingColorImageTest(realimage, imagefile, plotdiffs)
% BitsPlusImagingColorImageTest([realimage][, imagefile][, plotdiffs])
%
% Test use of Color++ mode with imaging pipeline...
%
% Renders a HDR color image, converts it to Bits++ Color++ format, once via
% Matlab Bits++ toolbox function BitsPlusPackColorImage and once via GPU
% based imaging pipeline. Reads back and compares results to assess
% correctness and accuracy of the new GPU implementation.
%
% realimage = 0 (Default) Create color gradient pattern with all 65535
% possible values for a systematic test.
%
% realimage = 1 Load an example HDR image of a "real" rendered scene for
% comparison.
%
% UNFINISHED BETA CODE!
%

% Define screen
whichScreen=max(Screen('Screens'));

if nargin < 3 || isempty(plotdiffs)
    plotdiffs = 0;
end

% If realimage == 0 or left out, create synthetic color test gradient:
if nargin < 1 || isempty(realimage) || realimage==0
    % No image spec'd. Genereate a synthetic grating that covers the whole
    % intensity range:
    theImage=zeros(256,256,3);
    theImage(:,:,1)=reshape(double(linspace(0, 2^16 - 1, 2^16)), 256, 256)' / (2^16 - 1);
    theImage(:,:,2)=theImage(:,:,1);
    theImage(:,:,3)=theImage(:,:,1);
else
    % Want a real image to display. Filename spec'd?
    if nargin < 2 || isempty(imagefile)
        % Load it in.  This is an RGB image with
        % data in range 0-1.
        fprintf('Loading a high dynamic range test image\n');
        load colorPlusTest
    else
        % Load HDR image file:
        theImage = HDRRead(imagefile);
    end
end

% Normalize to 0.0 - 1.0 intensity range:
%maxpixel = max(max(max(theImage)))
%theImage = theImage / maxpixel * 10;
%theImage = theImage.^0.4545;
% Convert input image via Bits++ toolbox routines (Matlab code):
fprintf('Converting to color++ format\n');
tic;
packedImage = BitsPlusPackColorImage(theImage,0);
toc

% Show the image
[m,n,p] = size(packedImage);
rect = [0 0 m n];
fprintf('Showing image\n');

% Open a double buffered fullscreen window with black background, configured for the Bits++
% Color++ Mode, i.e., with proper setup of imaging pipeline and conversion shaders:
[window,screenRect] = BitsPlusPlus('OpenWindowColor++', whichScreen, 0);

% Find out how big the window is:
[screenWidth, screenHeight]=Screen('WindowSize', window);

% Find the color values which correspond to white and black.  Though on OS
% X we currently only support true color and thus, for scalar color
% arguments,
% black is always 0 and white 255, this rule is not true on other platforms will
% not remain true on OS X after we add other color depth modes.
white=WhiteIndex(whichScreen);
black=BlackIndex(whichScreen);
gray=(white+black)/2;
if round(gray)==white
    gray=black;
end

% draw a black background on front and back buffers
Screen('FillRect',window, black);
Screen('Flip', window);
Screen('FillRect',window, black);

% Disable Bits++ Color++ output formatter:
Screen('HookFunction', window, 'Disable', 'FinalOutputFormattingBlit');

% Show input image in LDR mode first. This will just show it with 8bpc on a
% standard display, will probably create trash on the real Bits++
textureIndex= Screen('MakeTexture',window, theImage, [], [], 2);
dstRect = Screen('Rect', textureIndex);
Screen('DrawTexture',window,textureIndex, [], dstRect);

texpacked= Screen('MakeTexture', window, packedImage);
hdrtexIndex= Screen('MakeTexture',window, theImage, [], [], 2);

the8BPCImage = uint8(theImage * 255.0 + 0.5);
ldrtexIndex= Screen('MakeTexture',window, the8BPCImage, [], [], 0);

Screen(window,'Flip');
fprintf('Hit any character to continue\n');

KbWait;
while KbCheck; end;

% Repeated toggling between Bits++ box image and GPU image, until space key
% pressed twice.
while 1
    % Disable Bits++ Color++ output formatter:
    Screen('HookFunction', window, 'Disable', 'FinalOutputFormattingBlit');

    % Draw Color++ image as generated by Bits++ Matlab-toolbox:
    dstRect = Screen('Rect', texpacked);
    Screen('DrawTexture',window,texpacked, [], dstRect, [], 0);
    Screen('FillRect', window, [255 0 0], [400 400 450 450]);
    %DrawFormattedText(window, 'Bits++ Toolbox computed, drawn via DrawTexture', 'center', 'center');
    Screen(window,'Flip');

    fprintf('Hit any character to continue\n');
    KbWait;
    while KbCheck; end;

    % Enable Bits++ Color++ output formatter:
    Screen('HookFunction', window, 'Enable', 'FinalOutputFormattingBlit');

    % Draw Color++ image as generated by PTB GPU imaging pipeline:
    dstRect = Screen('Rect', hdrtexIndex);
    Screen('DrawTexture',window,hdrtexIndex, [], dstRect, [], 0);
    %DrawFormattedText(window, 'Imaging Pipeline generated.', 'center', 'center', [], 1);
    Screen('FillRect', window, [0 255 0], [400 400 450 450]);
    % Enforce finalizing the image before we can make a screenshot:
glFinish;
tic;
    Screen('DrawingFinished', window);
glFinish;
toc

    % Take screenshot of GPU converted image:
    convImage=Screen('GetImage', window, ScaleRect(dstRect, 2, 1),'backBuffer');

    % Show it.
    Screen(window,'Flip');

    KbWait;
    [d1 d2 keycode]=KbCheck;
    while KbCheck; end;

    if keycode(KbName('space'))
        break;
    end
    
    Screen('DrawTexture',window,ldrtexIndex, [], dstRect, [], 0);
    % Show it.
    Screen(window,'Flip');

    KbWait;
    [d1 d2 keycode]=KbCheck;
    while KbCheck; end;
    
end

while ~KbCheck
    Screen('DrawTexture',window,hdrtexIndex, [], dstRect, [], 0);
    Screen(window,'Flip', GetSecs+0.02);
    Screen('DrawTexture',window,ldrtexIndex, [], dstRect, [], 0);
    Screen(window,'Flip', GetSecs+0.02);
end

Screen('CloseAll');

% Compute difference images and such...
ww = Screen('WindowSize', whichScreen);
ww = min(ww / 2, size(packedImage,2));

diffred = abs(double(packedImage(:,1:ww,1)) - double(convImage(:,1:ww,1)));
diffgreen = abs(double(packedImage(:,1:ww,2)) - double(convImage(:,1:ww,2)));
diffblue = abs(double(packedImage(:,1:ww,3)) - double(convImage(:,1:ww,3)));

% Compute maximum deviation of framebuffer raw data:
mdr = max(max(diffred));
mdg = max(max(diffgreen));
mdb = max(max(diffblue));

fprintf('Maximum raw data difference: red= %f green = %f blue = %f\n', mdr, mdg, mdb);

% If there is a difference, show plotted difference images:
if (mdr>0 || mdg>0 || mdb>0) && plotdiffs
    close all;
    imagesc(diffred);
    figure;
    imagesc(diffgreen);
    figure;
    imagesc(diffblue);
end

% Now compute a more meaningful difference: The difference between the
% stimulus as the Bits++ box would see it (i.e. how much do the 16 bit
% intensity values of each color channel differ?):
convImage = double(convImage);
packedImage = double(packedImage);

% For each color channel do...
for c=1:3
    % Invert conversion: Compute 16 bpc color values from high/low byte
    % pixel data:
    deconvImage = (zeros(size(convImage,1), size(convImage,2)/2));
    deconvImage(:,:) = 256 * convImage(:, 1:2:end-1, c) + convImage(:, 2:2:end, c);

    depackImage = (zeros(size(packedImage,1), size(packedImage,2)/2));
    depackImage(:,:) = 256 * packedImage(:, 1:2:end-1, c) + packedImage(:, 2:2:end, c);

    % Difference image:
    diffImage = (deconvImage - depackImage);

    % Find locations where pixels differ:
    idxdiff = find(abs(diffImage) > 0);
    numdiff(c) = length(idxdiff);
    numtot(c) = size(diffImage,1)*size(diffImage,2);
    maxdiff(c) = max(max(abs(diffImage)));
    [row col] = ind2sub(size(diffImage), idxdiff);

    % Print out all pixels values which differ, and their difference:
    if plotdiffs
        for j=1:length(row)
            fprintf('Diff: %.2f Input Value: %.20f\n', diffImage(row(j), col(j)), theImage(row(j), col(j), c) * 65535);
        end
    end
end

for c=1:3
    % Summarize for this color channel:
    fprintf('Channel %i: %i out of %i pixels differ. The maximum absolute difference is %i.\n', c, numdiff(c), numtot(c), maxdiff(c));
end

% Done for now.
return;


% if the system only has one screen, set the LUT in Bits++ to a linear ramp
% if the system has two or more screens, then blank the screen.
if (whichScreen == 0)
    % restore the Bits++ LUT to a linear ramp
    linear_lut =  repmat(round(linspace(0, 2^16 -1, 256))', 1, 3);
    BitsPlusSetClut(window,linear_lut);

    % draw a gray background on front and back buffers to clear out any old LUTs
    Screen('FillRect',window, gray);
    Screen('Flip', window);
    Screen('FillRect',window, gray);
    Screen('Flip', window);

    % Close the window.
    Screen('CloseAll');
else
    % Blank the screen
    BitsPlusSetClut(window,zeros(256,3));

    % draw a black background on front and back buffers to clear out any old LUTs
    Screen('FillRect',window, black);
    Screen('Flip', window);
    Screen('FillRect',window, black);
    Screen('Flip', window);

    Screen('CloseAll');
end