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