Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68371474
ambcomp.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jun 27, 02:44
Size
20 KB
Mime Type
text/x-c
Expires
Sat, Jun 29, 02:44 (2 d)
Engine
blob
Format
Raw Data
Handle
18549866
Attached To
R10977 RADIANCE Photon Map
ambcomp.c
View Options
#ifndef lint
static
const
char
RCSid
[]
=
"$Id: ambcomp.c,v 2.86 2021/02/17 01:29:22 greg Exp $"
;
#endif
/*
* Routines to compute "ambient" values using Monte Carlo
*
* Hessian calculations based on "Practical Hessian-Based Error Control
* for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz
* from ACM SIGGRAPH Asia 2012 conference proceedings.
*
* Added book-keeping optimization to avoid calculations that would
* cancel due to traversal both directions on edges that are adjacent
* to same-valued triangles. This cuts about half of Hessian math.
*
* Declarations of external symbols in ambient.h
*/
#include "copyright.h"
#include "ray.h"
#include "ambient.h"
#include "random.h"
#ifndef MINADIV
#define MINADIV 7
/* minimum # divisions in each dimension */
#endif
extern
void
SDsquare2disk
(
double
ds
[
2
],
double
seedx
,
double
seedy
);
typedef
struct
{
COLOR
v
;
/* hemisphere sample value */
float
d
;
/* reciprocal distance */
FVECT
p
;
/* intersection point */
}
AMBSAMP
;
/* sample value */
typedef
struct
{
RAY
*
rp
;
/* originating ray sample */
int
ns
;
/* number of samples per axis */
int
sampOK
;
/* acquired full sample set? */
COLOR
acoef
;
/* division contribution coefficient */
double
acol
[
3
];
/* accumulated color */
FVECT
ux
,
uy
;
/* tangent axis unit vectors */
AMBSAMP
sa
[
1
];
/* sample array (extends struct) */
}
AMBHEMI
;
/* ambient sample hemisphere */
#define AI(h,i,j) ((i)*(h)->ns + (j))
#define ambsam(h,i,j) (h)->sa[AI(h,i,j)]
typedef
struct
{
FVECT
r_i
,
r_i1
,
e_i
,
rcp
,
rI2_eJ2
;
double
I1
,
I2
;
}
FFTRI
;
/* vectors and coefficients for Hessian calculation */
static
int
ambcollision
(
/* proposed direciton collides? */
AMBHEMI
*
hp
,
int
i
,
int
j
,
FVECT
dv
)
{
double
cos_thresh
;
int
ii
,
jj
;
/* min. spacing = 1/4th division */
cos_thresh
=
(
PI
/
4.
)
/
(
double
)
hp
->
ns
;
cos_thresh
=
1.
-
.5
*
cos_thresh
*
cos_thresh
;
/* check existing neighbors */
for
(
ii
=
i
-
1
;
ii
<=
i
+
1
;
ii
++
)
{
if
(
ii
<
0
)
continue
;
if
(
ii
>=
hp
->
ns
)
break
;
for
(
jj
=
j
-
1
;
jj
<=
j
+
1
;
jj
++
)
{
AMBSAMP
*
ap
;
FVECT
avec
;
double
dprod
;
if
(
jj
<
0
)
continue
;
if
(
jj
>=
hp
->
ns
)
break
;
if
((
ii
==
i
)
&
(
jj
==
j
))
continue
;
ap
=
&
ambsam
(
hp
,
ii
,
jj
);
if
(
ap
->
d
<=
.5
/
FHUGE
)
continue
;
/* no one home */
VSUB
(
avec
,
ap
->
p
,
hp
->
rp
->
rop
);
dprod
=
DOT
(
avec
,
dv
);
if
(
dprod
>=
cos_thresh
*
VLEN
(
avec
))
return
(
1
);
/* collision */
}
}
return
(
0
);
/* nothing to worry about */
}
static
int
ambsample
(
/* initial ambient division sample */
AMBHEMI
*
hp
,
int
i
,
int
j
,
int
n
)
{
AMBSAMP
*
ap
=
&
ambsam
(
hp
,
i
,
j
);
RAY
ar
;
int
hlist
[
3
],
ii
;
double
spt
[
2
],
zd
;
/* generate hemispherical sample */
/* ambient coefficient for weight */
if
(
ambacc
>
FTINY
)
setcolor
(
ar
.
rcoef
,
AVGREFL
,
AVGREFL
,
AVGREFL
);
else
copycolor
(
ar
.
rcoef
,
hp
->
acoef
);
if
(
rayorigin
(
&
ar
,
AMBIENT
,
hp
->
rp
,
ar
.
rcoef
)
<
0
)
return
(
0
);
if
(
ambacc
>
FTINY
)
{
multcolor
(
ar
.
rcoef
,
hp
->
acoef
);
scalecolor
(
ar
.
rcoef
,
1.
/
AVGREFL
);
}
hlist
[
0
]
=
hp
->
rp
->
rno
;
hlist
[
1
]
=
j
;
hlist
[
2
]
=
i
;
multisamp
(
spt
,
2
,
urand
(
ilhash
(
hlist
,
3
)
+
n
));
resample:
SDsquare2disk
(
spt
,
(
j
+
spt
[
1
])
/
hp
->
ns
,
(
i
+
spt
[
0
])
/
hp
->
ns
);
zd
=
sqrt
(
1.
-
spt
[
0
]
*
spt
[
0
]
-
spt
[
1
]
*
spt
[
1
]);
for
(
ii
=
3
;
ii
--
;
)
ar
.
rdir
[
ii
]
=
spt
[
0
]
*
hp
->
ux
[
ii
]
+
spt
[
1
]
*
hp
->
uy
[
ii
]
+
zd
*
hp
->
rp
->
ron
[
ii
];
checknorm
(
ar
.
rdir
);
/* avoid coincident samples */
if
(
!
n
&&
ambcollision
(
hp
,
i
,
j
,
ar
.
rdir
))
{
spt
[
0
]
=
frandom
();
spt
[
1
]
=
frandom
();
goto
resample
;
/* reject this sample */
}
dimlist
[
ndims
++
]
=
AI
(
hp
,
i
,
j
)
+
90171
;
rayvalue
(
&
ar
);
/* evaluate ray */
ndims
--
;
zd
=
raydistance
(
&
ar
);
if
(
zd
<=
FTINY
)
return
(
0
);
/* should never happen */
multcolor
(
ar
.
rcol
,
ar
.
rcoef
);
/* apply coefficient */
if
(
zd
*
ap
->
d
<
1.0
)
/* new/closer distance? */
ap
->
d
=
1.0
/
zd
;
if
(
!
n
)
{
/* record first vertex & value */
if
(
zd
>
10.0
*
thescene
.
cusize
+
1000.
)
zd
=
10.0
*
thescene
.
cusize
+
1000.
;
VSUM
(
ap
->
p
,
ar
.
rorg
,
ar
.
rdir
,
zd
);
copycolor
(
ap
->
v
,
ar
.
rcol
);
}
else
{
/* else update recorded value */
hp
->
acol
[
RED
]
-=
colval
(
ap
->
v
,
RED
);
hp
->
acol
[
GRN
]
-=
colval
(
ap
->
v
,
GRN
);
hp
->
acol
[
BLU
]
-=
colval
(
ap
->
v
,
BLU
);
zd
=
1.0
/
(
double
)(
n
+
1
);
scalecolor
(
ar
.
rcol
,
zd
);
zd
*=
(
double
)
n
;
scalecolor
(
ap
->
v
,
zd
);
addcolor
(
ap
->
v
,
ar
.
rcol
);
}
addcolor
(
hp
->
acol
,
ap
->
v
);
/* add to our sum */
return
(
1
);
}
/* Estimate variance based on ambient division differences */
static
float
*
getambdiffs
(
AMBHEMI
*
hp
)
{
const
double
normf
=
1.
/
bright
(
hp
->
acoef
);
float
*
earr
=
(
float
*
)
calloc
(
hp
->
ns
*
hp
->
ns
,
sizeof
(
float
));
float
*
ep
;
AMBSAMP
*
ap
;
double
b
,
b1
,
d2
;
int
i
,
j
;
if
(
earr
==
NULL
)
/* out of memory? */
return
(
NULL
);
/* sum squared neighbor diffs */
for
(
ap
=
hp
->
sa
,
ep
=
earr
,
i
=
0
;
i
<
hp
->
ns
;
i
++
)
for
(
j
=
0
;
j
<
hp
->
ns
;
j
++
,
ap
++
,
ep
++
)
{
b
=
bright
(
ap
[
0
].
v
);
if
(
i
)
{
/* from above */
b1
=
bright
(
ap
[
-
hp
->
ns
].
v
);
d2
=
b
-
b1
;
d2
*=
d2
*
normf
/
(
b
+
b1
);
ep
[
0
]
+=
d2
;
ep
[
-
hp
->
ns
]
+=
d2
;
}
if
(
!
j
)
continue
;
/* from behind */
b1
=
bright
(
ap
[
-
1
].
v
);
d2
=
b
-
b1
;
d2
*=
d2
*
normf
/
(
b
+
b1
);
ep
[
0
]
+=
d2
;
ep
[
-
1
]
+=
d2
;
if
(
!
i
)
continue
;
/* diagonal */
b1
=
bright
(
ap
[
-
hp
->
ns
-
1
].
v
);
d2
=
b
-
b1
;
d2
*=
d2
*
normf
/
(
b
+
b1
);
ep
[
0
]
+=
d2
;
ep
[
-
hp
->
ns
-
1
]
+=
d2
;
}
/* correct for number of neighbors */
earr
[
0
]
*=
8.
/
3.
;
earr
[
hp
->
ns
-
1
]
*=
8.
/
3.
;
earr
[(
hp
->
ns
-
1
)
*
hp
->
ns
]
*=
8.
/
3.
;
earr
[(
hp
->
ns
-
1
)
*
hp
->
ns
+
hp
->
ns
-
1
]
*=
8.
/
3.
;
for
(
i
=
1
;
i
<
hp
->
ns
-
1
;
i
++
)
{
earr
[
i
*
hp
->
ns
]
*=
8.
/
5.
;
earr
[
i
*
hp
->
ns
+
hp
->
ns
-
1
]
*=
8.
/
5.
;
}
for
(
j
=
1
;
j
<
hp
->
ns
-
1
;
j
++
)
{
earr
[
j
]
*=
8.
/
5.
;
earr
[(
hp
->
ns
-
1
)
*
hp
->
ns
+
j
]
*=
8.
/
5.
;
}
return
(
earr
);
}
/* Perform super-sampling on hemisphere (introduces bias) */
static
void
ambsupersamp
(
AMBHEMI
*
hp
,
int
cnt
)
{
float
*
earr
=
getambdiffs
(
hp
);
double
e2rem
=
0
;
float
*
ep
;
int
i
,
j
,
n
,
nss
;
if
(
earr
==
NULL
)
/* just skip calc. if no memory */
return
;
/* accumulate estimated variances */
for
(
ep
=
earr
+
hp
->
ns
*
hp
->
ns
;
ep
>
earr
;
)
e2rem
+=
*--
ep
;
ep
=
earr
;
/* perform super-sampling */
for
(
i
=
0
;
i
<
hp
->
ns
;
i
++
)
for
(
j
=
0
;
j
<
hp
->
ns
;
j
++
)
{
if
(
e2rem
<=
FTINY
)
goto
done
;
/* nothing left to do */
nss
=
*
ep
/
e2rem
*
cnt
+
frandom
();
for
(
n
=
1
;
n
<=
nss
&&
ambsample
(
hp
,
i
,
j
,
n
);
n
++
)
if
(
!--
cnt
)
goto
done
;
e2rem
-=
*
ep
++
;
/* update remainder */
}
done:
free
(
earr
);
}
static
AMBHEMI
*
samp_hemi
(
/* sample indirect hemisphere */
COLOR
rcol
,
RAY
*
r
,
double
wt
)
{
AMBHEMI
*
hp
;
double
d
;
int
n
,
i
,
j
;
/* insignificance check */
if
(
bright
(
rcol
)
<=
FTINY
)
return
(
NULL
);
/* set number of divisions */
if
(
ambacc
<=
FTINY
&&
wt
>
(
d
=
0.8
*
intens
(
rcol
)
*
r
->
rweight
/
(
ambdiv
*
minweight
)))
wt
=
d
;
/* avoid ray termination */
n
=
sqrt
(
ambdiv
*
wt
)
+
0.5
;
i
=
1
+
(
MINADIV
-
1
)
*
(
ambacc
>
FTINY
);
if
(
n
<
i
)
/* use minimum number of samples? */
n
=
i
;
/* allocate sampling array */
hp
=
(
AMBHEMI
*
)
malloc
(
sizeof
(
AMBHEMI
)
+
sizeof
(
AMBSAMP
)
*
(
n
*
n
-
1
));
if
(
hp
==
NULL
)
error
(
SYSTEM
,
"out of memory in samp_hemi"
);
hp
->
rp
=
r
;
hp
->
ns
=
n
;
hp
->
acol
[
RED
]
=
hp
->
acol
[
GRN
]
=
hp
->
acol
[
BLU
]
=
0.0
;
memset
(
hp
->
sa
,
0
,
sizeof
(
AMBSAMP
)
*
n
*
n
);
hp
->
sampOK
=
0
;
/* assign coefficient */
copycolor
(
hp
->
acoef
,
rcol
);
d
=
1.0
/
(
n
*
n
);
scalecolor
(
hp
->
acoef
,
d
);
/* make tangent plane axes */
if
(
!
getperpendicular
(
hp
->
ux
,
r
->
ron
,
1
))
error
(
CONSISTENCY
,
"bad ray direction in samp_hemi"
);
VCROSS
(
hp
->
uy
,
r
->
ron
,
hp
->
ux
);
/* sample divisions */
for
(
i
=
hp
->
ns
;
i
--
;
)
for
(
j
=
hp
->
ns
;
j
--
;
)
hp
->
sampOK
+=
ambsample
(
hp
,
i
,
j
,
0
);
copycolor
(
rcol
,
hp
->
acol
);
if
(
!
hp
->
sampOK
)
{
/* utter failure? */
free
(
hp
);
return
(
NULL
);
}
if
(
hp
->
sampOK
<
hp
->
ns
*
hp
->
ns
)
{
hp
->
sampOK
*=
-
1
;
/* soft failure */
return
(
hp
);
}
if
(
hp
->
sampOK
<=
MINADIV
*
MINADIV
)
return
(
hp
);
/* don't bother super-sampling */
n
=
ambssamp
*
wt
+
0.5
;
if
(
n
>
8
)
{
/* perform super-sampling? */
ambsupersamp
(
hp
,
n
);
copycolor
(
rcol
,
hp
->
acol
);
}
return
(
hp
);
/* all is well */
}
/* Return brightness of farthest ambient sample */
static
double
back_ambval
(
AMBHEMI
*
hp
,
const
int
n1
,
const
int
n2
,
const
int
n3
)
{
if
(
hp
->
sa
[
n1
].
d
<=
hp
->
sa
[
n2
].
d
)
{
if
(
hp
->
sa
[
n1
].
d
<=
hp
->
sa
[
n3
].
d
)
return
(
colval
(
hp
->
sa
[
n1
].
v
,
CIEY
));
return
(
colval
(
hp
->
sa
[
n3
].
v
,
CIEY
));
}
if
(
hp
->
sa
[
n2
].
d
<=
hp
->
sa
[
n3
].
d
)
return
(
colval
(
hp
->
sa
[
n2
].
v
,
CIEY
));
return
(
colval
(
hp
->
sa
[
n3
].
v
,
CIEY
));
}
/* Compute vectors and coefficients for Hessian/gradient calcs */
static
void
comp_fftri
(
FFTRI
*
ftp
,
AMBHEMI
*
hp
,
const
int
n0
,
const
int
n1
)
{
double
rdot_cp
,
dot_e
,
dot_er
,
rdot_r
,
rdot_r1
,
J2
;
int
ii
;
VSUB
(
ftp
->
r_i
,
hp
->
sa
[
n0
].
p
,
hp
->
rp
->
rop
);
VSUB
(
ftp
->
r_i1
,
hp
->
sa
[
n1
].
p
,
hp
->
rp
->
rop
);
VSUB
(
ftp
->
e_i
,
hp
->
sa
[
n1
].
p
,
hp
->
sa
[
n0
].
p
);
VCROSS
(
ftp
->
rcp
,
ftp
->
r_i
,
ftp
->
r_i1
);
rdot_cp
=
1.0
/
DOT
(
ftp
->
rcp
,
ftp
->
rcp
);
dot_e
=
DOT
(
ftp
->
e_i
,
ftp
->
e_i
);
dot_er
=
DOT
(
ftp
->
e_i
,
ftp
->
r_i
);
rdot_r
=
1.0
/
DOT
(
ftp
->
r_i
,
ftp
->
r_i
);
rdot_r1
=
1.0
/
DOT
(
ftp
->
r_i1
,
ftp
->
r_i1
);
ftp
->
I1
=
acos
(
DOT
(
ftp
->
r_i
,
ftp
->
r_i1
)
*
sqrt
(
rdot_r
*
rdot_r1
)
)
*
sqrt
(
rdot_cp
);
ftp
->
I2
=
(
DOT
(
ftp
->
e_i
,
ftp
->
r_i1
)
*
rdot_r1
-
dot_er
*
rdot_r
+
dot_e
*
ftp
->
I1
)
*
0.5
*
rdot_cp
;
J2
=
(
0.5
*
(
rdot_r
-
rdot_r1
)
-
dot_er
*
ftp
->
I2
)
/
dot_e
;
for
(
ii
=
3
;
ii
--
;
)
ftp
->
rI2_eJ2
[
ii
]
=
ftp
->
I2
*
ftp
->
r_i
[
ii
]
+
J2
*
ftp
->
e_i
[
ii
];
}
/* Compose 3x3 matrix from two vectors */
static
void
compose_matrix
(
FVECT
mat
[
3
],
FVECT
va
,
FVECT
vb
)
{
mat
[
0
][
0
]
=
2.0
*
va
[
0
]
*
vb
[
0
];
mat
[
1
][
1
]
=
2.0
*
va
[
1
]
*
vb
[
1
];
mat
[
2
][
2
]
=
2.0
*
va
[
2
]
*
vb
[
2
];
mat
[
0
][
1
]
=
mat
[
1
][
0
]
=
va
[
0
]
*
vb
[
1
]
+
va
[
1
]
*
vb
[
0
];
mat
[
0
][
2
]
=
mat
[
2
][
0
]
=
va
[
0
]
*
vb
[
2
]
+
va
[
2
]
*
vb
[
0
];
mat
[
1
][
2
]
=
mat
[
2
][
1
]
=
va
[
1
]
*
vb
[
2
]
+
va
[
2
]
*
vb
[
1
];
}
/* Compute partial 3x3 Hessian matrix for edge */
static
void
comp_hessian
(
FVECT
hess
[
3
],
FFTRI
*
ftp
,
FVECT
nrm
)
{
FVECT
ncp
;
FVECT
m1
[
3
],
m2
[
3
],
m3
[
3
],
m4
[
3
];
double
d1
,
d2
,
d3
,
d4
;
double
I3
,
J3
,
K3
;
int
i
,
j
;
/* compute intermediate coefficients */
d1
=
1.0
/
DOT
(
ftp
->
r_i
,
ftp
->
r_i
);
d2
=
1.0
/
DOT
(
ftp
->
r_i1
,
ftp
->
r_i1
);
d3
=
1.0
/
DOT
(
ftp
->
e_i
,
ftp
->
e_i
);
d4
=
DOT
(
ftp
->
e_i
,
ftp
->
r_i
);
I3
=
(
DOT
(
ftp
->
e_i
,
ftp
->
r_i1
)
*
d2
*
d2
-
d4
*
d1
*
d1
+
3.0
/
d3
*
ftp
->
I2
)
/
(
4.0
*
DOT
(
ftp
->
rcp
,
ftp
->
rcp
)
);
J3
=
0.25
*
d3
*
(
d1
*
d1
-
d2
*
d2
)
-
d4
*
d3
*
I3
;
K3
=
d3
*
(
ftp
->
I2
-
I3
/
d1
-
2.0
*
d4
*
J3
);
/* intermediate matrices */
VCROSS
(
ncp
,
nrm
,
ftp
->
e_i
);
compose_matrix
(
m1
,
ncp
,
ftp
->
rI2_eJ2
);
compose_matrix
(
m2
,
ftp
->
r_i
,
ftp
->
r_i
);
compose_matrix
(
m3
,
ftp
->
e_i
,
ftp
->
e_i
);
compose_matrix
(
m4
,
ftp
->
r_i
,
ftp
->
e_i
);
d1
=
DOT
(
nrm
,
ftp
->
rcp
);
d2
=
-
d1
*
ftp
->
I2
;
d1
*=
2.0
;
for
(
i
=
3
;
i
--
;
)
/* final matrix sum */
for
(
j
=
3
;
j
--
;
)
{
hess
[
i
][
j
]
=
m1
[
i
][
j
]
+
d1
*
(
I3
*
m2
[
i
][
j
]
+
K3
*
m3
[
i
][
j
]
+
2.0
*
J3
*
m4
[
i
][
j
]
);
hess
[
i
][
j
]
+=
d2
*
(
i
==
j
);
hess
[
i
][
j
]
*=
-
1.0
/
PI
;
}
}
/* Reverse hessian calculation result for edge in other direction */
static
void
rev_hessian
(
FVECT
hess
[
3
])
{
int
i
;
for
(
i
=
3
;
i
--
;
)
{
hess
[
i
][
0
]
=
-
hess
[
i
][
0
];
hess
[
i
][
1
]
=
-
hess
[
i
][
1
];
hess
[
i
][
2
]
=
-
hess
[
i
][
2
];
}
}
/* Add to radiometric Hessian from the given triangle */
static
void
add2hessian
(
FVECT
hess
[
3
],
FVECT
ehess1
[
3
],
FVECT
ehess2
[
3
],
FVECT
ehess3
[
3
],
double
v
)
{
int
i
,
j
;
for
(
i
=
3
;
i
--
;
)
for
(
j
=
3
;
j
--
;
)
hess
[
i
][
j
]
+=
v
*
(
ehess1
[
i
][
j
]
+
ehess2
[
i
][
j
]
+
ehess3
[
i
][
j
]
);
}
/* Compute partial displacement form factor gradient for edge */
static
void
comp_gradient
(
FVECT
grad
,
FFTRI
*
ftp
,
FVECT
nrm
)
{
FVECT
ncp
;
double
f1
;
int
i
;
f1
=
2.0
*
DOT
(
nrm
,
ftp
->
rcp
);
VCROSS
(
ncp
,
nrm
,
ftp
->
e_i
);
for
(
i
=
3
;
i
--
;
)
grad
[
i
]
=
(
0.5
/
PI
)
*
(
ftp
->
I1
*
ncp
[
i
]
+
f1
*
ftp
->
rI2_eJ2
[
i
]
);
}
/* Reverse gradient calculation result for edge in other direction */
static
void
rev_gradient
(
FVECT
grad
)
{
grad
[
0
]
=
-
grad
[
0
];
grad
[
1
]
=
-
grad
[
1
];
grad
[
2
]
=
-
grad
[
2
];
}
/* Add to displacement gradient from the given triangle */
static
void
add2gradient
(
FVECT
grad
,
FVECT
egrad1
,
FVECT
egrad2
,
FVECT
egrad3
,
double
v
)
{
int
i
;
for
(
i
=
3
;
i
--
;
)
grad
[
i
]
+=
v
*
(
egrad1
[
i
]
+
egrad2
[
i
]
+
egrad3
[
i
]
);
}
/* Compute anisotropic radii and eigenvector directions */
static
void
eigenvectors
(
FVECT
uv
[
2
],
float
ra
[
2
],
FVECT
hessian
[
3
])
{
double
hess2
[
2
][
2
];
FVECT
a
,
b
;
double
evalue
[
2
],
slope1
,
xmag1
;
int
i
;
/* project Hessian to sample plane */
for
(
i
=
3
;
i
--
;
)
{
a
[
i
]
=
DOT
(
hessian
[
i
],
uv
[
0
]);
b
[
i
]
=
DOT
(
hessian
[
i
],
uv
[
1
]);
}
hess2
[
0
][
0
]
=
DOT
(
uv
[
0
],
a
);
hess2
[
0
][
1
]
=
DOT
(
uv
[
0
],
b
);
hess2
[
1
][
0
]
=
DOT
(
uv
[
1
],
a
);
hess2
[
1
][
1
]
=
DOT
(
uv
[
1
],
b
);
/* compute eigenvalue(s) */
i
=
quadratic
(
evalue
,
1.0
,
-
hess2
[
0
][
0
]
-
hess2
[
1
][
1
],
hess2
[
0
][
0
]
*
hess2
[
1
][
1
]
-
hess2
[
0
][
1
]
*
hess2
[
1
][
0
]);
if
(
i
==
1
)
/* double-root (circle) */
evalue
[
1
]
=
evalue
[
0
];
if
(
!
i
||
((
evalue
[
0
]
=
fabs
(
evalue
[
0
]))
<=
FTINY
*
FTINY
)
|
((
evalue
[
1
]
=
fabs
(
evalue
[
1
]))
<=
FTINY
*
FTINY
)
)
{
ra
[
0
]
=
ra
[
1
]
=
maxarad
;
return
;
}
if
(
evalue
[
0
]
>
evalue
[
1
])
{
ra
[
0
]
=
sqrt
(
sqrt
(
4.0
/
evalue
[
0
]));
ra
[
1
]
=
sqrt
(
sqrt
(
4.0
/
evalue
[
1
]));
slope1
=
evalue
[
1
];
}
else
{
ra
[
0
]
=
sqrt
(
sqrt
(
4.0
/
evalue
[
1
]));
ra
[
1
]
=
sqrt
(
sqrt
(
4.0
/
evalue
[
0
]));
slope1
=
evalue
[
0
];
}
/* compute unit eigenvectors */
if
(
fabs
(
hess2
[
0
][
1
])
<=
FTINY
)
return
;
/* uv OK as is */
slope1
=
(
slope1
-
hess2
[
0
][
0
])
/
hess2
[
0
][
1
];
xmag1
=
sqrt
(
1.0
/
(
1.0
+
slope1
*
slope1
));
for
(
i
=
3
;
i
--
;
)
{
b
[
i
]
=
xmag1
*
uv
[
0
][
i
]
+
slope1
*
xmag1
*
uv
[
1
][
i
];
a
[
i
]
=
slope1
*
xmag1
*
uv
[
0
][
i
]
-
xmag1
*
uv
[
1
][
i
];
}
VCOPY
(
uv
[
0
],
a
);
VCOPY
(
uv
[
1
],
b
);
}
static
void
ambHessian
(
/* anisotropic radii & pos. gradient */
AMBHEMI
*
hp
,
FVECT
uv
[
2
],
/* returned */
float
ra
[
2
],
/* returned (optional) */
float
pg
[
2
]
/* returned (optional) */
)
{
static
char
memerrmsg
[]
=
"out of memory in ambHessian()"
;
FVECT
(
*
hessrow
)[
3
]
=
NULL
;
FVECT
*
gradrow
=
NULL
;
FVECT
hessian
[
3
];
FVECT
gradient
;
FFTRI
fftr
;
int
i
,
j
;
/* be sure to assign unit vectors */
VCOPY
(
uv
[
0
],
hp
->
ux
);
VCOPY
(
uv
[
1
],
hp
->
uy
);
/* clock-wise vertex traversal from sample POV */
if
(
ra
!=
NULL
)
{
/* initialize Hessian row buffer */
hessrow
=
(
FVECT
(
*
)[
3
])
malloc
(
sizeof
(
FVECT
)
*
3
*
(
hp
->
ns
-
1
));
if
(
hessrow
==
NULL
)
error
(
SYSTEM
,
memerrmsg
);
memset
(
hessian
,
0
,
sizeof
(
hessian
));
}
else
if
(
pg
==
NULL
)
/* bogus call? */
return
;
if
(
pg
!=
NULL
)
{
/* initialize form factor row buffer */
gradrow
=
(
FVECT
*
)
malloc
(
sizeof
(
FVECT
)
*
(
hp
->
ns
-
1
));
if
(
gradrow
==
NULL
)
error
(
SYSTEM
,
memerrmsg
);
memset
(
gradient
,
0
,
sizeof
(
gradient
));
}
/* compute first row of edges */
for
(
j
=
0
;
j
<
hp
->
ns
-
1
;
j
++
)
{
comp_fftri
(
&
fftr
,
hp
,
AI
(
hp
,
0
,
j
),
AI
(
hp
,
0
,
j
+
1
));
if
(
hessrow
!=
NULL
)
comp_hessian
(
hessrow
[
j
],
&
fftr
,
hp
->
rp
->
ron
);
if
(
gradrow
!=
NULL
)
comp_gradient
(
gradrow
[
j
],
&
fftr
,
hp
->
rp
->
ron
);
}
/* sum each row of triangles */
for
(
i
=
0
;
i
<
hp
->
ns
-
1
;
i
++
)
{
FVECT
hesscol
[
3
];
/* compute first vertical edge */
FVECT
gradcol
;
comp_fftri
(
&
fftr
,
hp
,
AI
(
hp
,
i
,
0
),
AI
(
hp
,
i
+
1
,
0
));
if
(
hessrow
!=
NULL
)
comp_hessian
(
hesscol
,
&
fftr
,
hp
->
rp
->
ron
);
if
(
gradrow
!=
NULL
)
comp_gradient
(
gradcol
,
&
fftr
,
hp
->
rp
->
ron
);
for
(
j
=
0
;
j
<
hp
->
ns
-
1
;
j
++
)
{
FVECT
hessdia
[
3
];
/* compute triangle contributions */
FVECT
graddia
;
double
backg
;
backg
=
back_ambval
(
hp
,
AI
(
hp
,
i
,
j
),
AI
(
hp
,
i
,
j
+
1
),
AI
(
hp
,
i
+
1
,
j
));
/* diagonal (inner) edge */
comp_fftri
(
&
fftr
,
hp
,
AI
(
hp
,
i
,
j
+
1
),
AI
(
hp
,
i
+
1
,
j
));
if
(
hessrow
!=
NULL
)
{
comp_hessian
(
hessdia
,
&
fftr
,
hp
->
rp
->
ron
);
rev_hessian
(
hesscol
);
add2hessian
(
hessian
,
hessrow
[
j
],
hessdia
,
hesscol
,
backg
);
}
if
(
gradrow
!=
NULL
)
{
comp_gradient
(
graddia
,
&
fftr
,
hp
->
rp
->
ron
);
rev_gradient
(
gradcol
);
add2gradient
(
gradient
,
gradrow
[
j
],
graddia
,
gradcol
,
backg
);
}
/* initialize edge in next row */
comp_fftri
(
&
fftr
,
hp
,
AI
(
hp
,
i
+
1
,
j
+
1
),
AI
(
hp
,
i
+
1
,
j
));
if
(
hessrow
!=
NULL
)
comp_hessian
(
hessrow
[
j
],
&
fftr
,
hp
->
rp
->
ron
);
if
(
gradrow
!=
NULL
)
comp_gradient
(
gradrow
[
j
],
&
fftr
,
hp
->
rp
->
ron
);
/* new column edge & paired triangle */
backg
=
back_ambval
(
hp
,
AI
(
hp
,
i
+
1
,
j
+
1
),
AI
(
hp
,
i
+
1
,
j
),
AI
(
hp
,
i
,
j
+
1
));
comp_fftri
(
&
fftr
,
hp
,
AI
(
hp
,
i
,
j
+
1
),
AI
(
hp
,
i
+
1
,
j
+
1
));
if
(
hessrow
!=
NULL
)
{
comp_hessian
(
hesscol
,
&
fftr
,
hp
->
rp
->
ron
);
rev_hessian
(
hessdia
);
add2hessian
(
hessian
,
hessrow
[
j
],
hessdia
,
hesscol
,
backg
);
if
(
i
<
hp
->
ns
-
2
)
rev_hessian
(
hessrow
[
j
]);
}
if
(
gradrow
!=
NULL
)
{
comp_gradient
(
gradcol
,
&
fftr
,
hp
->
rp
->
ron
);
rev_gradient
(
graddia
);
add2gradient
(
gradient
,
gradrow
[
j
],
graddia
,
gradcol
,
backg
);
if
(
i
<
hp
->
ns
-
2
)
rev_gradient
(
gradrow
[
j
]);
}
}
}
/* release row buffers */
if
(
hessrow
!=
NULL
)
free
(
hessrow
);
if
(
gradrow
!=
NULL
)
free
(
gradrow
);
if
(
ra
!=
NULL
)
/* extract eigenvectors & radii */
eigenvectors
(
uv
,
ra
,
hessian
);
if
(
pg
!=
NULL
)
{
/* tangential position gradient */
pg
[
0
]
=
DOT
(
gradient
,
uv
[
0
]);
pg
[
1
]
=
DOT
(
gradient
,
uv
[
1
]);
}
}
/* Compute direction gradient from a hemispherical sampling */
static
void
ambdirgrad
(
AMBHEMI
*
hp
,
FVECT
uv
[
2
],
float
dg
[
2
])
{
AMBSAMP
*
ap
;
double
dgsum
[
2
];
int
n
;
FVECT
vd
;
double
gfact
;
dgsum
[
0
]
=
dgsum
[
1
]
=
0.0
;
/* sum values times -tan(theta) */
for
(
ap
=
hp
->
sa
,
n
=
hp
->
ns
*
hp
->
ns
;
n
--
;
ap
++
)
{
/* use vector for azimuth + 90deg */
VSUB
(
vd
,
ap
->
p
,
hp
->
rp
->
rop
);
/* brightness over cosine factor */
gfact
=
colval
(
ap
->
v
,
CIEY
)
/
DOT
(
hp
->
rp
->
ron
,
vd
);
/* sine = proj_radius/vd_length */
dgsum
[
0
]
-=
DOT
(
uv
[
1
],
vd
)
*
gfact
;
dgsum
[
1
]
+=
DOT
(
uv
[
0
],
vd
)
*
gfact
;
}
dg
[
0
]
=
dgsum
[
0
]
/
(
hp
->
ns
*
hp
->
ns
);
dg
[
1
]
=
dgsum
[
1
]
/
(
hp
->
ns
*
hp
->
ns
);
}
/* Compute potential light leak direction flags for cache value */
static
uint32
ambcorral
(
AMBHEMI
*
hp
,
FVECT
uv
[
2
],
const
double
r0
,
const
double
r1
)
{
const
double
max_d
=
1.0
/
(
minarad
*
ambacc
+
0.001
);
const
double
ang_res
=
0.5
*
PI
/
hp
->
ns
;
const
double
ang_step
=
ang_res
/
((
int
)(
16
/
PI
*
ang_res
)
+
1.01
);
double
avg_d
=
0
;
uint32
flgs
=
0
;
FVECT
vec
;
double
u
,
v
;
double
ang
,
a1
;
int
i
,
j
;
/* don't bother for a few samples */
if
(
hp
->
ns
<
8
)
return
(
0
);
/* check distances overhead */
for
(
i
=
hp
->
ns
*
3
/
4
;
i
--
>
hp
->
ns
>>
2
;
)
for
(
j
=
hp
->
ns
*
3
/
4
;
j
--
>
hp
->
ns
>>
2
;
)
avg_d
+=
ambsam
(
hp
,
i
,
j
).
d
;
avg_d
*=
4.0
/
(
hp
->
ns
*
hp
->
ns
);
if
(
avg_d
*
r0
>=
1.0
)
/* ceiling too low for corral? */
return
(
0
);
if
(
avg_d
>=
max_d
)
/* insurance */
return
(
0
);
/* else circle around perimeter */
for
(
i
=
0
;
i
<
hp
->
ns
;
i
++
)
for
(
j
=
0
;
j
<
hp
->
ns
;
j
+=
!
i
|
(
i
==
hp
->
ns
-
1
)
?
1
:
hp
->
ns
-
1
)
{
AMBSAMP
*
ap
=
&
ambsam
(
hp
,
i
,
j
);
if
((
ap
->
d
<=
FTINY
)
|
(
ap
->
d
>=
max_d
))
continue
;
/* too far or too near */
VSUB
(
vec
,
ap
->
p
,
hp
->
rp
->
rop
);
u
=
DOT
(
vec
,
uv
[
0
]);
v
=
DOT
(
vec
,
uv
[
1
]);
if
((
r0
*
r0
*
u
*
u
+
r1
*
r1
*
v
*
v
)
*
ap
->
d
*
ap
->
d
<=
u
*
u
+
v
*
v
)
continue
;
/* occluder outside ellipse */
ang
=
atan2a
(
v
,
u
);
/* else set direction flags */
for
(
a1
=
ang
-
ang_res
;
a1
<=
ang
+
ang_res
;
a1
+=
ang_step
)
flgs
|=
1L
<<
(
int
)(
16
/
PI
*
(
a1
+
2.
*
PI
*
(
a1
<
0
)));
}
return
(
flgs
);
}
int
doambient
(
/* compute ambient component */
COLOR
rcol
,
/* input/output color */
RAY
*
r
,
double
wt
,
FVECT
uv
[
2
],
/* returned (optional) */
float
ra
[
2
],
/* returned (optional) */
float
pg
[
2
],
/* returned (optional) */
float
dg
[
2
],
/* returned (optional) */
uint32
*
crlp
/* returned (optional) */
)
{
AMBHEMI
*
hp
=
samp_hemi
(
rcol
,
r
,
wt
);
FVECT
my_uv
[
2
];
double
d
,
K
;
AMBSAMP
*
ap
;
int
i
;
/* clear return values */
if
(
uv
!=
NULL
)
memset
(
uv
,
0
,
sizeof
(
FVECT
)
*
2
);
if
(
ra
!=
NULL
)
ra
[
0
]
=
ra
[
1
]
=
0.0
;
if
(
pg
!=
NULL
)
pg
[
0
]
=
pg
[
1
]
=
0.0
;
if
(
dg
!=
NULL
)
dg
[
0
]
=
dg
[
1
]
=
0.0
;
if
(
crlp
!=
NULL
)
*
crlp
=
0
;
if
(
hp
==
NULL
)
/* sampling falure? */
return
(
0
);
if
((
ra
==
NULL
)
&
(
pg
==
NULL
)
&
(
dg
==
NULL
)
||
(
hp
->
sampOK
<
0
)
|
(
hp
->
ns
<
MINADIV
))
{
free
(
hp
);
/* Hessian not requested/possible */
return
(
-
1
);
/* value-only return value */
}
if
((
d
=
bright
(
rcol
))
>
FTINY
)
{
/* normalize Y values */
d
=
0.99
*
(
hp
->
ns
*
hp
->
ns
)
/
d
;
K
=
0.01
;
}
else
{
/* or fall back on geometric Hessian */
K
=
1.0
;
pg
=
NULL
;
dg
=
NULL
;
crlp
=
NULL
;
}
ap
=
hp
->
sa
;
/* relative Y channel from here on... */
for
(
i
=
hp
->
ns
*
hp
->
ns
;
i
--
;
ap
++
)
colval
(
ap
->
v
,
CIEY
)
=
bright
(
ap
->
v
)
*
d
+
K
;
if
(
uv
==
NULL
)
/* make sure we have axis pointers */
uv
=
my_uv
;
/* compute radii & pos. gradient */
ambHessian
(
hp
,
uv
,
ra
,
pg
);
if
(
dg
!=
NULL
)
/* compute direction gradient */
ambdirgrad
(
hp
,
uv
,
dg
);
if
(
ra
!=
NULL
)
{
/* scale/clamp radii */
if
(
pg
!=
NULL
)
{
if
(
ra
[
0
]
*
(
d
=
fabs
(
pg
[
0
]))
>
1.0
)
ra
[
0
]
=
1.0
/
d
;
if
(
ra
[
1
]
*
(
d
=
fabs
(
pg
[
1
]))
>
1.0
)
ra
[
1
]
=
1.0
/
d
;
if
(
ra
[
0
]
>
ra
[
1
])
ra
[
0
]
=
ra
[
1
];
}
if
(
ra
[
0
]
<
minarad
)
{
ra
[
0
]
=
minarad
;
if
(
ra
[
1
]
<
minarad
)
ra
[
1
]
=
minarad
;
}
ra
[
0
]
*=
d
=
1.0
/
sqrt
(
wt
);
if
((
ra
[
1
]
*=
d
)
>
2.0
*
ra
[
0
])
ra
[
1
]
=
2.0
*
ra
[
0
];
if
(
ra
[
1
]
>
maxarad
)
{
ra
[
1
]
=
maxarad
;
if
(
ra
[
0
]
>
maxarad
)
ra
[
0
]
=
maxarad
;
}
/* flag encroached directions */
if
(
crlp
!=
NULL
)
*
crlp
=
ambcorral
(
hp
,
uv
,
ra
[
0
]
*
ambacc
,
ra
[
1
]
*
ambacc
);
if
(
pg
!=
NULL
)
{
/* cap gradient if necessary */
d
=
pg
[
0
]
*
pg
[
0
]
*
ra
[
0
]
*
ra
[
0
]
+
pg
[
1
]
*
pg
[
1
]
*
ra
[
1
]
*
ra
[
1
];
if
(
d
>
1.0
)
{
d
=
1.0
/
sqrt
(
d
);
pg
[
0
]
*=
d
;
pg
[
1
]
*=
d
;
}
}
}
free
(
hp
);
/* clean up and return */
return
(
1
);
}
Event Timeline
Log In to Comment