--- ray/src/util/genBSDF.pl 2011/02/22 22:51:23 2.10 +++ ray/src/util/genBSDF.pl 2012/09/10 18:06:08 2.39 @@ -1,24 +1,30 @@ #!/usr/bin/perl -w -# RCSid $Id: genBSDF.pl,v 2.10 2011/02/22 22:51:23 greg Exp $ +# RCSid $Id: genBSDF.pl,v 2.39 2012/09/10 18:06:08 greg Exp $ # # Compute BSDF based on geometry and material description # # G. Ward # use strict; +use File::Temp qw/ :mktemp /; sub userror { - print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-dim xmin xmax ymin ymax zmin zmax][{+|-}f][{+|-}b][{+|-}mgf][{+|-}geom] [input ..]\n"; + print STDERR "Usage: genBSDF [-n Nproc][-c Nsamp][-t{3|4} Nlog2][-r \"ropts\"][-dim xmin xmax ymin ymax zmin zmax][{+|-}f][{+|-}b][{+|-}mgf][{+|-}geom units] [input ..]\n"; exit 1; } -my $td = `mktemp -d /tmp/genBSDF.XXXXXX`; +my $td = mkdtemp("/tmp/genBSDF.XXXXXX"); chomp $td; -my $nsamp = 1000; +my @savedARGV = @ARGV; +my $tensortree = 0; +my $ttlog2 = 4; +my $nsamp = 2000; 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 $pctcull = 90; +my $gunit = "Meter"; my @dim; # Get options while ($#ARGV >= 0) { @@ -29,10 +35,23 @@ while ($#ARGV >= 0) { shift @ARGV; } elsif ("$ARGV[0]" =~ /^[-+]g/) { $geout = ("$ARGV[0]" =~ /^\+/); + $gunit = $ARGV[1]; + if ($gunit !~ /^(?i)(meter|foot|inch|centimeter|millimeter)$/) { + die "Illegal geometry unit '$gunit': must be meter, foot, inch, centimeter, or millimeter\n"; + } + shift @ARGV; } elsif ("$ARGV[0]" =~ /^[-+]f/) { $doforw = ("$ARGV[0]" =~ /^\+/); } elsif ("$ARGV[0]" =~ /^[-+]b/) { $doback = ("$ARGV[0]" =~ /^\+/); + } elsif ("$ARGV[0]" eq "-t") { + # Use value < 0 for rttree_reduce bypass + $pctcull = $ARGV[1]; + shift @ARGV; + } elsif ("$ARGV[0]" =~ /^-t[34]$/) { + $tensortree = substr($ARGV[0], 2, 1); + $ttlog2 = $ARGV[1]; + shift @ARGV; } elsif ("$ARGV[0]" eq "-c") { $nsamp = $ARGV[1]; shift @ARGV; @@ -51,7 +70,7 @@ 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); +die "Must have at least one of +forward or +backward\n" if (!$doforw && !$doback); # Get scene description and dimensions my $radscn = "$td/device.rad"; my $mgfscn = "$td/device.mgf"; @@ -61,7 +80,7 @@ if ( $mgfin ) { die "Could not load MGF input\n" if ( $? ); system "mgf2rad $mgfscn > $radscn"; } else { - system "cat @ARGV | xform -e > $radscn"; + system "xform -e @ARGV > $radscn"; die "Could not load Radiance input\n" if ( $? ); system "rad2mgf $radscn > $mgfscn" if ( $geout ); } @@ -81,12 +100,244 @@ close RADSCN; # Generate octree system "oconv -w $radscn > $octree"; die "Could not compile scene\n" if ( $? ); -# Set up sampling of interior portal +# Output XML prologue +print +' + +'; +print "\n"; +print +'System + + + + Name + Manufacturer +'; +printf "\t\t%.3f\n", $dim[5] - $dim[4]; +printf "\t\t%.3f\n", $dim[1] - $dim[0]; +printf "\t\t%.3f\n", $dim[3] - $dim[2]; +print "\t\tIntegral\n"; +# Output MGF description if requested +if ( $geout ) { + print "\t\t\n"; + printf "xf -t %.6f %.6f 0\n", -($dim[0]+$dim[1])/2, -($dim[2]+$dim[3])/2; + open(MGFSCN, "< $mgfscn"); + while () { print $_; } + close MGFSCN; + print "xf\n"; + print "\t\t\n"; +} +print " \n"; +# Set up surface sampling +my $nx = int(sqrt($nsamp*($dim[1]-$dim[0])/($dim[3]-$dim[2])) + .5); +my $ny = int($nsamp/$nx + .5); +$nsamp = $nx * $ny; +my $ns = 2**$ttlog2; +my (@pdiv, $disk2sq, $sq2disk, $tcal, $kcal); +# Create data segments (all the work happens here) +if ( $tensortree ) { + do_tree_bsdf(); +} else { + do_matrix_bsdf(); +} +# Output XML epilogue +print +' + + +'; +# Clean up temporary files and exit +exec("rm -rf $td"); + +#-------------- End of main program segment --------------# + +#++++++++++++++ Tensor tree BSDF generation ++++++++++++++# +sub do_tree_bsdf { +# Shirley-Chiu mapping from unit square to disk +$sq2disk = ' +in_square_a = 2*in_square_x - 1; +in_square_b = 2*in_square_y - 1; +in_square_rgn = if(in_square_a + in_square_b, + if(in_square_a - in_square_b, 1, 2), + if(in_square_b - in_square_a, 3, 4)); +out_disk_r = .999995*select(in_square_rgn, in_square_a, in_square_b, + -in_square_a, -in_square_b); +out_disk_phi = PI/4 * select(in_square_rgn, + in_square_b/in_square_a, + 2 - in_square_a/in_square_b, + 4 + in_square_b/in_square_a, + if(in_square_b*in_square_b, + 6 - in_square_a/in_square_b, 0)); +Dx = out_disk_r*cos(out_disk_phi); +Dy = out_disk_r*sin(out_disk_phi); +Dz = sqrt(1 - out_disk_r*out_disk_r); +'; +# Shirley-Chiu mapping from unit disk to square +$disk2sq = ' +norm_radians(p) : if(-p - PI/4, p + 2*PI, p); +in_disk_r = .999995*sqrt(Dx*Dx + Dy*Dy); +in_disk_phi = norm_radians(atan2(Dy, Dx)); +in_disk_rgn = floor((in_disk_phi + PI/4)/(PI/2)) + 1; +out_square_a = select(in_disk_rgn, + in_disk_r, + (PI/2 - in_disk_phi)*in_disk_r/(PI/4), + -in_disk_r, + (in_disk_phi - 3*PI/2)*in_disk_r/(PI/4)); +out_square_b = select(in_disk_rgn, + in_disk_phi*in_disk_r/(PI/4), + in_disk_r, + (PI - in_disk_phi)*in_disk_r/(PI/4), + -in_disk_r); +out_square_x = (out_square_a + 1)/2; +out_square_y = (out_square_b + 1)/2; +'; +# Announce ourselves in XML output +print "\t\n"; +print "\t\tTensorTree$tensortree\n"; +print "\t\n"; + +# Start rcontrib processes for compute each side +do_tree_rtcontrib(0) if ( $doback ); +do_tree_rtcontrib(1) if ( $doforw ); + +} # end of sub do_tree_bsdf() + +# Run rcontrib process to generate tensor tree samples +sub do_tree_rtcontrib { + my $forw = shift; + my $matargs = "-m $bmodnm"; + if ( !$forw || !$doback || $tensortree==3 ) { $matargs .= " -m $fmodnm"; } + my $cmd = "rcontrib $rtargs -h -ff -fo -n $nproc -c $nsamp " . + "-e '$disk2sq' -bn '$ns*$ns' " . + "-b '$ns*floor(out_square_x*$ns)+floor(out_square_y*$ns)' " . + "-o $td/%s.flt $matargs $octree"; + if ( $tensortree == 3 ) { + # Isotropic BSDF + my $ns2 = $ns / 2; + $cmd = "cnt $ns2 $ny $nx " . + "| rcalc -e 'r1=rand(.8681*recno-.673892)' " . + "-e 'r2=rand(-5.37138*recno+67.1737811)' " . + "-e 'r3=rand(+3.17603772*recno+83.766771)' " . + "-e 'Dx=1-2*(\$1+r1)/$ns;Dy:0;Dz=sqrt(1-Dx*Dx)' " . + "-e 'xp=(\$3+r2)*(($dim[1]-$dim[0])/$nx)+$dim[0]' " . + "-e 'yp=(\$2+r3)*(($dim[3]-$dim[2])/$ny)+$dim[2]' " . + "-e 'zp=$dim[5-$forw]' -e 'myDz=Dz*($forw*2-1)' " . + "-e '\$1=xp-Dx;\$2=yp-Dy;\$3=zp-myDz' " . + "-e '\$4=Dx;\$5=Dy;\$6=myDz' -of " . + "| $cmd"; + } else { + # Anisotropic BSDF + # Sample area vertically to improve load balance, since + # shading systems usually have bilateral symmetry (L-R) + $cmd = "cnt $ns $ns $ny $nx " . + "| rcalc -e 'r1=rand(.8681*recno-.673892)' " . + "-e 'r2=rand(-5.37138*recno+67.1737811)' " . + "-e 'r3=rand(3.17603772*recno+83.766771)' " . + "-e 'r4=rand(-2.3857833*recno-964.72738)' " . + "-e 'in_square_x=(\$1+r1)/$ns' " . + "-e 'in_square_y=(\$2+r2)/$ns' -e '$sq2disk' " . + "-e 'xp=(\$4+r3)*(($dim[1]-$dim[0])/$nx)+$dim[0]' " . + "-e 'yp=(\$3+r4)*(($dim[3]-$dim[2])/$ny)+$dim[2]' " . + "-e 'zp=$dim[5-$forw]' -e 'myDz=Dz*($forw*2-1)' " . + "-e '\$1=xp-Dx;\$2=yp-Dy;\$3=zp-myDz' " . + "-e '\$4=Dx;\$5=Dy;\$6=myDz' -of " . + "| $cmd"; + } +# print STDERR "Starting: $cmd\n"; + system "$cmd" || die "Failure running rcontrib"; + ttree_out($forw); +} # end of do_tree_rtcontrib() + +# Simplify and output tensor tree results +sub ttree_out { + my $forw = shift; + my $side = ("Back","Front")[$forw]; + my $cmd; +# Only output one transmitted anisotropic distribution, preferring backwards +if ( !$forw || !$doback || $tensortree==3 ) { +print +' + System + Visible + CIE Illuminant D65 1nm.ssp + ASTM E308 1931 Y.dsp + +'; +print "\t\t\tTransmission $side\n"; +print +' LBNL/Shirley-Chiu + BTDF + +'; +$cmd = "rcalc -if3 -e 'Omega:PI/($ns*$ns)' " . + q{-e '$1=(0.265*$1+0.670*$2+0.065*$3)/Omega' }; +if ($pctcull >= 0) { + $cmd .= "-of $td/" . ($bmodnm,$fmodnm)[$forw] . ".flt " . + "| rttree_reduce -a -h -ff -t $pctcull -r $tensortree -g $ttlog2"; + system "$cmd" || die "Failure running rttree_reduce"; +} else { + $cmd .= "$td/" . ($bmodnm,$fmodnm)[$forw] . ".flt"; + print "{\n"; + system "$cmd" || die "Failure running rcalc"; + for (my $i = ($tensortree==3)*$ns*$ns*$ns/2; $i-- > 0; ) { + print "0\n"; + } + print "}\n"; +} +print +' + + +'; +} +# Output reflection +print +' + System + Visible + CIE Illuminant D65 1nm.ssp + ASTM E308 1931 Y.dsp + +'; +print "\t\t\tReflection $side\n"; +print +' LBNL/Shirley-Chiu + BRDF + +'; +$cmd = "rcalc -if3 -e 'Omega:PI/($ns*$ns)' " . + q{-e '$1=(0.265*$1+0.670*$2+0.065*$3)/Omega' }; +if ($pctcull >= 0) { + $cmd .= "-of $td/" . ($fmodnm,$bmodnm)[$forw] . ".flt " . + "| rttree_reduce -a -h -ff -t $pctcull -r $tensortree -g $ttlog2"; + system "$cmd" || die "Failure running rttree_reduce"; +} else { + $cmd .= "$td/" . ($fmodnm,$bmodnm)[$forw] . ".flt"; + print "{\n"; + system "$cmd" || die "Failure running rcalc"; + for (my $i = ($tensortree==3)*$ns*$ns*$ns/2; $i-- > 0; ) { + print "0\n"; + } + print "}\n"; +} +print +' + + +'; +} # end of ttree_out() + +#------------- End of do_tree_bsdf() & subroutines -------------# + +#+++++++++++++++ Klems matrix BSDF generation +++++++++++++++# +sub do_matrix_bsdf { +# Set up sampling of portal # Kbin to produce incident direction in full Klems basis with (x1,x2) randoms -my $tcal = ' +$tcal = ' DEGREE : PI/180; sq(x) : x*x; -Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90); +Kpola(r) : select(r+1, 0, 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); Kmax : Kaccum(Knaz(0)); @@ -104,11 +355,12 @@ KprojOmega = PI * if(Kbin-.5, 1 - sq(cos(Kpola(1)*DEGREE))); '; # Compute Klems bin from exiting ray direction (forward or backward) -my $kcal = ' +$kcal = ' DEGREE : PI/180; -Acos(x) : 1/DEGREE * if(x-1, 0, if(-1-x, 0, acos(x))); +abs(x) : if(x, x, -x); +Acos(x) : if(x-1, 0, if(-1-x, PI, acos(x))) / DEGREE; posangle(a) : if(-a, a + 2*PI, a); -Atan2(y,x) : 1/DEGREE * posangle(atan2(y,x)); +Atan2(y,x) : posangle(atan2(y,x)) / DEGREE; kpola(r) : select(r, 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); @@ -126,37 +378,36 @@ kbin2(pol,azi) = select(kfindrow(1, pol), kaccum(7) + kazn(azi,360/knaz(8)), kaccum(8) + kazn(azi,360/knaz(9)) ); -kbin = if(Dz, kbin2(Acos(Dz),Atan2(Dy,Dx)), 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 +# Compute scattering data using rcontrib my @tfarr; my @rfarr; my @tbarr; my @rbarr; my $cmd; -my $rtcmd = "rtcontrib -h -ff -fo -n $nproc -c $nsamp " . +my $rtcmd = "rcontrib $rtargs -h -ff -fo -n $nproc -c $nsamp " . "-e '$kcal' -b kbin -bn $ndiv " . - "-o '$td/%s.flt' -m $fmodnm -m $bmodnm $rtargs $octree"; + "-o '$td/%s.flt' -m $fmodnm -m $bmodnm $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'}; + 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(.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 '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(1.21*recno+2.75);x2=rand(-3.55*recno-7.57)' } . + 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 ( $? ); +chomp(@tfarr); @rfarr = `$rccmd $td/$bmodnm.flt`; die "Failure running: $rccmd $td/$bmodnm.flt\n" if ( $? ); +chomp(@rfarr); } if ( $doback ) { $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " . @@ -164,119 +415,100 @@ $cmd = "cnt $ndiv $ny $nx | rcalc -of -e '$tcal' " . "-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' } . + 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 ( $? ); +chomp(@tbarr); @rbarr = `$rccmd $td/$fmodnm.flt`; die "Failure running: $rccmd $td/$fmodnm.flt\n" if ( $? ); +chomp(@rbarr); } -# Output XML prologue +# Output angle basis print -' - - System - - - - Name - Manufacturer -'; -printf "\t\t\t%.3f\n", $dim[5] - $dim[4]; -printf "\t\t\t%.3f\n", $dim[1] - $dim[0]; -printf "\t\t\t%.3f\n", $dim[3] - $dim[2]; -print "\t\t\tIntegral\n"; -# Output MGF description if requested -if ( $geout ) { - print "\t\t\t\n"; - printf "xf -t %.6f %.6f 0\n", -($dim[0]+$dim[1])/2, -($dim[2]+$dim[3])/2; - system "cat $mgfscn"; - print "xf\n"; - print "\t\t\t\n"; -} -print ' - - Columns - +' + Columns + LBNL/Klems Full - - 0 - 1 - + + 0 + 1 + 0 5 - - - - 10 - 8 - - 5 - 15 - - - - 20 - 16 - - 15 - 25 - - - - 30 - 20 - - 25 - 35 - - - - 40 - 24 - - 35 - 45 - - - - 50 - 24 - - 45 - 55 - - - - 60 - 24 - - 55 - 65 - - - - 70 - 16 - - 65 - 75 - - - - 82.5 - 12 - - 75 - 90 - + + + 10 + 8 + + 5 + 15 + + + + 20 + 16 + + 15 + 25 + + + + 30 + 20 + + 25 + 35 + + + + 40 + 24 + + 35 + 45 + + + + 50 + 24 + + 45 + 55 + + + + 60 + 24 + + 55 + 65 + + + + 70 + 16 + + 65 + 75 + + + + 82.5 + 12 + + 75 + 90 + + '; if ( $doforw ) { -print ' +print +' System Visible CIE Illuminant D65 1nm.ssp @@ -291,13 +523,13 @@ print ' # 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 $tfarr[$ndiv*$id + $od], ",\n"; } print "\n"; } print -' - +' + System @@ -314,18 +546,19 @@ print # 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 $rfarr[$ndiv*$id + $od], ",\n"; } print "\n"; } print -' - +' + '; } if ( $doback ) { -print ' +print +' System Visible CIE Illuminant D65 1nm.ssp @@ -340,13 +573,13 @@ print ' # 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 $tbarr[$ndiv*$id + $od], ",\n"; } print "\n"; } print -' - +' + System @@ -363,20 +596,15 @@ print # 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 $rbarr[$ndiv*$id + $od], ",\n"; } print "\n"; } print -' - +' + '; } -# Output XML epilogue -print ' - - -'; -# Clean up temporary files -system "rm -rf $td"; +} +#------------- End of do_matrix_bsdf() --------------#