% hdrSimulate
%
% Creates a Matlab visualization of the image that would appear on a Brightside HDR display
% for the passed-in image data.  Output values are measured in cd/m^2.
%
% Usage:
%
% result_image = hdrSimulate( image, OPTIONS)
%
% Required arguments:
%   image: An 0-255 integer-valued 3-channel framebuffer image intended for the Brightside HDR display, with size
%          width x height x 3
%   
% Optional arguments:
%   Specify optional arguments using '<argument_name>',<argument_value> pairs.  
%   Possible argument names and the values expected for the argument are:
%
%   led_x: The x positions of the Brightside display backlight LEDs. 
%          If not specified, loaded from file ledpos.mat
%   led_y: The y positions of the Brightside display backlight LEDs.
%          If not specified, loaded from file ledpos.mat
%   led_ord: The mapping from control value location to LED index. 
%            If not specified, loaded from file ledpos.mat
%   led_pos: Specify a different file to read LED x,y,and ordering data from.  
%            If not specified, ledpos.mat is used.
%   lcd_resp: The LCD response curve, mapping input 0-255 control values 
%             to output 0-255 intensity values.
%             This should be a n by 2 matrix, which will be linearly 
%             interpolated to come up with a full mapping.
%             Or you can pass a filename from which to read the curve from.
%             If not specified, an empirically fit curve is used.
%   led_resp: The LED response curve, mapping input 0-255 control values to 
%             output 0-255 intensity values.
%             This should be a n by 2 matrix, which will be linearly 
%             interpolated to come up with a full mapping.
%             It can also be a 3D array, with the the ith (:,:,i) plane being
%             the i:th LED's response curve.
%             Or it can be a filename from which to read the curve from.
%             If not specified, an empirically fit curve is used.  
%   led_psf: The point-spread function of a single backlight LED.  
%            If not specified, assumed to be a sum of two gaussians
%            This should either be a vector mapping distance from LED center 
%            in pixels to screen intensity at that distance.
%            or a 2D image of the PSF, with the LED center at the image center.
%            or a filename from which the PSF can be loaded from.
%
% For example, to specify a different led_x vector, use:
%
% result_image = hdrSimulate(foo_image, 'led_x', new_led_x);

function result_image = hdrSimulate( image, varargin )

%% Check validity of input image

if nargin < 1
   disp 'Required arguments missing.  Usage:'
   help hdrSimulate
   return
end

if length(size(image)) ~= 3 | size(image,3) ~= 3
   disp 'Input image needs to contain three color planes.  Usage:'
   help hdrSimulate
   return
end

%% Set up default values

w = size(image,2);
h = size(image,1);

ledpos_file = 'ledpos.mat';
lcd_resp = [[0:10:255 255]',217*(([0:10:255 255]'/255).^2)+0.5];  % Matches sample LCD reasonably well
led_resp = [[0:10:255 255]',max(98*(exp(1)-exp(1-[0:10:255 255]'/255))-11,0)];  % Matches sample LED reasonably well.   

% PSF is the sum of two gaussians
psf_radius = 500;  % Increasing this too much will drastically increase run time.
ledrad1  = 42;
ledcoef1 = 0.75;
ledrad2  = 140;
ledcoef2 = 1-ledcoef1;
led_psf = exp(-[0:psf_radius].^2/(2*ledrad1^2))*ledcoef1 + exp(-[0:psf_radius].^2/(2*ledrad2^2))*ledcoef2; % Radial slice
led_psf = led_psf / sum(led_psf);  % Normalization

%% Process 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_x'
      led_x = value;
   case 'led_y'
      led_y = value;
   case 'led_ord'
      led_ord = value;
   case 'led_pos'
      ledpos_file = value;
   case 'lcd_resp'
      lcd_resp = value;
      if (ischar(lcd_resp))
         lcd_resp=load(lcd_resp);
      end
      if (size(lcd_resp,2) ~= 2)
         disp 'LCD response should be a N x 2 matrix'
         return
      end
      if (min(lcd_resp(:,1)) ~= 0 | max(lcd_resp(:,1)) ~= 255)
         disp 'Warning: LCD response curve domain isn''t [0,255]'
      end
   case 'led_resp'
      led_resp = value;
      if (ischar(led_resp))
         led_resp=load(led_resp);
      end
      if (size(led_resp,2) ~= 2)
         disp 'LED response should be a N x 2 (x n_leds) matrix'
         return
      end
      if (ndims(led_resp) > 2)
          disp 'Using a unique LED response curve for each LED'
      else
          if (min(led_resp(:,1)) ~= 0 | max(led_resp(:,1)) ~= 255)
              disp 'Warning: LED response curve domain isn''t [0,255]'
          end
      end
   case 'led_psf'
      led_psf = value;
      if ischar(led_psf)
         led_psf = load(led_psf);
      end
      psf_radius = length(led_psf) - 1;
   case 'dead_leds'
      dead_leds = value;
   otherwise
      disp(sprintf('Unknown argument %s',argname))
   end
