Skip to content
体验新版
项目
组织
正在加载...
登录
切换导航
打开侧边栏
imjiangjun
pbrt-v4
提交
8d7bcfa2
P
pbrt-v4
项目概览
imjiangjun
/
pbrt-v4
12 个月 前同步成功
通知
7
Star
0
Fork
0
代码
文件
提交
分支
Tags
贡献者
分支图
Diff
Issue
0
列表
看板
标记
里程碑
合并请求
0
DevOps
流水线
流水线任务
计划
Wiki
0
Wiki
分析
仓库
DevOps
项目成员
Pages
P
pbrt-v4
项目概览
项目概览
详情
发布
仓库
仓库
文件
提交
分支
标签
贡献者
分支图
比较
Issue
0
Issue
0
列表
看板
标记
里程碑
合并请求
0
合并请求
0
Pages
DevOps
DevOps
流水线
流水线任务
计划
分析
分析
仓库分析
DevOps
Wiki
0
Wiki
成员
成员
收起侧边栏
关闭侧边栏
动态
分支图
创建新Issue
流水线任务
提交
Issue看板
体验新版 GitCode,发现更多精彩内容 >>
提交
8d7bcfa2
编写于
11月 28, 2020
作者:
M
Matt Pharr
浏览文件
操作
浏览文件
下载
电子邮件补丁
差异文件
Update from book source. No meaningful changes to functionality.
上级
17297225
变更
5
隐藏空白更改
内联
并排
Showing
5 changed file
with
155 addition
and
164 deletion
+155
-164
src/pbrt/gpu/pathintegrator.cpp
src/pbrt/gpu/pathintegrator.cpp
+0
-1
src/pbrt/lights.cpp
src/pbrt/lights.cpp
+51
-56
src/pbrt/lights.h
src/pbrt/lights.h
+34
-30
src/pbrt/util/sampling.h
src/pbrt/util/sampling.h
+62
-69
src/pbrt/util/sampling_test.cpp
src/pbrt/util/sampling_test.cpp
+8
-8
未找到文件。
src/pbrt/gpu/pathintegrator.cpp
浏览文件 @
8d7bcfa2
...
...
@@ -526,7 +526,6 @@ void GPUPathIntegrator::TraceShadowRays(int depth) {
accel
->
IntersectShadowTr
(
maxQueueSize
,
shadowRayQueue
,
&
pixelSampleState
);
else
accel
->
IntersectShadow
(
maxQueueSize
,
shadowRayQueue
,
&
pixelSampleState
);
// Reset shadow ray queue
GPUDo
(
"Reset shadowRayQueue"
,
PBRT_GPU_LAMBDA
()
{
...
...
src/pbrt/lights.cpp
浏览文件 @
8d7bcfa2
...
...
@@ -968,7 +968,7 @@ std::string ImageInfiniteLight::ToString() const {
// PortalImageInfiniteLight Method Definitions
PortalImageInfiniteLight
::
PortalImageInfiniteLight
(
const
Transform
&
renderFromLight
,
Image
equ
i
AreaImage
,
const
Transform
&
renderFromLight
,
Image
equ
al
AreaImage
,
const
RGBColorSpace
*
imageColorSpace
,
Float
scale
,
const
std
::
string
&
filename
,
std
::
vector
<
Point3f
>
p
,
Allocator
alloc
)
:
LightBase
(
LightType
::
Infinite
,
renderFromLight
,
MediumInterface
()),
...
...
@@ -977,8 +977,7 @@ PortalImageInfiniteLight::PortalImageInfiniteLight(
scale
(
scale
),
filename
(
filename
),
distribution
(
alloc
)
{
// Initialize sampling PDFs for infinite area light
ImageChannelDesc
channelDesc
=
equiAreaImage
.
GetChannelDesc
({
"R"
,
"G"
,
"B"
});
ImageChannelDesc
channelDesc
=
equalAreaImage
.
GetChannelDesc
({
"R"
,
"G"
,
"B"
});
if
(
!
channelDesc
)
ErrorExit
(
"%s: image used for PortalImageInfiniteLight doesn't have R, "
"G, B channels."
,
...
...
@@ -986,17 +985,17 @@ PortalImageInfiniteLight::PortalImageInfiniteLight(
CHECK_EQ
(
3
,
channelDesc
.
size
());
CHECK
(
channelDesc
.
IsIdentity
());
if
(
equ
iAreaImage
.
Resolution
().
x
!=
equi
AreaImage
.
Resolution
().
y
)
if
(
equ
alAreaImage
.
Resolution
().
x
!=
equal
AreaImage
.
Resolution
().
y
)
ErrorExit
(
"%s: image resolution (%d, %d) is non-square. It's unlikely "
"this is an "
"equirect environment map."
,
filename
,
equiAreaImage
.
Resolution
().
x
,
equiAreaImage
.
Resolution
().
y
);
"this is an equal area environment map."
,
filename
,
equalAreaImage
.
Resolution
().
x
,
equalAreaImage
.
Resolution
().
y
);
if
(
p
.
size
()
!=
4
)
ErrorExit
(
"Expected 4 vertices for infinite light portal but given %d"
,
p
.
size
());
for
(
int
i
=
0
;
i
<
4
;
++
i
)
portal
[
i
]
=
p
[
i
];
// PortalImageInfiniteLight constructor conclusion
// Compute frame for portal coordinate system
Vector3f
p01
=
Normalize
(
portal
[
1
]
-
portal
[
0
]);
Vector3f
p12
=
Normalize
(
portal
[
2
]
-
portal
[
1
]);
...
...
@@ -1009,37 +1008,36 @@ PortalImageInfiniteLight::PortalImageInfiniteLight(
if
(
std
::
abs
(
Dot
(
p01
,
p12
))
>
.001
||
std
::
abs
(
Dot
(
p12
,
p32
))
>
.001
||
std
::
abs
(
Dot
(
p32
,
p03
))
>
.001
||
std
::
abs
(
Dot
(
p03
,
p01
))
>
.001
)
Error
(
"Infinite light portal isn't a planar quadrilateral"
);
portalFrame
=
Frame
::
FromXY
(
p0
1
,
p03
);
portalFrame
=
Frame
::
FromXY
(
p0
3
,
p01
);
// Resample environment map into rectified coordinates
// Resample the latlong map into rectified coordinates
image
=
Image
(
PixelFormat
::
Float
,
equiAreaImage
.
Resolution
(),
{
"R"
,
"G"
,
"B"
},
equiAreaImage
.
Encoding
(),
alloc
);
// Resample environment map into rectified image
image
=
Image
(
PixelFormat
::
Float
,
equalAreaImage
.
Resolution
(),
{
"R"
,
"G"
,
"B"
},
equalAreaImage
.
Encoding
(),
alloc
);
ParallelFor
(
0
,
image
.
Resolution
().
y
,
[
&
](
int
y
)
{
for
(
int
x
=
0
;
x
<
image
.
Resolution
().
x
;
++
x
)
{
// [0,1]^2 image coordinates
Point2f
st
((
x
+
0.5
f
)
/
image
.
Resolution
().
x
,
// Resample _equalAreaImage_ to compute rectified image pixel $(x,y)$
// Find $(u,v)$ coordinates in equal-area image for pixel
Point2f
uv
((
x
+
0.5
f
)
/
image
.
Resolution
().
x
,
(
y
+
0.5
f
)
/
image
.
Resolution
().
y
);
Vector3f
w
=
RenderFromImage
(
st
);
Vector3f
w
=
RenderFromImage
(
uv
);
w
=
Normalize
(
renderFromLight
.
ApplyInverse
(
w
));
Point2f
uvEqui
=
EqualAreaSphereToSquare
(
w
);
Point2f
stEqui
=
EqualAreaSphereToSquare
(
w
);
for
(
int
c
=
0
;
c
<
3
;
++
c
)
image
.
SetChannel
(
{
x
,
y
},
c
,
equiAreaImage
.
BilerpChannel
(
stEqui
,
c
,
WrapMode
::
OctahedralSphere
));
for
(
int
c
=
0
;
c
<
3
;
++
c
)
{
Float
v
=
equalAreaImage
.
BilerpChannel
(
uvEqui
,
c
,
WrapMode
::
OctahedralSphere
);
image
.
SetChannel
({
x
,
y
},
c
,
v
);
}
}
});
// Initialize sampling
PDFs for infinite area
light
auto
duvdw
=
[
&
](
const
Point2f
&
p
)
{
// Initialize sampling
distribution for portal image infinite
light
auto
duv
_
dw
=
[
&
](
const
Point2f
&
p
)
{
Float
duv_dw
;
(
void
)
RenderFromImage
(
p
,
&
duv_dw
);
return
duv_dw
;
};
Array2D
<
Float
>
d
=
image
.
GetSamplingDistribution
(
duvdw
);
Array2D
<
Float
>
d
=
image
.
GetSamplingDistribution
(
duv
_
dw
);
distribution
=
WindowedPiecewiseConstant2D
(
d
,
alloc
);
}
...
...
@@ -1070,64 +1068,61 @@ SampledSpectrum PortalImageInfiniteLight::Phi(const SampledWavelengths &lambda)
SampledSpectrum
PortalImageInfiniteLight
::
Le
(
const
Ray
&
ray
,
const
SampledWavelengths
&
lambda
)
const
{
// Ignore world to light...
Vector3f
w
=
Normalize
(
ray
.
d
);
Point2f
st
=
ImageFromRender
(
w
);
if
(
!
Inside
(
st
,
ImageBounds
(
ray
.
o
)))
pstd
::
optional
<
Point2f
>
uv
=
ImageFromRender
(
Normalize
(
ray
.
d
));
pstd
::
optional
<
Bounds2f
>
b
=
ImageBounds
(
ray
.
o
);
if
(
!
uv
||
!
b
||
!
Inside
(
*
uv
,
*
b
))
return
SampledSpectrum
(
0.
f
);
return
ImageLookup
(
st
,
lambda
);
return
ImageLookup
(
*
uv
,
lambda
);
}
SampledSpectrum
PortalImageInfiniteLight
::
ImageLookup
(
const
Point2f
&
st
,
const
SampledWavelengths
&
lambda
)
const
{
Point2f
uv
,
const
SampledWavelengths
&
lambda
)
const
{
RGB
rgb
;
for
(
int
c
=
0
;
c
<
3
;
++
c
)
rgb
[
c
]
=
image
.
LookupNearestChannel
(
st
,
c
);
return
scale
*
RGBIlluminantSpectrum
(
*
imageColorSpace
,
ClampZero
(
rgb
)).
Sample
(
lambda
);
rgb
[
c
]
=
image
.
LookupNearestChannel
(
uv
,
c
);
RGBIlluminantSpectrum
spec
(
*
imageColorSpace
,
ClampZero
(
rgb
));
return
scale
*
spec
.
Sample
(
lambda
);
}
pstd
::
optional
<
LightLiSample
>
PortalImageInfiniteLight
::
SampleLi
(
LightSampleContext
ctx
,
Point2f
u
,
SampledWavelengths
lambda
,
LightSamplingMode
mode
)
const
{
Bounds2f
b
=
ImageBounds
(
ctx
.
p
());
// Find $(u,v)$ sample coordinates in infinite light texture
// Sample $(u,v)$ in potentially-visible region of light image
pstd
::
optional
<
Bounds2f
>
b
=
ImageBounds
(
ctx
.
p
());
if
(
!
b
)
return
{};
Float
mapPDF
;
Point2f
uv
=
distribution
.
Sample
(
u
,
b
,
&
mapPDF
);
Point2f
uv
=
distribution
.
Sample
(
u
,
*
b
,
&
mapPDF
);
if
(
mapPDF
==
0
)
return
{};
// Convert infinite light sample point to direction
// Note: ignore WorldToLight since we already folded it in when we
// resampled...
// Convert portal image sample point to direction and compute PDF
Float
duv_dw
;
Vector3f
wi
=
RenderFromImage
(
uv
,
&
duv_dw
);
if
(
duv_dw
==
0
)
return
{};
// Compute PDF for sampled infinite light direction
Float
pdf
=
mapPDF
/
duv_dw
;
CHECK
(
!
IsInf
(
pdf
));
// Compute radiance for portal light sample and return _LightLiSample_
SampledSpectrum
L
=
ImageLookup
(
uv
,
lambda
);
return
LightLiSample
(
L
,
wi
,
pdf
,
Interaction
(
ctx
.
p
()
+
wi
*
(
2
*
sceneRadius
),
&
mediumInterface
));
Point3f
pl
=
ctx
.
p
()
+
2
*
sceneRadius
*
wi
;
return
LightLiSample
(
L
,
wi
,
pdf
,
Interaction
(
pl
,
&
mediumInterface
));
}
Float
PortalImageInfiniteLight
::
PDF_Li
(
LightSampleContext
ctx
,
Vector3f
w
,
LightSamplingMode
mode
)
const
{
// Note: ignore WorldToLight since we already folded it in when we
// resampled...
// Find image $(u,v)$ coordinates corresponding to direction _w_
Float
duv_dw
;
Point2f
st
=
ImageFromRender
(
w
,
&
duv_dw
);
if
(
duv_dw
==
0
)
pstd
::
optional
<
Point2f
>
uv
=
ImageFromRender
(
w
,
&
duv_dw
);
if
(
!
uv
||
duv_dw
==
0
)
return
0
;
Bounds2f
b
=
ImageBounds
(
ctx
.
p
());
Float
pdf
=
distribution
.
PDF
(
st
,
b
);
// Return PDF for sampling $(u,v)$ from reference point
pstd
::
optional
<
Bounds2f
>
b
=
ImageBounds
(
ctx
.
p
());
if
(
!
b
)
return
{};
Float
pdf
=
distribution
.
PDF
(
*
uv
,
*
b
);
return
pdf
/
duv_dw
;
}
...
...
@@ -1181,15 +1176,15 @@ void PortalImageInfiniteLight::PDF_Le(const Ray &ray, Float *pdfPos,
// TODO: negate here or???
Vector3f
w
=
-
Normalize
(
ray
.
d
);
Float
duv_dw
;
Point2f
st
=
ImageFromRender
(
w
,
&
duv_dw
);
pstd
::
optional
<
Point2f
>
uv
=
ImageFromRender
(
w
,
&
duv_dw
);
if
(
duv_dw
==
0
)
{
if
(
!
uv
||
duv_dw
==
0
)
{
*
pdfPos
=
*
pdfDir
=
0
;
return
;
}
Bounds2f
b
(
Point2f
(
0
,
0
),
Point2f
(
1
,
1
));
Float
pdf
=
distribution
.
PDF
(
st
,
b
);
Float
pdf
=
distribution
.
PDF
(
*
uv
,
b
);
#if 0
Normal3f n = Normal3f(portalFrame.z);
...
...
src/pbrt/lights.h
浏览文件 @
8d7bcfa2
...
...
@@ -706,42 +706,46 @@ class PortalImageInfiniteLight : public LightBase {
private:
// PortalImageInfiniteLight Private Methods
PBRT_CPU_GPU
SampledSpectrum
ImageLookup
(
const
Point2f
&
st
,
const
SampledWavelengths
&
lambda
)
const
;
SampledSpectrum
ImageLookup
(
Point2f
uv
,
const
SampledWavelengths
&
lambda
)
const
;
PBRT_CPU_GPU
Vector3f
RenderFromImage
(
const
Point2f
&
st
,
Float
*
duv_dw
=
nullptr
)
const
{
Float
alpha
=
-
Pi
/
2
+
st
.
x
*
Pi
,
beta
=
-
Pi
/
2
+
st
.
y
*
Pi
;
Float
x
=
std
::
tan
(
alpha
),
y
=
std
::
tan
(
beta
);
DCHECK
(
!
IsInf
(
x
)
&&
!
IsInf
(
y
));
Vector3f
w
=
Normalize
(
Vector3f
(
x
,
y
,
-
1
));
if
(
w
.
z
==
0
)
w
.
z
=
1e-5
;
if
(
duv_dw
)
*
duv_dw
=
Pi
*
Pi
*
std
::
abs
((
1
-
w
.
y
*
w
.
y
)
*
(
1
-
w
.
x
*
w
.
x
)
/
w
.
z
);
return
portalFrame
.
FromLocal
(
w
);
}
PBRT_CPU_GPU
Point2f
ImageFromRender
(
const
Vector3f
&
wRender
,
Float
*
duv_dw
=
nullptr
)
const
{
pstd
::
optional
<
Point2f
>
ImageFromRender
(
Vector3f
wRender
,
Float
*
duv_dw
=
nullptr
)
const
{
Vector3f
w
=
portalFrame
.
ToLocal
(
wRender
);
if
(
w
.
z
==
0
)
w
.
z
=
1e-5
;
if
(
w
.
z
<=
0
)
return
{};
// Compute Jacobian determinant of mapping $\roman{d}(u,v)/\roman{d}\omega$ if
// needed
if
(
duv_dw
)
*
duv_dw
=
Pi
*
Pi
*
std
::
abs
((
1
-
w
.
y
*
w
.
y
)
*
(
1
-
w
.
x
*
w
.
x
)
/
w
.
z
)
;
*
duv_dw
=
Sqr
(
Pi
)
*
(
1
-
Sqr
(
w
.
x
))
*
(
1
-
Sqr
(
w
.
y
))
/
w
.
z
;
Float
alpha
=
std
::
atan
(
w
.
x
/
-
w
.
z
),
beta
=
std
::
atan
(
w
.
y
/
-
w
.
z
);
Float
alpha
=
std
::
atan
2
(
w
.
x
,
w
.
z
),
beta
=
std
::
atan2
(
w
.
y
,
w
.
z
);
DCHECK
(
!
IsNaN
(
alpha
+
beta
));
return
Point2f
(
Clamp
((
alpha
+
Pi
/
2
)
/
Pi
,
0
,
1
),
Clamp
((
beta
+
Pi
/
2
)
/
Pi
,
0
,
1
));
}
PBRT_CPU_GPU
Bounds2f
ImageBounds
(
const
Point3f
&
p
)
const
{
Point2f
p0
=
ImageFromRender
(
Normalize
(
portal
[
0
]
-
p
));
Point2f
p1
=
ImageFromRender
(
Normalize
(
portal
[
2
]
-
p
));
return
Bounds2f
(
p0
,
p1
);
Vector3f
RenderFromImage
(
Point2f
uv
,
Float
*
duv_dw
=
nullptr
)
const
{
Float
alpha
=
-
Pi
/
2
+
uv
[
0
]
*
Pi
,
beta
=
-
Pi
/
2
+
uv
[
1
]
*
Pi
;
Float
x
=
std
::
tan
(
alpha
),
y
=
std
::
tan
(
beta
);
DCHECK
(
!
IsInf
(
x
)
&&
!
IsInf
(
y
));
Vector3f
w
=
Normalize
(
Vector3f
(
x
,
y
,
1
));
// Compute Jacobian determinant of mapping $\roman{d}(u,v)/\roman{d}\omega$ if
// needed
if
(
duv_dw
)
*
duv_dw
=
Sqr
(
Pi
)
*
(
1
-
Sqr
(
w
.
x
))
*
(
1
-
Sqr
(
w
.
y
))
/
w
.
z
;
return
portalFrame
.
FromLocal
(
w
);
}
PBRT_CPU_GPU
pstd
::
optional
<
Bounds2f
>
ImageBounds
(
const
Point3f
&
p
)
const
{
pstd
::
optional
<
Point2f
>
p0
=
ImageFromRender
(
Normalize
(
portal
[
0
]
-
p
));
pstd
::
optional
<
Point2f
>
p1
=
ImageFromRender
(
Normalize
(
portal
[
2
]
-
p
));
if
(
!
p0
||
!
p1
)
return
{};
return
Bounds2f
(
*
p0
,
*
p1
);
}
PBRT_CPU_GPU
...
...
@@ -750,15 +754,15 @@ class PortalImageInfiniteLight : public LightBase {
}
// PortalImageInfiniteLight Private Members
std
::
string
filename
;
pstd
::
array
<
Point3f
,
4
>
portal
;
Frame
portalFrame
;
Image
image
;
WindowedPiecewiseConstant2D
distribution
;
const
RGBColorSpace
*
imageColorSpace
;
Float
scale
;
Frame
portalFrame
;
pstd
::
array
<
Point3f
,
4
>
portal
;
WindowedPiecewiseConstant2D
distribution
;
Point3f
sceneCenter
;
Float
sceneRadius
;
std
::
string
filename
;
Point3f
sceneCenter
;
};
// SpotLight Definition
...
...
src/pbrt/util/sampling.h
浏览文件 @
8d7bcfa2
...
...
@@ -906,62 +906,51 @@ class SummedAreaTable {
// SummedAreaTable Public Methods
SummedAreaTable
(
Allocator
alloc
)
:
sum
(
alloc
)
{}
SummedAreaTable
(
const
Array2D
<
Float
>
&
values
,
Allocator
alloc
=
{})
:
sum
(
integrate
(
values
,
alloc
))
{}
:
sum
(
values
.
xSize
(),
values
.
ySize
(),
alloc
)
{
sum
(
0
,
0
)
=
values
(
0
,
0
);
// Compute sums along first scanline and column
for
(
int
x
=
1
;
x
<
sum
.
xSize
();
++
x
)
sum
(
x
,
0
)
=
values
(
x
,
0
)
+
sum
(
x
-
1
,
0
);
for
(
int
y
=
1
;
y
<
sum
.
ySize
();
++
y
)
sum
(
0
,
y
)
=
values
(
0
,
y
)
+
sum
(
0
,
y
-
1
);
PBRT_CPU_GPU
Float
Sum
(
const
Bounds2f
&
extent
)
const
{
double
s
=
((
lookup
(
extent
.
pMax
.
x
,
extent
.
pMax
.
y
)
-
lookup
(
extent
.
pMin
.
x
,
extent
.
pMax
.
y
))
+
(
lookup
(
extent
.
pMin
.
x
,
extent
.
pMin
.
y
)
-
lookup
(
extent
.
pMax
.
x
,
extent
.
pMin
.
y
)));
return
std
::
max
<
Float
>
(
s
,
0
);
// Compute sums for the remainder of the entries
for
(
int
y
=
1
;
y
<
sum
.
ySize
();
++
y
)
for
(
int
x
=
1
;
x
<
sum
.
xSize
();
++
x
)
sum
(
x
,
y
)
=
(
values
(
x
,
y
)
+
sum
(
x
-
1
,
y
)
+
sum
(
x
,
y
-
1
)
-
sum
(
x
-
1
,
y
-
1
));
}
PBRT_CPU_GPU
Float
Average
(
const
Bounds2f
&
extent
)
const
{
return
Sum
(
extent
)
/
extent
.
Area
();
}
Float
Integral
(
const
Bounds2f
&
extent
)
const
{
double
s
=
(((
double
)
Lookup
(
extent
.
pMax
.
x
,
extent
.
pMax
.
y
)
-
(
double
)
Lookup
(
extent
.
pMin
.
x
,
extent
.
pMax
.
y
))
+
((
double
)
Lookup
(
extent
.
pMin
.
x
,
extent
.
pMin
.
y
)
-
(
double
)
Lookup
(
extent
.
pMax
.
x
,
extent
.
pMin
.
y
)));
return
std
::
max
<
Float
>
(
s
/
(
sum
.
xSize
()
*
sum
.
ySize
()),
0
);
}
std
::
string
ToString
()
const
;
private:
// SummedAreaTable Private Methods
Array2D
<
double
>
integrate
(
const
Array2D
<
Float
>
&
values
,
Allocator
alloc
)
{
auto
f
=
[
&
values
](
int
x
,
int
y
)
{
return
values
(
x
,
y
)
/
(
values
.
xSize
()
*
values
.
ySize
());
};
Array2D
<
double
>
result
(
values
.
xSize
(),
values
.
ySize
(),
alloc
);
result
(
0
,
0
)
=
f
(
0
,
0
);
// Compute sums along first scanline and column
for
(
int
x
=
1
;
x
<
result
.
xSize
();
++
x
)
result
(
x
,
0
)
=
f
(
x
,
0
)
+
result
(
x
-
1
,
0
);
for
(
int
y
=
1
;
y
<
result
.
ySize
();
++
y
)
result
(
0
,
y
)
=
f
(
0
,
y
)
+
result
(
0
,
y
-
1
);
// Compute sums for the remainder of the entries
for
(
int
y
=
1
;
y
<
result
.
ySize
();
++
y
)
for
(
int
x
=
1
;
x
<
result
.
xSize
();
++
x
)
result
(
x
,
y
)
=
(
f
(
x
,
y
)
+
result
(
x
-
1
,
y
)
+
result
(
x
,
y
-
1
)
-
result
(
x
-
1
,
y
-
1
));
return
result
;
}
PBRT_CPU_GPU
double
l
ookup
(
Float
x
,
Float
y
)
const
{
Float
L
ookup
(
Float
x
,
Float
y
)
const
{
// Rescale $(x,y)$ to table resolution and compute integer coordinates
x
*=
sum
.
xSize
();
y
*=
sum
.
ySize
();
int
x0
=
(
int
)
x
,
y0
=
(
int
)
y
;
// Bilinearly interpolate between surrounding table values
Float
v00
=
lookup
(
x0
,
y0
),
v10
=
lookup
(
x0
+
1
,
y0
);
Float
v01
=
lookup
(
x0
,
y0
+
1
),
v11
=
lookup
(
x0
+
1
,
y0
+
1
);
Float
v00
=
LookupInt
(
x0
,
y0
),
v10
=
LookupInt
(
x0
+
1
,
y0
);
Float
v01
=
LookupInt
(
x0
,
y0
+
1
),
v11
=
LookupInt
(
x0
+
1
,
y0
+
1
);
Float
dx
=
x
-
int
(
x
),
dy
=
y
-
int
(
y
);
return
(
1
-
dx
)
*
(
1
-
dy
)
*
v00
+
(
1
-
dx
)
*
dy
*
v01
+
dx
*
(
1
-
dy
)
*
v10
+
dx
*
dy
*
v11
;
}
PBRT_CPU_GPU
double
lookup
(
int
x
,
int
y
)
const
{
Float
LookupInt
(
int
x
,
int
y
)
const
{
// Return zero at lower boundaries
if
(
x
==
0
||
y
==
0
)
return
0
;
...
...
@@ -987,80 +976,84 @@ class WindowedPiecewiseConstant2D {
PBRT_CPU_GPU
Point2f
Sample
(
const
Point2f
&
u
,
const
Bounds2f
&
b
,
Float
*
pdf
)
const
{
// Handle zero-valued function for windowed sampling
if
(
sat
.
Sum
(
b
)
==
0
)
{
if
(
sat
.
Integral
(
b
)
==
0
)
{
*
pdf
=
0
;
return
{};
}
//
Sample marginal windowed function in the first dimens
ion
Float
sumb
=
sat
.
Sum
(
b
);
//
Define lambda function _Px_ for marginal cumulative distribut
ion
Float
bInt
=
sat
.
Integral
(
b
);
auto
Px
=
[
&
,
this
](
Float
x
)
->
Float
{
Bounds2f
bx
=
b
;
bx
.
pMax
.
x
=
x
;
return
sat
.
Sum
(
bx
)
/
sumb
;
return
sat
.
Integral
(
bx
)
/
bInt
;
};
// Sample marginal windowed function in $x$
Point2f
p
;
p
.
x
=
SampleBisection
(
Px
,
u
[
0
],
b
.
pMin
.
x
,
b
.
pMax
.
x
,
func
.
xSize
());
// Sample conditional windowed function in $y$
// Compute 2D bounds _bCond_ for conditional sampling
int
nx
=
func
.
xSize
();
p
.
x
=
sample
(
Px
,
u
[
0
],
b
.
pMin
.
x
,
b
.
pMax
.
x
,
nx
);
// Sample conditional windowed function in the second dimension
Bounds2f
by
(
Point2f
(
std
::
floor
(
p
.
x
*
nx
)
/
nx
,
b
.
pMin
.
y
),
Point2f
(
std
::
ceil
(
p
.
x
*
nx
)
/
nx
,
b
.
pMax
.
y
));
if
(
by
.
pMin
.
x
==
by
.
pMax
.
x
)
by
.
pMax
.
x
+=
1.
f
/
nx
;
if
(
sat
.
Sum
(
by
)
==
0
)
{
// This can happen when we're provided a really narrow initial
// bounding box
Bounds2f
bCond
(
Point2f
(
std
::
floor
(
p
.
x
*
nx
)
/
nx
,
b
.
pMin
.
y
),
Point2f
(
std
::
ceil
(
p
.
x
*
nx
)
/
nx
,
b
.
pMax
.
y
));
if
(
bCond
.
pMin
.
x
==
bCond
.
pMax
.
x
)
bCond
.
pMax
.
x
+=
1.
f
/
nx
;
if
(
sat
.
Integral
(
bCond
)
==
0
)
{
*
pdf
=
0
;
return
{};
}
Float
sumby
=
sat
.
Sum
(
by
);
// Define lambda function for conditional distribution and sample $y$
Float
condIntegral
=
sat
.
Integral
(
bCond
);
auto
Py
=
[
&
,
this
](
Float
y
)
->
Float
{
Bounds2f
by
y
=
by
;
by
y
.
pMax
.
y
=
y
;
return
sat
.
Sum
(
byy
)
/
sumby
;
Bounds2f
by
=
bCond
;
by
.
pMax
.
y
=
y
;
return
sat
.
Integral
(
by
)
/
condIntegral
;
};
p
.
y
=
sample
(
Py
,
u
[
1
],
b
.
pMin
.
y
,
b
.
pMax
.
y
,
func
.
ySize
());
p
.
y
=
SampleBisection
(
Py
,
u
[
1
],
b
.
pMin
.
y
,
b
.
pMax
.
y
,
func
.
ySize
());
// Compute PDF and return point sampled from windowed function
*
pdf
=
PDF
(
p
,
b
)
;
*
pdf
=
Eval
(
p
)
/
bInt
;
return
p
;
}
PBRT_CPU_GPU
Float
PDF
(
const
Point2f
&
p
,
const
Bounds2f
&
b
)
const
{
if
(
sat
.
Sum
(
b
)
==
0
)
if
(
sat
.
Integral
(
b
)
==
0
)
return
0
;
return
Eval
(
p
)
/
sat
.
Sum
(
b
);
return
Eval
(
p
)
/
sat
.
Integral
(
b
);
}
private:
// WindowedPiecewiseConstant2D Private Methods
PBRT_CPU_GPU
Float
Eval
(
const
Point2f
&
p
)
const
{
Point2i
pi
(
std
::
min
<
int
>
(
p
[
0
]
*
func
.
xSize
(),
func
.
xSize
()
-
1
),
std
::
min
<
int
>
(
p
[
1
]
*
func
.
ySize
(),
func
.
ySize
()
-
1
));
return
func
[
pi
];
}
template
<
typename
F
>
PBRT_CPU_GPU
static
Float
sample
(
F
func
,
Float
u
,
Float
min
,
Float
max
,
int
n
)
{
template
<
typename
CDF
>
PBRT_CPU_GPU
static
Float
SampleBisection
(
CDF
P
,
Float
u
,
Float
min
,
Float
max
,
int
n
)
{
// Apply bisection to bracket _u_
while
(
std
::
ceil
(
n
*
max
)
-
std
::
floor
(
n
*
min
)
>
1
)
{
DCHECK_LE
(
func
(
min
),
u
);
DCHECK_GE
(
func
(
max
),
u
);
DCHECK_LE
(
P
(
min
),
u
);
DCHECK_GE
(
P
(
max
),
u
);
Float
mid
=
(
min
+
max
)
/
2
;
if
(
func
(
mid
)
>
u
)
if
(
P
(
mid
)
>
u
)
max
=
mid
;
else
min
=
mid
;
}
// Find sample by interpolating between _min_ and _max_
Float
t
=
(
u
-
func
(
min
))
/
(
func
(
max
)
-
func
(
min
));
Float
t
=
(
u
-
P
(
min
))
/
(
P
(
max
)
-
P
(
min
));
return
Clamp
(
Lerp
(
t
,
min
,
max
),
min
,
max
);
}
PBRT_CPU_GPU
Float
Eval
(
const
Point2f
&
p
)
const
{
Point2i
pi
(
std
::
min
<
int
>
(
p
[
0
]
*
func
.
xSize
(),
func
.
xSize
()
-
1
),
std
::
min
<
int
>
(
p
[
1
]
*
func
.
ySize
(),
func
.
ySize
()
-
1
));
return
func
[
pi
];
}
// WindowedPiecewiseConstant2D Private Members
SummedAreaTable
sat
;
Array2D
<
Float
>
func
;
...
...
src/pbrt/util/sampling_test.cpp
浏览文件 @
8d7bcfa2
...
...
@@ -1243,11 +1243,11 @@ TEST(SummedArea, Constant) {
SummedAreaTable
sat
(
v
);
EXPECT_EQ
(
1
,
sat
.
Sum
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
1
,
1
))));
EXPECT_EQ
(
0.5
,
sat
.
Sum
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
1
,
0.5
))));
EXPECT_EQ
(
0.5
,
sat
.
Sum
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
0.5
,
1
))));
EXPECT_EQ
(
3.
/
16.
,
sat
.
Sum
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
.25
,
.75
))));
EXPECT_EQ
(
3.
/
16.
,
sat
.
Sum
(
Bounds2f
(
Point2f
(
0.5
,
0.25
),
Point2f
(
0.75
,
1
))));
EXPECT_EQ
(
1
,
sat
.
Integral
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
1
,
1
))));
EXPECT_EQ
(
0.5
,
sat
.
Integral
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
1
,
0.5
))));
EXPECT_EQ
(
0.5
,
sat
.
Integral
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
0.5
,
1
))));
EXPECT_EQ
(
3.
/
16.
,
sat
.
Integral
(
Bounds2f
(
Point2f
(
0
,
0
),
Point2f
(
.25
,
.75
))));
EXPECT_EQ
(
3.
/
16.
,
sat
.
Integral
(
Bounds2f
(
Point2f
(
0.5
,
0.25
),
Point2f
(
0.75
,
1
))));
}
TEST
(
SummedArea
,
Rect
)
{
...
...
@@ -1271,7 +1271,7 @@ TEST(SummedArea, Rect) {
Bounds2f
b
(
Point2f
(
Float
(
x0
)
/
v
.
xSize
(),
Float
(
y0
)
/
v
.
ySize
()),
Point2f
(
Float
(
x1
)
/
v
.
xSize
(),
Float
(
y1
)
/
v
.
ySize
()));
EXPECT_EQ
(
mySum
/
(
v
.
xSize
()
*
v
.
ySize
()),
sat
.
Sum
(
b
));
EXPECT_EQ
(
mySum
/
(
v
.
xSize
()
*
v
.
ySize
()),
sat
.
Integral
(
b
));
}
}
...
...
@@ -1300,7 +1300,7 @@ TEST(SummedArea, Randoms) {
ref
+=
v
[
p
];
ref
/=
v
.
xSize
()
*
v
.
ySize
();
double
s
=
sat
.
Sum
(
bf
);
double
s
=
sat
.
Integral
(
bf
);
if
(
ref
!=
s
)
EXPECT_LT
(
std
::
abs
((
ref
-
s
)
/
ref
),
1e-3
f
)
<<
StringPrintf
(
"ref %f s %f"
,
ref
,
s
);
...
...
@@ -1332,7 +1332,7 @@ TEST(SummedArea, NonCellAligned) {
}
Float
sampledResult
=
sampledSum
*
b
.
Area
()
/
nSamples
;
double
s
=
sat
.
Sum
(
b
);
double
s
=
sat
.
Integral
(
b
);
if
(
sampledResult
!=
s
)
EXPECT_LT
(
std
::
abs
((
sampledResult
-
s
)
/
sampledResult
),
1e-3
f
)
<<
StringPrintf
(
"sampled %f s %f"
,
sampledResult
,
s
);
...
...
编辑
预览
Markdown
is supported
0%
请重试
或
添加新附件
.
添加附件
取消
You are about to add
0
people
to the discussion. Proceed with caution.
先完成此消息的编辑!
取消
想要评论请
注册
或
登录