ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genBSDF.pl
Revision: 2.10
Committed: Tue Feb 22 22:51:23 2011 UTC (13 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +14 -19 lines
Log Message:
Added -r option and fixed ray leakage (thanks to David G-M for help)

File Contents

# User Rev Content
1 greg 2.1 #!/usr/bin/perl -w
2 greg 2.10 # RCSid $Id: genBSDF.pl,v 2.9 2011/02/21 22:48:51 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.9 print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-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     Acos(x) : 1/DEGREE * if(x-1, 0, if(-1-x, 0, acos(x)));
110     posangle(a) : if(-a, a + 2*PI, a);
111     Atan2(y,x) : 1/DEGREE * posangle(atan2(y,x));
112     kpola(r) : select(r, 5, 15, 25, 35, 45, 55, 65, 75, 90);
113     knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
114     kaccum(r) : if(r-.5, knaz(r) + kaccum(r-1), 0);
115     kfindrow(r, pol) : if(r-kpola(0)+.5, r,
116     if(pol-kpola(r), kfindrow(r+1, pol), r) );
117     kazn(azi,inc) : if((360-.5*inc)-azi, floor((azi+.5*inc)/inc), 0);
118     kbin2(pol,azi) = select(kfindrow(1, pol),
119     kazn(azi,360/knaz(1)),
120     kaccum(1) + kazn(azi,360/knaz(2)),
121     kaccum(2) + kazn(azi,360/knaz(3)),
122     kaccum(3) + kazn(azi,360/knaz(4)),
123     kaccum(4) + kazn(azi,360/knaz(5)),
124     kaccum(5) + kazn(azi,360/knaz(6)),
125     kaccum(6) + kazn(azi,360/knaz(7)),
126     kaccum(7) + kazn(azi,360/knaz(8)),
127     kaccum(8) + kazn(azi,360/knaz(9))
128     );
129 greg 2.9 kbin = if(Dz, kbin2(Acos(Dz),Atan2(Dy,Dx)), kbin2(Acos(-Dz),Atan2(-Dy,-Dx)));
130 greg 2.1 ';
131     my $ndiv = 145;
132     my $nx = int(sqrt($nsamp*($dim[1]-$dim[0])/($dim[3]-$dim[2])) + .5);
133     my $ny = int($nsamp/$nx + .5);
134     $nsamp = $nx * $ny;
135 greg 2.3 # Compute scattering data using rtcontrib
136 greg 2.9 my @tfarr;
137     my @rfarr;
138     my @tbarr;
139     my @rbarr;
140     my $cmd;
141     my $rtcmd = "rtcontrib -h -ff -fo -n $nproc -c $nsamp " .
142     "-e '$kcal' -b kbin -bn $ndiv " .
143     "-o '$td/%s.flt' -m $fmodnm -m $bmodnm $rtargs $octree";
144     my $rccmd = "rcalc -e '$tcal' " .
145     "-e 'mod(n,d):n-floor(n/d)*d' -e 'Kbin=mod(recno-.999,$ndiv)' " .
146     q{-if3 -e '$1=(0.265*$1+0.670*$2+0.065*$3)/KprojOmega'};
147     if ( $doforw ) {
148     $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " .
149 greg 2.3 "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " .
150     "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " .
151 greg 2.10 "-e 'zp:$dim[4]' " .
152 greg 2.3 q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } .
153 greg 2.10 q{-e '$1=xp-Dx;$2=yp-Dy;$3=zp-Dz;$4=Dx;$5=Dy;$6=Dz' } .
154 greg 2.9 "| $rtcmd";
155     system "$cmd" || die "Failure running: $cmd\n";
156     @tfarr = `$rccmd $td/$fmodnm.flt`;
157     die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? );
158     @rfarr = `$rccmd $td/$bmodnm.flt`;
159     die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? );
160     }
161     if ( $doback ) {
162     $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " .
163     "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " .
164     "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " .
165 greg 2.10 "-e 'zp:$dim[5]' " .
166 greg 2.9 q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } .
167 greg 2.10 q{-e '$1=xp+Dx;$2=yp+Dy;$3=zp+Dz;$4=-Dx;$5=-Dy;$6=-Dz' } .
168 greg 2.9 "| $rtcmd";
169     system "$cmd" || die "Failure running: $cmd\n";
170     @tbarr = `$rccmd $td/$bmodnm.flt`;
171     die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? );
172     @rbarr = `$rccmd $td/$fmodnm.flt`;
173     die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? );
174     }
175 greg 2.1 # Output XML prologue
176     print
177     '<?xml version="1.0" encoding="UTF-8"?>
178     <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">
179     <WindowElementType>System</WindowElementType>
180     <Optical>
181     <Layer>
182     <Material>
183     <Name>Name</Name>
184     <Manufacturer>Manufacturer</Manufacturer>
185     ';
186     printf "\t\t\t<Thickness unit=\"Meter\">%.3f</Thickness>\n", $dim[5] - $dim[4];
187     printf "\t\t\t<Width unit=\"Meter\">%.3f</Width>\n", $dim[1] - $dim[0];
188     printf "\t\t\t<Height unit=\"Meter\">%.3f</Height>\n", $dim[3] - $dim[2];
189     print "\t\t\t<DeviceType>Integral</DeviceType>\n";
190     # Output MGF description if requested
191     if ( $geout ) {
192     print "\t\t\t<Geometry format=\"MGF\" unit=\"Meter\">\n";
193     printf "xf -t %.6f %.6f 0\n", -($dim[0]+$dim[1])/2, -($dim[2]+$dim[3])/2;
194     system "cat $mgfscn";
195     print "xf\n";
196     print "\t\t\t</Geometry>\n";
197     }
198     print ' </Material>
199     <DataDefinition>
200     <IncidentDataStructure>Columns</IncidentDataStructure>
201     <AngleBasis>
202     <AngleBasisName>LBNL/Klems Full</AngleBasisName>
203     <AngleBasisBlock>
204     <Theta>0</Theta>
205     <nPhis>1</nPhis>
206     <ThetaBounds>
207     <LowerTheta>0</LowerTheta>
208     <UpperTheta>5</UpperTheta>
209     </ThetaBounds>
210     </AngleBasisBlock>
211     <AngleBasisBlock>
212     <Theta>10</Theta>
213     <nPhis>8</nPhis>
214     <ThetaBounds>
215     <LowerTheta>5</LowerTheta>
216     <UpperTheta>15</UpperTheta>
217     </ThetaBounds>
218     </AngleBasisBlock>
219     <AngleBasisBlock>
220     <Theta>20</Theta>
221     <nPhis>16</nPhis>
222     <ThetaBounds>
223     <LowerTheta>15</LowerTheta>
224     <UpperTheta>25</UpperTheta>
225     </ThetaBounds>
226     </AngleBasisBlock>
227     <AngleBasisBlock>
228     <Theta>30</Theta>
229     <nPhis>20</nPhis>
230     <ThetaBounds>
231     <LowerTheta>25</LowerTheta>
232     <UpperTheta>35</UpperTheta>
233     </ThetaBounds>
234     </AngleBasisBlock>
235     <AngleBasisBlock>
236     <Theta>40</Theta>
237     <nPhis>24</nPhis>
238     <ThetaBounds>
239     <LowerTheta>35</LowerTheta>
240     <UpperTheta>45</UpperTheta>
241     </ThetaBounds>
242     </AngleBasisBlock>
243     <AngleBasisBlock>
244     <Theta>50</Theta>
245     <nPhis>24</nPhis>
246     <ThetaBounds>
247     <LowerTheta>45</LowerTheta>
248     <UpperTheta>55</UpperTheta>
249     </ThetaBounds>
250     </AngleBasisBlock>
251     <AngleBasisBlock>
252     <Theta>60</Theta>
253     <nPhis>24</nPhis>
254     <ThetaBounds>
255     <LowerTheta>55</LowerTheta>
256     <UpperTheta>65</UpperTheta>
257     </ThetaBounds>
258     </AngleBasisBlock>
259     <AngleBasisBlock>
260     <Theta>70</Theta>
261     <nPhis>16</nPhis>
262     <ThetaBounds>
263     <LowerTheta>65</LowerTheta>
264     <UpperTheta>75</UpperTheta>
265     </ThetaBounds>
266     </AngleBasisBlock>
267     <AngleBasisBlock>
268     <Theta>82.5</Theta>
269     <nPhis>12</nPhis>
270     <ThetaBounds>
271     <LowerTheta>75</LowerTheta>
272     <UpperTheta>90</UpperTheta>
273     </ThetaBounds>
274     </AngleBasisBlock>
275     </AngleBasis>
276     </DataDefinition>
277 greg 2.9 ';
278     if ( $doforw ) {
279     print ' <WavelengthData>
280     <LayerNumber>System</LayerNumber>
281     <Wavelength unit="Integral">Visible</Wavelength>
282     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
283     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
284     <WavelengthDataBlock>
285     <WavelengthDataDirection>Transmission Front</WavelengthDataDirection>
286     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
287     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
288     <ScatteringDataType>BTDF</ScatteringDataType>
289     <ScatteringData>
290     ';
291     # Output front transmission (transposed order)
292     for (my $od = 0; $od < $ndiv; $od++) {
293     for (my $id = 0; $id < $ndiv; $id++) {
294     print $tfarr[$ndiv*$id + $od];
295     }
296     print "\n";
297     }
298     print
299     ' </ScatteringData>
300     </WavelengthDataBlock>
301     </WavelengthData>
302 greg 2.1 <WavelengthData>
303     <LayerNumber>System</LayerNumber>
304     <Wavelength unit="Integral">Visible</Wavelength>
305     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
306     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
307     <WavelengthDataBlock>
308 greg 2.9 <WavelengthDataDirection>Reflection Front</WavelengthDataDirection>
309     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
310     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
311     <ScatteringDataType>BRDF</ScatteringDataType>
312     <ScatteringData>
313     ';
314     # Output front reflection (transposed order)
315     for (my $od = 0; $od < $ndiv; $od++) {
316     for (my $id = 0; $id < $ndiv; $id++) {
317     print $rfarr[$ndiv*$id + $od];
318     }
319     print "\n";
320     }
321     print
322     ' </ScatteringData>
323     </WavelengthDataBlock>
324     </WavelengthData>
325     ';
326     }
327     if ( $doback ) {
328     print ' <WavelengthData>
329     <LayerNumber>System</LayerNumber>
330     <Wavelength unit="Integral">Visible</Wavelength>
331     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
332     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
333     <WavelengthDataBlock>
334     <WavelengthDataDirection>Transmission Back</WavelengthDataDirection>
335 greg 2.1 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
336     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
337     <ScatteringDataType>BTDF</ScatteringDataType>
338     <ScatteringData>
339     ';
340 greg 2.9 # Output back transmission (transposed order)
341 greg 2.3 for (my $od = 0; $od < $ndiv; $od++) {
342     for (my $id = 0; $id < $ndiv; $id++) {
343 greg 2.9 print $tbarr[$ndiv*$id + $od];
344     }
345     print "\n";
346     }
347     print
348     ' </ScatteringData>
349     </WavelengthDataBlock>
350     </WavelengthData>
351     <WavelengthData>
352     <LayerNumber>System</LayerNumber>
353     <Wavelength unit="Integral">Visible</Wavelength>
354     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
355     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
356     <WavelengthDataBlock>
357     <WavelengthDataDirection>Reflection Back</WavelengthDataDirection>
358     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
359     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
360     <ScatteringDataType>BRDF</ScatteringDataType>
361     <ScatteringData>
362     ';
363     # Output back reflection (transposed order)
364     for (my $od = 0; $od < $ndiv; $od++) {
365     for (my $id = 0; $id < $ndiv; $id++) {
366     print $rbarr[$ndiv*$id + $od];
367 greg 2.3 }
368     print "\n";
369     }
370 greg 2.1 print
371     ' </ScatteringData>
372     </WavelengthDataBlock>
373     </WavelengthData>
374 greg 2.9 ';
375     }
376     # Output XML epilogue
377     print '</Layer>
378 greg 2.1 </Optical>
379     </WindowElement>
380     ';
381     # Clean up temporary files
382     system "rm -rf $td";