#!/bin/csh -fe
#
#
# Prepare an HDR image for display on new LED/LCD display
#
# August 2005	Greg Ward
#
# The $resdir directory must contain the following files:
#	descrip.txt	- Display description
#	lcdresp.dat	- LCD response data: input output
#	ledresp0.dat	- LED response data: input output
#	ledresp.dat	- LED interpolation data: # inp_min out_max
#	ledord.dat	- LED value order, increasing
#
if ( $#argv < 1 ) then
	goto userr
endif

#
# Script constants
#
# LCD display size
set dstsiz=(1280 1024)
# 2 Gaussians summed together, with following radii
set ledrad1=42
set ledcoef1=0.75
set ledrad2=140
set ledcoef2=`ev 1-$ledcoef1`
# Number of LEDs in backlight
set nleds=759
set tmpdir=/tmp/hdr$$
set origdir = "$cwd"

#
# Set display directory
#
set resdir="/Users/brian/bin/disp14"
set ledpow=.5
set ledmax=190
set sca=1
while ( "$argv[1]" =~ -* )
	switch ( $argv[1] )
	case -u:			# upper LED value
		shift argv
		set ledmax=$argv[1]
		if ($ledmax > 255 || $ledmax < 5) then
			echo "Illegal range for led max [5-255]"
			exit 1
		endif
		breaksw
	case -b:			# blackpoint insertion file
		shift argv
		set blackptf="$argv[1]"
		breaksw
	case -d:			# device directory
		shift argv
		set resdir="$argv[1]"
		if ( "$resdir" !~ /* ) then
			set resdir="$origdir/$resdir"
		endif
		if (! -d "$resdir" ) then
			echo "$resdir : no such directory"
			exit 1
		endif
		breaksw
	case -a:			# animation threshold
		shift argv
		set athresh = $argv[1]
		shift argv
		set dthresh = $argv[1]
		breaksw
	case -n:			# animation count (hack)
		shift argv
		set astep = $argv[1]
		set acount = 0
		breaksw
	case -s:			# spatial zoom factor
		shift argv
		set sca = $argv[1]
		breaksw
	case -p:			# LED exponent (fore/background)
		shift argv
		set ledpow = $argv[1]
		breaksw
	case -E:			# uncalibrated fallback
		set stabilize
	case -e:			# fixed exposure setting
		shift argv
		set exp="$argv[1]"
		if ( "$exp" =~ [+-]* ) then
			set exp=`ev "2^$exp"`
		endif
		breaksw
	case -m:			# multiplier on exposure
		shift argv
		set mul="$argv[1]"
		breaksw
	default:
		goto userr
		breaksw
	endsw
	shift argv
end
cat "$resdir/descrip.txt"
foreach i (lcdresp.dat ledresp0.dat ledresp.dat ledord.dat)
	if (! -r "$resdir/$i" ) then
		echo Missing essential file: "$resdir/$i"
		exit 1
	endif
end
if ( $?athresh ) then
	if ( $athresh < 1 || $athresh > 99 ) then
		echo "Animation threshold A must be between 1 and 99"
		exit 1
	endif
	if ( $dthresh < 1 || $dthresh > 99 ) then
		echo "Animation threshold B must be between 1 and 99"
		exit 1
	endif
	set athresh=`ev "floor($sca*$sca*$nleds*$athresh/100)"`
endif

#
# Work in a temporary directory
#
mkdir $tmpdir
onintr done
cd $tmpdir

#
# Create default black point insertion file if none given
#
if (! $?blackptf ) then
	set blackptf=blackpt.dat
	cnt $nleds | rcalc -e '$1=255' > $blackptf
endif

#
# Create diagonal matrix solver file
#
cat > contrib.cal << '_EOF_'
frac(x) : x - floor(x);
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 }
					{ neighbor values }
neigh(i) = select(i,
		gi(1,-1,1),
		gi(1,0,2),
		gi(1,1,1),
		gi(1,1,-1),
		gi(1,0,-2),
		gi(1,-1,-1)
	);
blval = gi(1) - NEIGH_FRAC*NORM_MULT*(neigh(1)+neigh(2)+neigh(3)+
					neigh(4)+neigh(5)+neigh(6));
indent_row = if(frac(y/2 + .25) - .5, 1, 0);
onbl = frac((x+1-indent_row)/2 + .25) - .5;
lo = if(onbl, blval, 0);
'_EOF_'

#
# Create response inversion file
#
set lcdmax=`sed -n 's/^255[ 	]//p' "$resdir/lcdresp.dat"`
rcalc -e '$1=$2/'$lcdmax -e '$2=$1/255' "$resdir/lcdresp.dat" \
	| tabfunc -i lcdinv > resp.cal
set ledmax2=`sed -n 's/^255[ 	]//p' "$resdir/ledresp0.dat"`
rcalc -e '$1=$2/'$ledmax2 -e '$2=$1' "$resdir/ledresp0.dat" \
	| tabfunc -i ledinv >> resp.cal
sort -n "$resdir/ledresp.dat" \
	| tabfunc ledimin ledomax >> resp.cal
cat << '_EOF_' >> resp.cal
min(a,b) : if(a-b, b, a);
ledndx(x,y) : (y*33 + x + 1)/2;
ledinv2(i,v) : if( v-ledomax(i), 255,
		if(v-.001,ledimin(i),0) + (255-ledimin(i))*v/ledomax(i) );
'_EOF_'

#
# Loop through input HDR images
#
foreach inp ($*)
if ( "$inp" !~ /* ) then
	set inp = "$origdir/$inp"
endif
if (! -r "$inp" ) then
	echo $inp cannot be read.
	continue
endif

#
# Create properly-sized HDR original
#
rm -f $tmpdir/resized.hdr
set siz=`getinfo -d < "$inp" | rcalc -i '-Y ${yr} +X ${xr}' -e '$1=xr;$2=yr'`
if ($siz[1] != $dstsiz[1] || $siz[2] != $dstsiz[2]) then
	pfilt -1 -r .6 -p 1 -x `ev "$sca*$dstsiz[1]"` \
			-y `ev "$sca*$dstsiz[2]"` "$inp" \
		| pcompos -x $dstsiz[1] -y $dstsiz[2] \=00 - \
			`ev "floor($dstsiz[1]/2)" "floor($dstsiz[2]/2)"` \
		> $tmpdir/resized.hdr
else
	ln -s "$inp" $tmpdir/resized.hdr
endif
set dest = "$inp:t"
set dest = "$origdir/$dest:r"

#
# Create "ideal" backlight image, square pixels
#
pcomb -e "lo=li(1)^$ledpow" resized.hdr \
	| pfilt -1 -b -x /12.929 -y /7.420 \
	| pfilt -1 -r 1 -x /3 -y /3 > backlight0.hdr

#
# Solve for best values on hex array with cross-talk
#
pcomb -f contrib.cal backlight0.hdr > backlight1.hdr



#
# Get maximum value to normalize output
#
set blmax=`pextrem < backlight1.hdr`
set blmax = $blmax[9]
set fgmax=`pextrem < resized.hdr`
set fgmax = $fgmax[9]
echo $fgmax
if ( $?exp ) then
	if ( $?stabilize ) then
		set blmax=`ev "$exp^-$ledpow"`
	else
		set blmax=`ev "$blmax/($exp*$fgmax)^$ledpow"`
	endif
	set fgmax=`ev "1/$exp"`
else
	if ( $?mul ) then
		set fgmax=`ev "$fgmax/$mul"`
	endif
	echo "Exposure set to -e" `ev "1/$fgmax"`
endif

#
# Output values as data file in proper order
#
pcomb -f resp.cal -e "lo=min(ledomax(ledndx(x,y)),gi(1)/$blmax)" \
		backlight1.hdr \
	> backlight1at.hdr
#
# See if we can reuse last backlight
#
unset computed_backlight
if ( $?athresh && -r backlight2.hdr ) then
	goto reuse_backlight
endif
if ( $?astep ) then
	if ( $acount <= 1 ) then
		set acount = $astep
	else
		@ acount--
		goto reuse_backlight
	endif
endif
compute_backlight:
mv backlight1at.hdr backlight1a.hdr
#
# Convert backlight drive values to 0-$ledmax using LED response
#
pvalue -h -H -pG backlight1a.hdr \
	| rcalc -f contrib.cal -e 'x=$1;y=$2;cond=onbl' -f resp.cal \
		-e "nf:$ledmax/255" \
		-e '$1=floor(nf*ledinv($3)+.5)' \
	| rlam "$resdir/ledord.dat" - $blackptf $blackptf \
	| sort -n \
	| sed -e 's/^[1-9][0-9]*$/0 255 255/' -e 's/^[1-9][0-9]*	//' \
	> led.out

#
# Compute actual backlight distribution using these values
# Multiply by final factor of 2 to compensate for blank pixels
#
pfilt -1 -c -r `ev $ledrad1/30.5` -x 112 -y 98 backlight1a.hdr \
	> backlight2a.hdr
pfilt -1 -c -r `ev $ledrad2/30.5` -x 112 -y 98 backlight1a.hdr \
	| pcomb -s $ledcoef1 backlight2a.hdr -s $ledcoef2 - \
	| pfilt -1 -e 2 -r 1 -x $dstsiz[1] -y $dstsiz[2] > backlight2.hdr



set computed_backlight
reuse_backlight:
#
# Check how many backlight pixels are incorrect
#
if (! $?computed_backlight ) then
	if ( $?athresh ) then
		# Difference old backlight with new version
		set npdiff=`pcomb -e 'clamp(x):if(x-1,1,x);pct(x,y):100*if(y,if(x-y,(x-y)/y,(y-x)/y),-x)' -e "lo=if(pct(clamp(gi(1)),clamp(gi(2)))-$dthresh,1,0)" backlight1a.hdr backlight1at.hdr | pvalue -h -H -d -pG | total | sed 's/\..*//'`
		if ( $npdiff >= $athresh ) then
			echo "Frame delta over threshold: recomputing backlight..."
			goto compute_backlight
		endif
	endif
	echo "Reusing backlight from previous frame"
endif
#
# Compute LCD front image and store as RGB PPM
# XXX color correction factors added for yellow LED cast
#
pcomb -f resp.cal -e "m=if(gi(2)-1e-10,1/$fgmax/gi(2),1)" \
		-e 'ro=lcdinv(min(m*ri(1),1)/1)' \
		-e 'go=lcdinv(min(m*gi(1),1)/1.15)' \
		-e 'bo=lcdinv(m*bi(1))' \
		resized.hdr backlight2.hdr \
	| pcompos -x $dstsiz[1] -y $dstsiz[2] - 0 0 \
 "\!pvalue -r -di -h -y 1 +x `wc -l < led.out` led.out" \
			2 `ev $dstsiz[2]-2` \
	| ra_bmp -g 1 \
	> "${dest}_hdr.bmp"

echo $dest done.
#
# Loop to next HDR input
#
end

#
# All done, clean up
#
done:
cd
rm -r $tmpdir
exit 0
userr:
echo "Usage: $0 [-d disp_dir][-u ledmax][-b blkpt][-e exp][-s sca][-p pwr] input.hdr .."
exit 1