|
1 #!/bin/csh -f |
|
2 # RCSid: $Id: vlpic.csh,v 3.4 2008/08/25 04:50:32 greg Exp $ |
|
3 # |
|
4 # Compute falsecolor image of visibility level |
|
5 # using the wacky formulas of Werner Adrian. |
|
6 # |
|
7 set age=40 # default age (years) |
|
8 set tim=.2 # default time (seconds) |
|
9 while ($#argv > 1) |
|
10 switch ($argv[1]) |
|
11 case -a: |
|
12 shift argv |
|
13 set age="$argv[1]" |
|
14 breaksw |
|
15 case -t: |
|
16 shift argv |
|
17 set tim="$argv[1]" |
|
18 breaksw |
|
19 default: |
|
20 echo bad option "'$argv[1]'" |
|
21 exit 1 |
|
22 endsw |
|
23 shift argv |
|
24 end |
|
25 set tc=`mktemp /tmp/vlcal.XXXXXXX` |
|
26 set tp1=`mktemp /tmp/vlpic.XXXXXX` |
|
27 set tp2=`mktemp /tmp/vlr2pic.XXXXXX` |
|
28 set tp4=`mktemp /tmp/vlr4pic.XXXXXX` |
|
29 set tp8=`mktemp /tmp/vlr8pic.XXXXXX` |
|
30 set tf=($tc $tp1 $tp2 $tp4 $tp8) |
|
31 set inpic=$argv[1] |
|
32 onintr quit |
|
33 set pr=(`getinfo -d < $inpic | sed 's/^-Y \([1-9][0-9]*\) +X \([1-9][0-9]*\)$/\2 \1/'`) |
|
34 # ( vwright V < $inpic ; cat ) > $tc << _EOF_ |
|
35 cat > $tc << _EOF_ |
|
36 { A : 3438 * sqrt(Vhn/xmax*Vvn/ymax); { pixel size (in minutes) } } |
|
37 PI : 3.14159265358979323846; |
|
38 A = 180*60/PI * sqrt(S(1)/PI); |
|
39 sq(x) : x*x; |
|
40 max(a,b) : if(a-b, a, b); |
|
41 Lv : 0; { veiling luminance {temporarily zero} } |
|
42 { find direction of maximum contrast } |
|
43 xoff(d) = select(d, 1, 1, 0, -1, -1, -1, 0, 1); |
|
44 yoff(d) = select(d, 0, 1, 1, 1, 0, -1, -1, -1); |
|
45 cf(lf, lb) = (lf - lb)/(lb + Lv); |
|
46 contrast(d) = cf(li(1,xoff(d),yoff(d)), li(1,-xoff(d),-yoff(d))); |
|
47 cwin(d1, d2) = if(contrast(d1) - contrast(d2), d1, d2); |
|
48 bestdir(d) = if(d-1.5, cwin(d, bestdir(d-1)), d); |
|
49 maxdir = bestdir(8); |
|
50 Lt = WE*li(1,xoff(maxdir),yoff(maxdir)) + Lv; |
|
51 La = WE*li(1,-xoff(maxdir),-yoff(maxdir)) + Lv; |
|
52 LLa = log10(La); |
|
53 F = if(La-.6, log10(4.1925) + LLa*.1556 + .1684*La^.5867, |
|
54 if(La-.00418, 10^(-.072 + LLa*(.3372 + LLa*.0866)), |
|
55 10^(.028 + .173*LLa))); |
|
56 L = if(La-.6, .05946*La^.466, |
|
57 if(La-.00418, 10^(-1.256 + .319*LLa), |
|
58 10^(-.891 + LLa*(.5275 + LLa*.0227)))); |
|
59 LAt = log10(A) + .523; |
|
60 AA = .36 - .0972*LAt*LAt/(LAt*(LAt - 2.513) + 2.7895); |
|
61 LLat = LLa + 6; |
|
62 AL = .355 - .1217*LLat*LLat/(LLat*(LLat - 10.4) + 52.28); |
|
63 AZ = sqrt(AA*AA + AL*AL)/2.1; |
|
64 DL1 = 2.6*sq(F/A + L); |
|
65 M = if(La-.004, 10^-(10^-(if(La-.1,.125,.075)*sq(LLa+1) + .0245)), 0); |
|
66 TGB = .6*La^-.1488; |
|
67 FCP = 1 - M*A^-TGB/(2.4*(DL1*(AZ+2)/2)); |
|
68 DL2 = DL1*(AZ + T)/T; |
|
69 FA = if(Age-64, sq(Age-56.6)/116.3 + 1.43, sq(Age-19)/2160 + .99); |
|
70 DL3 = DL2 * FA; |
|
71 DL4 = if(La-Lt, DL3*FCP, DL3); |
|
72 lo = (Lt - La)/DL4; { Output VL } |
|
73 _EOF_ |
|
74 pcomb -w -e "Age:$age;T:$tim" -f $tc -o $inpic > $tp1 |
|
75 pfilt -1 -x /2 -y /2 $inpic \ |
|
76 | pcomb -w -e "Age:$age;T:$tim" -f $tc -o - \ |
|
77 | pfilt -1 -r 1 -x $pr[1] -y $pr[2] > $tp2 |
|
78 pfilt -1 -x /4 -y /4 $inpic \ |
|
79 | pcomb -w -e "Age:$age;T:$tim" -f $tc -o - \ |
|
80 | pfilt -1 -r 1 -x $pr[1] -y $pr[2] > $tp4 |
|
81 pfilt -1 -x /8 -y /8 $inpic \ |
|
82 | pcomb -w -e "Age:$age;T:$tim" -f $tc -o - \ |
|
83 | pfilt -1 -r 1 -x $pr[1] -y $pr[2] > $tp8 |
|
84 pcomb -e 'max(a,b):if(a-b,a,b)' -e 'lo=max(gi(1),max(gi(2),max(gi(3),gi(4))))' \ |
|
85 $tp1 $tp2 $tp4 $tp8 \ |
|
86 | falsecolor -l VL -s 56 -m 1 -r 'if(v-1/8,1.6*8/7*(v-1/8)-.6,0)' \ |
|
87 -g 'if(v-1/8,if(v-.453125,1.6-1.6*8/7*(v-1/8),8/3*8/7*(v-1/8)),0)' \ |
|
88 -b 'if(v-1/8,1-8/3*8/7*(v-1/8),0)' |
|
89 quit: |
|
90 rm -f $tf |