Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86260600
pmaptkdt.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
Sat, Oct 5, 08:55
Size
14 KB
Mime Type
text/x-c
Expires
Mon, Oct 7, 08:55 (2 d)
Engine
blob
Format
Raw Data
Handle
21348587
Attached To
R10977 RADIANCE Photon Map
pmaptkdt.c
View Options
/*
======================================================================
In-core kd-tree for 4D TRANSIENT photon map. Based on pmapkdt.h with
extensions to handle photon path length (representing time of flight).
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Tokyo University of Science,
supported by the JSPS Grants-in-Aid for Scientific Research
(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
======================================================================
$Id$
*/
#include "pmapdata.h"
/* Includes pmaptkdt.h */
#include "source.h"
#include "otspecial.h"
#include "random.h"
/* PHOTON MAP BUILD ROUTINES ------------------------------------------ */
static
unsigned
long
kdT_TransMedianPartition
(
const
Photon
*
heap
,
unsigned
long
*
heapIdx
,
unsigned
long
*
heapXdi
,
unsigned
long
left
,
unsigned
long
right
,
unsigned
dim
)
/* Returns index to median in heap from indices left to right
(inclusive) in dimension dim [0..3], where dim=3 is path length
(corresponding to time of flight). The heap is partitioned relative to
median using a quicksort algorithm. The heap indices in heapIdx are
sorted rather than the heap itself. */
{
const
float
*
p
;
float
pLen
;
unsigned
long
l
,
r
,
lg2
,
n2
,
m
,
n
=
right
-
left
+
1
;
unsigned
d
;
/* Round down n to nearest power of 2 */
for
(
lg2
=
0
,
n2
=
n
;
n2
>
1
;
n2
>>=
1
,
++
lg2
);
n2
=
1
<<
lg2
;
/* Determine median position; this takes into account the fact that
only the last level in the heap can be partially empty, and that
it fills from left to right */
m
=
left
+
((
n
-
n2
)
>
(
n2
>>
1
)
-
1
?
n2
-
1
:
n
-
(
n2
>>
1
));
while
(
right
>
left
)
{
/* Pivot node */
p
=
heap
[
heapIdx
[
right
]].
pos
;
pLen
=
heap
[
heapIdx
[
right
]].
aux
.
pathLen
;
l
=
left
;
r
=
right
-
1
;
/* l & r converge, swapping elements out of order with respect to
pivot node. Identical keys are resolved by cycling through
dim. The convergence point is then the pivot's position. */
do
{
while
(
l
<=
r
)
{
d
=
dim
;
while
(
d
<
3
?
heap
[
heapIdx
[
l
]].
pos
[
d
]
==
p
[
d
]
:
heap
[
heapIdx
[
l
]].
aux
.
pathLen
==
pLen
)
{
d
=
(
d
+
1
)
%
4
;
if
(
d
==
dim
)
{
/* Ignore dupes? */
error
(
WARNING
,
"duplicate keys in photon heap"
);
l
++
;
break
;
}
}
if
(
d
<
3
?
heap
[
heapIdx
[
l
]].
pos
[
d
]
<
p
[
d
]
:
heap
[
heapIdx
[
l
]].
aux
.
pathLen
<
pLen
)
l
++
;
else
break
;
}
while
(
r
>
l
)
{
d
=
dim
;
while
(
d
<
3
?
heap
[
heapIdx
[
r
]].
pos
[
d
]
==
p
[
d
]
:
heap
[
heapIdx
[
r
]].
aux
.
pathLen
==
pLen
)
{
d
=
(
d
+
1
)
%
4
;
if
(
d
==
dim
)
{
/* Ignore dupes? */
error
(
WARNING
,
"duplicate keys in photon heap"
);
r
--
;
break
;
}
}
if
(
d
<
3
?
heap
[
heapIdx
[
r
]].
pos
[
d
]
>
p
[
d
]
:
heap
[
heapIdx
[
r
]].
aux
.
pathLen
>
pLen
)
r
--
;
else
break
;
}
/* Swap indices (not the nodes they point to) */
n2
=
heapIdx
[
l
];
heapIdx
[
l
]
=
heapIdx
[
r
];
heapIdx
[
r
]
=
n2
;
/* Update reverse indices */
heapXdi
[
heapIdx
[
l
]]
=
l
;
heapXdi
[
n2
]
=
r
;
}
while
(
l
<
r
);
/* Swap indices of convergence and pivot nodes */
heapIdx
[
r
]
=
heapIdx
[
l
];
heapIdx
[
l
]
=
heapIdx
[
right
];
heapIdx
[
right
]
=
n2
;
/* Update reverse indices */
heapXdi
[
heapIdx
[
r
]]
=
r
;
heapXdi
[
heapIdx
[
l
]]
=
l
;
heapXdi
[
n2
]
=
right
;
if
(
l
>=
m
)
right
=
l
-
1
;
if
(
l
<=
m
)
left
=
l
+
1
;
}
/* Once left & right have converged at m, we have found the median */
return
m
;
}
static
void
kdT_TransBuild
(
Photon
*
heap
,
unsigned
long
*
heapIdx
,
unsigned
long
*
heapXdi
,
const
float
min
[
4
],
const
float
max
[
4
],
unsigned
long
left
,
unsigned
long
right
,
unsigned
long
root
)
/* Recursive part of transBalancePhotons(..). Builds heap from subarray
defined by indices left and right. min and max are the minimum resp.
maximum photon positions and path lengths (corresponding to the photons'
time of flight) in the array. root is the index of the current subtree's
root, which corresponds to the median's 1-based index in the heap.
heapIdx are the balanced heap indices. The heap is accessed indirectly
through these. heapXdi are the reverse indices from the heap to heapIdx
so that heapXdi [heapIdx [i]] = i. */
{
float
maxLeft
[
4
],
minRight
[
4
],
maxDist
,
dDist
;
Photon
rootNode
;
unsigned
d
,
dim
;
unsigned
long
median
;
/* Choose median for dimension with largest spread and partition
accordingly */
for
(
maxDist
=
0
,
dim
=
0
,
d
=
0
;
d
<
4
;
d
++
)
if
((
dDist
=
max
[
d
]
-
min
[
d
])
>
maxDist
)
{
maxDist
=
dDist
;
dim
=
d
;
}
median
=
left
==
right
?
left
:
kdT_TransMedianPartition
(
heap
,
heapIdx
,
heapXdi
,
left
,
right
,
dim
);
/* Place median at root of current subtree. This consists of swapping
the median and the root nodes and updating the heap indices */
memcpy
(
&
rootNode
,
heap
+
heapIdx
[
median
],
sizeof
(
Photon
));
memcpy
(
heap
+
heapIdx
[
median
],
heap
+
root
-
1
,
sizeof
(
Photon
));
rootNode
.
discr
=
dim
;
memcpy
(
heap
+
root
-
1
,
&
rootNode
,
sizeof
(
Photon
));
heapIdx
[
heapXdi
[
root
-
1
]]
=
heapIdx
[
median
];
heapXdi
[
heapIdx
[
median
]]
=
heapXdi
[
root
-
1
];
heapIdx
[
median
]
=
root
-
1
;
heapXdi
[
root
-
1
]
=
median
;
/* Update bounds for left and right subtrees and recurse on them */
for
(
d
=
0
;
d
<
4
;
d
++
)
if
(
d
==
dim
)
maxLeft
[
d
]
=
minRight
[
d
]
=
d
<
3
?
rootNode
.
pos
[
d
]
:
rootNode
.
aux
.
pathLen
;
else
{
maxLeft
[
d
]
=
max
[
d
];
minRight
[
d
]
=
min
[
d
];
}
if
(
left
<
median
)
kdT_TransBuild
(
heap
,
heapIdx
,
heapXdi
,
min
,
maxLeft
,
left
,
median
-
1
,
root
<<
1
);
if
(
right
>
median
)
kdT_TransBuild
(
heap
,
heapIdx
,
heapXdi
,
minRight
,
max
,
median
+
1
,
right
,
(
root
<<
1
)
+
1
);
}
void
kdT_TransBuildPhotonMap
(
struct
PhotonMap
*
pmap
)
{
Photon
*
nodes
;
unsigned
long
i
;
unsigned
long
*
heapIdx
,
/* Photon index array */
*
heapXdi
;
/* Reverse index to heapIdx */
float
min
[
4
],
max
[
4
];
/* Allocate kd-tree nodes and load photons from heap file */
if
(
!
(
nodes
=
calloc
(
pmap
->
numPhotons
,
sizeof
(
Photon
))))
error
(
SYSTEM
,
"failed in-core heap allocation in kdT_TransBuildPhotonMap"
);
rewind
(
pmap
->
heap
);
i
=
fread
(
nodes
,
sizeof
(
Photon
),
pmap
->
numPhotons
,
pmap
->
heap
);
if
(
i
!=
pmap
->
numPhotons
)
error
(
SYSTEM
,
"failed loading photon heap in kdT_TransBuildPhotonMap"
);
pmap
->
store
.
nodes
=
nodes
;
heapIdx
=
calloc
(
pmap
->
numPhotons
,
sizeof
(
unsigned
long
));
heapXdi
=
calloc
(
pmap
->
numPhotons
,
sizeof
(
unsigned
long
));
if
(
!
heapIdx
||
!
heapXdi
)
error
(
SYSTEM
,
"failed heap index allocation in kdT_TransBuildPhotonMap"
);
/* Initialize index arrays */
for
(
i
=
0
;
i
<
pmap
->
numPhotons
;
i
++
)
heapXdi
[
i
]
=
heapIdx
[
i
]
=
i
;
/* Build kd-tree */
VCOPY
(
min
,
pmap
->
minPos
);
min
[
3
]
=
pmap
->
minPathLen
;
VCOPY
(
max
,
pmap
->
maxPos
);
max
[
3
]
=
pmap
->
maxPathLen
;
kdT_TransBuild
(
nodes
,
heapIdx
,
heapXdi
,
min
,
max
,
0
,
pmap
->
numPhotons
-
1
,
1
);
/* Cleanup */
free
(
heapIdx
);
free
(
heapXdi
);
}
/* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */
static
void
kdT_TransFindNearest
(
PhotonMap
*
pmap
,
const
float
pos
[
4
],
const
kdT_SearchFilterCallback
*
filter
,
const
kdT_SearchAttribCallback
*
attrib
,
unsigned
long
node
)
/* Recursive part of kdT_TransFindPhotons(). Locate pmap -> squeue.len
* nearest neighbours to 4D pos which pass the specified filt (if not
* NULL) and return in search queue starting at pmap -> squeue.node. */
{
Photon
*
p
=
(
Photon
*
)
pmap
->
store
.
nodes
+
node
-
1
;
/* Signed distance to current photon's splitting plane */
float
d
=
pos
[
p
->
discr
]
-
(
p
->
discr
<
3
?
p
->
pos
[
p
->
discr
]
:
p
->
aux
.
pathLen
),
d2
=
d
*
d
,
dv
[
4
];
/* Search subtree closer to pos first; exclude other subtree if the
distance to the splitting plane is greater than maxDist */
if
(
d
<
0
)
{
if
(
node
<<
1
<=
pmap
->
numPhotons
)
kdT_TransFindNearest
(
pmap
,
pos
,
filter
,
attrib
,
node
<<
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<
pmap
->
numPhotons
)
kdT_TransFindNearest
(
pmap
,
pos
,
filter
,
attrib
,
(
node
<<
1
)
+
1
);
}
else
{
if
(
node
<<
1
<
pmap
->
numPhotons
)
kdT_TransFindNearest
(
pmap
,
pos
,
filter
,
attrib
,
(
node
<<
1
)
+
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<=
pmap
->
numPhotons
)
kdT_TransFindNearest
(
pmap
,
pos
,
filter
,
attrib
,
node
<<
1
);
}
/* Bail out if photon is rejected by filter */
if
(
filter
&&
!
filter
->
filtFunc
(
p
,
filter
->
state
))
return
;
/* Squared distance to current photon (dist2() would require doubles) */
VSUB
(
dv
,
pos
,
p
->
pos
);
d2
=
DOT
(
dv
,
dv
);
dv
[
3
]
=
pos
[
3
]
-
p
->
aux
.
pathLen
;
d2
+=
dv
[
3
]
*
dv
[
3
];
/* Accept photon if closer than current max dist & add to search queue */
if
(
d2
<
pmap
->
maxDist2
)
{
/* Insert record in search queue if it passes filter and its SQUARED
* dist to key is below maxDist2 */
if
((
d2
=
kdT_PutNearest
(
&
pmap
->
squeue
,
d2
,
p
,
attrib
))
>=
0
)
/* Update maxDist2 if search queue is full */
pmap
->
maxDist2
=
d2
;
}
}
int
kdT_TransFindPhotons
(
struct
PhotonMap
*
pmap
,
const
FVECT
pos
,
const
FVECT
norm
)
{
kdT_SearchFilterCallback
filt
;
kdT_SearchFilterState
filtState
;
#ifdef PMAP_PATHFILT
kdT_SearchAttribCallback
paths
,
*
pathsPtr
=
NULL
;
#endif
float
p
[
4
],
n
[
3
];
/* 4th key dimension: photon path length at given time in units of
* scene geometry for consistency with pos */
const
float
pathLen
=
pmap
->
time
*
pmap
->
velocity
;
/* Photon pos & normal stored at lower precision; note path length is
* appended to position as 4th key dimension */
VCOPY
(
p
,
pos
);
p
[
3
]
=
pathLen
;
if
(
norm
)
VCOPY
(
n
,
norm
);
/* Set up filter callback */
filtState
.
pmap
=
pmap
;
filtState
.
pos
=
pos
;
filtState
.
norm
=
norm
?
n
:
NULL
;
filt
.
state
=
&
filtState
;
filt
.
filtFunc
=
filterPhoton
;
#ifdef PMAP_PATHFILT
/* Set up path ID callback to prune photons from duplicate paths */
paths
.
state
=
pmap
;
paths
.
findFunc
=
findPhotonPath
;
paths
.
addFunc
=
addPhotonPath
;
paths
.
delFunc
=
deletePhotonPath
;
paths
.
checkFunc
=
checkPhotonPaths
;
pathsPtr
=
&
paths
;
resetPhotonPaths
(
pmap
);
kdT_TransFindNearest
(
pmap
,
p
,
&
filt
,
pathsPtr
,
1
);
#else
kdT_TransFindNearest
(
pmap
,
p
,
&
filt
,
NULL
,
1
);
#endif
/* PMAP_PATHFILT */
/* Get (maximum distance)^2 from search queue head */
pmap
->
maxDist2
=
pmap
->
squeue
.
node
[
0
].
dist2
;
/* Return success or failure (empty queue => none found) */
return
pmap
->
squeue
.
tail
?
0
:
-
1
;
}
static
void
kdT_TransFind1Nearest
(
PhotonMap
*
pmap
,
const
float
pos
[
4
],
const
kdT_SearchFilterCallback
*
filter
,
Photon
**
photon
,
unsigned
long
node
)
/* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
* 4D pos with similar normal. Note that all heap and queue indices are
* 1-based, but accesses to the arrays are 0-based! */
{
Photon
*
p
=
(
Photon
*
)
pmap
->
store
.
nodes
+
node
-
1
;
/* Signed distance to current photon's splitting plane */
float
d
=
pos
[
p
->
discr
]
-
(
p
->
discr
<
3
?
p
->
pos
[
p
->
discr
]
:
p
->
aux
.
pathLen
),
d2
=
d
*
d
,
dv
[
4
];
/* Search subtree closer to pos first; exclude other subtree if the
distance to the splitting plane is greater than maxDist */
if
(
d
<
0
)
{
if
(
node
<<
1
<=
pmap
->
numPhotons
)
kdT_TransFind1Nearest
(
pmap
,
pos
,
filter
,
photon
,
node
<<
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<
pmap
->
numPhotons
)
kdT_TransFind1Nearest
(
pmap
,
pos
,
filter
,
photon
,
(
node
<<
1
)
+
1
);
}
else
{
if
(
node
<<
1
<
pmap
->
numPhotons
)
kdT_TransFind1Nearest
(
pmap
,
pos
,
filter
,
photon
,
(
node
<<
1
)
+
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<=
pmap
->
numPhotons
)
kdT_TransFind1Nearest
(
pmap
,
pos
,
filter
,
photon
,
node
<<
1
);
}
/* Squared distance to current photon (dist2() would require doubles) */
VSUB
(
dv
,
pos
,
p
->
pos
);
d2
=
DOT
(
dv
,
dv
);
dv
[
3
]
=
pos
[
3
]
-
p
->
aux
.
pathLen
;
d2
+=
dv
[
3
]
*
dv
[
3
];
if
(
d2
<
pmap
->
maxDist2
&&
(
!
filter
||
filter
->
filtFunc
(
p
,
filter
->
state
))
)
{
/* Closest photon so far that passes filter */
pmap
->
maxDist2
=
d2
;
*
photon
=
p
;
}
}
int
kdT_Find1TransPhoton
(
struct
PhotonMap
*
pmap
,
const
FVECT
pos
,
const
FVECT
norm
,
Photon
*
photon
)
{
kdT_SearchFilterCallback
filt
;
kdT_SearchFilterState
filtState
;
float
p
[
4
],
n
[
3
];
Photon
*
pnn
=
NULL
;
/* 4th key dimension: photon path length at given time in units of
* scene geometry for consistency with pos */
const
float
pathLen
=
pmap
->
time
*
pmap
->
velocity
;
/* Photon pos & normal stored at lower precision; note path length is
* appended to position as 4th key dimension */
VCOPY
(
p
,
pos
);
p
[
3
]
=
pathLen
;
if
(
norm
)
VCOPY
(
n
,
norm
);
/* Set up filter callback */
filtState
.
pmap
=
pmap
;
filtState
.
norm
=
norm
?
n
:
NULL
;
filtState
.
srcMod
=
NULL
;
filt
.
state
=
&
filtState
;
filt
.
filtFunc
=
filterPhoton
;
kdT_TransFind1Nearest
(
pmap
,
p
,
&
filt
,
&
pnn
,
1
);
if
(
!
pnn
)
/* No photon found => failed */
return
-
1
;
else
{
/* Copy found photon => successs */
memcpy
(
photon
,
pnn
,
sizeof
(
Photon
));
return
0
;
}
}
Event Timeline
Log In to Comment