orig/psquish.csh
author Axel Jacobs <axel@jaloxa.eu>
Mon, 14 Apr 2014 22:21:17 +0100
changeset 73 87fa3132a136
parent 0 0aa115157c9c
permissions -rw-r--r--
New IES test file with spaces in file name
     1 #!/bin/csh -f
     2 # RCSid: $Id: psquish.csh,v 3.5 2008/11/10 19:08:19 greg Exp $
     3 set Lmin=.0001		# minimum visible world luminance
     4 set Ldmin=1		# minimum display luminance
     5 set Ldmax=100		# maximum display luminance
     6 set nsteps=100		# number of steps in perceptual histogram
     7 set cvratio=0.05	# fraction of pixels to ignore in envelope clipping
     8 set td=/tmp
     9 set tf0=$td/tf$$
    10 set tf1=$td/hist$$
    11 set tf1b=$td/hist$$.diff
    12 set tf2=$td/cumt$$
    13 set tf3=$td/histeq$$.cal
    14 set tf4=$td/cf$$.cal
    15 set tf=($tf0 $tf1 $tf1b $tf2 $tf3 $tf4)
    16 if ( "$argv[1]" == "-a" ) then
    17 	set adaptive
    18 	shift argv
    19 endif
    20 if ( $#argv != 1 ) then
    21 	echo "Usage: $0 [-a] input.hdr > output.hdr"
    22 	exit 1
    23 endif
    24 set ifile=$1
    25 set ibase=$ifile:t
    26 if ( "$ibase" =~ *.hdr ) set ibase=$ibase:r
    27 set ibase=$ibase:t
    28 onintr quit
    29 pextrem -o $ifile > $tf0
    30 set Lmin=`rcalc -e 'L=179*(.265*$3+.67*$4+.065*$5)' -e 'cond=1.5-recno;$1=if('L-$Lmin,L,$Lmin')' $tf0`
    31 set Lmax=`rcalc -e 'L=179*(.265*$3+.67*$4+.065*$5)' -e 'cond=recno-1.5;$1=if('$Ldmax-L,$Ldmax,L')' $tf0`
    32 cat > $tf3 << _EOF_
    33 min(a,b) : if(a-b, b, a);
    34 WE : 179;			{ Radiance white luminous efficacy }
    35 Lmin : $Lmin ;			{ minimum visible luminance }
    36 Lmax : $Lmax ;			{ maximum picture luminance }
    37 Ldmin : $Ldmin ;		{ minimum output luminance }
    38 Ldmax : $Ldmax ;		{ maximum output luminance }
    39 Stepsiz : (Bl(Lmax)-Bl(Lmin))/ $nsteps ;	{ brightness step size }
    40 			{ Logarithmic brightness function }
    41 Bl(L) : log(L);
    42 Lb(B) : exp(B);
    43 BLw(Lw) : Bl(Ldmin) + (Bl(Ldmax)-Bl(Ldmin))*cf(Bl(Lw));
    44 			{ first derivative functions }
    45 Bl1(L) : 1/L;
    46 Lb1(B) : exp(B);
    47 			{ histogram equalization function }
    48 lin = li(1);
    49 Lw = WE/le(1) * lin;
    50 Lout = Lb(BLw(Lw));
    51 mult = if(Lw-Lmin, (Lout-Ldmin)/(Ldmax-Ldmin)/lin, 0) ;
    52 _EOF_
    53 if ( $?adaptive ) then
    54 	cat >> $tf3 << _EOF_
    55 			{ Ferwerda contrast sensitivity function }
    56 				{ log10 of cone threshold luminance }
    57 ltp(lLa) : if(-2.6 - lLa, -.72, if(lLa - 1.9, lLa - 1.255,
    58 		(.249*lLa + .65)^2.7 - .72));
    59 				{ log10 of rod threshold luminance }
    60 lts(lLa) : if(-3.94 - lLa, -2.86, if(lLa - -1.44, lLa - .395,
    61 		(.405*lLa + 1.6)^2.18 - 2.86));
    62 				{ threshold is minimum of rods and cones }
    63 ldL2(lLa) : min(ltp(lLa),lts(lLa));
    64 dL(La) : 10^ldL2(log10(La));
    65 			{ derivative clamping function }
    66 clamp2(L, bLw) : dL(Lb(bLw))/dL(L)/Lb1(bLw)/(Bl(Ldmax)-Bl(Ldmin))/Bl1(L);
    67 clamp(L) : clamp2(L, BLw(L));
    68 			{ shift direction for histogram }
    69 shiftdir(B) : if(cf(B) - (B - Bl(Lmin))/(Bl(Ldmax) - Bl(Lmin)), 1, -1);
    70 			{ Scotopic/Photopic color adjustment }
    71 sL(r,g,b) : .062*r + .608*g + .330*b;	{ approx. scotopic brightness }
    72 BotMesopic : 10^-2.25;		{ top of scotopic range }
    73 TopMesopic : 10^0.75;		{ bottom of photopic range }
    74 incolor = if(Lw-TopMesopic, 1, if(BotMesopic-Lw, 0,
    75 		(Lw-BotMesopic)/(TopMesopic-BotMesopic)));
    76 slf = (1 - incolor)*sL(ri(1),gi(1),bi(1));
    77 ro = mult*(incolor*ri(1) + slf);
    78 go = mult*(incolor*gi(1) + slf);
    79 bo = mult*(incolor*bi(1) + slf);
    80 _EOF_
    81 else
    82 	cat >> $tf3 << _EOF_
    83 			{ derivative clamping function }
    84 clamp2(L, bLw) : Lb(bLw)/L/Lb1(bLw)/(Bl(Ldmax)-Bl(Ldmin))/Bl1(L);
    85 clamp(L) : clamp2(L, BLw(L));
    86 			{ shift direction for histogram }
    87 shiftdir(B) : -1;
    88 ro = mult*ri(1);
    89 go = mult*gi(1);
    90 bo = mult*bi(1);
    91 _EOF_
    92 endif
    93 # Compute brightness histogram
    94 pfilt -1 -p 1 -x 128 -y 128 $ifile | pvalue -o -b -d -h -H \
    95 	| rcalc -f $tf3 -e 'Lw=WE*$1;$1=if(Lw-Lmin,Bl(Lw),-1)' \
    96 	| histo `ev "log($Lmin)" "log($Lmax)"` $nsteps > $tf1
    97 # Clamp frequency distribution
    98 set totcount=`sed 's/^.*[ 	]//' $tf1 | total`
    99 set margin=`ev "floor($totcount*$cvratio+.5)"`
   100 while ( 1 )
   101 	# Compute mapping function
   102 	sed 's/^.*[ 	]//' $tf1 | total -1 -r \
   103 		| rcalc -e '$1=$1/'$totcount | rlam $tf1 - \
   104 		| tabfunc -i 0 cf > $tf4
   105 	# Compute difference with visible envelope
   106 	rcalc -f $tf4 -f $tf3 -e "T:$totcount*Stepsiz" \
   107 			-e 'clfq=floor(T*clamp(Lb($1))+.5)' \
   108 			-e '$1=$2-clfq;$2=shiftdir($1)' $tf1 > $tf1b
   109 	if (`sed 's/[	 ].*$//' $tf1b | total` >= 0) then
   110 		# Nothing visible? -- just normalize
   111 		pfilt -1 -e `pextrem $ifile | rcalc -e 'cond=recno-1.5;$1=1/(.265*$3+.67*$4+.065*$5)'` $ifile
   112 		goto quit
   113 	endif
   114 	# Check to see if we're close enough
   115 	if (`rcalc -e '$1=if($1,$1,0)' $tf1b | total` <= $margin) break
   116 	# Squash frequency distribution
   117 	set diffs=(`sed 's/[ 	].*$//' $tf1b`)
   118 	set shftd=(`sed 's/^.*[ 	]//' $tf1b`)
   119 	while ( 1 )
   120 		set maxi=0
   121 		set maxd=0
   122 		set i=$nsteps
   123 		while ( $i > 0 )
   124 			if ( $diffs[$i] > $maxd ) then
   125 				set maxd=$diffs[$i]
   126 				set maxi=$i
   127 			endif
   128 			@ i--
   129 		end
   130 		if ( $maxd == 0 ) break
   131 		set i=0
   132 	tryagain:
   133 		set r=$shftd[$maxi]
   134 		while ( $i == 0 )
   135 			@ t= $maxi + $r
   136 			if ( $t < 1 || $t > $nsteps ) then
   137 				@ shftd[$maxi]= -($shftd[$maxi])
   138 				goto tryagain
   139 			endif
   140 			if ( $diffs[$t] < 0 ) set i=$t
   141 			@ r+= $shftd[$maxi]
   142 		end
   143 		if ( $diffs[$i] <= -$diffs[$maxi] ) then
   144 			@ diffs[$i]+= $diffs[$maxi]
   145 			set diffs[$maxi]=0
   146 		else
   147 			@ diffs[$maxi]+= $diffs[$i]
   148 			set diffs[$i]=0
   149 		endif
   150 	end
   151 	# Mung histogram
   152 	echo $diffs | tr ' ' '\012' | rlam $tf1 - \
   153 		| rcalc -f $tf4 -f $tf3 -e "T:$totcount*Stepsiz" \
   154 			-e 'clfq=floor(T*clamp(Lb($1))+.5)' \
   155 			-e '$1=$1;$2=$3+clfq' > $tf1b
   156 	mv -f $tf1b $tf1
   157 end
   158 # Plot the mapping function if we are in debug mode
   159 if ( $?DEBUG ) then
   160 	cat > ${ibase}_histo.plt << _EOF_
   161 include=curve.plt
   162 title="Brightness Frequency Distribution"
   163 subtitle= $ibase
   164 ymin=0
   165 xlabel="Perceptual Brightness B(Lw)"
   166 ylabel="Frequency Count"
   167 Alabel="Histogram"
   168 Alintype=0
   169 Blabel="Envelope"
   170 Bsymsize=0
   171 Adata=
   172 _EOF_
   173 	(cat $tf1; echo \;; echo Bdata=) >> ${ibase}_histo.plt
   174 	rcalc -f $tf4 -f $tf3 -e "T:$totcount*Stepsiz" \
   175 		-e '$1=$1;$2=T*clamp(Lb($1))' $tf1 \
   176 		>> ${ibase}_histo.plt
   177 	cat > ${ibase}_brmap.plt << _EOF_
   178 include=line.plt
   179 title="Brightness Mapping Function"
   180 subtitle= $ibase
   181 xlabel="World Luminance (log cd/m^2)"
   182 ylabel="Display Luminance (cd/m^2)"
   183 ymax= $Ldmax
   184 Adata=
   185 _EOF_
   186 	cnt 100 | rcalc -f $tf4 -f $tf3 -e '$1=lx;$2=Lb(BLw(10^lx))' \
   187 		-e 'lx=$1/99*(log10(Lmax)-log10(Lmin))+log10(Lmin)' \
   188 		>> ${ibase}_brmap.plt
   189 	if ( $?DISPLAY ) then
   190 		bgraph ${ibase}_histo.plt ${ibase}_brmap.plt | x11meta &
   191 	endif
   192 endif
   193 # Map our picture
   194 getinfo < $ifile | egrep '^((VIEW|PIXASPECT|PRIMARIES)=|[^ ]*(rpict|rview|pinterp) )'
   195 pcomb -f $tf4 -f $tf3 $ifile
   196 quit:
   197 rm -f $tf