ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genklemsamp.pl
Revision: 2.4
Committed: Fri Jun 19 06:49:42 2009 UTC (14 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +14 -10 lines
Log Message:
Improvements and bug fixes

File Contents

# User Rev Content
1 greg 2.1 #!/usr/bin/perl -w
2 greg 2.4 # RCSid $Id: genklemsamp.pl,v 2.3 2009/06/16 17:18:38 greg Exp $
3 greg 2.1 #
4     # Sample Klems (full) directions impinging on surface(s)
5     #
6 greg 2.2 # G. Ward
7     #
8 greg 2.3 use strict;
9 greg 2.1 if ($#ARGV < 0) {
10 greg 2.4 print STDERR "Usage: genklemsamp [-c N ][-f{a|f|d}] [view opts] [geom.rad ..]\n";
11 greg 2.1 exit 1;
12     }
13     my $td = `mktemp -d /tmp/genklemsamp.XXXXXX`;
14     chomp $td;
15     my $nsamp = 1000;
16     my $fmt = "a";
17     my @vopts="-vo 1e-5";
18     while ($#ARGV >= 0) {
19     if ("$ARGV[0]" eq "-c") {
20     shift @ARGV;
21     $nsamp = $ARGV[0];
22     } elsif ("$ARGV[0]" =~ /^-f[afd]$/) {
23     $fmt = substr($ARGV[0], 2, 1);
24     } elsif ("$ARGV[0]" =~ /^-v[pdu]$/) {
25     push @vopts, "@ARGV[0..3]";
26     shift @ARGV; shift @ARGV; shift @ARGV;
27     } elsif ("$ARGV[0]" =~ /^-v[fhvo]$/) {
28     push @vopts, "@ARGV[0..1]";
29     shift @ARGV;
30     } elsif ("$ARGV[0]" =~ /^-v./) {
31 greg 2.4 print STDERR "Unsupported view option: $ARGV[0]\n";
32 greg 2.1 exit 1;
33     } elsif ("$ARGV[0]" =~ /^-./) {
34 greg 2.4 print STDERR "Unknown option: $ARGV[0]\n";
35 greg 2.1 exit 1;
36     } else {
37     last;
38     }
39     shift @ARGV;
40     }
41     my $ofmt = "";
42     $ofmt = "-o$fmt" if ("$fmt" ne "a");
43     # Require parallel view and compile settings
44     push @vopts, "-vtl";
45     my $vwset = `vwright @vopts V`;
46     $_ = $vwset; s/^.*Vhn://; s/;.*$//;
47     my $width = $_;
48     $_ = $vwset; s/^.*Vvn://; s/;.*$//;
49     my $height = $_;
50     my $sca = sqrt($nsamp/($width*$height));
51     # Kbin is input to get direction in full Klems basis with (x1,x2) randoms
52     my $tcal = '
53     DEGREE : PI/180;
54     Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90);
55     Knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
56     Kaccum(r) : if(r-.5, Knaz(r) + Kaccum(r-1), 0);
57     Kmax : Kaccum(Knaz(0));
58     Kfindrow(r, rem) : if(rem-Knaz(r)-.5, Kfindrow(r+1, rem-Knaz(r)), r);
59     Krow = if(Kbin-(Kmax+.5), 0, Kfindrow(1, Kbin));
60     Kcol = Kbin - Kaccum(Krow-1) - 1;
61     Kazi = 360*DEGREE * (Kcol + .5 - x2) / Knaz(Krow);
62     Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1));
63     sin_kpol = sin(Kpol);
64     K0 = cos(Kazi)*sin_kpol;
65     K1 = sin(Kazi)*sin_kpol;
66     K2 = cos(Kpol);
67     ';
68     my $ndiv = 145;
69     # Do we have any Radiance input files?
70     if ($#ARGV >= 0) {
71     system "oconv -f @ARGV > $td/surf.oct";
72     # Set our own view center and size based on bounding cube
73     $_ = $vwset; s/^.*Vdx://; s/;.*$//;
74     my @vd = $_;
75     $_ = $vwset; s/^.*Vdy://; s/;.*$//;
76     push @vd, $_;
77     $_ = $vwset; s/^.*Vdz://; s/;.*$//;
78     push @vd, $_;
79     my @bcube = split /\s+/, `getinfo -d < $td/surf.oct`;
80     $width = $bcube[3]*sqrt(3);
81     $height = $width;
82     push @vopts, ("-vp", $bcube[0]+$bcube[3]/2-$width/2*$vd[0],
83     $bcube[1]+$bcube[3]/2-$width/2*$vd[1],
84     $bcube[2]+$bcube[3]/2-$width/2*$vd[2]);
85     push @vopts, ("-vh", $width, "-vh", $height);
86     $vwset = `vwright @vopts V`;
87     my $xres;
88     my $yres;
89 greg 2.4 my $ntot = 0;
90 greg 2.1 # This generally passes through the loop twice to get density right
91 greg 2.4 while ($ntot < $nsamp) {
92 greg 2.1 $xres = int($width*$sca) + 1;
93     $yres = int($height*$sca) + 1;
94     system "vwrays -ff -x $xres -y $yres -pa 0 -pj .7 @vopts " .
95     "| rtrace -w -h -ff -opLN $td/surf.oct " .
96     q{| rcalc -e 'cond=and(5e9-$4,nOK);$1=$1;$2=$2;$3=$3' } .
97     "-e 'and(a,b):if(a,b,a);sq(x):x*x' -e '$vwset' " .
98     q{-e 'nOK=sq(Vdx*$5+Vdy*$6+Vdz*$7)-.999' } .
99     "-if7 -of > $td/origins.flt";
100     $ntot = -s "$td/origins.flt";
101     $ntot /= 3*4;
102     if ($ntot == 0) {
103 greg 2.4 if ($sca < sqrt(199/($width*$height))) {
104     $sca = sqrt(200/($width*$height));
105     redo;
106     }
107     print STDERR "View direction does not correspond to any surfaces\n";
108 greg 2.1 exit 1;
109     }
110     $sca *= 1.05 * sqrt($nsamp/$ntot);
111 greg 2.4 }
112 greg 2.1 # All set to produce our samples
113 greg 2.4 for (my $k = 1; $k <= $ndiv; $k++) {
114 greg 2.1 my $rn = rand(10);
115     my $r1 = rand; my $r2 = rand;
116     # Chance of using = (number_still_needed)/(number_left_avail)
117     system "rcalc $ofmt -e '$tcal' -e '$vwset' " .
118     "-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
119     "-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
120     q{-e '$1=$1+Vo*Vdx; $2=$2+Vo*Vdy; $3=$3+Vo*Vdz' } .
121     q{-e '$4=K0*Vhx+K1*Vvx+K2*Vdx' } .
122     q{-e '$5=K0*Vhy+K1*Vvy+K2*Vdy' } .
123     q{-e '$6=K0*Vhz+K1*Vvz+K2*Vdz' } .
124     "-if3 $td/origins.flt";
125     }
126     } else {
127     # No Radiance input files, so just sample over parallel view
128     my $xres = int($width*$sca) + 1;
129     my $yres = int($height*$sca) + 1;
130     my $ntot = $xres * $yres;
131 greg 2.4 for (my $k = 1; $k <= $ndiv; $k++) {
132 greg 2.1 my $rn = rand(10);
133     my $r1 = rand; my $r2 = rand;
134     my $r3 = rand; my $r4 = rand;
135     # Chance of using = (number_still_needed)/(number_left_avail)
136     system "cnt $yres $xres | rcalc $ofmt -e '$tcal' -e '$vwset' " .
137     q{-e 'xc=$2;yc=$1' } .
138     "-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
139     "-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
140     "-e 'hpos=Vhn*((xc+rand($r3+recno))/$xres - .5)' " .
141     "-e 'vpos=Vvn*((yc+rand($r4+recno))/$yres - .5)' " .
142     q{-e '$1=Vpx+Vo*Vdx+hpos*Vhx+vpos*Vvx' } .
143     q{-e '$2=Vpy+Vo*Vdy+hpos*Vhy+vpos*Vvy' } .
144     q{-e '$3=Vpz+Vo*Vdz+hpos*Vhz+vpos*Vvz' } .
145     q{-e '$4=K0*Vhx+K1*Vvx+K2*Vdx' } .
146     q{-e '$5=K0*Vhy+K1*Vvy+K2*Vdy' } .
147     q{-e '$6=K0*Vhz+K1*Vvz+K2*Vdz' } ;
148     }
149     }
150     system "rm -rf $td";