ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genBSDF.pl
Revision: 2.14
Committed: Sat Apr 16 01:13:22 2011 UTC (13 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.13: +4 -2 lines
Log Message:
Eliminated Unix dependencies

File Contents

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