--- ray/src/util/genBSDF.pl 2010/09/02 02:29:24 2.1 +++ ray/src/util/genBSDF.pl 2011/02/24 20:27:00 2.12 @@ -1,5 +1,5 @@ #!/usr/bin/perl -w -# RCSid $Id: genBSDF.pl,v 2.1 2010/09/02 02:29:24 greg Exp $ +# RCSid $Id: genBSDF.pl,v 2.12 2011/02/24 20:27:00 greg Exp $ # # Compute BSDF based on geometry and material description # @@ -7,22 +7,32 @@ # use strict; sub userror { - print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-dim xmin xmax ymin ymax zmin zmax][{+|-}mgf][{+|-}geom] [input ..]"; + print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-r \"ropts\"][-dim xmin xmax ymin ymax zmin zmax][{+|-}f][{+|-}b][{+|-}mgf][{+|-}geom] [input ..]\n"; exit 1; } my $td = `mktemp -d /tmp/genBSDF.XXXXXX`; chomp $td; my $nsamp = 1000; +my $rtargs = "-w -ab 5 -ad 700 -lw 3e-6"; my $mgfin = 0; my $geout = 1; my $nproc = 1; +my $doforw = 0; +my $doback = 1; my @dim; # Get options while ($#ARGV >= 0) { if ("$ARGV[0]" =~ /^[-+]m/) { $mgfin = ("$ARGV[0]" =~ /^\+/); + } elsif ("$ARGV[0]" eq "-r") { + $rtargs = "$rtargs $ARGV[1]"; + shift @ARGV; } elsif ("$ARGV[0]" =~ /^[-+]g/) { $geout = ("$ARGV[0]" =~ /^\+/); + } elsif ("$ARGV[0]" =~ /^[-+]f/) { + $doforw = ("$ARGV[0]" =~ /^\+/); + } elsif ("$ARGV[0]" =~ /^[-+]b/) { + $doback = ("$ARGV[0]" =~ /^\+/); } elsif ("$ARGV[0]" eq "-c") { $nsamp = $ARGV[1]; shift @ARGV; @@ -31,7 +41,7 @@ while ($#ARGV >= 0) { shift @ARGV; } elsif ("$ARGV[0]" =~ /^-d/) { userror() if ($#ARGV < 6); - @dim = "@ARGV[1..6]"; + @dim = @ARGV[1..6]; shift @ARGV for (1..6); } elsif ("$ARGV[0]" =~ /^[-+]./) { userror(); @@ -40,6 +50,8 @@ while ($#ARGV >= 0) { } shift @ARGV; } +# Check that we're actually being asked to do something +die "Must have at least one of +forward or +backward" if (!$doforw && !$doback); # Get scene description and dimensions my $radscn = "$td/device.rad"; my $mgfscn = "$td/device.mgf"; @@ -54,27 +66,26 @@ if ( $mgfin ) { system "rad2mgf $radscn > $mgfscn" if ( $geout ); } if ($#dim != 5) { - @dim = split /\s+/, `getbbox -h $radscn`; - shift @dim; + @dim = split ' ', `getbbox -h $radscn`; } print STDERR "Warning: Device extends into room\n" if ($dim[5] > 1e-5); -# Add receiver surface (rectangle) -my $modnm="_receiver_black_"; +# Add receiver surfaces (rectangular) +my $fmodnm="receiver_face"; +my $bmodnm="receiver_behind"; open(RADSCN, ">> $radscn"); -print RADSCN "void glow $modnm\n0\n0\n4 0 0 0 0\n\n"; -print RADSCN "$modnm polygon _receiver_\n0\n0\n12\n"; -print RADSCN "\t",$dim[0],"\t",$dim[2],"\t",$dim[5]+1e-5,"\n"; -print RADSCN "\t",$dim[0],"\t",$dim[3],"\t",$dim[5]+1e-5,"\n"; -print RADSCN "\t",$dim[1],"\t",$dim[3],"\t",$dim[5]+1e-5,"\n"; -print RADSCN "\t",$dim[1],"\t",$dim[2],"\t",$dim[5]+1e-5,"\n"; +print RADSCN "void glow $fmodnm\n0\n0\n4 1 1 1 0\n\n"; +print RADSCN "$fmodnm source f_receiver\n0\n0\n4 0 0 1 180\n"; +print RADSCN "void glow $bmodnm\n0\n0\n4 1 1 1 0\n\n"; +print RADSCN "$bmodnm source b_receiver\n0\n0\n4 0 0 -1 180\n"; close RADSCN; # Generate octree system "oconv -w $radscn > $octree"; die "Could not compile scene\n" if ( $? ); -# Set up sampling +# Set up sampling of interior portal # Kbin to produce incident direction in full Klems basis with (x1,x2) randoms my $tcal = ' DEGREE : PI/180; +sq(x) : x*x; Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90); Knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12); Kaccum(r) : if(r-.5, Knaz(r) + Kaccum(r-1), 0); @@ -82,16 +93,20 @@ Kmax : Kaccum(Knaz(0)); Kfindrow(r, rem) : if(rem-Knaz(r)+.5, Kfindrow(r+1, rem-Knaz(r)), r); Krow = if(Kbin-(Kmax-.5), 0, Kfindrow(1, Kbin)); Kcol = Kbin - Kaccum(Krow-1); -Kazi = 360*DEGREE * (Kcol + .5 - x2) / Knaz(Krow); +Kazi = 360*DEGREE * (Kcol + (.5 - x2)) / Knaz(Krow); Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1)); sin_kpol = sin(Kpol); -Dx = -cos(Kazi)*sin_kpol; +Dx = cos(Kazi)*sin_kpol; Dy = sin(Kazi)*sin_kpol; Dz = sqrt(1 - sin_kpol*sin_kpol); +KprojOmega = PI * if(Kbin-.5, + (sq(cos(Kpola(Krow-1)*DEGREE)) - sq(cos(Kpola(Krow)*DEGREE)))/Knaz(Krow), + 1 - sq(cos(Kpola(1)*DEGREE))); '; -# Compute Klems bin from exiting ray direction +# Compute Klems bin from exiting ray direction (forward or backward) my $kcal = ' DEGREE : PI/180; +abs(x) : if(x, x, -x); Acos(x) : 1/DEGREE * if(x-1, 0, if(-1-x, 0, acos(x))); posangle(a) : if(-a, a + 2*PI, a); Atan2(y,x) : 1/DEGREE * posangle(atan2(y,x)); @@ -112,12 +127,52 @@ kbin2(pol,azi) = select(kfindrow(1, pol), kaccum(7) + kazn(azi,360/knaz(8)), kaccum(8) + kazn(azi,360/knaz(9)) ); -kbin = kbin2(Acos(Dz), Atan2(Dy, -Dx)); +kbin = kbin2(Acos(abs(Dz)),Atan2(Dy,Dx)); '; my $ndiv = 145; my $nx = int(sqrt($nsamp*($dim[1]-$dim[0])/($dim[3]-$dim[2])) + .5); my $ny = int($nsamp/$nx + .5); $nsamp = $nx * $ny; +# Compute scattering data using rtcontrib +my @tfarr; +my @rfarr; +my @tbarr; +my @rbarr; +my $cmd; +my $rtcmd = "rtcontrib -h -ff -fo -n $nproc -c $nsamp " . + "-e '$kcal' -b kbin -bn $ndiv " . + "-o '$td/%s.flt' -m $fmodnm -m $bmodnm $rtargs $octree"; +my $rccmd = "rcalc -e '$tcal' " . + "-e 'mod(n,d):n-floor(n/d)*d' -e 'Kbin=mod(recno-.999,$ndiv)' " . + q{-if3 -e '$1=(0.265*$1+0.670*$2+0.065*$3)/KprojOmega'}; +if ( $doforw ) { +$cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " . + "-e 'xp=(\$3+rand(.12*recno+288))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " . + "-e 'yp=(\$2+rand(.37*recno-44))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " . + "-e 'zp:$dim[4]' " . + q{-e 'Kbin=$1;x1=rand(2.75*recno+3.1);x2=rand(-2.01*recno-3.37)' } . + q{-e '$1=xp-Dx;$2=yp-Dy;$3=zp-Dz;$4=Dx;$5=Dy;$6=Dz' } . + "| $rtcmd"; +system "$cmd" || die "Failure running: $cmd\n"; +@tfarr = `$rccmd $td/$fmodnm.flt`; +die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? ); +@rfarr = `$rccmd $td/$bmodnm.flt`; +die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? ); +} +if ( $doback ) { +$cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " . + "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " . + "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " . + "-e 'zp:$dim[5]' " . + q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } . + q{-e '$1=xp-Dx;$2=yp-Dy;$3=zp+Dz;$4=Dx;$5=Dy;$6=-Dz' } . + "| $rtcmd"; +system "$cmd" || die "Failure running: $cmd\n"; +@tbarr = `$rccmd $td/$bmodnm.flt`; +die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? ); +@rbarr = `$rccmd $td/$fmodnm.flt`; +die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? ); +} # Output XML prologue print ' @@ -220,7 +275,9 @@ print ' - +'; +if ( $doforw ) { +print ' System Visible CIE Illuminant D65 1nm.ssp @@ -232,22 +289,93 @@ print ' BTDF '; -# Compute actual scattering data using rtcontrib -system "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " . - "-e 'xp=(\$3+rand(.35*recno-15))*(($dim[1]-$dim[0])/$nx)+$dim[0]' " . - "-e 'yp=(\$2+rand(.86*recno+11))*(($dim[3]-$dim[2])/$ny)+$dim[2]' " . - "-e 'zp:$dim[4]-1e-5' " . - q{-e 'Kbin=$1;x1=rand(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } . - q{-e '$1=xp;$2=yp;$3=zp;$4=Dx;$5=Dy;$6=Dz' } . - "| rtcontrib -h -ff -n $nproc -c $nsamp -e '$kcal' -b kbin -bn $ndiv " . - "-m $modnm -w -ab 4 -lw 1e-5 $octree " . - q{| rcalc -if3 -e '$1=0.265*$1+0.670*$2+0.065*$3'}; -# Output XML epilogue +# Output front transmission (transposed order) +for (my $od = 0; $od < $ndiv; $od++) { + for (my $id = 0; $id < $ndiv; $id++) { + print $tfarr[$ndiv*$id + $od]; + } + print "\n"; +} print ' - + + System + Visible + CIE Illuminant D65 1nm.ssp + ASTM E308 1931 Y.dsp + + Reflection Front + LBNL/Klems Full + LBNL/Klems Full + BRDF + +'; +# Output front reflection (transposed order) +for (my $od = 0; $od < $ndiv; $od++) { + for (my $id = 0; $id < $ndiv; $id++) { + print $rfarr[$ndiv*$id + $od]; + } + print "\n"; +} +print +' + + +'; +} +if ( $doback ) { +print ' + System + Visible + CIE Illuminant D65 1nm.ssp + ASTM E308 1931 Y.dsp + + Transmission Back + LBNL/Klems Full + LBNL/Klems Full + BTDF + +'; +# Output back transmission (transposed order) +for (my $od = 0; $od < $ndiv; $od++) { + for (my $id = 0; $id < $ndiv; $id++) { + print $tbarr[$ndiv*$id + $od]; + } + print "\n"; +} +print +' + + + + System + Visible + CIE Illuminant D65 1nm.ssp + ASTM E308 1931 Y.dsp + + Reflection Back + LBNL/Klems Full + LBNL/Klems Full + BRDF + +'; +# Output back reflection (transposed order) +for (my $od = 0; $od < $ndiv; $od++) { + for (my $id = 0; $id < $ndiv; $id++) { + print $rbarr[$ndiv*$id + $od]; + } + print "\n"; +} +print +' + + +'; +} +# Output XML epilogue +print ' ';