src/util/genklemsamp.pl

#!/usr/bin/perl -w
# RCSid $Id: genklemsamp.pl,v 2.9 2011/04/11 03:47:46 greg Exp $
#
# Sample Klems (full) directions impinging on surface(s)
#
#       G. Ward
#
use strict;
use File::Temp qw/ :mktemp  /;
if ($#ARGV < 0) {
        print STDERR "Usage: genklemsamp [-c N ][-f{a|f|d}] [view opts] [geom.rad ..]\n";
        exit 1;
}
my $td = mkdtemp("/tmp/genklemsamp.XXXXXX");
chomp $td;
my $nsamp = 1000;
my $fmt = "a";
my @vopts="-vo 1e-5";
while ($#ARGV >= 0) {
        if ("$ARGV[0]" eq "-c") {
                shift @ARGV;
                $nsamp = $ARGV[0];
        } elsif ("$ARGV[0]" =~ /^-f[afd]$/) {
                $fmt = substr($ARGV[0], 2, 1);
        } elsif ("$ARGV[0]" =~ /^-v[pdu]$/) {
                push @vopts, "@ARGV[0..3]";
                shift @ARGV; shift @ARGV; shift @ARGV;
        } elsif ("$ARGV[0]" =~ /^-v[fhvo]$/) {
                push @vopts, "@ARGV[0..1]";
                shift @ARGV;
        } elsif ("$ARGV[0]" =~ /^-v./) {
                print STDERR "Unsupported view option: $ARGV[0]\n";
                exit 1;
        } elsif ("$ARGV[0]" =~ /^-./) {
                print STDERR "Unknown option: $ARGV[0]\n";
                exit 1;
        } else {
                last;
        }
        shift @ARGV;
}
my $ofmt = "";
$ofmt = "-o$fmt" if ("$fmt" ne "a");
# Require parallel view and compile settings
push @vopts, "-vtl";
my $vwset = `vwright @vopts V`;
$_ = $vwset; s/^.*Vhn://; s/;.*$//;
my $width = $_;
$_ = $vwset; s/^.*Vvn://; s/;.*$//;
my $height = $_;
my $sca = sqrt($nsamp/($width*$height));
# Kbin is input to get direction in full Klems basis with (x1,x2) randoms
my $tcal = '
DEGREE : PI/180;
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);
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) - 1;
Kazi = 360*DEGREE * (Kcol + .5 - x2) / Knaz(Krow);
Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1));
sin_kpol = sin(Kpol);
K0 = cos(Kazi)*sin_kpol;
K1 = sin(Kazi)*sin_kpol;
K2 = sqrt(1 - sin_kpol*sin_kpol);
';
my $ndiv = 145;
# Do we have any Radiance input files?
if ($#ARGV >= 0) {
        system "oconv -f @ARGV > $td/surf.oct";
        # Set our own view center and size based on bounding cube
        $_ = $vwset; s/^.*Vdx://; s/;.*$//;
        my @vd = $_;
        $_ = $vwset; s/^.*Vdy://; s/;.*$//;
        push @vd, $_;
        $_ = $vwset; s/^.*Vdz://; s/;.*$//;
        push @vd, $_;
        my @bcube = split ' ', `getinfo -d < $td/surf.oct`;
        $width = $bcube[3]*sqrt(3);
        $height = $width;
        push @vopts, ("-vp", $bcube[0]+$bcube[3]/2-$width/2*$vd[0],
                        $bcube[1]+$bcube[3]/2-$width/2*$vd[1],
                        $bcube[2]+$bcube[3]/2-$width/2*$vd[2]);
        push @vopts, ("-vh", $width, "-vv", $height);
        $vwset = `vwright @vopts V`;
        $sca = sqrt($nsamp/($width*$height));
        my $xres;
        my $yres;
        my $ntot = 0;
        # This generally passes through the loop twice to get density right
        while ($ntot < $nsamp) {
                $xres = int($width*$sca) + 1;
                $yres = int($height*$sca) + 1;
                system "vwrays -ff -x $xres -y $yres -pa 0 -pj .7 @vopts " .
                        "| rtrace -w -h -ff -opLN $td/surf.oct " .
                        q{| rcalc -e 'cond=and(5e9-$4,nOK);$1=$1;$2=$2;$3=$3' } .
                                "-e 'and(a,b):if(a,b,a);sq(x):x*x' -e '$vwset' " .
                                q{-e 'nOK=sq(Vdx*$5+Vdy*$6+Vdz*$7)-.999' } .
                                "-if7 -of > $td/origins.flt";
                $ntot = -s "$td/origins.flt";
                $ntot /= 3*4;           # number of bytes per sample position
                if ($ntot == 0) {
                        if ($nsamp < 200) {
                                $sca = sqrt(200/($width*$height));
                                redo;
                        }
                        print STDERR "View direction does not correspond to any surfaces\n";
                        exit 1;
                }
                $sca *= 1.05 * sqrt($nsamp/$ntot);
        }
        # All set to produce our samples
        for (my $k = 1; $k <= $ndiv; $k++) {
                my $rn = rand(10);
                my $r1 = rand; my $r2 = rand;
                # Chance of using = (number_still_needed)/(number_left_avail)
                system "rcalc $ofmt -e '$tcal' -e '$vwset' " .
                        "-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
                        "-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
                        q{-e '$1=$1+Vo*Vdx; $2=$2+Vo*Vdy; $3=$3+Vo*Vdz' } .
                        q{-e '$4=-K0*Vhx-K1*Vvx+K2*Vdx' } .
                        q{-e '$5=-K0*Vhy-K1*Vvy+K2*Vdy' } .
                        q{-e '$6=-K0*Vhz-K1*Vvz+K2*Vdz' } .
                        "-if3 $td/origins.flt";
        }
} else {
        # No Radiance input files, so just sample over parallel view
        my $xres = int($width*$sca) + 1;
        my $yres = int($height*$sca) + 1;
        my $ntot = $xres * $yres;
        for (my $k = 1; $k <= $ndiv; $k++) {
                my $rn = rand(10);
                my $r1 = rand; my $r2 = rand;
                my $r3 = rand; my $r4 = rand;
                # Chance of using = (number_still_needed)/(number_left_avail)
                system "cnt $yres $xres | rcalc $ofmt -e '$tcal' -e '$vwset' " .
                        q{-e 'xc=$2;yc=$1' } .
                        "-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
                        "-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
                        "-e 'hpos=Vhn*((xc+rand($r3+recno))/$xres - .5)' " .
                        "-e 'vpos=Vvn*((yc+rand($r4+recno))/$yres - .5)' " .
                        q{-e '$1=Vpx+Vo*Vdx+hpos*Vhx+vpos*Vvx' } .
                        q{-e '$2=Vpy+Vo*Vdy+hpos*Vhy+vpos*Vvy' } .
                        q{-e '$3=Vpz+Vo*Vdz+hpos*Vhz+vpos*Vvz' } .
                        q{-e '$4=-K0*Vhx-K1*Vvx+K2*Vdx' } .
                        q{-e '$5=-K0*Vhy-K1*Vvy+K2*Vdy' } .
                        q{-e '$6=-K0*Vhz-K1*Vvz+K2*Vdz' } ;
        }
}
system "rm -rf $td";
Revision:	2.10
Committed:	Sat Apr 16 00:39:07 2011 UTC (14 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.9:	+3 -2 lines
Log Message:	Bug fix for 32-bit OS
#	User	Rev	Content
1	greg	2.1	#!/usr/bin/perl -w
2	greg	2.10	# RCSid $Id: genklemsamp.pl,v 2.9 2011/04/11 03:47:46 greg Exp $
3	greg	2.1	#
4			# Sample Klems (full) directions impinging on surface(s)
5			#
6	greg	2.2	# G. Ward
7			#
8	greg	2.3	use strict;
9	greg	2.10	use File::Temp qw/ :mktemp /;
10	greg	2.1	if ($#ARGV < 0) {
11	greg	2.4	print STDERR "Usage: genklemsamp [-c N ][-f{a\|f\|d}] [view opts] [geom.rad ..]\n";
12	greg	2.1	exit 1;
13			}
14	greg	2.10	my $td = mkdtemp("/tmp/genklemsamp.XXXXXX");
15	greg	2.1	chomp $td;
16			my $nsamp = 1000;
17			my $fmt = "a";
18			my @vopts="-vo 1e-5";
19			while ($#ARGV >= 0) {
20			if ("$ARGV[0]" eq "-c") {
21			shift @ARGV;
22			$nsamp = $ARGV[0];
23			} elsif ("$ARGV[0]" =~ /^-f[afd]$/) {
24			$fmt = substr($ARGV[0], 2, 1);
25			} elsif ("$ARGV[0]" =~ /^-v[pdu]$/) {
26			push @vopts, "@ARGV[0..3]";
27			shift @ARGV; shift @ARGV; shift @ARGV;
28			} elsif ("$ARGV[0]" =~ /^-v[fhvo]$/) {
29			push @vopts, "@ARGV[0..1]";
30			shift @ARGV;
31			} elsif ("$ARGV[0]" =~ /^-v./) {
32	greg	2.4	print STDERR "Unsupported view option: $ARGV[0]\n";
33	greg	2.1	exit 1;
34			} elsif ("$ARGV[0]" =~ /^-./) {
35	greg	2.4	print STDERR "Unknown option: $ARGV[0]\n";
36	greg	2.1	exit 1;
37			} else {
38			last;
39			}
40			shift @ARGV;
41			}
42			my $ofmt = "";
43			$ofmt = "-o$fmt" if ("$fmt" ne "a");
44			# Require parallel view and compile settings
45			push @vopts, "-vtl";
46			my $vwset = `vwright @vopts V`;
47			$_ = $vwset; s/^.Vhn://; s/;.$//;
48			my $width = $_;
49			$_ = $vwset; s/^.Vvn://; s/;.$//;
50			my $height = $_;
51			my $sca = sqrt($nsamp/($width*$height));
52			# Kbin is input to get direction in full Klems basis with (x1,x2) randoms
53			my $tcal = '
54			DEGREE : PI/180;
55			Kpola(r) : select(r+1, -5, 5, 15, 25, 35, 45, 55, 65, 75, 90);
56			Knaz(r) : select(r, 1, 8, 16, 20, 24, 24, 24, 16, 12);
57			Kaccum(r) : if(r-.5, Knaz(r) + Kaccum(r-1), 0);
58			Kmax : Kaccum(Knaz(0));
59			Kfindrow(r, rem) : if(rem-Knaz(r)-.5, Kfindrow(r+1, rem-Knaz(r)), r);
60			Krow = if(Kbin-(Kmax+.5), 0, Kfindrow(1, Kbin));
61			Kcol = Kbin - Kaccum(Krow-1) - 1;
62			Kazi = 360DEGREE (Kcol + .5 - x2) / Knaz(Krow);
63			Kpol = DEGREE * (x1Kpola(Krow) + (1-x1)Kpola(Krow-1));
64			sin_kpol = sin(Kpol);
65			K0 = cos(Kazi)*sin_kpol;
66			K1 = sin(Kazi)*sin_kpol;
67	greg	2.5	K2 = sqrt(1 - sin_kpol*sin_kpol);
68	greg	2.1	';
69			my $ndiv = 145;
70			# Do we have any Radiance input files?
71			if ($#ARGV >= 0) {
72			system "oconv -f @ARGV > $td/surf.oct";
73			# Set our own view center and size based on bounding cube
74			$_ = $vwset; s/^.Vdx://; s/;.$//;
75			my @vd = $_;
76			$_ = $vwset; s/^.Vdy://; s/;.$//;
77			push @vd, $_;
78			$_ = $vwset; s/^.Vdz://; s/;.$//;
79			push @vd, $_;
80	greg	2.8	my @bcube = split ' ', `getinfo -d < $td/surf.oct`;
81	greg	2.1	$width = $bcube[3]*sqrt(3);
82			$height = $width;
83			push @vopts, ("-vp", $bcube[0]+$bcube[3]/2-$width/2*$vd[0],
84			$bcube[1]+$bcube[3]/2-$width/2*$vd[1],
85			$bcube[2]+$bcube[3]/2-$width/2*$vd[2]);
86	greg	2.7	push @vopts, ("-vh", $width, "-vv", $height);
87	greg	2.1	$vwset = `vwright @vopts V`;
88	greg	2.6	$sca = sqrt($nsamp/($width*$height));
89	greg	2.1	my $xres;
90			my $yres;
91	greg	2.4	my $ntot = 0;
92	greg	2.1	# This generally passes through the loop twice to get density right
93	greg	2.4	while ($ntot < $nsamp) {
94	greg	2.1	$xres = int($width*$sca) + 1;
95			$yres = int($height*$sca) + 1;
96			system "vwrays -ff -x $xres -y $yres -pa 0 -pj .7 @vopts " .
97			"\| rtrace -w -h -ff -opLN $td/surf.oct " .
98			q{\| rcalc -e 'cond=and(5e9-$4,nOK);$1=$1;$2=$2;$3=$3' } .
99			"-e 'and(a,b):if(a,b,a);sq(x):x*x' -e '$vwset' " .
100			q{-e 'nOK=sq(Vdx$5+Vdy$6+Vdz*$7)-.999' } .
101			"-if7 -of > $td/origins.flt";
102			$ntot = -s "$td/origins.flt";
103	greg	2.7	$ntot /= 3*4; # number of bytes per sample position
104	greg	2.1	if ($ntot == 0) {
105	greg	2.6	if ($nsamp < 200) {
106	greg	2.4	$sca = sqrt(200/($width*$height));
107			redo;
108			}
109			print STDERR "View direction does not correspond to any surfaces\n";
110	greg	2.1	exit 1;
111			}
112			$sca = 1.05 sqrt($nsamp/$ntot);
113	greg	2.4	}
114	greg	2.1	# All set to produce our samples
115	greg	2.4	for (my $k = 1; $k <= $ndiv; $k++) {
116	greg	2.1	my $rn = rand(10);
117			my $r1 = rand; my $r2 = rand;
118			# Chance of using = (number_still_needed)/(number_left_avail)
119			system "rcalc $ofmt -e '$tcal' -e '$vwset' " .
120			"-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
121			"-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
122			q{-e '$1=$1+VoVdx; $2=$2+VoVdy; $3=$3+Vo*Vdz' } .
123	greg	2.9	q{-e '$4=-K0Vhx-K1Vvx+K2*Vdx' } .
124			q{-e '$5=-K0Vhy-K1Vvy+K2*Vdy' } .
125			q{-e '$6=-K0Vhz-K1Vvz+K2*Vdz' } .
126	greg	2.1	"-if3 $td/origins.flt";
127			}
128			} else {
129			# No Radiance input files, so just sample over parallel view
130			my $xres = int($width*$sca) + 1;
131			my $yres = int($height*$sca) + 1;
132			my $ntot = $xres * $yres;
133	greg	2.4	for (my $k = 1; $k <= $ndiv; $k++) {
134	greg	2.1	my $rn = rand(10);
135			my $r1 = rand; my $r2 = rand;
136			my $r3 = rand; my $r4 = rand;
137			# Chance of using = (number_still_needed)/(number_left_avail)
138			system "cnt $yres $xres \| rcalc $ofmt -e '$tcal' -e '$vwset' " .
139			q{-e 'xc=$2;yc=$1' } .
140			"-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
141			"-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
142			"-e 'hpos=Vhn*((xc+rand($r3+recno))/$xres - .5)' " .
143			"-e 'vpos=Vvn*((yc+rand($r4+recno))/$yres - .5)' " .
144			q{-e '$1=Vpx+VoVdx+hposVhx+vpos*Vvx' } .
145			q{-e '$2=Vpy+VoVdy+hposVhy+vpos*Vvy' } .
146			q{-e '$3=Vpz+VoVdz+hposVhz+vpos*Vvz' } .
147	greg	2.9	q{-e '$4=-K0Vhx-K1Vvx+K2*Vdx' } .
148			q{-e '$5=-K0Vhy-K1Vvy+K2*Vdy' } .
149			q{-e '$6=-K0Vhz-K1Vvz+K2*Vdz' } ;
150	greg	2.1	}
151			}
152			system "rm -rf $td";