ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genBSDF.pl
Revision: 2.12
Committed: Thu Feb 24 20:27:00 2011 UTC (13 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +2 -2 lines
Log Message:
Fixed usage message

File Contents

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