ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/genBSDF.pl
Revision: 2.5
Committed: Thu Dec 9 17:00:43 2010 UTC (13 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +7 -6 lines
Log Message:
Slight correction to averaged cosine (thanks to Jacob)

File Contents

# User Rev Content
1 greg 2.1 #!/usr/bin/perl -w
2 greg 2.5 # RCSid $Id: genBSDF.pl,v 2.4 2010/12/09 06:14:04 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.3 print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-dim xmin xmax ymin ymax zmin zmax][{+|-}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     my @dim;
20     # Get options
21     while ($#ARGV >= 0) {
22     if ("$ARGV[0]" =~ /^[-+]m/) {
23     $mgfin = ("$ARGV[0]" =~ /^\+/);
24     } elsif ("$ARGV[0]" =~ /^[-+]g/) {
25     $geout = ("$ARGV[0]" =~ /^\+/);
26     } elsif ("$ARGV[0]" eq "-c") {
27     $nsamp = $ARGV[1];
28     shift @ARGV;
29     } elsif ("$ARGV[0]" eq "-n") {
30     $nproc = $ARGV[1];
31     shift @ARGV;
32     } elsif ("$ARGV[0]" =~ /^-d/) {
33     userror() if ($#ARGV < 6);
34     @dim = "@ARGV[1..6]";
35     shift @ARGV for (1..6);
36     } elsif ("$ARGV[0]" =~ /^[-+]./) {
37     userror();
38     } else {
39     last;
40     }
41     shift @ARGV;
42     }
43     # Get scene description and dimensions
44     my $radscn = "$td/device.rad";
45     my $mgfscn = "$td/device.mgf";
46     my $octree = "$td/device.oct";
47     if ( $mgfin ) {
48     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";
49     die "Could not load MGF input\n" if ( $? );
50     system "mgf2rad $mgfscn > $radscn";
51     } else {
52     system "cat @ARGV | xform -e > $radscn";
53     die "Could not load Radiance input\n" if ( $? );
54     system "rad2mgf $radscn > $mgfscn" if ( $geout );
55     }
56     if ($#dim != 5) {
57     @dim = split /\s+/, `getbbox -h $radscn`;
58     shift @dim;
59     }
60     print STDERR "Warning: Device extends into room\n" if ($dim[5] > 1e-5);
61     # Add receiver surface (rectangle)
62     my $modnm="_receiver_black_";
63     open(RADSCN, ">> $radscn");
64     print RADSCN "void glow $modnm\n0\n0\n4 0 0 0 0\n\n";
65     print RADSCN "$modnm polygon _receiver_\n0\n0\n12\n";
66     print RADSCN "\t",$dim[0],"\t",$dim[2],"\t",$dim[5]+1e-5,"\n";
67     print RADSCN "\t",$dim[0],"\t",$dim[3],"\t",$dim[5]+1e-5,"\n";
68     print RADSCN "\t",$dim[1],"\t",$dim[3],"\t",$dim[5]+1e-5,"\n";
69     print RADSCN "\t",$dim[1],"\t",$dim[2],"\t",$dim[5]+1e-5,"\n";
70     close RADSCN;
71     # Generate octree
72     system "oconv -w $radscn > $octree";
73     die "Could not compile scene\n" if ( $? );
74     # Set up sampling
75     # Kbin to produce incident direction in full Klems basis with (x1,x2) randoms
76     my $tcal = '
77     DEGREE : PI/180;
78 greg 2.5 sq(x) : x*x;
79 greg 2.1 Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90);
80     Knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
81     Kaccum(r) : if(r-.5, Knaz(r) + Kaccum(r-1), 0);
82     Kmax : Kaccum(Knaz(0));
83     Kfindrow(r, rem) : if(rem-Knaz(r)+.5, Kfindrow(r+1, rem-Knaz(r)), r);
84     Krow = if(Kbin-(Kmax-.5), 0, Kfindrow(1, Kbin));
85     Kcol = Kbin - Kaccum(Krow-1);
86 greg 2.2 Kazi = 360*DEGREE * (Kcol + (.5 - x2)) / Knaz(Krow);
87 greg 2.1 Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1));
88     sin_kpol = sin(Kpol);
89     Dx = -cos(Kazi)*sin_kpol;
90     Dy = sin(Kazi)*sin_kpol;
91     Dz = sqrt(1 - sin_kpol*sin_kpol);
92 greg 2.5 KprojOmega = PI * if(Kbin-.5,
93     (sq(cos(Kpola(Krow-1)*DEGREE)) - sq(cos(Kpola(Krow)*DEGREE)))/Knaz(Krow),
94     1 - sq(cos(Kpola(1)*DEGREE)));
95 greg 2.1 ';
96     # Compute Klems bin from exiting ray direction
97     my $kcal = '
98     DEGREE : PI/180;
99     Acos(x) : 1/DEGREE * if(x-1, 0, if(-1-x, 0, acos(x)));
100     posangle(a) : if(-a, a + 2*PI, a);
101     Atan2(y,x) : 1/DEGREE * posangle(atan2(y,x));
102     kpola(r) : select(r, 5, 15, 25, 35, 45, 55, 65, 75, 90);
103     knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
104     kaccum(r) : if(r-.5, knaz(r) + kaccum(r-1), 0);
105     kfindrow(r, pol) : if(r-kpola(0)+.5, r,
106     if(pol-kpola(r), kfindrow(r+1, pol), r) );
107     kazn(azi,inc) : if((360-.5*inc)-azi, floor((azi+.5*inc)/inc), 0);
108     kbin2(pol,azi) = select(kfindrow(1, pol),
109     kazn(azi,360/knaz(1)),
110     kaccum(1) + kazn(azi,360/knaz(2)),
111     kaccum(2) + kazn(azi,360/knaz(3)),
112     kaccum(3) + kazn(azi,360/knaz(4)),
113     kaccum(4) + kazn(azi,360/knaz(5)),
114     kaccum(5) + kazn(azi,360/knaz(6)),
115     kaccum(6) + kazn(azi,360/knaz(7)),
116     kaccum(7) + kazn(azi,360/knaz(8)),
117     kaccum(8) + kazn(azi,360/knaz(9))
118     );
119     kbin = kbin2(Acos(Dz), Atan2(Dy, -Dx));
120     ';
121     my $ndiv = 145;
122     my $nx = int(sqrt($nsamp*($dim[1]-$dim[0])/($dim[3]-$dim[2])) + .5);
123     my $ny = int($nsamp/$nx + .5);
124     $nsamp = $nx * $ny;
125 greg 2.3 # Compute scattering data using rtcontrib
126     my $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " .
127     "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " .
128     "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " .
129     "-e 'zp:$dim[4]-1e-5' " .
130     q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } .
131     q{-e '$1=xp;$2=yp;$3=zp;$4=Dx;$5=Dy;$6=Dz' } .
132     "| rtcontrib -h -ff -n $nproc -c $nsamp -e '$kcal' -b kbin -bn $ndiv " .
133     "-m $modnm -w -ab 5 -ad 700 -lw 3e-6 $octree " .
134 greg 2.5 "| rcalc -e '$tcal' " .
135 greg 2.4 "-e 'mod(n,d):n-floor(n/d)*d' -e 'Kbin=mod(recno-1,$ndiv)' " .
136 greg 2.5 q{-if3 -e '$1=(0.265*$1+0.670*$2+0.065*$3)/KprojOmega'};
137 greg 2.3 my @darr = `$cmd`;
138     die "Failure running: $cmd\n" if ( $? );
139 greg 2.1 # Output XML prologue
140     print
141     '<?xml version="1.0" encoding="UTF-8"?>
142     <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">
143     <WindowElementType>System</WindowElementType>
144     <Optical>
145     <Layer>
146     <Material>
147     <Name>Name</Name>
148     <Manufacturer>Manufacturer</Manufacturer>
149     ';
150     printf "\t\t\t<Thickness unit=\"Meter\">%.3f</Thickness>\n", $dim[5] - $dim[4];
151     printf "\t\t\t<Width unit=\"Meter\">%.3f</Width>\n", $dim[1] - $dim[0];
152     printf "\t\t\t<Height unit=\"Meter\">%.3f</Height>\n", $dim[3] - $dim[2];
153     print "\t\t\t<DeviceType>Integral</DeviceType>\n";
154     # Output MGF description if requested
155     if ( $geout ) {
156     print "\t\t\t<Geometry format=\"MGF\" unit=\"Meter\">\n";
157     printf "xf -t %.6f %.6f 0\n", -($dim[0]+$dim[1])/2, -($dim[2]+$dim[3])/2;
158     system "cat $mgfscn";
159     print "xf\n";
160     print "\t\t\t</Geometry>\n";
161     }
162     print ' </Material>
163     <DataDefinition>
164     <IncidentDataStructure>Columns</IncidentDataStructure>
165     <AngleBasis>
166     <AngleBasisName>LBNL/Klems Full</AngleBasisName>
167     <AngleBasisBlock>
168     <Theta>0</Theta>
169     <nPhis>1</nPhis>
170     <ThetaBounds>
171     <LowerTheta>0</LowerTheta>
172     <UpperTheta>5</UpperTheta>
173     </ThetaBounds>
174     </AngleBasisBlock>
175     <AngleBasisBlock>
176     <Theta>10</Theta>
177     <nPhis>8</nPhis>
178     <ThetaBounds>
179     <LowerTheta>5</LowerTheta>
180     <UpperTheta>15</UpperTheta>
181     </ThetaBounds>
182     </AngleBasisBlock>
183     <AngleBasisBlock>
184     <Theta>20</Theta>
185     <nPhis>16</nPhis>
186     <ThetaBounds>
187     <LowerTheta>15</LowerTheta>
188     <UpperTheta>25</UpperTheta>
189     </ThetaBounds>
190     </AngleBasisBlock>
191     <AngleBasisBlock>
192     <Theta>30</Theta>
193     <nPhis>20</nPhis>
194     <ThetaBounds>
195     <LowerTheta>25</LowerTheta>
196     <UpperTheta>35</UpperTheta>
197     </ThetaBounds>
198     </AngleBasisBlock>
199     <AngleBasisBlock>
200     <Theta>40</Theta>
201     <nPhis>24</nPhis>
202     <ThetaBounds>
203     <LowerTheta>35</LowerTheta>
204     <UpperTheta>45</UpperTheta>
205     </ThetaBounds>
206     </AngleBasisBlock>
207     <AngleBasisBlock>
208     <Theta>50</Theta>
209     <nPhis>24</nPhis>
210     <ThetaBounds>
211     <LowerTheta>45</LowerTheta>
212     <UpperTheta>55</UpperTheta>
213     </ThetaBounds>
214     </AngleBasisBlock>
215     <AngleBasisBlock>
216     <Theta>60</Theta>
217     <nPhis>24</nPhis>
218     <ThetaBounds>
219     <LowerTheta>55</LowerTheta>
220     <UpperTheta>65</UpperTheta>
221     </ThetaBounds>
222     </AngleBasisBlock>
223     <AngleBasisBlock>
224     <Theta>70</Theta>
225     <nPhis>16</nPhis>
226     <ThetaBounds>
227     <LowerTheta>65</LowerTheta>
228     <UpperTheta>75</UpperTheta>
229     </ThetaBounds>
230     </AngleBasisBlock>
231     <AngleBasisBlock>
232     <Theta>82.5</Theta>
233     <nPhis>12</nPhis>
234     <ThetaBounds>
235     <LowerTheta>75</LowerTheta>
236     <UpperTheta>90</UpperTheta>
237     </ThetaBounds>
238     </AngleBasisBlock>
239     </AngleBasis>
240     </DataDefinition>
241     <WavelengthData>
242     <LayerNumber>System</LayerNumber>
243     <Wavelength unit="Integral">Visible</Wavelength>
244     <SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>
245     <DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>
246     <WavelengthDataBlock>
247     <WavelengthDataDirection>Transmission Front</WavelengthDataDirection>
248     <ColumnAngleBasis>LBNL/Klems Full</ColumnAngleBasis>
249     <RowAngleBasis>LBNL/Klems Full</RowAngleBasis>
250     <ScatteringDataType>BTDF</ScatteringDataType>
251     <ScatteringData>
252     ';
253 greg 2.3 # Output computed data (transposed order)
254     for (my $od = 0; $od < $ndiv; $od++) {
255     for (my $id = 0; $id < $ndiv; $id++) {
256     print $darr[$ndiv*$id + $od];
257     }
258     print "\n";
259     }
260 greg 2.1 # Output XML epilogue
261     print
262     ' </ScatteringData>
263     </WavelengthDataBlock>
264     </WavelengthData>
265     </Layer>
266     </Optical>
267     </WindowElement>
268     ';
269     # Clean up temporary files
270     system "rm -rf $td";