Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61808217
pmapdens.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, May 9, 02:34
Size
19 KB
Mime Type
text/x-c
Expires
Sat, May 11, 02:34 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17561144
Attached To
R10977 RADIANCE Photon Map
pmapdens.c
View Options
#ifndef lint
static
const
char
RCSid
[]
=
"$Id$"
;
#endif
/*
======================================================================
Photon map density estimation routines
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")
(c) Tokyo University of Science,
supported by the JSPS Grants-in-Aid for Scientific Research
(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
======================================================================
$Id$
*/
#include "pmapdens.h"
#include "pmapbias.h"
#include "otypes.h"
#include "random.h"
/* Photon map lookup functions per type */
void
(
*
pmapLookup
[
NUM_PMAP_TYPES
])(
PhotonMap
*
,
RAY
*
,
COLOR
)
=
{
photonDensity
,
photonPreCompDensity
,
photonDensity
,
volumePhotonDensity
,
photonDensity
,
photonDensity
,
#ifdef PMAP_PHOTONFLOW
lightFlowPhotonDensity
#endif
};
void
photonDensity
(
PhotonMap
*
pmap
,
RAY
*
ray
,
COLOR
irrad
)
/* Surface-bound photon density estimate for global and caustic photons.
Returns irradiance at ray -> rop. */
{
unsigned
i
;
float
r2
;
COLOR
flux
;
Photon
*
photon
;
const
PhotonSearchQueueNode
*
sqn
;
setcolor
(
irrad
,
0
,
0
,
0
);
if
(
!
pmap
->
maxGather
)
return
;
/* Ignore sources */
if
(
ray
->
ro
&&
islight
(
objptr
(
ray
->
ro
->
omod
)
->
otype
))
return
;
findPhotons
(
pmap
,
ray
);
/* Need at least 2 photons */
if
(
pmap
->
squeue
.
tail
<
2
)
{
#ifdef PMAP_NONEFOUND
sprintf
(
errmsg
,
"no photons found on %s at (%.3f, %.3f, %.3f)"
,
ray
->
ro
?
ray
->
ro
->
oname
:
"<null>"
,
ray
->
rop
[
0
],
ray
->
rop
[
1
],
ray
->
rop
[
2
]
);
error
(
WARNING
,
errmsg
);
#endif
return
;
}
if
(
pmap
->
minGather
==
pmap
->
maxGather
)
{
/* No bias compensation. Just do a plain vanilla estimate */
sqn
=
pmap
->
squeue
.
node
+
1
;
/* Average radius^2 between furthest two photons to improve accuracy */
r2
=
max
(
sqn
->
dist2
,
(
sqn
+
1
)
->
dist2
);
r2
=
0.25
*
(
pmap
->
maxDist2
+
r2
+
2
*
sqrt
(
pmap
->
maxDist2
*
r2
));
/* Skip the extra photon */
for
(
i
=
1
;
i
<
pmap
->
squeue
.
tail
;
i
++
,
sqn
++
)
{
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
getPhotonFlux
(
photon
,
flux
);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on photon dist */
scalecolor
(
flux
,
2
*
(
1
-
sqn
->
dist2
/
r2
));
#endif
addcolor
(
irrad
,
flux
);
}
/* Divide by search area PI * r^2, 1 / PI required as ambient
normalisation factor */
scalecolor
(
irrad
,
1
/
(
PI
*
PI
*
r2
));
return
;
}
else
/* Apply bias compensation to density estimate */
biasComp
(
pmap
,
irrad
);
}
void
photonPreCompDensity
(
PhotonMap
*
pmap
,
RAY
*
r
,
COLOR
irrad
)
/* Surface-bound photon density estimate for (single) precomputed photon.
Returns precomputed irradiance at ray -> rop. */
{
Photon
p
;
setcolor
(
irrad
,
0
,
0
,
0
);
/* Ignore sources */
if
(
r
->
ro
&&
islight
(
objptr
(
r
->
ro
->
omod
)
->
otype
))
return
;
if
(
find1Photon
(
preCompPmap
,
r
,
&
p
))
/* p contains a found photon, so get its irradiance, otherwise it
* remains zero under the assumption all photons are too distant
* to contribute significantly */
getPhotonFlux
(
&
p
,
irrad
);
}
void
volumePhotonDensity
(
PhotonMap
*
pmap
,
RAY
*
ray
,
COLOR
irrad
)
/* Volume photon density estimate using Henyey-Greenstein phase function.
Returns inscattered irradiance at ray -> rop. */
{
unsigned
i
;
float
r2
,
gecc2
,
ph
;
COLOR
flux
;
Photon
*
photon
;
const
PhotonSearchQueueNode
*
sqn
;
setcolor
(
irrad
,
0
,
0
,
0
);
if
(
!
pmap
->
maxGather
)
return
;
findPhotons
(
pmap
,
ray
);
/* Need at least 2 photons */
if
(
pmap
->
squeue
.
tail
<
2
)
return
;
gecc2
=
ray
->
gecc
*
ray
->
gecc
;
sqn
=
pmap
->
squeue
.
node
+
1
;
/* Average radius^2 between furthest two photons to improve accuracy */
r2
=
max
(
sqn
->
dist2
,
(
sqn
+
1
)
->
dist2
);
r2
=
0.25
*
(
pmap
->
maxDist2
+
r2
+
2
*
sqrt
(
pmap
->
maxDist2
*
r2
));
/* Skip the extra photon */
for
(
i
=
1
;
i
<
pmap
->
squeue
.
tail
;
i
++
,
sqn
++
)
{
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
/* Compute phase function for inscattering from photon */
if
(
gecc2
<=
FTINY
)
ph
=
1
;
else
{
ph
=
DOT
(
ray
->
rdir
,
photon
->
norm
)
/
127
;
ph
=
1
+
gecc2
-
2
*
ray
->
gecc
*
ph
;
ph
=
(
1
-
gecc2
)
/
(
ph
*
sqrt
(
ph
));
}
getPhotonFlux
(
photon
,
flux
);
scalecolor
(
flux
,
ph
);
addcolor
(
irrad
,
flux
);
}
/* Divide by search volume 4 / 3 * PI * r^3 and phase function
normalization factor 1 / (4 * PI) */
scalecolor
(
irrad
,
3
/
(
16
*
PI
*
PI
*
r2
*
sqrt
(
r2
)));
}
#ifdef PMAP_PHOTONFLOW
void
volumePhotonSphIrrad
(
PhotonMap
*
pmap
,
RAY
*
ray
,
COLOR
irrad
)
/* Evaluate mean spherical irradiance from volume photons inscattered at
* ray->rop, ignoring participating medium. This evaluates the scalar
* irradiance of a physical light field represented by "photon flow". */
{
unsigned
i
;
float
r2
;
COLOR
flux
;
Photon
*
photon
;
const
PhotonSearchQueueNode
*
sqn
;
setcolor
(
irrad
,
0
,
0
,
0
);
if
(
!
pmap
->
maxGather
)
return
;
findPhotons
(
pmap
,
ray
);
/* Need at least 2 photons */
if
(
pmap
->
squeue
.
tail
<
2
)
return
;
/* Issue warning if few too photons; this will particularly happen
if PMAP_PHOTONFILT is enabled, since photons from duplicate paths
are pruned during the search, and this is not accounted for by
the automatic search radius adjustment. As a workaround for now,
it's up to the user provide a sufficiently large fixed search
radius with the -ma option to rtrace */
if
(
pmap
->
squeue
.
tail
<
pmap
->
squeue
.
len
)
{
sprintf
(
errmsg
,
"short lookup in photon flow; consider -am %.3f or greater"
,
2
*
sqrt
(
pmap
->
maxDist2
)
);
error
(
WARNING
,
errmsg
);
}
sqn
=
pmap
->
squeue
.
node
+
1
;
/* Average radius^2 between furthest two photons to improve accuracy */
r2
=
max
(
sqn
->
dist2
,
(
sqn
+
1
)
->
dist2
);
r2
=
0.25
*
(
pmap
->
maxDist2
+
r2
+
2
*
sqrt
(
pmap
->
maxDist2
*
r2
));
/* Skip the extra photon */
for
(
i
=
1
;
i
<
pmap
->
squeue
.
tail
;
i
++
,
sqn
++
)
{
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
getPhotonFlux
(
photon
,
flux
);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on photon dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor
(
flux
,
1
-
sqn
->
dist2
/
r2
);
#endif
addcolor
(
irrad
,
flux
);
#ifdef PMAP_DEBUGPATHS
printf
(
"%f
\t
%f
\t
%f
\t
%f
\t
%f
\t
%d:%010d
\n
"
,
sqrt
(
sqn
->
dist2
),
photon
->
pos
[
0
],
photon
->
pos
[
1
],
photon
->
pos
[
2
],
flux
[
0
],
photon
->
proc
,
photon
->
primary
);
#endif
}
#ifdef PMAP_DEBUGPATHS
printf
(
"N = %d
\t
r = %f
\n
"
,
pmap
->
squeue
.
tail
,
sqrt
(
r2
));
#endif
#ifdef PMAP_PATHFILT
/* Divide accumulated flux by sphere surface area and ambient
normalisation factor PI */
scalecolor
(
irrad
,
1
/
(
4
*
PI
*
PI
*
r2
));
#else
/* Divide accumulated flux by sphere volume and ambient normalisation
factor PI */
scalecolor
(
irrad
,
3
/
(
4
*
PI
*
PI
*
r2
*
sqrt
(
r2
)));
#endif
}
void
volumePhotonCubicIllum
(
PhotonMap
*
pmap
,
RAY
*
ray
,
COLOR
illum
)
/* Evaluate cubic illuminance according to Cuttle from volume photons
* inscattered at ray->rop, ignoring participating medium. The evaluation
* interprets the volume photon map as a representation of the physical
* light field ("photon flow"). */
{
unsigned
ax
,
p
,
i
,
j
;
float
r2
,
cosNorm
,
volNorm
;
Photon
*
photon
;
const
PhotonSearchQueueNode
*
sqn
;
COLOR
photonFlux
,
E
[
6
],
Evec
[
3
],
Esymm
[
3
],
EvecLum
,
EsymmSum
=
{
0
,
0
,
0
};
FVECT
norm0
;
const
FVECT
norm
[
6
]
=
{
/* Cube face normals (INVERTED since volume photon directions
point TOWARDS their halfspace of incidence */
{
-
1
,
0
,
0
},
{
1
,
0
,
0
},
{
0
,
-
1
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
-
1
},
{
0
,
0
,
1
}
};
setcolor
(
illum
,
0
,
0
,
0
);
if
(
!
pmap
->
maxGather
)
return
;
/* Save ray normal */
VCOPY
(
norm0
,
ray
->
ron
);
#if 1
/* Perform separate lookup for each axis. This improves accuracy
particularly for low bandwitdhs, as the radius adjusts to the
density of incident photons along each axis. This comes at the
expense of 6 lookups.*/
for
(
ax
=
0
;
ax
<
6
;
ax
++
)
{
/* Set cube face normal for current axis */
VCOPY
(
ray
->
ron
,
norm
[
ax
]);
setcolor
(
E
[
ax
],
0
,
0
,
0
);
/* Find photons incident from this axis, filtering by normal */
findPhotons
(
pmap
,
ray
);
/* Need at least 2 photons, else skip */
if
(
pmap
->
squeue
.
tail
<
2
)
continue
;
if
(
pmap
->
squeue
.
tail
<
pmap
->
squeue
.
len
)
{
sprintf
(
errmsg
,
"short lookup in photon flow; consider -am %.3f or greater"
,
2
*
sqrt
(
pmap
->
maxDist2
)
);
error
(
WARNING
,
errmsg
);
}
sqn
=
pmap
->
squeue
.
node
+
1
;
/* Avg radius^2 between furthest two photons to improve accuracy */
r2
=
max
(
sqn
->
dist2
,
(
sqn
+
1
)
->
dist2
);
r2
=
0.25
*
(
pmap
->
maxDist2
+
r2
+
2
*
sqrt
(
pmap
->
maxDist2
*
r2
));
#ifdef _PMAP_DEBUGPATHS
printf
(
"N = %d
\t
r = %f
\n
"
,
pmap
->
squeue
.
tail
,
sqrt
(
r2
));
#endif
/* Skip the extra photon */
for
(
p
=
1
;
p
<
pmap
->
squeue
.
tail
;
p
++
,
sqn
++
)
{
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
getPhotonFlux
(
photon
,
photonFlux
);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor
(
photonFlux
,
1
-
sqn
->
dist2
/
r2
);
#endif
/* Project photon dir onto axis */
cosNorm
=
DOT
(
norm
[
ax
],
photon
->
norm
)
/
127
;
scalecolor
(
photonFlux
,
cosNorm
);
addcolor
(
E
[
ax
],
photonFlux
);
#ifdef PMAP_DEBUGPATHS
printf
(
"%f
\t
%f
\t
%f
\t
%f
\t
%f
\t
%d:%010d
\n
"
,
sqrt
(
sqn
->
dist2
),
photon
->
pos
[
0
],
photon
->
pos
[
1
],
photon
->
pos
[
2
],
photonFlux
[
0
],
photon
->
proc
,
photon
->
primary
);
#endif
}
/* Normalise illuminance for spherical volume */
scalecolor
(
E
[
ax
],
3
/
(
4
*
PI
*
r2
*
sqrt
(
r2
)));
printf
(
"E [%s%d] = %.1f
\n
"
,
ax
&
1
?
"-"
:
"+"
,
ax
>>
1
,
luminance
(
E
[
ax
])
);
}
#else
/* Perform 1 lookup for all axes, and distribute photon flux
according incident direction. While faster than separate lookups
per axis, the accuracy is poor along those axes with few incident
photons, particularly in conjunction with low bandwidths. */
findPhotons
(
pmap
,
ray
);
/* Need at least 2 photons */
if
(
pmap
->
squeue
.
tail
<
2
)
return
;
if
(
pmap
->
squeue
.
tail
<
pmap
->
squeue
.
len
)
{
sprintf
(
errmsg
,
"short lookup in photon flow; consider -am %.3f or greater"
,
2
*
sqrt
(
pmap
->
maxDist2
)
);
error
(
WARNING
,
errmsg
);
}
sqn
=
pmap
->
squeue
.
node
+
1
;
/* Avg radius^2 between furthest two photons to improve accuracy */
r2
=
max
(
sqn
->
dist2
,
(
sqn
+
1
)
->
dist2
);
r2
=
0.25
*
(
pmap
->
maxDist2
+
r2
+
2
*
sqrt
(
pmap
->
maxDist2
*
r2
));
/* Volumetric normalisation factor */
#ifdef PMAP_PATHFILT
volNorm
=
1
/
(
2
*
PI
*
r2
);
#else
# if 1
volNorm
=
3
/
(
4
*
PI
*
r2
*
sqrt
(
r2
));
#else
/* Hemispherical volume */
volNorm
=
3
/
(
2
*
PI
*
r2
*
sqrt
(
r2
));
#endif
#endif
#ifdef PMAP_DEBUGPATHS
printf
(
"N = %d
\t
r = %f
\n
"
,
pmap
->
squeue
.
tail
,
sqrt
(
r2
));
#endif
for
(
ax
=
0
;
ax
<
6
;
ax
++
)
{
/* Set cube face normal for current axis */
setcolor
(
E
[
ax
],
0
,
0
,
0
);
}
/* Skip the extra photon */
for
(
p
=
1
;
p
<
pmap
->
squeue
.
tail
;
p
++
,
sqn
++
)
{
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
getPhotonFlux
(
photon
,
photonFlux
);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor
(
photonFlux
,
1
-
sqn
->
dist2
/
r2
);
#endif
for
(
ax
=
0
;
ax
<
6
;
ax
++
)
{
cosNorm
=
DOT
(
norm
[
ax
],
photon
->
norm
)
/
127
;
if
(
cosNorm
>
0
)
{
/* Project photon dir onto axis */
scalecolor
(
photonFlux
,
cosNorm
);
addcolor
(
E
[
ax
],
photonFlux
);
}
}
#ifdef PMAP_DEBUGPATHS
printf
(
"%f
\t
%f
\t
%f
\t
%f
\t
%f
\t
%d:%010d
\n
"
,
sqrt
(
sqn
->
dist2
),
photon
->
pos
[
0
],
photon
->
pos
[
1
],
photon
->
pos
[
2
],
photonFlux
[
0
],
photon
->
proc
,
photon
->
primary
);
#endif
}
for
(
ax
=
0
;
ax
<
6
;
ax
++
)
{
/* Normalise illuminance for spherical volume */
scalecolor
(
E
[
ax
],
volNorm
);
printf
(
"E [%s%d] = %.1f
\n
"
,
ax
&
1
?
"-"
:
"+"
,
ax
>>
1
,
luminance
(
E
[
ax
])
);
}
#endif
/* Calculate illuminance vector and symmetric illuminance. Unlike
Cuttle's original formulation, we defer scalar reduction and
preserve RGB components in chromatic analysis is required. */
for
(
i
=
0
;
i
<
3
;
i
++
)
{
for
(
j
=
0
;
j
<
3
;
j
++
)
{
Evec
[
i
][
j
]
=
E
[
i
<<
1
][
j
]
-
E
[(
i
<<
1
)
+
1
][
j
];
Esymm
[
i
][
j
]
=
0.5
*
(
E
[
i
<<
1
][
j
]
+
E
[(
i
<<
1
)
+
1
][
j
]
-
fabs
(
Evec
[
i
][
j
])
);
}
/* Reduce Evec to luminance per axis, sum Esymm over axes */
EvecLum
[
i
]
=
luminance
(
Evec
[
i
]);
addcolor
(
EsymmSum
,
Esymm
[
i
]);
}
scalecolor
(
EsymmSum
,
1.0
/
3
);
printf
(
"Evec
\t
= [%.1f
\t
%.1f
\t
%.1f]
\n
"
,
EvecLum
[
0
],
EvecLum
[
1
],
EvecLum
[
2
]
);
printf
(
"Esym
\t
= [%.1f
\t
%.1f
\t
%.1f]
\n
"
,
luminance
(
Esymm
[
0
]),
luminance
(
Esymm
[
1
]),
luminance
(
Esymm
[
2
])
);
/* Reduce to scalars and return as triplet:
- illuminance vector magnitude |Evec|
- scalar symmetry |Esymm|
- scalar illuminance |Esr|. */
illum
[
0
]
=
sqrt
(
DOT
(
EvecLum
,
EvecLum
));
illum
[
1
]
=
luminance
(
EsymmSum
);
illum
[
2
]
=
0.25
*
illum
[
0
]
+
illum
[
1
];
/* Normalise by RADIANCE-specific irradiance factor PI */
scalecolor
(
illum
,
1
/
PI
);
/* Restore ray normal */
VCOPY
(
ray
->
ron
,
norm0
);
}
#endif
/* PMAP_PHOTONFLOW */
#ifdef PMAP_PHOTONFLOW
void
lightFlowPhotonDensity
(
PhotonMap
*
pmap
,
RAY
*
ray
,
COLOR
irrad
)
/* Evaluate irradiance from volume photons at ray->rop for plane defined
* by normal ray -> ron, ignoring participating medium. The evaluation
* interprets the volume photon map as a representation of the "flow"
* of light, i.e. a physical light field. */
{
unsigned
p
,
i
;
float
r2
,
cosNorm
;
Photon
*
photon
;
const
PhotonSearchQueueNode
*
sqn
;
COLOR
photonFlux
;
FVECT
photonDir
;
setcolor
(
irrad
,
0
,
0
,
0
);
if
(
!
pmap
->
maxGather
)
return
;
/* Flip normal since volume photon directions point towards plane of
incidence */
flipsurface
(
ray
);
/* Find photons with incident direction filtered by dot product with
normal (note that found photons may also be located _behind_ the
plane, as if having passed through it. */
findPhotons
(
pmap
,
ray
);
if
(
pmap
->
squeue
.
tail
<
pmap
->
squeue
.
len
)
{
sprintf
(
errmsg
,
"short lookup in photon flow; consider -am %.3f or greater"
,
2
*
sqrt
(
pmap
->
maxDist2
)
);
error
(
WARNING
,
errmsg
);
}
/* Need at least 2 photons, else bail out */
if
(
pmap
->
squeue
.
tail
<
2
)
return
;
sqn
=
pmap
->
squeue
.
node
+
1
;
/* Avg radius^2 between furthest two photons to improve accuracy */
r2
=
max
(
sqn
->
dist2
,
(
sqn
+
1
)
->
dist2
);
r2
=
0.25
*
(
pmap
->
maxDist2
+
r2
+
2
*
sqrt
(
pmap
->
maxDist2
*
r2
));
#ifdef PMAP_DEBUGPATHS
printf
(
"N = %d
\t
r = %f
\n
"
,
pmap
->
squeue
.
tail
,
sqrt
(
r2
));
#endif
/* Skip the extra photon */
for
(
p
=
1
;
p
<
pmap
->
squeue
.
tail
;
p
++
,
sqn
++
)
{
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
getPhotonFlux
(
photon
,
photonFlux
);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor
(
photonFlux
,
1
-
sqn
->
dist2
/
r2
);
#endif
/* Jitter photon direction since stored as 8-bit signed int */
for
(
i
=
0
;
i
<
3
;
i
++
)
{
photonDir
[
i
]
=
photon
->
norm
[
i
];
if
(
photon
->
norm
[
i
]
&&
photon
->
norm
[
i
]
&
127
<
127
)
photonDir
[
i
]
+=
photon
->
norm
[
i
]
&
0x80
?
frandom
()
:
-
frandom
();
}
/* Project photon direction onto normal */
cosNorm
=
DOT
(
ray
->
ron
,
photonDir
)
/
127
;
scalecolor
(
photonFlux
,
cosNorm
);
addcolor
(
irrad
,
photonFlux
);
#ifdef PMAP_DEBUGPATHS
printf
(
"%f
\t
%f
\t
%f
\t
%f
\t
%f
\t
%d:%010d
\n
"
,
sqrt
(
sqn
->
dist2
),
photon
->
pos
[
0
],
photon
->
pos
[
1
],
photon
->
pos
[
2
],
photonFlux
[
0
],
photon
->
proc
,
photon
->
primary
);
#endif
}
/* Normalise accumulated flux by spherical volume, include
RADIANCE-specific irradiance factor PI */
scalecolor
(
irrad
,
3
/
(
4
*
PI
*
PI
*
r2
*
sqrt
(
r2
)));
/* Restore ray normal */
flipsurface
(
ray
);
}
#endif
/* PMAP_PHOTONFLOW */
Event Timeline
Log In to Comment