Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F74239934
pmapkdt.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
Fri, Jul 26, 16:05
Size
23 KB
Mime Type
text/x-c
Expires
Sun, Jul 28, 16:05 (2 d)
Engine
blob
Format
Raw Data
Handle
19364794
Attached To
R10977 RADIANCE Photon Map
pmapkdt.c
View Options
/*
======================================================================
In-core kd-tree for photon map
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Fraunhofer Institute for Solar Energy Systems,
supported by the German Research Foundation
(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS))
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #147053, "Daylight Redirecting Components",
SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
(c) Tokyo University of Science,
supported by the JSPS Grants-in-Aid for Scientific Research
(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
======================================================================
$Id: pmapkdt.c,v 1.1 2020/11/10 01:10:57 u-no-hoo Exp u-no-hoo $
*/
#include "pmapdata.h"
/* Includes pmapkdt.h */
#include "pmapfilt.h"
#include "source.h"
#include "otspecial.h"
#include "random.h"
/* PHOTON MAP BUILD ROUTINES ------------------------------------------ */
void
kdT_Null
(
PhotonKdTree
*
kdt
)
{
kdt
->
nodes
=
NULL
;
}
static
unsigned
long
kdT_MedianPartition
(
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. 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
;
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
;
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
(
heap
[
heapIdx
[
l
]].
pos
[
d
]
==
p
[
d
])
{
d
=
(
d
+
1
)
%
3
;
if
(
d
==
dim
)
{
/* Ignore dupes? */
error
(
WARNING
,
"duplicate keys in photon heap"
);
l
++
;
break
;
}
}
if
(
heap
[
heapIdx
[
l
]].
pos
[
d
]
<
p
[
d
])
l
++
;
else
break
;
}
while
(
r
>
l
)
{
d
=
dim
;
while
(
heap
[
heapIdx
[
r
]].
pos
[
d
]
==
p
[
d
])
{
d
=
(
d
+
1
)
%
3
;
if
(
d
==
dim
)
{
/* Ignore dupes? */
error
(
WARNING
,
"duplicate keys in photon heap"
);
r
--
;
break
;
}
}
if
(
heap
[
heapIdx
[
r
]].
pos
[
d
]
>
p
[
d
])
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_Build
(
Photon
*
heap
,
unsigned
long
*
heapIdx
,
unsigned
long
*
heapXdi
,
const
float
min
[
3
],
const
float
max
[
3
],
unsigned
long
left
,
unsigned
long
right
,
unsigned
long
root
)
/* Recursive part of balancePhotons(..). Builds heap from subarray
defined by indices left and right. min and max are the minimum resp.
maximum photon positions 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
[
3
],
minRight
[
3
];
Photon
rootNode
;
unsigned
d
;
/* Choose median for dimension with largest spread and partition
accordingly */
const
float
d0
=
max
[
0
]
-
min
[
0
],
d1
=
max
[
1
]
-
min
[
1
],
d2
=
max
[
2
]
-
min
[
2
];
const
unsigned
char
dim
=
d0
>
d1
?
d0
>
d2
?
0
:
2
:
d1
>
d2
?
1
:
2
;
const
unsigned
long
median
=
left
==
right
?
left
:
kdT_MedianPartition
(
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
<=
2
;
d
++
)
if
(
d
==
dim
)
maxLeft
[
d
]
=
minRight
[
d
]
=
rootNode
.
pos
[
d
];
else
{
maxLeft
[
d
]
=
max
[
d
];
minRight
[
d
]
=
min
[
d
];
}
if
(
left
<
median
)
kdT_Build
(
heap
,
heapIdx
,
heapXdi
,
min
,
maxLeft
,
left
,
median
-
1
,
root
<<
1
);
if
(
right
>
median
)
kdT_Build
(
heap
,
heapIdx
,
heapXdi
,
minRight
,
max
,
median
+
1
,
right
,
(
root
<<
1
)
+
1
);
}
void
kdT_BuildPhotonMap
(
struct
PhotonMap
*
pmap
)
{
Photon
*
nodes
;
unsigned
long
i
;
unsigned
long
*
heapIdx
,
/* Photon index array */
*
heapXdi
;
/* Reverse index to heapIdx */
/* 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_BuildPhotonMap"
);
rewind
(
pmap
->
heap
);
i
=
fread
(
nodes
,
sizeof
(
Photon
),
pmap
->
numPhotons
,
pmap
->
heap
);
if
(
i
!=
pmap
->
numPhotons
)
error
(
SYSTEM
,
"failed loading photon heap in kdT_BuildPhotonMap"
);
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_BuildPhotonMap"
);
/* Initialize index arrays */
for
(
i
=
0
;
i
<
pmap
->
numPhotons
;
i
++
)
heapXdi
[
i
]
=
heapIdx
[
i
]
=
i
;
/* Build kd-tree */
kdT_Build
(
nodes
,
heapIdx
,
heapXdi
,
pmap
->
minPos
,
pmap
->
maxPos
,
0
,
pmap
->
numPhotons
-
1
,
1
);
/* Cleanup */
free
(
heapIdx
);
free
(
heapXdi
);
}
/* PHOTON MAP I/O ROUTINES ------------------------------------------ */
int
kdT_SavePhotons
(
const
struct
PhotonMap
*
pmap
,
FILE
*
out
)
{
unsigned
long
i
,
j
;
Photon
*
p
=
(
Photon
*
)
pmap
->
store
.
nodes
;
for
(
i
=
0
;
i
<
pmap
->
numPhotons
;
i
++
,
p
++
)
{
/* Write photon attributes */
for
(
j
=
0
;
j
<
3
;
j
++
)
putflt
(
p
->
pos
[
j
],
out
);
/* Bytewise dump otherwise we have portability probs */
for
(
j
=
0
;
j
<
3
;
j
++
)
putint
(
p
->
norm
[
j
],
1
,
out
);
#ifdef PMAP_FLOAT_FLUX
for
(
j
=
0
;
j
<
3
;
j
++
)
putflt
(
p
->
flux
[
j
],
out
);
#else
for
(
j
=
0
;
j
<
4
;
j
++
)
putint
(
p
->
flux
[
j
],
1
,
out
);
#endif
putbinary
(
&
p
->
aux
,
1
,
sizeof
(
p
->
aux
),
out
);
putint
(
p
->
flags
,
1
,
out
);
if
(
ferror
(
out
))
return
-
1
;
}
return
0
;
}
int
kdT_LoadPhotons
(
struct
PhotonMap
*
pmap
,
FILE
*
in
)
{
unsigned
long
i
,
j
;
Photon
*
p
;
/* Allocate kd-tree based on initialised pmap -> numPhotons */
pmap
->
store
.
nodes
=
calloc
(
sizeof
(
Photon
),
pmap
->
numPhotons
);
if
(
!
pmap
->
store
.
nodes
)
error
(
SYSTEM
,
"failed kd-tree allocation in kdT_LoadPhotons"
);
/* Get photon attributes */
for
(
i
=
0
,
p
=
pmap
->
store
.
nodes
;
i
<
pmap
->
numPhotons
;
i
++
,
p
++
)
{
for
(
j
=
0
;
j
<
3
;
j
++
)
p
->
pos
[
j
]
=
getflt
(
in
);
/* Bytewise grab otherwise we have portability probs */
for
(
j
=
0
;
j
<
3
;
j
++
)
p
->
norm
[
j
]
=
getint
(
1
,
in
);
#ifdef PMAP_FLOAT_FLUX
for
(
j
=
0
;
j
<
3
;
j
++
)
p
->
flux
[
j
]
=
getflt
(
in
);
#else
for
(
j
=
0
;
j
<
4
;
j
++
)
p
->
flux
[
j
]
=
getint
(
1
,
in
);
#endif
getbinary
(
&
p
->
aux
,
1
,
sizeof
(
p
->
aux
),
in
);
p
->
flags
=
getint
(
1
,
in
);
if
(
feof
(
in
))
return
-
1
;
}
return
0
;
}
/* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */
void
kdT_InitFindPhotons
(
struct
PhotonMap
*
pmap
)
{
pmap
->
squeue
.
len
=
pmap
->
maxGather
+
1
;
pmap
->
squeue
.
node
=
calloc
(
pmap
->
squeue
.
len
,
sizeof
(
PhotonSearchQueueNode
)
);
if
(
!
pmap
->
squeue
.
node
)
error
(
SYSTEM
,
"can't allocate photon search queue"
);
#ifdef PMAP_PATHFILT
initPhotonPaths
(
pmap
);
#endif
}
#ifdef DEBUG_KDT_NN
static
int
kdT_CheckSearchQueue
(
kdT_SearchQueue
*
squeue
,
const
float
pos
[
3
],
unsigned
root
)
/* Priority queue sanity check */
{
unsigned
kid
;
const
kdT_SearchQueueNode
*
sqn
=
squeue
->
node
;
const
Photon
*
photon
;
float
d2
,
dv
[
3
];
if
(
root
<
squeue
->
tail
)
{
photon
=
kdT_GetNearestPhoton
(
squeue
,
sqn
[
root
].
idx
);
VSUB
(
dv
,
pos
,
photon
->
pos
);
d2
=
DOT
(
dv
,
dv
);
if
(
sqn
[
root
].
dist2
!=
d2
)
return
-
1
;
if
((
kid
=
(
root
<<
1
)
+
1
)
<
squeue
->
tail
)
{
if
(
sqn
[
kid
].
dist2
>
sqn
[
root
].
dist2
)
return
-
1
;
else
return
kdT_CheckSearchQueue
(
squeue
,
pos
,
kid
);
}
if
(
++
kid
<
squeue
->
tail
)
{
if
(
sqn
[
kid
].
dist2
>
sqn
[
root
].
dist2
)
return
-
1
;
else
return
kdT_CheckSearchQueue
(
squeue
,
pos
,
kid
);
}
}
return
0
;
}
#endif
float
kdT_PutNearest
(
kdT_SearchQueue
*
squeue
,
float
d2
,
const
Photon
*
photon
,
const
kdT_SearchAttribCallback
*
attrib
)
/* Insert new photon with SQUARED distance to query point into the search
* priority queue, maintaining the most distant photon at the queue head.
* If the queue is full, the new photon is only inserted if it is closer
* than the queue head.
*
* The search priority queue is represented as a linearised binary tree with
* the root corresponding to the queue head, and the tail corresponding to
* the last leaf.
*
* The optional attrib callback maps photon attributes to their queue nodes
* by calling attrib->findFunc() if attrib != NULL. It enables finding a
* previously enqueued photon and replacing it if the new photon's SQUARED
* distance d2 is lower.
*
* This function returns the new maximum SQUARED distance at the head if the
* search queue is full. Otherwise it returns -1, indicating a maximum for
* the entire queue is as yet undefined.
*/
{
const
Photon
*
tPhoton
=
NULL
;
kdT_SearchQueueNode
*
sqn
=
squeue
->
node
,
**
attribSqn
=
NULL
,
**
tAttribSqn
=
NULL
;
unsigned
root
,
kid
,
kid2
;
/* Start with undefined max distance */
float
d2max
=
-
1
;
#ifdef PMAP_PATHFILT
/* Find previously enqueued photon with same attribute if callback
defined */
#ifdef PMAP_DEBUGPATHS
printf
(
"-------- Enqueue --------
\n
"
);
#endif
attribSqn
=
(
kdT_SearchQueueNode
**
)(
attrib
?
attrib
->
findFunc
(
photon
,
attrib
->
state
)
:
NULL
);
#ifdef PMAP_DEBUGPATHS
printf
(
attribSqn
?
"found
\n
"
:
"not found
\n
"
);
#endif
#endif
if
(
!
attribSqn
&&
squeue
->
tail
<
squeue
->
len
)
{
/* No duplicate attribute in queue, and queue not full;
insert photon at tail and resort towards head */
kid
=
squeue
->
tail
++
;
while
(
kid
)
{
root
=
(
kid
-
1
)
>>
1
;
/* Compare with parent and swap if necessary (bubble up), else
* terminate */
if
(
d2
>
sqn
[
root
].
dist2
)
{
sqn
[
kid
].
dist2
=
sqn
[
root
].
dist2
;
sqn
[
kid
].
idx
=
sqn
[
root
].
idx
;
#ifdef PMAP_PATHFILT
if
(
attrib
)
{
/* Update root's attribute map entry to refer to kid */
tPhoton
=
kdT_GetNearestPhoton
(
squeue
,
sqn
[
root
].
idx
);
tAttribSqn
=
(
kdT_SearchQueueNode
**
)(
attrib
->
findFunc
(
tPhoton
,
attrib
->
state
)
);
#ifdef PMAP_DEBUGPATHS
printf
(
"update
\n
"
);
#endif
if
(
!
tAttribSqn
)
{
/* Consistency error: a previously enqueued photon must
* have an entry in the attribute map! */
error
(
CONSISTENCY
,
"kdT_PutNearest: corrupt attribute map"
);
exit
(
-
1
);
}
*
tAttribSqn
=
&
sqn
[
kid
];
}
#endif
kid
=
root
;
}
else
break
;
}
/* Priority queue restored; assign new photon's queue entry */
sqn
[
kid
].
dist2
=
d2
;
sqn
[
kid
].
idx
=
(
PhotonIdx
)
photon
;
#ifdef PMAP_PATHFILT
if
(
attrib
)
/* Add new photon to attribute map */
attrib
->
addFunc
(
photon
,
&
sqn
[
kid
],
attrib
->
state
);
#endif
}
else
{
/* Search queue full or duplicate attribute in queue and new photon may
be closer; set root to replaced photon, otherwise to queue head */
root
=
attribSqn
?
*
attribSqn
-
sqn
:
0
;
if
(
d2
<
sqn
[
root
].
dist2
)
{
/* New photon closer than maximum at root; replace and resort towards
* leaves (=search queue tail) */
#ifdef PMAP_PATHFILT
if
(
attrib
&&
!
attribSqn
)
{
/* New photon will replace root and has unique attribute; delete
* root from the attribute map */
tPhoton
=
kdT_GetNearestPhoton
(
squeue
,
sqn
[
root
].
idx
);
attrib
->
delFunc
(
tPhoton
,
attrib
->
state
);
}
#endif
while
((
kid
=
(
root
<<
1
)
+
1
)
<
squeue
->
tail
)
{
/* Compare with larger kid & swap if necessary (bubble down),
* else terminate */
if
(
(
kid2
=
(
kid
+
1
))
<
squeue
->
tail
&&
sqn
[
kid2
].
dist2
>
sqn
[
kid
].
dist2
)
kid
=
kid2
;
if
(
d2
<
sqn
[
kid
].
dist2
)
{
sqn
[
root
].
dist2
=
sqn
[
kid
].
dist2
;
sqn
[
root
].
idx
=
sqn
[
kid
].
idx
;
#ifdef PMAP_PATHFILT
if
(
attrib
)
{
/* Update kid's attribute map entry to refer to root */
tPhoton
=
kdT_GetNearestPhoton
(
squeue
,
sqn
[
kid
].
idx
);
tAttribSqn
=
(
kdT_SearchQueueNode
**
)(
attrib
->
findFunc
(
tPhoton
,
attrib
->
state
)
);
#ifdef PMAP_DEBUGPATHS
printf
(
"update
\n
"
);
#endif
if
(
!
tAttribSqn
)
{
/* Consistency error: a previously enqueued photon must
* have an entry in the attribute map! */
error
(
CONSISTENCY
,
"kdT_PutNearest: corrupt attribute map"
);
exit
(
-
1
);
}
*
tAttribSqn
=
&
sqn
[
root
];
}
#endif
}
else
break
;
root
=
kid
;
}
/* Priority queue restored; reassign head or replaced node */
sqn
[
root
].
dist2
=
d2
;
sqn
[
root
].
idx
=
(
PhotonIdx
)
photon
;
#ifdef PMAP_PATHFILT
if
(
attrib
)
{
if
(
!
attribSqn
)
/* Add new photon to attribute map */
attrib
->
addFunc
(
photon
,
&
sqn
[
root
],
attrib
->
state
);
else
{
#ifdef PMAP_DEBUGPATHS
attrib
->
findFunc
(
photon
,
attrib
->
state
);
printf
(
"update
\n
"
);
#endif
/* Update new photon's existing attribute map entry */
*
attribSqn
=
&
sqn
[
root
];
}
}
#endif
}
}
if
(
squeue
->
tail
>=
squeue
->
len
)
/* Search queue full; update max distance^2 from head node */
d2max
=
sqn
[
0
].
dist2
;
#if defined(PMAP_PATHFILT) && defined (PMAP_DEBUGPATHS)
if
(
attrib
)
/* Run sanity check on attribute map */
attrib
->
checkFunc
(
attrib
->
state
);
#endif
return
d2max
;
}
static
void
kdT_FindNearest
(
PhotonMap
*
pmap
,
const
float
pos
[
3
],
const
kdT_SearchFilterCallback
*
filter
,
const
kdT_SearchAttribCallback
*
attrib
,
unsigned
long
node
)
/* Recursive part of kdT_FindPhotons(). Locate pmap -> squeue.len nearest
* neighbours 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
->
pos
[
p
->
discr
],
d2
=
d
*
d
,
dv
[
3
];
/* 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_FindNearest
(
pmap
,
pos
,
filter
,
attrib
,
node
<<
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<
pmap
->
numPhotons
)
kdT_FindNearest
(
pmap
,
pos
,
filter
,
attrib
,
(
node
<<
1
)
+
1
);
}
else
{
if
(
node
<<
1
<
pmap
->
numPhotons
)
kdT_FindNearest
(
pmap
,
pos
,
filter
,
attrib
,
(
node
<<
1
)
+
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<=
pmap
->
numPhotons
)
kdT_FindNearest
(
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 (note dist2() requires doubles) */
VSUB
(
dv
,
pos
,
p
->
pos
);
d2
=
DOT
(
dv
,
dv
);
/* 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
;
#ifdef DEBUG_KDT_NN
if
(
kdT_CheckSearchQueue
(
&
pmap
->
squeue
,
pos
,
0
))
error
(
CONSISTENCY
,
"kdT_PutNearest: corrupt search queue"
);
#endif
}
}
}
int
kdT_FindPhotons
(
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
[
3
],
n
[
3
];
/* Photon pos & normal stored at lower precision */
VCOPY
(
p
,
pos
);
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
;
/* No contribs with kd-tree */
filtState
.
srcMod
=
NULL
;
#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_FindNearest
(
pmap
,
p
,
&
filt
,
pathsPtr
,
1
);
#else
kdT_FindNearest
(
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_Find1Nearest
(
PhotonMap
*
pmap
,
const
float
pos
[
3
],
const
kdT_SearchFilterCallback
*
filter
,
Photon
**
photon
,
unsigned
long
node
)
/* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
* 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
->
pos
[
p
->
discr
],
d2
=
d
*
d
,
dv
[
3
];
/* 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_Find1Nearest
(
pmap
,
pos
,
filter
,
photon
,
node
<<
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<
pmap
->
numPhotons
)
kdT_Find1Nearest
(
pmap
,
pos
,
filter
,
photon
,
(
node
<<
1
)
+
1
);
}
else
{
if
(
node
<<
1
<
pmap
->
numPhotons
)
kdT_Find1Nearest
(
pmap
,
pos
,
filter
,
photon
,
(
node
<<
1
)
+
1
);
if
(
d2
<
pmap
->
maxDist2
&&
node
<<
1
<=
pmap
->
numPhotons
)
kdT_Find1Nearest
(
pmap
,
pos
,
filter
,
photon
,
node
<<
1
);
}
/* Squared distance to current photon */
VSUB
(
dv
,
pos
,
p
->
pos
);
d2
=
DOT
(
dv
,
dv
);
if
(
d2
<
pmap
->
maxDist2
&&
(
!
filter
||
filter
->
filtFunc
(
p
,
filter
->
state
))
)
{
/* Closest photon so far that passes filter */
pmap
->
maxDist2
=
d2
;
*
photon
=
p
;
}
}
int
kdT_Find1Photon
(
struct
PhotonMap
*
pmap
,
const
FVECT
pos
,
const
FVECT
norm
,
Photon
*
photon
)
{
kdT_SearchFilterCallback
filt
;
kdT_SearchFilterState
filtState
;
float
p
[
3
],
n
[
3
];
/* Photon pos & normal stored at lower precision */
VCOPY
(
p
,
pos
);
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
;
Photon
*
pnn
=
NULL
;
kdT_Find1Nearest
(
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
;
}
}
/* PHOTON ACCESS ROUTINES ------------------------------------------ */
int
kdT_GetPhoton
(
const
struct
PhotonMap
*
pmap
,
PhotonIdx
idx
,
Photon
*
photon
)
{
memcpy
(
photon
,
idx
,
sizeof
(
Photon
));
return
0
;
}
Photon
*
kdT_GetNearestPhoton
(
const
PhotonSearchQueue
*
squeue
,
PhotonIdx
idx
)
{
return
idx
;
}
PhotonIdx
kdT_FirstPhoton
(
const
struct
PhotonMap
*
pmap
)
{
return
pmap
->
store
.
nodes
;
}
void
kdT_Delete
(
PhotonKdTree
*
kdt
)
{
free
(
kdt
->
nodes
);
kdt
->
nodes
=
NULL
;
}
Event Timeline
Log In to Comment