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

File Contents

# User Rev Content
1 greg 2.1 #!/usr/bin/perl -w
2 greg 2.13 # RCSid $Id: genBSDF.pl,v 2.12 2011/02/24 20:27:00 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     system "cat $mgfscn";
197     print "xf\n";
198     print "\t\t\t</Geometry>\n";
199     }
200     print ' </Material>
201     <DataDefinition>
202     <IncidentDataStructure>Columns</IncidentDataStructure>
203     <AngleBasis>
204     <AngleBasisName>LBNL/Klems Full</AngleBasisName>
205     <AngleBasisBlock>
206     <Theta>0</Theta>
207     <nPhis>1</nPhis>
208     <ThetaBounds>
209     <LowerTheta>0</LowerTheta>
210     <UpperTheta>5</UpperTheta>
211     </ThetaBounds>
212     </AngleBasisBlock>
213     <AngleBasisBlock>
214     <Theta>10</Theta>
215     <nPhis>8</nPhis>
216     <ThetaBounds>
217     <LowerTheta>5</LowerTheta>
218     <UpperTheta>15</UpperTheta>
219     </ThetaBounds>
220     </AngleBasisBlock>
221     <AngleBasisBlock>
222     <Theta>20</Theta>
223     <nPhis>16</nPhis>
224     <ThetaBounds>
225     <LowerTheta>15</LowerTheta>
226     <UpperTheta>25</UpperTheta>
227     </ThetaBounds>
228     </AngleBasisBlock>
229     <AngleBasisBlock>
230     <Theta>30</Theta>
231     <nPhis>20</nPhis>
232     <ThetaBounds>
233     <LowerTheta>25</LowerTheta>
234     <UpperTheta>35</UpperTheta>
235     </ThetaBounds>
236     </AngleBasisBlock>
237     <AngleBasisBlock>
238     <Theta>40</Theta>
239     <nPhis>24</nPhis>
240     <ThetaBounds>
241     <LowerTheta>35</LowerTheta>
242     <UpperTheta>45</UpperTheta>
243     </ThetaBounds>
244     </AngleBasisBlock>
245     <AngleBasisBlock>
246     <Theta>50</Theta>
247     <nPhis>24</nPhis>
248     <ThetaBounds>
249     <LowerTheta>45</LowerTheta>
250     <UpperTheta>55</UpperTheta>
251     </ThetaBounds>
252     </AngleBasisBlock>
253     <AngleBasisBlock>
254     <Theta>60</Theta>
255     <nPhis>24</nPhis>
256     <ThetaBounds>
257     <LowerTheta>55</LowerTheta>
258     <UpperTheta>65</UpperTheta>
259     </ThetaBounds>
260     </AngleBasisBlock>
261     <AngleBasisBlock>
262     <Theta>70</Theta>
263     <nPhis>16</nPhis>
264     <ThetaBounds>
265     <LowerTheta>65</LowerTheta>
266     <UpperTheta>75</UpperTheta>
267     </ThetaBounds>
268     </AngleBasisBlock>
269     <AngleBasisBlock>
270     <Theta>82.5</Theta>
271     <nPhis>12</nPhis>
272     <ThetaBounds>
273     <LowerTheta>75</LowerTheta>
274     <UpperTheta>90</UpperTheta>
275     </ThetaBounds>
276     </AngleBasisBlock>
277     </AngleBasis>
278     </DataDefinition>
279 greg 2.9 ';
280     if ( $doforw ) {
281     print ' <WavelengthData>
282     <LayerNumber>System</LayerNumber>
283     <Wavelength unit="Integral">Visible</Wavelength>
284     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
285     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
286     <WavelengthDataBlock>
287     <WavelengthDataDirection>Transmission Front</WavelengthDataDirection>
288     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
289     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
290     <ScatteringDataType>BTDF</ScatteringDataType>
291     <ScatteringData>
292     ';
293     # Output front transmission (transposed order)
294     for (my $od = 0; $od < $ndiv; $od++) {
295     for (my $id = 0; $id < $ndiv; $id++) {
296     print $tfarr[$ndiv*$id + $od];
297     }
298     print "\n";
299     }
300     print
301     ' </ScatteringData>
302     </WavelengthDataBlock>
303     </WavelengthData>
304 greg 2.1 <WavelengthData>
305     <LayerNumber>System</LayerNumber>
306     <Wavelength unit="Integral">Visible</Wavelength>
307     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
308     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
309     <WavelengthDataBlock>
310 greg 2.9 <WavelengthDataDirection>Reflection Front</WavelengthDataDirection>
311     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
312     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
313     <ScatteringDataType>BRDF</ScatteringDataType>
314     <ScatteringData>
315     ';
316     # Output front reflection (transposed order)
317     for (my $od = 0; $od < $ndiv; $od++) {
318     for (my $id = 0; $id < $ndiv; $id++) {
319     print $rfarr[$ndiv*$id + $od];
320     }
321     print "\n";
322     }
323     print
324     ' </ScatteringData>
325     </WavelengthDataBlock>
326     </WavelengthData>
327     ';
328     }
329     if ( $doback ) {
330     print ' <WavelengthData>
331     <LayerNumber>System</LayerNumber>
332     <Wavelength unit="Integral">Visible</Wavelength>
333     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
334     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
335     <WavelengthDataBlock>
336     <WavelengthDataDirection>Transmission Back</WavelengthDataDirection>
337 greg 2.1 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
338     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
339     <ScatteringDataType>BTDF</ScatteringDataType>
340     <ScatteringData>
341     ';
342 greg 2.9 # Output back transmission (transposed order)
343 greg 2.3 for (my $od = 0; $od < $ndiv; $od++) {
344     for (my $id = 0; $id < $ndiv; $id++) {
345 greg 2.9 print $tbarr[$ndiv*$id + $od];
346     }
347     print "\n";
348     }
349     print
350     ' </ScatteringData>
351     </WavelengthDataBlock>
352     </WavelengthData>
353     <WavelengthData>
354     <LayerNumber>System</LayerNumber>
355     <Wavelength unit="Integral">Visible</Wavelength>
356     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
357     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
358     <WavelengthDataBlock>
359     <WavelengthDataDirection>Reflection Back</WavelengthDataDirection>
360     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
361     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
362     <ScatteringDataType>BRDF</ScatteringDataType>
363     <ScatteringData>
364     ';
365     # Output back reflection (transposed order)
366     for (my $od = 0; $od < $ndiv; $od++) {
367     for (my $id = 0; $id < $ndiv; $id++) {
368     print $rbarr[$ndiv*$id + $od];
369 greg 2.3 }
370     print "\n";
371     }
372 greg 2.1 print
373     ' </ScatteringData>
374     </WavelengthDataBlock>
375     </WavelengthData>
376 greg 2.9 ';
377     }
378     # Output XML epilogue
379     print '</Layer>
380 greg 2.1 </Optical>
381     </WindowElement>
382     ';
383     # Clean up temporary files
384     system "rm -rf $td";