diff --git a/README-photonflow.txt b/README-photonflow.txt new file mode 100644 index 0000000..56f9d14 --- /dev/null +++ b/README-photonflow.txt @@ -0,0 +1,197 @@ +This is README file for RADIANCE photon flow. This software was developed +as part of the project "Three-Dimensional Light Flow" hosted by Tokyo +University of Science, supported by the JSPS Grants-in-Aid for Scientific +Research (KAKENHI) under the grant number JP19KK0115. + +$Revision: 1.3 $ $Date: 2020/11/16 22:41:09 $ + + + +INTRODUCTION + +These source files implement a modified volume photon density estimate in +RADIANCE to evaluate the illuminance of a physical lightfield using photon +flow. This code is an experimental variant of the RADIANCE photon map extension +intended for testing prior to committing it to the RADIANCE CVS as part of the +official distribution. Whether this (admittedly specialised) functionality +should be made available to the RADIANCE community as a build option is +debatable and may be decided in due course. + +This code is essentially a C implementation of the prototype Python script +volPhotonIllum.py. As such, it is orders of magnitude faster, which comes to +bear with large photon maps and many sensor positions. One crucial difference +between the RADIANCE and Python implementation is that the former performs +photon lookups in terms of a fixed number of photons around the sensor +position (nearest neighbour search), while the latter does so in terms of a +fixed radius (range search). The RADIANCE implementation however offers an +additional maximum radius parameter for optimisation, which is infact +required to obtain the specified number of photons in conjunction with +photon flow (see USAGE below). + +The illuminance evaluated from photon flow approximates Cuttle's cubic +illuminance, which in turn quantifies a sensor position in a physical light +field. Photon flow constitues volume photons deposited by mkpmap in a +participating medium (mist); due to the interaction of a photon with the +medium, it will be deposited at multiple locations along its path. Each +photon path has a unique ID consisting of a prefix (mkpmap child process +number if multiprocessing) and a suffix (a 32-bit serial ray number); these +can be output by pmapdump with the -a and -P options. + +Evaluating photon flow at a given position in the lightfield requires pruning +photons with identical path IDs within the search radius, keeping only the +closest photon from each intercepted path. To this end, the search routines use +filtering routines which maintain a lookup table containing the intercepted +path IDs in order to identify duplicates. + +Some code was refactored to filter photons and their path IDs in a new module +pmapfilt.c shared by the two underlying data structures which store the photon +map: the (legacy) kd-tree and the out-of-core octree. The nearest-neighbour +search routines of both variants were modified to use these filter functions +via callbacks to facilitate maintenance. In this way, volume photons can be +filtered for duplicate path IDs "on the fly" as the search progresses, as +opposed to pruning these a posteriori, resulting in a lower photon count. In +addition, a new volume photon density estimate routine volumePhotonSphIrrad() +was added to pmapdens.c to evaluate the photon flow; this routine differs +from the already existing volume photon density estimate routine +volumePhotonDensity() in that it does not apply the Henyey-Greenstein phase +function to compute outscattering, and completely ignores the participating +medium. + +Because photon flow not only includes photons that have been reflected off one +or more surfaces (indirect contributions), but also photons that have just been +emitted from light sources (direct contributions), the new volume photon +density estimate is called by multiDirectPmap() in pmapsrc.c, which normally +calls the direct photon map density estimate (typically only used for +debugging). To prevent overcounting, the ambient calculation is therefore +disabled in photon flow mode (see ampPmap() in pmapamb.c). + + + +COMPILING + +The files in this package are intended to replace those from the current +RADIANCE CVS release; these are confined to the ray/src/rt subdirectory, +where the photon mapping code resides alongside RADIANCE's raytracing code. + +The photon flow code is built by defining PMAP_PHOTONFLOW during compilation +of the RADIANCE suite. This may be accelerated by only (re)compiling in +ray/src/rt if a RADIANCE build already exists: + + rmake OPT+=-DPMAP_PHOTONFLOW install + +The photon flow code contains debugging code that can be output for each +photon path lookup table operation. In addition, a built-in sanity check +of the lookup table will be run after each such operation. These options can +be enabled by defining PMAP_DEBUGPATHS in addition to PMAP_PHOTONFLOW: + + rmake OPT+="-DPMAP_PHOTONFLOW -DPMAP_DEBUGPATHS" install + +Note that these options are independent of the underlying data structure used +to store the photon map. The default options build the legacy kd-tree. If +a larger, out-of-core photon map is desired (to accommodate >1 billion photons, +for example), this build option must be enabled by defining PMAP_OOC: + + rmake OPT+="-DPMAP_OOC -DPMAP_PHOTONFLOW" install + +For details, see the RADIANCE Out-Of-Core Photon Map Technical Report at +http://dx.doi.org/10.13140/RG.2.1.2158.9363 + + + +USAGE + +A volume photon map representing the physical lightfield is generated with +mkpmap by defining a global participating medium on the command line that +does not absorb (albedo = 1.0) nor scatter (scattering eccentricity = 1), +and thus does not alter a photon's trajectory along its path as it traverses +the medium: + + mkpmap -apv photonFlow.vpm 1M [-apo port1 ... -apo portN] -apD 0.001 \ + -me .1 .1 .1 -ma 1 1 1 -mg 1 scene.oct + +The extinction -me is a free parameter here and determines the photon density +per path; this should be higher than typical for mist to capture sufficient +detail of the physical lightfield; a value of .1 is already a very dense +medium, but can be increased if more detail is desired. + +There is a caveat with using high values of -me, namely that the two-pass +photon distribution algorithm used by mkpmap to achieve the target number of +photons fails. The fraction of photons emitted in the 1st pass (prepass) is +specified with -apD (predistribution factor). Its default value of 0.25 is +tuned for "well behaved" geometry and materials using global or caustic +photons. This value is too high for volume photons in a dense participating +medium, and the photon distribution overflows in the prepass with a "photon +map overflow" error, meaning that the target photon count was already exceeded (technically this does not yield an incorrect solution, but it's not what the +user wanted). When using photon flow, this value should be significantly +reduced from its default to prevent an overflow in the prepass; in the example +above, an -apD of 0.001 accommodates the size of the photon map that results +with a relatively high -me of 0.1. + +Once the volume photon map has been generated, the physical lightfield it +represents can be evaluated using rtrace. Photon flow is enabled using a new +(experimental) option -apF, which disables the normal ambient calculation and +performs the new volume photon density estimate at the given sensor +position(s). Because the ambient calculation is entirely circumvented, the +-ab parameter is now irrelevant; the volume photon density is evaluated +immediately at each sensor position: + + rtrace -h -I -ap photonFlow.vpm 250 -apF -am 100 scene.oct < sensors.dat \ + rcalc -e '$1 = 179 * $1' + +In this example, the irradiance from photon flow is evaluated with a bandwidth +of 250 photons (from unique paths as described above) around each sensor +position, and converted to illuminance using rcalc with a standard luminous +efficacy for daylight of 179 lm/W. + +Here too, there is a caveat regarding the maximum search radius used by the +nearest neighbour search routines to locate photons. This radius is +automatically initialised from the average photon density to balance +performance against conservative pruning, and is optimised for "well behaved" +global and caustic photon maps. However, this fails with photon flow due to +the additional pruning action of the photon path filter, resulting in far +fewer photons for the density estimate than requested, which in turn increases +noise in the illuminance. In this case, rtrace issues a "short lookup in +photon flow" warning, and recommends using the -am parameter. + +The -am parameter manually sets the maximum search radius to fine-tune the +lookups; if set, the bandwidth becomes the _maximum_ number of photons that +will be sought within the given radius. This can be used, for example, to +limit bias by setting the radius to the sensor grid spacing (or a multiple +thereof, to allow for some overlap). Even if set, rtrace may still issue +occasional warnings recommending a higher radius; this is left to the user's +discretion, as increasing the radius also increases bias, and may indicate +that the size of the photon map overall may be too small for the given +sensor grid spacing. + +Finally, the photon flow distribution may be visualised by either dumping it +as a RADIANCE scene description, and rendering the photons as spheres, e.g. + + pmapdump -n 1m photonFlow.vpm | objview + +or by dumping it as a point list for use in a point cloud viewer: + + pmapdump -a -f -N -P -n 1m photonFlow.vpm > photonFlow.xyz + +In the latter example, pmapdump also outputs each photon's radiant flux +(in W), its direction (or normal for surface-bound global and caustic photons), +and its path ID for subsequent filtering. See the pmapdump man page or the +RADIANCE photon map manual at http://dx.doi.org/10.13140/RG.2.1.4330.8405/1 +for details. + + + +CONCLUSION + +This code is experimental and needs thorough testing; I welcome your feedback +in terms of the results but also its usage. The integration of photon flow +into the RADIANCE suite is currently very hacky and could be improved. In +particular, there may be a more elegant way to make this functionality +accessible than using an obscure parameter -apF. Once these issues have been +sorted out, we should consider whether the code be made public and uploaded +to the RADIANCE CVS repository. + +Many thanks in advance for your cooperation and feedback. Happy photon flowing! + + +--Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) + diff --git a/Rmakefile b/Rmakefile index 6091e3e..8438059 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,435 +1,463 @@ # RCSid: $Id: Rmakefile,v 2.85 2018/11/08 00:54:07 greg Exp $ # # Compiles for ray tracing programs. # OPT = -O MACH = -DBSD CFLAGS = -I../common -L../lib $(OPT) $(MACH) SPECIAL = CC = cc AR = ar MLIB = -lm LINT = lint LINTFLAGS = -DBSD # # The following are user-definable: # DESTDIR = . INSTDIR = /usr/local/bin INSTALL = cp # # The following paths must exist and be relative to root: # DEVDIR = $(INSTDIR)/dev LIBDIR = /usr/local/lib/ray # # Library routines: # RLIB = ../lib/libradiance.a RCLIB = ../lib/libraycalls.a LIBS = -lrtrad $(MLIB) # # Device drivers for rvu (see also devtable.c): # DOBJS = devtable.o devcomm.o editline.o x11.o x11twind.o \ colortab.o DSRC = devtable.c devcomm.c editline.c x11.c x11twind.c \ colortab.c DLIBS = -lX11 # # Standard object files: # RTOBJS = rtmain.o rtrace.o duphead.o persist.o RTSRC = rtmain.c rtrace.c duphead.c persist.c RPOBJS = rpmain.o rpict.o srcdraw.o duphead.o persist.o RPSRC = rpmain.c rpict.c srcdraw.c duphead.c persist.c RVOBJS = rvmain.o rview.o rv2.o rv3.o $(DOBJS) RVSRC = rvmain.c rview.c rv2.c rv3.c $(DSRC) RCOBJS = rcmain.o rcontrib.o rc2.o rc3.o RCSRC = rcmain.c rcontrib.c rc2.c rc3.c RLOBJS = raycalls.o raypcalls.o rayfifo.o RLSRC = raycalls.c raypcalls.c rayfifo.c ROBJS = $(RAYOBJS) $(SURFOBJS) $(MATOBJS) \ $(MODOBJS) $(SUPPOBJS) $(PMOBJS) RSRC = $(RAYSRC) $(SURFSRC) $(MATSRC) \ $(MODSRC) $(SUPPSRC) RAYOBJS = ambcomp.o ambient.o ambio.o freeobjmem.o initotypes.o \ preload.o raytrace.o renderopts.o RAYSRC = ambcomp.c ambient.c ambio.c freeobjmem.c initotypes.c \ preload.c raytrace.c renderopts.c SURFOBJS = source.o sphere.o srcobstr.o srcsupp.o srcsamp.o virtuals.o \ o_face.o o_cone.o o_instance.o o_mesh.o SURFSRC = sphere.c source.c srcobstr.c srcsupp.c srcsamp.c virtuals.c \ o_face.c srcsamp.c o_cone.c o_instance.c o_mesh.c MATOBJS = aniso.o normal.o dielectric.o m_clip.o glass.o m_brdf.o \ m_mirror.o m_direct.o m_mist.o fprism.o m_alias.o m_bsdf.o \ ashikhmin.o MATSRC = aniso.c normal.c dielectric.c m_clip.c glass.c m_brdf.c \ m_mirror.c m_direct.c m_mist.c fprism.c m_alias.c m_bsdf.c \ ashikhmin.c MODOBJS = p_func.o t_func.o p_data.o t_data.o text.o mx_func.o mx_data.o MODSRC = p_func.c t_func.c p_data.c t_data.c text.c mx_func.c mx_data.c SUPPOBJS = func.o noise3.o data.o SUPPSRC = func.c noise3.c data.c -PMOBJS = pmap.o pmapsrc.o pmapmat.o pmaprand.o pmapio.o pmapdata.o \ +PMOBJS = pmap.o pmapsrc.o pmapmat.o pmaprand.o pmapio.o pmapdata.o pmapdens.o \ pmapbias.o pmapparm.o pmapcontrib.o pmapamb.o pmapray.o pmapopt.o \ pmapdiag.o pmaptype.o oocmorton.o oococt.o oocsort.o oocbuild.o \ - oocnn.o ooccache.o pmutil.o pmcontrib2.o -PMSRC = pmap.c pmapsrc.c pmapmat.c pmaprand.c pmapio.c pmapdata.c \ + oocnn.o ooccache.o pmutil.o pmcontrib2.o pmapfilt.o +PMSRC = pmap.c pmapsrc.c pmapmat.c pmaprand.c pmapio.c pmapdata.c pmapdens.c \ pmapbias.c pmapparm.c pmapcontrib.c pmapamb.c pmapray.c pmapopt.c \ pmapdiag.c pmaptype.c pmapkdt.c pmapooc.c oocmorton.c oococt.c \ - oocsort.c oocbuild.c oocnn.c ooccache.c pmutil.c pmcontrib2.c + oocsort.c oocbuild.c oocnn.c ooccache.c pmutil.c pmcontrib2.c \ + pmapfilt.c HEADERS = ambient.h ray.h data.h otspecial.h source.h # # What this makefile produces: # PROGS = $(DESTDIR)/rtrace $(DESTDIR)/rpict $(DESTDIR)/rvu $(DESTDIR)/rcontrib \ $(DESTDIR)/lookamb $(DESTDIR)/mkpmap $(DESTDIR)/pmapdump all: $(PROGS) $(RCLIB) $(SPECIAL) install: all rayinit.cal $(INSTALL) $(PROGS) $(INSTDIR) cd $(INSTDIR) ; rm -f rview ; ln -s rvu rview cd $(INSTDIR) ; rm -f rtcontrib ; ln -s rcontrib rtcontrib rm -f $(LIBDIR)/rayinit.cal cp rayinit.cal $(LIBDIR) ogl: clean: set nonomatch; rm -f $(PROGS) *.o lint: $(RVSRC) $(LINT) $(LINTFLAGS) -DRVIEW $(RVSRC) $(LIBS) # # Links: # $(DESTDIR)/rtrace: $(RTOBJS) $(RCLIB) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rtrace $(RTOBJS) $(RCLIB) $(RLIB) $(LIBS) $(DESTDIR)/rpict: $(RPOBJS) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rpict $(RPOBJS) $(RLIB) $(LIBS) $(DESTDIR)/rvu: $(RVOBJS) $(RCLIB) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rvu $(RVOBJS) $(RCLIB) \ $(RLIB) $(LIBS) $(DLIBS) $(DESTDIR)/rcontrib: $(RCOBJS) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rcontrib $(RCOBJS) $(RLIB) $(LIBS) $(DESTDIR)/lookamb: lookamb.o ambio.o $(CC) $(CFLAGS) -o $(DESTDIR)/lookamb lookamb.o ambio.o $(LIBS) $(DESTDIR)/mkpmap: mkpmap.o $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/mkpmap mkpmap.o $(RLIB) $(LIBS) $(DESTDIR)/pmapdump: pmapdump.o pmaptype.o pmapparm.o $(CC) $(CFLAGS) -o pmapdump pmapdump.o pmaptype.o pmapparm.o $(LIBS) $(RLIB): $(ROBJS) Version.o rm -f $(RLIB) $(AR) rc $(RLIB) $(ROBJS) Version.o -ranlib $(RLIB) $(RCLIB): $(RLOBJS) rm -f $(RCLIB) $(AR) rc $(RCLIB) $(RLOBJS) -ranlib $(RCLIB) # # Uncomment the following to model dispersion: # dielectric.o: dielectric.c source.h $(CC) $(CFLAGS) -DDISPERSE -c dielectric.c # end of dispersion compiles. devcomm.o: devcomm.c $(CC) $(CFLAGS) -DDEVPATH=\"$(DEVDIR)\" -c devcomm.c # # Version module: # Version.c: VERSION $(RSRC) $(HEADERS) ( cat VERSION ; date ; whoami ; hostname ) > Version.c ed - Version.c < verscript.ed # # Include dependencies: # ambio.o colortab.o data.o devcomm.o \ devmain.o lookamb.o rview.o x11.o: ../common/color.h freeobjmem.o o_cone.o srcsupp.o: ../common/cone.h data.o freeobjmem.o m_brdf.o mx_data.o \ p_data.o raycalls.o t_data.o: data.h devcomm.o devmain.o devtable.o \ editline.o x11.o: driver.h freeobjmem.o o_face.o srcsupp.o: ../common/face.h ambient.o raytrace.o rpmain.o rtmain.o \ rtrace.o rvmain.o rv2.o rv3.o: ../common/octree.h o_instance.o: ../common/instance.h ambient.o aniso.o ashikhmin.o dielectric.o freeobjmem.o func.o glass.o \ initotypes.o m_brdf.o m_bsdf.o m_direct.o m_mirror.o normal.o o_cone.o \ preload.o raycalls.o raytrace.o rtrace.o rv2.o source.o sphere.o \ srcsupp.o text.o srcdraw.o srcobstr.o virtuals.o: ../common/otypes.h ambient.o ambcomp.o aniso.o ashikhmin.o normal.o raycalls.o raytrace.o \ rpict.o rvmain.o rtmain.o rpmain.o rcmain.o persist.o source.o rv3.o \ srcsamp.o virtuals.o: ../common/random.h ambcomp.o ambient.o aniso.o ashikhmin.o dielectric.o freeobjmem.o func.o \ glass.o m_bsdf.o m_brdf.o m_clip.o m_direct.o m_mirror.o m_mist.o mx_data.o \ o_mesh.o mx_func.o normal.o o_cone.o o_face.o o_instance.o p_data.o p_func.o \ raycalls.o raypcalls.o rayfifo.o raytrace.o rpict.o rtrace.o rv2.o rv3.o rview.o \ source.o sphere.o srcdraw.o srcobstr.o srcsamp.o srcsupp.o t_data.o t_func.o \ text.o rpmain.o rtmain.o rvmain.o virtuals.o m_alias.o rcmain.o \ rcontrib.o rc2.o rc3.o: ray.h \ ../common/standard.h ../common/rtmisc.h ../common/rtio.h ../common/rtmath.h \ ../common/rterror.h ../common/octree.h \ ../common/mat4.h ../common/fvect.h ../common/object.h ../common/color.h rv2.o rv3.o rview.o: rpaint.h driver.h ../common/view.h ../common/resolu.h m_direct.o m_mirror.o m_mist.o dielectric.o raycalls.o \ rpict.o rpmain.o rtmain.o rvmain.o source.o srcdraw.o \ srcobstr.o srcsamp.o srcsupp.o virtuals.o: source.h cone.o data.o devcomm.o initotypes.o fprism.o preload.o \ duphead.o octree.o: ../common/standard.h ../common/rtmisc.h \ ../common/rtio.h ../common/rtmath.h \ ../common/rterror.h ../common/mat4.h ../common/fvect.h ambio.o: ../common/fvect.h ambient.o initotypes.o m_alias.o pmapcontrib.o pmapkdt.o \ pmapooc.o pmapsrc.o pmcontrib2.o \ raytrace.o rtrace.o rv2.o rv3.o source.o \ srcdraw.o srcobstr.o virtuals.o: otspecial.h rpmain.o rtmain.o rvmain.o rpict.o \ srcdraw.o: ../common/view.h ../common/resolu.h rpict.o: ../common/hilbert.h x11.o x11twind.o: x11twind.h x11.o: x11icon.h ambient.o ambcomp.o ambio.o lookamb.o raycalls.o: ambient.h data.o rpmain.o rtmain.o rvmain.o rpict.o rtrace.o \ rv2.o: ../common/resolu.h aniso.o func.o m_brdf.o m_direct.o mx_data.o mx_func.o p_data.o \ p_func.o t_data.o t_func.o: func.h ../common/calcomp.h preload.o: data.h func.h ../common/object.h ../common/face.h \ ../common/cone.h ../common/instance.h ../common/mesh.h \ ../common/color.h ../common/bsdf.h ../common/otypes.h rtmain.o rpmain.o rvmain.o persist.o duphead.o \ renderopts.o rpict.o: ../common/paths.h freeobjmem.o raycalls.o text.o: ../common/font.h raypcalls.o: ../common/selcall.h o_mesh.o: ../common/mesh.h noise3.o: ../common/calcomp.h aniso.o ashikhmin.o dielectric.o freeobjmem.o glass.o initotypes.o \ m_alias.o m_brdf.o m_clip.o m_direct.o m_mirror.o m_mist.o \ mx_data.o mx_func.o normal.o o_cone.o o_face.o o_instance.o \ o_mesh.o p_data.o p_func.o source.o sphere.o t_data.o t_func.o \ srcobstr.o text.o: rtotypes.h m_bsdf.o: ambient.h source.h func.h \ ../common/calcomp.h ../common/bsdf.h ../common/random.h rcmain.o rcontrib.o rc2.o rc3.o: rcontrib.h \ ../common/platform.h ../common/paths.h ../common/lookup.h \ func.h ../common/calcomp.h ../common/rtprocess.h ambient.o rcmain.o: ambient.h rcmain.o: source.h rcontrib.o: source.h ../common/otypes.h rc2.o: ../common/resolu.h rc3.o: ../common/selcall.h # # Photon map include dependencies (via 'gcc -MM -I../common') # ambient.o: pmapparm.h pmaptype.h pmapamb.h pmapdata.h dielectric.o glass.o normal.o m_brdf.o m_bsdf.o ashikhmin.o aniso.o: \ pmapparm.h pmaptype.h pmapmat.h pmap.h pmapdata.h raycalls.o rpmain.o rcmain.o rtmain.o rvmain.o: \ pmapparm.h pmaptype.h pmapray.h rcmain.o: pmapparm.h pmaptype.h pmapray.h pmapcontrib.h pmapdata.h raytrace.o: pmapparm.h pmaptype.h pmap.h pmapdata.h renderopts.o: pmapparm.h pmaptype.h pmapopt.h rpict.o: pmapparm.h pmaptype.h pmapbias.h pmapdata.h pmapdiag.h source.o: pmapparm.h pmaptype.h pmap.h pmapdata.h pmapsrc.h -pmapamb.o: pmapamb.c pmapamb.h pmapdata.h ray.h ../common/standard.h \ - ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ - ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ - ../common/rterror.h ../common/octree.h ../common/object.h \ - ../common/color.h pmapparm.h pmaptype.h ../common/lookup.h pmap.h - pmapbias.o: pmapbias.c pmapbias.h pmapdata.h ray.h ../common/standard.h \ ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ ../common/rterror.h ../common/octree.h ../common/object.h \ ../common/color.h pmapparm.h pmaptype.h ../common/lookup.h pmap.h \ pmaprand.h -pmap.o: pmap.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ - ../common/standard.h ../common/copyright.h ../common/rtio.h \ - ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ - ../common/fvect.h ../common/rterror.h ../common/octree.h \ - ../common/object.h ../common/color.h ../common/paths.h \ - ../common/lookup.h pmapkdt.h pmapooc.h pmapmat.h pmapsrc.h source.h pmaprand.h \ - pmapio.h pmapbias.h pmapdiag.h ../common/platform.h ../common/otypes.h - pmapcontrib.o: pmapcontrib.c pmapcontrib.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ ../common/paths.h ../common/lookup.h pmapkdt.h pmapooc.h pmapmat.h pmap.h \ pmapsrc.h source.h pmaprand.h pmapio.h pmapdiag.h ../common/platform.h \ rcontrib.h ../common/rtprocess.h ../common/paths.h func.h \ ../common/calcomp.h ../common/otypes.h -pmapdata.o: pmapdata.c pmapdata.h ray.h ../common/standard.h \ - ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ - ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ - ../common/rterror.h ../common/octree.h ../common/object.h \ - ../common/color.h pmapparm.h pmaptype.h ../common/paths.h \ - ../common/lookup.h pmapkdt.h pmapooc.h pmaprand.h pmapmat.h pmap.h \ - ../common/otypes.h source.h rcontrib.h ../common/platform.h \ - ../common/rtprocess.h ../common/paths.h func.h ../common/calcomp.h \ - ../common/random.h pmapkdt.c pmapooc.c - pmapdiag.o: pmapdiag.c pmapdiag.h ../common/platform.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ ../common/lookup.h -pmapdump.o: pmapdump.c pmapio.h pmapdata.h ray.h ../common/standard.h \ - ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ - ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ - ../common/rterror.h ../common/octree.h ../common/object.h \ - ../common/color.h pmapparm.h pmaptype.h ../common/lookup.h \ - ../common/rtio.h ../common/resolu.h - pmapio.o: pmapio.c pmapio.h pmapdata.h ray.h ../common/standard.h \ ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ ../common/rterror.h ../common/octree.h ../common/object.h \ ../common/color.h pmapparm.h pmaptype.h ../common/paths.h \ ../common/lookup.h pmapkdt.h pmapooc.h pmapdiag.h ../common/platform.h \ ../common/resolu.h -pmapmat.o: pmapmat.c pmapmat.h pmap.h pmapparm.h pmaptype.h pmapdata.h \ +pmapray.o: pmapray.c pmapray.h ray.h ../common/standard.h \ + ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ + ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ + ../common/rterror.h ../common/octree.h ../common/object.h \ + ../common/color.h pmapparm.h pmaptype.h pmap.h pmapdata.h \ + ../common/lookup.h + +pmaptype.o: pmaptype.c pmaptype.h + +ooccache.o: ooccache.c ooccache.h + +oocsort.o: oocsort.c oocsort.h ../common/fvect.h oocmorton.h + +oocbuild.o: oocbuild.c oococt.h oocmorton.h ../common/fvect.h ooccache.h \ + oocsort.h + +oococt.o: oococt.c oococt.h oocmorton.h ../common/fvect.h ooccache.h \ + ../common/rtio.h + +oocnn.o: oocnn.c oocnn.h oococt.h oocmorton.h ../common/fvect.h \ + ooccache.h pmapfilt.h ../common/lookup.h oocsort.h + +pmapfilt.o: pmapfilt.c pmapfilt.h ../common/lookup.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ + ../common/paths.h pmapkdt.h ../common/random.h source.h otspecial.h + +pmapdens.o: pmapdens.c pmapdens.h pmap.h pmapparm.h pmaptype.h pmapdata.h \ ray.h ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ - ../common/object.h ../common/color.h ../common/lookup.h pmaprand.h \ - ../common/otypes.h data.h func.h ../common/calcomp.h ../common/bsdf.h \ - ../common/ccolor.h + ../common/object.h ../common/color.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmapbias.h ../common/otypes.h -pmapopt.o: pmapopt.c pmapparm.h pmaptype.h ../common/rtio.h \ - ../common/rterror.h +mkpmap.o: mkpmap.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h ../common/paths.h \ + ../common/lookup.h pmapmat.h pmapsrc.h source.h pmapcontrib.h \ + pmaprand.h ambient.h ../common/resolu.h pmapkdt.c pmapkdt.h \ + pmapooc.h pmapooc.c -pmapparm.o: pmapparm.c pmapparm.h pmaptype.h pmapdata.h ray.h \ +pmapdump.o: pmapdump.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h ../common/paths.h \ - ../common/lookup.h pmapkdt.h pmapooc.h + ../common/lookup.h pmapkdt.h pmapio.h ../common/rtio.h \ + ../common/resolu.h ../common/random.h -pmapray.o: pmapray.c pmapray.h ray.h ../common/standard.h \ +pmapdata.o: pmapdata.c pmapdata.h ray.h ../common/standard.h \ ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ ../common/rterror.h ../common/octree.h ../common/object.h \ - ../common/color.h pmapparm.h pmaptype.h pmap.h pmapdata.h \ - ../common/lookup.h + ../common/color.h pmapparm.h pmaptype.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmapdens.h pmap.h pmaprand.h pmapmat.h \ + ../common/otypes.h ../common/random.h pmapkdt.c source.h otspecial.h + +pmap.o: pmap.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmapmat.h pmapsrc.h source.h pmaprand.h \ + pmapio.h pmapdens.h pmapbias.h pmapdiag.h ../common/platform.h \ + ../common/otypes.h + +pmapamb.o: pmapamb.c pmapamb.h pmapdata.h ray.h ../common/standard.h \ + ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ + ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ + ../common/rterror.h ../common/octree.h ../common/object.h \ + ../common/color.h pmapparm.h pmaptype.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmap.h + +pmapmat.o: pmapmat.c pmapmat.h pmap.h pmapparm.h pmaptype.h pmapdata.h \ + ray.h ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmaprand.h ../common/otypes.h data.h func.h \ + ../common/calcomp.h ../common/bsdf.h ../common/ccolor.h \ + ../common/platform.h pmapsrc.o: pmapsrc.c pmapsrc.h ray.h ../common/standard.h \ ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ ../common/rterror.h ../common/octree.h ../common/object.h \ ../common/color.h pmapparm.h pmaptype.h source.h pmap.h pmapdata.h \ - ../common/lookup.h pmaprand.h ../common/otypes.h + ../common/paths.h ../common/lookup.h pmapkdt.h pmaprand.h \ + ../common/otypes.h otspecial.h -pmutil.o: pmap.h pmapio.h pmapdata.h pmapparm.h pmaptype.h pmapbias.h \ - ../common/otypes.h +pmapopt.o: pmapopt.c ray.h ../common/standard.h ../common/copyright.h \ + ../common/rtio.h ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h \ + ../common/mat4.h ../common/fvect.h ../common/rterror.h \ + ../common/octree.h ../common/object.h ../common/color.h pmapparm.h \ + pmaptype.h -pmaptype.o: pmaptype.c pmaptype.h - -mkpmap.o: mkpmap.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ +pmapparm.o: pmapparm.c pmapparm.h pmaptype.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h ../common/paths.h \ - ../common/lookup.h pmapkdt.h pmapooc.h pmapmat.h pmapcontrib.h pmaprand.h \ - ambient.h ../common/resolu.h source.h - -oococt.o: oococt.c oococt.h oocsort.h ooccache.h \ - ../common/rtio.h ../common/fvect.h - -ooccache.o: ooccache.c ooccache.h - -oocsort.o: oocsort.c oocsort.h ../common/fvect.h - -oocbuild.o: oocbuild.c oocbuild.h oocsort.h oococt.h - -oocnn.o: oocnn.c oocnn.h oococt.h oocsort.h + ../common/lookup.h pmapkdt.h + +pmutil.o: pmutil.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmapio.h diff --git a/mkpmap.c b/mkpmap.c index b32ed1f..70cb64b 100644 --- a/mkpmap.c +++ b/mkpmap.c @@ -1,649 +1,651 @@ #ifndef lint static const char RCSid[] = "$Id: mkpmap.c,v 2.10 2020/08/07 01:21:13 rschregle Exp $"; #endif /* ====================================================================== Photon map generator Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation (DFG) - under the FARESYS project. + 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). + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, - supported by the JSPS KAKENHI Grant Number JP19KK0115. + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: mkpmap.c,v 2.10 2020/08/07 01:21:13 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmapcontrib.h" #include "pmaprand.h" #include "paths.h" #include "ambient.h" #include "resolu.h" #include "source.h" #include #include /* Enable options for Ze Ekspertz only! */ #define PMAP_EKSPERTZ extern char VersionID []; char* progname; /* argv[0] */ int dimlist [MAXDIM]; /* sampling dimensions */ int ndims = 0; /* number of sampling dimenshunns */ char* octname = NULL; /* octree name */ CUBE thescene; /* scene top-level octree */ OBJECT nsceneobjs; /* number of objects in scene */ double srcsizerat = 0.01; /* source partition size ratio */ int backvis = 1; /* back face visibility */ int clobber = 0; /* overwrite output */ COLOR cextinction = BLKCOLOR; /* global extinction coefficient */ COLOR salbedo = BLKCOLOR; /* global scattering albedo */ double seccg = 0; /* global scattering eccentricity */ char *amblist [AMBLLEN + 1]; /* ambient include/exclude list */ int ambincl = -1; /* include == 1, exclude == 0 */ char *diagFile = NULL; /* diagnostics output file */ int rand_samp = 1; /* uncorrelated random sampling */ unsigned nproc = 1; /* number of parallel processes */ #ifdef EVALDRC_HACK char *angsrcfile = NULL; /* angular source file for EvalDRC */ #endif /* Dummies for linkage */ COLOR ambval = BLKCOLOR; double shadthresh = .05, ambacc = 0.2, shadcert = .5, minweight = 5e-3, ssampdist = 0, dstrsrc = 0.0, specthresh = 0.15, specjitter = 1.0, avgrefl = 0.5; int ambvwt = 0, ambssamp = 0, ambres = 32, ambounce = 0, directrelay = 1, directvis = 1, samplendx, do_irrad = 0, ambdiv = 128, vspretest = 512, maxdepth = 6, contrib = 0; char *shm_boundary = NULL, *ambfile = NULL, *RCCONTEXT = NULL; void (*trace)() = NULL, (*addobjnotify [])() = {ambnotify, NULL}; void printdefaults() /* print default values to stdout */ { #ifdef EVALDRC_HACK /* EvalDRC support */ puts("-A\t\t\t\t# angular source file"); #endif puts("-ae mod\t\t\t\t# exclude modifier"); puts("-aE file\t\t\t\t# exclude modifiers from file"); puts("-ai mod\t\t\t\t# include modifier"); puts("-aI file\t\t\t\t# include modifiers from file"); #ifdef PMAP_EKSPERTZ puts("-api xmin ymin zmin xmax ymax zmax\t# region of interest"); #endif puts("-apg file nPhotons\t\t\t# global photon map"); puts("-apc file nPhotons\t\t\t# caustic photon map"); puts("-apd file nPhotons\t\t\t# direct photon map"); puts("-app file nPhotons bwidth\t\t# precomputed global photon map"); puts("-apv file nPhotons\t\t\t# volume photon map"); puts("-apC file nPhotons\t\t\t# contribution photon map"); printf("-apD %f\t\t\t\t# predistribution factor\n", preDistrib); printf("-apM %d\t\t\t\t\t# max predistrib passes\n", maxPreDistrib); #if 1 /* Kept for backwards compat, will be gradually phased out by -ld, -lr */ printf("-apm %ld\t\t\t\t# limit photon bounces\n", photonMaxBounce); #endif puts("-apo+ mod\t\t\t\t# photon port modifier"); puts("-apO+ file\t\t\t\t# photon ports from file"); printf("-apP %f\t\t\t\t# precomputation factor\n", finalGather); printf("-apr %d\t\t\t\t\t# random seed\n", randSeed); puts("-aps mod\t\t\t\t# antimatter sensor modifier"); puts("-apS file\t\t\t\t# antimatter sensors from file"); printf(backvis ? "-bv+\t\t\t\t\t# back face visibility on\n" : "-bv-\t\t\t\t\t# back face visibility off\n"); printf("-dp %.1f\t\t\t\t# PDF samples / sr\n", pdfSamples); printf("-ds %f\t\t\t\t# source partition size ratio\n", srcsizerat); printf("-e %s\t\t\t\t# diagnostics output file\n", diagFile); printf(clobber ? "-fo+\t\t\t\t\t# force overwrite\n" : "-fo-\t\t\t\t\t# do not overwrite\n"); #ifdef PMAP_EKSPERTZ /* (Formerly) NU STUFF for Ze Exspertz! */ printf("-ld %.1f\t\t\t\t\t# limit photon distance\n", photonMaxDist); printf("-lr %ld\t\t\t\t# limit photon bounces\n", photonMaxBounce); #endif printf("-ma %.2f %.2f %.2f\t\t\t# scattering albedo\n", colval(salbedo,RED), colval(salbedo,GRN), colval(salbedo,BLU)); printf("-me %.2e %.2e %.2e\t\t# extinction coefficient\n", colval(cextinction,RED), colval(cextinction,GRN), colval(cextinction,BLU)); printf("-mg %.2f\t\t\t\t# scattering eccentricity\n", seccg); #if NIX /* Multiprocessing on NIX only; so tuff luck, Windoze Weenies! */ printf("-n %d\t\t\t\t\t# number of parallel processes\n", nproc); #endif printf("-t %-9d\t\t\t\t# time between reports\n", photonRepTime); printf(verbose ? "-v+\t\t\t\t\t# verbose console output\n" : "-v-\t\t\t\t\t# terse console output\n"); } int main (int argc, char* argv []) { #define check(ol, al) if ( \ argv [i][ol] || badarg(argc - i - 1,argv + i + 1, al) \ ) goto badopt /* Evaluate boolean option, setting var accordingly */ #define check_bool(olen, var) switch (argv [i][olen]) { \ case '\0': \ var = !var; break; \ case 'y': case 'Y': case 't': case 'T': case '+': case '1': \ var = 1; break; \ case 'n': case 'N': case 'f': case 'F': case '-': case '0': \ var = 0; break; \ default: \ goto badopt; \ } /* Evaluate trinary option, setting bits v1 and v2 in var accordingly */ #define check_tri(olen, v1, v2, var) switch (argv [i][olen]) { \ case '\0': case '+': \ var = v1; break; \ case '-': \ var = v2; break;\ case '0': \ var = v1 | v2; break; \ default: \ goto badopt; \ } int loadflags = IO_CHECK | IO_SCENE | IO_TREE | IO_BOUNDS, rval, i, j, n; char **portLp = photonPortList, **sensLp = photonSensorList, **amblp = NULL, sbuf [MAXSTR], portFlags [2] = "\0\0"; struct stat pmstat; /* Global program name */ progname = fixargv0(argv [0]); /* Initialize object types */ initotypes(); /* Parse options */ for (i = 1; i < argc; i++) { /* Eggs-pand arguments */ while ((rval = expandarg(&argc, &argv, i))) if (rval < 0) { sprintf(errmsg, "cannot eggs-pand '%s'", argv [i]); error(SYSTEM, errmsg); } if (argv[i] == NULL) break; if (!strcmp(argv [i], "-version")) { puts(VersionID); quit(0); } if (!strcmp(argv [i], "-defaults") || !strcmp(argv [i], "-help")) { printdefaults(); quit(0); } /* Get octree */ if (i == argc - 1) { octname = argv [i]; break; } switch (argv [i][1]) { case 'a': /* Ambient */ switch (argv [i][2]) { case 'i': /* Ambient include */ case 'I': check(3, "s"); if (ambincl != 1) { ambincl = 1; amblp = amblist; } if (argv [i][2] == 'I') { /* Add modifiers from file */ rval = wordfile(amblp, AMBLLEN - (amblp - amblist), getpath(argv [++i], getrlibpath(), R_OK)); if (rval < 0) { sprintf(errmsg, "cannot open ambient include file \"%s\"", argv [i]); error(SYSTEM, errmsg); } amblp += rval; } else { /* Add modifier from next arg */ *amblp++ = savqstr(argv [++i]); *amblp = NULL; } break; case 'e': /* Ambient exclude */ case 'E': check(3, "s"); if (ambincl != 0) { ambincl = 0; amblp = amblist; } if (argv [i][2] == 'E') { /* Add modifiers from file */ rval = wordfile(amblp, AMBLLEN - (amblp - amblist), getpath(argv [++i], getrlibpath(), R_OK)); if (rval < 0) { sprintf(errmsg, "cannot open ambient exclude file \"%s\"", argv [i]); error(SYSTEM, errmsg); } amblp += rval; } else { /* Add modifier from next arg */ *amblp++ = savqstr(argv [++i]); *amblp = NULL; } break; case 'p': /* Pmap-specific */ switch (argv [i][3]) { case 'g': /* Global photon map */ check(4, "ss"); globalPmapParams.fileName = argv [++i]; globalPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!globalPmapParams.distribTarget) goto badopt; globalPmapParams.minGather = globalPmapParams.maxGather = 0; break; case 'p': /* Precomputed global photon map */ check(4, "ssi"); preCompPmapParams.fileName = argv [++i]; preCompPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!preCompPmapParams.distribTarget) goto badopt; preCompPmapParams.minGather = preCompPmapParams.maxGather = atoi(argv [++i]); if (!preCompPmapParams.maxGather) goto badopt; break; case 'c': /* Caustic photon map */ check(4, "ss"); causticPmapParams.fileName = argv [++i]; causticPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!causticPmapParams.distribTarget) goto badopt; break; case 'v': /* Volume photon map */ check(4, "ss"); volumePmapParams.fileName = argv [++i]; volumePmapParams.distribTarget = parseMultiplier(argv [++i]); if (!volumePmapParams.distribTarget) goto badopt; break; case 'd': /* Direct photon map */ check(4, "ss"); directPmapParams.fileName = argv [++i]; directPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!directPmapParams.distribTarget) goto badopt; break; case 'C': /* Contribution photon map */ check(4, "ss"); contribPmapParams.fileName = argv [++i]; contribPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!contribPmapParams.distribTarget) goto badopt; break; case 'D': /* Predistribution factor */ check(4, "f"); preDistrib = atof(argv [++i]); if (preDistrib <= 0) error(USER, "predistrib factor must be > 0"); break; case 'M': /* Max predistribution passes */ check(4, "i"); maxPreDistrib = atoi(argv [++i]); if (maxPreDistrib <= 0) error(USER, "max predistrib passes must be > 0"); break; #if 1 /* Kept for backwards compat, to be phased out by -lr */ case 'm': /* Max photon bounces */ check(4, "i"); photonMaxBounce = atol(argv [++i]); if (photonMaxBounce <= 0) error(USER, "max photon bounces must be > 0"); break; #endif #ifdef PMAP_EKSPERTZ case 'i': /* Add region of interest */ check(4, "ffffff"); n = pmapNumROI; pmapROI = realloc(pmapROI, ++pmapNumROI * sizeof(PhotonMapROI)); if (!pmapROI) error(SYSTEM, "failed to allocate ROI"); pmapROI [n].min [0] = atof(argv [++i]); pmapROI [n].min [1] = atof(argv [++i]); pmapROI [n].min [2] = atof(argv [++i]); pmapROI [n].max [0] = atof(argv [++i]); pmapROI [n].max [1] = atof(argv [++i]); pmapROI [n].max [2] = atof(argv [++i]); for (j = 0; j < 3; j++) if (pmapROI [n].min [j] >= pmapROI [n].max [j]) error(USER, "invalid region of interest " "(swapped min/max?)"); break; #endif case 'P': /* Global photon precomp ratio */ check(4, "f"); finalGather = atof(argv [++i]); if (finalGather <= 0 || finalGather > 1) error(USER, "global photon precomputation ratio " "must be in range ]0, 1]"); break; case 'o': /* Photon port */ case 'O': /* Check for bad arg and length, taking into account * default forward orientation if none specified, in * order to maintain previous behaviour */ check(argv [i][4] ? 5 : 4, "s"); /* Get port orientation flags */ check_tri( 4, PMAP_PORTFWD, PMAP_PORTBWD, portFlags [0] ); if (argv [i][3] == 'O') { /* Add port modifiers from file */ rval = wordfile( portLp, MAXSET - (portLp - photonPortList), getpath(argv [++i], getrlibpath(), R_OK) ); if (rval < 0) { sprintf( errmsg, "cannot open photon port file %s", argv [i] ); error(SYSTEM, errmsg); } /* HACK: append port orientation flags to every * modifier; note this requires reallocation */ for (; rval--; portLp++) { j = strlen(*portLp); if (!(*portLp = realloc(*portLp, j + 2))) { sprintf( errmsg, "cannot allocate photon port modifiers" " from file %s", argv [i] ); error(SYSTEM, errmsg); } strcat(*portLp, portFlags); } } else { /* Append port flags to port modifier, add to * port list and mark of end list with NULL */ strcpy(sbuf, argv [++i]); strcat(sbuf, portFlags); *portLp++ = savqstr(sbuf); *portLp = NULL; } break; case 'r': /* Random seed */ check(4, "i"); randSeed = atoi(argv [++i]); break; case 's': /* Antimatter sensor */ case 'S': check(4, "s"); if (argv[i][3] == 'S') { /* Add sensor modifiers from file */ rval = wordfile(sensLp, MAXSET - (sensLp - photonSensorList), getpath(argv [++i], getrlibpath(), R_OK)); if (rval < 0) { sprintf(errmsg, "cannot open antimatter sensor file %s", argv [i]); error(SYSTEM, errmsg); } sensLp += rval; } else { /* Append modifier to sensor list, mark end with * NULL */ *sensLp++ = savqstr(argv [++i]); *sensLp = NULL; } break; default: goto badopt; } break; default: goto badopt; } break; case 'b': /* Back face visibility */ if (argv [i][2] == 'v') { check_bool(3, backvis); } else goto badopt; break; case 'd': /* Direct */ switch (argv [i][2]) { case 'p': /* PDF samples */ check(3, "f"); pdfSamples = atof(argv [++i]); break; case 's': /* Source partition size ratio */ check(3, "f"); srcsizerat = atof(argv [++i]); break; default: goto badopt; } break; case 'e': /* Diagnostics file */ check(2, "s"); diagFile = argv [++i]; break; case 'f': /* Force overwrite */ if (argv [i][2] == 'o') { check_bool(3, clobber); } else goto badopt; break; #ifdef PMAP_EKSPERTZ case 'l': /* Limits */ switch (argv [i][2]) { case 'd': /* Limit photon path distance */ check(3, "f"); photonMaxDist = atof(argv [++i]); if (photonMaxDist <= 0) error(USER, "max photon distance must be > 0"); break; case 'r': /* Limit photon bounces */ check(3, "i"); photonMaxBounce = atol(argv [++i]); if (photonMaxBounce <= 0) error(USER, "max photon bounces must be > 0"); break; default: goto badopt; } break; #endif case 'm': /* Medium */ switch (argv[i][2]) { case 'e': /* Eggs-tinction coefficient */ check(3, "fff"); setcolor(cextinction, atof(argv [i + 1]), atof(argv [i + 2]), atof(argv [i + 3])); i += 3; break; case 'a': /* Albedo */ check(3, "fff"); setcolor(salbedo, atof(argv [i + 1]), atof(argv [i + 2]), atof(argv [i + 3])); i += 3; break; case 'g': /* Scattering eccentricity */ check(3, "f"); seccg = atof(argv [++i]); break; default: goto badopt; } break; #if NIX case 'n': /* Num parallel processes (NIX only) */ check(2, "i"); nproc = atoi(argv [++i]); if (nproc > PMAP_MAXPROC) { nproc = PMAP_MAXPROC; sprintf(errmsg, "too many parallel processes, clamping to " "%d\n", nproc); error(WARNING, errmsg); } break; #endif case 't': /* Timer */ check(2, "i"); photonRepTime = atoi(argv [++i]); break; case 'v': /* Verbosity */ check_bool(2, verbose); break; #ifdef EVALDRC_HACK case 'A': /* Angular source file */ check(2,"s"); angsrcfile = argv[++i]; break; #endif default: goto badopt; } } /* Open diagnostics file */ if (diagFile) { if (!freopen(diagFile, "a", stderr)) quit(2); fprintf(stderr, "**************\n*** PID %5d: ", getpid()); printargs(argc, argv, stderr); putc('\n', stderr); fflush(stderr); } #ifdef NICE /* Lower priority */ nice(NICE); #endif if (octname == NULL) error(USER, "missing octree argument"); /* Allocate photon maps and set parameters */ for (i = 0; i < NUM_PMAP_TYPES; i++) { setPmapParam(photonMaps + i, pmapParams + i); /* Don't overwrite existing photon map unless clobbering enabled */ if (photonMaps [i] && !stat(photonMaps [i] -> fileName, &pmstat) && !clobber) { sprintf(errmsg, "photon map file %s exists, not overwritten", photonMaps [i] -> fileName); error(USER, errmsg); } } for (i = 0; i < NUM_PMAP_TYPES && !photonMaps [i]; i++); if (i >= NUM_PMAP_TYPES) error(USER, "no photon maps specified"); readoct(octname, loadflags, &thescene, NULL); #ifdef EVALDRC_HACK if (angsrcfile) readobj(angsrcfile); /* load angular sources */ #endif nsceneobjs = nobjects; /* Get sources */ marksources(); /* Do forward pass and build photon maps */ if (contribPmap) /* Just build contrib pmap, ignore others */ distribPhotonContrib(contribPmap, nproc); else distribPhotons(photonMaps, nproc); /* Save photon maps; no idea why GCC needs an explicit cast here... */ savePmaps((const PhotonMap**)photonMaps, argc, argv); cleanUpPmaps(photonMaps); quit(0); badopt: sprintf(errmsg, "command line error at '%s'", argv[i]); error(USER, errmsg); #undef check #undef check_bool return 0; } diff --git a/oocnn.c b/oocnn.c index e6388fb..31c8e25 100644 --- a/oocnn.c +++ b/oocnn.c @@ -1,328 +1,476 @@ #ifndef lint -static const char RCSid[] = "$Id: oocnn.c,v 2.2 2017/08/14 21:12:10 rschregle Exp $"; +static const char RCSid[] = "$Id: oocnn.c,v 1.2 2020/11/10 01:09:47 u-no-hoo Exp u-no-hoo $"; #endif /* ========================================================================= k-nearest neighbour lookup routines for out-of-core octree data structure Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation (SNSF, #147053) + 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: oocnn.c,v 2.2 2017/08/14 21:12:10 rschregle Exp $ + $Id: oocnn.c,v 1.2 2020/11/10 01:09:47 u-no-hoo Exp u-no-hoo $ */ -#if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) +#if !defined(_WIN32) && !defined(_WIN64) && defined(PMAP_OOC) /* No Windoze support for now */ #include "oocnn.h" #include "oocsort.h" #include #include #include +/* ========================= SEARCH QUEUE STUFF ========================= */ + #ifdef DEBUG_OOC_NN -static int OOC_SearchQueueCheck (OOC_SearchQueue *queue, const FVECT key, - RREAL *(*keyFunc)(const void*), - unsigned root) -/* Priority queue sanity check */ -{ - unsigned kid; - const OOC_SearchQueueNode *qn = queue -> node; - void *rec; - float d2; - - if (root < queue -> tail) { - rec = OOC_GetNearest(queue, qn [root].idx); - d2 = dist2(keyFunc(rec), key); + static int OOC_CheckSearchQueue ( + OOC_SearchQueue *squeue, const FVECT key, + RREAL *(*keyFunc)(const void*), unsigned root + ) + /* Priority queue sanity check */ + { + unsigned kid; + const OOC_SearchQueueNode *sqn = squeue -> node; + void *rec; + float d2; - if (qn [root].dist2 != d2) - return -1; - - if ((kid = (root << 1) + 1) < queue -> tail) { - if (qn [kid].dist2 > qn [root].dist2) - return -1; - else return OOC_SearchQueueCheck(queue, key, keyFunc, kid); - } + if (root < squeue -> tail) { + rec = OOC_GetNearest(squeue, sqn [root].idx); + d2 = dist2(keyFunc(rec), key); - if (++kid < queue -> tail) { - if (qn [kid].dist2 > qn [root].dist2) + if (sqn [root].dist2 != d2) return -1; - else return OOC_SearchQueueCheck(queue, key, keyFunc, kid); + + if ((kid = (root << 1) + 1) < squeue -> tail) { + if (sqn [kid].dist2 > sqn [root].dist2) + return -1; + else return OOC_CheckSearchQueue(squeue, key, keyFunc, kid); + } + + if (++kid < squeue -> tail) { + if (sqn [kid].dist2 > sqn [root].dist2) + return -1; + else return OOC_CheckSearchQueue(squeue, key, keyFunc, kid); + } } + + return 0; } - - return 0; -} #endif -static float OOC_PutNearest (OOC_SearchQueue *queue, float d2, void *rec) +static float OOC_PutNearest ( + OOC_SearchQueue *squeue, float d2, void *rec, + const OOC_SearchAttribCallback *attrib +) /* Insert data record with SQUARED distance to query point into search * priority queue, maintaining the most distant record at the queue head. * If the queue is full, the new record is only inserted if it is closer * than the queue head. - * Returns the new maximum SQUARED distance at the head if the queue is - * full. Otherwise returns -1, indicating a maximum for the entire queue is - * as yet undefined - * The data record is copied into the queue's local record buffa for - * post-search retrieval to minimise redundant disk access. Note that it - * suffices to only rearrange the corresponding indices in the queue nodes - * when restoring the priority queue after every insertion, rather than - * moving the actual records. */ + * + * 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 record attributes to their queue nodes + * by calling attrib->findFunc() if attrib != NULL. It enables finding a + * previously queued data record and replacing it by rec if the latter's + * SQUARED distance d2 is lower. + * + * The record is copied into the queue's local record buffa for post-search + * retrieval to minimise redundant disk access. Note that it suffices to + * only rearrange the corresponding indices in the queue nodes when + * restoring the priority queue after every insertion, rather than moving + * the actual records. + * + * 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. + */ { - OOC_SearchQueueNode *qn = queue -> node; + void *tRec = NULL; + OOC_SearchQueueNode *sqn = squeue -> node, + **attribSqn = NULL, **tAttribSqn = NULL; unsigned root, kid, kid2, rootIdx; - float d2max = -1; /* Undefined max distance ^2 */ + /* Start with undefined max distance */ + float d2max = -1; - /* The queue is represented as a linearised binary tree with the root - * corresponding to the queue head, and the tail corresponding to the - * last leaf */ - if (queue -> tail < queue -> len) { - /* Enlarge queue if not full, insert at tail and resort towards head */ - kid = queue -> tail++; +#ifdef PMAP_PHOTONFLOW + /* Find previously enqueued record with same attribute if callback + defined */ + #ifdef PMAP_DEBUGPATHS + printf("-------- Enqueue --------\n"); + #endif + attribSqn = (OOC_SearchQueueNode**)( + attrib ? attrib -> findFunc(rec, 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 record at tail and resort towards head */ + kid = squeue -> tail++; while (kid) { root = (kid - 1) >> 1; - /* Compare with parent and swap if necessary, else terminate */ - if (d2 > qn [root].dist2) { - qn [kid].dist2 = qn [root].dist2; - qn [kid].idx = qn [root].idx; + /* 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_PHOTONFLOW + if (attrib) { + /* Update root's attribute map entry to refer to kid */ + tRec = OOC_GetNearest(squeue, sqn [root].idx); + tAttribSqn = (OOC_SearchQueueNode**)( + attrib -> findFunc(tRec, attrib -> state) + ); + #ifdef PMAP_DEBUGPATHS + printf("update\n"); + #endif + if (!tAttribSqn) { + /* Consistency error: a previously enqueued record must + * have an entry in the attribute map! */ + fputs("OOC_PutNearest: corrupt attribute map\n", stderr); + exit(-1); + } + + *tAttribSqn = &sqn [kid]; + } +#endif kid = root; } else break; } - /* Assign tail position as linear index into record buffa - * queue -> nnRec and append record */ - qn [kid].dist2 = d2; - qn [kid].idx = queue -> tail - 1; - memcpy(OOC_GetNearest(queue, qn [kid].idx), rec, queue -> recSize); + /* Priority queue restored; assign tail position as linear index into + * record buffa squeue -> nnRec and append record */ + sqn [kid].dist2 = d2; + sqn [kid].idx = squeue -> tail - 1; + memcpy(OOC_GetNearest(squeue, sqn [kid].idx), rec, squeue -> recSize); + +#ifdef PMAP_PHOTONFLOW + if (attrib) + /* Add new record to attribute map */ + attrib -> addFunc(rec, &sqn [kid], attrib -> state); +#endif } - else if (d2 < qn [0].dist2) { - /* Queue full and new record closer than maximum at head; replace head - * and resort towards tail */ - root = 0; - rootIdx = qn [root].idx; + else { + /* Search queue full or duplicate attribute in queue and new node may + be closer; set root to replaced node, otherwise to queue head */ + root = attribSqn ? *attribSqn - sqn : 0; - while ((kid = (root << 1) + 1) < queue -> tail) { - /* Compare with larger kid & swap if necessary, else terminate */ - if ((kid2 = (kid + 1)) < queue -> tail && - qn [kid2].dist2 > qn [kid].dist2) - kid = kid2; - - if (d2 < qn [kid].dist2) { - qn [root].dist2 = qn [kid].dist2; - qn [root].idx = qn [kid].idx; + if (d2 < sqn [root].dist2) { + /* New record closer than maximum at root; replace and resort + towards leaves (=search queue tail) */ + rootIdx = sqn [root].idx; + +#ifdef PMAP_PHOTONFLOW + if (attrib && !attribSqn) { + /* New record will replace root and has unique attribute; delete + * root from the attribute map */ + tRec = OOC_GetNearest(squeue, sqn [root].idx); + attrib -> delFunc(tRec, 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_PHOTONFLOW + if (attrib) { + /* Update kid's attribute map entry to refer to root */ + tRec = OOC_GetNearest(squeue, sqn [kid].idx); + tAttribSqn = (OOC_SearchQueueNode**)( + attrib -> findFunc(tRec, attrib -> state) + ); + #ifdef PMAP_DEBUGPATHS + printf("update\n"); + #endif + if (!tAttribSqn) { + /* Consistency error: a previously enqueued record must + * have an entry in the attribute map! */ + fputs("OOC_PutNearest: corrupt attribute map\n", stderr); + exit(-1); + } + + *tAttribSqn = &sqn [root]; + } +#endif + } + else break; + + root = kid; } - else break; - root = kid; + /* Reassign head's / replaced node's previous buffa index and + overwrite corresponding record */ + sqn [root].dist2 = d2; + sqn [root].idx = rootIdx; + memcpy(OOC_GetNearest(squeue, sqn[root].idx), rec, squeue->recSize); + +#ifdef PMAP_PHOTONFLOW + if (attrib) { + if (!attribSqn) + /* Add new record to attribute map */ + attrib -> addFunc(rec, &sqn [root], attrib -> state); + else { + #ifdef PMAP_DEBUGPATHS + attrib -> findFunc(rec, attrib -> state); + printf("update\n"); + #endif + /* Update new record's existing attribute map entry */ + *attribSqn = &sqn [root]; + } + } +#endif } - - /* Reassign head's previous buffa index and overwrite corresponding - * record */ - qn [root].dist2 = d2; - qn [root].idx = rootIdx; - memcpy(OOC_GetNearest(queue, qn [root].idx), rec, queue -> recSize); - - /* Update SQUARED maximum distance from head node */ - d2max = qn [0].dist2; } + if (squeue -> tail >= squeue -> len) + /* Search queue full; update max distance^2 from head node */ + d2max = sqn [0].dist2; + +#if defined(PMAP_PHOTONFLOW) && defined(PMAP_DEBUGPATHS) + if (attrib) + /* Run sanity check on attribute map */ + attrib -> checkFunc(attrib -> state); +#endif + return d2max; } -int OOC_InitNearest (OOC_SearchQueue *squeue, - unsigned len, unsigned recSize) +int OOC_InitNearest (OOC_SearchQueue *squeue, unsigned len, unsigned recSize) { squeue -> len = len; squeue -> recSize = recSize; squeue -> node = calloc(len, sizeof(OOC_SearchQueueNode)); squeue -> nnRec = calloc(len, recSize); if (!squeue -> node || !squeue -> nnRec) { perror("OOC_InitNearest: failed search queue allocation"); return -1; } return 0; } void *OOC_GetNearest (const OOC_SearchQueue *squeue, unsigned idx) { return squeue -> nnRec + idx * squeue -> recSize; } +/* ============================ OCTREE STUFF ============================ */ + static float OOC_BBoxDist2 (const FVECT bbOrg, RREAL bbSiz, const FVECT key) /* Return minimum *SQUARED* distance between key and bounding box defined by * bbOrg and bbSiz; a distance of 0 implies the key lies INSIDE the bbox */ { /* Explicit comparison with bbox corners */ int i; FVECT d; for (i = 0; i < 3; i++) { d [i] = key [i] - bbOrg [i]; d [i] = d [i] < 0 ? -d [i] : d [i] - bbSiz; /* Set to 0 if inside bbox */ if (d [i] < 0) d [i] = 0; } return DOT(d, d); } -float OOC_FindNearest (OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, - const FVECT org, float size, const FVECT key, - const OOC_SearchFilter *filter, - OOC_SearchQueue *queue, float maxDist2) +float OOC_FindNearest ( + OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, + const FVECT org, float size, const FVECT key, + const OOC_SearchFilterCallback *filter, + const OOC_SearchAttribCallback *attrib, + OOC_SearchQueue *squeue, float maxDist2 +) { const float kidSize = size * 0.5; unsigned i, kid, kid0; float d2; char rec [oct -> recSize]; FVECT kidOrg; OOC_DataIdx kidDataIdx, recIdx; OOC_Node *kidNode; /* Start with suboctant closest to key */ for (kid0 = 0, i = 0; i < 3; i++) kid0 |= (key [i] > org [i] + kidSize) << i; for (i = 0; i < 7; i++) { kid = kid0 ^ i; kidNode = node; kidDataIdx = dataIdx + OOC_GetKid(oct, &kidNode, kid); /* Prune empty suboctant */ - if ((!kidNode && !OOC_ISLEAF(node)) || - (OOC_ISLEAF(node) && !node -> leaf.num [kid])) + if ( + (!kidNode && !OOC_ISLEAF(node)) || + (OOC_ISLEAF(node) && !node -> leaf.num [kid]) + ) continue; /* Set up suboctant */ VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, kid, kidSize); /* Prune suboctant if not overlapped by maxDist2 */ if (OOC_BBoxDist2(kidOrg, kidSize, key) > maxDist2) continue; if (kidNode) { /* Internal node; recurse into non-empty suboctant */ - maxDist2 = OOC_FindNearest(oct, kidNode, kidDataIdx, kidOrg, - kidSize, key, filter, queue, maxDist2); + maxDist2 = OOC_FindNearest( + oct, kidNode, kidDataIdx, kidOrg, kidSize, key, + filter, attrib, squeue, maxDist2 + ); if (maxDist2 < 0) /* Bail out on error */ break; } - else if (OOC_ISLEAF(node)) + else if (OOC_ISLEAF(node)) { /* Leaf node; do linear check of all records in suboctant */ - for (recIdx = kidDataIdx; - recIdx < kidDataIdx + node -> leaf.num [kid]; recIdx++) { + for ( + recIdx = kidDataIdx; + recIdx < kidDataIdx + node -> leaf.num [kid]; + recIdx++ + ) { if (OOC_GetData(oct, recIdx, rec)) return -1; - if (!filter || filter -> func(rec, filter -> data)) - /* Insert record in search queue SQUARED dist to key below - * maxDist2 and passes filter */ + if (!filter || filter -> filtFunc(rec, filter -> state)) { + /* Insert record in search queue if it passes filter + * and its SQUARED dist to key is below maxDist2 */ if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) { - if ((d2 = OOC_PutNearest(queue, d2, rec)) >= 0) - /* Update maxDist2 if queue is full */ + if ((d2 = OOC_PutNearest(squeue, d2, rec, attrib)) >= 0) + /* Update maxDist2 if search queue is full */ maxDist2 = d2; #ifdef DEBUG_OOC_NN - if (OOC_SearchQueueCheck(queue, key, oct -> key, 0)) { - fprintf(stderr, "OOC_SearchPush: priority queue " - "inconsistency\n"); - return -1; - } + if (OOC_CheckSearchQueue(squeue, key, oct->key, 0)) { + fprintf( + stderr, "OOC_PutNearest: corrupt search queue\n" + ); + return -1; + } #endif } + } } + } } return maxDist2; } -float OOC_Find1Nearest (OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, - const FVECT org, float size, const FVECT key, - const OOC_SearchFilter *filter, void *nnRec, - float maxDist2) +float OOC_Find1Nearest ( + OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, + const FVECT org, float size, const FVECT key, + const OOC_SearchFilterCallback *filter, void *nnRec, float maxDist2 +) { const float kidSize = size * 0.5; unsigned i, kid, kid0; float d2; char rec [oct -> recSize]; FVECT kidOrg; OOC_DataIdx kidDataIdx, recIdx; OOC_Node *kidNode; /* Start with suboctant closest to key */ for (kid0 = 0, i = 0; i < 3; i++) kid0 |= (key [i] > org [i] + kidSize) << i; for (i = 0; i < 7; i++) { kid = kid0 ^ i; kidNode = node; kidDataIdx = dataIdx + OOC_GetKid(oct, &kidNode, kid); /* Prune empty suboctant */ - if ((!kidNode && !OOC_ISLEAF(node)) || - (OOC_ISLEAF(node) && !node -> leaf.num [kid])) + if ( + (!kidNode && !OOC_ISLEAF(node)) || + (OOC_ISLEAF(node) && !node -> leaf.num [kid]) + ) continue; /* Set up suboctant */ VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, kid, kidSize); /* Prune suboctant if not overlapped by maxDist2 */ if (OOC_BBoxDist2(kidOrg, kidSize, key) > maxDist2) continue; if (kidNode) { /* Internal node; recurse into non-empty suboctant */ - maxDist2 = OOC_Find1Nearest(oct, kidNode, kidDataIdx, kidOrg, - kidSize, key, filter, nnRec, maxDist2); + maxDist2 = OOC_Find1Nearest( + oct, kidNode, kidDataIdx, kidOrg, kidSize, + key, filter, nnRec, maxDist2 + ); if (maxDist2 < 0) /* Bail out on error */ break; } else if (OOC_ISLEAF(node)) /* Leaf node; do linear check of all records in suboctant */ - for (recIdx = kidDataIdx; - recIdx < kidDataIdx + node -> leaf.num [kid]; recIdx++) { + for ( + recIdx = kidDataIdx; + recIdx < kidDataIdx + node -> leaf.num [kid]; + recIdx++ + ) { if (OOC_GetData(oct, recIdx, rec)) return -1; - if (!filter || filter -> func(rec, filter -> data)) + if (!filter || filter -> filtFunc(rec, filter -> state)) /* Update closest record and max SQUARED dist to key if it * passes filter */ if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) { memcpy(nnRec, rec, oct -> recSize); maxDist2 = d2; } } } return maxDist2; } #endif /* NIX / PMAP_OOC */ + diff --git a/oocnn.h b/oocnn.h index 55108c4..04ce3e8 100644 --- a/oocnn.h +++ b/oocnn.h @@ -1,87 +1,120 @@ /* ========================================================================= k-nearest neighbour lookup routines for out-of-core octree data structure Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation (SNSF, #147053) + 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: oocnn.h,v 2.1 2016/05/17 17:39:47 rschregle Exp $ + $Id: oocnn.h,v 1.2 2020/11/10 01:15:09 u-no-hoo Exp u-no-hoo $ */ #ifndef OOC_SEARCH_H #define OOC_SEARCH_H #include "oococt.h" + /* Priority queue node for octree lookups */ typedef struct { - float dist2; /* Record's *SQUARED* distance to query point */ - unsigned idx; /* Record's index in queue buffer */ + float dist2; /* Record's distance^2 to query point */ + unsigned idx; /* Record's index in queue buffer */ } OOC_SearchQueueNode; - - + /* Priority queue for octree lookups. Note that we explicitly store the * NN candidates in a local in-core buffer nnRec for later retrieval - * without incurring additional disk I/O. */ + * without incurring additional disk I/O */ typedef struct { OOC_SearchQueueNode *node; unsigned len, tail, recSize; void *nnRec; } OOC_SearchQueue; - - - /* Filter for k-NN search; records are accepted for which func(rec, data) - * returns nonzero */ - typedef struct { - int (*func)(void *rec, void *data); - void *data; - } OOC_SearchFilter; - - - float OOC_FindNearest (OOC_Octree *oct, OOC_Node *node, - OOC_DataIdx dataIdx, const FVECT org, float size, - const FVECT key, const OOC_SearchFilter *filter, - OOC_SearchQueue *queue, float maxDist2); + /* Requires above typedefs */ + #include "pmapfilt.h" + + + float OOC_FindNearest ( + OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, + const FVECT org, float size, const FVECT key, + const OOC_SearchFilterCallback *filter, + const OOC_SearchAttribCallback *attrib, + OOC_SearchQueue *queue, float maxDist2 + ); /* Find k nearest neighbours (where k = queue -> len) within a maximum * SQUARED distance of maxDist2 around key. Returns the corresponding * data records with their SQUARED distances in the search queue * (queue -> node [0] .. queue -> node [queue -> tail - 1]), with the * furthest neighbour at the queue head and queue->tail <= queue->len. * * This is a recursive function requiring params for the current node: * DataIdx is the data offset for records in the current node, which is * necessary for implied addressing into the leaf file. Org and size are * the origin and size of the current node's bounding box. * - * An optional filter may be specified; if filter != NULL, filter->func() - * is called for each potential nearest neighbour along with additional - * data provided by the caller; see OOC_SearchFilter typedef above. + * This function accepts two optional callbacks for selective search: + * /SHALL MODE ON * - * Return value is the SQUARED distance to furthest neighbour on success, - * else -1 on failure. */ + * -- If filter != NULL, filter->filtFunc() is called for each candidate + * nearest neighbour along with the encapsulated state data provided + * by the caller; records SHALL be accepted if filter->filtFunc() + * returns nonzero. + * + * -- If attrib != NULL, attrib->findFunc() is called for each candidate + * nearest neighbour that has been accepted by filter->filtFunc() and + * will therefore be added to the search queue. + * This function SHALL map a record attribute to its search queue node + * via a lookup table ("attribute map") accessible via attrib->state, + * and SHALL return a modifiable pointer to the corresponding search + * queue node, or NULL if the attribute is absent. + * This callback SHALL support replacing a previously enqueued record + * with the candidate if they share the same attribute and the latter + * is closer to the key; this avoids accumulating records with + * duplicate attributes in the queue, which may be desirable in some + * applications. Or not. Maybe not. + * This callback SHALL also support updating the attribute map entry + * for each resorted record when restoring the priority queue property + * during insertions; the search queue nodes referenced in the + * attribute map are directly modified by dereferencing the returned + * pointers. + * In addition, the callbacks attrib->addFunc() and attrib->delFunc() + * SHALL insert new entries into the attribute map, or delete evicted + * ones once the search queue is full. + * + * /SHALL MODE OFF + * We use the word "SHALL" a lot here 'cos Lars likes it and it sounds + * very imposing and professional and stuff. + * + * This function returns the SQUARED distance to furthest neighbour on + * success, else -1 on failure. */ - float OOC_Find1Nearest (OOC_Octree *oct, OOC_Node *node, - OOC_DataIdx dataIdx, const FVECT org, float size, - const FVECT key, const OOC_SearchFilter *filter, - void *nnRec, float maxDist2); - /* Find single nearest neighbour within max SQUARED distance maxDist2 - * around key and copy corresponding record in nnRec. This is an - * optimised version of OOC_FindNearest() without a search queue */ + float OOC_Find1Nearest ( + OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, + const FVECT org, float size, const FVECT key, + const OOC_SearchFilterCallback *filter, void *nnRec, float maxDist2 + ); + /* Find single nearest neighbour within max SQUARED distance maxDist2 + * around key and copy corresponding record in nnRec. This is an + * optimised version of OOC_FindNearest() without a search queue */ - int OOC_InitNearest (OOC_SearchQueue *squeue, - unsigned len, unsigned recSize); - /* Initialise NN search queue of length len and local buffa for records - * of size recSize; precondition to calling OOC_FindNearest() */ + int OOC_InitNearest ( + OOC_SearchQueue *squeue, unsigned len, unsigned recSize + ); + /* Initialise NN search queue of length len and local buffa for records + * of size recSize; precondition to calling OOC_FindNearest() */ - void *OOC_GetNearest (const OOC_SearchQueue *queue, unsigned idx); - /* Return pointer to record at pos idx in search queue after calling - * OOC_FindNearest() */ + void *OOC_GetNearest (const OOC_SearchQueue *queue, unsigned idx); + /* Return pointer to record at pos idx in search queue after calling + * OOC_FindNearest() */ #endif + diff --git a/pmap.c b/pmap.c index b79e8a7..bdc517e 100644 --- a/pmap.c +++ b/pmap.c @@ -1,795 +1,816 @@ #ifndef lint static const char RCSid[] = "$Id: pmap.c,v 2.17 2018/12/18 22:14:04 rschregle Exp $"; #endif /* ====================================================================== Photon map main module 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) + 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: pmap.c,v 2.17 2018/12/18 22:14:04 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmapio.h" +#include "pmapdens.h" #include "pmapbias.h" #include "pmapdiag.h" #include "otypes.h" #include #if NIX #include #include #include #endif void savePmaps (const PhotonMap **pmaps, int argc, char **argv) { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv); } } static int photonParticipate (RAY *ray) /* Trace photon through participating medium. Returns 1 if passed through, or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ { int i; RREAL xi1, cosTheta, phi, du, dv; const float cext = colorAvg(ray -> cext), albedo = colorAvg(ray -> albedo), gecc = ray -> gecc, gecc2 = sqr(gecc); FVECT u, v; COLOR cvext; /* Mean free distance until interaction with medium */ ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { if (!incube(&thescene, ray -> rop)) { /* Terminate photon if it has leaked from the scene */ #ifdef DEBUG_PMAP fprintf(stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", ray -> rop [0], ray -> rop [1], ray -> rop [2]); #endif return 0; } setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]), exp(-ray -> rmax * ray -> cext [1]), exp(-ray -> rmax * ray -> cext [2])); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); VCOPY(ray -> rorg, ray -> rop); #if 0 if (albedo > FTINY && ray -> rlvl > 0) #else /* Store volume photons unconditionally in mist to also account for direct inscattering from sources */ if (albedo > FTINY) #endif /* Add to volume photon map */ newPhoton(volumePmap, ray); /* Absorbed? */ if (pmapRandom(rouletteState) > albedo) return 0; /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); scalecolor(ray -> rcol, 1 / albedo); /* Scatter photon */ xi1 = pmapRandom(scatterState); cosTheta = ray -> gecc <= FTINY ? 2 * xi1 - 1 : 0.5 / gecc * (1 + gecc2 - sqr((1 - gecc2) / (1 + gecc * (2 * xi1 - 1)))); phi = 2 * PI * pmapRandom(scatterState); du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */ du *= cos(phi); dv *= sin(phi); /* Get axes u & v perpendicular to photon direction */ i = 0; do { v [0] = v [1] = v [2] = 0; v [i++] = 1; fcross(u, v, ray -> rdir); } while (normalize(u) < FTINY); fcross(v, ray -> rdir, u); for (i = 0; i < 3; i++) ray -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * ray -> rdir [i]; ray -> rlvl++; ray -> rmax = -log(pmapRandom(mediumState)) / cext; + /* fprintf(stderr, "%.3f\n", ray -> rmax); */ } /* Passed through medium until intersecting local object */ setcolor(cvext, exp(-ray -> rot * ray -> cext [0]), exp(-ray -> rot * ray -> cext [1]), exp(-ray -> rot * ray -> cext [2])); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); return 1; } void tracePhoton (RAY *ray) /* Follow photon as it bounces around... */ { long mod; OBJREC *mat, *port = NULL; if (!ray -> parent) { /* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for * !!! primary ray from ray -> ro, then reset the latter to NULL so * !!! as not to interfere with localhit() */ port = ray -> ro; ray -> ro = NULL; } if (ray -> rlvl > photonMaxBounce) { #ifdef PMAP_RUNAWAY_WARN error(WARNING, "runaway photon!"); #endif return; } if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray)) return; if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; if (port && ray -> ro != port) { /* !!! PHOTON PORT REJECTION SAMPLING HACK !!! * Terminate photon if emitted from port without intersecting it; * this can happen when the port's partitions extend beyond its * actual geometry, e.g. with polygons. Since the total flux * relayed by the port is based on the (in this case) larger * partition area, it is overestimated; terminating these photons * constitutes rejection sampling and thereby compensates any bias * incurred by the overestimated flux. */ #ifdef PMAP_PORTREJECT_WARN sprintf(errmsg, "photon outside port %s", ray -> ro -> oname); error(WARNING, errmsg); #endif return; } if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) { /* Transfer ray if modifier is VOID or clipped within antimatta */ RAY tray; photonRay(ray, &tray, PMAP_XFER, NULL); tracePhoton(&tray); } else { /* Scatter for modifier material */ mat = objptr(mod); photonScatter [mat -> otype] (mat, ray); } } } static void preComputeGlobal (PhotonMap *pmap) /* Precompute irradiance from global photons for final gathering for a random subset of finalGather * pmap -> numPhotons photons, and builds the photon map, discarding the original photons. */ /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { unsigned long i, numPreComp; unsigned j; PhotonIdx pIdx; Photon photon; RAY ray; PhotonMap nuPmap; repComplete = numPreComp = finalGather * pmap -> numPhotons; if (verbose) { sprintf(errmsg, "\nPrecomputing irradiance for %ld global photons\n", numPreComp); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Copy photon map for precomputed photons */ memcpy(&nuPmap, pmap, sizeof(PhotonMap)); /* Zero counters, init new heap and extents */ nuPmap.numPhotons = 0; initPhotonHeap(&nuPmap); for (j = 0; j < 3; j++) { nuPmap.minPos [j] = FHUGE; nuPmap.maxPos [j] = -FHUGE; } /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { /* Get random photon from stratified distribution in source heap to * avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)((i + pmapRandom(pmap -> randState)) / finalGather); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection at photon position */ VCOPY(ray.rop, photon.pos); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; /* Get density estimate at photon position */ photonDensity(pmap, &ray, ray.rcol); /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray); /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Flush heap */ flushPhotonHeap(&nuPmap); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Trash original pmap, replace with precomputed one */ deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); if (verbose) { eputs("\nRebuilding precomputed photon map\n"); #if NIX fflush(stderr); #endif } /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); } typedef struct { unsigned long numPhotons [NUM_PMAP_TYPES], numEmitted, numComplete; } PhotonCnt; void distribPhotons (PhotonMap **pmaps, unsigned numProc) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned t, srcIdx, proc; double totalFlux = 0; int shmFile, stat, pid; PhotonMap *pm; PhotonCnt *photonCnt; for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++); if (t >= NUM_PMAP_TYPES) error(USER, "no photon maps defined in distribPhotons"); if (!nsources) error(USER, "no light sources in distribPhotons"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(INTERNAL, "can't allocate source partitions in distribPhotons"); /* Initialise all defined photon maps */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { initPhotonMap(pmaps [t], t); /* Open photon heapfile */ initPhotonHeap(pmaps [t]); /* Per-subprocess target count */ pmaps [t] -> distribTarget /= numProc; if (!pmaps [t] -> distribTarget) error(INTERNAL, "no photons to distribute in distribPhotons"); } initPhotonEmissionFuncs(); initPhotonScatterFuncs(); /* Get photon ports from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); #if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0) error(SYSTEM, "failed shared mem init in distribPhotons"); photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0); if (photonCnt == MAP_FAILED) error(SYSTEM, "failed mapping shared memory in distribPhotons"); #else /* Allocate photon counters statically on Windoze */ if (!(photonCnt = malloc(sizeof(PhotonCnt)))) error(SYSTEM, "failed trivial malloc in distribPhotons"); photonCnt -> numEmitted = photonCnt -> numComplete = 0; #endif /* NIX */ if (verbose) { sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); if (photonPorts) { sprintf(errmsg2, " via %d ports", numPhotonPorts); strcat(errmsg, errmsg2); } strcat(errmsg, "\n"); eputs(errmsg); } /* =================================================================== * FLUX INTEGRATION - Get total photon flux from light sources * =================================================================== */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; - do { /* Need at least one iteration if no ports! */ + do { + /* Set next photon port if defined; note this loop iterates once if + * no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose) { sprintf(errmsg, "\tIntegrating flux from source %s ", source [srcIdx].so -> oname); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } - for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; - emap.partitionCnt++) { + for ( + emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; + emap.partitionCnt++ + ) { initPhotonEmission(&emap, pdfSamples); totalFlux += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); } if (totalFlux < FTINY) error(USER, "zero flux from light sources"); /* Record start time for progress reports */ repStartTime = time(NULL); if (verbose) { sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); eputs(errmsg); } /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { #if NIX if (!(pid = fork())) { /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */ #else if (1) { /* No subprocess under Windoze */ #endif /* Local photon counters for this subprocess */ unsigned passCnt = 0, prePassCnt = 0; unsigned long lastNumPhotons [NUM_PMAP_TYPES]; unsigned long localNumEmitted = 0; /* Num photons emitted by this subprocess alone */ /* Seed RNGs from PID for decorellated photon distribution */ pmapSeed(randSeed + proc, partState); pmapSeed(randSeed + (proc + 1) % numProc, emitState); pmapSeed(randSeed + (proc + 2) % numProc, cntState); pmapSeed(randSeed + (proc + 3) % numProc, mediumState); pmapSeed(randSeed + (proc + 4) % numProc, scatterState); pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); #ifdef DEBUG_PMAP /* Output child process PID after random delay to prevent corrupted * console output due to race condition */ usleep(1e6 * pmapRandom(rouletteState)); fprintf(stderr, "Proc %d: PID = %d " "(waiting 10 sec to attach debugger...)\n", proc, getpid()); /* Allow time for debugger to attach to child process */ sleep(10); #endif for (t = 0; t < NUM_PMAP_TYPES; t++) lastNumPhotons [t] = 0; /* ============================================================= * 2-PASS PHOTON DISTRIBUTION * Pass 1 (pre): emit fraction of target photon count * Pass 2 (main): based on outcome of pass 1, estimate remaining * number of photons to emit to approximate target * count * ============================================================= */ do { double numEmit; if (!passCnt) { /* INIT PASS 1 */ /* Skip if no photons contributed after sufficient * iterations; make it clear to user which photon maps are * missing so (s)he can check geometry and materials */ if (++prePassCnt > maxPreDistrib) { sprintf(errmsg, "proc %d: too many prepasses", proc); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { sprintf(errmsg2, ", no %s photons stored", pmapName [t]); strcat(errmsg, errmsg2); } error(USER, errmsg); break; } /* Num to emit is fraction of minimum target count */ numEmit = FHUGE; for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) numEmit = min(pmaps [t] -> distribTarget, numEmit); numEmit *= preDistrib; } else { /* INIT PASS 2 */ /* Based on the outcome of the predistribution we can now * estimate how many more photons we have to emit for each * photon map to meet its respective target count. This * value is clamped to 0 in case the target has already been * exceeded in the pass 1. */ double maxDistribRatio = 0; /* Set the distribution ratio for each map; this indicates * how many photons of each respective type are stored per * emitted photon, and is used as probability for storing a * photon by newPhoton(). Since this biases the photon * density, newPhoton() promotes the flux of stored photons * to compensate. */ for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { pm -> distribRatio = (double)pm -> distribTarget / pm -> numPhotons - 1; /* Check if photon map "overflowed", i.e. exceeded its * target count in the prepass; correcting the photon * flux via the distribution ratio is no longer * possible, as no more photons of this type will be * stored, so notify the user rather than deliver * incorrect results. In future we should handle this * more intelligently by using the photonFlux in each * photon map to individually correct the flux after * distribution. */ if (pm -> distribRatio <= FTINY) { sprintf(errmsg, "%s photon map overflow in " "prepass, reduce -apD", pmapName [t]); error(INTERNAL, errmsg); } maxDistribRatio = max(pm -> distribRatio, maxDistribRatio); } /* Normalise distribution ratios and calculate number of * photons to emit in main pass */ for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) pm -> distribRatio /= maxDistribRatio; if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY) /* No photons left to distribute in main pass */ break; } /* Update shared completion counter for progreport by parent */ photonCnt -> numComplete += numEmit; /* PHOTON DISTRIBUTION LOOP */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; - do { /* Need at least one iteration if no ports! */ + do { + /* Set next photon port if defined; note this loop iterates + * once if no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose && !proc) { /* Output from subproc 0 only to avoid race condition * on console I/O */ if (!passCnt) sprintf(errmsg, "\tPREPASS %d on source %s ", prePassCnt, source [srcIdx].so -> oname); else sprintf(errmsg, "\tMAIN PASS on source %s ", source [srcIdx].so -> oname); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } - for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; - emap.partitionCnt++) { + for ( + emap.partitionCnt = 0; + emap.partitionCnt < emap.numPartitions; + emap.partitionCnt++ + ) { double partNumEmit; unsigned long partEmitCnt; /* Get photon origin within current source partishunn * and build emission map */ photonOrigin [emap.src -> so -> otype] (&emap); initPhotonEmission(&emap, pdfSamples); /* Number of photons to emit from ziss partishunn -- * proportional to flux; photon ray weight and scalar * flux are uniform (latter only varying in RGB). */ partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux; partEmitCnt = (unsigned long)partNumEmit; /* Probabilistically account for fractional photons */ if (pmapRandom(cntState) < partNumEmit - partEmitCnt) partEmitCnt++; /* Update local and shared (global) emission counter */ photonCnt -> numEmitted += partEmitCnt; localNumEmitted += partEmitCnt; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon based on PDF and trace through scene * until absorbed/leaked */ emitPhoton(&emap, &photonRay); #if 1 if (emap.port) /* !!! PHOTON PORT REJECTION SAMPLING HACK: set * !!! photon port as fake hit object for * !!! primary ray to check for intersection in * !!! tracePhoton() */ + /* TODO: Move this into emitPhoton()? */ photonRay.ro = emap.port -> so; #endif + /* Set subprocess ID in photonRay; extends path ID. + Used by newPhoton() as photon attributes. */ + PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); } /* Update shared global photon count for each pmap */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { photonCnt -> numPhotons [t] += pmaps [t] -> numPhotons - lastNumPhotons [t]; lastNumPhotons [t] = pmaps [t] -> numPhotons; } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) { repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; pmapDistribReport(); } #endif } portCnt++; } while (portCnt < numPhotonPorts); } for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { /* Double preDistrib in case a photon map is empty and * redo pass 1 --> possibility of infinite loop for * pathological scenes (e.g. absorbing materials) */ preDistrib *= 2; break; } if (t >= NUM_PMAP_TYPES) /* No empty photon maps found; now do pass 2 */ passCnt++; } while (passCnt < 2); /* Flush heap buffa for every pmap one final time; * avoids potential data corruption! */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { flushPhotonHeap(pmaps [t]); /* Heap file closed automatically on exit fclose(pmaps [t] -> heap); */ #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: total %ld photons\n", proc, pmaps [t] -> numPhotons); eputs(errmsg); #endif } #if NIX /* Terminate subprocess */ exit(0); #endif } else if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotons"); } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) error(USER, "failed photon distribution"); --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; repProgress = repComplete = 0; for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { /* Get global photon count from shmem updated by subprocs */ repProgress += pm -> numPhotons = photonCnt -> numPhotons [t]; repComplete += pm -> distribTarget; } repComplete *= numProc; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* =================================================================== * POST-DISTRIBUTION - Set photon flux and build data struct for photon * storage, etc. * =================================================================== */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); /* Set photon flux */ totalFlux /= photonCnt -> numEmitted; #if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); #else free(photonCnt); #endif if (verbose) eputs("\n"); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { if (verbose) { sprintf(errmsg, "Building %s photon map\n", pmapName [t]); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc); } /* Precompute photon irradiance if necessary */ if (preCompPmap) { if (verbose) eputs("\n"); preComputeGlobal(preCompPmap); } if (verbose) eputs("\n"); } diff --git a/pmap.h b/pmap.h index cda4c26..0089e7e 100644 --- a/pmap.h +++ b/pmap.h @@ -1,91 +1,105 @@ /* RCSid $Id: pmap.h,v 2.9 2018/01/24 19:39:05 rschregle Exp $ */ /* ====================================================================== Photon map main header 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) + 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: pmap.h,v 2.9 2018/01/24 19:39:05 rschregle Exp $ + + TODO: + -- Precomp. contrib photons */ #ifndef PMAP_H #define PMAP_H - #ifndef NIX - #if defined(_WIN32) || defined(_WIN64) - #define NIX 0 - #else - #define NIX 1 - #endif + #ifndef NIX + #if defined(_WIN32) || defined(_WIN64) + #define NIX 0 + #else + #define NIX 1 + #endif #endif #include "pmapparm.h" #include "pmapdata.h" #ifndef min #define min(a, b) ((a) < (b) ? (a) : (b)) #endif #ifndef max #define max(a, b) ((a) > (b) ? (a) : (b)) #endif #define sqr(a) ((a) * (a)) /* Average over colour channels */ #define colorAvg(col) ((col [0] + col [1] + col [2]) / 3) /* Macros to test for enabled photon maps */ #define photonMapping (globalPmap || preCompPmap || \ causticPmap || contribPmap) - #define causticPhotonMapping (causticPmap != NULL) - #define directPhotonMapping (directPmap != NULL) - #define volumePhotonMapping (volumePmap != NULL) + #define causticPhotonMapping (causticPmap) + #ifdef PMAP_PHOTONFLOW + /* Photon flow also accounts for direct component */ + #define directPhotonMapping (directPmap || photonFlow) + #else + #define directPhotonMapping (directPmap) + #endif + #define volumePhotonMapping (volumePmap) /* #define contribPhotonMapping (contribPmap && contribPmap -> srcContrib) */ #define contribPhotonMapping (contribPmap) extern void (*pmapLookup [])(PhotonMap*, RAY*, COLOR); /* Photon map lookup functions per type */ void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *params); /* Load photon and set their respective parameters, checking timestamps * relative to octree for possible staleness */ void savePmaps (const PhotonMap **pmaps, int argc, char **argv); /* Save all defined photon maps with specified command line */ void cleanUpPmaps (PhotonMap **pmaps); /* Trash all photon maps after processing is complete */ void distribPhotons (PhotonMap **pmaps, unsigned numProc); /* Emit photons from light sources and build photon maps for non-NULL * entries in photon map array */ void tracePhoton (RAY*); /* Follow photon as it bounces around the scene. Analogon to * raytrace(). */ void photonDensity (PhotonMap*, RAY*, COLOR irrad); /* Perform surface density estimate from incoming photon flux at ray's intersection point. Returns irradiance from found photons. */ void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad); /* Returns precomputed photon density estimate at ray -> rop. */ void volumePhotonDensity (PhotonMap*, RAY*, COLOR); /* Perform volume density estimate from incoming photon flux at ray's intersection point. Returns irradiance. */ void colorNorm (COLOR); /* Normalise colour channels to average of 1 */ #endif diff --git a/pmapamb.c b/pmapamb.c index 5dbc240..82bb3f6 100644 --- a/pmapamb.c +++ b/pmapamb.c @@ -1,87 +1,102 @@ #ifndef lint static const char RCSid[] = "$Id: pmapamb.c,v 2.8 2016/05/17 17:39:47 rschregle Exp $"; #endif /* ================================================================== Photon map interface to RADIANCE ambient calculation 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) + 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: pmapamb.c,v 2.8 2016/05/17 17:39:47 rschregle Exp $ */ #include "pmapamb.h" #include "pmap.h" int ambPmap (COLOR aval, RAY *r, int rdepth) /* Factor irradiance from photon map into ambient coefficient aval; * return 1 on success, else 0 (with aval unmodified) */ { COLOR rcoef, photonIrrad; /* Handle precedence in case of multiple photon maps: * contrib > precomp > global */ PhotonMap *pmap = contribPhotonMapping ? contribPmap : preCompPmap ? preCompPmap : globalPmap; - + +#if PMAP_PHOTONFLOW + if (photonFlow) + /* Disable ambient component in photon flow mode; direct and indirect + * irradiance is instead evaluated from volume photons in + * multDirectPmap() (see pmapsrc.c). Here we simply return the + * unmodified irradiance and signal success. */ + return 1; +#endif + /* Get photon irradiance either via 1 ambient bounce (final * gather) if ambounce > 0, or directly if ambounce < 0. */ if (pmap && (rdepth || ambounce < 0)) { /* Temporarily factor ambient value into ray coefficient * (required for contribution photon map) */ copycolor(rcoef, r -> rcoef); multcolor(r -> rcoef, aval); /* Get photon irradiance via callback */ pmap -> lookupFlags = 0; (pmap -> lookup)(pmap, r, photonIrrad); /* Factor irradiance into ambient value and restore ray coeficient */ multcolor(aval, photonIrrad); copycolor(r -> rcoef, rcoef); return 1; } return 0; } int ambPmapCaustic (COLOR aval, RAY *r, int rdepth) /* Factor specular-diffuse (caustic) irradiance from photon map into ambient * coeff aval; return 1 if successful, else 0 (with aval set to zero) */ { COLOR rcoef, photonIrrad; /* Handle precedence in case of multiple photon maps: contrib > caustic */ PhotonMap *pmap = contribPhotonMapping ? contribPmap : causticPmap; /* Get caustic photon density estimate at primary rays or when * filling in ambient rays that have no global photon map to use */ if (pmap && (!rdepth || !globalPmap & !contribPmap & !preCompPmap)) { /* Temporarily factor ambient value into ray coefficient * (required for contribution photon map) */ copycolor(rcoef, r -> rcoef); multcolor(r -> rcoef, aval); /* Set caustic flag and get photon irradiance via callback */ pmap -> lookupCaustic = 1; (pmap -> lookup)(pmap, r, photonIrrad); /* Factor irradiance into ambient value and restore ray coeficient */ multcolor(aval, photonIrrad); copycolor(r -> rcoef, rcoef); return 1; } setcolor(aval, 0, 0, 0); return 0; } diff --git a/pmapdata.c b/pmapdata.c index 3a63d5a..32b1715 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,699 +1,736 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $"; #endif /* ========================================================================= - Photon map types and interface to nearest neighbour lookups in underlying - point cloud data structure. + Photon map types and high level interface to nearest neighbour lookups in + underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). 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) + 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: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $ */ #include "pmapdata.h" +#include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" #include "otypes.h" -#include "source.h" -#include "rcontrib.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #endif /* Ambient include/exclude set (from ambient.c) */ #ifndef MAXASET #define MAXASET 4095 #endif extern OBJECT ambset [MAXASET+1]; + /* Callback to print photon attributes acc. to user defined format */ int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); +/* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ + void initPhotonMap (PhotonMap *pmap, PhotonMapType t) /* Init photon map 'n' stuff... */ { if (!pmap) return; pmap -> numPhotons = 0; pmap -> biasCompHist = NULL; pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; pmap -> gatherTolerance = gatherTolerance; pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; pmap -> numDensity = 0; pmap -> distribRatio = 1; pmap -> type = t; pmap -> squeue.node = NULL; pmap -> squeue.len = 0; +#ifdef PMAP_PHOTONFLOW + /* Init path lookup table and its key buffer */ + pmap -> pathLUT = NULL; + pmap -> pathLUTKeys = NULL; + pmap -> numPathLUTKeys = 0; +#endif + /* Init local RNG state */ pmap -> randState [0] = 10243; pmap -> randState [1] = 39829; pmap -> randState [2] = 9433; pmapSeed(randSeed, pmap -> randState); /* Set up type-specific photon lookup callback */ pmap -> lookup = pmapLookup [t]; /* Mark primary photon ray as unused */ pmap -> lastPrimary.srcIdx = -1; pmap -> numPrimary = 0; pmap -> primaries = NULL; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; #ifdef PMAP_OOC OOC_Null(&pmap -> store); #else kdT_Null(&pmap -> store); #endif } void initPhotonHeap (PhotonMap *pmap) { int fdFlags; if (!pmap) error(INTERNAL, "undefined photon map in initPhotonHeap"); if (!pmap -> heap) { /* Open heap file */ mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) error(SYSTEM, "failed opening heap file in initPhotonHeap"); -#ifdef F_SETFL /* XXX is there an alternate needed for Windows? */ +#ifdef F_SETFL /* XXX is there an alternative needed for Windows? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif/* ftruncate(fileno(pmap -> heap), 0); */ } } void flushPhotonHeap (PhotonMap *pmap) { int fd; const unsigned long len = pmap -> heapBufLen * sizeof(Photon); if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined heap in flushPhotonHeap"); */ return; } /* Atomically seek and write block to heap */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ fd = fileno(pmap -> heap); #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon)); eputs(errmsg); #endif /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) error(SYSTEM, "failed append to heap file in flushPhotonHeap"); #if NIX if (fsync(fd)) error(SYSTEM, "failed fsync in flushPhotonHeap"); #endif pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file)) { dup = 1; for (i = 0; i <= 2; i++) { if (p.pos [i] < thescene.cuorg [i] || p.pos [i] > thescene.cuorg [i] + thescene.cusize) { sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2]); error(WARNING, errmsg); } dup &= p.pos [i] == lastp.pos [i]; } if (dup) { sprintf(errmsg, "consecutive duplicate photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2]); error(WARNING, errmsg); } } return 0; } #endif int newPhoton (PhotonMap* pmap, const RAY* ray) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray -> ro -> omod)) return -1; if (pmapNumROI && pmapROI) { unsigned inROI = 0; /* Store photon if within a region of interest (for ze Ecksperts!) */ for (i = 0; !inROI && i < pmapNumROI; i++) inROI = (ray -> rop [0] >= pmapROI [i].min [0] && ray -> rop [0] <= pmapROI [i].max [0] && ray -> rop [1] >= pmapROI [i].min [1] && ray -> rop [1] <= pmapROI [i].max [1] && ray -> rop [2] >= pmapROI [i].min [2] && ray -> rop [2] <= pmapROI [i].max [2]); if (!inROI) return -1; } /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ scalecolor(photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1)); #else scalecolor(photonFlux, 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1)); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set contrib photon's primary ray and subprocess index (the latter * to linearise the primary ray indices after photon distribution is * complete). Also set primary ray's source index, thereby marking it * as used. */ if (isContribPmap(pmap)) { photon.primary = pmap -> numPrimary; photon.proc = PMAP_GETRAYPROC(ray); pmap -> lastPrimary.srcIdx = ray -> rsrc; } - else photon.primary = 0; + /* Set photon's path ID from ray and subprocess */ + else { + photon.primary = ray -> rno; + photon.proc = PMAP_GETRAYPROC(ray); + } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i] : ray -> ron [i]); if (!pmap -> heapBuf) { /* Lazily allocate heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "failed heap buffer allocation in newPhoton"); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Heap buffa full, flush to heap file */ flushPhotonHeap(pmap); pmap -> numPhotons++; /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); return 0; } void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonPrimaryIdx *primaryOfs, unsigned nproc) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); if (!pmap -> numPhotons) error(INTERNAL, "empty photon map in buildPhotonMap"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) error(SYSTEM, "failed to allocate postprocessed photon heap in" "buildPhotonMap"); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) error(SYSTEM, "failed to open postprocessed photon heap in " "buildPhotonMap"); rewind(pmap -> heap); #ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { #ifdef DEBUG_PMAP printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap)); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap); #ifdef DEBUG_PMAP printf("Got %lu\n", pmap -> heapBufLen); #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; else if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (primaryOfs) /* Linearise photon primary index from subprocess index using the * per-subprocess offsets in primaryOfs */ p -> primary += primaryOfs [p -> proc]; /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { scalecolor(flux, photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0]); setPhotonFlux(p, flux); } /* Update average photon flux; need a double here */ addcolor(avgFlux, flux); } /* Write modified photons to new heap */ fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); if (ferror(nuHeap)) error(SYSTEM, "failed postprocessing photon flux in " "buildPhotonMap"); nCheck += pmap -> heapBufLen; } #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } +/* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ + /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; /* Search position is ray -> rorg for volume photons, since we have no intersection point. Normals are ignored -- these are incident directions). */ - /* NOTE: status returned by XXX_FindPhotons() is currently ignored; + /* XXX/NOTE: status returned by XXX_FindPhotons() is currently ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ if (isVolumePmap(pmap)) { #ifdef PMAP_OOC OOC_FindPhotons(pmap, ray -> rorg, NULL); #else kdT_FindPhotons(pmap, ray -> rorg, NULL); #endif } else { #ifdef PMAP_OOC OOC_FindPhotons(pmap, ray -> rop, ray -> ron); #else kdT_FindPhotons(pmap, ray -> rop, ray -> ron); #endif } #ifdef PMAP_LOOKUP_INFO fprintf(stderr, "%d/%d %s photons found within radius %.3f " "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : ""); #endif if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) { /* Short lookup; too few photons found */ if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { /* Ignore short lookups which return fewer than * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there * really are no photons in the vicinity, and increasing the max * search radius therefore won't help */ #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : ""); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { - /* Bail out after warning if maxDist is fixed */ + /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; redo = 0; } } while (redo); } Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon) { /* Init (squared) search radius to avg photon dist to centre of gravity */ float maxDist2_0 = pmap -> CoGdist; int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } +/* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ + void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; + +#ifdef PMAP_PHOTONFLOW + if (pmap -> pathLUT) { + lu_done(pmap -> pathLUT); + free(pmap -> pathLUT); + } + if (pmap -> pathLUTKeys) { + unsigned k; + + for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); + free(pmap -> pathLUTKeys); + } +#endif } + diff --git a/pmapdata.h b/pmapdata.h index a677bad..137b826 100644 --- a/pmapdata.h +++ b/pmapdata.h @@ -1,338 +1,348 @@ /* RCSid $Id: pmapdata.h,v 2.14 2020/04/08 15:14:21 rschregle Exp $ */ /* ========================================================================= Photon map types and interface to nearest neighbour lookups in underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). Defining PMAP_FLOAT_FLUX stores photon flux as floats rather than packed RGBE for greater precision; this may be necessary when the flux differs significantly in individual colour channels, e.g. with highly saturated colours. 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) + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components") ========================================================================= $Id: pmapdata.h,v 2.14 2020/04/08 15:14:21 rschregle Exp $ */ #ifndef PMAPDATA_H #define PMAPDATA_H #ifndef NIX #if defined(_WIN32) || defined(_WIN64) #define NIX 0 #else #define NIX 1 #endif #endif #if (defined(PMAP_OOC) && !NIX) #error "OOC currently only supported on NIX -- tuff luck." #endif #ifdef PMAP_CBDM /* Enable photon primary hitpoints and incident directions (see struct * PhotonPrimary below). Note this will increase the size of photon * primaries 9-fold (10-fold after alignment)!!! */ #define PMAP_PRIMARYPOS #define PMAP_PRIMARYDIR #endif #include "ray.h" #include "pmaptype.h" #include "paths.h" #include "lookup.h" #include /* Primary photon ray for light source contributions */ typedef struct { int16 srcIdx; /* Index of emitting light source */ /* !!! REDUCED FROM 32 BITS !!! */ #ifdef PMAP_PRIMARYDIR int32 dir; /* Encoded ray direction */ #endif #ifdef PMAP_PRIMARYPOS float pos [3]; /* Hit point */ #endif } PhotonPrimary; #define photonSrcIdx(pm, p) ((pm)->primaries[(p)->primary].srcIdx) /* Photon primary ray index type and limit */ typedef uint32 PhotonPrimaryIdx; #define PMAP_MAXPRIMARY UINT32_MAX /* Macros for photon's generating subprocess field */ #ifdef PMAP_OOC #define PMAP_PROCBITS 7 #else #define PMAP_PROCBITS 5 #endif #define PMAP_MAXPROC (1 << PMAP_PROCBITS) #define PMAP_GETRAYPROC(r) ((r) -> crtype >> 8) #define PMAP_SETRAYPROC(r,p) ((r) -> crtype |= p << 8) /* Tolerance for photon normal check during lookups */ #define PMAP_NORM_TOL 0.02 typedef struct { float pos [3]; /* Photon position */ signed char norm [3]; /* Surface normal at pos (incident direction for volume photons) */ union { struct { #ifndef PMAP_OOC unsigned char discr : 2; /* kd-tree discriminator axis */ #endif unsigned char caustic : 1; /* Specularly scattered (=caustic) */ /* Photon's generating subprocess index, used for primary ray * index linearisation when building contrib pmaps; note this is * reduced for kd-tree to accommodate the discriminator field */ unsigned char proc : PMAP_PROCBITS; }; unsigned char flags; }; #ifdef PMAP_FLOAT_FLUX COLOR flux; #else COLR flux; #endif PhotonPrimaryIdx primary; /* Index to primary ray */ } Photon; /* Define PMAP_FLOAT_FLUX to store photon flux as floats instead of * compact RGBE, which was found to improve accuracy in analytical * validation. */ #ifdef PMAP_FLOAT_FLUX #define setPhotonFlux(p,f) copycolor((p) -> flux, f) #define getPhotonFlux(p,f) copycolor(f, (p) -> flux) #else #define setPhotonFlux(p,f) setcolr((p)->flux, (f)[0], (f)[1], (f)[2]) #define getPhotonFlux(p,f) colr_color(f, (p) -> flux) #endif - + /* Bias compensation history node */ typedef struct { COLOR irrad; float weight; } PhotonBiasCompNode; - /* Forward declaration */ struct PhotonMap; - - + + /* Define search queue and underlying data struct types */ #ifdef PMAP_OOC #include "pmapooc.h" #else #include "pmapkdt.h" #endif /* Mean size of heapfile write buffer, in number of photons */ #define PMAP_HEAPBUFSIZE 1e6 /* Mean idle time between heap locking attempts, in usec */ #define PMAP_HEAPBUFSLEEP 2e6 /* Temporary filename for photon heaps */ #define PMAP_TMPFNAME TEMPLATE #define PMAP_TMPFNLEN (TEMPLEN + 1) typedef struct PhotonMap { PhotonMapType type; /* See pmaptype.h */ char *fileName; /* Photon map file */ /* ================================================================ * PRE/POST-BUILD STORAGE * ================================================================ */ FILE *heap; /* Unsorted photon heap prior to construction of store */ char heapFname [sizeof(PMAP_TMPFNAME)]; Photon *heapBuf; /* Write buffer for above */ unsigned long heapBufLen, /* Current & max size of heapBuf */ heapBufSize; PhotonStorage store; /* Photon storage in space subdividing data struct */ /* ================================================================ * PHOTON DISTRIBUTION STUFF * ================================================================ */ unsigned long distribTarget, /* Num stored specified by user */ numPhotons; /* Num actually stored */ float distribRatio; /* Probability of photon storage */ COLOR photonFlux; /* Average photon flux */ unsigned short randState [3]; /* Local RNG state */ /* ================================================================ * PHOTON LOOKUP STUFF * ================================================================ */ union { /* Flags passed to findPhotons() */ char lookupCaustic : 1; char lookupFlags; }; PhotonSearchQueue squeue; /* Search queue for photon lookups */ unsigned minGather, /* Specified min/max photons per */ maxGather; /* density estimate */ /* NOTE: All distances are SQUARED */ float maxDist2, /* Max search radius */ maxDist0, /* Initial value for above */ maxDist2Limit, /* Hard limit for above */ gatherTolerance; /* Fractional deviation from minGather/ maxGather for short lookup */ void (*lookup)(struct PhotonMap*, RAY*, COLOR); /* Callback for type-specific photon - * lookup (usually density estimate) */ - + * lookup (usually density estimate) */ + +#ifdef PMAP_PHOTONFLOW + LUTAB *pathLUT; /* Photon path lookup table to filter + volume photons */ + char **pathLUTKeys; /* Preallocated buffer to store keys + for path lookup table */ + unsigned numPathLUTKeys; /* Num keys currently in key buffer + (= next free entry at tail) */ +#endif /* ================================================================ * BIAS COMPENSATION STUFF * ================================================================ */ PhotonBiasCompNode *biasCompHist; /* Bias compensation history */ /* ================================================================ * STATISTIX * ================================================================ */ unsigned long totalGathered, /* Total photons gathered */ numDensity, /* Num density estimates */ numLookups, /* Counters for short photon lookups */ numShortLookups; unsigned minGathered, /* Min/max photons actually gathered */ maxGathered, /* per density estimate */ shortLookupPct; /* % of short lookups for stats */ float minError, /* Min/max/rms density estimate error */ maxError, rmsError, CoGdist, /* Avg distance to centre of gravity */ maxPos [3], /* Max & min photon positions */ minPos [3]; FVECT CoG; /* Centre of gravity (avg photon pos) */ /* ================================================================ * PHOTON CONTRIB/COEFF STUFF * ================================================================ */ PhotonPrimary *primaries, /* Photon primary array for rendering */ lastPrimary; /* Current primary for photon distrib */ PhotonPrimaryIdx numPrimary; /* Number of primary rays */ LUTAB *srcContrib; /* Lookup table for source contribs */ } PhotonMap; /* Photon maps by type (see PhotonMapType) */ extern PhotonMap *photonMaps []; /* Macros for specific photon map types */ #define globalPmap (photonMaps [PMAP_TYPE_GLOBAL]) #define preCompPmap (photonMaps [PMAP_TYPE_PRECOMP]) #define causticPmap (photonMaps [PMAP_TYPE_CAUSTIC]) #define volumePmap (photonMaps [PMAP_TYPE_VOLUME]) #define directPmap (photonMaps [PMAP_TYPE_DIRECT]) #define contribPmap (photonMaps [PMAP_TYPE_CONTRIB]) /* Photon map type tests */ #define isGlobalPmap(p) ((p) -> type == PMAP_TYPE_GLOBAL) #define isCausticPmap(p) ((p) -> type == PMAP_TYPE_CAUSTIC) #define isContribPmap(p) ((p) -> type == PMAP_TYPE_CONTRIB) #define isVolumePmap(p) ((p) -> type == PMAP_TYPE_VOLUME) void initPhotonMap (PhotonMap *pmap, PhotonMapType t); /* Initialise empty photon map of specified type */ int newPhoton (PhotonMap *pmap, const RAY *ray); /* Create new photon with ray's direction, intersection point, and flux, and append to unsorted photon heap pmap -> heap. The photon is accepted with probability pmap -> distribRatio for global density control; if the photon is rejected, -1 is returned, else 0. The flux is scaled by ray -> rweight and 1 / pmap -> distribRatio. */ void initPhotonHeap (PhotonMap *pmap); /* Open photon heap file */ void flushPhotonHeap (PhotonMap *pmap); /* Flush photon heap buffa pmap -> heapBuf to heap file pmap -> heap; * used by newPhoton() and to finalise heap in distribPhotons(). */ void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonPrimaryIdx *primaryOfs, unsigned nproc); /* Postpress unsorted photon heap pmap -> heap and build underlying data * structure pmap -> store. This is prerequisite to photon lookups with * findPhotons(). */ /* PhotonFlux is the flux per photon averaged over RGB; this is * multiplied with each photon's flux during the postprocess. In the * case of a contribution photon map, this is an array with a separate * flux specific to each light source due to non-uniform photon emission; * Otherwise it is referenced as a scalar value. Flux is not scaled if * photonFlux == NULL. */ /* Photon map construction may be parallelised if nproc > 1, if * supported. The heap is destroyed on return. */ /* PrimaryOfs is an array of index offsets for the primaries in pmap -> * primaries generated by each of the nproc subprocesses during contrib * photon distribution (see distribPhotonContrib()). These offsets are * used to linearise the photon primary indices in the postprocess. This * linearisation is skipped if primaryOfs == NULL. */ void findPhotons (PhotonMap* pmap, const RAY *ray); /* Find pmap -> squeue.len closest photons to ray -> rop with similar normal. For volume photons ray -> rorg is used and the normal is ignored (being the incident direction in this case). Found photons are placed search queue starting with the furthest photon at pmap -> squeue.node, and pmap -> squeue.tail being the number actually found. */ Photon *find1Photon (PhotonMap *pmap, const RAY *ray, Photon *photon); /* Find single closest photon to ray -> rop with similar normal. Return NULL if none found, else the supplied Photon* buffer, indicating that it contains a valid photon. */ void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon); /* Retrieve photon referenced by idx from pmap -> store */ Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx); /* Retrieve photon from NN search queue after calling findPhotons() */ PhotonIdx firstPhoton (const PhotonMap *pmap); /* Index to first photon, to be passed to getPhoton(). Indices to * subsequent photons can be optained via increment operator (++) */ void deletePhotons (PhotonMap*); /* Free dem mammaries... */ #endif diff --git a/pmapdens.c b/pmapdens.c new file mode 100644 index 0000000..99a1aef --- /dev/null +++ b/pmapdens.c @@ -0,0 +1,251 @@ +#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" + + + +/* Photon map lookup functions per type */ +void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = { + photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity, + photonDensity, photonDensity +}; + + + +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 : "", + 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))); + return; +} + + + +#ifdef PMAP_PHOTONFLOW + void volumePhotonSphIrrad (PhotonMap *pmap, RAY *ray, COLOR irrad) + /* Evaluate 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". + * The found photons in the search queue must originate from disparate + * paths for an unbiased solution. This requires enabling the photon + * path ID filter callback in findPhotons() to remove duplicate paths. + */ + { + unsigned i; + float r2; + COLOR flux; + Photon *photon; + const PhotonSearchQueueNode *sqn; + + setcolor(irrad, 0, 0, 0); + + if (!pmap -> maxGather) + return; + + findPhotons(pmap, ray); + + /* Issue warning if few too photons, as the photon path filter almost + * invariably results in short lookups due to the reduced effective + * photon density after dropping duplicates, thereby invalidating the + * automatic search radius adjustment. As a workaround for now, it's + * up to the user to the user to compensate this with a sufficiently + * large fixed search radius. */ + 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 */ + scalecolor(flux, 2 * (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("r = %f\n", sqrt(r2)); +#endif + + /* Divide accumulated flux by sphere surface area and ambient + normalisation factor PI */ + scalecolor(irrad, 1 / (4 * PI * PI * r2)); + return; + } +#endif diff --git a/pmapdens.h b/pmapdens.h new file mode 100644 index 0000000..703978a --- /dev/null +++ b/pmapdens.h @@ -0,0 +1,55 @@ +/* RCSid $Id$ */ + +/* + ====================================================================== + 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$ +*/ + + + +#ifndef PMAPDENS_H + #define PMAPDENS_H + + #include "pmap.h" + + /* Photon map lookup functions per type */ + void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR); + + void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad); + /* Surface-bound photon density estimate for global and caustic photons. + Returns irradiance at ray -> rop. */ + + void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad); + /* Surface-bound photon density estimate for (single) precomputed photon. + Returns precomputed irradiance at ray -> rop. */ + + void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad); + /* Volume photon density estimate using Henyey-Greenstein phase function. + Returns inscattered irradiance at ray -> rop. */ + + #ifdef PMAP_PHOTONFLOW + void volumePhotonSphIrrad (PhotonMap *pmap, RAY *ray, COLOR irrad); + /* Evaluate 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". + /* The found photons in the search queue must originate from disparate + /* paths for an unbiased solution. This requires enabling the photon + /* path ID filter callback in findPhotons() to remove duplicate paths. + */ + #endif +#endif + diff --git a/pmapdump.c b/pmapdump.c index 7dc7565..a7b54f9 100644 --- a/pmapdump.c +++ b/pmapdump.c @@ -1,362 +1,394 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdump.c,v 2.17 2020/08/07 01:21:13 rschregle Exp $"; #endif /* ====================================================================== Dump photon maps as RADIANCE scene description or ASCII point list to stdout Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation (DFG) - under the FARESYS project. + 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). + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, - supported by the JSPS KAKENHI Grant Number JP19KK0115. + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapdump.c,v 2.17 2020/08/07 01:21:13 rschregle Exp $ */ #include "pmap.h" #include "pmapio.h" #include "rtio.h" #include "resolu.h" #include "random.h" #include "math.h" /* Defaults */ /* Sphere radius as fraction of avg. intersphere dist */ /* Relative scale for sphere radius (fudge factor) */ /* Number of spheres */ #define RADCOEFF 0.05 #define RADSCALE 1.0 #define NSPHERES 10000 -/* Format for optional ASCII output as XYZ RGB points */ -#define POINTFMT "%g\t%g\t%g\t%g\t%g\t%g\n" +/* Formats for ASCII output as XYZ RGB points plus optional normals + (incident directions for volume photons) and path IDs */ +#define POINTFMT "%g\t%g\t%g\t%g\t%g\t%g" +#define NORMFMT "\t%f\t%f\t%f" +#define PATHFMT "\t%u:%lu" /* RADIANCE material and object defs for each photon type */ typedef struct { char *mat, *obj; } RadianceDef; const RadianceDef radDefs [] = { { "void glow mat.global\n0\n0\n4 %g %g %g 0\n", "mat.global sphere obj.global\n0\n0\n4 %g %g %g %g\n" }, { "void glow mat.pglobal\n0\n0\n4 %g %g %g 0\n", "mat.pglobal sphere obj.pglobal\n0\n0\n4 %g %g %g %g\n" }, { "void glow mat.caustic\n0\n0\n4 %g %g %g 0\n", "mat.caustic sphere obj.caustic\n0\n0\n4 %g %g %g %g\n" }, { "void glow mat.volume\n0\n0\n4 %g %g %g 0\n", "mat.volume sphere obj.volume\n0\n0\n4 %g %g %g %g\n" }, { "void glow mat.direct\n0\n0\n4 %g %g %g 0\n", "mat.direct sphere obj.direct\n0\n0\n4 %g %g %g %g\n" }, { "void glow mat.contrib\n0\n0\n4 %g %g %g 0\n", "mat.contrib sphere obj.contrib\n0\n0\n4 %g %g %g %g\n" } }; /* Default colour codes are as follows: global = blue precomp global = cyan caustic = red volume = green direct = magenta contrib = yellow */ const COLOR colDefs [] = { {0.25, 0.25, 2}, {0.1, 1, 1}, {1, 0.1, 0.1}, {0.1, 1, 0.1}, {1, 0.1, 1}, {1, 1, 0.1} }; static int setBool(char *str, unsigned pos, unsigned *var) { switch ((str) [pos]) { case '\0': *var = !*var; break; case 'y': case 'Y': case 't': case 'T': case '+': case '1': *var = 1; break; case 'n': case 'N': case 'f': case 'F': case '-': case '0': *var = 0; break; default: return 0; } return 1; } int main (int argc, char** argv) { - char format [MAXFMTLEN]; - RREAL rad, radScale = RADSCALE, extent, dumpRatio; - unsigned arg, j, ptype, dim, fluxCol = 0, points = 0; - long numSpheres = NSPHERES; - COLOR col = {0, 0, 0}; - FILE *pmapFile; - PhotonMap pm; - PhotonPrimary pri; - Photon p; + char format [MAXFMTLEN]; + RREAL rad, radScale = RADSCALE, extent, dumpRatio; + unsigned arg, j, ptype, dim, + fluxCol = 0, points = 0, normals = 0, paths = 0; + long numSpheres = NSPHERES; + COLOR col = {0, 0, 0}; + FILE *pmapFile; + PhotonMap pm; + PhotonPrimary pri; + Photon p; #ifdef PMAP_OOC - char leafFname [1024]; + char leafFname [1024]; #endif if (argc < 2) { puts("Dump photon maps as RADIANCE scene description " - "or ASCII point list\n"); - printf("Usage: %s " - "[-a] [-r radscale1] [-n num1] " - "[-f | -c rcol1 gcol1 bcol1] pmap1 " - "[-a] [-r radscale2] [-n num2] " - "[-f | -c rcol2 gcol2 bcol2] pmap2 " - "...\n", argv [0]); + "or ASCII point list\n" + ); + printf("Usage:\n" + "%s [-a [-N][-P]] [-r radscale1] [-n num1] [-f | -c r1 g1 b1] pmap1\n" + "\t[-a [-N][-P]] [-r radscale2] [-n num2] [-f | -c r2 g2 b2] pmap2 " + " ...\n", + argv [0] + ); return 1; } for (arg = 1; arg < argc; arg++) { /* Parse options */ if (argv [arg][0] == '-') { switch (argv [arg][1]) { case 'a': if (!setBool(argv [arg], 2, &points)) error(USER, "invalid option syntax at -a"); break; + case 'r': if ((radScale = atof(argv [++arg])) <= 0) error(USER, "invalid radius scale"); break; case 'n': if ((numSpheres = parseMultiplier(argv [++arg])) <= 0) error(USER, "invalid number of points/spheres"); break; + case 'N': + if (!setBool(argv [arg], 2, &normals)) + error(USER, "invalid option syntax at -N"); + break; + + case 'P': + if (!setBool(argv [arg], 2, &paths)) + error(USER, "invalid option syntax at -P"); + break; + case 'c': if (fluxCol) error(USER, "-f and -c are mutually exclusive"); if (badarg(argc - arg - 1, &argv [arg + 1], "fff")) error(USER, "invalid RGB colour"); for (j = 0; j < 3; j++) col [j] = atof(argv [++arg]); break; case 'f': if (intens(col) > 0) error(USER, "-f and -c are mutually exclusive"); if (!setBool(argv [arg], 2, &fluxCol)) error(USER, "invalid option syntax at -f"); break; default: sprintf(errmsg, "unknown option %s", argv [arg]); error(USER, errmsg); return -1; } continue; } + /* Check for compatible options */ + if ((normals || paths) && !points) + error(USER, "option -N or -P only valid in conjuction with -a"); + /* Open next photon map file */ if (!(pmapFile = fopen(argv [arg], "rb"))) { sprintf(errmsg, "can't open %s", argv [arg]); error(SYSTEM, errmsg); } /* Get format string */ strcpy(format, PMAP_FORMAT_GLOB); if (checkheader(pmapFile, format, NULL) != 1) { - sprintf(errmsg, "photon map file %s has unknown format %s", - argv [arg], format); + sprintf(errmsg, + "photon map file %s has unknown format %s", argv [arg], format + ); error(USER, errmsg); } /* Identify photon map type from format string */ for (ptype = 0; - ptype < NUM_PMAP_TYPES && strcmp(pmapFormat [ptype], format); - ptype++); + ptype < NUM_PMAP_TYPES && strcmp(pmapFormat [ptype], format); + ptype++ + ); if (!validPmapType(ptype)) { sprintf(errmsg, "file %s contains an unknown photon map type", - argv [arg]); + argv [arg] + ); error(USER, errmsg); } /* Get file format version and check for compatibility */ if (strcmp(getstr(format, pmapFile), PMAP_FILEVER)) error(USER, "incompatible photon map file format"); if (!points) { /* Dump command line as comment */ fputs("# ", stdout); printargs(argc, argv, stdout); fputc('\n', stdout); } /* Set point/sphere colour if independent of photon flux, output RADIANCE material def if required */ if (!fluxCol) { if (intens(col) <= 0) copycolor(col, colDefs [ptype]); if (!points) { printf(radDefs [ptype].mat, col [0], col [1], col [2]); fputc('\n', stdout); } } /* Get number of photons */ pm.numPhotons = getint(sizeof(pm.numPhotons), pmapFile); /* Skip avg photon flux */ for (j = 0; j < 3; j++) getflt(pmapFile); /* Get distribution extent (min & max photon positions) */ for (j = 0; j < 3; j++) { pm.minPos [j] = getflt(pmapFile); pm.maxPos [j] = getflt(pmapFile); } /* Skip centre of gravity, and avg photon dist to it */ for (j = 0; j < 4; j++) getflt(pmapFile); /* Sphere radius based on avg intersphere dist depending on whether the distribution occupies a 1D line (!), a 2D plane, or 3D volume (= sphere distrib density ^-1/d, where d is the dimensionality of the distribution) */ for (j = 0, extent = 1.0, dim = 0; j < 3; j++) { rad = pm.maxPos [j] - pm.minPos [j]; if (rad > FTINY) { dim++; extent *= rad; } } rad = radScale * RADCOEFF * pow(extent / numSpheres, 1./dim); /* Photon dump probability to satisfy target sphere count */ dumpRatio = min(1, (float)numSpheres / pm.numPhotons); - /* Skip primary rays */ + /* Skip primary rays (currently not loaded; -N only dumps index) */ pm.numPrimary = getint(sizeof(pm.numPrimary), pmapFile); while (pm.numPrimary-- > 0) { /* Skip source index & incident dir */ getint(sizeof(pri.srcIdx), pmapFile); #ifdef PMAP_PRIMARYDIR /* Skip primary incident dir */ getint(sizeof(pri.dir), pmapFile); #endif #ifdef PMAP_PRIMARYPOS /* Skip primary hitpoint */ for (j = 0; j < 3; j++) getflt(pmapFile); #endif } #ifdef PMAP_OOC /* Open leaf file with filename derived from pmap, replace pmapFile * (which is currently the node file) */ strncpy(leafFname, argv [arg], sizeof(leafFname) - 1); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - 1); fclose(pmapFile); if (!(pmapFile = fopen(leafFname, "rb"))) { sprintf(errmsg, "cannot open leaf file %s", leafFname); error(SYSTEM, errmsg); } #endif /* Read photons */ while (pm.numPhotons-- > 0) { #ifdef PMAP_OOC /* Get entire photon record from ooC octree leaf file !!! OOC PMAP FILES CURRENTLY DON'T USE PORTABLE I/O !!! */ if (!fread(&p, sizeof(p), 1, pmapFile)) { sprintf(errmsg, "error reading OOC leaf file %s", leafFname); error(SYSTEM, errmsg); } #else /* kd-tree */ /* Get photon position */ for (j = 0; j < 3; j++) p.pos [j] = getflt(pmapFile); - /* Get photon normal (currently not used) */ + /* Get photon normal (dumped with -N) */ for (j = 0; j < 3; j++) p.norm [j] = getint(1, pmapFile); /* Get photon flux */ #ifdef PMAP_FLOAT_FLUX for (j = 0; j < 3; j++) p.flux [j] = getflt(pmapFile); #else for (j = 0; j < 4; j++) p.flux [j] = getint(1, pmapFile); #endif - - - /* Skip primary ray index */ - getint(sizeof(p.primary), pmapFile); + /* Get primary ray / path index (dumped with -P) */ + p.primary = getint(sizeof(p.primary), pmapFile); - /* Skip flags */ - getint(sizeof(p.flags), pmapFile); + /* Get photons flags (includes subprocess prefix for path ID */ + p.flags = getint(sizeof(p.flags), pmapFile); #endif /* Dump photon probabilistically acc. to target sphere count */ if (frandom() <= dumpRatio) { if (fluxCol) { /* Get photon flux */ getPhotonFlux(&p, col); /* Scale by dumpRatio for energy conservation */ scalecolor(col, 1.0 / dumpRatio); } if (!points) { if (fluxCol) { /* Dump material def if variable (depends on flux) */ printf(radDefs [ptype].mat, col [0], col [1], col [2]); fputc('\n', stdout); } printf(radDefs [ptype].obj, p.pos [0], p.pos [1], p.pos [2], rad); fputc('\n', stdout); } - else /* Dump as XYZ RGB point */ + else { + /* Dump as XYZ RGB point */ printf(POINTFMT, p.pos [0], p.pos [1], p.pos [2], col [0], col [1] ,col [2]); + if (normals) + printf(NORMFMT, p.norm [0] / 127., p.norm [1] / 127., + p.norm [2] / 127. + ); + if (paths) + printf(PATHFMT, p.proc, (unsigned long)p.primary); + fputc('\n', stdout); + } } if (ferror(pmapFile) || feof(pmapFile)) { sprintf(errmsg, "error reading %s", argv [arg]); error(USER, errmsg); } } fclose(pmapFile); /* Reset defaults for next dump */ radScale = RADSCALE; numSpheres = NSPHERES; col [0] = col [1] = col [2] = 0; fluxCol = points = 0; } return 0; } diff --git a/pmapfilt.c b/pmapfilt.c new file mode 100644 index 0000000..d77b587 --- /dev/null +++ b/pmapfilt.c @@ -0,0 +1,349 @@ +/* + ====================================================================== + Photon map filter callbacks for nearest neighbour search + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) + (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: pmapfilt.c,v 1.2 2020/11/16 22:28:23 u-no-hoo Exp u-no-hoo $ +*/ + + + +#include "pmapfilt.h" +#include "pmapdata.h" +#include "random.h" +#include "source.h" +#include "otspecial.h" + + + +int filterPhoton (const void *rec, const void *state) +/* Filter callback for photon kNN search */ +{ + const Photon *photon = rec; + const PhotonSearchFilterState *filtState = state; + const PhotonMap *pmap = filtState -> pmap; + + /* Reject photon if normal faces away (ignored for volume photons) with + * tolerance to account for perturbation; note photon normal is coded + * in range [-127,127], hence we factor this in */ + if ( + filtState -> norm && + DOT(filtState->norm, photon->norm) <= PMAP_NORM_TOL * 127 * frandom() + ) + return 0; + + if (isContribPmap(pmap)) { + /* Lookup in contribution photon map; filter according to emitting + * light source if contrib list set, else accept all */ + + if (pmap -> srcContrib) { + OBJREC *srcMod; + const int srcIdx = photonSrcIdx(pmap, photon); + + if (srcIdx < 0 || srcIdx >= nsources) + error(INTERNAL, "invalid light source index in photon map"); + + srcMod = findmaterial(source [srcIdx].so); + + /* Reject photon if contributions from light source which emitted + * it are not sought */ + if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data) + return 0; + } + + /* Reject non-caustic photon if lookup for caustic contribs */ + if (pmap -> lookupCaustic && !photon -> caustic) + return 0; + } + + /* Accept photon */ + return 1; +} + + + +#ifdef PMAP_PHOTONFLOW + /* Length of key strings for photon path LUT (plus terminator) */ + #define PMAP_PATHKEYLEN ((sizeof(PhotonPrimaryIdx) + 1 << 1) + 1) + + void initPhotonPaths (void *state) + /* Lazily allocate and initialise photon path ID lookup table and + * associate key buffer */ + { + PhotonMap *pmap = state; + + if (!pmap -> pathLUT) { + /* Allocate and init photon path ID lookup table of same size as + * search queue */ + LUTAB initLUT = LU_SINIT(NULL, NULL); + + if (!(pmap -> pathLUT = malloc(sizeof(LUTAB)))) + error(SYSTEM, "failed photon path lookup table allocation"); + memcpy(pmap -> pathLUT, &initLUT, sizeof(initLUT)); + if (!lu_init(pmap -> pathLUT, pmap -> squeue.len)) + error(INTERNAL, "failed photon path lookup table init"); + } + + if (!pmap -> pathLUTKeys) { + /* Allocate key buffer of same size as search queue (+1 for + * scratchpad at tail) */ + unsigned k, maxKeys = pmap -> squeue.len + 1; + + pmap -> pathLUTKeys = (char**)calloc(maxKeys, sizeof(char*)); + for (k = 0; k < maxKeys; k++) + if ( + !pmap -> pathLUTKeys || + !(pmap -> pathLUTKeys [k] = (char*)malloc(PMAP_PATHKEYLEN)) + ) + error(SYSTEM, "failed photon path LUT key buffer allocation"); + } + + /* Reset photon path ID key buffer */ + pmap -> numPathLUTKeys = 0; + } + + + + static char *photonPathKey (const Photon *photon, char *key) + { + /* Generate LUT key from photon path ID. Writes transposed hex string + into specified key buffer. */ + PhotonPrimaryIdx pathId1 = photon -> primary; + unsigned char pathId2 = photon -> proc; + char *kp = key; + unsigned i; + + for ( + i = sizeof(PhotonPrimaryIdx) << 1; + i; + *kp++ = (pathId1 & 0x0f) + 'A', pathId1 >>= 4, i-- + ); + for ( + i = sizeof(unsigned char) << 1; + i; + *kp++ = (pathId2 & 0x0f) + 'A', pathId2 >>= 4, i-- + ); + + *kp++ = '\0'; + return key; + } + + + + int checkPhotonPaths (const void *state) + /* Run sanity check on photon path ID lookup table */ + { + const PhotonMap *pmap = state; + const PhotonSearchQueue *squeue = &pmap -> squeue; + const PhotonSearchQueueNode *squeueNode; + const Photon *photon; + const LUENT *lutEntry; + unsigned squeueIdx; + char key [PMAP_PATHKEYLEN]; + + if (!pmap -> pathLUT || !pmap -> pathLUTKeys) + error(INTERNAL, "photon path lookup table not initialised"); + + for ( + squeueIdx = 0, squeueNode = squeue -> node; + squeueIdx < squeue -> tail; + squeueIdx++, squeueNode++ + ) { + photon = getNearestPhoton(squeue, squeueNode -> idx); + photonPathKey(photon, key); + lutEntry = lu_find(pmap -> pathLUT, key); + if (!lutEntry || !lutEntry -> key || strcmp(lutEntry -> key, key)) + error(CONSISTENCY, "invalid key in photon path LUT check"); + + if ( + !lutEntry -> data || + (PhotonSearchQueueNode*)lutEntry -> data != squeueNode + ) + error(CONSISTENCY, "invalid data in photon path LUT check"); + } + + return 0; + } + + + + /* TODO: REPLACE pathLUTKeys SCRATCHPAD WITH LOCAL ARRAYS??? */ + + void **findPhotonPath (const void *rec, const void *state) + /* Photon path callback for k-NN search to map photon path IDs to search + * queue nodes via a lookup table. State points to the associated photon + * map containing the path lookup table. Returns modifiable pointer to data + * field of found LUT entry, which in turn points to the search queue node + * for the corresponding record. If no entry exists for the path ID, NULL + * is returned. */ + { + const Photon *photon = rec; + const PhotonMap *pmap = state; + char *key; + LUENT *lutEntry; + + if (!pmap -> pathLUT || !pmap -> pathLUTKeys) + error(INTERNAL, "photon path lookup table not initialised"); + + /* Generate key for photon's path and write into next free slot in key + buffer (using it as scratchpad) */ + key = photonPathKey( + photon, pmap -> pathLUTKeys [pmap -> numPathLUTKeys] + ); + + #ifdef PMAP_DEBUGPATHS + printf("%s ", key); + #endif + + /* Look up attribute for generated key, set last LUT entry in state */ + if (!(lutEntry = lu_find(pmap -> pathLUT, key))) + /* Outta mammaries */ + error(SYSTEM, "failed photon path lookup entry allocation"); + + return (void**)(lutEntry -> data ? &lutEntry -> data : NULL); + } + + + + void addPhotonPath (const void *rec, void *data, void *state) + /* Photon path callback for k-NN search to add photon path IDs to lookup + table entries. State points to the associated photon map containing the + path lookup table. The new lookup table entry is initialised with the + specified data pointer */ + { + const Photon *photon = rec; + PhotonMap *pmap = state; + char *key; + LUENT *lutEntry; + + if (!pmap -> pathLUT || !pmap -> pathLUTKeys) + error(INTERNAL, "photon path lookup table not initialised"); + + /* Generate key for photon's path and write into next free slot in key + buffer */ + key = photonPathKey( + photon, pmap -> pathLUTKeys [pmap -> numPathLUTKeys] + ); + + /* Find free LUT entry for generated key */ + if (!(lutEntry = lu_find(pmap -> pathLUT, key))) + /* Outta mammaries */ + error(SYSTEM, "failed photon path lookup entry allocation"); + + if (lutEntry -> data) + /* Occupied */ + error(INTERNAL, "attempt to overwrite photon path lookup entry"); + + #ifdef PMAP_DEBUGPATHS + printf("%s add\n", key); + #endif + + /* Set data and key fields, advance key buffer tail */ + lutEntry -> data = data; + lutEntry -> key = pmap -> pathLUTKeys [pmap -> numPathLUTKeys]; + + if (pmap -> numPathLUTKeys > pmap -> squeue.tail) + /* Number of keys (and distinct path IDs) can never exceed the number + of photons currently in the search queue (+1 for scratchpad) */ + error(INTERNAL, "photon path lookup table key overflow"); + else + pmap -> numPathLUTKeys++; + } + + + + static char emptyKey = '\0'; + + void deletePhotonPath (const void *rec, void *state) + /* Photon path callback for k-NN search to delete photon path IDs from + lookup table. State points to the associated photon map containing the + path ID lookup table. */ + { + const Photon *photon = rec; + PhotonMap *pmap = state; + LUENT *lutEntry; + char *delKey, *relocKey; + + if (!pmap -> pathLUT || !pmap -> pathLUTKeys) + error(INTERNAL, "photon path lookup table not initialised"); + + /* Generate key for photon's path and write into scratchpad at key buffer + tail (will be overwritten) */ + delKey = photonPathKey( + photon, pmap -> pathLUTKeys [pmap -> numPathLUTKeys] + ); + + #ifdef PMAP_DEBUGPATHS + printf("%s del\n", delKey); + #endif + + /* Find photon's LUT entry, check if valid, and get its key buffer slot */ + lutEntry = lu_find(pmap -> pathLUT, delKey); + if (!lutEntry || !lutEntry -> key || !lutEntry -> data) + error(CONSISTENCY, "attempt to delete empty photon path LUT entry"); + + /* Delete photon LUT entry */ + lu_delete(pmap -> pathLUT, delKey); + delKey = lutEntry -> key; + lutEntry -> key = &emptyKey; /* ??? */ + + /* Decrement key buffer tail and get last valid key */ + relocKey = pmap -> pathLUTKeys [--pmap -> numPathLUTKeys]; + if (delKey < relocKey) { + /* Reclaim deleted photon's key buffer slot by overwriting with last + valid key */ + memcpy(delKey, relocKey, PMAP_PATHKEYLEN); + relocKey = delKey; + #ifdef PMAP_DEBUGPATHS + printf("%s reloc\n", relocKey); + #endif + + /* Update the LUT entry owning the relocated key to point to its new + position in the overwritten slot, again checking for validity */ + lutEntry = lu_find(pmap -> pathLUT, relocKey); + if (!lutEntry || !lutEntry -> key || !lutEntry -> data) + error(CONSISTENCY, "attempt to delete empty photon path LUT entry"); + lutEntry -> key = relocKey; + } + } + + + + int clearPhotonPath (const LUENT *lutEntry, void *lut) + /* Clear entry in photon path ID lookup table */ + { + #ifdef PMAP_DEBUGPATHS + printf("%s clr\n", lutEntry -> key); + #endif + lu_delete(lut, lutEntry -> key); + /* Apologies to Don Gregorio for abusing his LUT here... */ + ((LUENT*)lutEntry) -> key = &emptyKey; + + return 0; + } + + + + void resetPhotonPaths (void *state) + /* Reset (clear) all entries in photon path ID lookup table */ + { + PhotonMap *pmap = state; + + if (!pmap -> pathLUT) + error(INTERNAL, "photon path lookup table not initialised"); + + lu_doall(pmap -> pathLUT, clearPhotonPath, pmap -> pathLUT); + + /* Reset photon path ID key buffer */ + pmap -> numPathLUTKeys = 0; + } +#endif /* PMAP_PHOTONFLOW */ + diff --git a/pmapfilt.h b/pmapfilt.h new file mode 100644 index 0000000..43e44b6 --- /dev/null +++ b/pmapfilt.h @@ -0,0 +1,112 @@ +/* + ====================================================================== + Photon map filter callbacks for nearest neighbour search + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) + (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$ +*/ + + +#ifndef PMAP_FILT_H + #define PMAP_FILT_H + + #include "lookup.h" + + /* Forward declarations */ + struct PhotonMap; + + /* Generic callback function container with context/state data to filter + photons during k-NN search */ + typedef struct { + int (*filtFunc)(const void *rec, const void *state); + void *state; + } PhotonSearchFilterCallback; + + /* State data for PhotonSearchFilterCallback, used by filterPhoton() */ + typedef struct { + const struct PhotonMap *pmap; + const float *norm; + } PhotonSearchFilterState; + + #ifdef PMAP_OOC + /* Interface to out-of-core photon map routines */ + typedef PhotonSearchFilterCallback OOC_SearchFilterCallback; + typedef PhotonSearchFilterState OOC_SearchFilterState; + #else + /* Interface to in-core kd-tree photon map routines */ + typedef PhotonSearchFilterCallback kdT_SearchFilterCallback; + typedef PhotonSearchFilterState kdT_SearchFilterState; + #endif + + + int filterPhoton (const void *rec, const void *state); + /* Filter callback for photon kNN search */ + + + #ifdef PMAP_PHOTONFLOW + /* Generic callback function container with context/state data to + maintain table of photon attributes and exclude possible duplicates + during k-NN search */ + typedef struct { + void **(*findFunc)(const void *rec, const void *state); + void (*addFunc)(const void *rec, void *data, void *state); + void (*delFunc)(const void *rec, void *state); + int (*checkFunc)(const void *state); + void *state; + } PhotonSearchAttribCallback; + + #ifdef PMAP_OOC + /* Interface to out-of-core photon map routines */ + typedef PhotonSearchAttribCallback OOC_SearchAttribCallback; + #else + /* Interface to in-core kd-tree photon map routines */ + typedef PhotonSearchAttribCallback kdT_SearchAttribCallback; + #endif + + void initPhotonPaths (void *state); + /* Lazily allocate and initialise photon path ID lookup table and + * associate key buffer */ + + int checkPhotonPaths (const void *state); + /* Run sanity check on photon path ID lookup table */ + + void **findPhotonPath (const void *rec, const void *state); + /* Photon path callback for k-NN search to map photon path IDs to + * search queue nodes via a lookup table. State points to the + * associated photon map containing the path lookup table. Returns + * modifiable pointer to data field of found LUT entry, which in turn + * points to the search queue node for the corresponding record. If + * no entry exists for the path ID, NULL is returned. */ + + void addPhotonPath (const void *rec, void *data, void *state); + /* Photon path callback for k-NN search to add photon path IDs to + lookup table entries. State points to the associated photon map + containing the path lookup table. The new lookup table entry is + initialised with the specified data pointer */ + + void deletePhotonPath (const void *rec, void *state); + /* Photon path callback for k-NN search to delete photon path IDs from + lookup table. State points to the associated photon map containing + the path ID lookup table. */ + + int clearPhotonPath (const LUENT *lutEntry, void *lut); + /* Clear entry in photon path ID lookup table */ + + void resetPhotonPaths (void *state); + /* Reset (clear) all entries in photon path ID lookup table */ + #else + #ifdef PMAP_OOC + typedef void OOC_SearchAttribCallback; + #else + typedef void kdT_SearchAttribCallback; + #endif + #endif /* PMAP_PHOTONFLOW */ +#endif diff --git a/pmapkdt.c b/pmapkdt.c index 04ca423..5cd6b96 100644 --- a/pmapkdt.c +++ b/pmapkdt.c @@ -1,534 +1,776 @@ /* ====================================================================== 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) + 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: pmapkdt.c,v 1.6 2020/04/08 15:14:21 rschregle Exp $ + $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) +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) +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); + 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); + 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); + 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); + 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; + 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 putint(p -> primary, sizeof(p -> primary), 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 p -> primary = getint(sizeof(p -> primary), 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)); + pmap -> squeue.node = calloc( + pmap -> squeue.len, sizeof(PhotonSearchQueueNode) + ); if (!pmap -> squeue.node) error(SYSTEM, "can't allocate photon search queue"); + +#ifdef PMAP_PHOTONFLOW + initPhotonPaths(pmap); +#endif } -static void kdT_FindNearest (PhotonMap *pmap, const float pos [3], - const float norm [3], unsigned long node) -/* Recursive part of kdT_FindPhotons(). Locate pmap -> squeue.len nearest - * neighbours to pos with similar normal and return in search queue starting - * at pmap -> squeue.node. Note that all heap and queue indices are - * 1-based, but accesses to the arrays are 0-based! */ +#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 + + + +static 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_PHOTONFLOW + /* 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_PHOTONFLOW + 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_PHOTONFLOW + 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_PHOTONFLOW + 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_PHOTONFLOW + 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_PHOTONFLOW + 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_PHOTONFLOW) && 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; - unsigned i, j; + 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]; - PhotonSearchQueueNode* sq = pmap -> squeue.node; - const unsigned sqSize = pmap -> squeue.len; + 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, norm, node << 1); + kdT_FindNearest(pmap, pos, filter, attrib, node << 1); if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons) - kdT_FindNearest(pmap, pos, norm, (node << 1) + 1); + kdT_FindNearest(pmap, pos, filter, attrib, (node << 1) + 1); } else { if (node << 1 < pmap -> numPhotons) - kdT_FindNearest(pmap, pos, norm, (node << 1) + 1); + kdT_FindNearest(pmap, pos, filter, attrib, (node << 1) + 1); if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons) - kdT_FindNearest(pmap, pos, norm, node << 1); + kdT_FindNearest(pmap, pos, filter, attrib, node << 1); } - /* Reject photon if normal faces away (ignored for volume photons) with - * tolerance to account for perturbation; note photon normal is coded - * in range [-127,127], hence we factor this in */ - if (norm && DOT(norm, p -> norm) <= PMAP_NORM_TOL * 127 * frandom()) - return; - - if (isContribPmap(pmap)) { - /* Lookup in contribution photon map; filter according to emitting - * light source if contrib list set, else accept all */ - - if (pmap -> srcContrib) { - OBJREC *srcMod; - const int srcIdx = photonSrcIdx(pmap, p); - - if (srcIdx < 0 || srcIdx >= nsources) - error(INTERNAL, "invalid light source index in photon map"); - - srcMod = findmaterial(source [srcIdx].so); - - /* Reject photon if contributions from light source which emitted it - * are not sought */ - if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data) - return; - } - - /* Reject non-caustic photon if lookup for caustic contribs */ - if (pmap -> lookupCaustic & !p -> caustic) - return; - } + /* 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 priority queue */ + /* Accept photon if closer than current max dist & add to search queue */ if (d2 < pmap -> maxDist2) { - if (pmap -> squeue.tail < sqSize) { - /* Priority queue not full; append photon and restore heap */ - i = ++pmap -> squeue.tail; + /* 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; - while (i > 1 && sq [(i >> 1) - 1].dist2 <= d2) { - sq [i - 1].idx = sq [(i >> 1) - 1].idx; - sq [i - 1].dist2 = sq [(i >> 1) - 1].dist2; - i >>= 1; - } - - sq [--i].idx = (PhotonIdx)p; - sq [i].dist2 = d2; - /* Update maxDist if we've just filled the queue */ - if (pmap -> squeue.tail >= pmap -> squeue.len) - pmap -> maxDist2 = sq [0].dist2; - } - else { - /* Priority queue full; replace maximum, restore heap, and - update maxDist */ - i = 1; - - while (i <= sqSize >> 1) { - j = i << 1; - if (j < sqSize && sq [j - 1].dist2 < sq [j].dist2) - j++; - if (d2 >= sq [j - 1].dist2) - break; - sq [i - 1].idx = sq [j - 1].idx; - sq [i - 1].dist2 = sq [j - 1].dist2; - i = j; - } - - sq [--i].idx = (PhotonIdx)p; - sq [i].dist2 = d2; - pmap -> maxDist2 = sq [0].dist2; +#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) +int kdT_FindPhotons ( + struct PhotonMap *pmap, const FVECT pos, const FVECT norm +) { - float p [3], n [3]; - + kdT_SearchFilterCallback filt; + kdT_SearchFilterState filtState; +#ifdef PMAP_PHOTONFLOW + 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); - kdT_FindNearest(pmap, p, norm ? n : NULL, 1); + + /* Set up filter callback */ + filtState.pmap = pmap; + filtState.norm = norm ? n : NULL; + filt.state = &filtState; + filt.filtFunc = filterPhoton; + +#ifdef PMAP_PHOTONFLOW + if (isVolumePmap(pmap)) { + /* Set up path ID callback */ + 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 /* 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 float norm [3], Photon **photon, - unsigned long node) +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, norm, photon, node << 1); + kdT_Find1Nearest(pmap, pos, filter, photon, node << 1); if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons) - kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1); + kdT_Find1Nearest(pmap, pos, filter, photon, (node << 1) + 1); } else { if (node << 1 < pmap -> numPhotons) - kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1); + kdT_Find1Nearest(pmap, pos, filter, photon, (node << 1) + 1); if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons) - kdT_Find1Nearest(pmap, pos, norm, photon, node << 1); + 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 && - (!norm || DOT(norm, p -> norm) > PMAP_NORM_TOL * 127 * frandom())) { - /* Closest photon so far with similar normal. We allow for tolerance - * to account for perturbation in the latter; note the photon normal - * is coded in the range [-127,127], hence we factor this in */ + 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) +int kdT_Find1Photon ( + struct PhotonMap *pmap, const FVECT pos, const FVECT norm, Photon *photon +) { - float p [3], n [3]; - Photon *pnn = NULL; - + 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); - kdT_Find1Nearest(pmap, p, norm ? n : NULL, &pnn, 1); + VCOPY(n, norm); + + /* Set up filter callback */ + filtState.pmap = pmap; + filtState.norm = norm ? n : 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; } } -int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx, - Photon *photon) +/* 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; } diff --git a/pmapkdt.h b/pmapkdt.h index a2f9041..3b75440 100644 --- a/pmapkdt.h +++ b/pmapkdt.h @@ -1,95 +1,104 @@ /* ================================================================== 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) + 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: pmapkdt.h,v 1.2 2020/04/08 15:14:21 rschregle Exp $ */ #ifndef PMAPKDT_H #define PMAPKDT_H -/* #include - #include "fvect.h" */ #include "pmapdata.h" /* Forward declarations to break dependency loop with pmapdata.h */ struct PhotonMap; /* k-d tree as linear array of photons */ typedef struct { Photon *nodes; /* k-d tree as linear array */ } PhotonKdTree; typedef PhotonKdTree PhotonStorage; typedef Photon* PhotonIdx; /* Priority queue node for kd-tree lookups */ typedef struct { PhotonIdx idx; /* Pointer to photon */ float dist2; /* Photon's SQUARED distance to query point */ } kdT_SearchQueueNode; typedef struct { kdT_SearchQueueNode *node; unsigned len, tail; } kdT_SearchQueue; typedef kdT_SearchQueueNode PhotonSearchQueueNode; typedef kdT_SearchQueue PhotonSearchQueue; void kdT_Null (PhotonKdTree *kdt); /* Initialise kd-tree prior to storing photons */ void kdT_BuildPhotonMap (struct PhotonMap *pmap); /* Build a balanced kd-tree pmap -> store from photons in unsorted * heapfile pmap -> heap to guarantee logarithmic search times. The heap * is destroyed on return. */ int kdT_SavePhotons (const struct PhotonMap *pmap, FILE *out); /* Save photons in kd-tree to file. Return -1 on error, else 0 */ int kdT_LoadPhotons (struct PhotonMap *pmap, FILE *in); /* Load photons in kd-tree from file. Return -1 on error, else 0 */ void kdT_InitFindPhotons (struct PhotonMap *pmap); /* Initialise NN search queue prior to calling kdT_FindPhotons() */ - int kdT_FindPhotons (struct PhotonMap* pmap, const FVECT pos, - const FVECT norm); + int kdT_FindPhotons ( + struct PhotonMap* pmap, const FVECT pos, const FVECT norm + ); /* Locate pmap -> squeue.len nearest photons to pos with similar normal * (NULL for volume photons) and return in search queue pmap -> squeue, * starting with the further photon at pmap -> squeue.node. Return -1 * if none found, else 0. */ - int kdT_Find1Photon (struct PhotonMap* pmap, const FVECT pos, - const FVECT norm, Photon *photon); + int kdT_Find1Photon ( + struct PhotonMap* pmap, const FVECT pos, const FVECT norm, + Photon *photon + ); /* Locate single nearest photon to pos with similar normal. Return -1 * if none found, else 0. */ - int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx, - Photon *photon); + int kdT_GetPhoton ( + const struct PhotonMap *pmap, PhotonIdx idx,Photon *photon + ); /* Retrieve photon referenced by idx from kd-tree and return -1 on * error, else 0. */ - Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, - PhotonIdx idx); + Photon *kdT_GetNearestPhoton ( + const kdT_SearchQueue *squeue, PhotonIdx idx + ); /* Retrieve photon from NN search queue after OOC_FindPhotons() */ PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap); /* Return pointer to first photon in kd-Tree */ void kdT_Delete (PhotonKdTree *kdt); /* Self-destruct */ #endif diff --git a/pmapooc.c b/pmapooc.c index ad16ead..7884157 100644 --- a/pmapooc.c +++ b/pmapooc.c @@ -1,338 +1,331 @@ /* ====================================================================== Photon map interface to out-of-core octree Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation (SNSF, #147053) + 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: pmapooc.c,v 1.6 2020/04/08 15:14:21 rschregle Exp $ + $Id: pmapooc.c,v 1.2 2020/11/10 01:08:56 u-no-hoo Exp u-no-hoo $ */ #include "pmapdata.h" /* Includes pmapooc.h */ #include "source.h" #include "otspecial.h" #include "oocsort.h" #include "oocbuild.h" +/* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ + /* Returns photon position as sorting key for OOC_Octree & friends (notably * for Morton code generation). - * !!! Uses type conversion from float via TEMPORARY storage; + * !!! XXX: Uses type conversion from float via TEMPORARY storage; * !!! THIS IS NOT THREAD SAFE! * !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */ RREAL *OOC_PhotonKey (const void *p) { static FVECT photonPos; /* Temp storage for type conversion */ VCOPY(photonPos, ((Photon*)p) -> pos); return photonPos; } #ifdef DEBUG_OOC -static int OOC_CheckKeys (FILE *file, const OOC_Octree *oct) -{ - Photon p, lastp; - RREAL *key; - OOC_MortonIdx zIdx, lastzIdx = 0; - - rewind(file); - memset(&lastp, 0, sizeof(lastp)); - - while (fread(&p, sizeof(p), 1, file) > 0) { - key = OOC_PhotonKey(&p); - zIdx = OOC_KEY2MORTON(key, oct); + static void OOC_CheckKeys (FILE *file, const OOC_Octree *oct) + { + Photon p, lastp; + RREAL *key; + OOC_MortonIdx zIdx, lastzIdx = 0; - if (zIdx < lastzIdx) { - error(INTERNAL, "photons not sorted"); - return -1; - } + rewind(file); + memset(&lastp, 0, sizeof(lastp)); - if (zIdx == lastzIdx) { - sprintf(errmsg, "identical key %021ld " - "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]", - zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], - p.pos [0], p.pos [1], p.pos [2]); - error(WARNING, errmsg); - } - - lastzIdx = zIdx; - memcpy(&lastp, &p, sizeof(p)); + while (fread(&p, sizeof(p), 1, file) > 0) { + key = OOC_PhotonKey(&p); + zIdx = OOC_KEY2MORTON(key, oct); + + if (zIdx < lastzIdx) + error(INTERNAL, "photons not sorted"); + + if (zIdx == lastzIdx) { + sprintf(errmsg, "identical key %021ld " + "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]", + zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], + p.pos [0], p.pos [1], p.pos [2]); + error(WARNING, errmsg); + } + + lastzIdx = zIdx; + memcpy(&lastp, &p, sizeof(p)); + } } - - return 0; -} #endif void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc) { FILE *leafFile; char leafFname [1024]; FVECT d, octOrg; int i; RREAL octSize = 0; /* Determine octree size and origin from pmap extent and init octree */ VCOPY(octOrg, pmap -> minPos); VSUB(d, pmap -> maxPos, pmap -> minPos); for (i = 0; i < 3; i++) if (octSize < d [i]) octSize = d [i]; if (octSize < FTINY) error(INTERNAL, "zero octree size in OOC_BuildPhotonMap"); /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1); if (!(leafFile = fopen(leafFname, "w+b"))) error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Sorting photons by Morton code...\n"); #endif /* Sort photons in heapfile by Morton code and write to leaf file */ if (OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey)) error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap"); /* Init and build octree */ OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey, leafFile); #ifdef DEBUG_OOC eputs("Checking leaf file consistency...\n"); OOC_CheckKeys(leafFile, &pmap -> store); eputs("Building out-of-core octree...\n"); #endif if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH)) error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Checking out-of-core octree consistency...\n"); if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), octOrg, octSize, 0)) error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri"); #endif } +/* PHOTON MAP I/O ROUTINES ------------------------------------------ */ + int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out) { return OOC_SaveOctree(&pmap -> store, out); } int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *nodeFile) { FILE *leafFile; char leafFname [1024]; /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1); if (!(leafFile = fopen(leafFname, "r"))) error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons"); if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile)) return -1; return 0; } +/* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ + void OOC_InitFindPhotons (struct PhotonMap *pmap) { - if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1, - sizeof(Photon))) + if (OOC_InitNearest( + &pmap -> squeue, pmap -> maxGather + 1, sizeof(Photon) + )) error(SYSTEM, "can't allocate photon search queue"); + +#ifdef PMAP_PHOTONFLOW + if (isVolumePmap(pmap) && photonFlow) + initPhotonPaths(pmap); +#endif } static void OOC_InitPhotonCache (struct PhotonMap *pmap) /* Initialise OOC photon cache */ { static char warn = 1; if (!pmap -> store.cache && !pmap -> numDensity) { if (pmapCacheSize > 0) { const unsigned pageSize = pmapCachePageSize * pmap -> maxGather, numPages = pmapCacheSize / pageSize; /* Allocate & init I/O cache in octree */ pmap -> store.cache = malloc(sizeof(OOC_Cache)); if (!pmap -> store.cache || - OOC_CacheInit(pmap -> store.cache, numPages, pageSize, - sizeof(Photon))) { + OOC_CacheInit( + pmap -> store.cache, numPages, pageSize, sizeof(Photon) + )) { error(SYSTEM, "failed OOC photon map cache init"); } } else if (warn) { error(WARNING, "OOC photon map cache DISABLED"); warn = 0; } } } -typedef struct { - const PhotonMap *pmap; - const float *norm; -} OOC_FilterData; - - - -int OOC_FilterPhoton (void *p, void *fd) -/* Filter callback for photon kNN search, used by OOC_FindNearest() */ +int OOC_FindPhotons ( + struct PhotonMap *pmap, const FVECT pos, const FVECT norm +) { - const Photon *photon = p; - const OOC_FilterData *filtData = fd; - const PhotonMap *pmap = filtData -> pmap; - - /* Reject photon if normal faces away (ignored for volume photons) with - * tolerance to account for perturbation; note photon normal is coded - * in range [-127,127], hence we factor this in */ - if (filtData -> norm && - DOT(filtData->norm, photon->norm) <= PMAP_NORM_TOL * 127 * frandom()) - return 0; - - if (isContribPmap(pmap)) { - /* Lookup in contribution photon map; filter according to emitting - * light source if contrib list set, else accept all */ - - if (pmap -> srcContrib) { - OBJREC *srcMod; - const int srcIdx = photonSrcIdx(pmap, photon); - - if (srcIdx < 0 || srcIdx >= nsources) - error(INTERNAL, "invalid light source index in photon map"); - - srcMod = findmaterial(source [srcIdx].so); - - /* Reject photon if contributions from light source which emitted - * it are not sought */ - if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data) - return 0; - } - - /* Reject non-caustic photon if lookup for caustic contribs */ - if (pmap -> lookupCaustic && !photon -> caustic) - return 0; - } - - /* Accept photon */ - return 1; -} - - - -int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm) -{ - OOC_SearchFilter filt; - OOC_FilterData filtData; - float n [3]; + OOC_SearchFilterCallback filt; + OOC_SearchFilterState filtState; +#ifdef PMAP_PHOTONFLOW + OOC_SearchAttribCallback paths, *pathsPtr = NULL; +#endif + float n [3]; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ - filtData.pmap = pmap; if (norm) - VCOPY(n, norm); - filtData.norm = norm ? n : NULL; - filt.data = &filtData; - filt.func = OOC_FilterPhoton; + VCOPY(n, norm); + + filtState.pmap = pmap; + filtState.norm = norm ? n : NULL; + filt.state = &filtState; + filt.filtFunc = filterPhoton; + +#ifdef PMAP_PHOTONFLOW + if (isVolumePmap(pmap) && photonFlow) { + /* Set up path ID callback */ + paths.state = pmap; + paths.findFunc = findPhotonPath; + paths.addFunc = addPhotonPath; + paths.delFunc = deletePhotonPath; + paths.checkFunc = checkPhotonPaths; + pathsPtr = &paths; + resetPhotonPaths(pmap); + } - pmap -> maxDist2 = OOC_FindNearest(&pmap -> store, - OOC_ROOT(&pmap -> store), 0, - pmap -> store.org, pmap -> store.size, - pos, &filt, &pmap -> squeue, - pmap -> maxDist2); + pmap -> maxDist2 = OOC_FindNearest( + &pmap -> store, OOC_ROOT(&pmap -> store), 0, + pmap -> store.org, pmap -> store.size, pos, &filt, pathsPtr, + &pmap -> squeue, pmap -> maxDist2 + ); +#else + pmap -> maxDist2 = OOC_FindNearest( + &pmap -> store, OOC_ROOT(&pmap -> store), 0, + pmap -> store.org, pmap -> store.size, pos, &filt, NULL, + &pmap -> squeue, pmap -> maxDist2 + ); +#endif /* PMAP_PHOTONFLOW */ if (pmap -> maxDist2 < 0) error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons"); - + /* Return success or failure (empty queue => none found) */ return pmap -> squeue.tail ? 0 : -1; } -int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, - const FVECT norm, Photon *photon) +int OOC_Find1Photon ( + struct PhotonMap* pmap, const FVECT pos, const FVECT norm, Photon *photon +) { - OOC_SearchFilter filt; - OOC_FilterData filtData; - float n [3], maxDist2; + OOC_SearchFilterCallback filt; + OOC_SearchFilterState filtState; + float n [3], maxDist2; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ - filtData.pmap = pmap; if (norm) VCOPY(n, norm); - filtData.norm = norm ? n : NULL; - filt.data = &filtData; - filt.func = OOC_FilterPhoton; + + filtState.pmap = pmap; + filtState.norm = norm ? n : NULL; + filt.state = &filtState; + filt.filtFunc = filterPhoton; - maxDist2 = OOC_Find1Nearest(&pmap -> store, - OOC_ROOT(&pmap -> store), 0, - pmap -> store.org, pmap -> store.size, - pos, &filt, photon, pmap -> maxDist2); + maxDist2 = OOC_Find1Nearest( + &pmap -> store, OOC_ROOT(&pmap -> store), 0, + pmap -> store.org, pmap -> store.size,pos, + &filt, photon, pmap -> maxDist2 + ); if (maxDist2 < 0) error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon"); - + if (maxDist2 >= pmap -> maxDist2) /* No photon found => failed */ return -1; else { /* Set photon distance => success */ pmap -> maxDist2 = maxDist2; return 0; } } +/* PHOTON ACCESS ROUTINES ------------------------------------------ */ + int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon) { return OOC_GetData(&pmap -> store, idx, photon); } Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { return OOC_GetNearest(squeue, idx); } PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap) { return 0; } + diff --git a/pmapooc.h b/pmapooc.h index be1a20d..328b3bd 100644 --- a/pmapooc.h +++ b/pmapooc.h @@ -1,85 +1,90 @@ /* ================================================================== Photon map interface to out-of-core octree Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation (SNSF, #147053) + 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: pmapooc.h,v 1.2 2020/04/08 15:14:21 rschregle Exp $ + $Id: pmapooc.h,v 1.1 2020/10/29 16:20:11 u-no-hoo Exp u-no-hoo $ */ #ifndef PMAPOOC_H #define PMAPOOC_H #include "oocnn.h" - /* Suffixes for octree filenames */ /* #define PMAP_OOC_NODESUFFIX ".node" #define PMAP_OOC_HEAPSUFFIX ".heap" */ #define PMAP_OOC_LEAFSUFFIX ".leaf" /* Out-of-core octree constants */ #define PMAP_OOC_NUMBLK 32 /* Num blocks for external sort */ #define PMAP_OOC_BLKSIZE 1e8 /* Block size for external sort */ #define PMAP_OOC_LEAFMAX (OOC_OCTCNT_MAX) /* Max photons per leaf */ #define PMAP_OOC_MAXDEPTH (OOC_MORTON_BITS) /* Max octree depth */ - typedef OOC_SearchQueueNode PhotonSearchQueueNode; typedef OOC_SearchQueue PhotonSearchQueue; typedef OOC_Octree PhotonStorage; typedef unsigned PhotonIdx; /* Forward declarations to break dependency loop with pmapdata.h */ struct PhotonMap; - void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc); /* Build out-of-core octree pmap -> store from photons in unsorted * heapfile pmap -> heap and generate nodes and leaf file with prefix * pmap -> fileName. Photon map construction may be parallelised if * numProc > 1, if supported. The heap is destroyed on return. */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out); /* Save photons in out-of-core octree to file. Return -1 on error, * else 0 */ int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *in); /* Load photons in out-of-core octree from file. Return -1 on error, * else 0 */ void OOC_InitFindPhotons (struct PhotonMap *pmap); /* Initialise NN search queue prior to calling kdT_FindPhotons() */ - int OOC_FindPhotons (struct PhotonMap* pmap, const FVECT pos, - const FVECT norm); + int OOC_FindPhotons ( + struct PhotonMap* pmap, const FVECT pos, const FVECT norm + ); /* Locate pmap -> squeue.len nearest photons to pos with similar normal * (NULL for volume photons) and return in search queue pmap -> squeue, * starting with the further photon at pmap -> squeue.node. Return -1 * if none found, else 0. */ - int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, - const FVECT norm, Photon *photon); + int OOC_Find1Photon ( + struct PhotonMap* pmap, const FVECT pos, const FVECT norm, + Photon *photon + ); /* Locate single nearest photon to pos with similar normal. Return -1 * if none found, else 0. */ - int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, - Photon *photon); + int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon); /* Retrieve photon referenced by idx from leaf file and return -1 on * error, else 0. */ - Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, - PhotonIdx idx); + Photon *OOC_GetNearestPhoton ( + const PhotonSearchQueue *squeue, PhotonIdx idx + ); /* Retrieve photon from NN search queue after OOC_FindPhotons() */ PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap); /* Return index to first photon in octree */ + #endif diff --git a/pmapopt.c b/pmapopt.c index 6333a55..2179e01 100644 --- a/pmapopt.c +++ b/pmapopt.c @@ -1,121 +1,157 @@ #ifndef lint static const char RCSid[] = "$Id: pmapopt.c,v 2.10 2020/06/15 22:18:57 rschregle Exp $"; #endif /* ================================================================== Photon map interface to RADIANCE render options 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) + 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: pmapopt.c,v 2.10 2020/06/15 22:18:57 rschregle Exp $ */ #include "ray.h" #include "pmapparm.h" int getPmapRenderOpt (int ac, char *av []) /* Parse next render option for photon map; interface to getrenderopt(); * return -1 if parsing failed, else number of parameters consumed */ { #define check(ol,al) (av[0][ol] || badarg(ac-1,av+1,al)) + /* Evaluate boolean option, setting var accordingly */ + #define check_bool(olen, var) switch (av [0][olen]) { \ + case '\0': \ + var = !var; break; \ + case 'y': case 'Y': case 't': case 'T': case '+': case '1': \ + var = 1; break; \ + case 'n': case 'N': case 'f': case 'F': case '-': case '0': \ + var = 0; break; \ + default: \ + return -1; \ + } + static int t = -1; /* pmap parameter index */ if (ac < 1 || !av [0] || av [0][0] != '-') return -1; switch (av [0][1]) { case 'a': switch (av [0][2]) { case 'p': /* photon map */ - /* Asking for photon map, ergo ambounce != 0 */ + /* Photon map bumps up ambounce from its default 0; SHOULD + THIS REMAIN THE DEFAULT BEHAVIOUR? In some cases it's + desirable to render photons directly via ab -1 */ ambounce += (ambounce == 0); - + +#ifdef PMAP_PHOTONFLOW + if (av [0][3] == 'F') { + /* Evaluate immediate irrad. from volume photon map (if + * specified) as a lightfield, ignoring part. medium */ + /* FOR ZE EKSPERTZ ONLY! */ + check_bool(4, photonFlow); + return 0; + } +#endif + if (!check(3, "s")) { /* File -> assume bwidth = 1 or precomputed pmap */ if (++t >= NUM_PMAP_TYPES) error(USER, "too many photon maps"); pmapParams [t].fileName = savqstr(av [1]); pmapParams [t].minGather = pmapParams [t].maxGather = defaultGather; } else return -1; if (!check(3, "si")) { /* File + bandwidth */ pmapParams [t].minGather = pmapParams [t].maxGather = atoi(av [2]); if (!pmapParams [t].minGather) return -1; } else { sprintf(errmsg, "missing photon lookup bandwidth, " "defaulting to %d", defaultGather); error(WARNING, errmsg); return 1; } if (!check(3, "sii")) { /* File + min bwidth + max bwidth -> bias compensation */ pmapParams [t].maxGather = atoi(av [3]); if (pmapParams [t].minGather >= pmapParams [t].maxGather) return -1; return 3; } else return 2; case 'm': /* Fixed max photon search radius */ if (check(3, "f") || (maxDistFix = atof(av [1])) <= 0) error(USER, "invalid max photon search radius"); return 1; #ifdef PMAP_OOC case 'c': /* OOC pmap cache page size ratio */ if (check(3, "f") || (pmapCachePageSize = atof(av [1])) <= 0) error(USER, "invalid photon cache page size ratio"); return 1; case 'C': /* OOC pmap cache size (in photons); 0 = no caching */ if (check(3, "s")) error(USER, "invalid number of cached photons"); /* Parsing failure sets to zero and disables caching */ pmapCacheSize = parseMultiplier(av [1]); return 1; #endif } } #undef check /* Unknown option */ return -1; } void printPmapDefaults () /* Print defaults for photon map render options */ { printf("-am %.1f\t\t\t\t# max photon search radius\n", maxDistFix); #ifdef PMAP_OOC printf("-ac %.1f\t\t\t\t# photon cache page size ratio\n", pmapCachePageSize); printf("-aC %ld\t\t\t# num cached photons\n", pmapCacheSize); #endif +#ifdef PMAP_PHOTONFLOW + printf( + photonFlow ? "-apF+\t\t\t\t# photon flow from volume pmap on\n" + : "-apF-\t\t\t\t# photon flow from volume pmap off\n" + ); +#endif } diff --git a/pmapparm.c b/pmapparm.c index 21a51ad..6508eff 100644 --- a/pmapparm.c +++ b/pmapparm.c @@ -1,110 +1,120 @@ #ifndef lint static const char RCSid[] = "$Id: pmapparm.c,v 2.9 2018/02/02 19:47:55 rschregle Exp $"; #endif /* ====================================================================== Parameters for photon map generation and rendering; used by mkpmap and rpict/rvu/rtrace. 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) + 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: pmapparm.c,v 2.9 2018/02/02 19:47:55 rschregle Exp $ */ #include "pmapparm.h" #include "pmapdata.h" #include float pdfSamples = 1000, /* PDF samples per steradian */ finalGather = 0.25, /* fraction of global photons for irradiance precomputation */ preDistrib = 0.25, /* fraction of num photons for distribution prepass */ gatherTolerance = 0.5, /* Photon map lookup tolerance; lookups returning fewer than this fraction of minGather/maxGather are restarted with a larger search radius */ maxDistFix = 0, /* Static maximum photon search radius (radius is adaptive if this is zero) */ photonMaxDist = 0; /* Maximum cumulative distance of photon path */ #ifdef PMAP_OOC float pmapCachePageSize = 8; /* OOC cache pagesize as multiple * of maxGather */ unsigned long pmapCacheSize = 1e6; /* OOC cache size in photons */ #endif - /* Regions of interest */ -unsigned pmapNumROI = 0; -PhotonMapROI *pmapROI = NULL; - +unsigned pmapNumROI = 0; +PhotonMapROI *pmapROI = NULL; -unsigned verbose = 0; /* Verbose console output */ +unsigned verbose = 0; /* Verbose console output */ unsigned long photonMaxBounce = 5000; /* Runaway photon bounce limit */ unsigned photonRepTime = 0, /* Seconds between reports */ maxPreDistrib = 4, /* Max predistrib passes */ defaultGather = 40; /* Default num photons for lookup */ + +#ifdef PMAP_PHOTONFLOW +unsigned photonFlow = 0; /* Evaluate immediate irradiance from + volume photons as lightfield, + ignoring participating medium */ +#endif /* Per photon map params */ PhotonMapParams pmapParams [NUM_PMAP_TYPES] = { {NULL, 0, 0, 0}, {NULL, 0, 0, 0}, {NULL, 0, 0, 0}, {NULL, 0, 0, 0}, {NULL, 0, 0, 0}, {NULL, 0, 0, 0} }; int setPmapParam (PhotonMap** pm, const PhotonMapParams* parm) { if (parm && parm -> fileName) { if (!(*pm = (PhotonMap*)malloc(sizeof(PhotonMap)))) error(INTERNAL, "failed to allocate photon map"); (*pm) -> fileName = parm -> fileName; (*pm) -> minGather = parm -> minGather; (*pm) -> maxGather = parm -> maxGather; (*pm) -> distribTarget = parm -> distribTarget; (*pm) -> maxDist0 = FHUGE; (*pm) -> srcContrib = NULL; return 1; } return 0; } unsigned long parseMultiplier (const char *num) { unsigned long mult = 1; const int strEnd = strlen(num) - 1; if (strEnd <= 0) return 0; if (!isdigit(num [strEnd])) switch (toupper(num [strEnd])) { case 'G': mult *= 1000; case 'M': mult *= 1000; case 'K': mult *= 1000; break; default : error(USER, "unknown multiplier"); } return (unsigned long)(mult * atof(num)); } diff --git a/pmapparm.h b/pmapparm.h index ad0c775..2ee7c6f 100644 --- a/pmapparm.h +++ b/pmapparm.h @@ -1,72 +1,82 @@ /* RCSid $Id: pmapparm.h,v 2.9 2018/02/02 19:47:55 rschregle Exp $ */ /* ====================================================================== Parameters for photon map generation and rendering; used by mkpmap and rpict/rvu/rtrace. 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) + 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: pmapparm.h,v 2.9 2018/02/02 19:47:55 rschregle Exp $ */ #ifndef PMAPPARAMS_H #define PMAPPARAMS_H #include "pmaptype.h" /* Struct for passing params per photon map from rpict/rtrace/rvu */ typedef struct { char *fileName; /* Photon map file */ unsigned minGather, maxGather; /* Num photons to gather */ unsigned long distribTarget; /* Num photons to store */ } PhotonMapParams; /* Bounding box for region of interest */ typedef struct { float min [3], max [3]; } PhotonMapROI; extern PhotonMapParams pmapParams [NUM_PMAP_TYPES]; /* Macros for type specific photon map parameters */ #define globalPmapParams (pmapParams [PMAP_TYPE_GLOBAL]) #define preCompPmapParams (pmapParams [PMAP_TYPE_PRECOMP]) #define causticPmapParams (pmapParams [PMAP_TYPE_CAUSTIC]) #define volumePmapParams (pmapParams [PMAP_TYPE_VOLUME]) #define directPmapParams (pmapParams [PMAP_TYPE_DIRECT]) #define contribPmapParams (pmapParams [PMAP_TYPE_CONTRIB]) extern float pdfSamples, preDistrib, finalGather, gatherTolerance, maxDistFix, pmapMaxDist, photonMaxDist; extern unsigned long photonHeapSizeInc, photonMaxBounce; extern unsigned photonRepTime, maxPreDistrib, defaultGather, verbose; extern unsigned pmapNumROI; extern PhotonMapROI *pmapROI; #ifdef PMAP_OOC extern float pmapCachePageSize; extern unsigned long pmapCacheSize; #endif +#ifdef PMAP_PHOTONFLOW + extern unsigned photonFlow; +#endif + struct PhotonMap; int setPmapParam (struct PhotonMap **pm, const PhotonMapParams *parm); /* Allocate photon map and set its parameters from parm */ unsigned long parseMultiplier (const char *num); /* Evaluate numeric parameter string with optional multiplier suffix (G = 10^9, M = 10^6, K = 10^3). Returns 0 if parsing fails. */ #endif diff --git a/pmapsrc.c b/pmapsrc.c index bf0bdce..b42f51b 100644 --- a/pmapsrc.c +++ b/pmapsrc.c @@ -1,962 +1,977 @@ #ifndef lint static const char RCSid[] = "$Id: pmapsrc.c,v 2.19 2020/08/10 19:51:20 rschregle Exp $"; #endif /* ====================================================================== Photon map support routines for emission from light sources Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation (DFG) - under the FARESYS project. + 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). + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, - supported by the JSPS KAKENHI Grant Number JP19KK0115. + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapsrc.c,v 2.19 2020/08/10 19:51:20 rschregle Exp $" */ #include "pmapsrc.h" #include "pmap.h" #include "pmaprand.h" #include "otypes.h" #include "otspecial.h" +#ifdef PMAP_PHOTONFLOW + #include "pmapdens.h" +#endif /* List of photon port modifier names */ char *photonPortList [MAXSET + 1] = {NULL}; /* Photon port objects (with modifiers in photonPortMods) */ SRCREC *photonPorts = NULL; unsigned numPhotonPorts = 0; void (*photonPartition [NUMOTYPE]) (EmissionMap*); void (*photonOrigin [NUMOTYPE]) (EmissionMap*); /* PHOTON PORT SUPPORT ROUTINES ------------------------------------------ */ /* Get/set photon port orientation flags from/in source flags. * HACK: the port orientation flags are embedded in the source flags and * shifted so they won't clobber the latter, since these are interpreted * by the *PhotonPartition() and *PhotonOrigin() routines! */ #define PMAP_SETPORTFLAGS(portdir) ((int)(portdir) << 14) #define PMAP_GETPORTFLAGS(sflags) ((sflags) >> 14) /* Set number of source partitions. * HACK: this is doubled if the source acts as a bidirectionally * emitting photon port, resulting in alternating front/backside partitions, * although essentially each partition is just used twice with opposing * normals. */ #define PMAP_SETNUMPARTITIONS(emap) ( \ (emap) -> numPartitions <<= ( \ (emap) -> port && \ PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \ ) \ ) /* Get current source partition and numer of partitions * HACK: halve the number partitions if the source acts as a bidrectionally * emitting photon port, since each partition is used twice with opposing * normals. */ #define PMAP_GETNUMPARTITIONS(emap) (\ (emap) -> numPartitions >> ( \ (emap) -> port && \ PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \ ) \ ) #define PMAP_GETPARTITION(emap) ( \ (emap) -> partitionCnt >> ( \ (emap) -> port && \ PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \ ) \ ) void getPhotonPorts (char **portList) /* Find geometry declared as photon ports from modifiers in portList */ { OBJECT i; OBJREC *obj, *mat; int mLen; char **lp; /* Init photon port objects */ photonPorts = NULL; if (!portList [0]) return; for (i = numPhotonPorts = 0; i < nobjects; i++) { obj = objptr(i); mat = findmaterial(obj); /* Check if object is a surface and NOT a light source (duh) and * resolve its material (if any) via any aliases, then check for * inclusion in modifier list; note use of strncmp() to ignore port * flags */ if (issurface(obj -> otype) && mat && !islight(mat -> otype)) { mLen = strlen(mat -> oname); for (lp = portList; *lp && strncmp(mat -> oname, *lp, mLen); lp++); if (*lp) { /* Add photon port */ photonPorts = (SRCREC*)realloc( photonPorts, (numPhotonPorts + 1) * sizeof(SRCREC) ); if (!photonPorts) error(USER, "can't allocate photon ports"); photonPorts [numPhotonPorts].so = obj; /* Extract port orientation flags and embed in source flags. * Note setsource() combines (i.e. ORs) these with the actual * source flags below. */ photonPorts [numPhotonPorts].sflags = PMAP_SETPORTFLAGS((*lp) [mLen]); if (!sfun [obj -> otype].of || !sfun[obj -> otype].of -> setsrc) objerror(obj, USER, "illegal photon port"); setsource(photonPorts + numPhotonPorts, obj); numPhotonPorts++; } } } if (!numPhotonPorts) error(USER, "no valid photon ports found"); } static void setPhotonPortNormal (EmissionMap* emap) /* Set normal for current photon port partition (if defined) based on its * orientation */ { int i, portFlags; if (emap -> port) { /* Extract photon port orientation flags, set surface normal as follows: -- Port oriented forwards --> flip surface normal to point outwards, since normal points inwards per mkillum convention) -- Port oriented backwards --> surface normal is NOT flipped, since it already points inwards. -- Port is bidirectionally/bilaterally oriented --> flip normal based on the parity of the current partition emap -> partitionCnt. In this case, photon emission alternates between port front/back faces for consecutive partitions. */ portFlags = PMAP_GETPORTFLAGS(emap -> port -> sflags); if ( portFlags == PMAP_PORTFWD || portFlags == PMAP_PORTBI && !(emap -> partitionCnt & 1) ) for (i = 0; i < 3; i++) emap -> ws [i] = -emap -> ws [i]; } } /* SOURCE / PHOTON PORT PARTITIONING ROUTINES----------------------------- */ static int flatPhotonPartition2 ( EmissionMap* emap, unsigned long mp, FVECT cent, FVECT u, FVECT v, double du2, double dv2 ) /* Recursive part of flatPhotonPartition(..) */ { FVECT newct, newax; unsigned long npl, npu; if (mp > emap -> maxPartitions) { /* Enlarge partition array */ emap -> maxPartitions <<= 1; emap -> partitions = (unsigned char*)realloc( emap -> partitions, emap -> maxPartitions >> 1 ); if (!emap -> partitions) error(USER, "can't allocate source partitions"); memset( emap -> partitions + (emap -> maxPartitions >> 2), 0, emap -> maxPartitions >> 2 ); } if (du2 * dv2 <= 1) { /* hit limit? */ setpart(emap -> partitions, emap -> partitionCnt, S0); emap -> partitionCnt++; return 1; } if (du2 > dv2) { /* subdivide in U */ setpart(emap -> partitions, emap -> partitionCnt, SU); emap -> partitionCnt++; newax [0] = 0.5 * u [0]; newax [1] = 0.5 * u [1]; newax [2] = 0.5 * u [2]; u = newax; du2 *= 0.25; } else { /* subdivide in V */ setpart(emap -> partitions, emap -> partitionCnt, SV); emap -> partitionCnt++; newax [0] = 0.5 * v [0]; newax [1] = 0.5 * v [1]; newax [2] = 0.5 * v [2]; v = newax; dv2 *= 0.25; } /* lower half */ newct [0] = cent [0] - newax [0]; newct [1] = cent [1] - newax [1]; newct [2] = cent [2] - newax [2]; npl = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2); /* upper half */ newct [0] = cent [0] + newax [0]; newct [1] = cent [1] + newax [1]; newct [2] = cent [2] + newax [2]; npu = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2); /* return total */ return npl + npu; } static void flatPhotonPartition (EmissionMap* emap) /* Partition flat source for photon emission */ { RREAL *vp; double du2, dv2; memset(emap -> partitions, 0, emap -> maxPartitions >> 1); emap -> partArea = srcsizerat * thescene.cusize; emap -> partArea *= emap -> partArea; vp = emap -> src -> ss [SU]; du2 = DOT(vp, vp) / emap -> partArea; vp = emap -> src -> ss [SV]; dv2 = DOT(vp, vp) / emap -> partArea; emap -> partitionCnt = 0; emap -> numPartitions = flatPhotonPartition2( emap, 1, emap -> src -> sloc, emap -> src -> ss [SU], emap -> src -> ss [SV], du2, dv2 ); emap -> partitionCnt = 0; emap -> partArea = emap -> src -> ss2 / emap -> numPartitions; } static void sourcePhotonPartition (EmissionMap* emap) /* Partition scene cube faces or photon port for photon emission from distant source */ { if (emap -> port) { /* Relay partitioning to photon port */ SRCREC *src = emap -> src; emap -> src = emap -> port; photonPartition [emap -> src -> so -> otype] (emap); PMAP_SETNUMPARTITIONS(emap); emap -> src = src; } else { /* No photon ports defined; partition scene cube faces (SUBOPTIMAL) */ memset(emap -> partitions, 0, emap -> maxPartitions >> 1); setpart(emap -> partitions, 0, S0); emap -> partitionCnt = 0; emap -> numPartitions = 1 / srcsizerat; emap -> numPartitions *= emap -> numPartitions; emap -> partArea = sqr(thescene.cusize) / emap -> numPartitions; emap -> numPartitions *= 6; } } static void spherePhotonPartition (EmissionMap* emap) /* Partition spherical source into equal solid angles using uniform mapping for photon emission */ { unsigned numTheta, numPhi; memset(emap -> partitions, 0, emap -> maxPartitions >> 1); setpart(emap -> partitions, 0, S0); emap -> partArea = 4 * PI * sqr(emap -> src -> srad); emap -> numPartitions = emap -> partArea / sqr(srcsizerat * thescene.cusize); numTheta = max(sqrt(2 * emap -> numPartitions / PI) + 0.5, 1); numPhi = 0.5 * PI * numTheta + 0.5; emap -> numPartitions = (unsigned long)numTheta * numPhi; emap -> partitionCnt = 0; emap -> partArea /= emap -> numPartitions; } static int cylPhotonPartition2 ( EmissionMap* emap, unsigned long mp, FVECT cent, FVECT axis, double d2 ) /* Recursive part of cyPhotonPartition(..) */ { FVECT newct, newax; unsigned long npl, npu; if (mp > emap -> maxPartitions) { /* Enlarge partition array */ emap -> maxPartitions <<= 1; emap -> partitions = (unsigned char*)realloc( emap -> partitions, emap -> maxPartitions >> 1 ); if (!emap -> partitions) error(USER, "can't allocate source partitions"); memset( emap -> partitions + (emap -> maxPartitions >> 2), 0, emap -> maxPartitions >> 2 ); } if (d2 <= 1) { /* hit limit? */ setpart(emap -> partitions, emap -> partitionCnt, S0); emap -> partitionCnt++; return 1; } /* subdivide */ setpart(emap -> partitions, emap -> partitionCnt, SU); emap -> partitionCnt++; newax [0] = 0.5 * axis [0]; newax [1] = 0.5 * axis [1]; newax [2] = 0.5 * axis [2]; d2 *= 0.25; /* lower half */ newct [0] = cent [0] - newax [0]; newct [1] = cent [1] - newax [1]; newct [2] = cent [2] - newax [2]; npl = cylPhotonPartition2(emap, mp << 1, newct, newax, d2); /* upper half */ newct [0] = cent [0] + newax [0]; newct [1] = cent [1] + newax [1]; newct [2] = cent [2] + newax [2]; npu = cylPhotonPartition2(emap, mp << 1, newct, newax, d2); /* return total */ return npl + npu; } static void cylPhotonPartition (EmissionMap* emap) /* Partition cylindrical source for photon emission */ { double d2; memset(emap -> partitions, 0, emap -> maxPartitions >> 1); d2 = srcsizerat * thescene.cusize; d2 = PI * emap -> src -> ss2 / (2 * emap -> src -> srad * sqr(d2)); d2 *= d2 * DOT(emap -> src -> ss [SU], emap -> src -> ss [SU]); emap -> partitionCnt = 0; emap -> numPartitions = cylPhotonPartition2( emap, 1, emap -> src -> sloc, emap -> src -> ss [SU], d2 ); emap -> partitionCnt = 0; emap -> partArea = PI * emap -> src -> ss2 / emap -> numPartitions; } /* PHOTON ORIGIN ROUTINES ------------------------------------------------ */ static void flatPhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on flat light source surface. Also sets source aperture and sampling hemisphere axes at origin */ { int i, cent[3], size[3], parr[2]; FVECT vpos; cent [0] = cent [1] = cent [2] = 0; size [0] = size [1] = size [2] = emap -> maxPartitions; parr [0] = 0; parr [1] = PMAP_GETPARTITION(emap); if (!skipparts(cent, size, parr, emap -> partitions)) error(CONSISTENCY, "bad source partition in flatPhotonOrigin"); vpos [0] = (1 - 2 * pmapRandom(partState)) * size [0] / emap -> maxPartitions; vpos [1] = (1 - 2 * pmapRandom(partState)) * size [1] / emap -> maxPartitions; vpos [2] = 0; for (i = 0; i < 3; i++) vpos [i] += (double)cent [i] / emap -> maxPartitions; /* Get origin */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> src -> sloc [i] + vpos [SU] * emap -> src -> ss [SU][i] + vpos [SV] * emap -> src -> ss [SV][i] + vpos [SW] * emap -> src -> ss [SW][i]; /* Get surface axes */ VCOPY(emap -> us, emap -> src -> ss [SU]); normalize(emap -> us); VCOPY(emap -> ws, emap -> src -> ss [SW]); /* Flip normal emap -> ws if port and required by its orientation */ setPhotonPortNormal(emap); fcross(emap -> vs, emap -> ws, emap -> us); /* Get hemisphere axes & aperture */ if (emap -> src -> sflags & SSPOT) { VCOPY(emap -> wh, emap -> src -> sl.s -> aim); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI); } else { VCOPY(emap -> uh, emap -> us); VCOPY(emap -> vh, emap -> vs); VCOPY(emap -> wh, emap -> ws); emap -> cosThetaMax = 0; } } static void spherePhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on spherical light source. Also sets source aperture and sampling hemisphere axes at origin */ { int i = 0; unsigned numTheta, numPhi, t, p; RREAL cosTheta, sinTheta, phi; /* Get current partition */ numTheta = max(sqrt(2 * PMAP_GETNUMPARTITIONS(emap) / PI) + 0.5, 1); numPhi = 0.5 * PI * numTheta + 0.5; t = PMAP_GETPARTITION(emap) / numPhi; p = PMAP_GETPARTITION(emap) - t * numPhi; emap -> ws [2] = cosTheta = 1 - 2 * (t + pmapRandom(partState)) / numTheta; sinTheta = sqrt(1 - sqr(cosTheta)); phi = 2 * PI * (p + pmapRandom(partState)) / numPhi; emap -> ws [0] = cos(phi) * sinTheta; emap -> ws [1] = sin(phi) * sinTheta; /* Flip normal emap -> ws if port and required by its orientation */ setPhotonPortNormal(emap); /* Get surface axes us & vs perpendicular to ws */ do { emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0; emap -> vs [i++] = 1; fcross(emap -> us, emap -> vs, emap -> ws); } while (normalize(emap -> us) < FTINY); fcross(emap -> vs, emap -> ws, emap -> us); /* Get origin */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> src -> sloc [i] + emap -> src -> srad * emap -> ws [i]; /* Get hemisphere axes & aperture */ if (emap -> src -> sflags & SSPOT) { VCOPY(emap -> wh, emap -> src -> sl.s -> aim); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI); } else { VCOPY(emap -> uh, emap -> us); VCOPY(emap -> vh, emap -> vs); VCOPY(emap -> wh, emap -> ws); emap -> cosThetaMax = 0; } } static void sourcePhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on scene cube face for distant light source. Also sets source aperture (solid angle) and sampling hemisphere axes at origin */ { unsigned long i, partsPerDim, partsPerFace; unsigned face; RREAL du, dv; if (emap -> port) { /* Relay to photon port; get origin on its surface */ SRCREC *src = emap -> src; emap -> src = emap -> port; photonOrigin [emap -> src -> so -> otype] (emap); emap -> src = src; } else { /* No ports defined, so get origin on scene cube face (SUBOPTIMAL) */ /* Get current face from partition number */ partsPerDim = 1 / srcsizerat; partsPerFace = sqr(partsPerDim); face = emap -> partitionCnt / partsPerFace; if (!(emap -> partitionCnt % partsPerFace)) { /* Skipped to a new face; get its normal */ emap -> ws [0] = emap -> ws [1] = emap -> ws [2] = 0; emap -> ws [face >> 1] = face & 1 ? 1 : -1; /* Get surface axes us & vs perpendicular to ws */ face >>= 1; emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0; emap -> vs [(face + (emap -> ws [face] > 0 ? 2 : 1)) % 3] = 1; fcross(emap -> us, emap -> vs, emap -> ws); } /* Get jittered offsets within face from partition number (in range [-0.5, 0.5]) */ i = emap -> partitionCnt % partsPerFace; du = (i / partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5; dv = (i % partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5; /* Jittered destination point within partition */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = thescene.cuorg [i] + thescene.cusize * ( 0.5 + du * emap -> us [i] + dv * emap -> vs [i] + 0.5 * emap -> ws [i] ); } /* Get hemisphere axes & aperture */ VCOPY(emap -> wh, emap -> src -> sloc); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); /* Get aperture */ emap -> cosThetaMax = 1 - emap -> src -> ss2 / (2 * PI); emap -> cosThetaMax = min(1, max(-1, emap -> cosThetaMax)); } static void cylPhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on cylindrical light source surface. Also sets source aperture and sampling hemisphere axes at origin */ { int i, cent[3], size[3], parr[2]; FVECT v; cent [0] = cent [1] = cent [2] = 0; size [0] = size [1] = size [2] = emap -> maxPartitions; parr [0] = 0; parr [1] = PMAP_GETPARTITION(emap); if (!skipparts(cent, size, parr, emap -> partitions)) error(CONSISTENCY, "bad source partition in cylPhotonOrigin"); v [SU] = 0; v [SV] = (1 - 2 * pmapRandom(partState)) * (double)size [1] / emap -> maxPartitions; v [SW] = (1 - 2 * pmapRandom(partState)) * (double)size [2] / emap -> maxPartitions; normalize(v); v [SU] = (1 - 2 * pmapRandom(partState)) * (double)size [1] / emap -> maxPartitions; for (i = 0; i < 3; i++) v [i] += (double)cent [i] / emap -> maxPartitions; /* Get surface axes */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> ws [i] = ( v [SV] * emap -> src -> ss [SV][i] + v [SW] * emap -> src -> ss [SW][i] ) / 0.8559; /* Flip normal emap -> ws if port and required by its orientation */ setPhotonPortNormal(emap); normalize(emap -> ws); VCOPY(emap -> us, emap -> src -> ss [SU]); normalize(emap -> us); fcross(emap -> vs, emap -> ws, emap -> us); /* Get origin */ for (i = 0; i < 3; i++) emap -> photonOrg [i] += v [SU] * emap -> src -> ss [SU][i] + emap -> src -> sloc [i]; /* Get hemisphere axes & aperture */ if (emap -> src -> sflags & SSPOT) { VCOPY(emap -> wh, emap -> src -> sl.s -> aim); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI); } else { VCOPY(emap -> uh, emap -> us); VCOPY(emap -> vh, emap -> vs); VCOPY(emap -> wh, emap -> ws); emap -> cosThetaMax = 0; } } /* PHOTON EMISSION ROUTINES ---------------------------------------------- */ static void defaultEmissionFunc (EmissionMap* emap) /* Default behaviour when no emission funcs defined for this source type */ { objerror(emap -> src -> so, INTERNAL, "undefined photon emission function"); } void initPhotonEmissionFuncs () /* Init photonPartition[] and photonOrigin[] dispatch tables */ { int i; for (i = 0; i < NUMOTYPE; i++) photonPartition [i] = photonOrigin [i] = defaultEmissionFunc; photonPartition [OBJ_FACE] = photonPartition [OBJ_RING] = flatPhotonPartition; photonPartition [OBJ_SOURCE] = sourcePhotonPartition; photonPartition [OBJ_SPHERE] = spherePhotonPartition; photonPartition [OBJ_CYLINDER] = cylPhotonPartition; photonOrigin [OBJ_FACE] = photonOrigin [OBJ_RING] = flatPhotonOrigin; photonOrigin [OBJ_SOURCE] = sourcePhotonOrigin; photonOrigin [OBJ_SPHERE] = spherePhotonOrigin; photonOrigin [OBJ_CYLINDER] = cylPhotonOrigin; } void initPhotonEmission (EmissionMap *emap, float numPdfSamples) /* Initialise photon emission from partitioned light source emap -> src; * this involves integrating the flux emitted from the current photon * origin emap -> photonOrg and setting up a PDF to sample the emission * distribution with numPdfSamples samples */ { unsigned i, t, p; double phi, cosTheta, sinTheta, du, dv, dOmega, thetaScale; EmissionSample* sample; const OBJREC* mod = findmaterial(emap -> src -> so); static RAY r; photonOrigin [emap -> src -> so -> otype] (emap); setcolor(emap -> partFlux, 0, 0, 0); emap -> cdf = 0; emap -> numSamples = 0; cosTheta = DOT(emap -> ws, emap -> wh); if (cosTheta <= 0 && sqrt(1-sqr(cosTheta)) <= emap -> cosThetaMax + FTINY) /* Aperture completely below surface; no emission from current origin */ return; if ( mod -> omod == OVOID && !emap -> port && ( cosTheta >= 1 - FTINY || ( emap -> src -> sflags & SDISTANT && acos(cosTheta) + acos(emap -> cosThetaMax) <= 0.5 * PI ) ) ) { /* Source is unmodified and has no port (which requires testing for occlusion), and is either local with normal aligned aperture or distant with aperture above surface --> get flux & PDF via analytical solution */ setcolor( emap -> partFlux, mod -> oargs.farg [0], mod -> oargs.farg [1], mod -> oargs.farg [2] ); /* Multiply radiance by projected Omega * dA to get flux */ scalecolor( emap -> partFlux, PI * cosTheta * (1 - sqr(max(emap -> cosThetaMax, 0))) * emap -> partArea ); } else { /* Source is either modified, has a port, is local with off-normal aperture, or distant with aperture partly below surface --> get flux via numerical integration */ thetaScale = (1 - emap -> cosThetaMax); /* Figga out numba of aperture samples for integration; numTheta / numPhi ratio is 1 / PI */ du = sqrt(pdfSamples * 2 * thetaScale); emap -> numTheta = max(du + 0.5, 1); emap -> numPhi = max(PI * du + 0.5, 1); dOmega = 2 * PI * thetaScale / (emap -> numTheta * emap -> numPhi); thetaScale /= emap -> numTheta; /* Allocate PDF, baby */ sample = emap -> samples = (EmissionSample*)realloc( emap -> samples, sizeof(EmissionSample) * emap -> numTheta * emap -> numPhi ); if (!emap -> samples) error(USER, "can't allocate emission PDF"); VCOPY(r.rorg, emap -> photonOrg); VCOPY(r.rop, emap -> photonOrg); r.rmax = 0; for (t = 0; t < emap -> numTheta; t++) { for (p = 0; p < emap -> numPhi; p++) { /* This uniform mapping handles 0 <= thetaMax <= PI */ cosTheta = 1 - (t + pmapRandom(emitState)) * thetaScale; sinTheta = sqrt(1 - sqr(cosTheta)); phi = 2 * PI * (p + pmapRandom(emitState)) / emap -> numPhi; du = cos(phi) * sinTheta; dv = sin(phi) * sinTheta; rayorigin(&r, PRIMARY, NULL, NULL); for (i = 0; i < 3; i++) r.rdir [i] = ( du * emap -> uh [i] + dv * emap -> vh [i] + cosTheta * emap -> wh [i] ); /* Sample behind surface? */ VCOPY(r.ron, emap -> ws); if ((r.rod = DOT(r.rdir, r.ron)) <= 0) continue; /* Get radiance emitted in this direction; to get flux we multiply by cos(theta_surface), dOmega, and dA. Ray is directed towards light source for raytexture(). */ if (!(emap -> src -> sflags & SDISTANT)) for (i = 0; i < 3; i++) r.rdir [i] = -r.rdir [i]; /* Port occluded in this direction? */ if (emap -> port && localhit(&r, &thescene)) continue; raytexture(&r, mod -> omod); setcolor( r.rcol, mod -> oargs.farg [0], mod -> oargs.farg [1], mod -> oargs.farg [2] ); multcolor(r.rcol, r.pcol); /* Multiply by cos(theta_surface) */ scalecolor(r.rcol, r.rod); /* Add PDF sample if nonzero; importance info for photon emission * could go here... */ if (colorAvg(r.rcol)) { copycolor(sample -> pdf, r.rcol); sample -> cdf = emap -> cdf += colorAvg(sample -> pdf); sample -> theta = t; sample++ -> phi = p; emap -> numSamples++; addcolor(emap -> partFlux, r.rcol); } } } /* Multiply by dOmega * dA */ scalecolor(emap -> partFlux, dOmega * emap -> partArea); } } #define vomitPhoton emitPhoton #define bluarrrghPhoton vomitPhoton void emitPhoton (const EmissionMap* emap, RAY* ray) /* Emit photon from current partition emap -> partitionCnt based on emission distribution. Returns new photon ray. */ { unsigned long i, lo, hi; const EmissionSample* sample = emap -> samples; RREAL du, dv, cosTheta, cosThetaSqr, sinTheta, phi; const OBJREC* mod = findmaterial(emap -> src -> so); /* Choose a new origin within current partition for every emitted photon to break up clustering artifacts */ photonOrigin [emap -> src -> so -> otype] ((EmissionMap*)emap); /* If we have a local glow source with a maximum radius, then restrict our photon to the specified distance, otherwise we set the limit imposed by photonMaxDist (or no limit if 0) */ if ( mod -> otype == MAT_GLOW && !(emap -> src -> sflags & SDISTANT) && mod -> oargs.farg[3] > FTINY ) ray -> rmax = mod -> oargs.farg[3]; else ray -> rmax = photonMaxDist; rayorigin(ray, PRIMARY, NULL, NULL); if (!emap -> numSamples) { /* Source is unmodified and has no port, and either local with normal aligned aperture, or distant with aperture above surface --> use cosine weighted distribution */ cosThetaSqr = (1 - pmapRandom(emitState) * (1 - sqr(max(emap -> cosThetaMax, 0))) ); cosTheta = sqrt(cosThetaSqr); sinTheta = sqrt(1 - cosThetaSqr); phi = 2 * PI * pmapRandom(emitState); setcolor( ray -> rcol, mod -> oargs.farg [0], mod -> oargs.farg [1], mod -> oargs.farg [2] ); } else { /* Source is either modified, has a port, is local with off-normal aperture, or distant with aperture partly below surface --> choose direction from constructed cumulative distribution function with Monte Carlo inversion using binary search. */ du = pmapRandom(emitState) * emap -> cdf; lo = 1; hi = emap -> numSamples; while (hi > lo) { i = (lo + hi) >> 1; sample = emap -> samples + i - 1; if (sample -> cdf >= du) hi = i; if (sample -> cdf < du) lo = i + 1; } /* This is a uniform mapping, mon */ cosTheta = (1 - (sample -> theta + pmapRandom(emitState)) * (1 - emap -> cosThetaMax) / emap -> numTheta ); sinTheta = sqrt(1 - sqr(cosTheta)); phi = 2 * PI * (sample -> phi + pmapRandom(emitState)) / emap -> numPhi; copycolor(ray -> rcol, sample -> pdf); } /* Normalize photon flux so that average over RGB is 1 */ colorNorm(ray -> rcol); VCOPY(ray -> rorg, emap -> photonOrg); du = cos(phi) * sinTheta; dv = sin(phi) * sinTheta; for (i = 0; i < 3; i++) ray -> rdir [i] = du * emap -> uh [i] + dv * emap -> vh [i] + cosTheta * emap -> wh [i]; if (emap -> src -> sflags & SDISTANT) /* Distant source; reverse ray direction to point into the scene. */ for (i = 0; i < 3; i++) ray -> rdir [i] = -ray -> rdir [i]; if (emap -> port) /* Photon emitted from port; move origin just behind port so it will be scattered */ for (i = 0; i < 3; i++) ray -> rorg [i] -= 2 * FTINY * ray -> rdir [i]; /* Assign emitting light source index */ ray -> rsrc = emap -> src - source; } /* SOURCE CONTRIBS FROM DIRECT / VOLUME PHOTONS -------------------------- */ void multDirectPmap (RAY *r) /* Factor irradiance from direct photons into r -> rcol; interface to * direct() */ { COLOR photonIrrad; - - /* Lookup direct photon irradiance */ - (directPmap -> lookup)(directPmap, r, photonIrrad); + +#ifdef PMAP_PHOTONFLOW + if (photonFlow) + if (volumePhotonMapping) + /* Photon flow mode: evaluate volume photon map as lightfield + * (including direct component), ignoring participating medium. */ + volumePhotonSphIrrad(volumePmap, r, photonIrrad); + else + error(USER, "photon flow mode requires a volume photon map"); + else +#endif + /* Lookup direct photon irradiance */ + (directPmap -> lookup)(directPmap, r, photonIrrad); /* Multiply with coefficient in ray */ multcolor(r -> rcol, photonIrrad); return; } void inscatterVolumePmap (RAY *r, COLOR inscatter) /* Add inscattering from volume photon map; interface to srcscatter() */ { /* Add ambient in-scattering via lookup callback */ (volumePmap -> lookup)(volumePmap, r, inscatter); } diff --git a/pmapsrc.h b/pmapsrc.h index 79cecc3..f26ab7f 100644 --- a/pmapsrc.h +++ b/pmapsrc.h @@ -1,105 +1,107 @@ /* RCSid $Id: pmapsrc.h,v 2.7 2020/08/07 01:21:13 rschregle Exp $ */ /* ====================================================================== Photon map support routines for emission from light sources Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation (DFG) - under the FARESYS project. + 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). + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, - supported by the JSPS KAKENHI Grant Number JP19KK0115. + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapsrc.h,v 2.7 2020/08/07 01:21:13 rschregle Exp $ */ #ifndef PMAPSRC_H #define PMAPSRC_H #include "ray.h" #include "source.h" /* Data structures for photon emission */ typedef struct { unsigned theta, phi; /* Direction indices */ COLOR pdf; /* Probability of emission in this direction per RGB */ float cdf; /* Cumulative probability up to this sample */ } EmissionSample; typedef struct { unsigned numTheta, numPhi; /* Num divisions in theta & phi */ RREAL cosThetaMax; /* cos(source aperture) */ FVECT uh, vh, wh; /* Emission aperture axes at origin, w is normal*/ FVECT us, vs, ws; /* Source surface axes at origin, w is normal */ FVECT photonOrg; /* Current photon origin */ EmissionSample *samples; /* PDF samples for photon emission directions */ unsigned long numPartitions; /* Number of surface partitions */ RREAL partArea; /* Area covered by each partition */ SRCREC *src, *port; /* Current source and port */ unsigned long partitionCnt, /* Current surface partition */ maxPartitions, /* Max allocated partitions */ numSamples; /* Number of PDF samples */ unsigned char* partitions; /* Source surface partitions */ COLOR partFlux; /* Total flux emitted by partition */ float cdf; /* Cumulative probability */ } EmissionMap; /* Photon port flags (orientation relative to surface normal): * Forward, backward, both (bidirectional). */ #define PMAP_PORTFWD 1 #define PMAP_PORTBWD 2 #define PMAP_PORTBI (PMAP_PORTFWD | PMAP_PORTBWD) /* Photon ports for emission from geometry en lieu of light sources */ extern char *photonPortList [MAXSET + 1]; extern SRCREC *photonPorts; extern unsigned numPhotonPorts; /* Dispatch tables for partitioning light source surfaces and generating an origin for photon emission */ extern void (*photonPartition []) (EmissionMap*); extern void (*photonOrigin []) (EmissionMap*); void getPhotonPorts (char **portList); /* Find geometry declared as photon ports from modifiers in portList */ void initPhotonEmissionFuncs (); /* Init photonPartition[] and photonOrigin[] dispatch tables with source specific emission routines */ void initPhotonEmission (EmissionMap *emap, float numPdfSamples); /* Initialize photon emission from partitioned light source emap -> src; * this involves integrating the flux emitted from the current photon * origin emap -> photonOrg and setting up a PDF to sample the emission * distribution with numPdfSamples samples */ void emitPhoton (const EmissionMap*, RAY* ray); /* Emit photon from current source partition based on emission distribution and return new photon ray */ void multDirectPmap (RAY *r); /* Factor irradiance from direct photons into r -> rcol; interface to * direct() */ void inscatterVolumePmap (RAY *r, COLOR inscatter); /* Add inscattering from volume photon map; interface to srcscatter() */ #endif diff --git a/pmutil.c b/pmutil.c index 468444f..938a6ba 100644 --- a/pmutil.c +++ b/pmutil.c @@ -1,277 +1,115 @@ #ifndef lint static const char RCSid[] = "$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $"; #endif /* ====================================================================== - Photon map utilities + Auxiliary photon map utilities Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF, #147053) ====================================================================== $Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $ */ #include "pmap.h" #include "pmapio.h" -#include "pmapbias.h" -#include "otypes.h" #include extern char *octname; -/* Photon map lookup functions per type */ -void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = { - photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity, - photonDensity, photonDensity -}; - - - void colorNorm (COLOR c) /* Normalise colour channels to average of 1 */ { const float avg = colorAvg(c); if (!avg) return; c [0] /= avg; c [1] /= avg; c [2] /= avg; } - void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm) { unsigned t; struct stat octstat, pmstat; PhotonMap *pm; PhotonMapType type; for (t = 0; t < NUM_PMAP_TYPES; t++) if (setPmapParam(&pm, parm + t)) { /* Check if photon map newer than octree */ if (pm -> fileName && octname && !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) && octstat.st_mtime > pmstat.st_mtime) { sprintf(errmsg, "photon map in file %s may be stale", pm -> fileName); error(USER, errmsg); } /* Load photon map from file and get its type */ if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE) error(USER, "failed loading photon map"); /* Assign to appropriate photon map type (deleting previously * loaded photon map of same type if necessary) */ if (pmaps [type]) { sprintf(errmsg, "multiple %s photon maps, dropping previous", pmapName [type]); error(WARNING, errmsg); deletePhotons(pmaps [type]); free(pmaps [type]); } pmaps [type] = pm; /* Check for valid density estimate bandwidths */ if ((pm -> minGather > 1 || pm -> maxGather > 1) && (type == PMAP_TYPE_PRECOMP)) { /* Force bwidth to 1 for precomputed pmap */ error(WARNING, "ignoring bandwidth for precomp photon map"); pm -> minGather = pm -> maxGather = 1; } if ((pm -> maxGather > pm -> minGather) && (type == PMAP_TYPE_VOLUME)) { /* Biascomp for volume pmaps (see volumePhotonDensity() below) is considered redundant, and there's probably no point in recovering by using the lower bandwidth, since it's probably not what the user wants, so bail out. */ sprintf(errmsg, "bias compensation is not available with %s photon maps", pmapName [type]); error(USER, errmsg); } if (pm -> maxGather > pm -> numPhotons) { error(WARNING, "adjusting density estimate bandwidth"); pm -> minGather = pm -> maxGather = pm -> numPhotons; } } } void cleanUpPmaps (PhotonMap **pmaps) { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) { deletePhotons(pmaps [t]); free(pmaps [t]); } } } - - - -void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) -/* Photon density estimate. 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 : "", - 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) -/* Returns precomputed photon density estimate 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) -/* Photon volume density estimate. Returns 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; - -#if 0 - /* Volume biascomp disabled (probably redundant) */ - if (pmap -> minGather == pmap -> maxGather) -#endif - { - /* No bias compensation. Just do a plain vanilla estimate */ - 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))); - return; - } -#if 0 - else - /* Apply bias compensation to density estimate */ - volumeBiasComp(pmap, ray, irrad); -#endif -}