diff --git a/Rmakefile b/Rmakefile index 2cedd97..d5cbe20 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,470 +1,471 @@ # RCSid: $Id: Rmakefile,v 2.87 2021/01/19 23:31:47 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 pmapdens.o \ pmapbias.o pmapparm.o pmapcontrib.o pmcontrib2.o pmapamb.o pmapray.o \ pmapopt.o pmapdiag.o pmaptype.o morton.o oococt.o oocsort.o \ - oocbuild.o oocnn.o ooccache.o pmutil.o pmapfilt.o mrgbe.o + oocbuild.o oocnn.o ooccache.o pmutil.o pmaproi.o pmapfilt.o mrgbe.o PMSRC = pmap.c pmapsrc.c pmapmat.c pmaprand.c pmapio.c pmapdata.c pmapdens.c \ pmapbias.c pmapparm.c pmapcontrib.c pmcontrib2.c pmapamb.c pmapray.c \ pmapopt.c pmapdiag.c pmaptype.c morton.c oococt.c oocsort.c \ - oocbuild.c oocnn.c ooccache.c pmutil.c pmapfilt.c pmapkdt.c \ - pmaptkdt.c pmapooc.c mrgbe.c + oocbuild.c oocnn.c ooccache.c pmutil.c pmapfilt.c pmaproi.c \ + pmapkdt.c pmaptkdt.c pmapooc.c mrgbe.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) 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 || rm Version.c # # 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 pmapdata.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') +# TODO: Update after adding pmaproi module! # 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 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 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 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 \ pmapkdt.h pmapkdt.c pmapooc.h 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 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 pmapdiag.h ../common/platform.h ../common/resolu.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c 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 morton.h oocbuild.o: oocbuild.c oococt.h morton.h ../common/fvect.h ooccache.h \ oocsort.h oococt.o: oococt.c oococt.h morton.h ../common/fvect.h ooccache.h \ ../common/rtio.h oocnn.o: oocnn.c oocnn.h oococt.h morton.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 ../common/random.h source.h otspecial.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c 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/paths.h \ ../common/lookup.h pmapbias.h ../common/otypes.h \ pmapkdt.h pmapkdt.c pmaptkdt.h pmaptkdt.c pmapooc.h pmapooc.c mkpmap.o: 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.h pmapkdt.c pmaptkdt.h pmaptkdt.c pmapooc.h pmapooc.c 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 pmapio.h ../common/rtio.h \ ../common/resolu.h ../common/random.h \ pmapkdt.h pmapkdt.c pmaptkdt.h pmaptkdt.c pmapooc.h pmapooc.c 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 pmapdens.h pmap.h pmaprand.h pmapmat.h \ ../common/otypes.h ../common/random.h source.h otspecial.h \ pmapkdt.h pmapkdt.c pmaptkdt.h pmaptkdt.c pmapooc.h pmapooc.c 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 pmapmat.h pmapsrc.h source.h pmaprand.h pmapio.h \ pmapdens.h pmapbias.h pmapdiag.h ../common/platform.h ../common/otypes.h \ pmapkdt.h pmapkdt.c pmaptkdt.h pmaptkdt.c pmapooc.h pmapooc.c 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 pmap.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c 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 pmaprand.h ../common/otypes.h data.h func.h \ ../common/calcomp.h ../common/bsdf.h ../common/ccolor.h \ ../common/platform.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c 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/paths.h ../common/lookup.h pmaprand.h \ ../common/otypes.h otspecial.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c 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 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 pmapkdt.c pmapooc.h pmapooc.c -pmutil.o: pmutil.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ +pmutil.o: pmutil.h 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 pmapio.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c diff --git a/Version.c b/Version.c index 18eb02c..e57da59 100644 --- a/Version.c +++ b/Version.c @@ -1,5 +1,5 @@ /* * This file was created automatically during make. */ -char VersionID[]="RADIANCE 5.4a lastmod Wed 30 Jun 01:11:05 CEST 2021 by u-no-hoo on ovo"; +char VersionID[]="RADIANCE 5.4a lastmod Tue 31 Aug 02:36:52 CEST 2021 by u-no-hoo on ovo"; diff --git a/mkpmap.c b/mkpmap.c index feb899f..8f1b374 100644 --- a/mkpmap.c +++ b/mkpmap.c @@ -1,947 +1,980 @@ #ifndef lint static const char RCSid[] = "$Id: mkpmap.c,v 2.11 2021/04/14 11:26:25 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 LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: mkpmap.c,v 2.11 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmapcontrib.h" #include "pmaprand.h" +#include "pmaproi.h" +#include "pmapio.h" #include "paths.h" #include "ambient.h" #include "resolu.h" #include "source.h" #include #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 dimensions */ 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 */ int nproc = 1; /* number of parallel processes */ #ifdef EVALDRC_HACK char *angsrcfile = NULL; /* angular source file for EvalDRC */ #endif /* Contribution photon map stuff: light source modifier lookup table and associated cleanup func */ static void chickenMcFreeNuggetz(void *contrib) { epfree((*(MODCONT *)contrib).binv); free(contrib); } LUTAB contribTab = LU_SINIT(NULL, chickenMcFreeNuggetz); #if 0 char srcMod [MAXMODLIST]; #endif unsigned numContribs = 0; /* The variables below interface with the RADIANCE suite merely as dummies to resolve external references and make the linker happy; most of them are unused. */ 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; char *shm_boundary = NULL, *ambfile = NULL, RCCONTEXT [] = "RC."; /* Evaluation context for contrib pmap */ 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("-aph mod\t\t\t\t# polyhedral region of interest"); + puts("-apH file\t\t\t\t# polyhedral region of interest"); puts("-api xmin ymin zmin xmax ymax zmax" "\t# rectangular region of interest" ); puts("-apI xpos ypos zpos radius\t\t# spherical region of interest"); #endif puts("-apc file nPhotons\t\t\t# caustic photon map"); puts("-apC file nPhotons bwidth compression" "\t# precomputed contribution photon map" ); puts("-apd file nPhotons\t\t\t# direct photon map"); puts("-apg file nPhotons\t\t\t# global photon map"); puts("-app file nPhotons bwidth\t\t# precomputed global photon map"); /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC puts("-apt file nPhotons velocity\t\t# transient photon map"); #endif #ifdef PMAP_PHOTONFLOW /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC puts("-apT file nPhotons velocity\t\t# transient light flow photon map"); #endif puts("-apV file nPhotons\t\t\t# light flow photon map"); #endif puts("-apv file nPhotons\t\t\t# volume 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"); + puts("-bn\t\t\t\t\t# contribution binning resolution"); 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); + puts("-e expr\t\t\t\t\t# init expression for contrib binning"); + printf("-ef %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", photonMaxPathDist); 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 + puts("-p\t\t\t\t\t# per-modifier contrib binning params"); 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" ); #if 0 /* TODO: Do contrib/coeffs here or in rcontrib? */ printf( "-V%c\t\t\t\t# precompute %s\n", contrib ? '+' : '-', contrib ? "contributions" : "coefficients" ); #endif } #if 0 /* Disabled for now; will attempt to extract bins from photon map header */ static char *appendParams (char *accParams, char **params, int n) /* Accumulate n next parameter strings from params to string accParams. Returns new accmulated string. */ { int i; if (params && n) for (i = 0; i < n; i++) if (params [i] && strlen(params [i])) /* Reallocate and factor in delimiter (space) and terminator */ if (accParams = realloc( accParams, strlen(accParams) + strlen(params [i]) + 2 )) { strcat(accParams, " "); strcat(accParams, params [i]); } else error( SYSTEM, "cannot append binning parameters for contrib photon map" ); return accParams; } #endif 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; \ } /* Check target number of photons with optional multiplier suffix begins with a digit */ #define checkNumPhotons(i) if (!isdigit(argv [i][0])) \ if (argv [i][0] == '-') \ goto badopt; \ else { \ sprintf(errmsg, "invalid number of photons %s", argv [i]); \ error(USER, errmsg); \ } int loadflags = IO_CHECK | IO_SCENE | IO_TREE | IO_BOUNDS, - rval, i, j, n, binCnt = 0; + rval, i, j, n, binRes = 0; char **portLp = photonPortList, **sensLp = photonSensorList, - **amblp = NULL, sbuf [MAXSTR], portFlags [2] = "\0\0", - sensFlags [2] = "\0\0", *binVal = NULL, *binParm = NULL; + **roiModLp = pmapROImodList, **amblp = NULL, sbuf [MAXSTR], + portFlags [2] = "\0\0", sensFlags [2] = "\0\0", *binParm = NULL; struct stat pmstat; /* Global program name */ progname = fixargv0(argv [0]); /* Initialize object types */ initotypes(); /* Initialize cal function routines for contrib photon binning */ initfunc(); setcontext(RCCONTEXT); + + /* Compile default orientation variables for contribution binning */ + scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); #if defined(_WIN32) || defined(_WIN64) /* Increase file limit to maximum */ /* XXX: DO WE NEED THIS FOR CONTRIB PHOTONS? _setmaxstdio(2048); */ #endif /* 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 (isupper(argv [i][2])) { /* 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 (isupper(argv [i][2])) { /* 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 */ checkNumPhotons(i + 2); 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 */ checkNumPhotons(i + 2); 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 */ checkNumPhotons(i + 2); check(4, "ss"); causticPmapParams.fileName = argv [++i]; causticPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!causticPmapParams.distribTarget) goto badopt; break; case 'v': /* Volume photon map */ checkNumPhotons(i + 2); check(4, "ss"); volumePmapParams.fileName = argv [++i]; volumePmapParams.distribTarget = parseMultiplier(argv [++i]); if (!volumePmapParams.distribTarget) goto badopt; break; case 'd': /* Direct photon map */ checkNumPhotons(i + 2); check(4, "ss"); directPmapParams.fileName = argv [++i]; directPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!directPmapParams.distribTarget) goto badopt; break; case 'C': /* Precomputed contribution photon map */ checkNumPhotons(i + 2); check(4, "ssif"); contribPmapParams.fileName = argv [++i]; contribPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!contribPmapParams.distribTarget) goto badopt; contribPmapParams.minGather = contribPmapParams.maxGather = atoi(argv [++i]); if (!contribPmapParams.maxGather) goto badopt; compressRatio = atof(argv [++i]); if (compressRatio < FTINY || compressRatio > 1 - FTINY ) error(USER, "contribution photon compression ratio must " "be in range ]0, 1[" ); break; /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC case 't': /* Transient photon map */ checkNumPhotons(i + 2); check(4, "ssf"); transientPmapParams.fileName = argv [++i]; transientPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!transientPmapParams.distribTarget) goto badopt; photonVelocity = atof(argv [++i]); if (photonVelocity < FTINY) error(USER, "can't halt or reverse time! " "We're not Steven bloody Hawking, you know!" ); break; #endif #ifdef PMAP_PHOTONFLOW /* Light flow is a variant of volume photon map */ case 'V': /* Light flow photon map */ checkNumPhotons(i + 2); check(4, "ss"); lightFlowParams.fileName = argv [++i]; lightFlowParams.distribTarget = parseMultiplier(argv [++i]); if (!lightFlowParams.distribTarget) goto badopt; /* Set zero absorption and forward scattering for global mist; extinction up to user of local mist defined in octree */ /* setcolor(cextinction, 1, 1, 1); */ setcolor(salbedo, 1, 1, 1); seccg = 1; break; /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC case 'T': /* Transient light flow photon map */ checkNumPhotons(i + 2); check(4, "ssf"); transLightFlowParams.fileName = argv [++i]; transLightFlowParams.distribTarget = parseMultiplier(argv [++i]); if (!transLightFlowParams.distribTarget) goto badopt; photonVelocity = atof(argv [++i]); if (photonVelocity < FTINY) error(USER, "can't halt or reverse time! " "We're not Steven bloody Hawking, you know!" ); /* Set zero absorption and forward scattering for global mist; extinction is left up to user or defined as local mist boundary in octree */ /* setcolor(cextinction, 1, 1, 1); */ setcolor(salbedo, 1, 1, 1); seccg = 1; break; #endif #endif 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 'h': + case 'H': + /* Add polyhedral region(s) of interest by + modifier(s) */ + check(4, "s"); + + if (isupper(argv[i][3])) { + /* Add ROI modifiers from file */ + rval = wordfile(roiModLp, + MAXSET - (roiModLp - pmapROImodList), + getpath(argv [++i], getrlibpath(), R_OK) + ); + if (rval < 0) { + sprintf(errmsg, + "cannot open ROI modifier file %s", + argv [i] + ); + error(SYSTEM, errmsg); + } + } + else { + /* Add modifier string to list, mark end of list + * with NULL */ + *roiModLp++ = savqstr(argv [++i]); + *roiModLp = NULL; + } + break; + case 'i': /* Add rectangular region of interest */ case 'I': /* Add spherical region of interest */ check(4, isupper(argv [j=i][3]) ? "ffff" : "ffffff"); n = pmapNumROI; pmapROI = realloc(pmapROI, ++pmapNumROI * sizeof(PhotonMapROI) ); if (!pmapROI) error(SYSTEM, "failed to allocate ROI"); pmapROI [n].pos [0] = atof(argv [++i]); pmapROI [n].pos [1] = atof(argv [++i]); pmapROI [n].pos [2] = atof(argv [++i]); pmapROI [n].siz [0] = atof(argv [++i]); if (isupper(argv [j][3])) { /* Spherical ROI; radius^2 */ pmapROI [n].siz [0] *= pmapROI [n].siz [0]; PMAP_ROI_SETSPHERE(pmapROI + n); if (pmapROI [n].siz [0] <= FTINY) error(USER, "region of interest has invalid radius" ); } else { /* Rectangular ROI */ pmapROI [n].siz [1] = atof(argv [++i]); pmapROI [n].siz [2] = atof(argv [++i]); for (j = 0; j < 3; j++) { /* Pos at rectangle centre, siz symmetric */ pmapROI [n].pos [j] = 0.5 * ( pmapROI [n].pos [j] + pmapROI [n].siz [j] ); pmapROI [n].siz [j] = fabs( pmapROI [n].siz [j] - pmapROI [n].pos [j] ); if (pmapROI [n].siz [j] <= FTINY) error(USER, "region of interest has invalid size" ); } } break; #endif case 'P': /* Photon precomputation ratio */ check(4, "f"); finalGather = atof(argv [++i]); if (finalGather <= 0 || finalGather > 1) error(USER, "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 (isupper(argv [i][3])) { /* 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 orientation flags to all modifier * strings; 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 { /* HACK: append flags to modifier string, add to - * port list and mark of end list with NULL */ + * port list and mark end of 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 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_SENSFWD, PMAP_SENSBWD, sensFlags [0] ); if (isupper(argv[i][3])) { /* 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); } /* HACK: append orientation flags to all modifier * strings; note this requires reallocation */ for (; rval--; sensLp++) { j = strlen(*sensLp); if (!(*sensLp = realloc(*sensLp, j + 2))) { sprintf(errmsg, "cannot allocate antimatter sensor " "modifiers from file %s", argv [i] ); error(SYSTEM, errmsg); } strcat(*sensLp, sensFlags); } } else { /* HACK: append flags to modifier string, add to * sensor list and mark of end list with NULL */ strcpy(sbuf, argv [++i]); strcat(sbuf, sensFlags); *sensLp++ = savqstr(sbuf); *sensLp = NULL; } break; default: goto badopt; } break; default: goto badopt; } break; case 'b': switch (argv [i][2]) { case 'v': /* Back face visibility */ check_bool(3, backvis); break; case 'n': - /* Number of bins for precomputed contrib. photon map. */ - check(3, "s"); - /* binParams = appendParams(binParams, argv + i, 2); */ - binCnt = (int)(eval(argv [++i]) + .5); + /* Bin resolution for precomp contrib pmap */ + check(3, "i"); + binRes = atoi(argv [++i]); + if (binRes <= 1 || binRes > PMAP_CONTRIB_MAXBINRES) + goto badopt; break; default: 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': switch (argv [i][2]) { case 'f': /* Diagnostics file */ - check(2, "s"); + check(3, "s"); diagFile = argv [++i]; break; default: - /* Functional expression for precomputed contrib. pmap */ + /* Functional expression for precomputed contrib. pmap, + e.g. to set constants for disk2square.cal */ check(2, "s"); - /* binParams = appendParams(binParams, argv + i, 2); */ scompile(argv [++i], NULL, 0); } break; case 'f': switch (argv [i][2]) { case 'o': /* Force overwrite */ check_bool(3, clobber); break; default: - /* Function file for precomputed contribution photon map */ - check(2, "s"); - /* binParams = appendParams(binParams, argv + i, 2); */ - loadfunc(argv[++i]); + goto badopt; } break; #ifdef PMAP_EKSPERTZ case 'l': /* Limits */ switch (argv [i][2]) { case 'd': /* Limit photon path distance */ check(3, "f"); photonMaxPathDist = atof(argv [++i]); if (photonMaxPathDist <= 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': switch (argv[i][2]) { case 'e': /* Mist 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': /* Mist albedo */ check(3, "fff"); setcolor(salbedo, atof(argv [i + 1]), atof(argv [i + 2]), atof(argv [i + 3]) ); i += 3; break; case 'g': /* Mist scattering eccentricity */ check(3, "f"); seccg = atof(argv [++i]); break; default: /* Modifier name for precomputed contrib. photon map */ check(2, "s"); - /* binParams = appendParams(binParams, argv + i, 2); */ addContribModifier(&contribTab, &numContribs, - argv [++i], binParm, binVal, binCnt + argv [++i], binParm, NULL, PMAP_CONTRIB_SCBINS(binRes) ); break; } break; case 'M': /* Modifier file for precomputed contribution photon map */ check(2, "s"); - /* binParams = appendParams(binParams, argv + i, 2); */ addContribModfile(&contribTab, &numContribs, - argv [++i], binParm, binVal, binCnt + argv [++i], binParm, NULL, PMAP_CONTRIB_SCBINS(binRes) ); 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 'p': /* Parameter setting(s) for precomputed contrib. photon map */ check(2, "s"); - /* binParams = appendParams(binParams, argv + i, 2); */ set_eparams(binParm = argv [++i]); break; case 't': /* Timer */ check(2, "i"); photonRepTime = atoi(argv [++i]); break; case 'v': /* Verbosity */ check_bool(2, verbose); break; #ifdef EVALDRC_HACK case 'A': /* Capt. B's (now historic) angular source file hack */ 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, &contribTab, numContribs, nproc); else distribPhotons(photonMaps, nproc); /* TODO: Where do we save binParams and wavelet coeffs? */ /* 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/pmap.c b/pmap.c index bf4483d..79ef24d 100644 --- a/pmap.c +++ b/pmap.c @@ -1,934 +1,927 @@ #ifndef lint static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 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, "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.18 2021/02/19 02:10:35 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" +#include "pmaproi.h" #include "pmapio.h" #include "pmapdens.h" #include "pmapbias.h" #include "pmapdiag.h" +#include "pmutil.h" #include "otypes.h" #include "otspecial.h" #include #if NIX #include #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, rmax0; const RREAL cext = colorAvg(ray -> cext), invCext = 1 / cext, albedo = colorAvg(ray -> albedo), invAlbedo = 1 / albedo, gecc = ray -> gecc, gecc2 = sqr(gecc); FVECT u, v; COLOR cvext; /* Save incoming ray's max distance; this is required since rmax is set to the mean free distance for each path segment (see below) */ rmax0 = ray -> rmax; /* Mean free distance until interaction with medium */ ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { /* No geometry hit within mean free dist ray->rmax; photon interacts with medium */ if (!incube(&thescene, ray -> rop)) { /* Terminate photon if it has leaked from the scene */ #ifdef DEBUG_PMAP_LEAKED fprintf( stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); #endif return 0; } /* RGB attenuation over mean free distance */ 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); /* Set ray origin for next interaction */ VCOPY(ray -> rorg, ray -> rop); /* Update ray's path length; terminate photon if maximum reached. See also photonRay(). */ rmax0 -= ray -> rot; if (photonMaxPathDist > 0 && rmax0 < 0) return 0; /* Store volume photon if nonzero albedo; this also accounts for direct inscattering from light sources */ if (albedo > FTINY) { /* Set ray's path length to store with photon */ ray -> rmax = rmax0; #ifdef PMAP_PHOTONFLOW if (lightFlowPhotonMapping) { /* Temporarily correct normalised flux for lightflow photons based on extinction (= mean num photons deposited per unit path length). */ scalecolor(ray -> rcol, invCext); newPhoton(lightFlowPmap, ray); /* Undo correction for next iteration after storing photon */ scalecolor(ray -> rcol, cext); } else #endif newPhoton(volumePmap, ray); } /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); scalecolor(ray -> rcol, invAlbedo); if (pmapRandom(rouletteState) < 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); /* Restore outgoing ray's max distance */ ray -> rmax = rmax0; 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; /* NOTE: localhit() limits intersections to the aft plane's distance ray->rmax. This mechanism is (ab)used here to limit the photon path length by initialising ray->rmax to photonMaxPathDist; see emitPhoton(). With unlimited path length (photonMaxPathDist=0), ray->rmax becomes negative (see photonRay()), and the aft plane has no effect. */ if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; if (port && ray -> ro != port) { /* Terminate photon if emitted from port without intersecting it */ #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; /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ #define PMAP_CNTUPDATE 0xffL 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, lastPhotonCnt; pid_t procPids [PMAP_MAXPROC]; 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); + /* Get polyhedral regions of interest */ + getPolyROIs(pmapROImodList); + #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 */ /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); 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 { /* 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++ ) { 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; } /* PHOTON DISTRIBUTION LOOP */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; 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++ ) { 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++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon based on PDF, skip if rejected. */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ if (!emitPhoton(&emap, &photonRay)) continue; /* Set subprocess ID in photonRay; extends path ID. Used by newPhoton() as photon attributes. */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter*/ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update global counters shared with siblings in intervals to reduce overhead and contention */ lastPhotonCnt.numEmitted = photonCnt -> numEmitted; photonCnt -> numEmitted += PMAP_CNTUPDATE; lastPhotonCnt.numComplete = photonCnt -> numComplete; photonCnt -> numComplete += PMAP_CNTUPDATE; /* Check for overflow using previous values */ if (photonCnt -> numEmitted < lastPhotonCnt.numEmitted || photonCnt -> numComplete < lastPhotonCnt.numComplete ) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt.numEmitted, photonCnt -> numEmitted ); error(INTERNAL, errmsg); } /* Update global photon count for each pmap */ for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) { lastPhotonCnt.numPhotons [t] = photonCnt -> numPhotons [t]; /* Differential increment using local counter lastNumPhotons */ photonCnt -> numPhotons [t] += pmaps [t] -> numPhotons - lastNumPhotons [t]; lastNumPhotons [t] = pmaps [t] -> numPhotons; /* Check for photon counter overflow (this * could only happen before an emission * counter overflow if the scene has an * absurdly high albedo and/or very dense * geometry) */ if (photonCnt -> numPhotons [t] < lastPhotonCnt.numPhotons [t] ) { sprintf(errmsg, "photon counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt.numPhotons [t], photonCnt -> numPhotons [t] ); error(INTERNAL, errmsg); } } } } } #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 { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotons"); else /* Save child process IDs */ procPids [proc] = pid; } } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses to complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); 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 1e85bb3..c3d6f19 100644 --- a/pmap.h +++ b/pmap.h @@ -1,112 +1,106 @@ /* 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, "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: Combine type-specific properties of pmap types as binary flags, e.g. PMAP_TYPE_TRANSLIGHTFLOW = PMAP_TYPE_VOLUME | PMAP_TYPE_TRANSIENT | PMAP_TYPE_LIGHTFLOW ??? */ #ifndef PMAP_H #define PMAP_H #ifndef NIX #if defined(_WIN32) || defined(_WIN64) #define NIX 0 #else #define NIX 1 #endif #endif #include "pmapparm.h" #include "pmapdata.h" + #include "pmutil.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 map (also see pmapdata.h) */ #define photonMapping (\ globalPmap || preCompPmap || causticPmap || contribPmap || \ transientPmap \ ) #define causticPhotonMapping (causticPmap) #define directPhotonMapping (directPmap || transientPhotonMapping) #define volumePhotonMapping (volumePmap) /* #define contribPhotonMapping (contribPmap && contribPmap -> srcContrib) */ #define contribPhotonMapping (contribPmap) #define transientPhotonMapping (transientPmap) #ifdef PMAP_PHOTONFLOW #define lightFlowPhotonMapping (lightFlowPmap) #endif extern void (*pmapLookup [])(PhotonMap*, RAY*, COLOR); /* Photon map lookup functions per type */ - + +#if 0 /* Moved to pmapio */ 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 */ +#endif +#if 0 /* Moved to pmutil */ void cleanUpPmaps (PhotonMap **pmaps); /* Trash all photon maps after processing is complete */ +#endif 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. */ +#if 0 /* Moved to pmutils */ void colorNorm (COLOR); /* Normalise colour channels to average of 1 */ +#endif #endif diff --git a/pmapcontrib.c b/pmapcontrib.c index 7976373..d29c682 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1150 +1,1208 @@ #ifndef lint static const char RCSid[] = "$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to mkpmap. 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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $ */ #include "pmapcontrib.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" +#include "pmaproi.h" #include "pmapio.h" #include "pmapdiag.h" +#include "wavelet2.h" #include "otypes.h" #include "otspecial.h" #if NIX #include #include #endif +extern void SDdisk2square(double sq[2], double diskx, double disky); -MODCONT *addContribModifier(LUTAB *contribTab, unsigned *numContribs, + +MODCONT *addContribModifier ( + LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, char *binVal, int binCnt ) /* Add light source modifier mod to contribution lookup table contribTab, and update numContribs. Return initialised contribution data for this modifier. */ { LUENT *lutEntry = lu_find(contribTab, mod); MODCONT *contrib; EPNODE *eBinVal; if (lutEntry -> data) { /* Reject duplicate modifiers */ sprintf(errmsg, "duplicate light source modifier %s", mod); error(USER, errmsg); } if (*numContribs >= MAXMODLIST) { sprintf(errmsg, "too many modifiers (max. %d)", MAXMODLIST); error(INTERNAL, errmsg); } lutEntry -> key = mod; if (!binVal) /* Fall back to single bin if no value specified */ binVal = "0"; eBinVal = eparse(binVal); if (eBinVal -> type == NUM) { /* Bin value is constant */ binCnt = (int)(evalue(eBinVal) + 1.5); if (binCnt != 1) { sprintf(errmsg, "invalid non-zero constant for bin %s", binVal); error(USER, errmsg); } } else if (binCnt <= 0) { sprintf(errmsg, "undefined/invalid bin count for modifier %s", mod); error(USER, errmsg); } /* Allocate and init contributions */ contrib = (MODCONT *)malloc( sizeof(MODCONT) + sizeof(DCOLOR) * (binCnt - 1) ); if (!contrib) error(SYSTEM, "out of memory in addContribModifier()"); contrib -> outspec = NULL; contrib -> modname = mod; contrib -> params = binParm; contrib -> nbins = binCnt; contrib -> binv = eBinVal; contrib -> bin0 = 0; memset(contrib -> cbin, 0, sizeof(DCOLOR) * binCnt); lutEntry -> data = lutEntry -> data = (char *)contrib; ++(*numContribs); return(contrib); } -void addContribModfile(LUTAB *contribTab, unsigned *numContribs, +void addContribModfile ( + LUTAB *contribTab, unsigned *numContribs, char *modFile, char *binParm, char *binVal, int binCnt ) /* Add light source modifiers from file modFile to contribution lookup table * contribTab, and update numContribs. * NOTE: This code is adapted from rcontrib.c */ { - char *mods [MAXMODLIST]; - int i; + char *mods [MAXMODLIST]; + int i; /* Find file and store strings */ i = wordfile(mods, MAXMODLIST, getpath(modFile, getrlibpath(), R_OK)); if (i < 0) { sprintf(errmsg, "can't open modifier file %s", modFile); error(SYSTEM, errmsg); } if (*numContribs + i >= MAXMODLIST - 1) { sprintf(errmsg, "too many modifiers (max. %d) in file %s", MAXMODLIST - 1, modFile ); error(INTERNAL, errmsg); } for (i = 0; mods [i]; i++) /* Add each modifier */ addContribModifier(contribTab, numContribs, mods [i], binParm, binVal, binCnt ); } +static unsigned logDim (unsigned size) +/* Return log2(sqrt(size)) = l, where size = (2^l)(2^l). + Is there a faster (and portable!) way? */ +{ + unsigned i, dim, sz; + + if (!size) + return 0; + + for (i = 0, dim = 1, sz = size; sz >>= 1; dim++); + return dim >> 1; +} + + + +static int serialContribSourceBin (unsigned l, int x, int y) +/* Serialise 2D coordinates to 1D for array of size (2^l) x (2^l). + Returns -1 if coordinates invalid */ +{ + return x < 0 || y < 0 + ? -1 + : (x << l) + y; +} + + + +static void deserialContribSourceBin (unsigned l, int bin, int *x, int *y) +/* Deserialise 1D bin to 2D coordinates for array of size (2^l) x (2^l). + Returns -1 if bin invalid */ +{ + /* Obviously this is faster than integer division/modulo */ + *x = bin < 0 ? -1 : bin >> l; + *y = bin < 0 ? -1 : bin & (1 << l) - 1; +} + + + static int contribSourceBin (LUTAB *contribs, const RAY *ray) -/* Map contribution source ray to its bin for light source ray -> rsrc. +/* Map contribution source ray to its bin for light source ray -> rsrc, + using Shirley-Chiu disk-to-square mapping. Return invalid bin -1 if mapping failed. */ { const SRCREC *src; const OBJREC *srcMod; const MODCONT *srcCont; - double binReal; - int bin; + int i, bin, scRHS; + FVECT scNorm, scUp; + RAY srcRay; + double diskx, disky, scCoord [2]; + unsigned l, scDim, scBin [2]; /* Check we have a valid ray and contribution LUT */ if (!ray || !contribs) return -1; src = &source [ray -> rsrc]; srcMod = findmaterial(src -> so); srcCont = (MODCONT*)lu_find(contribs, srcMod -> oname) -> data; if (!srcCont) /* We're not interested in this source (modifier not in contrib LUT) */ return -1; - if (srcCont -> binv -> type != NUM) { - /* Binning function defined; set up shadow ray pointing to source - for evaluation */ - RAY srcRay; - int i; - - rayorigin(&srcRay, SHADOW, NULL, NULL); - srcRay.rsrc = ray -> rsrc; - VCOPY(srcRay.rorg, ray -> rop); - for (i = 0; i < 3; i++) - srcRay.rdir [i] = -ray -> rdir [i]; - - if (!(src -> sflags & SDISTANT - ? sourcehit(&srcRay) - : (*ofun[srcMod -> otype].funp)(srcMod, &srcRay) - )) - /* (Redundant?) sanity check for valid source ray? */ - return -1; - - worldfunc(RCCONTEXT, &srcRay); - set_eparams((char*)srcCont -> params); - } - - if ((binReal = evalue(srcCont -> binv)) < -.5) + /* Set up shadow ray pointing to source for disk2square mapping */ + rayorigin(&srcRay, SHADOW, NULL, NULL); + srcRay.rsrc = ray -> rsrc; + VCOPY(srcRay.rorg, ray -> rop); + for (i = 0; i < 3; i++) + srcRay.rdir [i] = -ray -> rdir [i]; + + if (!(src -> sflags & SDISTANT + ? sourcehit(&srcRay) + : (*ofun[srcMod -> otype].funp)(srcMod, &srcRay) + )) + /* (Redundant?) sanity check for valid source ray? */ return -1; - if ((bin = binReal + .5) >= srcCont -> nbins) { + worldfunc(RCCONTEXT, &srcRay); + set_eparams((char*)srcCont -> params); + + /* Get Shirley-Chiu mapping orientation from function parameters */ + scRHS = varvalue(PMAP_CONTRIB_SCRHS); + scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX); + scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY); + scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ); + scUp [0] = varvalue(PMAP_CONTRIB_SCUPX); + scUp [1] = varvalue(PMAP_CONTRIB_SCUPY); + scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ); + + /* Apply Shirley-Chiu mapping of disk to square */ + SDdisk2square(scCoord, diskx, disky); + l = logDim(srcCont -> nbins); + scDim = 1 << l; + /* Map Shirley-Chiu square coords to 2D bins and serialise */ + scBin [0] = scCoord [0] * scDim + 0.5; + scBin [1] = scCoord [1] * scDim + 0.5; + + if ((bin = serialContribSourceBin(l, scBin [0], scBin [1])) < 0) { error(WARNING, "Ignoring invalid bin in contribSourceBin()"); return -1; } return bin; } static PhotonContribSourceIdx newPhotonContribSource ( PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap ) /* Add contribution source ray for emitted contribution photon and save * light source index and binned direction. The current contribution source * is stored in pmap -> lastContribSrc. If the previous contribution source * spawned photons (i.e. has srcIdx >= 0), it's appended to contribSrcHeap. * If contribSrcRay == NULL, the current contribution source is still * flushed, but no new source is set. Returns updated contribution source * counter pmap -> numContribSrc. */ { if (!pmap || !contribSrcHeap) return 0; /* Check if last contribution source has spawned photons (srcIdx >= 0, * see newPhoton()), in which case we save it to the heap file before * clobbering it. (Note this is short-term storage, so we doan' need * da portable I/O stuff here). */ if (pmap -> lastContribSrc.srcIdx >= 0) { if (!fwrite(&pmap -> lastContribSrc, sizeof(PhotonContribSource), 1, contribSrcHeap )) error(SYSTEM, "failed writing photon contrib source in newPhotonContribSource()" ); pmap -> numContribSrc++; if (!pmap -> numContribSrc || pmap -> numContribSrc > PMAP_MAXCONTRIBSRC ); } /* Mark this contribution source unused with a negative source index until its path spawns a photon (see newPhoton()) */ pmap -> lastContribSrc.srcIdx = -1; /* Map ray to bin in anticipation that this contribution source will be used, since the ray will be lost once a photon is spawned */ pmap -> lastContribSrc.srcBin = contribSourceBin( pmap -> contribTab, contribSrcRay ); return pmap -> numContribSrc; } static PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ) /* Consolidate per-subprocess contribution sources heaps into array * pmap -> contribSrc. Returns offset for contribution source index * linearisation in pmap -> numContribSrc. The heap files in contribSrcHeap * are closed on return. */ { PhotonContribSourceIdx heapLen; unsigned heap; if (!pmap || !contribSrcHeap || !contribSrcOfs || !numHeaps) return 0; pmap -> numContribSrc = 0; for (heap = 0; heap < numHeaps; heap++) { contribSrcOfs [heap] = pmap -> numContribSrc; if (fseek(contribSrcHeap [heap], 0, SEEK_END) < 0) error(SYSTEM, "failed photon contribution source seek in buildContribSources()" ); pmap -> numContribSrc += heapLen = ftell(contribSrcHeap [heap]) / sizeof(PhotonContribSource); if (!(pmap -> contribSrc = realloc(pmap -> contribSrc, pmap -> numContribSrc * sizeof(PhotonContribSource) ))) error(SYSTEM, "failed photon contrib source alloc in buildContribSources()" ); rewind(contribSrcHeap [heap]); if (fread(pmap -> contribSrc + contribSrcOfs [heap], sizeof(PhotonContribSource), heapLen, contribSrcHeap [heap] ) != heapLen) error(SYSTEM, "failed reading photon contrib source in buildContribSources()" ); fclose(contribSrcHeap [heap]); unlink(contribSrcHeapFname [heap]); } return pmap -> numContribSrc; } static MODCONT *photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Locate photons near ray -> rop which originated from the light source modifier indexed by ray -> rsrc, and accumulate their contributions in pmap -> srcContrib. Return total contributions in irrad and pointer to binned contributions (or NULL if none were found). */ { const OBJREC *srcMod = findmaterial(source [ray -> rsrc].so); MODCONT *contrib = (MODCONT*)lu_find( pmap -> contribTab, srcMod -> oname ) -> data; Photon *photon; COLOR flux; PhotonSearchQueueNode *sqn; float r2, invArea; unsigned i; /* Zero bins for this source modifier */ memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins); setcolor(irrad, 0, 0, 0); if (!contrib || !pmap -> maxGather) /* Modifier not in LUT or zero bandwidth */ return NULL; /* Lookup photons */ pmap -> squeue.tail = 0; /* Pass light source index to filter in findPhotons() via pmap -> lastContribSrc */ pmap -> lastContribSrc.srcIdx = ray -> rsrc; 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 NULL; } /* Average radius^2 between furthest two photons to improve accuracy and get inverse search area 1 / (PI * r^2). */ /* XXX: OMIT extra normalisation factor 1/PI for ambient calculation? */ sqn = pmap -> squeue.node + 1; r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); invArea = 1 / (PI * PI * r2); /* Skip the extra photon */ for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); scalecolor(flux, invArea); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on photon distance */ scalecolor(flux, 2 * (1 - sqn -> dist2 / r2)); #endif addcolor(irrad, flux); addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } static void preComputeContrib (PhotonMap *pmap) /* Precompute contributions for a random subset of finalGather * pmap -> numPhotons photons, and build the precomputed contribution photon map, discarding the original photons. */ -/* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ +/* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { unsigned long i, numPreComp; unsigned j; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; PhotonMap nuPmap; repComplete = numPreComp = finalGather * pmap -> numPhotons; if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); -#if NIX +#if NIX fflush(stderr); -#endif +#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++) { do { /* 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 and normal at photon position. Set emitting light source index from origin. */ VCOPY(ray.rop, photon.pos); ray.rsrc = photonSrcIdx(pmap, &photon); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; /* Get contributions at photon position; retry with another photon * if no contribs found */ binnedContribs = photonContrib(pmap, &ray, ray.rcol); } while (!binnedContribs); /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray); /* TODO: Apply wavelet xform to binned contribs, apply lossy compression, encode to RGBE */ #if 0 encodeContrib(binnedContribs, pmap -> compRatio); #endif /* 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 contribution photon map\n"); #if NIX fflush(stderr); #endif } /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); /* TODO: Save wavelet coeffs here? */ } /* Defs for photon emission counter array passed by sub-processes to parent * via shared memory */ typedef unsigned long PhotonContribCnt; /* Indices for photon emission counter array: num photons stored and num * emitted per source */ #define PHOTONCNT_NUMPHOT 0 #define PHOTONCNT_NUMEMIT(n) (1 + n) /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ #define PMAP_CNTUPDATE 0xffL void distribPhotonContrib (PhotonMap* pmap, LUTAB *contribTab, unsigned numContribs, unsigned numProc ) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned srcIdx, proc; int shmFile, stat, pid; double *srcFlux, /* Emitted flux per light source */ srcDistribTarget; /* Target photon count per source */ PhotonContribCnt *photonCnt, /* Photon emission counter array */ lastPhotonCnt [PHOTONCNT_NUMEMIT(nsources)]; unsigned photonCntSize = (sizeof(PhotonContribCnt) * PHOTONCNT_NUMEMIT(nsources) ); FILE **contribSrcHeap = NULL; char **contribSrcHeapFname = NULL; PhotonContribSourceIdx *contribSrcOfs = NULL; pid_t procPids [PMAP_MAXPROC]; if (!pmap) error(USER, "no contribution photon map specified"); if (!nsources) error(USER, "no light sources"); if (nsources > MAXMODLIST) error(USER, "too many light sources"); if (!contribTab || !numContribs) error(USER, "no modifiers specified for contribution photon map"); /* Allocate photon flux per light source; this differs for every * source as all sources contribute the same number of distributed * photons (srcDistribTarget), hence the number of photons emitted per * source does not correlate with its emitted flux. The resulting flux * per photon is therefore adjusted individually for each source. */ if (!(srcFlux = calloc(nsources, sizeof(double)))) error(SYSTEM, "can't allocate source flux in distribPhotonContrib"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.src = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(USER, "can't allocate source partitions in distribPhotonContrib"); /* Initialise contrib photon map */ initPhotonMap(pmap, PMAP_TYPE_CONTRIB); initPmapContrib(contribTab); initPhotonHeap(pmap); initPhotonEmissionFuncs(); initPhotonScatterFuncs(); /* Per-subprocess / per-source target counts */ pmap -> distribTarget /= numProc; srcDistribTarget = nsources ? (double)pmap -> distribTarget / nsources : 0; if (!pmap -> distribTarget) error(INTERNAL, "no photons to distribute in distribPhotonContrib"); /* Get photon ports from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); + + /* Get polyhedral regions of interest */ + getPolyROIs(pmapROImodList); #if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0) error(SYSTEM, "failed shared mem init in distribPhotonContrib"); photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0 ); if (photonCnt == MAP_FAILED) error(SYSTEM, "failed shsared mem mapping in distribPhotonContrib"); #else /* Allocate photon counters statically on Windoze */ if (!(photonCnt = malloc(photonCntSize))) error(SYSTEM, "failed trivial malloc in distribPhotonContrib"); for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++) photonCnt [srcIdx] = 0; #endif /* NIX */ /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); 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 flux emitted from sources/ports * ============================================================= */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; const OBJREC *srcMod = findmaterial(source [srcIdx].so); srcFlux [srcIdx] = 0; /* Skip this source if its contributions are not sought */ if (!lu_find(pmap -> contribTab, srcMod -> oname) -> data) { sprintf(errmsg, "ignoring contributions from source %s", source [srcIdx].so -> oname ); error(WARNING, errmsg); continue; } emap.src = source + srcIdx; do { /* Need at least one iteration if no ports! */ 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++ ) { initPhotonEmission(&emap, pdfSamples); srcFlux [srcIdx] += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); if (srcFlux [srcIdx] < FTINY) { sprintf(errmsg, "source %s has zero emission", source [srcIdx].so -> oname ); error(WARNING, errmsg); } } /* Allocate & init per-subprocess contribution source heap files */ contribSrcHeap = calloc(numProc, sizeof(FILE*)); contribSrcHeapFname = calloc(numProc, sizeof(char*)); contribSrcOfs = calloc(numProc, sizeof(PhotonContribSourceIdx)); if (!contribSrcHeap || !contribSrcHeapFname || !contribSrcOfs) error(SYSTEM, "failed contribution source heap allocation " "in distribPhotonContrib()" ); for (proc = 0; proc < numProc; proc++) { contribSrcHeapFname [proc] = malloc(PMAP_TMPFNLEN); if (!contribSrcHeapFname [proc]) error(SYSTEM, "failed contribution source heap file allocation " "in distribPhotonContrib()" ); mktemp(strcpy(contribSrcHeapFname [proc], PMAP_TMPFNAME)); if (!(contribSrcHeap [proc] = fopen(contribSrcHeapFname [proc], "w+b"))) error(SYSTEM, "failed opening contribution source heap file " "in distribPhotonContrib()" ); } /* 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; opened and mmapped files inherited */ #else if (1) { /* No subprocess under Windoze */ #endif /* Local photon counters for this subprocess */ unsigned long lastNumPhotons = 0, localNumEmitted = 0; double photonFluxSum = 0; /* Accum. photon flux */ /* 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 PMAP_SIGUSR { double partNumEmit; unsigned long partEmitCnt; double srcPhotonFlux, avgPhotonFlux; unsigned portCnt, passCnt, prePassCnt; float srcPreDistrib; double srcNumEmit; /* # to emit from source */ unsigned long srcNumDistrib; /* # stored */ void sigUsrDiags() /* Loop diags via SIGUSR1 */ { sprintf(errmsg, "********************* Proc %d Diags *********************\n" "srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n" "srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n" "partNumEmit = %f\npartEmitCnt = %lu\n\n", proc, srcIdx, findmaterial(source [srcIdx].so) -> oname, portCnt, photonPorts [portCnt].so -> oname, passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux, partNumEmit, partEmitCnt ); eputs(errmsg); fflush(stderr); } } #endif #ifdef PMAP_SIGUSR signal(SIGUSR1, sigUsrDiags); #endif #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 /* ============================================================= * 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 * ============================================================= */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { const unsigned numEmitIdx = PHOTONCNT_NUMEMIT(srcIdx); #ifndef PMAP_SIGUSR unsigned portCnt, passCnt = 0, prePassCnt = 0; float srcPreDistrib = preDistrib; double srcNumEmit = 0; /* # to emit from source */ unsigned long srcNumDistrib = pmap -> numPhotons; /* #stored */ #else passCnt = prePassCnt = 0; srcPreDistrib = preDistrib; srcNumEmit = 0; /* # to emit from source */ srcNumDistrib = pmap -> numPhotons; /* # stored */ #endif if (srcFlux [srcIdx] < FTINY) /* Source has zero emission or was skipped in prepass because its contributions are not sought */ continue; while (passCnt < 2) { if (!passCnt) { /* INIT PASS 1 */ if (++prePassCnt > maxPreDistrib) { /* Warn if no photons contributed after sufficient * iterations; only output from subprocess 0 to reduce * console clutter */ if (!proc) { sprintf(errmsg, "source %s: too many prepasses, skipped", source [srcIdx].so -> oname ); error(WARNING, errmsg); } break; } /* Num to emit is fraction of target count */ srcNumEmit = srcPreDistrib * srcDistribTarget; } else { /* INIT PASS 2 */ #ifndef PMAP_SIGUSR double srcPhotonFlux, avgPhotonFlux; #endif /* Based on the outcome of the predistribution we can now * figure out how many more photons we have to emit from * the current source to meet the target count, * srcDistribTarget. This value is clamped to 0 in case * the target has already been exceeded in pass 1. * srcNumEmit and srcNumDistrib is the number of photons * emitted and distributed (stored) from the current * source in pass 1, respectively. */ srcNumDistrib = pmap -> numPhotons - srcNumDistrib; srcNumEmit *= srcNumDistrib ? max(srcDistribTarget/srcNumDistrib, 1) - 1 : 0; if (!srcNumEmit) /* No photons left to distribute in main pass */ break; srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit; avgPhotonFlux = photonFluxSum / (srcIdx + 1); if (avgPhotonFlux > FTINY && srcPhotonFlux / avgPhotonFlux < FTINY ) { /* Skip source if its photon flux is grossly below the * running average, indicating negligible contributions * at the expense of excessive distribution time; only * output from subproc 0 to reduce console clutter */ if (!proc) { sprintf(errmsg, "source %s: itsy bitsy photon flux, skipped", source [srcIdx].so -> oname ); error(WARNING, errmsg); } srcNumEmit = 0; /* Or just break??? */ } /* Update sum of photon flux per light source */ photonFluxSum += srcPhotonFlux; } portCnt = 0; do { /* Need at least one iteration if no ports! */ emap.src = source + srcIdx; 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++ ) { #ifndef PMAP_SIGUSR double partNumEmit; unsigned long partEmitCnt; #endif /* 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; * scale according to its normalised contribushunn to * the emitted source flux */ partNumEmit = (srcNumEmit * colorAvg(emap.partFlux) / srcFlux [srcIdx] ); partEmitCnt = (unsigned long)partNumEmit; /* Probabilistically account for fractional photons */ if (pmapRandom(cntState) < partNumEmit - partEmitCnt) partEmitCnt++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon according to PDF (if any). If * accepted, allocate associated contribution * origin, and trace through scene until * absorbed/leaked; emitPhoton() sets the emitting * light source index in photonRay */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ if (!emitPhoton(&emap, &photonRay)) continue; newPhotonContribSource(pmap, &photonRay, contribSrcHeap [proc] ); /* Set subprocess index in photonRay for post- * distrib contribution source index linearisation; * this is propagated with the contrib source index * in photonRay and set for photon hits by * newPhoton() */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter */ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update global counters shared with siblings in intervals to reduce overhead/contention */ lastPhotonCnt [numEmitIdx] = photonCnt [numEmitIdx]; photonCnt [numEmitIdx] += PMAP_CNTUPDATE; /* Check for overflow using previous values */ if (photonCnt [numEmitIdx] < lastPhotonCnt [numEmitIdx] ) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt [numEmitIdx], photonCnt [numEmitIdx] ); error(INTERNAL, errmsg); } /* Differentially increment photon counter */ lastNumPhotons = pmap -> numPhotons; photonCnt [PHOTONCNT_NUMPHOT] += pmap -> numPhotons - lastNumPhotons; /* Check for photon counter overflow (this could only happen before an emission counter overflow if the scene has an absurdly high albedo and/or very dense geometry) */ if (photonCnt [PHOTONCNT_NUMPHOT] < lastNumPhotons ) { sprintf(errmsg, "photon counter overlow " "(was: %ld, is: %ld)", lastNumPhotons, photonCnt [PHOTONCNT_NUMPHOT] ); error(INTERNAL, errmsg); } } } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime ) { unsigned s; repComplete = pmap -> distribTarget * numProc; repProgress = photonCnt [PHOTONCNT_NUMPHOT]; for (repEmitted = 0, s = 0; s < nsources; s++) repEmitted += photonCnt [numEmitIdx]; pmapDistribReport(); } #endif } portCnt++; } while (portCnt < numPhotonPorts); if (pmap -> numPhotons == srcNumDistrib) { /* Double predistrib factor in case no photons were stored * for this source and redo pass 1 */ srcPreDistrib *= 2; } else { /* Now do pass 2 */ passCnt++; } } } /* Flush heap buffa one final time to prevent data corruption */ flushPhotonHeap(pmap); /* Flush last contribution origin to origin heap file */ newPhotonContribSource(pmap, NULL, contribSrcHeap [proc]); /* Heap files closed automatically on exit fclose(pmap -> heap); fclose(orgHeap [proc]); */ #ifdef DEBUG_PMAP sprintf( errmsg, "Proc %d total %ld photons\n", proc, pmap -> numPhotons ); eputs(errmsg); fflush(stderr); #endif #ifdef PMAP_SIGUSR signal(SIGUSR1, SIG_DFL); #endif #if NIX /* Terminate subprocess */ exit(0); #endif } else { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotonContrib()" ); else /* Saves child process IDs */ procPids [proc] = pid; } } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses to complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); error(USER, "failed photon distribution"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ repComplete = pmap -> distribTarget * numProc; repProgress = photonCnt [PHOTONCNT_NUMPHOT]; for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++) repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; /* Get global photon count from shmem updated by subprocs */ pmap -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT]; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* ================================================================ * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc. * ================================================================ */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); if (!pmap -> numPhotons) error(USER, "empty contribution photon map"); /* Load per-subprocess contribution sources into pmap -> contribSrc */ /* Dumb compilers apparently need the char** cast */ pmap -> numContribSrc = buildContribSources(pmap, contribSrcHeap, (char**)contribSrcHeapFname, contribSrcOfs, numProc ); if (!pmap -> numContribSrc) error(INTERNAL, "no contribution sources in contribution photon map"); /* Set photon flux per source */ /* TODO: HOW/WHERE DO WE HANDLE COEFFICIENT MODE??? */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; #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("\nBuilding contribution photon map...\n"); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmap, srcFlux, contribSrcOfs, numProc); /* Precompute binned photon contributions */ if (verbose) eputs("\n"); preComputeContrib(contribPmap); /* Free per-subprocess origin heap files */ for (proc = 0; proc < numProc; proc++) free(contribSrcHeapFname [proc]); free(contribSrcHeapFname); free(contribSrcHeap); free(contribSrcOfs); if (verbose) eputs("\n"); } diff --git a/pmapcontrib.h b/pmapcontrib.h index 2423862..032c3fa 100644 --- a/pmapcontrib.h +++ b/pmapcontrib.h @@ -1,81 +1,111 @@ /* RCSid $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to mkpmap and rcontrib. 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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef _PMAPCONTRIB_H #define _PMAPCONTRIB_H #include "pmapdata.h" + #include "mrgbe.h" #ifndef MAXPROCESS #include "rcontrib.h" #endif #ifndef MAXMODLIST /* Max number of contributing sources */ #define MAXMODLIST 1024 #endif /* Filename templates for precomputed contrib photon map: Wavelet compressed contributions (per modifier) */ #define PMAP_CONTRIB_WAVELETSUFFIX "%s-%s.wvt" + + /* The following variables can be specified to override the orientation + of the Shirley-Chiu mapping (see also disk2square.cal). + We use the built-in functions in disk2square.c for efficiency and + in order to not depend on an external function file. These variables + merely mimic the function file interace. + + RHS : right-hand-coordinate system (-1 for left-hand) + rNx, rNy, rNz : surface normal + Ux, Uy, Uz : up vector (defines phi = 0) */ + #define PMAP_CONTRIB_SCRHS "RHS" + #define PMAP_CONTRIB_SCNORMX "rNx" + #define PMAP_CONTRIB_SCNORMY "rNy" + #define PMAP_CONTRIB_SCNORMZ "rNz" + #define PMAP_CONTRIB_SCUPX "Ux" + #define PMAP_CONTRIB_SCUPY "Uy" + #define PMAP_CONTRIB_SCUPZ "Uz" + #define PMAP_CONTRIB_SCDEFAULTS ( \ + "RHS=1; rNx=0; rNy=0; rNz=1; Ux=0; Uy=1; Uz=0;" \ + ) + + /* Maximum Shirley-Chiu binning resolution l*/ + #define PMAP_CONTRIB_MAXBINRES (mRGBE_DATABITS >> 1) + + /* Shirley-Chiu square dimensions and number of bins for resolution l */ + #define PMAP_CONTRIB_SCDIM(l) (1 << l) + #define PMAP_CONTRIB_SCBINS(l) (1 << (l << 1)) + + typedef struct { /* Contribution tagged by associated bin; expanded representation of below */ DCOLOR contrib; unsigned bin; } BinnedContrib; MODCONT *addContribModifier( LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, char *binVal, int binCnt ); /* Add light source modifier mod to contribution lookup table contribsTab, and update numContribs. Return initialised contribution data for this modifier. */ void addContribModfile( LUTAB *contribTab, unsigned *numContribs, char *modFile, char *binParm, char *binVal, int binCnt ); /* Add light source modifiers from file modFile to contribution lookup * table contribTab, and update numContribs. */ void initPmapContrib (LUTAB *contribTab); /* Set up photon map contributions (interface to rcmain.c) */ void distribPhotonContrib ( PhotonMap *pmap, LUTAB *contribTab, unsigned numContribs, unsigned numProc ); /* Emit photons with binned light source contributions, precompute their contributions and build photon map */ void preCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib); /* Fetch and decode precomputed light source contributions from single closest precomputed contribution photon at ray -> rop, and accumulate them in pmap -> contribTab. Also returns total precomputed contribs. */ #endif diff --git a/pmapdata.c b/pmapdata.c index 8c3e64f..0b25be6 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,836 +1,821 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif /* ========================================================================= 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, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" +#include "pmaproi.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL #ifdef PMAP_PHOTONFLOW , NULL, NULL #endif }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #include "pmaptkdt.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_PATHFILT /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init transient photon mapping stuff */ pmap -> velocity = photonVelocity; pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0; /* 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 photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; pmap -> numContribSrc = 0; pmap -> contribSrc = 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 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()"); /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); } #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; - FVECT photonDist; - - /* Store photon if within a region of interest (for ze Ecksperts!) - Note size of spherical ROI is squared. */ - for (i = 0; !inROI && i < pmapNumROI; i++) { - VSUB(photonDist, ray -> rop, pmapROI [i].pos); - - inROI = (PMAP_ROI_ISSPHERE(pmapROI + i) - ? DOT(photonDist, photonDist) <= pmapROI [i].siz [0] - : fabs(photonDist [0]) <= pmapROI [i].siz [0] && - fabs(photonDist [1]) <= pmapROI [i].siz [1] && - fabs(photonDist [2]) <= pmapROI [i].siz [2] - ); - } - if (!inROI) - return -1; - } + /* Store photon if within a region of interest (for ze Ecksperts!) */ + if (!photonInROI(ray)) + 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 photon's subprocess index */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { case PMAP_TYPE_CONTRIB: /* Set contrib photon's origin and subprocess index (the latter to * linearise the origin array after photon distribution completes). * Also set contribution source index, thereby marking it as used. * Note the contribution source bin has already been set by * newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; pmap -> lastContribSrc.srcIdx = ray -> rsrc; break; case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; break; default: /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? 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, PhotonContribSourceIdx *contribSrcOfs, 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, avgPathLen = 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 (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; /* 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()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } 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; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= 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); if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } 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 if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); 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; const RREAL *norm = ray -> ron, *pos = ray -> rop; 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; if (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } /* 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. */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #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 { /* 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_PATHFILT 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 26de791..36cd913 100644 --- a/pmapdata.h +++ b/pmapdata.h @@ -1,412 +1,412 @@ /* 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, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: 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 #include "ray.h" #include "pmaptype.h" #include "paths.h" #include "lookup.h" #include /* Source of a contribution photon. This consists of the emitting light source and binned direction. These records are only used to precompute contribution photons. They are referenced by contribution photons (see contribIdx field in struct Photon below) in a surjective mapping, since multiple photons may share the same emitting source and direction if they lie along its associated path. For this reason it is more efficient to factor this data out of the photons themselves and consolidate it here until the photons have been precomputed, after which it is no longer needed. */ typedef struct { int16 srcIdx, /* Index of emitting light source */ srcBin; /* Binned incident direction */ } PhotonContribSource; typedef uint32 PhotonPathID; typedef uint32 PhotonContribSourceIdx; #define PMAP_MAXCONTRIBSRC UINT32_MAX #define photonSrcIdx(pm, p) ((pm) -> contribSrc \ ? (pm) -> contribSrc [(p) -> aux.contribSrc].srcIdx \ : (p) -> aux.pathID\ ) #define photonSrcBin(pm, p) ( \ (pm) -> contribSrc [(p) -> aux.contribSrc].srcBin \ ) #define photonSrcMod(pm, p) findmaterial(source [photonSrcIdx(pm, p)].so) /* Multipurpose auxiliary photon attribute type */ typedef union { - /* Distance travelled by photon (= path length / time of flight) + /* Photon's propagation distance (= path length / time of flight) for temporal light flow */ float pathLen; /* Index into contribution photon's emitting source and binned direction; see struct PhotonContribSource above */ PhotonContribSourceIdx contribSrc; /* Unique path ID for all other photon types */ PhotonPathID pathID; } PhotonAuxAttrib; /* 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) typedef struct { float pos [3]; /* Photon position */ signed char norm [3]; /* Encoded normal / incident direction [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; }; /* Photon flux in watts or lumen / photon contribution [contrib photons] / average wavelet coefficient [precomputed contrib photons] */ #ifdef PMAP_FLOAT_FLUX COLOR flux; #else COLR flux; #endif /* Auxiliary field; this is a multipurpose, type-specific field used by the following photon types (as identified by enum PhotonMapType in pmaptype.h): PMAP_TYPE_CONTRIB: Index into photon map's contrib origin array. PMAP_TYPE_TEMPLIGHTFLOW: Distance travelled by photon / time of flight All others: Photon path ID. */ PhotonAuxAttrib aux; } 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" #include "pmaptkdt.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) */ /* ================================================================ * TRANSIENT PHOTON STUFF * ================================================================ */ double velocity, /* Speed of light in units of scene geometry [1/sec] */ time, /* Photons' time of flight for transient lookups */ minPathLen, maxPathLen, avgPathLen; /* Min/max/avg path length */ /* ================================================================ * CONTRIBUTION PHOTON STUFF * ================================================================ */ PhotonContribSource *contribSrc, /* Contribution source array */ lastContribSrc; /* Current contrib source */ PhotonContribSourceIdx numContribSrc; /* Number of contrib sources */ LUTAB *contribTab; /* LUT for binned contribs */ /* ================================================================ * 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) */ #ifdef PMAP_PATHFILT /* ================================================================ * PHOTON PATH FILTERING STUFF * ================================================================ */ 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 } 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 directPmap (photonMaps [PMAP_TYPE_DIRECT]) #define contribPmap (photonMaps [PMAP_TYPE_CONTRIB]) #define volumePmap (photonMaps [PMAP_TYPE_VOLUME]) #define transientPmap (photonMaps [PMAP_TYPE_TRANSIENT]) #ifdef PMAP_PHOTONFLOW /* Transient lightflow has precedence */ #define lightFlowPmap (transLightFlowPmap \ ? transLightFlowPmap : photonMaps [PMAP_TYPE_LIGHTFLOW] \ ) #define transLightFlowPmap (photonMaps [PMAP_TYPE_TRANSLIGHTFLOW]) #endif /* 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) #define isTransientPmap(p) ((p) -> type == PMAP_TYPE_TRANSIENT) #ifdef PMAP_PHOTONFLOW /* lightflow also implies transient lightflow */ #define isLightFlowPmap(p) ( \ (p) -> type == PMAP_TYPE_LIGHTFLOW || isTransLightFlowPmap(p) \ ) #define isTransLightFlowPmap(p) ( \ (p) -> type == PMAP_TYPE_TRANSLIGHTFLOW \ ) #endif 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, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ); /* Postprocess 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. */ /* OriginOfs is an array of index offsets for the contribution photon * origins in pmap->contribOrg generated by each of the nproc subprocesses * during contrib photon distribution (see distribPhotonContrib()). These * offsets are used to linearise the photon origin indices in the * postprocess. This linearisation is skipped if originOfs == NULL, * e.g. when building a global/caustic/volume photon map, where the * origins are serial path IDs. */ 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/pmapio.c b/pmapio.c index 4f12698..8f90bd8 100644 --- a/pmapio.c +++ b/pmapio.c @@ -1,221 +1,318 @@ #ifndef lint static const char RCSid[] = "$Id: pmapio.c,v 2.13 2021/02/18 17:08:50 rschregle Exp $"; #endif /* ========================================================================= Photon map portable file I/O Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS") (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: pmapio.c,v 2.13 2021/02/18 17:08:50 rschregle Exp $ */ #include "pmapio.h" #include "pmapdiag.h" #include "resolu.h" +#include + + + +extern char *octname; void savePhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ) { unsigned long i, j; FILE *file; if (!pmap || !pmap -> numPhotons || !validPmapType(pmap -> type)) { error(INTERNAL, "attempt to save empty or invalid photon map"); return; } if (verbose) { sprintf(errmsg, "Saving %s (%ld photons)\n", fname, pmap -> numPhotons ); eputs(errmsg); fflush(stderr); } if (!(file = fopen(fname, "wb"))) { sprintf(errmsg, "can't open photon map file %s", fname); error(SYSTEM, errmsg); } /* Write header */ newheader("RADIANCE", file); /* Write command line */ printargs(argc, argv, file); /* Include statistics in info text */ fprintf(file, "NumPhotons\t= %ld\n" "AvgFlux\t\t= [%.2e, %.2e, %.2e]\n" "Bbox\t\t= [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n" "CoG\t\t= [%.3f, %.3f, %.3f]\n" "MaxDist^2\t= %.3f\n" "Velocity\t= %g / sec\n" "Timespan\t= [%e, %e] sec\n" "AvgPathLen\t= %.3f\n", pmap -> numPhotons, pmap -> photonFlux [0], pmap -> photonFlux [1], pmap -> photonFlux [2], pmap -> minPos [0], pmap -> minPos [1], pmap -> minPos [2], pmap -> maxPos [0], pmap -> maxPos [1], pmap -> maxPos [2], pmap -> CoG [0], pmap -> CoG [1], pmap -> CoG [2], pmap -> CoGdist , pmap -> velocity, pmap -> minPathLen / pmap -> velocity, pmap -> maxPathLen / pmap -> velocity, pmap -> avgPathLen ); /* Write format, including human-readable file version */ fputformat((char *)pmapFormat [pmap -> type], file); fprintf(file, "VERSION=%s\n", PMAP_FILEVER); /* Empty line = end of header */ putc('\n', file); /* Write machine-readable file format version */ putstr(PMAP_FILEVER, file); /* Write number of photons as fixed size, which possibly results in * padding of MSB with 0 on some platforms. Unlike sizeof() however, * this ensures portability since this value may span 32 or 64 bits * depending on platform. */ putint(pmap -> numPhotons, PMAP_LONGSIZE, file); /* Write average photon flux */ for (j = 0; j < 3; j++) putflt(pmap -> photonFlux [j], file); /* Write max and min photon positions */ for (j = 0; j < 3; j++) { putflt(pmap -> minPos [j], file); putflt(pmap -> maxPos [j], file); } /* Write centre of gravity */ for (j = 0; j < 3; j++) putflt(pmap -> CoG [j], file); /* Write avg distance to centre of gravity */ putflt(pmap -> CoGdist, file); /* Save photon's velocity (speed of light in units of scene geometry), min/max/avg photon path length (= timespan and temporal CoG) */ putflt(pmap -> velocity, file); putflt(pmap -> minPathLen, file); putflt(pmap -> maxPathLen, file); putflt(pmap -> avgPathLen, file); /* Save photon storage */ #ifdef PMAP_OOC if (OOC_SavePhotons(pmap, file)) { #else if (kdT_SavePhotons(pmap, file)) { #endif sprintf(errmsg, "error writing photon map file %s", fname); error(SYSTEM, errmsg); } fclose(file); } PhotonMapType loadPhotonMap (PhotonMap *pmap, const char *fname) { PhotonMapType ptype = PMAP_TYPE_NONE; char format [MAXFMTLEN]; unsigned long i, j; FILE *file; if (!pmap) return PMAP_TYPE_NONE; if ((file = fopen(fname, "rb")) == NULL) { sprintf(errmsg, "can't open photon map file %s", fname); error(SYSTEM, errmsg); } /* Get format string */ strcpy(format, PMAP_FORMAT_GLOB); if (checkheader(file, format, NULL) != 1) { sprintf(errmsg, "photon map file %s has an unknown format", fname); error(USER, errmsg); } /* Identify photon map type from format string */ for (ptype = 0; ptype < NUM_PMAP_TYPES && strcmp(pmapFormat [ptype], format); ptype++ ); if (!validPmapType(ptype)) { sprintf(errmsg, "file %s contains an unknown photon map type", fname); error(USER, errmsg); } initPhotonMap(pmap, ptype); /* Get file format version and check for compatibility */ if (strcmp(getstr(format, file), PMAP_FILEVER)) error(USER, "incompatible photon map file format"); /* Get number of photons as fixed size, which possibly results in * padding of MSB with 0 on some platforms. Unlike sizeof() however, * this ensures portability since this value may span 32 or 64 bits * depending on platform. */ pmap -> numPhotons = getint(PMAP_LONGSIZE, file); /* Get average photon flux */ for (j = 0; j < 3; j++) pmap -> photonFlux [j] = getflt(file); /* Get max and min photon positions */ for (j = 0; j < 3; j++) { pmap -> minPos [j] = getflt(file); pmap -> maxPos [j] = getflt(file); } /* Get centre of gravity */ for (j = 0; j < 3; j++) pmap -> CoG [j] = getflt(file); /* Get avg distance to centre of gravity */ pmap -> CoGdist = getflt(file); /* Get photon's velocity (speed of light in units of scene geometry), min/max/avg photon path length (= timespan) and temporal CoG */ pmap -> velocity = getflt(file); pmap -> minPathLen = getflt(file); pmap -> maxPathLen = getflt(file); pmap -> avgPathLen = getflt(file); /* Load photon storage */ #ifdef PMAP_OOC if (OOC_LoadPhotons(pmap, file)) { #else if (kdT_LoadPhotons(pmap, file)) { #endif sprintf(errmsg, "error reading photon map file %s", fname); error(SYSTEM, errmsg); } fclose(file); return ptype; } + + +void savePmaps (const PhotonMap **pmaps, int argc, char **argv) +/* Wrapper to save all defined photon maps with specified command line */ +{ + unsigned t; + + for (t = 0; t < NUM_PMAP_TYPES; t++) { + if (pmaps [t]) + savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv); + } +} + + + +void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm) +/* Wrapper to save all defined photon maps with specified command line */ +{ + 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) { + /* Clamp lookup bandwidth to total number of photons (minus one, + since density estimate gets extra photon to obtain averaged + radius) */ + sprintf( + errmsg, "clamping density estimate bandwidth to %ld", + pm -> numPhotons + ); + error(WARNING, errmsg); + pm -> minGather = pm -> maxGather = pm -> numPhotons - 1; + } + } +} diff --git a/pmapio.h b/pmapio.h index ec12afb..1dbbf79 100644 --- a/pmapio.h +++ b/pmapio.h @@ -1,62 +1,68 @@ /* RCSid $Id: pmapio.h,v 2.8 2021/02/18 17:08:50 rschregle Exp $ */ /* ====================================================================== Photon map portable file I/O Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS") (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ====================================================================== $Id: pmapio.h,v 2.8 2021/02/18 17:08:50 rschregle Exp $ */ #ifndef PMAPIO_H #define PMAPIO_H #include "pmapdata.h" /* File format version with feature suffix */ #define PMAP_FILEVER_MAJ "5.0" #ifdef PMAP_OOC #define PMAP_FILEVER_TYP "o" #else #define PMAP_FILEVER_TYP "k" #endif #ifdef PMAP_FLOAT_FLUX #define PMAP_FILEVER_FLX "f" #else #define PMAP_FILEVER_FLX "" #endif #define PMAP_FILEVER (PMAP_FILEVER_MAJ PMAP_FILEVER_TYP PMAP_FILEVER_FLX) /* Maximum portable size of a long int in photon map file */ #define PMAP_LONGSIZE 8 /* Maximum portable size of a long int in photon map file */ #define PMAP_LONGSIZE 8 void savePhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ); /* Save portable photon map of specified type to filename, along with the corresponding command line. */ PhotonMapType loadPhotonMap (PhotonMap *pmap, const char* fname); /* Load portable photon map from specified filename, and return photon * map type as identified from file header. */ + void savePmaps (const PhotonMap **pmaps, int argc, char **argv); + /* Wrapper to save all defined photon maps with specified command line */ + + void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm); + /* Wrapper to save all defined photon maps with specified command line */ + #endif diff --git a/pmapmat.c b/pmapmat.c index 43762db..e72d5f8 100644 --- a/pmapmat.c +++ b/pmapmat.c @@ -1,2072 +1,2075 @@ #ifndef lint static const char RCSid[] = "$Id: pmapmat.c,v 2.24 2021/02/22 13:27:49 rschregle Exp $"; #endif /* ====================================================================== Photon map support routines for scattering by materials. 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") ====================================================================== */ #include "pmapmat.h" #include "pmapdata.h" #include "pmaprand.h" #include "otypes.h" #include "data.h" #include "func.h" #include "bsdf.h" #include /* Stuff ripped off from material modules */ #define MAXITER 10 #define SP_REFL 01 #define SP_TRAN 02 #define SP_PURE 04 #define SP_FLAT 010 #define SP_BADU 040 #define MLAMBDA 500 #define RINDEX 1.52 #define FRESNE(ci) (exp(-5.85*(ci)) - 0.00287989916) /* Get antimatter sensor orientation flags. HACK: these replace the last character in the object name! */ #define PMAP_GETSENSFLAGS(obj) ((obj) -> oname [strlen((obj) -> oname) - 1]) typedef struct { OBJREC *mp; RAY *rp; short specfl; COLOR mcolor, scolor; FVECT vrefl, prdir, pnorm; double alpha2, rdiff, rspec, trans, tdiff, tspec, pdot; } NORMDAT; typedef struct { OBJREC *mp; RAY *rp; short specfl; COLOR mcolor, scolor; FVECT vrefl, prdir, u, v, pnorm; double u_alpha, v_alpha, rdiff, rspec, trans, tdiff, tspec, pdot; } ANISODAT; typedef struct { OBJREC *mp; RAY *pr; DATARRAY *dp; COLOR mcolor; COLOR rdiff; COLOR tdiff; double rspec; double trans; double tspec; FVECT pnorm; double pdot; } BRDFDAT; typedef struct { OBJREC *mp; RAY *pr; FVECT pnorm; FVECT vray; double sr_vpsa [2]; RREAL toloc [3][3]; RREAL fromloc [3][3]; double thick; SDData *sd; COLOR runsamp; COLOR rdiff; COLOR tunsamp; COLOR tdiff; } BSDFDAT; extern const SDCDst SDemptyCD; /* Per-material scattering function dispatch table; return value is usually * zero, indicating photon termination */ int (*photonScatter [NUMOTYPE]) (OBJREC*, RAY*); /* List of antimatter sensor modifier names and associated object set */ char *photonSensorList [MAXSET + 1] = {NULL}; static OBJECT photonSensorSet [MAXSET + 1] = {0}; /* ================ General support routines ================ */ void photonRay (const RAY *rayIn, RAY *rayOut, int rayOutType, COLOR fluxAtten ) /* Spawn a new photon ray from a previous one; this is effectively a * customised rayorigin(). * A SPECULAR rayOutType flags this photon as _caustic_ for subsequent hits. * It is preserved for transferred rays (of type PMAP_XFER). * fluxAtten specifies the RGB attenuation of the photon flux effected by * the scattering material. The outgoing flux is then normalised to maintain * a uniform average of 1 over RGB. If fluxAtten == NULL, the flux remains * unchanged for the outgoing photon. fluxAtten is ignored for transferred * rays. * The ray direction is preserved for transferred rays, and undefined for * scattered rays and must be subsequently set by the caller. */ { rayorigin(rayOut, rayOutType, rayIn, NULL); if (rayIn) { /* Transfer flux */ copycolor(rayOut -> rcol, rayIn -> rcol); /* Copy caustic flag & direction for transferred rays */ if (rayOutType == PMAP_XFER) { /* rayOut -> rtype |= rayIn -> rtype & SPECULAR; */ rayOut -> rtype |= rayIn -> rtype; VCOPY(rayOut -> rdir, rayIn -> rdir); } else if (fluxAtten) { /* Attenuate and normalise flux for scattered rays */ multcolor(rayOut -> rcol, fluxAtten); colorNorm(rayOut -> rcol); } /* Propagate index of emitting light source */ rayOut -> rsrc = rayIn -> rsrc; /* Update maximum photon path distance; this is passed as aft plane distance to localhit() in tracePhoton(); it has no effect if negative, but limits the ray distance if positive. */ rayOut -> rmax = rayIn -> rmax - rayIn -> rot; } } static void addPhotons (const RAY *r) /* Insert photon hits, where applicable */ { if (!r -> rlvl) /* Add direct photon at primary hitpoint */ newPhoton(directPmap, r); else { /* Add global or precomputed photon at indirect hitpoint */ newPhoton(preCompPmap ? preCompPmap : globalPmap, r); /* Store caustic photon if specular flag set */ if (PMAP_CAUSTICRAY(r)) newPhoton(causticPmap, r); /* Store in contribution photon map */ newPhoton(contribPmap, r); } /* Add transient photon at all hitpoints */ newPhoton(transientPmap, r); } void getPhotonSensors (char **sensorList) /* Find antimatter geometry declared as photon sensors */ { OBJECT i; OBJREC *obj; char **lp; /* Init sensor set */ photonSensorSet [0] = 0; if (!sensorList [0]) return; for (i = 0; i < nobjects; i++) { obj = objptr(i); /* Insert object in sensor set if it's in the specified sensor list * and of type antimatter. Note the string comparison omits the last * character, since this contains the orientation flags. */ for (lp = sensorList; *lp; lp++) { if (!strncmp(obj -> oname, *lp, strlen(*lp) - 1)) { if (obj -> otype != MAT_CLIP) { sprintf(errmsg, "photon sensor modifier %s is not antimatter", obj -> oname ); error(USER, errmsg); } if (photonSensorSet [0] >= AMBLLEN) error(USER, "too many photon sensor modifiers"); /* Set orientation flags in object */ PMAP_GETSENSFLAGS(obj) = (*lp) [strlen(*lp) - 1]; insertelem(photonSensorSet, i); } } } if (!photonSensorSet [0]) error(USER, "no photon sensors found"); } + + + /* ================ Material specific scattering routines ================ */ static int isoSpecPhotonScatter (NORMDAT *nd, RAY *rayOut) /* Generate direction for isotropically specularly reflected or transmitted ray. Returns 1 if successful. */ { FVECT u, v, h; RAY *rayIn = nd -> rp; double d, d2, sinp, cosp; int niter, i = 0; /* Set up sample coordinates */ getperpendicular(u, nd -> pnorm, 1); fcross(v, nd -> pnorm, u); if (nd -> specfl & SP_REFL) { /* Specular reflection; make MAXITER attempts at getting a ray */ for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d); sinp = sin(d); d2 = pmapRandom(scatterState); d = d2 <= FTINY ? 1 : sqrt(nd -> alpha2 * -log(d2)); for (i = 0; i < 3; i++) h [i] = nd -> pnorm [i] + d * (cosp * u [i] + sinp * v [i]); d = -2 * DOT(h, rayIn -> rdir) / (1 + d * d); VSUM(rayOut -> rdir, rayIn -> rdir, h, d); if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY) return 1; } return 0; } else { /* Specular transmission; make MAXITER attempts at getting a ray */ for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d); sinp = sin(d); d2 = pmapRandom(scatterState); d = d2 <= FTINY ? 1 : sqrt(-log(d2) * nd -> alpha2); for (i = 0; i < 3; i++) rayOut -> rdir [i] = nd -> prdir [i] + d * (cosp * u [i] + sinp * v [i]); if (DOT(rayOut -> rdir, rayIn -> ron) < -FTINY) { normalize(rayOut -> rdir); return 1; } } return 0; } } static void diffPhotonScatter (FVECT normal, RAY* rayOut) /* Generate cosine-weighted direction for diffuse ray */ { const RREAL cosThetaSqr = pmapRandom(scatterState), cosTheta = sqrt(cosThetaSqr), sinTheta = sqrt(1 - cosThetaSqr), phi = 2 * PI * pmapRandom(scatterState), du = cos(phi) * sinTheta, dv = sin(phi) * sinTheta; FVECT u, v; int i = 0; /* Set up sample coordinates */ getperpendicular(u, normal, 1); fcross(v, normal, u); /* Convert theta & phi to cartesian */ for (i = 0; i < 3; i++) rayOut -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * normal [i]; normalize(rayOut -> rdir); } static int normalPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for isotropic material and recurse */ { NORMDAT nd; int i, hastexture; float xi, albedo, prdiff, ptdiff, prspec, ptspec; double d, fresnel; RAY rayOut; if (mat -> oargs.nfargs != (mat -> otype == MAT_TRANS ? 7 : 5)) objerror(mat, USER, "bad number of arguments"); /* Check for back side; reorient if back is visible */ if (rayIn -> rod < 0) if (!backvis && mat -> otype != MAT_TRANS) return 0; else { /* Get modifiers */ raytexture(rayIn, mat -> omod); flipsurface(rayIn); } else raytexture(rayIn, mat -> omod); nd.mp = mat; nd.rp = rayIn; /* Get material color */ copycolor(nd.mcolor, mat -> oargs.farg); /* Get roughness */ nd.specfl = 0; nd.alpha2 = mat -> oargs.farg [4]; if ((nd.alpha2 *= nd.alpha2) <= FTINY) nd.specfl |= SP_PURE; if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype)) nd.specfl |= SP_FLAT; /* Perturb normal */ if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)) )) nd.pdot = raynormal(nd.pnorm, rayIn); else { VCOPY(nd.pnorm, rayIn -> ron); nd.pdot = rayIn -> rod; } nd.pdot = max(nd.pdot, .001); /* Modify material color */ multcolor(nd.mcolor, rayIn -> pcol); nd.rspec = mat -> oargs.farg [3]; /* Approximate Fresnel term */ if (nd.specfl & SP_PURE && nd.rspec > FTINY) { fresnel = FRESNE(rayIn -> rod); nd.rspec += fresnel * (1 - nd.rspec); } else fresnel = 0; /* Transmission params */ if (mat -> otype == MAT_TRANS) { nd.trans = mat -> oargs.farg [5] * (1 - nd.rspec); nd.tspec = nd.trans * mat -> oargs.farg [6]; nd.tdiff = nd.trans - nd.tspec; } else nd.tdiff = nd.tspec = nd.trans = 0; /* Specular reflection params */ if (nd.rspec > FTINY) { /* Specular color */ if (mat -> otype != MAT_METAL) setcolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec); else if (fresnel > FTINY) { d = nd.rspec * (1 - fresnel); for (i = 0; i < 3; i++) nd.scolor [i] = fresnel + nd.mcolor [i] * d; } else { copycolor(nd.scolor, nd.mcolor); scalecolor(nd.scolor, nd.rspec); } } else setcolor(nd.scolor, 0, 0, 0); /* Diffuse reflection params */ nd.rdiff = 1 - nd.trans - nd.rspec; /* Set up probabilities */ prdiff = ptdiff = ptspec = colorAvg(nd.mcolor); prdiff *= nd.rdiff; ptdiff *= nd.tdiff; prspec = colorAvg(nd.scolor); ptspec *= nd.tspec; albedo = prdiff + ptdiff + prspec + ptspec; /* Insert direct and indirect photon hits if diffuse component */ if (prdiff > FTINY || ptdiff > FTINY) addPhotons(rayIn); xi = pmapRandom(rouletteState); if (xi > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prspec)) { /* Specular reflection */ nd.specfl |= SP_REFL; if (nd.specfl & SP_PURE) { /* Perfect specular reflection */ for (i = 0; i < 3; i++) { /* Reflected ray */ nd.vrefl [i] = rayIn -> rdir [i] + 2 * nd.pdot * nd.pnorm [i]; } /* Penetration? */ if (hastexture && DOT(nd.vrefl, rayIn -> ron) <= FTINY) for (i = 0; i < 3; i++) { /* Safety measure */ nd.vrefl [i] = rayIn -> rdir [i] + 2 * rayIn -> rod * rayIn -> ron [i]; } VCOPY(rayOut.rdir, nd.vrefl); } else if (!isoSpecPhotonScatter(&nd, &rayOut)) return 0; photonRay(rayIn, &rayOut, PMAP_SPECREFL, nd.scolor); } else if (xi > (albedo -= ptspec)) { /* Specular transmission */ nd.specfl |= SP_TRAN; if (hastexture) { /* Perturb */ for (i = 0; i < 3; i++) nd.prdir [i] = rayIn -> rdir [i] - rayIn -> pert [i]; if (DOT(nd.prdir, rayIn -> ron) < -FTINY) normalize(nd.prdir); else VCOPY(nd.prdir, rayIn -> rdir); } else VCOPY(nd.prdir, rayIn -> rdir); if ((nd.specfl & (SP_TRAN | SP_PURE)) == (SP_TRAN | SP_PURE)) /* Perfect specular transmission */ VCOPY(rayOut.rdir, nd.prdir); else if (!isoSpecPhotonScatter(&nd, &rayOut)) return 0; photonRay(rayIn, &rayOut, PMAP_SPECTRANS, nd.mcolor); } else if (xi > (albedo -= prdiff)) { /* Diffuse reflection */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor); diffPhotonScatter(hastexture ? nd.pnorm : rayIn -> ron, &rayOut); } else { /* Diffuse transmission */ flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor); if (hastexture) { FVECT bnorm; bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; diffPhotonScatter(bnorm, &rayOut); } else diffPhotonScatter(rayIn -> ron, &rayOut); } tracePhoton(&rayOut); return 0; } static void getacoords (ANISODAT *nd) /* Set up coordinate system for anisotropic sampling; cloned from aniso.c */ { MFUNC *mf; int i; mf = getfunc(nd -> mp, 3, 0x7, 1); setfunc(nd -> mp, nd -> rp); errno = 0; for (i = 0; i < 3; i++) nd -> u [i] = evalue(mf -> ep [i]); if (errno == EDOM || errno == ERANGE) nd -> u [0] = nd -> u [1] = nd -> u [2] = 0.0; if (mf -> fxp != &unitxf) multv3(nd -> u, nd -> u, mf -> fxp -> xfm); fcross(nd -> v, nd -> pnorm, nd -> u); if (normalize(nd -> v) == 0.0) { if (fabs(nd -> u_alpha - nd -> v_alpha) > 0.001) objerror(nd -> mp, WARNING, "illegal orientation vector"); getperpendicular(nd -> u, nd -> pnorm, 1); fcross(nd -> v, nd -> pnorm, nd -> u); nd -> u_alpha = nd -> v_alpha = sqrt(0.5 * (sqr(nd -> u_alpha) + sqr(nd -> v_alpha))); } else fcross(nd -> u, nd -> v, nd -> pnorm); } static int anisoSpecPhotonScatter (ANISODAT *nd, RAY *rayOut) /* Generate direction for anisotropically specularly reflected or transmitted ray. Returns 1 if successful. */ { FVECT h; double d, d2, sinp, cosp; int niter, i; RAY *rayIn = nd -> rp; if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype)) nd -> specfl |= SP_FLAT; /* set up coordinates */ getacoords(nd); if (rayOut -> rtype & TRANS) { /* Specular transmission */ if (DOT(rayIn -> pert, rayIn -> pert) <= sqr(FTINY)) VCOPY(nd -> prdir, rayIn -> rdir); else { /* perturb */ for (i = 0; i < 3; i++) nd -> prdir [i] = rayIn -> rdir [i] - rayIn -> pert [i]; if (DOT(nd -> prdir, rayIn -> ron) < -FTINY) normalize(nd -> prdir); else VCOPY(nd -> prdir, rayIn -> rdir); } /* Make MAXITER attempts at getting a ray */ for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d) * nd -> u_alpha; sinp = sin(d) * nd -> v_alpha; d = sqrt(sqr(cosp) + sqr(sinp)); cosp /= d; sinp /= d; d2 = pmapRandom(scatterState); d = d2 <= FTINY ? 1 : sqrt(-log(d2) / (sqr(cosp) / sqr(nd -> u_alpha) + sqr(sinp) / (nd -> v_alpha * nd -> u_alpha) ) ); for (i = 0; i < 3; i++) rayOut -> rdir [i] = nd -> prdir [i] + d * (cosp * nd -> u [i] + sinp * nd -> v [i]); if (DOT(rayOut -> rdir, rayIn -> ron) < -FTINY) { normalize(rayOut -> rdir); return 1; } } return 0; } else { /* Specular reflection */ /* Make MAXITER attempts at getting a ray */ for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d) * nd -> u_alpha; sinp = sin(d) * nd -> v_alpha; d = sqrt(sqr(cosp) + sqr(sinp)); cosp /= d; sinp /= d; d2 = pmapRandom(scatterState); d = d2 <= FTINY ? 1 : sqrt(-log(d2) / (sqr(cosp) / sqr(nd -> u_alpha) + sqr(sinp) / (nd->v_alpha * nd->v_alpha) ) ); for (i = 0; i < 3; i++) h [i] = nd -> pnorm [i] + d * (cosp * nd -> u [i] + sinp * nd -> v [i]); d = -2 * DOT(h, rayIn -> rdir) / (1 + d * d); VSUM(rayOut -> rdir, rayIn -> rdir, h, d); if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY) return 1; } return 0; } } static int anisoPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for anisotropic material and recurse */ { ANISODAT nd; float xi, albedo, prdiff, ptdiff, prspec, ptspec; RAY rayOut; if (mat -> oargs.nfargs != (mat -> otype == MAT_TRANS2 ? 8 : 6)) objerror(mat, USER, "bad number of real arguments"); nd.mp = mat; nd.rp = rayIn; /* get material color */ copycolor(nd.mcolor, mat -> oargs.farg); /* get roughness */ nd.specfl = 0; nd.u_alpha = mat -> oargs.farg [4]; nd.v_alpha = mat -> oargs.farg [5]; if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY) objerror(mat, USER, "roughness too small"); /* check for back side; reorient if back is visible */ if (rayIn -> rod < 0) if (!backvis && mat -> otype != MAT_TRANS2) return 0; else { /* get modifiers */ raytexture(rayIn, mat -> omod); flipsurface(rayIn); } else raytexture(rayIn, mat -> omod); /* perturb normal */ nd.pdot = max(raynormal(nd.pnorm, rayIn), .001); /* modify material color */ multcolor(nd.mcolor, rayIn -> pcol); nd.rspec = mat -> oargs.farg [3]; /* transmission params */ if (mat -> otype == MAT_TRANS2) { nd.trans = mat -> oargs.farg [6] * (1 - nd.rspec); nd.tspec = nd.trans * mat -> oargs.farg [7]; nd.tdiff = nd.trans - nd.tspec; if (nd.tspec > FTINY) nd.specfl |= SP_TRAN; } else nd.tdiff = nd.tspec = nd.trans = 0; /* specular reflection params */ if (nd.rspec > FTINY) { nd.specfl |= SP_REFL; /* compute specular color */ if (mat -> otype == MAT_METAL2) copycolor(nd.scolor, nd.mcolor); else setcolor(nd.scolor, 1, 1, 1); scalecolor(nd.scolor, nd.rspec); } else setcolor(nd.scolor, 0, 0, 0); /* diffuse reflection params */ nd.rdiff = 1 - nd.trans - nd.rspec; /* Set up probabilities */ prdiff = ptdiff = ptspec = colorAvg(nd.mcolor); prdiff *= nd.rdiff; ptdiff *= nd.tdiff; prspec = colorAvg(nd.scolor); ptspec *= nd.tspec; albedo = prdiff + ptdiff + prspec + ptspec; /* Insert direct and indirect photon hits if diffuse component */ if (prdiff > FTINY || ptdiff > FTINY) addPhotons(rayIn); xi = pmapRandom(rouletteState); if (xi > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prspec)) /* Specular reflection */ if (!(nd.specfl & SP_BADU)) { photonRay(rayIn, &rayOut, PMAP_SPECREFL, nd.scolor); if (!anisoSpecPhotonScatter(&nd, &rayOut)) return 0; } else return 0; else if (xi > (albedo -= ptspec)) /* Specular transmission */ if (!(nd.specfl & SP_BADU)) { /* Specular transmission */ photonRay(rayIn, &rayOut, PMAP_SPECTRANS, nd.mcolor); if (!anisoSpecPhotonScatter(&nd, &rayOut)) return 0; } else return 0; else if (xi > (albedo -= prdiff)) { /* Diffuse reflection */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor); diffPhotonScatter(nd.pnorm, &rayOut); } else { /* Diffuse transmission */ FVECT bnorm; flipsurface(rayIn); bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor); diffPhotonScatter(bnorm, &rayOut); } tracePhoton(&rayOut); return 0; } static double mylog (double x) /* special log for extinction coefficients; cloned from dielectric.c */ { if (x < 1e-40) return(-100.); if (x >= 1.) return(0.); return(log(x)); } static int dielectricPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for dielectric material and recurse */ { double cos1, cos2, nratio, d1, d2, refl; COLOR ctrans, talb; FVECT dnorm; int hastexture, i; RAY rayOut; if (mat -> oargs.nfargs != (mat -> otype == MAT_DIELECTRIC ? 5 : 8)) objerror(mat, USER, "bad arguments"); /* get modifiers */ raytexture(rayIn, mat -> omod); if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)))) /* Perturb normal */ cos1 = raynormal(dnorm, rayIn); else { VCOPY(dnorm, rayIn -> ron); cos1 = rayIn -> rod; } /* index of refraction */ nratio = mat -> otype == MAT_DIELECTRIC ? mat->oargs.farg[3] + mat->oargs.farg[4] / MLAMBDA : mat->oargs.farg[3] / mat->oargs.farg[7]; if (cos1 < 0) { /* inside */ hastexture = -hastexture; cos1 = -cos1; dnorm [0] = -dnorm [0]; dnorm [1] = -dnorm [1]; dnorm [2] = -dnorm [2]; setcolor(rayIn -> cext, -mylog(mat -> oargs.farg [0] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [1] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [2] * rayIn -> pcol [2]) ); setcolor(rayIn -> albedo, 0, 0, 0); rayIn -> gecc = 0; if (mat -> otype == MAT_INTERFACE) { setcolor(ctrans, -mylog(mat -> oargs.farg [4] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [5] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [6] * rayIn -> pcol [2]) ); setcolor(talb, 0, 0, 0); } else { copycolor(ctrans, cextinction); copycolor(talb, salbedo); } } else { /* outside */ nratio = 1.0 / nratio; setcolor(ctrans, -mylog(mat -> oargs.farg [0] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [1] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [2] * rayIn -> pcol [2]) ); setcolor(talb, 0, 0, 0); if (mat -> otype == MAT_INTERFACE) { setcolor(rayIn -> cext, -mylog(mat -> oargs.farg [4] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [5] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [6] * rayIn -> pcol [2]) ); setcolor(rayIn -> albedo, 0, 0, 0); rayIn -> gecc = 0; } } /* compute cos theta2 */ d2 = 1 - sqr(nratio) * (1 - sqr(cos1)); if (d2 < FTINY) { /* Total reflection */ refl = cos2 = 1.0; } else { /* Refraction, compute Fresnel's equations */ cos2 = sqrt(d2); d1 = cos1; d2 = nratio * cos2; d1 = (d1 - d2) / (d1 + d2); refl = sqr(d1); d1 = 1 / cos1; d2 = nratio / cos2; d1 = (d1 - d2) / (d1 + d2); refl += sqr(d1); refl *= 0.5; } if (pmapRandom(rouletteState) > refl) { /* Refraction */ photonRay(rayIn, &rayOut, PMAP_REFRACT, NULL); d1 = nratio * cos1 - cos2; for (i = 0; i < 3; i++) rayOut.rdir [i] = nratio * rayIn -> rdir [i] + d1 * dnorm [i]; if (hastexture && DOT(rayOut.rdir, rayIn->ron)*hastexture >= -FTINY) { d1 *= hastexture; for (i = 0; i < 3; i++) rayOut.rdir [i] = nratio * rayIn -> rdir [i] + d1 * rayIn -> ron [i]; normalize(rayOut.rdir); } copycolor(rayOut.cext, ctrans); copycolor(rayOut.albedo, talb); } else { /* Reflection */ photonRay(rayIn, &rayOut, PMAP_SPECREFL, NULL); VSUM(rayOut.rdir, rayIn -> rdir, dnorm, 2 * cos1); if (hastexture && DOT(rayOut.rdir, rayIn->ron) * hastexture <= FTINY) for (i = 0; i < 3; i++) rayOut.rdir [i] = rayIn -> rdir [i] + 2 * rayIn -> rod * rayIn -> ron [i]; } /* Ray is modified by medium defined by cext and albedo in * photonParticipate() */ tracePhoton(&rayOut); return 0; } static int glassPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for glass material and recurse */ { float albedo, xi, ptrans; COLOR mcolor, refl, trans; double pdot, cos2, d, r1e, r1m, rindex = 0.0; FVECT pnorm, pdir; int hastexture, i; RAY rayOut; /* check arguments */ if (mat -> oargs.nfargs == 3) rindex = RINDEX; else if (mat -> oargs.nfargs == 4) rindex = mat -> oargs.farg [3]; else objerror(mat, USER, "bad arguments"); copycolor(mcolor, mat -> oargs.farg); /* get modifiers */ raytexture(rayIn, mat -> omod); /* reorient if necessary */ if (rayIn -> rod < 0) flipsurface(rayIn); if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)))) pdot = raynormal(pnorm, rayIn); else { VCOPY(pnorm, rayIn -> ron); pdot = rayIn -> rod; } /* Modify material color */ multcolor(mcolor, rayIn -> pcol); /* angular transmission */ cos2 = sqrt((1 - 1 / sqr(rindex)) + sqr(pdot / rindex)); setcolor(mcolor, pow(mcolor [0], 1 / cos2), pow(mcolor [1], 1 / cos2), pow(mcolor [2], 1 / cos2) ); /* compute reflection */ r1e = (pdot - rindex * cos2) / (pdot + rindex * cos2); r1e *= r1e; r1m = (1 / pdot - rindex / cos2) / (1 / pdot + rindex / cos2); r1m *= r1m; for (i = 0; i < 3; i++) { double r1ed2, r1md2, d2; d = mcolor [i]; d2 = sqr(d); r1ed2 = sqr(r1e) * d2; r1md2 = sqr(r1m) * d2; /* compute transmittance */ trans [i] = 0.5 * d * (sqr(1 - r1e) / (1 - r1ed2) + sqr(1 - r1m) / (1 - r1md2) ); /* compute reflectance */ refl [i] = 0.5 * (r1e * (1 + (1 - 2 * r1e) * d2) / (1 - r1ed2) + r1m * (1 + (1 - 2 * r1m) * d2) / (1 - r1md2) ); } /* Set up probabilities */ ptrans = colorAvg(trans); albedo = colorAvg(refl) + ptrans; xi = pmapRandom(rouletteState); if (xi > albedo) /* Absorbed */ return 0; if (xi > (albedo -= ptrans)) { /* Transmitted */ if (hastexture) { /* perturb direction */ VSUM(pdir, rayIn -> rdir, rayIn -> pert, 2 * (1 - rindex)); if (normalize(pdir) == 0) { objerror(mat, WARNING, "bad perturbation"); VCOPY(pdir, rayIn -> rdir); } } else VCOPY(pdir, rayIn -> rdir); VCOPY(rayOut.rdir, pdir); photonRay(rayIn, &rayOut, PMAP_SPECTRANS, mcolor); } else { /* reflected ray */ VSUM(rayOut.rdir, rayIn -> rdir, pnorm, 2 * pdot); photonRay(rayIn, &rayOut, PMAP_SPECREFL, mcolor); } tracePhoton(&rayOut); return 0; } static int aliasPhotonScatter (OBJREC *mat, RAY *rayIn) /* Transfer photon scattering to alias target */ { OBJECT aliasObj; OBJREC aliasRec, *aliasPtr; /* Straight replacement? */ if (!mat -> oargs.nsargs) { /* Skip void modifier! */ if (mat -> omod != OVOID) { mat = objptr(mat -> omod); photonScatter [mat -> otype] (mat, rayIn); } return 0; } /* Else replace alias */ if (mat -> oargs.nsargs != 1) objerror(mat, INTERNAL, "bad # string arguments"); aliasPtr = mat; aliasObj = objndx(aliasPtr); /* Follow alias trail */ do { aliasObj = aliasPtr -> oargs.nsargs == 1 ? lastmod(aliasObj, aliasPtr -> oargs.sarg [0]) : aliasPtr -> omod; if (aliasObj < 0) objerror(aliasPtr, USER, "bad reference"); aliasPtr = objptr(aliasObj); } while (aliasPtr -> otype == MOD_ALIAS); /* Copy alias object */ aliasRec = *aliasPtr; /* Substitute modifier */ aliasRec.omod = mat -> omod; /* Replacement scattering routine */ photonScatter [aliasRec.otype] (&aliasRec, rayIn); /* Avoid potential memory leak? */ if (aliasRec.os != aliasPtr -> os) { if (aliasPtr -> os) free_os(aliasPtr); aliasPtr -> os = aliasRec.os; } return 0; } static int clipPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for antimatter material and recurse */ { OBJECT obj = objndx(mat), mod, cset [MAXSET + 1], *modset; int entering, inside = 0, i; const RAY *rp; RAY rayOut; if ((modset = (OBJECT*)mat -> os) == NULL) { if (mat -> oargs.nsargs < 1 || mat -> oargs.nsargs > MAXSET) objerror(mat, USER, "bad # arguments"); modset = (OBJECT*)malloc((mat -> oargs.nsargs + 1) * sizeof(OBJECT)); if (modset == NULL) error(SYSTEM, "out of memory in clipPhotonScatter"); modset [0] = 0; for (i = 0; i < mat -> oargs.nsargs; i++) { if (!strcmp(mat -> oargs.sarg [i], VOIDID)) continue; if ((mod = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "unknown modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, WARNING, errmsg); continue; } if (inset(modset, mod)) { objerror(mat, WARNING, "duplicate modifier"); continue; } insertelem(modset, mod); } mat -> os = (char*)modset; } if (rayIn -> clipset != NULL) setcopy(cset, rayIn -> clipset); else cset [0] = 0; entering = rayIn -> rod > 0; if (inset(photonSensorSet, obj)) { /* Material defined as virtual sensor; store photon if incident * from side(s) corresponding to orientation flags. Note PMAPSENSBI * subsumes both PMAP_SENSFWD and PMAP_SENSBWD. */ const char sensFlags = PMAP_GETSENSFLAGS(mat); if (entering && sensFlags & PMAP_SENSFWD || !entering && sensFlags & PMAP_SENSBWD ) if (entering) addPhotons(rayIn); else { /* Flip normal if incident from back, then flip back after saving photon */ flipsurface(rayIn); addPhotons(rayIn); flipsurface(rayIn); } } for (i = modset [0]; i > 0; i--) { if (entering) { if (!inset(cset, modset [i])) { if (cset [0] >= MAXSET) error(INTERNAL, "set overflow in clipPhotonScatter"); insertelem(cset, modset [i]); } } else if (inset(cset, modset [i])) deletelem(cset, modset [i]); } rayIn -> newcset = cset; if (strcmp(mat -> oargs.sarg [0], VOIDID)) { for (rp = rayIn; rp -> parent != NULL; rp = rp -> parent) { if ( !(rp -> rtype & RAYREFL) && rp->parent->ro != NULL && inset(modset, rp -> parent -> ro -> omod)) { if (rp -> parent -> rod > 0) inside++; else inside--; } } if (inside > 0) { flipsurface(rayIn); mat = objptr(lastmod(obj, mat -> oargs.sarg [0])); photonScatter [mat -> otype] (mat, rayIn); return 0; } } /* Else transfer ray */ photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); return 0; } static int mirrorPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for mirror material and recurse */ { RAY rayOut; int rpure = 1, i; FVECT pnorm; double pdot; float albedo; COLOR mcolor; /* check arguments */ if (mat -> oargs.nfargs != 3 || mat -> oargs.nsargs > 1) objerror(mat, USER, "bad number of arguments"); /* back is black */ if (rayIn -> rod < 0) return 0; /* get modifiers */ raytexture(rayIn, mat -> omod); /* assign material color */ copycolor(mcolor, mat -> oargs.farg); multcolor(mcolor, rayIn -> pcol); /* Set up probabilities */ albedo = colorAvg(mcolor); if (pmapRandom(rouletteState) > albedo) /* Absorbed */ return 0; /* compute reflected ray */ photonRay(rayIn, &rayOut, PMAP_SPECREFL, mcolor); if (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)) { /* use textures */ pdot = raynormal(pnorm, rayIn); for (i = 0; i < 3; i++) rayOut.rdir [i] = rayIn -> rdir [i] + 2 * pdot * pnorm [i]; rpure = 0; } /* Check for penetration */ if (rpure || DOT(rayOut.rdir, rayIn -> ron) <= FTINY) for (i = 0; i < 3; i++) rayOut.rdir [i] = rayIn -> rdir [i] + 2 * rayIn -> rod * rayIn -> ron [i]; tracePhoton(&rayOut); return 0; } static int mistPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray within mist and recurse */ { COLOR mext; RREAL re, ge, be; RAY rayOut; /* check arguments */ if (mat -> oargs.nfargs > 7) objerror(mat, USER, "bad arguments"); if (mat -> oargs.nfargs > 2) { /* compute extinction */ copycolor(mext, mat -> oargs.farg); /* get modifiers */ raytexture(rayIn, mat -> omod); multcolor(mext, rayIn -> pcol); } else setcolor(mext, 0, 0, 0); photonRay(rayIn, &rayOut, PMAP_XFER, NULL); if (rayIn -> rod > 0) { /* entering ray */ addcolor(rayOut.cext, mext); if (mat -> oargs.nfargs > 5) copycolor(rayOut.albedo, mat -> oargs.farg + 3); if (mat -> oargs.nfargs > 6) rayOut.gecc = mat -> oargs.farg [6]; } else { /* leaving ray */ re = max(rayIn -> cext [0] - mext [0], cextinction [0]); ge = max(rayIn -> cext [1] - mext [1], cextinction [1]); be = max(rayIn -> cext [2] - mext [2], cextinction [2]); setcolor(rayOut.cext, re, ge, be); if (mat -> oargs.nfargs > 5) copycolor(rayOut.albedo, salbedo); if (mat -> oargs.nfargs > 6) rayOut.gecc = seccg; } tracePhoton(&rayOut); return 0; } static int mx_dataPhotonScatter (OBJREC *mat, RAY *rayIn) /* Pass photon on to materials selected by mixture data */ { OBJECT obj; double coef, pt [MAXDIM]; DATARRAY *dp; OBJECT mod [2]; MFUNC *mf; int i; if (mat -> oargs.nsargs < 6) objerror(mat, USER, "bad # arguments"); obj = objndx(mat); for (i = 0; i < 2; i++) if (!strcmp(mat -> oargs.sarg [i], VOIDID)) mod [i] = OVOID; else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, USER, errmsg); } dp = getdata(mat -> oargs.sarg [3]); i = (1 << dp -> nd) - 1; mf = getfunc(mat, 4, i << 5, 0); setfunc(mat, rayIn); errno = 0; for (i = 0; i < dp -> nd; i++) { pt [i] = evalue(mf -> ep [i]); if (errno) { objerror(mat, WARNING, "compute error"); return 0; } } coef = datavalue(dp, pt); errno = 0; coef = funvalue(mat -> oargs.sarg [2], 1, &coef); if (errno) objerror(mat, WARNING, "compute error"); else { OBJECT mxMod = mod [pmapRandom(rouletteState) < coef ? 0 : 1]; if (mxMod != OVOID) { mat = objptr(mxMod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ RAY rayOut; photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } } return 0; } static int mx_pdataPhotonScatter (OBJREC *mat, RAY *rayIn) /* Pass photon on to materials selected by mixture picture */ { OBJECT obj; double col [3], coef, pt [MAXDIM]; DATARRAY *dp; OBJECT mod [2]; MFUNC *mf; int i; if (mat -> oargs.nsargs < 7) objerror(mat, USER, "bad # arguments"); obj = objndx(mat); for (i = 0; i < 2; i++) if (!strcmp(mat -> oargs.sarg [i], VOIDID)) mod [i] = OVOID; else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, USER, errmsg); } dp = getpict(mat -> oargs.sarg [3]); mf = getfunc(mat, 4, 0x3 << 5, 0); setfunc(mat, rayIn); errno = 0; pt [1] = evalue(mf -> ep [0]); pt [0] = evalue(mf -> ep [1]); if (errno) { objerror(mat, WARNING, "compute error"); return 0; } for (i = 0; i < 3; i++) col [i] = datavalue(dp + i, pt); errno = 0; coef = funvalue(mat -> oargs.sarg [2], 3, col); if (errno) objerror(mat, WARNING, "compute error"); else { OBJECT mxMod = mod [pmapRandom(rouletteState) < coef ? 0 : 1]; if (mxMod != OVOID) { mat = objptr(mxMod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ RAY rayOut; photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } } return 0; } static int mx_funcPhotonScatter (OBJREC *mat, RAY *rayIn) /* Pass photon on to materials selected by mixture function */ { OBJECT obj, mod [2]; int i; double coef; MFUNC *mf; if (mat -> oargs.nsargs < 4) objerror(mat, USER, "bad # arguments"); obj = objndx(mat); for (i = 0; i < 2; i++) if (!strcmp(mat -> oargs.sarg [i], VOIDID)) mod [i] = OVOID; else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, USER, errmsg); } mf = getfunc(mat, 3, 0x4, 0); setfunc(mat, rayIn); errno = 0; /* bound coefficient */ coef = min(1, max(0, evalue(mf -> ep [0]))); if (errno) objerror(mat, WARNING, "compute error"); else { OBJECT mxMod = mod [pmapRandom(rouletteState) < coef ? 0 : 1]; if (mxMod != OVOID) { mat = objptr(mxMod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ RAY rayOut; photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } } return 0; } static int pattexPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for pattern or texture modifier and recurse. This code is brought to you by Henkel! :^) */ { RAY rayOut; /* Get pattern */ ofun [mat -> otype].funp(mat, rayIn); if (mat -> omod != OVOID) { /* Scatter using modifier (if any) */ mat = objptr(mat -> omod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } return 0; } static int setbrdfunc(BRDFDAT *bd) /* Set up brdf function and variables; ripped off from m_brdf.c */ { FVECT v; if (setfunc(bd -> mp, bd -> pr) == 0) return 0; /* (Re)Assign func variables */ multv3(v, bd -> pnorm, funcxf.xfm); varset("NxP", '=', v [0] / funcxf.sca); varset("NyP", '=', v [1] / funcxf.sca); varset("NzP", '=', v [2] / funcxf.sca); varset("RdotP", '=', bd -> pdot <= -1. ? -1. : bd -> pdot >= 1. ? 1. : bd -> pdot); varset("CrP", '=', colval(bd -> mcolor, RED)); varset("CgP", '=', colval(bd -> mcolor, GRN)); varset("CbP", '=', colval(bd -> mcolor, BLU)); return 1; } static int brdfPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for BRTDfunc material and recurse. Only ideal reflection and transmission are sampled for the specular componentent. */ { int hitfront = 1, hastexture, i; BRDFDAT nd; RAY rayOut; COLOR rspecCol, tspecCol; double prDiff, ptDiff, prSpec, ptSpec, albedo, xi; MFUNC *mf; FVECT bnorm; /* Check argz */ if (mat -> oargs.nsargs < 10 || mat -> oargs.nfargs < 9) objerror(mat, USER, "bad # arguments"); nd.mp = mat; nd.pr = rayIn; /* Dummiez */ nd.rspec = nd.tspec = 1.0; nd.trans = 0.5; /* Diffuz reflektanz */ if (rayIn -> rod > 0.0) setcolor(nd.rdiff, mat -> oargs.farg[0], mat -> oargs.farg [1], mat -> oargs.farg [2]); else setcolor(nd.rdiff, mat-> oargs.farg [3], mat -> oargs.farg [4], mat -> oargs.farg [5]); /* Diffuz tranzmittanz */ setcolor(nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7], mat -> oargs.farg [8]); /* Get modz */ raytexture(rayIn, mat -> omod); hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)); if (hastexture) { /* Perturb normal */ nd.pdot = raynormal(nd.pnorm, rayIn); } else { VCOPY(nd.pnorm, rayIn -> ron); nd.pdot = rayIn -> rod; } if (rayIn -> rod < 0.0) { /* Orient perturbed valuz */ nd.pdot = -nd.pdot; for (i = 0; i < 3; i++) { nd.pnorm [i] = -nd.pnorm [i]; rayIn -> pert [i] = -rayIn -> pert [i]; } hitfront = 0; } /* Get pattern kolour, modify diffuz valuz */ copycolor(nd.mcolor, rayIn -> pcol); multcolor(nd.rdiff, nd.mcolor); multcolor(nd.tdiff, nd.mcolor); /* Load cal file, evaluate spekula refl/tranz varz */ nd.dp = NULL; mf = getfunc(mat, 9, 0x3f, 0); setbrdfunc(&nd); errno = 0; setcolor(rspecCol, evalue(mf->ep[0]), evalue(mf->ep[1]), evalue(mf->ep[2])); setcolor(tspecCol, evalue(mf->ep[3]), evalue(mf->ep[4]), evalue(mf->ep[5])); if (errno == EDOM || errno == ERANGE) objerror(mat, WARNING, "compute error"); else { /* Set up probz */ prDiff = colorAvg(nd.rdiff); ptDiff = colorAvg(nd.tdiff); prSpec = colorAvg(rspecCol); ptSpec = colorAvg(tspecCol); albedo = prDiff + ptDiff + prSpec + ptSpec; } /* Insert direct and indirect photon hitz if diffuz komponent */ if (prDiff > FTINY || ptDiff > FTINY) addPhotons(rayIn); /* Stochastically sample absorption or scattering evenz */ if ((xi = pmapRandom(rouletteState)) > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prSpec)) { /* Ideal spekula reflekzion */ photonRay(rayIn, &rayOut, PMAP_SPECREFL, rspecCol); VSUM(rayOut.rdir, rayIn -> rdir, nd.pnorm, 2 * nd.pdot); checknorm(rayOut.rdir); } else if (xi > (albedo -= ptSpec)) { /* Ideal spekula tranzmission */ photonRay(rayIn, &rayOut, PMAP_SPECTRANS, tspecCol); if (hastexture) { /* Perturb direkzion */ VSUB(rayOut.rdir, rayIn -> rdir, rayIn -> pert); if (normalize(rayOut.rdir) == 0.0) { objerror(mat, WARNING, "illegal perturbation"); VCOPY(rayOut.rdir, rayIn -> rdir); } } else VCOPY(rayOut.rdir, rayIn -> rdir); } else if (xi > (albedo -= prDiff)) { /* Diffuz reflekzion */ if (!hitfront) flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor); diffPhotonScatter(nd.pnorm, &rayOut); } else { /* Diffuz tranzmission */ if (hitfront) flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor); bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; diffPhotonScatter(bnorm, &rayOut); } tracePhoton(&rayOut); return 0; } int brdf2PhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for procedural or data driven BRDF material and recurse. Only diffuse reflection and transmission are sampled. */ { BRDFDAT nd; RAY rayOut; double dtmp, prDiff, ptDiff, albedo, xi; MFUNC *mf; FVECT bnorm; /* Check argz */ if (mat -> oargs.nsargs < (hasdata(mat -> otype) ? 4 : 2) || mat -> oargs.nfargs < (mat -> otype == MAT_TFUNC || mat -> otype == MAT_TDATA ? 6 : 4)) objerror(mat, USER, "bad # arguments"); if (rayIn -> rod < 0.0) { /* Hit backside; reorient if visible, else transfer photon */ if (!backvis) { photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); return 0; } raytexture(rayIn, mat -> omod); flipsurface(rayIn); } else raytexture(rayIn, mat -> omod); nd.mp = mat; nd.pr = rayIn; /* Material kolour */ setcolor(nd.mcolor, mat -> oargs.farg [0], mat -> oargs.farg [1], mat -> oargs.farg [2]); /* Spekula komponent */ nd.rspec = mat -> oargs.farg [3]; /* Tranzmittanz */ if (mat -> otype == MAT_TFUNC || mat -> otype == MAT_TDATA) { nd.trans = mat -> oargs.farg [4] * (1.0 - nd.rspec); nd.tspec = nd.trans * mat -> oargs.farg [5]; dtmp = nd.trans - nd.tspec; setcolor(nd.tdiff, dtmp, dtmp, dtmp); } else { nd.tspec = nd.trans = 0.0; setcolor(nd.tdiff, 0.0, 0.0, 0.0); } /* Reflektanz */ dtmp = 1.0 - nd.trans - nd.rspec; setcolor(nd.rdiff, dtmp, dtmp, dtmp); /* Perturb normal */ nd.pdot = raynormal(nd.pnorm, rayIn); /* Modify material kolour */ multcolor(nd.mcolor, rayIn -> pcol); multcolor(nd.rdiff, nd.mcolor); multcolor(nd.tdiff, nd.mcolor); /* Load auxiliary filez */ if (hasdata(mat -> otype)) { nd.dp = getdata(mat -> oargs.sarg [1]); getfunc(mat, 2, 0, 0); } else { nd.dp = NULL; getfunc(mat, 1, 0, 0); } /* Set up probz */ prDiff = colorAvg(nd.rdiff); ptDiff = colorAvg(nd.tdiff); albedo = prDiff + ptDiff; /* Insert direct and indirect photon hitz if diffuz komponent */ if (prDiff > FTINY || ptDiff > FTINY) addPhotons(rayIn); /* Stochastically sample absorption or scattering evenz */ if ((xi = pmapRandom(rouletteState)) > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prDiff)) { /* Diffuz reflekzion */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff); diffPhotonScatter(nd.pnorm, &rayOut); } else { /* Diffuz tranzmission */ flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; diffPhotonScatter(bnorm, &rayOut); } tracePhoton(&rayOut); return 0; } /* ====================================================================== The following code is (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") ====================================================================== */ static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for BSDF modifier and recurse. */ { int hasthick = (mat->otype == MAT_BSDF); int hitFront; SDError err; SDValue bsdfVal; FVECT upvec; MFUNC *mf; BSDFDAT nd; RAY rayOut; COLOR bsdfRGB; int transmitted; double prDiff, ptDiff, prDiffSD, ptDiffSD, prSpecSD, ptSpecSD, albedo, xi; const double patAlb = bright(rayIn -> pcol); /* Following code adapted from m_bsdf() */ /* Check arguments */ if ( mat -> oargs.nsargs < hasthick+5 || mat -> oargs.nfargs > 9 || mat -> oargs.nfargs % 3 ) objerror(mat, USER, "bad # arguments"); hitFront = (rayIn -> rod > 0); /* Load cal file */ mf = hasthick ? getfunc(mat, 5, 0x1d, 1) : getfunc(mat, 4, 0xe, 1); /* Get thickness */ nd.thick = 0; if (hasthick) { nd.thick = evalue(mf -> ep [0]); if ((-FTINY <= nd.thick) & (nd.thick <= FTINY)) nd.thick = .0; } /* Get BSDF data */ nd.sd = loadBSDF(mat -> oargs.sarg [hasthick]); /* Extra diffuse reflectance from material def */ if (hitFront) { if (mat -> oargs.nfargs < 3) setcolor(nd.rdiff, .0, .0, .0); else setcolor( nd.rdiff, mat -> oargs.farg [0], mat -> oargs.farg [1], mat -> oargs.farg [2] ); } else if (mat -> oargs.nfargs < 6) { /* Check for absorbing backside */ if (!backvis && !nd.sd -> rb && !nd.sd -> tf) { SDfreeCache(nd.sd); return 0; } setcolor(nd.rdiff, .0, .0, .0); } else setcolor( nd.rdiff, mat -> oargs.farg [3], mat -> oargs.farg [4], mat -> oargs.farg [5] ); /* Extra diffuse transmittance from material def */ if (mat -> oargs.nfargs < 9) setcolor(nd.tdiff, .0, .0, .0); else setcolor( nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7], mat -> oargs.farg [8] ); nd.mp = mat; nd.pr = rayIn; /* Get modifiers */ raytexture(rayIn, mat -> omod); /* Modify diffuse values */ multcolor(nd.rdiff, rayIn -> pcol); multcolor(nd.tdiff, rayIn -> pcol); /* Get up vector & xform to world coords */ upvec [0] = evalue(mf -> ep [hasthick+0]); upvec [1] = evalue(mf -> ep [hasthick+1]); upvec [2] = evalue(mf -> ep [hasthick+2]); if (mf -> fxp != &unitxf) { multv3(upvec, upvec, mf -> fxp -> xfm); nd.thick *= mf -> fxp -> sca; } if (rayIn -> rox) { multv3(upvec, upvec, rayIn -> rox -> f.xfm); nd.thick *= rayIn -> rox -> f.sca; } /* Perturb normal */ raynormal(nd.pnorm, rayIn); /* Xform incident dir to local BSDF coords */ err = SDcompXform(nd.toloc, nd.pnorm, upvec); if (!err) { nd.vray [0] = -rayIn -> rdir [0]; nd.vray [1] = -rayIn -> rdir [1]; nd.vray [2] = -rayIn -> rdir [2]; err = SDmapDir(nd.vray, nd.toloc, nd.vray); } if (!err) err = SDinvXform(nd.fromloc, nd.toloc); if (err) { objerror(mat, WARNING, "Illegal orientation vector"); return 0; } /* Determine BSDF resolution */ err = SDsizeBSDF( nd.sr_vpsa, nd.vray, NULL, SDqueryMin + SDqueryMax, nd.sd ); if (err) objerror(mat, USER, transSDError(err)); nd.sr_vpsa [0] = sqrt(nd.sr_vpsa [0]); nd.sr_vpsa [1] = sqrt(nd.sr_vpsa [1]); /* Orient perturbed normal towards incident side */ if (!hitFront) { nd.pnorm [0] = -nd.pnorm [0]; nd.pnorm [1] = -nd.pnorm [1]; nd.pnorm [2] = -nd.pnorm [2]; } /* Get scatter probabilities (weighted by pattern except for spec refl) * prDiff, ptDiff: extra diffuse component in material def * prDiffSD, ptDiffSD: diffuse (constant) component in SDF * prSpecSD, ptSpecSD: non-diffuse ("specular") component in SDF * albedo: sum of above, inverse absorption probability */ prDiff = colorAvg(nd.rdiff); ptDiff = colorAvg(nd.tdiff); prDiffSD = patAlb * SDdirectHemi(nd.vray, SDsampDf | SDsampR, nd.sd); ptDiffSD = patAlb * SDdirectHemi(nd.vray, SDsampDf | SDsampT, nd.sd); prSpecSD = SDdirectHemi(nd.vray, SDsampSp | SDsampR, nd.sd); ptSpecSD = patAlb * SDdirectHemi(nd.vray, SDsampSp | SDsampT, nd.sd); albedo = prDiff + ptDiff + prDiffSD + ptDiffSD + prSpecSD + ptSpecSD; /* if (albedo > 1) objerror(mat, WARNING, "Invalid albedo"); */ /* Insert direct and indirect photon hits if diffuse component */ if (prDiff + ptDiff + prDiffSD + ptDiffSD > FTINY) addPhotons(rayIn); xi = pmapRandom(rouletteState); if (xi > albedo) /* Absorbtion */ return 0; transmitted = 0; if ((xi -= prDiff) <= 0) { /* Diffuse reflection (extra component in material def) */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff); diffPhotonScatter(nd.pnorm, &rayOut); } else if ((xi -= ptDiff) <= 0) { /* Diffuse transmission (extra component in material def) */ photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); diffPhotonScatter(nd.pnorm, &rayOut); transmitted = 1; } else { /* Sample SDF */ if ((xi -= prDiffSD) <= 0) { /* Diffuse SDF reflection (constant component) */ if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampDf | SDsampR, nd.sd ))) objerror(mat, USER, transSDError(err)); /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); multcolor(bsdfRGB, rayIn -> pcol); photonRay(rayIn, &rayOut, PMAP_DIFFREFL, bsdfRGB); } else if ((xi -= ptDiffSD) <= 0) { /* Diffuse SDF transmission (constant component) */ if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampDf | SDsampT, nd.sd ))) objerror(mat, USER, transSDError(err)); /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); multcolor(bsdfRGB, rayIn -> pcol); addcolor(bsdfRGB, nd.tdiff); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, bsdfRGB); transmitted = 1; } else if ((xi -= prSpecSD) <= 0) { /* Non-diffuse ("specular") SDF reflection */ if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampSp | SDsampR, nd.sd ))) objerror(mat, USER, transSDError(err)); ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); photonRay(rayIn, &rayOut, PMAP_SPECREFL, bsdfRGB); } else { /* Non-diffuse ("specular") SDF transmission */ if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampSp | SDsampT, nd.sd ))) objerror(mat, USER, transSDError(err)); /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); multcolor(bsdfRGB, rayIn -> pcol); photonRay(rayIn, &rayOut, PMAP_SPECTRANS, bsdfRGB); transmitted = 1; } /* Xform outgoing dir to world coords */ if ((err = SDmapDir(rayOut.rdir, nd.fromloc, nd.vray))) { objerror(mat, USER, transSDError(err)); return 0; } } /* Clean up */ SDfreeCache(nd.sd); /* Offset outgoing photon origin by thickness to bypass proxy geometry */ if (transmitted && nd.thick != 0) VSUM(rayOut.rorg, rayOut.rorg, rayIn -> ron, -nd.thick); tracePhoton(&rayOut); return 0; } static int lightPhotonScatter (OBJREC* mat, RAY* ray) /* Light sources doan' reflect, mang */ { return 0; } void initPhotonScatterFuncs () /* Init photonScatter[] dispatch table */ { int i; /* Catch-all for inconsistencies */ for (i = 0; i < NUMOTYPE; i++) photonScatter [i] = o_default; photonScatter [MAT_LIGHT] = photonScatter [MAT_ILLUM] = photonScatter [MAT_GLOW] = photonScatter [MAT_SPOT] = lightPhotonScatter; photonScatter [MAT_PLASTIC] = photonScatter [MAT_METAL] = photonScatter [MAT_TRANS] = normalPhotonScatter; photonScatter [MAT_PLASTIC2] = photonScatter [MAT_METAL2] = photonScatter [MAT_TRANS2] = anisoPhotonScatter; photonScatter [MAT_DIELECTRIC] = photonScatter [MAT_INTERFACE] = dielectricPhotonScatter; photonScatter [MAT_MIST] = mistPhotonScatter; photonScatter [MAT_GLASS] = glassPhotonScatter; photonScatter [MAT_CLIP] = clipPhotonScatter; photonScatter [MAT_MIRROR] = mirrorPhotonScatter; photonScatter [MIX_FUNC] = mx_funcPhotonScatter; photonScatter [MIX_DATA] = mx_dataPhotonScatter; photonScatter [MIX_PICT]= mx_pdataPhotonScatter; photonScatter [PAT_BDATA] = photonScatter [PAT_CDATA] = photonScatter [PAT_BFUNC] = photonScatter [PAT_CFUNC] = photonScatter [PAT_CPICT] = photonScatter [TEX_FUNC] = photonScatter [TEX_DATA] = pattexPhotonScatter; photonScatter [MOD_ALIAS] = aliasPhotonScatter; photonScatter [MAT_BRTDF] = brdfPhotonScatter; photonScatter [MAT_PFUNC] = photonScatter [MAT_MFUNC] = photonScatter [MAT_PDATA] = photonScatter [MAT_MDATA] = photonScatter [MAT_TFUNC] = photonScatter [MAT_TDATA] = brdf2PhotonScatter; photonScatter [MAT_BSDF] = photonScatter [MAT_ABSDF] = bsdfPhotonScatter; } diff --git a/pmapparm.c b/pmapparm.c index 9dd7789..20c8203 100644 --- a/pmapparm.c +++ b/pmapparm.c @@ -1,131 +1,128 @@ #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, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapparm.c,v 2.9 2018/02/02 19:47:55 rschregle Exp $ */ #include "pmapparm.h" #include "pmapdata.h" +#include "otspecial.h" +#include "otypes.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) */ compressRatio = 0.9, /* Compression ratio for precomputed contribution photon map */ photonMaxPathDist = 0; /* Maximum cumulative distance of photon path */ double photonVelocity = 299792458.0; /* Photon velocity in units of scene geometry (default m/s) */ #ifdef PMAP_OOC float pmapCachePageSize = 8; /* OOC cache pagesize as multiple * of maxGather */ unsigned long pmapCacheSize = 1e6; /* OOC cache size in photons */ #endif #ifdef PMAP_PHOTONFLOW int sphericalIrrad = 0, /* Toggle for spherical / planar irrad from lightflow pmap */ hemiLightFlow = 0; /* Toggle hemispherical/spherical lookups for lightflow pmap */ #endif -/* Regions of interest */ -unsigned pmapNumROI = 0; -PhotonMapROI *pmapROI = NULL; - 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 */ /* Per photon map lookup params */ PhotonMapParams pmapParams [NUM_PMAP_TYPES] = { {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.} #ifdef PMAP_PHOTONFLOW , {NULL, 0, 0, 0}, {NULL, 0, 0, 0} #endif }; 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) -> contribTab = NULL; (*pm) -> time = parm -> time; 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 5bcd48a..690ca57 100644 --- a/pmapparm.h +++ b/pmapparm.h @@ -1,95 +1,85 @@ /* RCSid $Id: pmapparm.h,v 2.10 2021/04/14 11:26:25 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, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapparm.h,v 2.10 2021/04/14 11:26:25 rschregle Exp $ */ #ifndef PMAPPARAMS_H #define PMAPPARAMS_H - #include "pmaptype.h" + #include "pmaptype.h" /* Struct for passing lookup params per photon map for rpict/rtrace/rvu */ typedef struct { char *fileName; /* Photon map file */ unsigned minGather, maxGather; /* Num photons to gather */ unsigned long distribTarget; /* Num photons to store */ double time; /* Transient photon time */ } PhotonMapParams; - /* Region of interest */ - typedef struct { - /* siz [1], siz [2] <= 0 --> sphere, else rectangle */ - float pos [3], siz [3]; - } PhotonMapROI; - - #define PMAP_ROI_ISSPHERE(roi) ((roi)->siz[1] <= 0 && (roi)->siz[2] <= 0) - #define PMAP_ROI_SETSPHERE(roi) ((roi)->siz[1] = (roi)->siz[2] = -1) + #include "ray.h" /* Requires PhotonMapParams */ 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]) #define transientPmapParams (pmapParams [PMAP_TYPE_TRANSIENT]) #ifdef PMAP_PHOTONFLOW #define lightFlowParams (pmapParams [PMAP_TYPE_LIGHTFLOW]) #define transLightFlowParams (pmapParams [PMAP_TYPE_TRANSLIGHTFLOW]) #endif extern float pdfSamples, preDistrib, finalGather, gatherTolerance, maxDistFix, pmapMaxDist, photonMaxPathDist, compressRatio; extern double photonVelocity; 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 int sphericalIrrad, hemiLightFlow; #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/pmapray.c b/pmapray.c index f9a4b23..ef035cc 100644 --- a/pmapray.c +++ b/pmapray.c @@ -1,79 +1,80 @@ #ifndef lint static const char RCSid[] = "$Id: pmapray.c,v 2.7 2016/11/02 22:09:14 greg Exp $"; #endif /* ================================================================== Photon map interface to RADIANCE raycalls 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: pmapray.c,v 2.7 2016/11/02 22:09:14 greg Exp $ */ #include "pmapray.h" -#include "pmap.h" +#include "pmapio.h" +#include "pmutil.h" void ray_init_pmap () /* Interface to ray_init(); init & load photon maps */ { loadPmaps(photonMaps, pmapParams); } void ray_done_pmap () /* Interface to ray_done(); free photon maps */ { cleanUpPmaps(photonMaps); } void ray_save_pmap (RAYPARAMS *rp) /* Interface to ray_save(); save photon map params */ { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmapParams [t].fileName) rp -> pmapParams [t].fileName = savqstr(pmapParams [t].fileName); else rp -> pmapParams [t].fileName = NULL; rp -> pmapParams [t].minGather = pmapParams [t].minGather; rp -> pmapParams [t].maxGather = pmapParams [t].maxGather; rp -> pmapParams [t].distribTarget = pmapParams [t].distribTarget; } } void ray_restore_pmap (RAYPARAMS *rp) /* Interface to ray_restore(); restore photon mapping params */ { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { pmapParams [t].fileName = rp -> pmapParams [t].fileName; pmapParams [t].minGather = rp -> pmapParams [t].minGather; pmapParams [t].maxGather = rp -> pmapParams [t].maxGather; pmapParams [t].distribTarget = rp -> pmapParams [t].distribTarget; } } void ray_defaults_pmap (RAYPARAMS *rp) /* Interface to ray_defaults(); set photon mapping defaults */ { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { rp -> pmapParams [t].fileName = NULL; rp -> pmapParams [t].minGather = 0; rp -> pmapParams [t].maxGather = 0; rp -> pmapParams [t].distribTarget = 0; } } diff --git a/pmaproi.c b/pmaproi.c new file mode 100644 index 0000000..f736f9c --- /dev/null +++ b/pmaproi.c @@ -0,0 +1,136 @@ +#ifndef lint +static const char RCSid[] = "$Id$"; +#endif + +/* + ====================================================================== + Regions of interest used by mkpmap to constrain location of photons + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com} + (c) Fraunhofer Institute for Solar Energy Systems, + supported by the German Research Foundation + (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components", + SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") + (c) Tokyo University of Science, + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") + ====================================================================== + + $Id$ +*/ + + +#include "pmaproi.h" +#include "otspecial.h" +#include "otypes.h" + + +/* Rectangular / spherical regions of interest */ +unsigned pmapNumROI = 0; +PhotonMapROI *pmapROI = NULL; + +/* Polyhedral regions of interest by modifier */ +char *pmapROImodList [MAXSET + 1] = {NULL}; +OBJECT polyROIset [MAXSET + 1] = {0}; + + + +void getPolyROIs (char **roiModList) +/* Find geometry comprising polyhedral regions of interest from modifiers in + * roiModList */ +{ + OBJECT i; + OBJREC *obj, *mat; + char **lp; + unsigned lastNumPolys = 0; + + /* Init poly ROI set */ + polyROIset [0] = 0; + + if (!roiModList [0]) + return; + + for (lp = roiModList; *lp; lp++) { + lastNumPolys = polyROIset [0]; + + for (i = 0; i < nobjects; i++) { + obj = objptr(i); + mat = findmaterial(obj); + + /* Check if object is geometry and its modifier is in list */ + if (issurface(obj -> otype) && mat && !strcmp(mat -> oname, *lp)) { + if (polyROIset [0] >= PMAP_MAXROI) + error(USER, "too many ROI objects"); + + insertelem(polyROIset, i); + } + } + + if (polyROIset [0] == lastNumPolys) { + /* No geometry found for this modifier */ + sprintf(errmsg, "no ROI geometry found for modifier %s", *lp); + error(USER, errmsg); + } + } + + if (!polyROIset [0]) + /* Made redundant by check above? */ + error(USER, "no ROI modifiers found"); +} + + + +int photonInROI (const RAY *photonRay) +/* Returns nonzero if photon defined by photonRay lies in any defined region + * of interest */ +{ + unsigned inROI = 0, i; + + if (pmapNumROI && pmapROI) { + /* Check rectangular / spherical ROIs */ + FVECT photonDist; + + for (i = 0; i < pmapNumROI; i++) { + VSUB(photonDist, photonRay -> rop, pmapROI [i].pos); + /* NOTE: Size of spherical ROI is _squared_ */ + inROI = (PMAP_ROI_ISSPHERE(pmapROI + i) + ? DOT(photonDist, photonDist) <= pmapROI [i].siz [0] + : fabs(photonDist [0]) <= pmapROI [i].siz [0] && + fabs(photonDist [1]) <= pmapROI [i].siz [1] && + fabs(photonDist [2]) <= pmapROI [i].siz [2] + ); + if (inROI) + /* Bail out early, skip remaining ROIs */ + return 1; + } + } + + if (polyROIset [0]) { + /* Check polyhedral ROIs; send ray and count number of crossed + * polygons in ROI set */ + RAY testRay; + + /* Init inclusion test ray from photon origin and incident direction */ + rayorigin(&testRay, PRIMARY, NULL, NULL); + /* Offset origin along reversed photon direction in case we have a + * surface-bound (i.e. global, caustic) photon, thereby avoiding + * potential boundary issues. */ + VSUM(testRay.rorg, photonRay -> rop, photonRay -> rdir, -FTINY); + VCOPY(testRay.rdir, photonRay -> rdir); + testRay.rmax = FHUGE; + + /* Follow inclusion test ray until it leaves the scene */ + while (localhit(&testRay, &thescene)) { + /* Toggle odd/even number of crossings if intersected object is in + ROI set (note this works since inset() returns 0 or 1!) */ + inROI ^= inset(polyROIset, testRay.robj); + /* Update ray origin with intersection for next iteration */ + VCOPY(testRay.rorg, testRay.rop); + } + } + + return inROI; +} diff --git a/pmaproi.h b/pmaproi.h new file mode 100644 index 0000000..f6e5189 --- /dev/null +++ b/pmaproi.h @@ -0,0 +1,53 @@ +/* + ====================================================================== + Regions of interest used by mkpmap to constrain location of photons + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com} + (c) Fraunhofer Institute for Solar Energy Systems, + supported by the German Research Foundation + (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components", + SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") + (c) Tokyo University of Science, + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") + ====================================================================== + + $Id$ +*/ + +#ifndef _PMAP_ROI_H + #define _PMAP_ROI_H + + #include "ray.h" + + #define PMAP_MAXROI 1024 + #define PMAP_ROI_ISSPHERE(roi) ((roi)->siz[1] <= 0 && (roi)->siz[2] <= 0) + #define PMAP_ROI_SETSPHERE(roi) ((roi)->siz[1] = (roi)->siz[2] = -1) + + + /* Region of interest defined by geometry */ + typedef struct { + /* siz [1], siz [2] <= 0 --> sphere, else rectangle */ + float pos [3], siz [3]; + } PhotonMapROI; + + + /* Rectangular / spherical regions of interest */ + extern unsigned pmapNumROI; + extern PhotonMapROI *pmapROI; + + /* Polyhedral regions of interest by modifier */ + extern char *pmapROImodList []; + extern OBJECT polyROIset []; + + + void getPolyROIs (char **roiModList); + /* Find geometry of polyhedral regions of interest from modifier list */ + + int photonInROI (const RAY *photonRay); + /* Returns nonzero if photon defined by photonRay lies in any defined + * region of interest */ +#endif diff --git a/pmcontrib2.c b/pmcontrib2.c index 9920e60..6782c03 100644 --- a/pmcontrib2.c +++ b/pmcontrib2.c @@ -1,190 +1,190 @@ #ifndef lint static const char RCSid[] = "$Id: pmcontrib2.c,v 2.5 2018/11/08 00:54:07 greg Exp $"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to rcontrib. 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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id: pmcontrib2.c,v 2.5 2018/11/08 00:54:07 greg Exp $ */ #include "pmapcontrib.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmapio.h" #include "pmapdiag.h" #include "otypes.h" #include "otspecial.h" static void setPmapContribParams (PhotonMap *pmap, LUTAB *contribTab) /* Set parameters for light source contributions */ { /* Set light source modifier list and appropriate callback to extract * their contributions from the photon map */ if (pmap) { pmap -> contribTab = contribTab; pmap -> lookup = preCompPhotonContrib; /* XXX: Need tighter tolerance for contribution photon lookups? */ pmap -> gatherTolerance = 1.0; } } static int checkSourceContrib (const LUENT *lutEntry, void *p) /* Check lookup table entry for light source contributions; Returns 0 on success. */ { char *modName = (char *)((MODCONT *)lutEntry -> data) -> modname; OBJECT modObj = modifier(modName); if (modObj == OVOID) { sprintf(errmsg, "invalid modifier %s", modName); error(USER, errmsg); } else if (!islight(objptr(modObj) -> otype)) { sprintf(errmsg, "%s is not a light source modifier", modName); error(USER, errmsg); } else return 0; } - + void initPmapContrib (LUTAB *srcContrib) { unsigned t; - + /* XXX: Redundant check? */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) { sprintf( errmsg, "%s photon map does not support contributions", pmapName [t] ); error(USER, errmsg); } /* Check for valid modifiers in lookup table; this doesn't return on error. */ lu_doall(srcContrib, checkSourceContrib, NULL); /* Get params */ setPmapContribParams(contribPmap, srcContrib); #if 0 if (contribPhotonMapping) { if (contribPmap -> maxGather < numSrcContrib) { /* Warn if density estimate bandwidth is lower than modifier * count, since some contributions will be missing */ error(WARNING, "photon density estimate bandwidth too low," " contributions may be underestimated"); /* Sanity check */ checkPmapContribs(contribPmap, srcContrib); } #endif } void preCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib) /* Fetch and decode precomputed light source contributions from single closest precomputed contribution photon at ray -> rop, and accumulate them in pmap -> srcContrib. Also returns total precomputed contributions. */ { unsigned b; RREAL rayCoeff [3]; Photon photon; /* Get cumulative path coefficient up to photon lookup point */ raycontrib(rayCoeff, ray, PRIMARY); setcolor(totalContrib, 0, 0, 0); /* Ignore sources */ if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) return; if (find1Photon(pmap, ray, &photon)) { /* Get average wavelet coeff as total contributions */ getPhotonFlux(&photon, totalContrib); if (pmap -> contribTab) { /* TODO: fetch decoded binned contribs here */ const OBJREC *mod = photonSrcMod(pmap, &photon); MODCONT *contrib = (MODCONT*)lu_find( pmap -> contribTab, mod -> oname ) -> data; COLOR binContrib; if (!contrib) { /* Contribs from this light source not sought; this indicates * the photon map was precomputed with different rcontrib * parameters */ sprintf(errmsg, "missing LUT entry for light source modifier %s; " "mismatched contribution photon map parameters?", mod -> oname ); error(CONSISTENCY, errmsg); } #if 0 /* TODO */ decodeBins(&photon, &srcBins); if (srcBins -> numBins != contrib -> nbins) { /* Number of precomputed and allocated bins don't match */ sprintf( errmsg, "wrong number of bins for light source modifier %s; " "mismatched contribution photon map parameters?", srcMod -> oname ); error(CONSISTENCY, errmsg); } for (b = 0; b <= srcBins -> numBins; b++) { /* Accumulate contribs from decoded bin */ copycolor(binContrib, srcBins [b].contrib); if (!contrib) { /* Ray coefficient mode; normalise by light source radiance * after applying pattern */ int i; raytexture(ray, srcMod -> omod); setcolor( ray -> rcol, srcMod -> oargs.farg [0], srcMod -> oargs.farg [1], srcMod -> oargs.farg [2] ); multcolor(ray -> rcol, ray -> pcol); for (i = 0; i < 3; i++) binContrib [i] = ray -> rcol [i] ? binContrib [i] / ray -> rcol [i] : 0; } multcolor(contrib, rayCoeff); addcolor(srcContrib -> cbin [b], binContrib); } #endif } } } diff --git a/pmutil.c b/pmutil.c index b92b0d1..3696f79 100644 --- a/pmutil.c +++ b/pmutil.c @@ -1,132 +1,47 @@ #ifndef lint static const char RCSid[] = "$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $"; #endif /* ====================================================================== 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 - - -extern char *octname; +#include "pmutil.h" 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) { - /* Clamp lookup bandwidth to total number of photons (minus one, - since density estimate gets extra photon to obtain averaged - radius) */ - sprintf( - errmsg, "clamping density estimate bandwidth to %ld", - pm -> numPhotons - ); - error(WARNING, errmsg); - pm -> minGather = pm -> maxGather = pm -> numPhotons - 1; - } - } -} - - - void cleanUpPmaps (PhotonMap **pmaps) { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) { deletePhotons(pmaps [t]); free(pmaps [t]); } } } - diff --git a/pmutil.h b/pmutil.h new file mode 100644 index 0000000..f9ee4d8 --- /dev/null +++ b/pmutil.h @@ -0,0 +1,49 @@ +/* RCSid $Id$ */ + +/* + ====================================================================== + Auxiliary photon map utilities; stuff that doesn't fit anywhere else! + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) + (c) Fraunhofer Institute for Solar Energy Systems, + supported by the German Research Foundation + (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components", + SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") + (c) Tokyo University of Science, + supported by the JSPS Grants-in-Aid for Scientific Research + (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") + ====================================================================== + + $Id$ +*/ + + +#ifndef _PMAP_UTIL_H + #define _PMAP_UTIL_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) + + + void colorNorm (COLOR c); + /* Normalise colour channels to average of 1 */ + + void cleanUpPmaps (PhotonMap **pmaps); + /* Wrapper to delete all photon maps */ +#endif