ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/psquish.csh
Revision: 3.1
Committed: Fri Sep 20 16:44:00 1996 UTC (27 years, 7 months ago) by greg
Content type: application/x-csh
Branch: MAIN
Log Message:
Initial revision

File Contents

# User Rev Content
1 greg 3.1 #!/bin/csh -f
2     # SCCSid "$SunId$ LBL"
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=/usr/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.pic > output.pic"
22     exit 1
23     endif
24     set ifile=$1
25     set ibase=$ifile:t
26     if ( "$ibase" =~ *.pic ) 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 | lam $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' | lam $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