ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genBSDF.pl
Revision: 2.1
Committed: Thu Sep 2 02:29:24 2010 UTC (13 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial version of genBSDF, probably buggy

File Contents

# Content
1 #!/usr/bin/perl -w
2 # RCSid $Id$
3 #
4 # Compute BSDF based on geometry and material description
5 #
6 # G. Ward
7 #
8 use strict;
9 sub userror {
10 print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-dim xmin xmax ymin ymax zmin zmax][{+|-}mgf][{+|-}geom] [input ..]";
11 exit 1;
12 }
13 my $td = `mktemp -d /tmp/genBSDF.XXXXXX`;
14 chomp $td;
15 my $nsamp = 1000;
16 my $mgfin = 0;
17 my $geout = 1;
18 my $nproc = 1;
19 my @dim;
20 # Get options
21 while ($#ARGV >= 0) {
22 if ("$ARGV[0]" =~ /^[-+]m/) {
23 $mgfin = ("$ARGV[0]" =~ /^\+/);
24 } elsif ("$ARGV[0]" =~ /^[-+]g/) {
25 $geout = ("$ARGV[0]" =~ /^\+/);
26 } elsif ("$ARGV[0]" eq "-c") {
27 $nsamp = $ARGV[1];
28 shift @ARGV;
29 } elsif ("$ARGV[0]" eq "-n") {
30 $nproc = $ARGV[1];
31 shift @ARGV;
32 } elsif ("$ARGV[0]" =~ /^-d/) {
33 userror() if ($#ARGV < 6);
34 @dim = "@ARGV[1..6]";
35 shift @ARGV for (1..6);
36 } elsif ("$ARGV[0]" =~ /^[-+]./) {
37 userror();
38 } else {
39 last;
40 }
41 shift @ARGV;
42 }
43 # Get scene description and dimensions
44 my $radscn = "$td/device.rad";
45 my $mgfscn = "$td/device.mgf";
46 my $octree = "$td/device.oct";
47 if ( $mgfin ) {
48 system "mgfilt '#,o,xf,c,cxy,cspec,cmix,m,sides,rd,td,rs,ts,ir,v,p,n,f,fh,sph,cyl,cone,prism,ring,torus' @ARGV > $mgfscn";
49 die "Could not load MGF input\n" if ( $? );
50 system "mgf2rad $mgfscn > $radscn";
51 } else {
52 system "cat @ARGV | xform -e > $radscn";
53 die "Could not load Radiance input\n" if ( $? );
54 system "rad2mgf $radscn > $mgfscn" if ( $geout );
55 }
56 if ($#dim != 5) {
57 @dim = split /\s+/, `getbbox -h $radscn`;
58 shift @dim;
59 }
60 print STDERR "Warning: Device extends into room\n" if ($dim[5] > 1e-5);
61 # Add receiver surface (rectangle)
62 my $modnm="_receiver_black_";
63 open(RADSCN, ">> $radscn");
64 print RADSCN "void glow $modnm\n0\n0\n4 0 0 0 0\n\n";
65 print RADSCN "$modnm polygon _receiver_\n0\n0\n12\n";
66 print RADSCN "\t",$dim[0],"\t",$dim[2],"\t",$dim[5]+1e-5,"\n";
67 print RADSCN "\t",$dim[0],"\t",$dim[3],"\t",$dim[5]+1e-5,"\n";
68 print RADSCN "\t",$dim[1],"\t",$dim[3],"\t",$dim[5]+1e-5,"\n";
69 print RADSCN "\t",$dim[1],"\t",$dim[2],"\t",$dim[5]+1e-5,"\n";
70 close RADSCN;
71 # Generate octree
72 system "oconv -w $radscn > $octree";
73 die "Could not compile scene\n" if ( $? );
74 # Set up sampling
75 # Kbin to produce incident direction in full Klems basis with (x1,x2) randoms
76 my $tcal = '
77 DEGREE : PI/180;
78 Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90);
79 Knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
80 Kaccum(r) : if(r-.5, Knaz(r) + Kaccum(r-1), 0);
81 Kmax : Kaccum(Knaz(0));
82 Kfindrow(r, rem) : if(rem-Knaz(r)+.5, Kfindrow(r+1, rem-Knaz(r)), r);
83 Krow = if(Kbin-(Kmax-.5), 0, Kfindrow(1, Kbin));
84 Kcol = Kbin - Kaccum(Krow-1);
85 Kazi = 360*DEGREE * (Kcol + .5 - x2) / Knaz(Krow);
86 Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1));
87 sin_kpol = sin(Kpol);
88 Dx = -cos(Kazi)*sin_kpol;
89 Dy = sin(Kazi)*sin_kpol;
90 Dz = sqrt(1 - sin_kpol*sin_kpol);
91 ';
92 # Compute Klems bin from exiting ray direction
93 my $kcal = '
94 DEGREE : PI/180;
95 Acos(x) : 1/DEGREE * if(x-1, 0, if(-1-x, 0, acos(x)));
96 posangle(a) : if(-a, a + 2*PI, a);
97 Atan2(y,x) : 1/DEGREE * posangle(atan2(y,x));
98 kpola(r) : select(r, 5, 15, 25, 35, 45, 55, 65, 75, 90);
99 knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
100 kaccum(r) : if(r-.5, knaz(r) + kaccum(r-1), 0);
101 kfindrow(r, pol) : if(r-kpola(0)+.5, r,
102 if(pol-kpola(r), kfindrow(r+1, pol), r) );
103 kazn(azi,inc) : if((360-.5*inc)-azi, floor((azi+.5*inc)/inc), 0);
104 kbin2(pol,azi) = select(kfindrow(1, pol),
105 kazn(azi,360/knaz(1)),
106 kaccum(1) + kazn(azi,360/knaz(2)),
107 kaccum(2) + kazn(azi,360/knaz(3)),
108 kaccum(3) + kazn(azi,360/knaz(4)),
109 kaccum(4) + kazn(azi,360/knaz(5)),
110 kaccum(5) + kazn(azi,360/knaz(6)),
111 kaccum(6) + kazn(azi,360/knaz(7)),
112 kaccum(7) + kazn(azi,360/knaz(8)),
113 kaccum(8) + kazn(azi,360/knaz(9))
114 );
115 kbin = kbin2(Acos(Dz), Atan2(Dy, -Dx));
116 ';
117 my $ndiv = 145;
118 my $nx = int(sqrt($nsamp*($dim[1]-$dim[0])/($dim[3]-$dim[2])) + .5);
119 my $ny = int($nsamp/$nx + .5);
120 $nsamp = $nx * $ny;
121 # Output XML prologue
122 print
123 '<?xml version="1.0" encoding="UTF-8"?>
124 <WindowElement xmlns="http://windows.lbl.gov" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://windows.lbl.gov/BSDF-v1.4.xsd">
125 <WindowElementType>System</WindowElementType>
126 <Optical>
127 <Layer>
128 <Material>
129 <Name>Name</Name>
130 <Manufacturer>Manufacturer</Manufacturer>
131 ';
132 printf "\t\t\t<Thickness unit=\"Meter\">%.3f</Thickness>\n", $dim[5] - $dim[4];
133 printf "\t\t\t<Width unit=\"Meter\">%.3f</Width>\n", $dim[1] - $dim[0];
134 printf "\t\t\t<Height unit=\"Meter\">%.3f</Height>\n", $dim[3] - $dim[2];
135 print "\t\t\t<DeviceType>Integral</DeviceType>\n";
136 # Output MGF description if requested
137 if ( $geout ) {
138 print "\t\t\t<Geometry format=\"MGF\" unit=\"Meter\">\n";
139 printf "xf -t %.6f %.6f 0\n", -($dim[0]+$dim[1])/2, -($dim[2]+$dim[3])/2;
140 system "cat $mgfscn";
141 print "xf\n";
142 print "\t\t\t</Geometry>\n";
143 }
144 print ' </Material>
145 <DataDefinition>
146 <IncidentDataStructure>Columns</IncidentDataStructure>
147 <AngleBasis>
148 <AngleBasisName>LBNL/Klems Full</AngleBasisName>
149 <AngleBasisBlock>
150 <Theta>0</Theta>
151 <nPhis>1</nPhis>
152 <ThetaBounds>
153 <LowerTheta>0</LowerTheta>
154 <UpperTheta>5</UpperTheta>
155 </ThetaBounds>
156 </AngleBasisBlock>
157 <AngleBasisBlock>
158 <Theta>10</Theta>
159 <nPhis>8</nPhis>
160 <ThetaBounds>
161 <LowerTheta>5</LowerTheta>
162 <UpperTheta>15</UpperTheta>
163 </ThetaBounds>
164 </AngleBasisBlock>
165 <AngleBasisBlock>
166 <Theta>20</Theta>
167 <nPhis>16</nPhis>
168 <ThetaBounds>
169 <LowerTheta>15</LowerTheta>
170 <UpperTheta>25</UpperTheta>
171 </ThetaBounds>
172 </AngleBasisBlock>
173 <AngleBasisBlock>
174 <Theta>30</Theta>
175 <nPhis>20</nPhis>
176 <ThetaBounds>
177 <LowerTheta>25</LowerTheta>
178 <UpperTheta>35</UpperTheta>
179 </ThetaBounds>
180 </AngleBasisBlock>
181 <AngleBasisBlock>
182 <Theta>40</Theta>
183 <nPhis>24</nPhis>
184 <ThetaBounds>
185 <LowerTheta>35</LowerTheta>
186 <UpperTheta>45</UpperTheta>
187 </ThetaBounds>
188 </AngleBasisBlock>
189 <AngleBasisBlock>
190 <Theta>50</Theta>
191 <nPhis>24</nPhis>
192 <ThetaBounds>
193 <LowerTheta>45</LowerTheta>
194 <UpperTheta>55</UpperTheta>
195 </ThetaBounds>
196 </AngleBasisBlock>
197 <AngleBasisBlock>
198 <Theta>60</Theta>
199 <nPhis>24</nPhis>
200 <ThetaBounds>
201 <LowerTheta>55</LowerTheta>
202 <UpperTheta>65</UpperTheta>
203 </ThetaBounds>
204 </AngleBasisBlock>
205 <AngleBasisBlock>
206 <Theta>70</Theta>
207 <nPhis>16</nPhis>
208 <ThetaBounds>
209 <LowerTheta>65</LowerTheta>
210 <UpperTheta>75</UpperTheta>
211 </ThetaBounds>
212 </AngleBasisBlock>
213 <AngleBasisBlock>
214 <Theta>82.5</Theta>
215 <nPhis>12</nPhis>
216 <ThetaBounds>
217 <LowerTheta>75</LowerTheta>
218 <UpperTheta>90</UpperTheta>
219 </ThetaBounds>
220 </AngleBasisBlock>
221 </AngleBasis>
222 </DataDefinition>
223 <WavelengthData>
224 <LayerNumber>System</LayerNumber>
225 <Wavelength unit="Integral">Visible</Wavelength>
226 <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
227 <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
228 <WavelengthDataBlock>
229 <WavelengthDataDirection>Transmission Front</WavelengthDataDirection>
230 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
231 <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
232 <ScatteringDataType>BTDF</ScatteringDataType>
233 <ScatteringData>
234 ';
235 # Compute actual scattering data using rtcontrib
236 system "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " .
237 "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " .
238 "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " .
239 "-e 'zp:$dim[4]-1e-5' " .
240 q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } .
241 q{-e '$1=xp;$2=yp;$3=zp;$4=Dx;$5=Dy;$6=Dz' } .
242 "| rtcontrib -h -ff -n $nproc -c $nsamp -e '$kcal' -b kbin -bn $ndiv " .
243 "-m $modnm -w -ab 4 -lw 1e-5 $octree " .
244 q{| rcalc -if3 -e '$1=0.265*$1+0.670*$2+0.065*$3'};
245 # Output XML epilogue
246 print
247 ' </ScatteringData>
248 </WavelengthDataBlock>
249 </WavelengthData>
250 </Layer>
251 </Optical>
252 </WindowElement>
253 ';
254 # Clean up temporary files
255 system "rm -rf $td";