end

% Read in LED position values if needed

if ~exist('led_x','var')
   ledpos = load(ledpos_file);
   led_x = ledpos.led_x;
end

if ~exist('led_y','var')
   if ~exist('ledpos','var'); ledpos = load(ledpos_file); end
   led_y = ledpos.led_y;
end

if ~exist('led_ord','var')
   if ~exist('ledpos','var'); ledpos = load(ledpos_file); end
   led_ord = ledpos.ledord;
end

n_leds = length(led_x);
if n_leds ~= length(led_y) | n_leds ~= length(led_ord)
   disp 'led_x, led_y, and led_ord do not have matching length.'
   return
end

if size(led_resp,3) > 1
    if size(led_resp,3) ~= n_leds
        disp 'More than one LED respose curve specified, but not the right number'
        disp(sprintf('Have %d curves, need %d curves',size(led_resp,3),n_leds))
        return
    end
end

%% Convert image to double  

image = double(image);

%% Extract LED control values.  

led_control = image(2,2+led_ord(1:n_leds),1); % Red channel contains the drive value.

%% Mimic HDR display intensity curves

% Convert LCD/LED response arrays to lookup tables.  

lcd_lookup = interp1(lcd_resp(:,1),lcd_resp(:,2),[0:255]);
led_lookup = interp1(led_resp(:,1),led_resp(:,2),[0:255]);

% Convert control values to output intensity

led_drive = led_lookup(led_control+1);  
lcd_drive = reshape(lcd_lookup(image(:)+1),h,w,3);

%% Create backlight output light distribution

% Create 2D PSF of radial LED PSF

if length(led_psf) ~= length(led_psf(:))   % Test for 2D-ness
   led_psf_matrix = led_psf       % Passed in a 2D PSF already, use it
else
   % This is much faster than the straightforward for-loop
   x_map = repmat([0:psf_radius*2],1,psf_radius*2+1);
   y_map = repmat([0:psf_radius*2],psf_radius*2+1,1);
   y_map = y_map(:)';
   psf_val = interp1(led_psf,min(sqrt( (x_map-psf_radius) .^ 2 + (y_map - psf_radius) .^ 2)+1,psf_radius));
   led_psf_matrix = reshape(psf_val,psf_radius*2+1,psf_radius*2+1);
end

% Create backlight LED image as point lights

result_image = zeros(h,w);  
for i=1:n_leds
	result_image(led_y(i)+1,led_x(i)+1) = led_drive(i);
end

% Convolve with PSF to create real light distribution patterns

result_image = conv2(result_image, led_psf_matrix, 'same');

%% Convert backlight grayscale intensity to full-color - apply yellow tint. 

% Green should be 1.15 times red and blue, based on disp14enc script.  But keep overall output power the same

red_tint = (1/1.15)^(1/3);
green_tint = 1.15 * red_tint;
blue_tint = red_tint;  % red*green*blue = 1, green = 1.15*blue and red

result_image(:,:,2) = result_image(:,:) * green_tint;    
result_image(:,:,3) = result_image(:,:,1) * blue_tint;
result_image(:,:,1) = result_image(:,:,1) * red_tint; 

%% Modulate backlight with LCD image

result_image = result_image(1:h,1:w,:) .* lcd_drive;  % Need to crop result_image, due to LED positions outside the LCD

%% Normalize so that the brightness at the center of a single LED set to 150 is 490 cd/m^2

white_luminance = 256;  % cd/m^2
white_LED_drive = 150;  % LED drive which created the white_luminance value
max_LED_drive = 255;
max_white_luminance = white_luminance * (led_lookup(max_LED_drive)/led_lookup(white_LED_drive));
max_single_led_drive = max(led_psf) * max(lcd_lookup) * max(led_lookup); % Normalization factor 

% Assumes all color channels have equal luminance output at max.  Very
% likely not true.  
max_red_luminance = white_luminance / 3;
max_green_luminance = white_luminance / 3;
max_blue_luminance = white_luminance / 3;

result_image(:,:,1) = result_image(:,:,1) ./ max_single_led_drive * max_red_luminance;  % Result is in cd/m^2
result_image(:,:,2) = result_image(:,:,2) ./ max_single_led_drive * max_green_luminance;  % Result is in cd/m^2
result_image(:,:,3) = result_image(:,:,3) ./ max_single_led_drive * max_blue_luminance;  % Result is in cd/m^2

%% Done

return