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

# User Rev Content
1 greg 2.1 #!/usr/bin/perl -w
2 greg 2.9 # RCSid $Id: genBSDF.pl,v 2.8 2011/01/29 16:44:01 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     my $mgfin = 0;
17     my $geout = 1;
18     my $nproc = 1;
19 greg 2.9 my $doforw = 0;
20     my $doback = 1;
21 greg 2.1 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 greg 2.9 } elsif ("$ARGV[0]" =~ /^[-+]f/) {
29     $doforw = ("$ARGV[0]" =~ /^\+/);
30     } elsif ("$ARGV[0]" =~ /^[-+]b/) {
31     $doback = ("$ARGV[0]" =~ /^\+/);
32 greg 2.1 } 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 greg 2.8 @dim = @ARGV[1..6];
41 greg 2.1 shift @ARGV for (1..6);
42     } elsif ("$ARGV[0]" =~ /^[-+]./) {
43     userror();
44     } else {
45     last;
46     }
47     shift @ARGV;
48     }
49 greg 2.9 # 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 greg 2.1 # 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 greg 2.7 @dim = split ' ', `getbbox -h $radscn`;
66 greg 2.1 }
67     print STDERR "Warning: Device extends into room\n" if ($dim[5] > 1e-5);
68 greg 2.9 # Add receiver surfaces (rectangular)
69     my $bmodnm="receiver_behind";
70     my $fmodnm="receiver_face";
71 greg 2.1 open(RADSCN, ">> $radscn");
72 greg 2.9 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 greg 2.1 close RADSCN;
85     # Generate octree
86     system "oconv -w $radscn > $octree";
87     die "Could not compile scene\n" if ( $? );
88 greg 2.9 # Set up sampling of interior portal
89 greg 2.1 # Kbin to produce incident direction in full Klems basis with (x1,x2) randoms
90     my $tcal = '
91     DEGREE : PI/180;
92 greg 2.5 sq(x) : x*x;
93 greg 2.1 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 greg 2.2 Kazi = 360*DEGREE * (Kcol + (.5 - x2)) / Knaz(Krow);
101 greg 2.1 Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1));
102     sin_kpol = sin(Kpol);
103 greg 2.9 Dx = cos(Kazi)*sin_kpol;
104 greg 2.1 Dy = sin(Kazi)*sin_kpol;
105     Dz = sqrt(1 - sin_kpol*sin_kpol);
106 greg 2.5 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 greg 2.1 ';
110 greg 2.9 # Compute Klems bin from exiting ray direction (forward or backward)
111 greg 2.1 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 greg 2.9 kbin = if(Dz, kbin2(Acos(Dz),Atan2(Dy,Dx)), kbin2(Acos(-Dz),Atan2(-Dy,-Dx)));
134 greg 2.1 ';
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 greg 2.3 # Compute scattering data using rtcontrib
140 greg 2.9 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 greg 2.3 "-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 greg 2.9 "| $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 greg 2.1 # 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 greg 2.9 ';
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 greg 2.1 <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 greg 2.9 <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 greg 2.1 <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
341     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
342     <ScatteringDataType>BTDF</ScatteringDataType>
343     <ScatteringData>
344     ';
345 greg 2.9 # Output back transmission (transposed order)
346 greg 2.3 for (my $od = 0; $od < $ndiv; $od++) {
347     for (my $id = 0; $id < $ndiv; $id++) {
348 greg 2.9 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 greg 2.3 }
373     print "\n";
374     }
375 greg 2.1 print
376     ' </ScatteringData>
377     </WavelengthDataBlock>
378     </WavelengthData>
379 greg 2.9 ';
380     }
381     # Output XML epilogue
382     print '</Layer>
383 greg 2.1 </Optical>
384     </WindowElement>
385     ';
386     # Clean up temporary files
387     system "rm -rf $td";