ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genBSDF.pl
Revision: 2.9
Committed: Mon Feb 21 22:48:51 2011 UTC (13 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.8: +145 -27 lines
Log Message:
Added BRDF computation to genBSDF

File Contents

# Content
1 #!/usr/bin/perl -w
2 # RCSid $Id: genBSDF.pl,v 2.8 2011/01/29 16:44:01 greg Exp $
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][{+|-}f][{+|-}b][{+|-}mgf][{+|-}geom] [input ..]\n";
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 $doforw = 0;
20 my $doback = 1;
21 my @dim;
22 # Get options
23 while ($#ARGV >= 0) {
24 if ("$ARGV[0]" =~ /^[-+]m/) {
25 $mgfin = ("$ARGV[0]" =~ /^\+/);
26 } elsif ("$ARGV[0]" =~ /^[-+]g/) {
27 $geout = ("$ARGV[0]" =~ /^\+/);
28 } elsif ("$ARGV[0]" =~ /^[-+]f/) {
29 $doforw = ("$ARGV[0]" =~ /^\+/);
30 } elsif ("$ARGV[0]" =~ /^[-+]b/) {
31 $doback = ("$ARGV[0]" =~ /^\+/);
32 } elsif ("$ARGV[0]" eq "-c") {
33 $nsamp = $ARGV[1];
34 shift @ARGV;
35 } elsif ("$ARGV[0]" eq "-n") {
36 $nproc = $ARGV[1];
37 shift @ARGV;
38 } elsif ("$ARGV[0]" =~ /^-d/) {
39 userror() if ($#ARGV < 6);
40 @dim = @ARGV[1..6];
41 shift @ARGV for (1..6);
42 } elsif ("$ARGV[0]" =~ /^[-+]./) {
43 userror();
44 } else {
45 last;
46 }
47 shift @ARGV;
48 }
49 # Check that we're actually being asked to do something
50 die "Must have at least one of +forward or +backward" if (!$doforw && !$doback);
51 # Get scene description and dimensions
52 my $radscn = "$td/device.rad";
53 my $mgfscn = "$td/device.mgf";
54 my $octree = "$td/device.oct";
55 if ( $mgfin ) {
56 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";
57 die "Could not load MGF input\n" if ( $? );
58 system "mgf2rad $mgfscn > $radscn";
59 } else {
60 system "cat @ARGV | xform -e > $radscn";
61 die "Could not load Radiance input\n" if ( $? );
62 system "rad2mgf $radscn > $mgfscn" if ( $geout );
63 }
64 if ($#dim != 5) {
65 @dim = split ' ', `getbbox -h $radscn`;
66 }
67 print STDERR "Warning: Device extends into room\n" if ($dim[5] > 1e-5);
68 # Add receiver surfaces (rectangular)
69 my $bmodnm="receiver_behind";
70 my $fmodnm="receiver_face";
71 open(RADSCN, ">> $radscn");
72 print RADSCN "void glow $fmodnm\n0\n0\n4 0 0 0 0\n\n";
73 print RADSCN "$fmodnm polygon f_receiver\n0\n0\n12\n";
74 print RADSCN "\t",$dim[0],"\t",$dim[2],"\t",$dim[5]+2e-5,"\n";
75 print RADSCN "\t",$dim[0],"\t",$dim[3],"\t",$dim[5]+2e-5,"\n";
76 print RADSCN "\t",$dim[1],"\t",$dim[3],"\t",$dim[5]+2e-5,"\n";
77 print RADSCN "\t",$dim[1],"\t",$dim[2],"\t",$dim[5]+2e-5,"\n";
78 print RADSCN "void glow $bmodnm\n0\n0\n4 0 0 0 0\n\n";
79 print RADSCN "$bmodnm polygon b_receiver\n0\n0\n12\n";
80 print RADSCN "\t",$dim[1],"\t",$dim[2],"\t",$dim[4]-2e-5,"\n";
81 print RADSCN "\t",$dim[1],"\t",$dim[3],"\t",$dim[4]-2e-5,"\n";
82 print RADSCN "\t",$dim[0],"\t",$dim[3],"\t",$dim[4]-2e-5,"\n";
83 print RADSCN "\t",$dim[0],"\t",$dim[2],"\t",$dim[4]-2e-5,"\n";
84 close RADSCN;
85 # Generate octree
86 system "oconv -w $radscn > $octree";
87 die "Could not compile scene\n" if ( $? );
88 # Set up sampling of interior portal
89 # Kbin to produce incident direction in full Klems basis with (x1,x2) randoms
90 my $tcal = '
91 DEGREE : PI/180;
92 sq(x) : x*x;
93 Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90);
94 Knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
95 Kaccum(r) : if(r-.5, Knaz(r) + Kaccum(r-1), 0);
96 Kmax : Kaccum(Knaz(0));
97 Kfindrow(r, rem) : if(rem-Knaz(r)+.5, Kfindrow(r+1, rem-Knaz(r)), r);
98 Krow = if(Kbin-(Kmax-.5), 0, Kfindrow(1, Kbin));
99 Kcol = Kbin - Kaccum(Krow-1);
100 Kazi = 360*DEGREE * (Kcol + (.5 - x2)) / Knaz(Krow);
101 Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1));
102 sin_kpol = sin(Kpol);
103 Dx = cos(Kazi)*sin_kpol;
104 Dy = sin(Kazi)*sin_kpol;
105 Dz = sqrt(1 - sin_kpol*sin_kpol);
106 KprojOmega = PI * if(Kbin-.5,
107 (sq(cos(Kpola(Krow-1)*DEGREE)) - sq(cos(Kpola(Krow)*DEGREE)))/Knaz(Krow),
108 1 - sq(cos(Kpola(1)*DEGREE)));
109 ';
110 # Compute Klems bin from exiting ray direction (forward or backward)
111 my $kcal = '
112 DEGREE : PI/180;
113 Acos(x) : 1/DEGREE * if(x-1, 0, if(-1-x, 0, acos(x)));
114 posangle(a) : if(-a, a + 2*PI, a);
115 Atan2(y,x) : 1/DEGREE * posangle(atan2(y,x));
116 kpola(r) : select(r, 5, 15, 25, 35, 45, 55, 65, 75, 90);
117 knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
118 kaccum(r) : if(r-.5, knaz(r) + kaccum(r-1), 0);
119 kfindrow(r, pol) : if(r-kpola(0)+.5, r,
120 if(pol-kpola(r), kfindrow(r+1, pol), r) );
121 kazn(azi,inc) : if((360-.5*inc)-azi, floor((azi+.5*inc)/inc), 0);
122 kbin2(pol,azi) = select(kfindrow(1, pol),
123 kazn(azi,360/knaz(1)),
124 kaccum(1) + kazn(azi,360/knaz(2)),
125 kaccum(2) + kazn(azi,360/knaz(3)),
126 kaccum(3) + kazn(azi,360/knaz(4)),
127 kaccum(4) + kazn(azi,360/knaz(5)),
128 kaccum(5) + kazn(azi,360/knaz(6)),
129 kaccum(6) + kazn(azi,360/knaz(7)),
130 kaccum(7) + kazn(azi,360/knaz(8)),
131 kaccum(8) + kazn(azi,360/knaz(9))
132 );
133 kbin = if(Dz, kbin2(Acos(Dz),Atan2(Dy,Dx)), kbin2(Acos(-Dz),Atan2(-Dy,-Dx)));
134 ';
135 my $ndiv = 145;
136 my $nx = int(sqrt($nsamp*($dim[1]-$dim[0])/($dim[3]-$dim[2])) + .5);
137 my $ny = int($nsamp/$nx + .5);
138 $nsamp = $nx * $ny;
139 # Compute scattering data using rtcontrib
140 my @tfarr;
141 my @rfarr;
142 my @tbarr;
143 my @rbarr;
144 my $cmd;
145 my $rtargs = "-w -ab 5 -ad 700 -lw 3e-6";
146 my $rtcmd = "rtcontrib -h -ff -fo -n $nproc -c $nsamp " .
147 "-e '$kcal' -b kbin -bn $ndiv " .
148 "-o '$td/%s.flt' -m $fmodnm -m $bmodnm $rtargs $octree";
149 my $rccmd = "rcalc -e '$tcal' " .
150 "-e 'mod(n,d):n-floor(n/d)*d' -e 'Kbin=mod(recno-.999,$ndiv)' " .
151 q{-if3 -e '$1=(0.265*$1+0.670*$2+0.065*$3)/KprojOmega'};
152 if ( $doforw ) {
153 $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " .
154 "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " .
155 "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " .
156 "-e 'zp:$dim[4]-1e-5' " .
157 q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } .
158 q{-e '$1=xp;$2=yp;$3=zp;$4=Dx;$5=Dy;$6=Dz' } .
159 "| $rtcmd";
160 system "$cmd" || die "Failure running: $cmd\n";
161 @tfarr = `$rccmd $td/$fmodnm.flt`;
162 die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? );
163 @rfarr = `$rccmd $td/$bmodnm.flt`;
164 die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? );
165 }
166 if ( $doback ) {
167 $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " .
168 "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " .
169 "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " .
170 "-e 'zp:$dim[5]+1e-5' " .
171 q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } .
172 q{-e '$1=xp;$2=yp;$3=zp;$4=-Dx;$5=-Dy;$6=-Dz' } .
173 "| $rtcmd";
174 system "$cmd" || die "Failure running: $cmd\n";
175 @tbarr = `$rccmd $td/$bmodnm.flt`;
176 die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? );
177 @rbarr = `$rccmd $td/$fmodnm.flt`;
178 die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? );
179 }
180 # Output XML prologue
181 print
182 '<?xml version="1.0" encoding="UTF-8"?>
183 <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">
184 <WindowElementType>System</WindowElementType>
185 <Optical>
186 <Layer>
187 <Material>
188 <Name>Name</Name>
189 <Manufacturer>Manufacturer</Manufacturer>
190 ';
191 printf "\t\t\t<Thickness unit=\"Meter\">%.3f</Thickness>\n", $dim[5] - $dim[4];
192 printf "\t\t\t<Width unit=\"Meter\">%.3f</Width>\n", $dim[1] - $dim[0];
193 printf "\t\t\t<Height unit=\"Meter\">%.3f</Height>\n", $dim[3] - $dim[2];
194 print "\t\t\t<DeviceType>Integral</DeviceType>\n";
195 # Output MGF description if requested
196 if ( $geout ) {
197 print "\t\t\t<Geometry format=\"MGF\" unit=\"Meter\">\n";
198 printf "xf -t %.6f %.6f 0\n", -($dim[0]+$dim[1])/2, -($dim[2]+$dim[3])/2;
199 system "cat $mgfscn";
200 print "xf\n";
201 print "\t\t\t</Geometry>\n";
202 }
203 print ' </Material>
204 <DataDefinition>
205 <IncidentDataStructure>Columns</IncidentDataStructure>
206 <AngleBasis>
207 <AngleBasisName>LBNL/Klems Full</AngleBasisName>
208 <AngleBasisBlock>
209 <Theta>0</Theta>
210 <nPhis>1</nPhis>
211 <ThetaBounds>
212 <LowerTheta>0</LowerTheta>
213 <UpperTheta>5</UpperTheta>
214 </ThetaBounds>
215 </AngleBasisBlock>
216 <AngleBasisBlock>
217 <Theta>10</Theta>
218 <nPhis>8</nPhis>
219 <ThetaBounds>
220 <LowerTheta>5</LowerTheta>
221 <UpperTheta>15</UpperTheta>
222 </ThetaBounds>
223 </AngleBasisBlock>
224 <AngleBasisBlock>
225 <Theta>20</Theta>
226 <nPhis>16</nPhis>
227 <ThetaBounds>
228 <LowerTheta>15</LowerTheta>
229 <UpperTheta>25</UpperTheta>
230 </ThetaBounds>
231 </AngleBasisBlock>
232 <AngleBasisBlock>
233 <Theta>30</Theta>
234 <nPhis>20</nPhis>
235 <ThetaBounds>
236 <LowerTheta>25</LowerTheta>
237 <UpperTheta>35</UpperTheta>
238 </ThetaBounds>
239 </AngleBasisBlock>
240 <AngleBasisBlock>
241 <Theta>40</Theta>
242 <nPhis>24</nPhis>
243 <ThetaBounds>
244 <LowerTheta>35</LowerTheta>
245 <UpperTheta>45</UpperTheta>
246 </ThetaBounds>
247 </AngleBasisBlock>
248 <AngleBasisBlock>
249 <Theta>50</Theta>
250 <nPhis>24</nPhis>
251 <ThetaBounds>
252 <LowerTheta>45</LowerTheta>
253 <UpperTheta>55</UpperTheta>
254 </ThetaBounds>
255 </AngleBasisBlock>
256 <AngleBasisBlock>
257 <Theta>60</Theta>
258 <nPhis>24</nPhis>
259 <ThetaBounds>
260 <LowerTheta>55</LowerTheta>
261 <UpperTheta>65</UpperTheta>
262 </ThetaBounds>
263 </AngleBasisBlock>
264 <AngleBasisBlock>
265 <Theta>70</Theta>
266 <nPhis>16</nPhis>
267 <ThetaBounds>
268 <LowerTheta>65</LowerTheta>
269 <UpperTheta>75</UpperTheta>
270 </ThetaBounds>
271 </AngleBasisBlock>
272 <AngleBasisBlock>
273 <Theta>82.5</Theta>
274 <nPhis>12</nPhis>
275 <ThetaBounds>
276 <LowerTheta>75</LowerTheta>
277 <UpperTheta>90</UpperTheta>
278 </ThetaBounds>
279 </AngleBasisBlock>
280 </AngleBasis>
281 </DataDefinition>
282 ';
283 if ( $doforw ) {
284 print ' <WavelengthData>
285 <LayerNumber>System</LayerNumber>
286 <Wavelength unit="Integral">Visible</Wavelength>
287 <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
288 <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
289 <WavelengthDataBlock>
290 <WavelengthDataDirection>Transmission Front</WavelengthDataDirection>
291 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
292 <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
293 <ScatteringDataType>BTDF</ScatteringDataType>
294 <ScatteringData>
295 ';
296 # Output front transmission (transposed order)
297 for (my $od = 0; $od < $ndiv; $od++) {
298 for (my $id = 0; $id < $ndiv; $id++) {
299 print $tfarr[$ndiv*$id + $od];
300 }
301 print "\n";
302 }
303 print
304 ' </ScatteringData>
305 </WavelengthDataBlock>
306 </WavelengthData>
307 <WavelengthData>
308 <LayerNumber>System</LayerNumber>
309 <Wavelength unit="Integral">Visible</Wavelength>
310 <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
311 <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
312 <WavelengthDataBlock>
313 <WavelengthDataDirection>Reflection Front</WavelengthDataDirection>
314 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
315 <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
316 <ScatteringDataType>BRDF</ScatteringDataType>
317 <ScatteringData>
318 ';
319 # Output front reflection (transposed order)
320 for (my $od = 0; $od < $ndiv; $od++) {
321 for (my $id = 0; $id < $ndiv; $id++) {
322 print $rfarr[$ndiv*$id + $od];
323 }
324 print "\n";
325 }
326 print
327 ' </ScatteringData>
328 </WavelengthDataBlock>
329 </WavelengthData>
330 ';
331 }
332 if ( $doback ) {
333 print ' <WavelengthData>
334 <LayerNumber>System</LayerNumber>
335 <Wavelength unit="Integral">Visible</Wavelength>
336 <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
337 <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
338 <WavelengthDataBlock>
339 <WavelengthDataDirection>Transmission Back</WavelengthDataDirection>
340 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
341 <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
342 <ScatteringDataType>BTDF</ScatteringDataType>
343 <ScatteringData>
344 ';
345 # Output back transmission (transposed order)
346 for (my $od = 0; $od < $ndiv; $od++) {
347 for (my $id = 0; $id < $ndiv; $id++) {
348 print $tbarr[$ndiv*$id + $od];
349 }
350 print "\n";
351 }
352 print
353 ' </ScatteringData>
354 </WavelengthDataBlock>
355 </WavelengthData>
356 <WavelengthData>
357 <LayerNumber>System</LayerNumber>
358 <Wavelength unit="Integral">Visible</Wavelength>
359 <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
360 <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
361 <WavelengthDataBlock>
362 <WavelengthDataDirection>Reflection Back</WavelengthDataDirection>
363 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
364 <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
365 <ScatteringDataType>BRDF</ScatteringDataType>
366 <ScatteringData>
367 ';
368 # Output back reflection (transposed order)
369 for (my $od = 0; $od < $ndiv; $od++) {
370 for (my $id = 0; $id < $ndiv; $id++) {
371 print $rbarr[$ndiv*$id + $od];
372 }
373 print "\n";
374 }
375 print
376 ' </ScatteringData>
377 </WavelengthDataBlock>
378 </WavelengthData>
379 ';
380 }
381 # Output XML epilogue
382 print '</Layer>
383 </Optical>
384 </WindowElement>
385 ';
386 # Clean up temporary files
387 system "rm -rf $td";