diff --git a/Rmakefile b/Rmakefile index 2f72ccc..0b58d48 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,499 +1,543 @@ # 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 pmapamb.o pmapray.o \ pmapopt.o pmapdiag.o pmaptype.o morton.o oococt.o oocsort.o \ oocbuild.o oocnn.o ooccache.o pmutil.o pmaproi.o pmapfilt.o \ wavelet2.o mrgbe.o pmapcontrib.o pmcontrib2.o pmcontrib3.o \ - pmcontrib4.o + pmcontrib4.o pmcontribsort.o PMSRC = pmap.c pmapsrc.c pmapmat.c pmaprand.c pmapio.c pmapdata.c \ pmapdens.c pmapbias.c pmapparm.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 pmaproi.c \ pmapkdt.c pmaptkdt.c pmapooc.c wavelet2.c mrgbe.c \ - pmapcontrib.c pmcontrib2.c pmcontrib3.c pmcontrib4.c + pmapcontrib.c pmcontrib2.c pmcontrib3.c pmcontrib4.c \ + pmcontribsort.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 \ wavelet-test mrgbe-test pmapcontrib-test 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) # Pmap unit tests wavelet-test: wavelet2.c wavelet2.h $(CC) $(CFLAGS) -DWAVELET_TEST -DWAVELET_DBG \ -o wavelet-test -g wavelet2.c mrgbe-test: mrgbe.c mrgbe.h $(CC) $(CFLAGS) -DmRGBE_TEST -o mrgbe-test -g mrgbe.c $(LIBS) pmapcontrib-test: pmapcontrib.c pmapcontrib.h $(CC) $(CFLAGS) -DPMAP_CONTRIB_TEST -o pmapcontrib-test \ -g pmapcontrib.c $(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 \ pmcontrib4.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 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.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 pmaproi.o: pmaproi.c pmaproi.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 otspecial.h ../common/otypes.h mkpmap.o: mkpmap.c pmap.h pmapparm.h pmaptype.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 pmapdata.h ../common/paths.h \ ../common/lookup.h pmapkdt.h pmapfilt.h ../common/fvect.h pmaptkdt.h \ pmutil.h pmapmat.h pmapsrc.h source.h pmapcontrib.h mrgbe.h rcontrib.h \ ../common/platform.h ../common/rtprocess.h ../common/paths.h func.h \ ../common/calcomp.h pmaprand.h pmaproi.h pmapio.h ambient.h \ ../common/resolu.h mrgbe.o: mrgbe.c mrgbe.h wavelet2.o: wavelet2.c wavelet2.h + +pmapcontrib.o: pmapcontrib.c pmapcontrib.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ + ../common/paths.h ../common/lookup.h pmapkdt.h pmapfilt.h \ + ../common/fvect.h pmaptkdt.h wavelet2.h mrgbe.h rcontrib.h \ + ../common/platform.h ../common/rtprocess.h ../common/paths.h func.h \ + ../common/calcomp.h pmapdiag.h pmaprand.h pmapmat.h pmap.h pmutil.h \ + pmapsrc.h source.h otspecial.h ../common/otypes.h + +pmcontrib2.o: pmcontrib2.c pmapcontrib.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ + ../common/paths.h ../common/lookup.h pmapkdt.h pmapfilt.h \ + ../common/fvect.h pmaptkdt.h wavelet2.h mrgbe.h rcontrib.h \ + ../common/platform.h ../common/rtprocess.h ../common/paths.h func.h \ + ../common/calcomp.h pmapdiag.h pmaprand.h pmapmat.h pmap.h pmutil.h \ + pmaproi.h pmapsrc.h source.h pmapio.h otspecial.h + +pmcontrib3.o: pmcontrib3.c pmapcontrib.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ + ../common/paths.h ../common/lookup.h pmapkdt.h pmapfilt.h \ + ../common/fvect.h pmaptkdt.h wavelet2.h mrgbe.h rcontrib.h \ + ../common/platform.h ../common/rtprocess.h ../common/paths.h func.h \ + ../common/calcomp.h pmapdiag.h pmapio.h + +pmcontrib4.o: pmcontrib4.c pmapcontrib.h pmapdata.h ray.h \ + ../common/standard.h ../common/copyright.h ../common/rtio.h \ + ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ + ../common/fvect.h ../common/rterror.h ../common/octree.h \ + ../common/object.h ../common/color.h pmapparm.h pmaptype.h \ + ../common/paths.h ../common/lookup.h pmapkdt.h pmapfilt.h \ + ../common/fvect.h pmaptkdt.h wavelet2.h mrgbe.h rcontrib.h \ + ../common/platform.h ../common/rtprocess.h ../common/paths.h func.h \ + ../common/calcomp.h pmapmat.h pmap.h pmutil.h pmapsrc.h source.h \ + pmaprand.h pmapio.h pmapdiag.h ../common/otypes.h otspecial.h + +pmcontribsort.o: pmcontribsort.c oocsort.h ../common/fvect.h morton.h \ + pmapcontrib.h pmapdata.h ray.h ../common/standard.h \ + ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ + ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ + ../common/rterror.h ../common/octree.h ../common/object.h \ + ../common/color.h pmapparm.h pmaptype.h ../common/paths.h \ + ../common/lookup.h pmapkdt.h pmapfilt.h pmaptkdt.h wavelet2.h mrgbe.h \ + rcontrib.h ../common/platform.h ../common/rtprocess.h ../common/paths.h \ + func.h ../common/calcomp.h diff --git a/oocsort.c b/oocsort.c index 8f2ea8f..58fbf0f 100644 --- a/oocsort.c +++ b/oocsort.c @@ -1,597 +1,582 @@ #ifndef lint static const char RCSid[] = "$Id: oocsort.c,v 2.5 2021/02/10 21:48:50 rschregle Exp $"; #endif /* ========================================================================= N-way out-of-core merge sort for records with 3D keys. Recursively subdivides input into N blocks until these are of sufficiently small size to be sorted in-core according to Z-order (Morton code), then N-way merging them out-of-core using a priority queue as the stack unwinds. 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) ========================================================================== $Id: oocsort.c,v 2.5 2021/02/10 21:48:50 rschregle Exp $ */ #if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) -/* No Windoze support for now */ + /* No Windoze support for now */ -#include "oocsort.h" -#include "morton.h" -#include -#include -#include -#include -#include + #include "oocsort.h" + #include + #include + #include + #include + #include -/* Priority queue node */ -typedef struct { - MortonIdx pri; /* Record's priority (sort key) */ - unsigned blk; /* Block containing record */ -} OOC_SortQueueNode; + /* We resort to instancing the key data as a global variable instead of + passing it to the compare func via qsort_r(), since the latter is a + non-portable GNU extension. */ + static OOC_KeyData keyData; -/* Priority queue */ -typedef struct { - OOC_SortQueueNode *node; - unsigned len, tail; -} OOC_SortQueue; + static int OOC_KeyCompare (const void *rec1, const void *rec2) + /* Comparison function for in-core quicksort */ + { + MortonIdx pri1, pri2; + + if (!rec1 || !rec2) + return 0; + + pri1 = Key2Morton3D(keyData.key(rec1), keyData.bbOrg, + keyData.mortonScale + ); + pri2 = Key2Morton3D(keyData.key(rec2), keyData.bbOrg, + keyData.mortonScale + ); + + if (pri1 < pri2) + return -1; + else if (pri1 > pri2) + return 1; + else + return 0; + } -/* Additional data for qsort() compare function. We resort to instancing - * this as a global variable instead of passing it to the compare func via - * qsort_r(), since the latter is a non-portable GNU extension. */ -typedef struct { - RREAL *(*key)(const void*); /* Callback to access 3D key in record */ - FVECT bbOrg; /* Origin of bbox containing keys */ - RREAL mortonScale; /* Scale for generating Morton code */ -} OOC_KeyData; - -static OOC_KeyData keyData; + int OOC_SortRead (FILE *file, unsigned recSize, char *rec) + /* Read next record from file; return 0 and record in rec on success, + * else -1 */ + { + if (!rec || feof(file) || !fread(rec, recSize, 1, file)) + return -1; -static int OOC_KeyCompare (const void *rec1, const void *rec2) -/* Comparison function for in-core quicksort */ -{ - MortonIdx pri1, pri2; - - if (!rec1 || !rec2) return 0; - - pri1 = Key2Morton3D(keyData.key(rec1), keyData.bbOrg, - keyData.mortonScale - ); - pri2 = Key2Morton3D(keyData.key(rec2), keyData.bbOrg, - keyData.mortonScale - ); - - if (pri1 < pri2) - return -1; - else if (pri1 > pri2) - return 1; - else - return 0; -} + } -static int OOC_SortRead (FILE *file, unsigned recSize, char *rec) -/* Read next record from file; return 0 and record in rec on success, - * else -1 */ -{ - if (!rec || feof(file) || !fread(rec, recSize, 1, file)) - return -1; + int OOC_SortPeek (FILE *file, unsigned recSize, char *rec) + /* Return next record from file WITHOUT advancing file position; + * return 0 and record in rec on success, else -1 */ + { + const unsigned long filePos = ftell(file); - return 0; -} + if (OOC_SortRead(file, recSize, rec)) + return -1; + + /* Workaround; fseek(file, -recSize, SEEK_CUR) causes subsequent + * fread()'s to fail until rewind() */ + rewind(file); + if (fseek(file, filePos, SEEK_SET) < 0) + return -1; + return 0; + } -static int OOC_SortPeek (FILE *file, unsigned recSize, char *rec) -/* Return next record from file WITHOUT advancing file position; - * return 0 and record in rec on success, else -1 */ -{ - const unsigned long filePos = ftell(file); - if (OOC_SortRead(file, recSize, rec)) - return -1; + int OOC_SortWrite (FILE *file, unsigned recSize, const char *rec) + /* Output record to file; return 0 on success, else -1 */ + { + if (!rec || !fwrite(rec, recSize, 1, file)) + return -1; - /* Workaround; fseek(file, -recSize, SEEK_CUR) causes subsequent - * fread()'s to fail until rewind() */ - rewind(file); - if (fseek(file, filePos, SEEK_SET) < 0) - return -1; - - return 0; -} + return 0; + } -static int OOC_SortWrite (FILE *file, unsigned recSize, const char *rec) -/* Output record to file; return 0 on success, else -1 */ -{ - if (!rec || !fwrite(rec, recSize, 1, file)) - return -1; - - return 0; -} + #ifdef DEBUG_OOC_SORT + static int OOC_SortQueueCheck (OOC_SortQueue *pq, unsigned root) + /* Priority queue sanity check */ + { + unsigned kid; + + if (root < pq -> tail) { + if ((kid = (root << 1) + 1) < pq -> tail) { + if (pq -> head [kid].pri < pq -> head [root].pri) + return -1; + else return OOC_SortQueueCheck(pq, kid); + } + + if ((kid = kid + 1) < pq -> tail) { + if (pq -> head [kid].pri < pq -> head [root].pri) + return -1; + else return OOC_SortQueueCheck(pq, kid); + } + } + + return 0; + } + #endif -#ifdef DEBUG_OOC_SORT - static int OOC_SortQueueCheck (OOC_SortQueue *pq, unsigned root) - /* Priority queue sanity check */ + int OOC_SortPush (OOC_SortQueue *pq, MortonIdx pri, unsigned blk) + /* Insert specified block index into priority queue; return block index + * or -1 if queue is full */ { - unsigned kid; + OOC_SortQueueNode *pqn = pq -> head; + unsigned kid, root; - if (root < pq -> tail) { - if ((kid = (root << 1) + 1) < pq -> tail) { - if (pq -> node [kid].pri < pq -> node [root].pri) - return -1; - else return OOC_SortQueueCheck(pq, kid); - } - - if ((kid = kid + 1) < pq -> tail) { - if (pq -> node [kid].pri < pq -> node [root].pri) - return -1; - else return OOC_SortQueueCheck(pq, kid); + if (pq -> tail >= pq -> len) + /* Queue full */ + return -1; + + /* Append node at tail and re-sort */ + kid = pq -> tail++; + + while (kid) { + root = (kid - 1) >> 1; + + /* Compare with parent node and swap if necessary, else terminate */ + if (pri < pqn [root].pri) { + pqn [kid].pri = pqn [root].pri; + pqn [kid].blk = pqn [root].blk; + kid = root; } + else break; + } + + pqn [kid].pri = pri; + pqn [kid].blk = blk; + + #ifdef DEBUG_OOC_SORT + if (OOC_SortQueueCheck(pq, 0) < 0) { + fprintf(stderr, "OOC_Sort: priority queue inconsistency\n"); + return -1; } + #endif - return 0; + return blk; } -#endif -static int OOC_SortPush (OOC_SortQueue *pq, MortonIdx pri, unsigned blk) -/* Insert specified block index into priority queue; return block index - * or -1 if queue is full */ -{ - OOC_SortQueueNode *pqn = pq -> node; - unsigned kid, root; - - if (pq -> tail >= pq -> len) - /* Queue full */ - return -1; - - /* Append node at tail and re-sort */ - kid = pq -> tail++; - - while (kid) { - root = (kid - 1) >> 1; + int OOC_SortPop (OOC_SortQueue *pq) + /* Remove head of priority queue and return its block index + or -1 if queue empty */ + { + OOC_SortQueueNode *pqn = pq -> head; + MortonIdx pri; + unsigned kid, kid2, root = 0, blk, res; - /* Compare with parent node and swap if necessary, else terminate */ - if (pri < pqn [root].pri) { - pqn [kid].pri = pqn [root].pri; - pqn [kid].blk = pqn [root].blk; - kid = root; + if (!pq -> tail) + /* Queue empty */ + return -1; + + res = pqn -> blk; + pri = pqn [--pq -> tail].pri; + blk = pqn [pq -> tail].blk; + + /* Replace head with tail node and re-sort */ + while ((kid = (root << 1) + 1) < pq -> tail) { + /* Compare with smaller kid & swap if necessary, else terminate */ + if ((kid2 = (kid + 1)) < pq -> tail && + pqn [kid2].pri < pqn [kid].pri + ) + kid = kid2; + + if (pri > pqn [kid].pri) { + pqn [root].pri = pqn [kid].pri; + pqn [root].blk = pqn [kid].blk; + } + else break; + + root = kid; } - else break; - } - pqn [kid].pri = pri; - pqn [kid].blk = blk; + pqn [root].pri = pri; + pqn [root].blk = blk; -#ifdef DEBUG_OOC_SORT - if (OOC_SortQueueCheck(pq, 0) < 0) { - fprintf(stderr, "OOC_Sort: priority queue inconsistency\n"); - return -1; + #ifdef DEBUG_OOC_SORT + if (OOC_SortQueueCheck(pq, 0) < 0) { + fprintf(stderr, "OOC_Sort: priority queue inconsistency\n"); + return -1; + } + #endif + + return res; } -#endif - - return blk; -} -static int OOC_SortPop (OOC_SortQueue *pq) -/* Remove head of priority queue and return its block index - or -1 if queue empty */ -{ - OOC_SortQueueNode *pqn = pq -> node; - MortonIdx pri; - unsigned kid, kid2, root = 0, blk, res; - - if (!pq -> tail) - /* Queue empty */ - return -1; - - res = pqn -> blk; - pri = pqn [--pq -> tail].pri; - blk = pqn [pq -> tail].blk; - - /* Replace head with tail node and re-sort */ - while ((kid = (root << 1) + 1) < pq -> tail) { - /* Compare with smaller kid and swap if necessary, else terminate */ - if ((kid2 = (kid + 1)) < pq -> tail && pqn [kid2].pri < pqn [kid].pri) - kid = kid2; - - if (pri > pqn [kid].pri) { - pqn [root].pri = pqn [kid].pri; - pqn [root].blk = pqn [kid].blk; - } - else break; + /* Active subprocess counter updated by parent process; must persist when + * recursing into OOC_SortRecurse(), hence global */ + static unsigned procCnt = 0; + + static int OOC_SortRecurse (FILE *in, unsigned long blkLo, + unsigned long blkHi, FILE *out, unsigned numBlk, + unsigned long maxBlkSize, unsigned numProc, unsigned recSize, + char *sortBuf, OOC_SortQueue *pqueue, const OOC_KeyData *keyData + ) + /* Recursive part of OOC_Sort(). Reads block of records from input file + * within the interval [blkLo, blkHi] and divides it into numBlk blocks + * until the size (in bytes) does not exceed maxBlkSize, then quicksorts + * each block into a temporary file. These files are then mergesorted + * via a priority queue to the output file as the stack unwinds. + * NOTE: Critical sections are prepended with '!!!' + * + * Parameters are as follows: + * in Input file containing unsorted records (assumed to be open) + * blkLo Start of current block in input file, in number of records + * blkHi End of current block in input file, in number of records + * out Output file containing sorted records (assumed to be open) + * numBlk Number of blocks to divide into / merge from + * maxBlkSize Max block size and size of in-core sort buffer, in bytes + * numProc Number of parallel processes for in-core sort + * recSize Size of input records in bytes + * sortBuf Preallocated in-core sort buffer of size maxBlkSize + * pqueue Preallocated priority queue of length numBlk for block merge + * keyData Aggregate data for Morton code generation and comparison + */ + { + const unsigned long blkLen = blkHi - blkLo + 1, + blkSize = blkLen * recSize; + int stat; + + if (!blkLen || blkHi < blkLo) + return 0; - root = kid; - } - - pqn [root].pri = pri; - pqn [root].blk = blk; + if (blkSize > maxBlkSize) { + unsigned long blkLo2 = blkLo, blkHi2 = blkLo, blkSize2; + const double blkLen2 = (double)blkLen / numBlk; + FILE *blkFile [numBlk]; /* Violates C89! */ + char rec [recSize]; /* Ditto */ + MortonIdx pri; + int b, pid; + #ifdef DEBUG_OOC_SORT + unsigned long pqCnt = 0; + #endif + + /* ====================================================== + * Block too large for in-core sort -> divide into numBlk + * subblocks and recurse + * ====================================================== */ + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "OOC_Sort: splitting block [%lu - %lu]\n", + blkLo, blkHi + ); + #endif -#ifdef DEBUG_OOC_SORT - if (OOC_SortQueueCheck(pq, 0) < 0) { - fprintf(stderr, "OOC_Sort: priority queue inconsistency\n"); - return -1; - } -#endif - - return res; -} - - - -/* Active subprocess counter updated by parent process; must persist when - * recursing into OOC_SortRecurse(), hence global */ -static unsigned procCnt = 0; - -static int OOC_SortRecurse (FILE *in, unsigned long blkLo, - unsigned long blkHi, FILE *out, unsigned numBlk, unsigned long maxBlkSize, - unsigned numProc, unsigned recSize, char *sortBuf, OOC_SortQueue *pqueue, - const OOC_KeyData *keyData -) -/* Recursive part of OOC_Sort(). Reads block of records from input file - * within the interval [blkLo, blkHi] and divides it into numBlk blocks - * until the size (in bytes) does not exceed maxBlkSize, then quicksorts - * each block into a temporary file. These files are then mergesorted via a - * priority queue to the output file as the stack unwinds. NOTE: Critical - * sections are prepended with '!!!' - * - * Parameters are as follows: - * in Input file containing unsorted records (assumed to be open) - * blkLo Start of current block in input file, in number of records - * blkHi End of current block in input file, in number of records - * out Output file containing sorted records (assumed to be open) - * numBlk Number of blocks to divide into / merge from - * maxBlkSize Max block size and size of in-core sort buffer, in bytes - * numProc Number of parallel processes for in-core sort - * recSize Size of input records in bytes - * sortBuf Preallocated in-core sort buffer of size maxBlkSize - * pqueue Preallocated priority queue of length numBlk for block merge - * keyData Aggregate data for Morton code generation and comparison - */ -{ - const unsigned long blkLen = blkHi - blkLo + 1, - blkSize = blkLen * recSize; - int stat; - - if (!blkLen || blkHi < blkLo) - return 0; - - if (blkSize > maxBlkSize) { - unsigned long blkLo2 = blkLo, blkHi2 = blkLo, blkSize2; - const double blkLen2 = (double)blkLen / numBlk; - FILE *blkFile [numBlk]; /* Violates C89! */ - char rec [recSize]; /* Ditto */ - MortonIdx pri; - int b, pid; -#ifdef DEBUG_OOC_SORT - unsigned long pqCnt = 0; -#endif - - /* ====================================================== - * Block too large for in-core sort -> divide into numBlk - * subblocks and recurse - * ====================================================== */ - -#ifdef DEBUG_OOC_SORT - fprintf(stderr, "OOC_Sort: splitting block [%lu - %lu]\n", - blkLo, blkHi - ); -#endif + for (b = 0; b < numBlk; b++) { + /* Open temporary file as output for subblock */ + if (!(blkFile [b] = tmpfile())) { + perror("OOC_Sort: failed opening temporary block file"); + return -1; + } + + /* Subblock interval [blkLo2, blkHi2] of size blkSize2 */ + blkHi2 = blkLo + (b + 1) * blkLen2 - 1; + blkSize2 = (blkHi2 - blkLo2 + 1) * recSize; + + if (blkSize2 <= maxBlkSize) { + /* !!! Will be in-core sorted on recursion -> fork kid process, + * !!! but don't fork more than numProc kids, else wait */ + while (procCnt >= numProc && wait(&stat) >= 0) { + if (!WIFEXITED(stat) || WEXITSTATUS(stat)) + return -1; + + procCnt--; + } - for (b = 0; b < numBlk; b++) { - /* Open temporary file as output for subblock */ - if (!(blkFile [b] = tmpfile())) { - perror("OOC_Sort: failed opening temporary block file"); - return -1; - } - - /* Subblock interval [blkLo2, blkHi2] of size blkSize2 */ - blkHi2 = blkLo + (b + 1) * blkLen2 - 1; - blkSize2 = (blkHi2 - blkLo2 + 1) * recSize; - - if (blkSize2 <= maxBlkSize) { - /* !!! Will be in-core sorted on recursion -> fork kid process, - * !!! but don't fork more than numProc kids; wait if necessary */ - while (procCnt >= numProc && wait(&stat) >= 0) { - if (!WIFEXITED(stat) || WEXITSTATUS(stat)) + /* Now fork kid process */ + if (!(pid = fork())) { + /* Recurse on subblocks with new output filehandle */ + if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk, + maxBlkSize, numProc, recSize, sortBuf, pqueue, keyData + )) + exit(-1); + + /* !!! Apparently the parent's tmpfile isn't deleted when + * !!! the child exits (which is what we want), though + * !!! some sources suggest using _Exit() instead; is this + * !!! implementation specific? */ + exit(0); + } + else if (pid < 0) { + perror("OOC_Sort: failed to fork subprocess"); return -1; + } - procCnt--; + #ifdef DEBUG_OOC_FORK + fprintf(stderr, "OOC_Sort: forking kid %d for block %d\n", + procCnt, b + ); + #endif + + /* Parent continues here */ + procCnt++; } - - /* Now fork kid process */ - if (!(pid = fork())) { - /* Recurse on subblocks with new input filehandle; */ + else { + /* No in-core sort on recursion; don't fork */ if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk, maxBlkSize, numProc, recSize, sortBuf, pqueue, keyData )) - exit(-1); - - /* !!! Apparently the parent's tmpfile isn't deleted when the - * !!! child exits (which is what we want), though some - * !!! sources suggest using _Exit() instead; is this - * !!! implementation specific? */ - exit(0); - } - else if (pid < 0) { - perror("OOC_Sort: failed to fork subprocess"); - return -1; + return -1; } -#ifdef DEBUG_OOC_FORK - fprintf(stderr, "OOC_Sort: forking kid %d for block %d\n", - procCnt, b - ); -#endif - - /* Parent continues here */ - procCnt++; + /* Prepare for next block */ + blkLo2 = blkHi2 + 1; } - else { - /* Recurse on subblock; without forking */ - if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk, - maxBlkSize, numProc, recSize, sortBuf, pqueue, keyData - )) - return -1; - } - - /* Prepare for next block */ - blkLo2 = blkHi2 + 1; - } - - /* !!! Wait for any forked processes to terminate */ - while (procCnt && wait(&stat) >= 0) { - if (!WIFEXITED(stat) || WEXITSTATUS(stat)) - return -1; - procCnt--; - } - - /* =============================================================== - * Subblocks are now sorted; prepare priority queue by peeking and - * enqueueing first record from corresponding temporary file - * =============================================================== */ - -#ifdef DEBUG_OOC_SORT - fprintf(stderr, "OOC_Sort: merging block [%lu - %lu]\n", blkLo, blkHi); -#endif + /* !!! Wait for any forked processes to terminate */ + while (procCnt && wait(&stat) >= 0) { + if (!WIFEXITED(stat) || WEXITSTATUS(stat)) + return -1; - for (b = 0; b < numBlk; b++) { -#ifdef DEBUG_OOC_SORT - fseek(blkFile [b], 0, SEEK_END); - if (ftell(blkFile [b]) % recSize) { - fprintf(stderr, "OOC_Sort: truncated record in tmp block " - "file %d\n", b - ); - return -1; + procCnt--; } - - fprintf(stderr, "OOC_Sort: tmp block file %d contains %ld rec\n", - b, ftell(blkFile [b]) / recSize - ); -#endif - rewind(blkFile [b]); + /* =============================================================== + * Subblocks are now sorted; prepare priority queue by peeking and + * enqueueing first record from corresponding temporary file + * =============================================================== + */ - if (OOC_SortPeek(blkFile [b], recSize, rec)) { - fprintf(stderr, "OOC_Sort: error reading from block file\n"); - return -1; - } - - /* Enqueue record along with its Morton code as priority */ - pri = Key2Morton3D(keyData -> key(rec), keyData -> bbOrg, - keyData -> mortonScale + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "OOC_Sort: merging block [%lu - %lu]\n", + blkLo, blkHi ); - - if (OOC_SortPush(pqueue, pri, b) < 0) { - fprintf(stderr, "OOC_Sort: failed priority queue insertion\n"); - return -1; - } - } - - /* ========================================================== - * Subblocks now sorted and priority queue filled; merge from - * temporary files - * ========================================================== */ - - do { - /* Get next node in priority queue, read next record in corresponding - * block, and send to output */ - b = OOC_SortPop(pqueue); - - if (b >= 0) { - /* Priority queue non-empty */ - if (OOC_SortRead(blkFile [b], recSize, rec)) { - /* Corresponding record should still be in the input block */ - fprintf(stderr, - "OOC_Sort: unexpected end reading block file\n" + #endif + + for (b = 0; b < numBlk; b++) { + #ifdef DEBUG_OOC_SORT + fseek(blkFile [b], 0, SEEK_END); + if (ftell(blkFile [b]) % recSize) { + fprintf(stderr, "OOC_Sort: truncated record in tmp block " + "file %d\n", b ); return -1; } + + fprintf(stderr, "OOC_Sort: tmp block file %d contains %ld rec\n", + b, ftell(blkFile [b]) / recSize + ); + #endif + + rewind(blkFile [b]); - if (OOC_SortWrite(out, recSize, rec)) { - fprintf(stderr, "OOC_Sort; error writing output file\n"); + if (OOC_SortPeek(blkFile [b], recSize, rec)) { + fprintf(stderr, "OOC_Sort: error reading from block file\n"); + return -1; + } + + /* Enqueue record along with its Morton code as priority */ + pri = Key2Morton3D(keyData -> key(rec), keyData -> bbOrg, + keyData -> mortonScale + ); + + if (OOC_SortPush(pqueue, pri, b) < 0) { + fprintf(stderr, "OOC_Sort: failed priority queue insertion\n"); return -1; } + } + + /* ========================================================== + * Subblocks now sorted and priority queue filled; merge from + * temporary files + * ========================================================== */ + + do { + /* Get next node in priority queue, read next record in + * corresponding block, and send to output */ + b = OOC_SortPop(pqueue); + + if (b >= 0) { + /* Priority queue non-empty */ + if (OOC_SortRead(blkFile [b], recSize, rec)) { + /* Corresponding record should still be in input block */ + fprintf(stderr, + "OOC_Sort: unexpected end reading block file\n" + ); + return -1; + } -#ifdef DEBUG_OOC_SORT - pqCnt++; -#endif - - /* Peek next record from same block and insert in priority queue */ - if (!OOC_SortPeek(blkFile [b], recSize, rec)) { - /* Block not exhausted */ - pri = Key2Morton3D(keyData -> key(rec), keyData -> bbOrg, - keyData -> mortonScale - ); - - if (OOC_SortPush(pqueue, pri, b) < 0) { - fprintf(stderr, "OOC_Sort: failed priority queue insert\n"); + if (OOC_SortWrite(out, recSize, rec)) { + fprintf(stderr, "OOC_Sort; error writing output file\n"); return -1; } + + #ifdef DEBUG_OOC_SORT + pqCnt++; + #endif + + /* Peek next record from same block and enqueue */ + if (!OOC_SortPeek(blkFile [b], recSize, rec)) { + /* Block not exhausted */ + pri = Key2Morton3D(keyData -> key(rec), keyData -> bbOrg, + keyData -> mortonScale + ); + + if (OOC_SortPush(pqueue, pri, b) < 0) { + fprintf(stderr, + "OOC_Sort: failed priority queue insert\n" + ); + return -1; + } + } } + } while (b >= 0); + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "OOC_Sort: dequeued %lu rec\n", pqCnt); + fprintf(stderr, "OOC_Sort: merged file contains %lu rec\n", + ftell(out) / recSize + ); + #endif + + /* Priority queue now empty -> done; close temporary subblock files, + * causing them to be automatically deleted. */ + for (b = 0; b < numBlk; b++) { + fclose(blkFile [b]); } - } while (b >= 0); - -#ifdef DEBUG_OOC_SORT - fprintf(stderr, "OOC_Sort: dequeued %lu rec\n", pqCnt); - fprintf(stderr, "OOC_Sort: merged file contains %lu rec\n", - ftell(out) / recSize - ); -#endif - - /* Priority queue now empty -> done; close temporary subblock files, - * causing them to be automatically deleted. */ - for (b = 0; b < numBlk; b++) { - fclose(blkFile [b]); + + /* !!! Commit output file to disk before caller reads it; omitting + * !!! this apparently leads to corrupted files (race condition?) + * !!! when the caller tries to read them! */ + fflush(out); + fsync(fileno(out)); } - - /* !!! Commit output file to disk before caller reads it; omitting - * !!! this apparently leads to corrupted files (race condition?) when - * !!! the caller tries to read them! */ - fflush(out); - fsync(fileno(out)); - } - else { - /* ====================================== - * Block is small enough for in-core sort - * ====================================== */ - int ifd = fileno(in), ofd = fileno(out); + else { + /* ====================================== + * Block is small enough for in-core sort + * ====================================== */ + int ifd = fileno(in), ofd = fileno(out); -#ifdef DEBUG_OOC_SORT - fprintf(stderr, - "OOC_Sort: Proc %d (%d/%d) sorting block [%lu - %lu]\n", - getpid(), procCnt, numProc - 1, blkLo, blkHi - ); + #ifdef DEBUG_OOC_SORT + fprintf(stderr, + "OOC_Sort: Proc %d (%d/%d) sorting block [%lu - %lu]\n", + getpid(), procCnt, numProc - 1, blkLo, blkHi + ); + #endif + + /* Atomically seek and read block into in-core sort buffer */ + /* !!! Unbuffered I/O via pread() avoids potential race conditions + * !!! and buffer corruption which can occur with lseek()/fseek() + * !!! followed by read()/fread(). */ + if (pread(ifd, sortBuf, blkSize, blkLo * recSize) != blkSize) { + perror("OOC_Sort: error reading from input file"); + return -1; + } -#endif + /* Quicksort block in-core and write to output file */ + qsort(sortBuf, blkLen, recSize, OOC_KeyCompare); - /* Atomically seek and read block into in-core sort buffer */ - /* !!! Unbuffered I/O via pread() avoids potential race conditions - * !!! and buffer corruption which can occur with lseek()/fseek() - * !!! followed by read()/fread(). */ - if (pread(ifd, sortBuf, blkSize, blkLo * recSize) != blkSize) { - perror("OOC_Sort: error reading from input file"); - return -1; + if (write(ofd, sortBuf, blkSize) != blkSize) { + perror("OOC_Sort: error writing to block file"); + return -1; + } + + fsync(ofd); + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "OOC_Sort: proc %d wrote %ld records\n", + getpid(), lseek(ofd, 0, SEEK_END) / recSize + ); + #endif } - /* Quicksort block in-core and write to output file */ - qsort(sortBuf, blkLen, recSize, OOC_KeyCompare); + return 0; + } + + + + int OOC_Sort (FILE *in, FILE *out, unsigned numBlk, + unsigned long blkSize, unsigned numProc, unsigned recSize, + FVECT bbOrg, RREAL bbSize, RREAL *(*key)(const void*) + ) + /* Sort records in inFile and append to outFile by subdividing inFile into + * small blocks, sorting these in-core, and merging them out-of-core via a + * priority queue. + * + * This is implemented as a recursive (numBlk)-way sort; the input is + * successively split into numBlk smaller blocks until these are of size <= + * blkSize, i.e. small enough for in-core sorting, then merging the sorted + * blocks as the stack unwinds. The in-core sort is parallelised over + * numProx processes. + * + * Parameters are as follows: + * inFile Opened input file containing unsorted records + * outFile Opened output file containing sorted records + * numBlk Number of blocks to divide into / merge from + * blkSize Max block size and size of in-core sort buffer, in bytes + * numProc Number of parallel processes for in-core sort + * recSize Size of input records in bytes + * bbOrg Origin of bounding box containing keys for Morton code + * bbSize Extent of bounding box containing keys for Morton code + * key Callback to access 3D coords from records for Morton code + */ + { + unsigned long numRec; + OOC_SortQueue pqueue; + char *sortBuf = NULL; + int stat; + + if (numBlk < 1) + numBlk = 1; - if (write(ofd, sortBuf, blkSize) != blkSize) { - perror("OOC_Sort: error writing to block file"); + /* Open input file & get size in number of records */ + if (fseek(in, 0, SEEK_END) < 0) { + fputs("OOC_Sort: failed seek in input file\n", stderr); return -1; } - fsync(ofd); - -#ifdef DEBUG_OOC_SORT - fprintf(stderr, "OOC_Sort: proc %d wrote %ld records\n", - getpid(), lseek(ofd, 0, SEEK_END) / recSize - ); -#endif - } - - return 0; -} - - - -int OOC_Sort (FILE *in, FILE *out, unsigned numBlk, unsigned long blkSize, - unsigned numProc, unsigned recSize, FVECT bbOrg, RREAL bbSize, - RREAL *(*key)(const void*) -) -/* Sort records in inFile and append to outFile by subdividing inFile into - * small blocks, sorting these in-core, and merging them out-of-core via a - * priority queue. - * - * This is implemented as a recursive (numBlk)-way sort; the input is - * successively split into numBlk smaller blocks until these are of size <= - * blkSize, i.e. small enough for in-core sorting, then merging the sorted - * blocks as the stack unwinds. The in-core sort is parallelised over - * numProx processes. - * - * Parameters are as follows: - * inFile Opened input file containing unsorted records - * outFile Opened output file containing sorted records - * numBlk Number of blocks to divide into / merge from - * blkSize Max block size and size of in-core sort buffer, in bytes - * numProc Number of parallel processes for in-core sort - * recSize Size of input records in bytes - * bbOrg Origin of bounding box containing record keys for Morton code - * bbSize Extent of bounding box containing record keys for Morton code - * key Callback to access 3D coords from records for Morton code - */ -{ - unsigned long numRec; - OOC_SortQueue pqueue; - char *sortBuf = NULL; - int stat; - - if (numBlk < 1) - numBlk = 1; + if (!(numRec = ftell(in) / recSize)) { + fputs("OOC_Sort: empty input file\n", stderr); + return -1; + } - /* Open input file & get size in number of records */ - if (fseek(in, 0, SEEK_END) < 0) { - fputs("OOC_Sort: failed seek in input file\n", stderr); - return -1; - } - - if (!(numRec = ftell(in) / recSize)) { - fputs("OOC_Sort: empty input file\n", stderr); - return -1; - } + /* Allocate & init priority queue */ + if (!(pqueue.head = calloc(numBlk, sizeof(OOC_SortQueueNode)))) { + fputs("OOC_Sort: failure allocating priority queue\n", stderr); + return -1; + } + pqueue.tail = 0; + pqueue.len = numBlk; - /* Allocate & init priority queue */ - if (!(pqueue.node = calloc(numBlk, sizeof(OOC_SortQueueNode)))) { - fputs("OOC_Sort: failure allocating priority queue\n", stderr); - return -1; - } - pqueue.tail = 0; - pqueue.len = numBlk; + /* Allocate in-core sort buffa */ + if (!(sortBuf = malloc(blkSize))) { + fprintf(stderr, "OOC_Sort: failure allocating in-core sort buffer"); + return -1; + } + + /* Set up key data to pass to qsort()'s comparison func */ + keyData.key = key; + keyData.mortonScale = MORTON3D_MAX / bbSize; + VCOPY(keyData.bbOrg, bbOrg); + + stat = OOC_SortRecurse(in, 0, numRec - 1, out, numBlk, blkSize, + numProc, recSize, sortBuf, &pqueue, &keyData + ); - /* Allocate in-core sort buffa */ - if (!(sortBuf = malloc(blkSize))) { - fprintf(stderr, "OOC_Sort: failure allocating in-core sort buffer"); - return -1; + /* Cleanup */ + free(pqueue.head); + free(sortBuf); + + return stat; } - - /* Set up key data to pass to qsort()'s comparison func */ - keyData.key = key; - keyData.mortonScale = MORTON3D_MAX / bbSize; - VCOPY(keyData.bbOrg, bbOrg); - - stat = OOC_SortRecurse(in, 0, numRec - 1, out, numBlk, blkSize, numProc, - recSize, sortBuf, &pqueue, &keyData - ); - - /* Cleanup */ - free(pqueue.node); - free(sortBuf); - - return stat; -} #endif /* NIX / PMAP_OOC */ diff --git a/oocsort.h b/oocsort.h index ef0c475..bc7a336 100644 --- a/oocsort.h +++ b/oocsort.h @@ -1,44 +1,88 @@ /* ====================================================================== Header for N-way out-of-core merge sort for records with 3D keys. 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) ====================================================================== $Id: oocsort.h,v 2.3 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef OOC_SORT_H #define OOC_SORT_H #include "fvect.h" - #include - + #include "morton.h" + #include + + + + /* Priority queue node */ + typedef struct { + MortonIdx pri; /* Record's priority (sort key) */ + unsigned blk; /* Block containing record */ + } OOC_SortQueueNode; + + /* Priority queue */ + typedef struct { + OOC_SortQueueNode *head; + unsigned len, tail; + } OOC_SortQueue; + + /* Additional data for key comparison function.*/ + typedef struct { + RREAL *(*key)(const void*); /* Callback to access 3D key in record */ + FVECT bbOrg; /* Origin of bbox containing keys */ + RREAL mortonScale; /* Scale for generating Morton code */ + } OOC_KeyData; + + + + /* Read next record from file; return 0 and record in rec on success, + * else -1 */ + int OOC_SortRead (FILE *file, unsigned recSize, char *rec); + + /* Return next record from file WITHOUT advancing file position; + * return 0 and record in rec on success, else -1 */ + int OOC_SortPeek (FILE *file, unsigned recSize, char *rec); + + /* Output record to file; return 0 on success, else -1 */ + int OOC_SortWrite (FILE *file, unsigned recSize, const char *rec); + + /* Insert specified block index into priority queue; return block index + * or -1 if queue is full */ + int OOC_SortPush (OOC_SortQueue *pq, MortonIdx pri, unsigned blk); + + /* Remove head of priority queue and return its block index or -1 + if queue empty */ + int OOC_SortPop (OOC_SortQueue *pq); /* Sort records in inFile and append to outFile by subdividing inFile * into small blocks, sorting these in-core, and merging them out-of-core * via a priority queue. * * This is implemented as a recursive (numBlk)-way sort; the input is * successively split into numBlk smaller blocks until these are of size * <= blkSize, i.e. small enough for in-core sorting, then merging the * sorted blocks as the stack unwinds. The in-core sort is parallelised * over numProc processes. * * Parameters are as follows: * inFile Opened input file containing unsorted records * outFile Opened output file containing sorted records * numBlk Number of blocks to divide into / merge from * blkSize Max block size and size of in-core sort buffer, in bytes * numProc Number of parallel processes for in-core sort * recSize Size of input records in bytes * bbOrg Origin of bounding box containing record keys for Morton code * bbSize Extent of bounding box containing record keys for Morton code * key Callback to access 3D coords from records for Morton code - */ - int OOC_Sort (FILE *inFile, FILE *outFile, unsigned numBlk, - unsigned long blkSize, unsigned numProc, unsigned recSize, - FVECT bbOrg, RREAL bbSize, RREAL *(*key)(const void*)); + */ + int OOC_Sort (FILE *inFile, FILE *outFile, + unsigned numBlk, unsigned long blkSize, unsigned numProc, + unsigned recSize, FVECT bbOrg, RREAL bbSize, + RREAL *(*key)(const void*)); #endif + diff --git a/pmapdata.c b/pmapdata.c index 122f883..6847de0 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,851 +1,852 @@ #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 stats */ setcolor(pmap -> photonFlux, 0, 0, 0); pmap -> CoG [0] = pmap -> CoG [1] = pmap -> CoG [2] = 0; pmap -> CoGdist = 0; pmap -> maxDist2 = 0; /* 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]; /* Init precomputed contrib photon stuff */ pmap -> preCompContribTab = NULL; pmap -> contribHeap = NULL; + pmap -> waveletFile = NULL; /* 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; /* 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 (to consolidate contribution sources after photon distribution completes) */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { case PMAP_TYPE_CONTRIB: if (pmap -> lastContribSrc.srcIdx < -1) { /* HACK: Contrib photon being precomputed; set contribution source from index passed in ray */ photon.aux.contribSrc = ray -> rsrc; } else { /* HACK: Contrib photon before precomputation (i.e. in forward pass); set contribution source from last index in contrib source array. Note the contribution source bin has already been set by newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; /* Photon will be stored; set contribution source index, * thereby marking it as having spawned photon(s) */ 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]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* 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 /* Contrib stuff */ if (isContribPmap(pmap) && pmap -> preCompContribTab) { lu_done(pmap -> preCompContribTab); free(pmap -> preCompContribTab); pmap -> preCompContribTab = NULL; } } diff --git a/pmapdata.h b/pmapdata.h index 8e54d0e..3c35578 100644 --- a/pmapdata.h +++ b/pmapdata.h @@ -1,432 +1,434 @@ /* 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 { /* 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 /* 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) /* Bias compensation history node */ typedef struct { COLOR irrad; float weight; } PhotonBiasCompNode; /* Forward declaration */ struct PhotonMap; 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 */ LUTAB *preCompContribTab; /* LUT for per-modifier precomp. contributions including constituent photon maps */ - FILE *contribHeap; /* Out-of-core heap containing + FILE *contribHeap, /* Out-of-core heap containing unsorted precomputed contribution photon bins prior to construction of store */ + *waveletFile; /* Sorted precomputed contribution + photon bins after construction */ char contribHeapFname [sizeof(PMAP_TMPFNAME)]; /* ================================================================ * 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]) #else #define lightflowPmap NULL #define transLightFlowPmap NULL #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/pmapooc.c b/pmapooc.c index 4438275..92880b5 100644 --- a/pmapooc.c +++ b/pmapooc.c @@ -1,349 +1,358 @@ /* ====================================================================== Photon map interface to out-of-core octree Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "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: pmapooc.c,v 1.7 2021/03/22 23:00:00 rschregle Exp $ */ #include "pmapdata.h" /* Includes pmapooc.h */ #include "source.h" #include "otspecial.h" #include "oocsort.h" #include "oocbuild.h" /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ /* Returns photon position as sorting key for OOC_Octree & friends (notably * for Morton code generation). * !!! XXX: Uses type conversion from float via TEMPORARY storage; * !!! THIS IS NOT THREAD SAFE! * !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */ RREAL *OOC_PhotonKey (const void *p) { static FVECT photonPos; /* Temp storage for type conversion */ VCOPY(photonPos, ((Photon*)p) -> pos); return photonPos; } #ifdef DEBUG_OOC static void OOC_CheckKeys (FILE *file, const OOC_Octree *oct) { Photon p, lastp; RREAL *key; MortonIdx zIdx, lastzIdx = 0; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file) > 0) { key = OOC_PhotonKey(&p); zIdx = OOC_KEY2MORTON(key, oct); if (zIdx < lastzIdx) error(INTERNAL, "photons not sorted"); if (zIdx == lastzIdx) { sprintf(errmsg, "identical key %021ld " "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]", zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } lastzIdx = zIdx; memcpy(&lastp, &p, sizeof(p)); } } #endif void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc) { FILE *leafFile; char leafFname [1024]; FVECT d, octOrg; int i; RREAL octSize = 0; /* Determine octree size and origin from pmap extent and init octree */ VCOPY(octOrg, pmap -> minPos); VSUB(d, pmap -> maxPos, pmap -> minPos); for (i = 0; i < 3; i++) if (octSize < d [i]) octSize = d [i]; if (octSize < FTINY) error(INTERNAL, "zero octree size in OOC_BuildPhotonMap"); /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1 ); if (!(leafFile = fopen(leafFname, "w+b"))) error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Sorting photons by Morton code...\n"); #endif - - /* Sort photons in heapfile by Morton code and write to leaf file */ - if (OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, - numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey) + + /* Sort photons in heapfile by Morton code and write to leaf file; + redirect to contrib sorting routine for contribution photons! */ + if (isContribPmap(pmap) + ? contribPhotonSort(pmap -> heap, pmap -> contribHeap, + leafFile, pmap -> waveletFile, pmap -> numPhotons, + pmap -> preCompContrib -> nCompressedBins * sizeof(mRGBE), + octOrg, octSize, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, + numProc, OOC_PhotonKey + ) + : OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, + numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey + ) ) error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap"); /* Init and build octree */ OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey, leafFile ); #ifdef DEBUG_OOC eputs("Checking leaf file consistency...\n"); OOC_CheckKeys(leafFile, &pmap -> store); eputs("Building out-of-core octree...\n"); #endif if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH)) error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Checking out-of-core octree consistency...\n"); if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), octOrg, octSize, 0 )) error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri"); #endif } /* PHOTON MAP I/O ROUTINES ------------------------------------------ */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out) { return OOC_SaveOctree(&pmap -> store, out); } int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *nodeFile) { FILE *leafFile; char leafFname [1024]; /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1 ); if (!(leafFile = fopen(leafFname, "r"))) error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons"); if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile)) return -1; #ifdef DEBUG_OOC /* Check octree for consistency */ if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), pmap -> store.org, pmap -> store.size, 0 )) return -1; #endif return 0; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ void OOC_InitFindPhotons (struct PhotonMap *pmap) { if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1, sizeof(Photon) )) error(SYSTEM, "can't allocate photon search queue"); #ifdef PMAP_PATHFILT initPhotonPaths(pmap); #endif } static void OOC_InitPhotonCache (struct PhotonMap *pmap) /* Initialise OOC photon cache */ { static char warn = 1; if (!pmap -> store.cache && !pmap -> numDensity) { if (pmapCacheSize > 0) { const unsigned pageSize = pmapCachePageSize * pmap -> maxGather, numPages = pmapCacheSize / pageSize; /* Allocate & init I/O cache in octree */ pmap -> store.cache = malloc(sizeof(OOC_Cache)); if (!pmap -> store.cache || OOC_CacheInit(pmap -> store.cache, numPages, pageSize, sizeof(Photon) )) { error(SYSTEM, "failed OOC photon map cache init"); } } else if (warn) { error(WARNING, "OOC photon map cache DISABLED"); warn = 0; } } } int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; #ifdef PMAP_PATHFILT OOC_SearchAttribCallback paths, *pathsPtr = NULL; #endif float res, n [3]; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.pos = pos; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; /* Get sought light source modifier from pmap -> lastContribSrc if precomputing contribution photon */ filtState.srcMod = isContribPmap(pmap) ? findmaterial(source [pmap -> lastContribSrc.srcIdx].so) : NULL; #ifdef PMAP_PATHFILT /* Set up path ID callback to filter for photons from unique paths */ paths.state = pmap; paths.findFunc = findPhotonPath; paths.addFunc = addPhotonPath; paths.delFunc = deletePhotonPath; paths.checkFunc = checkPhotonPaths; pathsPtr = &paths; resetPhotonPaths(pmap); res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, pathsPtr, &pmap -> squeue, pmap -> maxDist2 ); #else res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, NULL, &pmap -> squeue, pmap -> maxDist2 ); #endif /* PMAP_PATHFILT */ if (res < 0) error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons"); /* Get (maximum distance)^2 from search queue head */ pmap -> maxDist2 = pmap -> squeue.node [0].dist2; /* Return success or failure (empty queue => none found) */ return pmap -> squeue.tail ? 0 : -1; } int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, const FVECT norm, Photon *photon ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; float n [3], maxDist2; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filtState.srcMod = NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; maxDist2 = OOC_Find1Nearest(&pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size,pos, &filt, photon, pmap -> maxDist2 ); if (maxDist2 < 0) error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon"); if (maxDist2 >= pmap -> maxDist2) /* No photon found => failed */ return -1; else { /* Set photon distance => success */ pmap -> maxDist2 = maxDist2; return 0; } } /* PHOTON ACCESS ROUTINES ------------------------------------------ */ int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon) { return OOC_GetData(&pmap -> store, idx, photon); } Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { return OOC_GetNearest(squeue, idx); } PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap) { return 0; } diff --git a/pmcontrib3.c b/pmcontrib3.c index cb9d069..a2193bd 100644 --- a/pmcontrib3.c +++ b/pmcontrib3.c @@ -1,192 +1,205 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines build and save precomputed contribution photon maps, and are used by 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$ */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmapio.h" #include #if NIX #define __USE_XOPEN_EXTENDED /* Enables nftw() under Linux */ #include #endif #if NIX static int ftwFunc (const char *path, const struct stat *statData, int typeFlags, struct FTW *ftwData ) /* Callback for nftw() to delete a directory entry */ { return remove(path); /* = unlink()/rmdir() depending on type */ } #endif int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p) /* Build per-modifier precomputed photon map. Returns 0 on success. */ { PreComputedContrib *preCompContrib = (PreComputedContrib*)( preCompContribNode -> data ); PhotonMap *preCompContribPmap = preCompContrib -> pmap; char dirName [1024]; struct stat dirStat; /* Flush photon/contrib heaps */ flushPhotonHeap(preCompContribPmap); fflush(preCompContribPmap -> contribHeap); if (verbose) { sprintf(errmsg, "Building precomputed contribution photon map " "for modifier %s\n", preCompContribNode -> key ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Set output directory for filename */ snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, preCompContribPmap -> fileName ); /* Allocate & set output filenames to subdirectory and modifier */ preCompContribPmap -> fileName = malloc(strlen(dirName) + strlen(preCompContribNode -> key) + strlen(PMAP_CONTRIB_FILE) ); preCompContrib -> waveletFname = malloc(strlen(dirName) + strlen(preCompContribNode -> key) + strlen(PMAP_CONTRIB_WAVELETFILE) ); if (!preCompContribPmap -> fileName || !preCompContrib -> waveletFname) error(SYSTEM, "out of memory allocating filename in " "buildPreCompContribPmap()" ); snprintf(preCompContribPmap -> fileName, 1028, PMAP_CONTRIB_FILE, dirName, preCompContribNode -> key ); snprintf(preCompContrib -> waveletFname, 1029, PMAP_CONTRIB_WAVELETFILE, dirName, preCompContribNode -> key ); /* Cleanup previous directory contents if necessary */ if (!stat(dirName, &dirStat)) { /* File/dir exists */ if (S_ISREG(dirStat.st_mode)) /* Found regular file; delete and skip rest */ unlink(dirName); else if (S_ISDIR(dirStat.st_mode)) { /* Found subdirectory; delete its contents, skipping symlinks */ #if NIX if (nftw(dirName, ftwFunc, 1, FTW_DEPTH | FTW_MOUNT | FTW_PHYS) < 0 ) { sprintf(errmsg, "failed cleanup of output directory %s", dirName ); error(SYSTEM, errmsg); } #else /* Apparently there's no nftw() under Wind0z, so just skip * cleanup. Tuff luck, Wind0z Weenie! There's probably some * replacement for this, but we couldn't be buggered... */ #endif } else { /* Found neither regular file nor directory; whatever it is, * just stuff it and quit! */ sprintf(errmsg, "cannot remove existing output file %s", dirName ); error(SYSTEM, errmsg); } } /* Create new directory for per-modifier contribution photon maps */ if (mkdir(dirName, S_IFDIR | S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH) < 0 ) { sprintf(errmsg, "error creating output directory %s", dirName); error(SYSTEM, errmsg); } - /* Rebuild underlying data structure, destroying heap */ - buildPhotonMap(preCompContribPmap, NULL, NULL, 1); + /* Open wavelet coefficient file; this is where the sorted contribs + are written to by contribPhotonSort() (via buildPhotonMap()) */ + if (!(preCompContribPmap -> waveletFile = fopen( + preCompContrib -> waveletFname, "wb" + )) + ) { + sprintf(errmsg, "can't open wavelet coefficient file %s", + preCompContrib -> waveletFname + ); + error(SYSTEM, errmsg); + } + + /* Rebuild underlying data structure, destroying heap, + and sort contribs into wavelet coeff file */ + buildPhotonMap(preCompContribPmap, NULL, NULL, 1); /* Free primary and transposed wavelet matrices; no longer needed */ free(preCompContrib -> waveletMatrix); freeWaveletMatrix(preCompContrib -> tWaveletMatrix, preCompContrib -> l); /* Free thresholded coefficients; no longer needed */ free(preCompContrib -> threshCoeffs); preCompContrib -> waveletMatrix = preCompContrib -> tWaveletMatrix = NULL; preCompContrib -> threshCoeffs = NULL; return 0; } static int savePreCompContrib (const LUENT *preCompContribNode, void *p) { PreComputedContrib *preCompContrib = (PreComputedContrib*)( preCompContribNode -> data ); PhotonMapArgz *argz = (PhotonMapArgz*)p; PhotonMap *preCompContribPmap = preCompContrib -> pmap; FILE *waveletFile; /* Save contrib photon map for current modifier; savePhotonMap will not * recurse and call saveContribPhotonMap() again because * preCompContribPmap -> preCompContribTab == NULL */ savePhotonMap(preCompContribPmap, preCompContribPmap -> fileName, argz -> argc, argz -> argv ); #if 0 /* Save wavelet coefficients */ if (!(waveletFile = fopen(preCompContrib -> waveletFname, "wb"))) { sprintf(errmsg, "can't open wavelet coefficient file %s", preCompContrib -> waveletFname ); error(SYSTEM, errmsg); } #endif return 0; } void saveContribPhotonMap (const PhotonMap *pmap, int argc, char **argv) { PhotonMapArgz argz = {argc, argv}; lu_doall(pmap -> preCompContribTab, savePreCompContrib, &argz); } diff --git a/pmcontribsort.c b/pmcontribsort.c new file mode 100644 index 0000000..67fc9ff --- /dev/null +++ b/pmcontribsort.c @@ -0,0 +1,511 @@ +#ifndef lint +static const char RCSid[] = "$Id$"; +#endif + + +/* + ========================================================================= + Support routines for N-way out-of-core merge sort of precomputed photon + contributions. These routines are used adapted from OOC_Sort() to sort + binned contributions in the same order as their associated precomputed + photons. + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation + (SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") + ========================================================================== + + $Id$ +*/ + + +#if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) + /* No Windoze 'cos no fork() */ + + #include "oocsort.h" + #include "pmutil.h" + #include "pmapcontrib.h" + #include + #include + #include + #include + #include + + #define PHOTONSIZE sizeof(Photon) + #define PMAPCONTRIB_KEY(p, kdat) Key2Morton3D( \ + kdat -> key(p), kdat -> bbOrg, kdat -> mortonScale \ + ) + + + extern RREAL *OOC_PhotonKey (const void *p); + + + static void contribPhotonQuickSort (Photon *photons, char *contribs, + unsigned contribSize, const OOC_KeyData *keyData, + unsigned long left, unsigned long right + ) + /* Partition photons and their binned contributions using a quicksort + * algorithm. Returns index to pivot's position after partitioning. */ + { + const unsigned long m = (right - left + 1) >> 1; /* Pivot in mid */ + unsigned long l = left, r = right; + const MortonIdx mKey = PMAPCONTRIB_KEY(photons + m, keyData); + char tmp [max(PHOTONSIZE, contribSize)]; + + if (left < right) { + /* l & r converge, swapping photons (and their corresponding + * binned contributions) out of order with respect to pivot at + * photons [m]. The convergence point l = r is the pivot's final + * position */ + do { + while (l <= r && PMAPCONTRIB_KEY(photons + l, keyData) < mKey) + l++; + + while (r > l && PMAPCONTRIB_KEY(photons + r, keyData) >= mKey) + r--; + + /* Photons out of order; swap them and their corresponding + * binned contributions */ + /* TODO: swap photon indices instead of photons themselves? */ + memcpy(tmp, photons + l, PHOTONSIZE); + memcpy(photons + l, photons + r, PHOTONSIZE); + memcpy(photons + r, tmp, PHOTONSIZE); + memcpy(tmp, contribs + l, contribSize); + memcpy(contribs + l, contribs + r, contribSize); + memcpy(contribs + r, tmp, contribSize); + } while (l < r); + + /* Now l = r; swap this photon with pivot, ditto for contribs */ + memcpy(tmp, photons + r, PHOTONSIZE); + memcpy(photons + r, photons + m, PHOTONSIZE); + memcpy(photons + m, tmp, PHOTONSIZE); + memcpy(tmp, contribs + r, contribSize); + memcpy(contribs + r, contribs + m, contribSize); + memcpy(contribs + m, tmp, contribSize); + + /* Recurse in left & right partitions */ + contribPhotonQuickSort(photons, contribs, contribSize, keyData, + left, m - 1 + ); + contribPhotonQuickSort(photons, contribs, contribSize, keyData, + m + 1, right + ); + } + } + + + + /* Active subprocess counter updated by parent process; must persist when + * recursing into contribPhotonSortRecurse(), hence global */ + static unsigned procCnt = 0; + + static int contribPhotonSortRecurse ( + FILE *photonIn, FILE *contribIn, FILE* photonOut, FILE* contribOut, + unsigned contribSize, unsigned numBlk, unsigned numProc, + Photon *photonBuf, char *contribBuf, unsigned long bufLen, + const OOC_KeyData *keyData, unsigned long blkLo, unsigned long blkHi + ) + /* Recursive part of contribPhotonSort(). Reads block of photons and + * their binned contributions from input files photonIn resp. + * contribOut, within the interval [blkLo, blkHi], and divides them into + * numBlk blocks until the number of photons/contributions in a block + * does not exceed bufLen, then quicksorts each block in-core into + * temporary files. These files are then mergesorted via a priority + * queue to the specified output files as the stack unwinds. + * NOTE: Critical sections are prepended with '!!!' + * + * Parameters are as follows: + * photonIn Opened primary input file containing unsorted photons; + * their Morton codes determine the ordering + * contribIn Opened secondary input file containing unsorted binned + * contributions size contribSize + * photonOut Opened primary output file for sorted photons + * contribOut Opened secondary output file for sorted binned + * contributions ordered by corresponding photons + * contribSize Size of binned contributions per photon, in bytes + * numBlk Number of blocks to divide into / merge from + * numProc Number of parallel processes for in-core sort + * photonBuf Preallocated in-core sort buffer for maxRec photons + * contribBuf Preallocated in-core sort buffer for maxRec contribs + * bufLen In-core sort buffer length, in number of photons/contribs + * keyData Aggregate data for Morton code generation as sort keys + * blkLo Start of current block in input files, in number of photons + * blkHi End of current block in input files, in number of photons + */ + { + const unsigned long blkLen = blkHi - blkLo + 1; + + if (!blkLen || blkHi < blkLo) + return 0; + + if (blkLen > bufLen) { + unsigned long blkLo2 = blkLo, blkHi2 = blkLo; + const double blkLen2 = (double)blkLen / numBlk; + FILE *pBlkFile [numBlk], *cBlkFile [numBlk]; + Photon photon; + char contrib [contribSize]; + MortonIdx pri; + int b, pid, stat; + OOC_SortQueueNode pqueueNodes [numBlk]; + OOC_SortQueue pqueue; + #ifdef DEBUG_OOC_SORT + unsigned long pqCnt = 0; + #endif + + /* ====================================================== + * Block too large for in-core sort -> divide into numBlk + * subblocks and recurse + * ====================================================== */ + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, + "*DBG* contribPhotonSort: splitting block [%lu - %lu]\n", + blkLo, blkHi + ); + #endif + + for (b = 0; b < numBlk; b++) { + /* Open temporary file as output for subblock */ + if (!(pBlkFile [b] = tmpfile()) || !(cBlkFile [b] = tmpfile())) { + perror("contribPhotonSort: failed to open temp block file"); + return -1; + } + + /* Set new subblock interval [blkLo2, blkHi2] */ + blkHi2 = blkLo + (b + 1) * blkLen2 - 1; + + if (blkHi2 - blkLo2 + 1 <= bufLen) { + /* !!! Subblock [blkLo2, blkHi2] will be sorted in-core on + * !!! recursion, so fork kid process. But don't fork more + * !!! than numProc kids; wait if necessary */ + while (procCnt >= numProc && wait(&stat) >= 0) { + /* Kid process terminated; check status, update counter */ + if (!WIFEXITED(stat) || WEXITSTATUS(stat)) + return -1; + procCnt--; + } + + /* Now fork kid process */ + if (!(pid = fork())) { + /* Recurse on subblock with temp output files */ + if (contribPhotonSortRecurse(photonIn, contribIn, + pBlkFile [b], cBlkFile [b], contribSize, + numBlk, numProc, photonBuf, contribBuf, bufLen, + keyData, blkLo2, blkHi2 + )) + exit(-1); + + /* !!! Kid exits here. Apparently the parent's tmpfile + * !!! isn't deleted (which is what we want), though + * !!! some sources suggest using _Exit() instead; + * !!! is this implementation specific? Anyway, we + * !!! leave cleanup to the parent (see below) */ + exit(0); + } + else if (pid < 0) { + perror("contribPhotonSort: failed to fork subprocess"); + return -1; + } + + #ifdef DEBUG_OOC_FORK + fprintf(stderr, + "*DBG* contribPhotonSort: forking kid %d for block %d\n", + procCnt, b + ); + #endif + + /* Parent continues here */ + procCnt++; + } + else { + /* Subblock will NOT be sorted in-core on recursion, + so DON'T fork */ + if (contribPhotonSortRecurse(photonIn, contribIn, + pBlkFile [b], cBlkFile [b], contribSize, + numBlk, numProc, photonBuf, contribBuf, bufLen, + keyData, blkLo2, blkHi2 + )) + return -1; + } + + /* Update offset for next subblock */ + blkLo2 = blkHi2 + 1; + } + + /* !!! Wait for any forked processes to terminate */ + while (procCnt && wait(&stat) >= 0) { + /* Check status, bail out on error or update process counter */ + if (!WIFEXITED(stat) || WEXITSTATUS(stat)) + return -1; + procCnt--; + } + + /* ============================================================= + * Subblocks are now sorted; prepare priority queue + * ============================================================= */ + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, + "*DBG* contribPhotonSort: merging block [%lu - %lu]\n", + blkLo, blkHi + ); + #endif + + /* Init priority queue */ + pqueue.head = pqueueNodes; + pqueue.tail = 0; + pqueue.len = numBlk; + + /* Fill priority queue by peeking and enqueueing first photon + from each subblock */ + for (b = 0; b < numBlk; b++) { + #ifdef DEBUG_OOC_SORT + fseek(pBlkFile [b], 0, SEEK_END); + fseek(cBlkFile [b], 0, SEEK_END); + + if (ftell(pBlkFile [b]) % PHOTONSIZE || + ftell(cBlkFile [b]) % contribSize + ) { + fprintf(stderr, "contribPhotonSort: truncated record " + "in temp block file %d\n", b + ); + return -1; + } + + fprintf(stderr, "contribPhotonSort: temp block file %d " + "contains %lu photons, %lu contribs\n", b, + ftell(pBlkFile [b]) / PHOTONSIZE, + ftell(cBlkFile [b]) / contribSize + ); + #endif + + rewind(pBlkFile [b]); + rewind(cBlkFile [b]); + + if (OOC_SortPeek(pBlkFile [b], PHOTONSIZE, (char*)&photon) || + OOC_SortPeek(cBlkFile [b], contribSize, contrib) + ) { + perror("contribPhotonSort: error reading block file"); + return -1; + } + + /* Enqueue photon's Morton code as priority */ + pri = PMAPCONTRIB_KEY(&photon, keyData); + + if (OOC_SortPush(&pqueue, pri, b) < 0) { + perror("contribPhotonSort: failed priority queue insertion"); + return -1; + } + } + + /* ========================================================== + * Subblocks now sorted and priority queue filled; + * merge subblocks from temporary files + * ========================================================== */ + + do { + /* Get head of priority queue, read next photon/contrib in + * corresponding subblock, and send to output */ + if ((b = OOC_SortPop(&pqueue)) >= 0) { + /* Priority queue not empty */ + if (OOC_SortRead(pBlkFile [b], PHOTONSIZE, (char*)&photon) || + OOC_SortRead(cBlkFile [b], contribSize, contrib) + ) { + /* Corresponding photon/contrib should be in subblock */ + perror("contribPhotonSort: unexpected EOF in block file"); + return -1; + } + + if (OOC_SortWrite(photonOut, PHOTONSIZE, (char*)&photon) || + OOC_SortWrite(contribOut, contribSize, contrib) + ) { + perror("contribPhotonSort: error writing output file"); + return -1; + } + + #ifdef DEBUG_OOC_SORT + pqCnt++; + #endif + /* Peek next photon/contrib from same subblock and enqueue + for next iteration */ + if (!OOC_SortPeek(pBlkFile [b], PHOTONSIZE, (char*)&photon) && + !OOC_SortPeek(cBlkFile [b], contribSize, contrib) + ) { + /* Subblock not exhausted */ + pri = PMAPCONTRIB_KEY(&photon, keyData); + + if (OOC_SortPush(&pqueue, pri, b) < 0) { + perror("contribPhotonSort: failed priority enqueue"); + return -1; + } + } + } + } while (b >= 0); + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "*DBG* contribPhotonSort: dequeued %lu rec\n", + pqCnt + ); + fprintf(stderr, "*DBG* contribPhotonSort: merged file contains " + "%lu photons, %lu contribs\n", + ftell(photonOut) / PHOTONSIZE, ftell(contribOut) / contribSize + ); + #endif + + /* Priority queue now empty -> done; close temporary subblock files, + * which deletes them (see note in kid kode above) */ + for (b = 0; b < numBlk; b++) { + fclose(pBlkFile [b]); + fclose(cBlkFile [b]); + } + + /* !!! Commit output file to disk before caller reads it; omitting + * !!! this apparently leads to corrupted files (race condition?) + * !!! when the caller tries to read them! */ + fflush(photonOut); + fflush(contribOut); + fsync(fileno(photonOut)); + fsync(fileno(contribOut)); + } + + else { + /* ========================================================== + * Block is small enough for in-core sort + * ========================================================== */ + + const unsigned long pBlkSize = blkLen * PHOTONSIZE, + cBlkSize = blkLen * contribSize; + int pIfd = fileno(photonIn), + pOfd = fileno(photonOut), + cIfd = fileno(contribIn), + cOfd = fileno(contribOut); + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "*DBG* contribPhotonSort: Proc %d (%d/%d) " + "sorting block [%lu - %lu]\n", + getpid(), procCnt, numProc - 1, blkLo, blkHi + ); + #endif + + /* Atomically seek and read block into in-core sort buffer */ + /* !!! Unbuffered I/O via pread() avoids potential race conditions + * !!! and buffer corruption which can occur with lseek()/fseek() + * !!! followed by read()/fread(). */ + if (pread(pIfd, photonBuf, pBlkSize, blkLo*PHOTONSIZE) != pBlkSize || + pread(cIfd, contribBuf, cBlkSize, blkLo*contribSize) != cBlkSize + ) { + perror("contribPhotonSort: error reading block for sort"); + return -1; + } + + /* Quicksort block in-core and write to output file */ + contribPhotonQuickSort(photonBuf, contribBuf, contribSize, + keyData, blkLo, blkHi + ); + + if (write(pOfd, photonBuf, pBlkSize) != pBlkSize || + write(cOfd, contribBuf, cBlkSize) != cBlkSize + ) { + perror("contribPhotonSort: error writing sorted block"); + return -1; + } + + fsync(pOfd); + fsync(cOfd); + + #ifdef DEBUG_OOC_SORT + fprintf(stderr, "*DBG* contribPhotonSort: proc %d wrote " + "%lu photons, %lu contribs\n", getpid(), + lseek(pOfd, 0, SEEK_END) / PHOTONSIZE, + lseek(cOfd, 0, SEEK_END) / contribSize + ); + #endif + } + + return 0; + } + + + + int contribPhotonSort ( + FILE *photonIn, FILE *contribIn, FILE *photonOut, FILE *contribOut, + unsigned long numPhotons, unsigned contribSize, + FVECT bbOrg, RREAL bbSize, unsigned numBlk, unsigned long blkSize, + unsigned numProc, RREAL *(*keyFunc)(const void*) + ) + /* Sort photons and their corresponding binned contributions by + * subdividing input files into small blocks, sorting these in-core, and + * merging them out-of-core via a priority queue. This code is adapted + * from oocsort.{h,c} + * + * The sort is implemented as a recursive (numBlk)-way sort; the input + * files are successively split into numBlk smaller blocks until these + * are of size <= blkSize, i.e. small enough for in-core sorting, then + * merging the sorted blocks as the stack unwinds. The in-core sort is + * parallelised over numProc processes. + * + * Parameters are as follows: + * photonIn Opened primary input file containing unsorted photons; + * their Morton codes determine the ordering + * contribIn Opened secondary input file containing unsorted binned + * contributions size contribSize + * photonOut Opened primary output file for sorted photons + * contribOut Opened secondary output file for sorted binned + * contributions ordered by corresponding photons + * numPhotons Number of photon/contrib records in input files + * contribSize Size of binned contributions per photon, in bytes + * numBlk Number of blocks to divide into / merge from + * blkSize Max block size and size of in-core sort buffer, in bytes + * numProc Number of parallel processes for in-core sort + * bbOrg Origin of bounding box containing photons + * bbSize Extent of bounding box containing photons + * keyFunc Callback to access photon coords to generate Morton code + */ + { + Photon *photonBuf = NULL; + char *contribBuf = NULL; + unsigned long bufLen; + int stat; + OOC_KeyData keyData; + + if (numBlk < 1) + numBlk = 1; + + if (!numPhotons || !photonIn || !contribIn || !contribSize || + !numBlk || !blkSize + ) + error(INTERNAL, "contribPhotonSort: nothing to sort"); + + /* Figure out maximum buffer length (number of photons/contribs we + can fit into the maximum block size) and allocate in-core sort + buffers accordingly. Quit if block size too small */ + if (!(bufLen = blkSize / max(PHOTONSIZE, contribSize))) + error(INTERNAL, "contribPhotonSort: block size too small"); + if (!(photonBuf = calloc(bufLen, PHOTONSIZE)) || + !(contribBuf = calloc(bufLen, contribSize)) + ) + error(SYSTEM, + "contribPhotonSort: cannot allocate in-core sort buffer" + ); + + /* Set key callback and data; reduces params on stack in recursion */ + keyData.key = keyFunc; + keyData.mortonScale = MORTON3D_MAX / bbSize; + VCOPY(keyData.bbOrg, bbOrg); + + stat = contribPhotonSortRecurse( + photonIn, contribIn, photonOut, contribOut, + contribSize, numBlk, numProc, photonBuf, contribBuf, bufLen, + &keyData, 0, numPhotons - 1 + ); + + /* Cleanup */ + free(photonBuf); + free(contribBuf); + + return stat; + } + +#endif /* NIX / PMAP_OOC */ + diff --git a/pmcontribsort.h b/pmcontribsort.h new file mode 100644 index 0000000..1e4b5b3 --- /dev/null +++ b/pmcontribsort.h @@ -0,0 +1,56 @@ +/* + ====================================================================== + Header for N-way out-of-core merge sort of precomputed photon + contributions. + + Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation + (SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") + ====================================================================== + + $Id$ +*/ + + +#ifndef _PMAPCONTRIB_SORT_H + #define _PMAPCONTRIB_SORT_H + + #include "fvect.h" + + + int contribPhotonSort ( + FILE *photonIn, FILE *contribIn, FILE *photonOut, FILE *contribOut, + unsigned long numPhotons, unsigned contribSize, + FVECT bbOrg, RREAL bbSize, unsigned numBlk, unsigned long blkSize, + unsigned numProc, RREAL *(*keyFunc)(const void*) + ); + /* Sort photons and their corresponding binned contributions by + * subdividing input files into small blocks, sorting these in-core, and + * merging them out-of-core via a priority queue. This code is adapted + * from oocsort.{h,c} + * + * The sort is implemented as a recursive (numBlk)-way sort; the input + * files are successively split into numBlk smaller blocks until these + * are of size <= blkSize, i.e. small enough for in-core sorting, then + * merging the sorted blocks as the stack unwinds. The in-core sort is + * parallelised over numProc processes. + * + * Parameters are as follows: + * photonIn Opened primary input file containing unsorted photons; + * their Morton codes determine the ordering + * contribIn Opened secondary input file containing unsorted binned + * contributions size contribSize + * photonOut Opened primary output file for sorted photons + * contribOut Opened secondary output file for sorted binned + * contributions ordered by corresponding photons + * numPhotons Number of photon/contrib records in input files + * contribSize Size of binned contributions per photon, in bytes + * numBlk Number of blocks to divide into / merge from + * blkSize Max block size and size of in-core sort buffer, in bytes + * numProc Number of parallel processes for in-core sort + * bbOrg Origin of bounding box containing photons + * bbSize Extent of bounding box containing photons + * keyFunc Callback to access photon coords to generate Morton code + */ +#endif