% hdrEnc
%
% Takes in an image array and converts it to a frame-buffer ready to be
% sent to the Brightside HDR display. for the passed-in image data.
% Duplicates the functionality of the disp14enc script from Brightside as
% closely as possible.  
%
% Also contains a convex optimization solver that results in more accurate 
% results at the cost of 10x runtime (CVX toolbox needed for convex optimization) 
%
% Usage:
%
% result_image = hdrEnc( images, OPTIONS )
%
% Required arguments:
%   image: An color image/image stack intended for the Brightside HDR display, with size
%          height x width x 3 x (# images).  
%
% Optional arguments:
%   Specify optional arguments using '<argument_name>',<argument_value> pairs.  
%   Possible argument names and the values expected for the argument are:
%
%  led_max:  Scalar value in the range [5,255] signifying the maximum
%            allowed LED brightness.  Defaults to 190. (-u in disp14enc)
%  blackpoint: Either a vector containing n_leds values, or the the black
%              point insertion file name.  Defaults to all 255 (-b)
%  dev_dir:  Directory from which to read device parameter files.  Defaults
%            to 'disp14' (-d)
%  anim_thres: A vector of size two, listing the two thresholds for
%            recalculating the backlight when encoding a sequence of frames.  The
%            second value is the percentage change threshold for individual LEDs, and
%            the first is the number of LEDs that must be over threshold before the
%            backlight will be recalculated. (-a)
%  anim_count: The animation count (hack) (-n)
%  scale:  Image magnification factor. (-s)
%  led_exp:  Exponent for approximating LED response curve.  Defaults to
%            0.5. (-p)
%  exp_auto: Set to '1' to use an uncalibrated fallback exposure.  (-E)
%  exp:  Set a fixed exposure setting.  Use a scalar for absolute exposure,
%         and a string like '+5' or '-5' to use powers-of-two exposures
%         (f-stops).  (-e)
%  exp_mult: Exposure multiplier (-m)
%  cvx: Use convex optimization solver instead of simple crosstalk solver.
%       1 - Turn on convex optimization solver (requires the CVX toolbox)
%       2 - Reuse last cvx solution (use only if the input image the last
%       time the convex optimization solver was used is the same as now)

function result_image = hdrEnc(images, varargin)

t_start = clock;

%% Script constants

dstsiz=[1024 1280];    % Display panel resolution

ledrad1  = 42;         % These define the PSF of the LEDs, as a sum of 
ledcoef1 = 0.75;       % two gaussians.  Radiai are in pixels.
ledrad2  = 140;
ledcoef2 = 1-ledcoef1;

nleds=759;

resdir = 'disp14';          % Default device data directory
ledpow=.5;             % LED gamma exponent
ledmax=190;            % Maximum allowed LED drive value
sca=1;                 % Default magnification factor
stabilize = 0;         % Autoexposure setting

use_cvx  = 0;           % Use basic crosstalk solver, not convex optimization

%% Read in optional arguments

if mod(length(varargin),2) == 1
   disp 'Optional argument missing value'
   help hdrSimulate
   return
end

for i=1:length(varargin)/2
   argname = lower(varargin{i*2-1});
   value = varargin{i*2};
   switch argname
       case 'led_max' % Maximum LED drive value
           if (ledmax > 255 || ledmax < 5)
               disp 'Illegal range for led_max [5-255]'
               return
           end
           ledmax = value;
       case 'blackpoint' % Not sure what this is for, actually
           if (ischar(value) )
               blackptf = load(value);
           else
               blackptf = value;
           end
       case 'dev_dir' % Device parameter directory
           if ( ~exist(value,'dir') )
               disp(sprintf('%s : no such directory',value))
               return
           end
           resdir = value;
       case 'anim_thres' % Thresholds for backlight recalculation
           athresh = value(1);
           dthresh = value(2);
       case 'anim_count' 
           astep = value;  % I don't think these do anything, but
           acount = 0;     % I copied their implementation anyway
       case 'scale' % Magnification factor
           sca = value;
       case 'led_exp' % LED gamma exponent
           ledpow = value;
       case 'exp_auto' % Turn on/off exposure stabilization 
           stabilize = value;
       case 'exp' % Set fixed exposure instead of automatically selecting one
           if (ischar(value))
               exposure = 2^eval(value);
           else
               exposure = value;
           end
       case 'exp_mult' % Multiply exposure by a fixed value
           exp_mul=value;
       case 'cvx' % Convex optimizer options
           use_cvx = value;
       otherwise
           disp 'Unknown option'
           help hdrEnc
           return
   end
end

%% Check for existance of important files

type([resdir '/descrip.txt'])
keyfiles = { 'lcdresp.dat', 'ledresp0.dat', 'ledresp.dat', 'ledord.dat' };
for i=1:length(keyfiles)
    if ~exist([resdir '/' keyfiles{i}],'file')
        disp(['Missing essential file ' resdir '/' keyfiles{i}])
        return
    end
end

%% Set up animation constants

if exist('athresh','var')
    if athresh < 1 || athresh > 99
        disp 'Animation threshold A must be between 1 and 99'
        return
    end
    if dthresh < 1 || dthresh > 99
        disp 'Animation threshold B must be between 1 and 99'
        return
    end
    athresh = floor(sca*sca*nleds*athresh/100);
end

% Create default black point insertion file if none given
%   (What is the black point insertion file? I don't know, but its values
%   are output as the green and blue channels at the LED driving pixels,
%   where red is the LED drive value)

if (~exist('blackptf','var')) 
    blackptf = ones(nleds,1)*255;
end

%% Create response inversion lookups

% Load row-order -> frame buffer order map
ledord=load([resdir '/ledord.dat']);

% Load LCD response curve, calculate normalized inverse curve
lcdresp=load([resdir '/lcdresp.dat']);
lcdmax = lcdresp(find(lcdresp(:,1)==255),2);
lcdinv=zeros(size(lcdresp));
lcdinv(:,1) = lcdresp(:,2)/lcdmax;
lcdinv(:,2) = lcdresp(:,1)/255;
% Load LED response curve, calculate normalized inverse curve
ledresp0=load([resdir '/ledresp0.dat']);
ledmax2 = ledresp0(find(ledresp0(:,1)==255),2);
ledinv=zeros(size(ledresp0));
ledinv(:,1) = ledresp0(:,2)/ledmax2;
ledinv(:,2) = ledresp0(:,1);
% Load LED maximum/minimum values
ledresp=load([resdir '/ledresp.dat']);
ledresp = sortrows(ledresp,1);
ledimin = [ledresp(:,1) ledresp(:,2)];
ledomax = [ledresp(:,1) ledresp(:,3)];

% Define backlight image pixel -> LED index mapping function
ledndx = inline( '(y*33 + x + 1)/2' , 'x' ,'y');
% ledinv2 defined at the end of the file

%% Loop though input HDR images

for index = 1:size(images,4)

    image = images(:,:,:,index);

%% Create properly-sized HDR original

    resized = hdrResize(image);
    
%% Create "ideal" backlight image, square pixels
    
    % Ideal backlight image is genererated by two consecutive
    % filtering/downsampling steps, first with an averaging (box) filter,
    % then with a gaussian filter
    box_filter_size = [7 12];  % disp14enc may use a non-integer-sized filter here
    intermediate_siz = [138 99]; % First downsampling destination size
    gaussian_radius = 2.1; % disp14enc uses 3, but 2.1 seems to match what disp14enc produces
    backlight_siz =[46 33]; % Size of backlight array if it were filled in to become rectangular
    
    backlight0 = im2bright(resized) .^ ledpow;
    backlight0 = imfilter(backlight0,fspecial('average',box_filter_size),'replicate'); 
    backlight0 = imresize(backlight0,intermediate_siz,'bilinear',0);
    backlight0 = imfilter(backlight0,fspecial('gaussian',15,gaussian_radius),'replicate'); 
    backlight0 = imresize(backlight0,backlight_siz,'bilinear',0);

    switch use_cvx
%% Solve for best values on hex array with cross-talk
        case 0
            backlight1 = contrib(backlight0);  % contrib defined near end of file
%% Solve using convex optimizer            
        case 1
            disp 'Starting convex optimizer'
            y = backlight0(:);
            load('backlightMap.mat');
            cvx_begin
               variable x(759,1)
               minimize (norm(y-backlightMap*x,2))
               subject to
                 x>=zeros(759,1);
            cvx_end
            save 'hdrEnc_lastCVXsoln.mat' x
            backlight1 = x;
            disp 'Convex optimizer completed'
%% Solve by loading old convex optimizer solution
        case 2
            disp 'Reusing last convex optimizer solution.'
            load('hdrEnc_lastCVXsoln.mat');
            load('backlightMap.mat');
            backlight1 = x;
    end
    
%% Get maximum value to normalize output

    blmax = max(backlight1(:));         % Maximum backlight brightness
    fgmax = max(max(resized(:,:,2)));   % Maximum green channel value from source
    
    % Calculate exposure normalizations
    if exist('exposure','var')
        if (stabilize == 1)
            blmax = exposure^(-ledpow);
        else
            blmax = blmax/(exposure*fgmax)^ledpow;
        end
        fgmax = 1/exposure;
    else
        if exist('exp_mul','var')
            fgmax = fgmax/exp_mul;
        end
        disp(sprintf('Exposure set to -e %0.5g',1/fgmax));
    end

    switch use_cvx

%% ------------------
%% Crosstalk solver backlight operations
        case 0
%% Output values in proper order        
            x=repmat(0:size(backlight1,2)-1, size(backlight1,1),1);
            y=repmat( (0:size(backlight1,1)-1)', 1, size(backlight1,2));
            backlight1at = min(interp1(ledomax(:,1),ledomax(:,2),ledndx(x,y)),backlight1/blmax);
   
%% See if we can reuse last backlight
            ok_to_reuse = 0;
            if exist('athresh','var') && exist('backlight2','var')
                ok_to_reuse = reuse_backlight(backlight1a, backlight1at, athresh, dthresh);
            elseif exist('astep','var')
                if acount <= 1
                    acount = astep;
                else
                    acount = acount - 1;
                    ok_to_reuse = reuse_backlight(backlight1a, backlight1at);
                end
            end
   
            if (ok_to_reuse == 0)
%% Convert backlight drive values to 0-ledmax using LED response

                backlight1a = backlight1at';
                onbl_val = onbl(backlight1a);
                led_out_t = backlight1a(onbl_val > 0 );   % Only grab values that are actual LEDs in the hex pattern
                led_out_t = floor(ledmax/255*interp1(ledinv(:,1),ledinv(:,2),led_out_t)+0.5);  % Use inverse LED response to map to framebuffer values
                led_out = repmat([0 255 255],length(ledord),1);
                led_out(ledord(1:nleds),:) = [led_out_t blackptf blackptf]; % Reorganize LED drive values

%% Computer actual backlight distribution using these values.  Multiply by
%% final factor of 2 to compensate for blank pixels.

                 % Filter constants
                 backlight_siz = [98 112]; % Intermediate size for backlight distribution creation
                 scalef = dstsiz(2)/backlight_siz(2); % LED PSF defined at full resolution, calculate scaling factor for filtering at lower res
                 fudge_factor = 0.75; % Empirical value to match disp14enc and Matlab output
                 gaussrad1 = ledrad1/scalef*fudge_factor; % Width of 1st gaussian filter, adjusted for scale and filter differences
                 gaussrad2 = ledrad2/scalef*fudge_factor; % Width of 2nd gaussian filter, adjusted for scale and filter differences
                 gausswidth1 = round(ledrad1/scalef*3)*2+3; % Width of filtering window, just needs to be big enough not to cut off gaussian
                 gausswidth2 = round(ledrad2/scalef*3)*2+3; % Width of filtering window, just needs to be big enough not to cut off gaussian
                 
                 % Resize point light image to filtering size
                 backlight2a = imresize(backlight1a',backlight_siz,'bicubic');
                 % Filter with the two gaussians used to define LED PSF
                 backlight2a_gauss1 = imfilter(backlight2a,fspecial('gaussian',gausswidth1,gaussrad1),'replicate');
                 backlight2a_gauss2 = imfilter(backlight2a,fspecial('gaussian',gausswidth2,gaussrad2),'replicate');
                 % Combine gaussians
                 backlight2a = backlight2a_gauss1 * ledcoef1 + ...
                     backlight2a_gauss2 * ledcoef2;
                 % Resize backlight to full resolution
                 backlight2 = 2*imresize(backlight2a,dstsiz,'bicubic');
                 
            end

%% ------------------
%% Convex optimization solver backlight operations

        case {1,2}

%% Output values in proper order

           backlight1at = min(ledomax(:,2),backlight1/blmax);

%% Convert backlight drive values to 0-ledmax using LED response

           led_out_t = backlight1at;   
           led_out_t = floor(ledmax/255*interp1(ledinv(:,1),ledinv(:,2),led_out_t)+0.5);  % Use inverse LED response to map to framebuffer values
           led_out = repmat([0 255 255],length(ledord),1);
           led_out(ledord(1:nleds),:) = [led_out_t blackptf blackptf]; % Reorganize LED drive values

%% Computer actual backlight distribution using these values.
           backlight1a = reshape(backlightMap*x,46,33);

           backlight2 = imresize(backlight1a,dstsiz,'bicubic');
           backlight2 = backlight2 / max(backlight2(:));
    end
%% ------------------
    
%% Computer LCD front image. XXX color correction factors added for yellow LED cast

    % Compute correction factors
    m = zeros(dstsiz);
    m(backlight2 > 1e-10) = 1/fgmax./backlight2(backlight2 > 1e-10);
    m(backlight2 <= 1e-10) = 1;
    % LED tint constants
    red_tint_correction = 1;
    green_tint_correction = 1.15;
    blue_tint_correction = 1;
    % LED maximum cutoffs
    red_cutoff = 1;
    green_cutoff = 1;
    blue_cutoff = (use_cvx == 0) * 1e6 + 1;  % Huge for crosstalk, 1 for CVX
    % Apply tint, LCD response curve
    %   (Why does the blue channel not have a min() with it? I don't know, but
    %    that's how disp14enc has it)
    result_image =        interp1(lcdinv(:,1), lcdinv(:,2),  min(m.*resized(:,:,1),red_cutoff)/red_tint_correction    ,[],'extrap');
    result_image(:,:,2) = interp1(lcdinv(:,1), lcdinv(:,2),  min(m.*resized(:,:,2),green_cutoff)/green_tint_correction,[],'extrap');
    result_image(:,:,3) = interp1(lcdinv(:,1), lcdinv(:,2),  min(m.*resized(:,:,3),blue_cutoff)/blue_tint_correction  ,[],'extrap');
    % Normalize on 0-255
    result_image = min(max(result_image,0),1)*255; 
    % Composite LED values onto the LCD image
    result_image(2,3:2+size(led_out,1),:) = reshape(led_out,1,size(led_out,1),size(led_out,2));
    % Convert to uint8
    result_image = uint8(result_image);

    if ~exist('result_images','var')
        result_images = result_image;
    else
        result_images(:,:,:,index) = result_image;
    end
end

%% Done - display elapsed time

t_end = clock;

disp(sprintf('Done - elapsed time: %.2f seconds',etime(t_end,t_start)));

end

%% Check how many backlight pixels are incorrect
function ok_to_reuse = reuse_backlight(old, new, athresh, dthresh)
    if exist('athresh','var')
        npdiff = sum(100*abs(old' - new)./new > dthresh);
        if npdiff>=athresh
            disp 'Frame delta over threshold: recomputing backlight...'
            ok_to_reuse = 0;
            return;
        end
    end
    disp 'Reusing backlight from previous frame'
    ok_to_reuse = 1;
    
end

%% Utility functions
% This is defined by not actually used by disp14enc.  Included for
% completeness
function val = ledinv2(v,ledomax,ledimin)
    if v > ledomax
        val = 255;
    else
        if v > 0.001
            val = ledimin;
        else 
            val = 0;
        end
        val = val + (255 - ledimin)*v/ledomax;
    end
end

% This function returns a matrix that indicates whether a rectangular grid
% position actually contains an LED, which are laid out on a hex grid.
function val = onbl(image)
    
    frac = inline('x-floor(x)','x');
    x = repmat(0:(size(image,2)-1), size(image,1), 1);
    y = repmat((0:(size(image,1)-1))', 1, size(image,2));
    
    indent_row = ( (frac(y/2 + 0.25) - 0.5) > 0);
    val = -(frac((x+1-indent_row)/2 + 0.25) - 0.5); % neg.sign seems backwards, but output matches script this way.
end

%% contrib.cal equivalent - diagonal matrix solve
% Adjusts LED drive values to account for crosstalk from neighboring LEDs
function hexbacklight = contrib(image)
   
    NEIGH_FRAC = 0.712 / 3.07;    % average neighbor contribution
    NEIGH_CONT = 6 * NEIGH_FRAC;  % approximate neighborhood total
    NORM_MULT = 1/(1+NEIGH_CONT); % uniform normalization factor
   
    % Crosstalk matrix for convolution
    conv_kernel = [         0              -NEIGH_FRAC*NORM_MULT            0; ...
                   -NEIGH_FRAC*NORM_MULT            0              -NEIGH_FRAC*NORM_MULT; ...
                            0                       1                       0; ...
                   -NEIGH_FRAC*NORM_MULT            0              -NEIGH_FRAC*NORM_MULT; ...
                            0              -NEIGH_FRAC*NORM_MULT            0];

    % disp14enc assumes image edge values continue to infinity, so pad for the
    % convolution
    im_pad = [image(1,1) image(1,1) image(1,:) image(1,end) image(1,end); ...
              image(1,1) image(1,1) image(1,:) image(1,end) image(1,end); ...
              image(:,1) image(:,1) image      image(:,end) image(:,end); ...
              image(end,1) image(end,1) image(end,:) image(end,end) image(end,end); ...
              image(end,1) image(end,1) image(end,:) image(end,end) image(end,end)];
    % Convolve, and clamp negative values to 0 
    blval = max(conv2(im_pad, conv_kernel,'same'),0);
    % Remove padding regions
    blval = blval(3:end-2,3:end-2);
    % Only keep values at actual LED positions
    hexbacklight = (onbl(image) > 0) .* blval;  
        
end

%% Color->luminance image converter, based on Radiance source code.
function brightimage = im2bright(color_image)
    % Following definitions are from the Radiance ray-tracer source code,
    % /src/common/color.h
    
    CIE_x_r = 0.640;
    CIE_y_r	= 0.330;
    CIE_x_g	= 0.290;
    CIE_y_g	= 0.600;
    CIE_x_b	= 0.150;
    CIE_y_b	= 0.060;
    CIE_x_w	= 0.3333;
    CIE_y_w	= 0.3333;

    STDPRIMS =	[CIE_x_r,CIE_y_r;...
                 CIE_x_g,CIE_y_g;...
 				CIE_x_b,CIE_y_b];

    CIE_D	  =	(CIE_x_r*(CIE_y_g - CIE_y_b) + ... 
				CIE_x_g*(CIE_y_b - CIE_y_r) +  ...
				CIE_x_b*(CIE_y_r - CIE_y_g)	);
    CIE_C_rD  = ( (1./CIE_y_w) * ...
				( CIE_x_w*(CIE_y_g - CIE_y_b) - ...
				  CIE_y_w*(CIE_x_g - CIE_x_b) + ...
				  CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g	) );
              
    CIE_C_gD  = ( (1./CIE_y_w) * ...
				( CIE_x_w*(CIE_y_b - CIE_y_r) - ...
				  CIE_y_w*(CIE_x_b - CIE_x_r) - ...
				  CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r	) );
              
    CIE_C_bD  = ( (1./CIE_y_w) * ...
        ( CIE_x_w*(CIE_y_r - CIE_y_g) - ...
        CIE_y_w*(CIE_x_r - CIE_x_g) + ...
        CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r	) );

    CIE_rf    = (CIE_y_r*CIE_C_rD/CIE_D);
    CIE_gf    = (CIE_y_g*CIE_C_gD/CIE_D);
    CIE_bf    = (CIE_y_b*CIE_C_bD/CIE_D);

    % As of 9-94, CIE_rf=.265074126, CIE_gf=.670114631 and CIE_bf=.064811243

    % The following definitions are valid for RGB colors only...

    brightimage	= max(CIE_rf*color_image(:,:,1)+CIE_gf*color_image(:,:,2)+CIE_bf*color_image(:,:,3),0);

end

       