## This file is a part of the supplemental metarial for the paper:
##
## Perception-motivated High Dynamic Range Video Encoding
## Rafal Mantiuk, Grzegorz Krawczyk, Karol Myszkowski, Hans-Peter Seidel
##
## @file percept_quant_lum.m
##
## Find the perceptual quantized luminance space using Ferwerda's t.v.i.
## function (refer to Section 2.1 in the paper).
##
## The script stores a look up table of the mapping in vector psi and a
## quality coefficient f in scalar f . The mapping can be plotted using:
## plot( l, log10(psi))
##
## The code can be run as a GNU Octave2.0 script
## (http://www.octave.org/).

1; # This is a script file

function T = tvi_ferwerda( L_adapt )
  ## Threshold versus Intensity function
  ## Proposed by Ferwerda'96
  ## L_adapt - adapting luminance
  ## T - absolute perceivable threshold
  logL_adapt = log10( L_adapt );
  logT = \
      (logL_adapt < -3.94) * -2.86 + \
      (logL_adapt >= -3.94) .* (logL_adapt < -1.44) .* \
         ((0.405*logL_adapt+1.6).^2.18-2.86) + \
      (logL_adapt >= -1.44) .* (logL_adapt < -0.0184) .* \
         (logL_adapt - 0.395) + \
      (logL_adapt >= -0.0184) .* (logL_adapt < 1.9) .* \
         ((0.249*logL_adapt+0.65).^2.7 - 0.72) + \
      (logL_adapt >= 1.9) .* (logL_adapt - 1.255);
  T = 10.^logT;
endfunction

global f;
 
function [f_out psi] = tpbp_shooting( thresholdFunc, l, y_min, y_max )
  ## Solve Two Point Boundary Problem using shooting method
  ## thresholdFunc - name of the function that returns
  ##     a right hand side of the differential equation
  ## l - vector that defines points in Lp space
  ## y_min - value of l(0); left boundary
  ## y_max - value of 1(last); right boundary
  ## f_out - a parameter in thresholdFunc
  ## psi - ODE solution; y values that fulfills equation
  global f;
  fl = 1;
  fr = length( l );
  while( fr-fl > 0.001 )        # Binary search for solution that
                                # fulfills boundary conditions
    f = (fl+fr)/2;
    psi = lsode( thresholdFunc, y_min, l ); # Solve ODE
    if( psi(length(l)) > y_max ) 
      fl = f;
    else
      fr = f;
    endif
  endwhile
  f_out = f;
endfunction

function xdot = tvi_psi( y, l )
  ## Right hand side of the differential equation
  global f;
  xdot = 2*tvi_ferwerda(y)/f;
endfunction

bits=11
l=1:2^bits;

## Find a perceptually quantized mapping function for
## 11-bit integer space 'l'
[f psi] = tpbp_shooting( "tvi_psi", l, 10^-4, 10^8 );