ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genklemsamp.pl
Revision: 2.10
Committed: Sat Apr 16 00:39:07 2011 UTC (11 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +3 -2 lines
Log Message:
Bug fix for 32-bit OS

File Contents

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