diff --git a/Rmakefile b/Rmakefile index 91a72c2..bea526f 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,469 +1,469 @@ # RCSid: $Id: Rmakefile,v 2.87 2021/01/19 23:31:47 greg Exp $ # # Compiles for ray tracing programs. # OPT = -O MACH = -DBSD CFLAGS = -I../common -L../lib $(OPT) $(MACH) SPECIAL = CC = cc AR = ar MLIB = -lm LINT = lint LINTFLAGS = -DBSD # # The following are user-definable: # DESTDIR = . INSTDIR = /usr/local/bin INSTALL = cp # # The following paths must exist and be relative to root: # DEVDIR = $(INSTDIR)/dev LIBDIR = /usr/local/lib/ray # # Library routines: # RLIB = ../lib/libradiance.a RCLIB = ../lib/libraycalls.a LIBS = -lrtrad $(MLIB) # # Device drivers for rvu (see also devtable.c): # DOBJS = devtable.o devcomm.o editline.o x11.o x11twind.o \ colortab.o DSRC = devtable.c devcomm.c editline.c x11.c x11twind.c \ colortab.c DLIBS = -lX11 # # Standard object files: # RTOBJS = rtmain.o rtrace.o duphead.o persist.o RTSRC = rtmain.c rtrace.c duphead.c persist.c RPOBJS = rpmain.o rpict.o srcdraw.o duphead.o persist.o RPSRC = rpmain.c rpict.c srcdraw.c duphead.c persist.c RVOBJS = rvmain.o rview.o rv2.o rv3.o $(DOBJS) RVSRC = rvmain.c rview.c rv2.c rv3.c $(DSRC) RCOBJS = rcmain.o rcontrib.o rc2.o rc3.o RCSRC = rcmain.c rcontrib.c rc2.c rc3.c RLOBJS = raycalls.o raypcalls.o rayfifo.o RLSRC = raycalls.c raypcalls.c rayfifo.c ROBJS = $(RAYOBJS) $(SURFOBJS) $(MATOBJS) \ $(MODOBJS) $(SUPPOBJS) $(PMOBJS) RSRC = $(RAYSRC) $(SURFSRC) $(MATSRC) \ $(MODSRC) $(SUPPSRC) RAYOBJS = ambcomp.o ambient.o ambio.o freeobjmem.o initotypes.o \ preload.o raytrace.o renderopts.o RAYSRC = ambcomp.c ambient.c ambio.c freeobjmem.c initotypes.c \ preload.c raytrace.c renderopts.c SURFOBJS = source.o sphere.o srcobstr.o srcsupp.o srcsamp.o virtuals.o \ o_face.o o_cone.o o_instance.o o_mesh.o SURFSRC = sphere.c source.c srcobstr.c srcsupp.c srcsamp.c virtuals.c \ o_face.c srcsamp.c o_cone.c o_instance.c o_mesh.c MATOBJS = aniso.o normal.o dielectric.o m_clip.o glass.o m_brdf.o \ m_mirror.o m_direct.o m_mist.o fprism.o m_alias.o m_bsdf.o \ ashikhmin.o MATSRC = aniso.c normal.c dielectric.c m_clip.c glass.c m_brdf.c \ m_mirror.c m_direct.c m_mist.c fprism.c m_alias.c m_bsdf.c \ ashikhmin.c MODOBJS = p_func.o t_func.o p_data.o t_data.o text.o mx_func.o mx_data.o MODSRC = p_func.c t_func.c p_data.c t_data.c text.c mx_func.c mx_data.c SUPPOBJS = func.o noise3.o data.o SUPPSRC = func.c noise3.c data.c PMOBJS = pmap.o pmapsrc.o pmapmat.o pmaprand.o pmapio.o pmapdata.o pmapdens.o \ pmapbias.o pmapparm.o pmapcontrib.o pmcontrib2.o pmapamb.o pmapray.o \ - pmapopt.o pmapdiag.o pmaptype.o oocmorton.o oococt.o oocsort.o \ + pmapopt.o pmapdiag.o pmaptype.o morton.o oococt.o oocsort.o \ oocbuild.o oocnn.o ooccache.o pmutil.o pmapfilt.o PMSRC = pmap.c pmapsrc.c pmapmat.c pmaprand.c pmapio.c pmapdata.c pmapdens.c \ pmapbias.c pmapparm.c pmapcontrib.c pmcontrib2.c pmapamb.c pmapray.c \ - pmapopt.c pmapdiag.c pmaptype.c oocmorton.c oococt.c oocsort.c \ + pmapopt.c pmapdiag.c pmaptype.c morton.c oococt.c oocsort.c \ oocbuild.c oocnn.c ooccache.c pmutil.c pmapfilt.c pmapkdt.c pmapooc.c HEADERS = ambient.h ray.h data.h otspecial.h source.h # # What this makefile produces: # PROGS = $(DESTDIR)/rtrace $(DESTDIR)/rpict $(DESTDIR)/rvu $(DESTDIR)/rcontrib \ $(DESTDIR)/lookamb $(DESTDIR)/mkpmap $(DESTDIR)/pmapdump all: $(PROGS) $(RCLIB) $(SPECIAL) install: all rayinit.cal $(INSTALL) $(PROGS) $(INSTDIR) rm -f $(LIBDIR)/rayinit.cal cp rayinit.cal $(LIBDIR) ogl: clean: set nonomatch; rm -f $(PROGS) *.o lint: $(RVSRC) $(LINT) $(LINTFLAGS) -DRVIEW $(RVSRC) $(LIBS) # # Links: # $(DESTDIR)/rtrace: $(RTOBJS) $(RCLIB) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rtrace $(RTOBJS) $(RCLIB) $(RLIB) $(LIBS) $(DESTDIR)/rpict: $(RPOBJS) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rpict $(RPOBJS) $(RLIB) $(LIBS) $(DESTDIR)/rvu: $(RVOBJS) $(RCLIB) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rvu $(RVOBJS) $(RCLIB) \ $(RLIB) $(LIBS) $(DLIBS) $(DESTDIR)/rcontrib: $(RCOBJS) $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/rcontrib $(RCOBJS) $(RLIB) $(LIBS) $(DESTDIR)/lookamb: lookamb.o ambio.o $(CC) $(CFLAGS) -o $(DESTDIR)/lookamb lookamb.o ambio.o $(LIBS) $(DESTDIR)/mkpmap: mkpmap.o $(RLIB) $(CC) $(CFLAGS) -o $(DESTDIR)/mkpmap mkpmap.o $(RLIB) $(LIBS) $(DESTDIR)/pmapdump: pmapdump.o pmaptype.o pmapparm.o $(CC) $(CFLAGS) -o pmapdump pmapdump.o pmaptype.o pmapparm.o $(LIBS) $(RLIB): $(ROBJS) Version.o rm -f $(RLIB) $(AR) rc $(RLIB) $(ROBJS) Version.o -ranlib $(RLIB) $(RCLIB): $(RLOBJS) rm -f $(RCLIB) $(AR) rc $(RCLIB) $(RLOBJS) -ranlib $(RCLIB) # # Uncomment the following to model dispersion: # dielectric.o: dielectric.c source.h $(CC) $(CFLAGS) -DDISPERSE -c dielectric.c # end of dispersion compiles. devcomm.o: devcomm.c $(CC) $(CFLAGS) -DDEVPATH=\"$(DEVDIR)\" -c devcomm.c # # Version module: # Version.c: VERSION $(RSRC) $(HEADERS) ( cat VERSION ; date ; whoami ; hostname ) > Version.c ed - Version.c < verscript.ed || rm Version.c # # Include dependencies: # ambio.o colortab.o data.o devcomm.o \ devmain.o lookamb.o rview.o x11.o: ../common/color.h freeobjmem.o o_cone.o srcsupp.o: ../common/cone.h data.o freeobjmem.o m_brdf.o mx_data.o \ p_data.o raycalls.o t_data.o: data.h devcomm.o devmain.o devtable.o \ editline.o x11.o: driver.h freeobjmem.o o_face.o srcsupp.o: ../common/face.h ambient.o raytrace.o rpmain.o rtmain.o \ rtrace.o rvmain.o rv2.o rv3.o: ../common/octree.h o_instance.o: ../common/instance.h ambient.o aniso.o ashikhmin.o dielectric.o freeobjmem.o func.o glass.o \ initotypes.o m_brdf.o m_bsdf.o m_direct.o m_mirror.o normal.o o_cone.o \ preload.o raycalls.o raytrace.o rtrace.o rv2.o source.o sphere.o \ srcsupp.o text.o srcdraw.o srcobstr.o virtuals.o: ../common/otypes.h ambient.o ambcomp.o aniso.o ashikhmin.o normal.o raycalls.o raytrace.o \ rpict.o rvmain.o rtmain.o rpmain.o rcmain.o persist.o source.o rv3.o \ srcsamp.o virtuals.o: ../common/random.h ambcomp.o ambient.o aniso.o ashikhmin.o dielectric.o freeobjmem.o func.o \ glass.o m_bsdf.o m_brdf.o m_clip.o m_direct.o m_mirror.o m_mist.o mx_data.o \ o_mesh.o mx_func.o normal.o o_cone.o o_face.o o_instance.o p_data.o p_func.o \ raycalls.o raypcalls.o rayfifo.o raytrace.o rpict.o rtrace.o rv2.o rv3.o rview.o \ source.o sphere.o srcdraw.o srcobstr.o srcsamp.o srcsupp.o t_data.o t_func.o \ text.o rpmain.o rtmain.o rvmain.o virtuals.o m_alias.o rcmain.o \ rcontrib.o rc2.o rc3.o: ray.h \ ../common/standard.h ../common/rtmisc.h ../common/rtio.h ../common/rtmath.h \ ../common/rterror.h ../common/octree.h \ ../common/mat4.h ../common/fvect.h ../common/object.h ../common/color.h rv2.o rv3.o rview.o: rpaint.h driver.h ../common/view.h ../common/resolu.h m_direct.o m_mirror.o m_mist.o dielectric.o raycalls.o \ rpict.o rpmain.o rtmain.o rvmain.o source.o srcdraw.o \ srcobstr.o srcsamp.o srcsupp.o virtuals.o: source.h cone.o data.o devcomm.o initotypes.o fprism.o preload.o \ duphead.o octree.o: ../common/standard.h ../common/rtmisc.h \ ../common/rtio.h ../common/rtmath.h \ ../common/rterror.h ../common/mat4.h ../common/fvect.h ambio.o: ../common/fvect.h ambient.o initotypes.o m_alias.o pmapcontrib.o pmapdata.o pmapsrc.o \ pmcontrib2.o raytrace.o rtrace.o rv2.o rv3.o source.o \ srcdraw.o srcobstr.o virtuals.o: otspecial.h rpmain.o rtmain.o rvmain.o rpict.o \ srcdraw.o: ../common/view.h ../common/resolu.h rpict.o: ../common/hilbert.h x11.o x11twind.o: x11twind.h x11.o: x11icon.h ambient.o ambcomp.o ambio.o lookamb.o raycalls.o: ambient.h data.o rpmain.o rtmain.o rvmain.o rpict.o rtrace.o \ rv2.o: ../common/resolu.h aniso.o func.o m_brdf.o m_direct.o mx_data.o mx_func.o p_data.o \ p_func.o t_data.o t_func.o: func.h ../common/calcomp.h preload.o: data.h func.h ../common/object.h ../common/face.h \ ../common/cone.h ../common/instance.h ../common/mesh.h \ ../common/color.h ../common/bsdf.h ../common/otypes.h rtmain.o rpmain.o rvmain.o persist.o duphead.o \ renderopts.o rpict.o: ../common/paths.h freeobjmem.o raycalls.o text.o: ../common/font.h raypcalls.o: ../common/selcall.h o_mesh.o: ../common/mesh.h noise3.o: ../common/calcomp.h aniso.o ashikhmin.o dielectric.o freeobjmem.o glass.o initotypes.o \ m_alias.o m_brdf.o m_clip.o m_direct.o m_mirror.o m_mist.o \ mx_data.o mx_func.o normal.o o_cone.o o_face.o o_instance.o \ o_mesh.o p_data.o p_func.o source.o sphere.o t_data.o t_func.o \ srcobstr.o text.o: rtotypes.h m_bsdf.o: ambient.h source.h func.h \ ../common/calcomp.h ../common/bsdf.h ../common/random.h rcmain.o rcontrib.o rc2.o rc3.o: rcontrib.h \ ../common/platform.h ../common/paths.h ../common/lookup.h \ func.h ../common/calcomp.h ../common/rtprocess.h ambient.o rcmain.o: ambient.h rcmain.o: source.h rcontrib.o: source.h ../common/otypes.h rc2.o: ../common/resolu.h rc3.o: ../common/selcall.h # # Photon map include dependencies (via 'gcc -MM -I../common') # 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 oocmorton.h +oocsort.o: oocsort.c oocsort.h ../common/fvect.h morton.h -oocbuild.o: oocbuild.c oococt.h oocmorton.h ../common/fvect.h ooccache.h \ +oocbuild.o: oocbuild.c oococt.h morton.h ../common/fvect.h ooccache.h \ oocsort.h -oococt.o: oococt.c oococt.h oocmorton.h ../common/fvect.h ooccache.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 oocmorton.h ../common/fvect.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 pmapooc.h pmapooc.c mkpmap.o: pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h ../common/paths.h \ ../common/lookup.h pmapmat.h pmapsrc.h source.h pmapcontrib.h \ pmaprand.h ambient.h ../common/resolu.h \ pmapkdt.h pmapkdt.c 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 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 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 pmapooc.h pmapooc.c pmapamb.o: pmapamb.c pmapamb.h pmapdata.h ray.h ../common/standard.h \ ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ ../common/rterror.h ../common/octree.h ../common/object.h \ ../common/color.h pmapparm.h pmaptype.h ../common/paths.h \ ../common/lookup.h pmap.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c pmapmat.o: pmapmat.c pmapmat.h pmap.h pmapparm.h pmaptype.h pmapdata.h \ ray.h ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h ../common/paths.h \ ../common/lookup.h pmaprand.h ../common/otypes.h data.h func.h \ ../common/calcomp.h ../common/bsdf.h ../common/ccolor.h \ ../common/platform.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c pmapsrc.o: pmapsrc.c pmapsrc.h ray.h ../common/standard.h \ ../common/copyright.h ../common/rtio.h ../common/rtmisc.h \ ../common/rtmath.h ../common/tiff.h ../common/mat4.h ../common/fvect.h \ ../common/rterror.h ../common/octree.h ../common/object.h \ ../common/color.h pmapparm.h pmaptype.h source.h pmap.h pmapdata.h \ ../common/paths.h ../common/lookup.h pmaprand.h \ ../common/otypes.h otspecial.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c pmapopt.o: pmapopt.c ray.h ../common/standard.h ../common/copyright.h \ ../common/rtio.h ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h \ ../common/mat4.h ../common/fvect.h ../common/rterror.h \ ../common/octree.h ../common/object.h ../common/color.h pmapparm.h \ pmaptype.h pmapparm.o: pmapparm.c pmapparm.h pmaptype.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h ../common/paths.h ../common/lookup.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c pmutil.o: pmutil.c pmap.h pmapparm.h pmaptype.h pmapdata.h ray.h \ ../common/standard.h ../common/copyright.h ../common/rtio.h \ ../common/rtmisc.h ../common/rtmath.h ../common/tiff.h ../common/mat4.h \ ../common/fvect.h ../common/rterror.h ../common/octree.h \ ../common/object.h ../common/color.h ../common/paths.h \ ../common/lookup.h pmapio.h \ pmapkdt.h pmapkdt.c pmapooc.h pmapooc.c diff --git a/morton.c b/morton.c new file mode 100644 index 0000000..9ee5fbc --- /dev/null +++ b/morton.c @@ -0,0 +1,67 @@ +#ifndef lint +static const char RCSid[] = "$Id$"; +#endif + + +/* + ========================================================================= + Routines to generate and compare Morton Codes, i.e. indices on space + filling Z-curves. + + 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$ +*/ + + +#if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) +/* No Windoze support for now */ + +#include "morton.h" + + + +MortonIdx Key2Morton3D (const FVECT key, const FVECT org, RREAL scale) +/* Compute Morton code (Z-curve index) of length 3 * MORTON3D_BITS bits + for 3D key within bounding box defined by org and scaled to maximum + index with scale */ +{ + unsigned i; + MortonIdx k [3]; + + /* Normalise key and map each dim to int of MORTON3D_BITS */ + for (i = 0; i < 3; i++) + k [i] = scale * (key [i] - org [i]); + + /* Interleave each dim with zeros and merge */ + return (Morton3DBitInterleave(k [0]) | + Morton3DBitInterleave(k [1]) << 1 | + Morton3DBitInterleave(k [2]) << 2 + ); +} + + + +MortonIdx Key2Morton2D (const RREAL key [2], const RREAL org [2], RREAL scale) +/* Compute Morton code (Z-curve index) of length 2 * MORTON2D_BITS bits + for 2D key within bounding box defined by org and scaled to maximum + index with scale */ +{ + unsigned i; + MortonIdx k [2]; + + /* Normalise key and map each dim to int of MORTON2D_BITS */ + for (i = 0; i < 2; i++) + k [i] = scale * (key [i] - org [i]); + + /* Interleave each dim with zeros and merge */ + return (Morton2DBitInterleave(k [0]) | Morton2DBitInterleave(k [1]) << 1); +} + +#endif /* NIX / PMAP_OOC */ + diff --git a/morton.h b/morton.h new file mode 100644 index 0000000..06e6be4 --- /dev/null +++ b/morton.h @@ -0,0 +1,74 @@ +/* + ====================================================================== + Routines to generate and compare Morton Codes, i.e. indices on space + filling Z-curves. + + 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$ +*/ + +#ifndef _MORTON_H + #define _MORTON_H + + #include "fvect.h" + #include + + + /* Type to represent Morton codes (Z-curve indices) */ + typedef uint64_t MortonIdx; + + + /* Number of bits per dimension for Morton code (Z-curve index) in 2D + and 3D. This corresponds to the maximum number of bounding box + subdivisions which can be resolved per axis; further subdivisions + will map to the same Morton code. */ + #define MORTON3D_BITS 21 + #define MORTON3D_MAX (((MortonIdx)1 << MORTON3D_BITS) - 1) + + #define MORTON2D_BITS 32 + #define MORTON2D_MAX (((MortonIdx)1 << MORTON2D_BITS) - 1) + + + /* Interleave lower MORTONxx_BITS bits of k with 00 for 3D (resp. 0 for + 2D), resulting in 3*MORTON3D_BITS (resp. 2*MORTON2D_BITS) bits. This + is an optimised "bitmask hack" version code taken from: + http://www.forceflow.be/2013/10/07/ + morton-encodingdecoding-through-bit-interleaving-implementations/ + and (for 32-bit codes): + https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/ + */ + #define Morton3DBitInterleave(k) ((k) &= MORTON3D_MAX, \ + (k) = ((k) | (k) << 32) & 0x001f00000000ffff, \ + (k) = ((k) | (k) << 16) & 0x001f0000ff0000ff, \ + (k) = ((k) | (k) << 8) & 0x100f00f00f00f00f, \ + (k) = ((k) | (k) << 4) & 0x10c30c30c30c30c3, \ + (k) = ((k) | (k) << 2) & 0x1249249249249249 \ + ) + + #define Morton2DBitInterleave(k) ((k) &= MORTON2D_MAX, \ + (k) = ((k) | (k) << 16) & 0x0000ffff0000ffff, \ + (k) = ((k) | (k) << 8) & 0x00ff00ff00ff00ff, \ + (k) = ((k) | (k) << 4) & 0x0f0f0f0f0f0f0f0f, \ + (k) = ((k) | (k) << 2) & 0x3333333333333333, \ + (k) = ((k) | (k) << 1) & 0x5555555555555555 \ + ) + + + /* Compute Morton code (Z-curve index) of length 3 * MORTON3D_BITS bits + for 3D key within bounding box defined by org and scaled to maximum + index with scale */ + MortonIdx Key2Morton3D (const FVECT key, const FVECT org, RREAL scale); + + /* As above, but compute Morton code of length 2 * MORTON2D_BITS bits + for 2D key */ + MortonIdx Key2Morton2D (const RREAL key [2], const RREAL org [2], + RREAL scale + ); +#endif /* _MORTON_H */ + diff --git a/oocbuild.c b/oocbuild.c index 04530b5..66bf6fb 100644 --- a/oocbuild.c +++ b/oocbuild.c @@ -1,359 +1,367 @@ #ifndef lint static const char RCSid[] = "$Id: oocbuild.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $"; #endif /* ======================================================================= Routines for building out-of-core octree data structure Adapted from: Kontkanen J., Tabellion E. and Overbeck R.S., "Coherent Out-of-Core Point-Based Global Illumination", EGSR 2011. 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: oocbuild.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $ */ #if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) /* No Windoze support for now */ #include "oococt.h" #include "oocsort.h" #include #include /* Test for empty/full input queue, return pointer to head/tail */ #define QueueFull(q) ((q) -> len == (q) -> cap) #define QueueEmpty(q) (!(q) -> len) #define QueueHead(q) ((q) -> data + (q) -> head * (q) -> recSize) #define QueueTail(q) ((q) -> data + \ ((q)->head + (q)->len-1) % (q)->cap * (q)->recSize) /* Input queue for bottom-up octree construction */ typedef struct { void *data; unsigned head, len, cap, recSize; /* Queue head, length (from head), capacity and record size */ FILE *in; /* Input file for data records */ } OOC_BuildQueue; static OOC_BuildQueue *QueueInit (OOC_BuildQueue *q, unsigned recSize, unsigned capacity) /* Initialise queue of #capacity records of size recSize each; returns queue * pointer or NULL if failed. */ { if (!(q && (q -> data = calloc(capacity, recSize)))) return NULL; q -> cap = capacity; q -> recSize = recSize; q -> head = q -> len = 0; return q; } static int QueuePush (OOC_BuildQueue *q, const void *rec) /* Append record to queue tail; return new queue length or -1 on failure */ { if (!q || !rec || QueueFull(q)) return -1; ++q->len; memcpy(QueueTail(q), rec, q -> recSize); return q -> len; } static int QueuePop (OOC_BuildQueue *q, void *rec) /* Remove record from queue head and return in rec if not NULL; return new * queue length or -1 on failure */ { if (!q || QueueEmpty(q)) return -1; /* Return head if rec != NULL */ if (rec) memcpy(rec, QueueHead(q), q -> recSize); q -> head = (q -> head + 1) % q -> cap; return --q -> len; } static int QueueFill (OOC_BuildQueue *q) /* Read records from q -> in until the queue is full; return queue * length or -1 on failure */ { static void *rec = NULL; if (!rec && !(rec = malloc(q -> recSize))) return -1; while (!QueueFull(q) && !feof(q -> in) && fread(rec, q -> recSize, 1, q -> in)) QueuePush(q, rec); return q -> len; } static OOC_DataIdx OOC_BuildRecurse (OOC_Octree *oct, OOC_Node* node, FVECT org, RREAL size, unsigned depth, OOC_BuildQueue *queue) /* Recursive part of OOC_Build(); insert records from input queue into * octree node and subdivide into kids if necessary. Returns number of * records in subtree or OOC_DATAIDX_ERR if failed. */ { int k; const RREAL kidSize = size * 0.5; if (!oct || !node) return OOC_DATAIDX_ERR; if (QueueEmpty(queue) || !OOC_InBBox(oct, org, size, oct -> key(QueueHead(queue)))) /* Input exhausted or queue head outside node */ return 0; if (QueueFull(queue) && depth < oct -> maxDepth && OOC_InBBox(oct, org, size, oct -> key(QueueTail(queue)))) { /*************************** SUBDIVIDE NODE ************************* * At least leafMax + 1 records (since the queue is full) lie inside * the current node's bounding box, and maxDepth has not been reached * ==> subdivide this node. * (Note it suffices to only test the queue tail against the bounding * box, as the records are in Z-order) ********************************************************************/ OOC_Node kid [8]; OOC_DataIdx dataCnt; FVECT kidOrg; #ifdef DEBUG_OOC_BUILD -FVECT key; -unsigned k2 = 0; + FVECT key; + unsigned k2 = 0; #endif /* We recurse on the nonempty kids first, then finalise their nodes so * they are ordered consecutively, since the parent only indexes the 1st kid */ for (k = 0; k < 8; k++) { /* Clear kid node and get its octant origin */ OOC_CLEARNODE(kid + k); VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, k, kidSize); /* Recurse on kid and check for early bailout */ if (OOC_BuildRecurse(oct, kid + k, kidOrg, kidSize, depth + 1, queue) == OOC_DATAIDX_ERR) return OOC_DATAIDX_ERR; #ifdef DEBUG_OOC_BUILD -if (!QueueEmpty(queue)) { -VCOPY(key, oct -> key(QueueHead(queue))); -k2 = OOC_Branch(oct, org, kidSize, key, NULL); -if (OOC_InBBox(oct, org, size, key) && k2 < k) { - fprintf(stderr, - "OOC_BuildRecurse, node subdiv: unsorted key [%f, %f, %f] with " - "octant %d (last %d with bbox [%f-%f, %f-%f, %f-%f])\n", - key [0], key [1], key [2], k2, k, - kidOrg [0], kidOrg [0] + kidSize, kidOrg [1], kidOrg [1] + kidSize, - kidOrg [2], kidOrg [2] + kidSize); -}} + if (!QueueEmpty(queue)) { + VCOPY(key, oct -> key(QueueHead(queue))); + k2 = OOC_Branch(oct, org, kidSize, key, NULL); + if (OOC_InBBox(oct, org, size, key) && k2 < k) { + fprintf(stderr, + "OOC_BuildRecurse, node subdiv: unsorted key [%f, %f, %f] " + "with octant %d (last %d with bbox [%f-%f, %f-%f, %f-%f])\n", + key [0], key [1], key [2], k2, k, + kidOrg [0], kidOrg [0] + kidSize, + kidOrg [1], kidOrg [1] + kidSize, + kidOrg [2], kidOrg [2] + kidSize + ); + } + } #endif } /* Now finalise consecutive kid nodes, skipping empty ones */ for (k = 0; k < 8; k++) if ((dataCnt = OOC_DATACNT(kid + k))) { /* Nonzero kid ==> allocate and set node */ if (!NewNode(oct)) { fputs("OOC_BuildRecurse: failed to allocate new node\n", stderr); return OOC_DATAIDX_ERR; } OOC_SETROOT(oct, kid + k); /* Sum kid's data count to parent's and check for overflow */ if ((dataCnt += node -> node.num) <= OOC_DATAIDX_MAX) node -> node.num = dataCnt; else { fputs("OOC_BuildRecurse: data count overflow in node\n", stderr); return OOC_DATAIDX_ERR; } /* Set kid index in parent (if first kid) and corresponding * branch bit. The kid is the most recent node and thus at the * end of the node array, which coincides with the current * subtree root */ if (!node -> node.branch) node -> node.kid = OOC_ROOTIDX(oct); OOC_SETBRANCH(node, k); } } else { /****************************** MAKE LEAF **************************** * Queue contains no more than leafMax records, queue tail lies * outside node's bounding box, or maxDepth reached * ==> make this node a leaf. *********************************************************************/ RREAL *key; #ifdef DEBUG_OOC_BUILD -OOC_MortonIdx zIdx, lastzIdx = 0; -FVECT /* key, */ - lastKey, kidOrg; -unsigned lastk = 0; + MortonIdx zIdx, lastzIdx = 0; + FVECT /* key, */ + lastKey, kidOrg; + unsigned lastk = 0; #endif /* Mark as leaf (note it's been cleared by the parent call) */ OOC_SETLEAF(node); while (!QueueEmpty(queue) && OOC_InBBox(oct, org, size, (key = oct->key(QueueHead(queue))))) { /* Record lies inside leaf; increment data counter for octant * containing record. */ if ((k = OOC_Branch(oct, org, kidSize, key, NULL)) < 0) { /* Shouldn't happen, as key tested within bbox above */ fprintf(stderr, "OOC_BuildRecurse: buggered Morton code, " "disruption in space-time continuum?\n"); return OOC_DATAIDX_ERR; } if (node -> leaf.num [k] == OOC_OCTCNT_MAX) { /* Currently we're buggered here; merge records instead? */ fprintf(stderr, "OOC_BuildRecurse: data count overflow in " "leaf: depth = %d, count = %d\n", depth, node -> leaf.num [k]); return OOC_DATAIDX_ERR; } ++node -> leaf.num [k]; - + #ifdef DEBUG_OOC_BUILD -/* VCOPY(key, oct -> key(QueueHead(queue))); */ -if ((zIdx = OOC_KEY2MORTON(key, oct)) < lastzIdx) { - fprintf(stderr, "OOC_BuildRecurse, make leaf: unsorted zIdx %020ld for " - "key [%f, %f, %f] (previous zIdx %020ld for " - "key [%f, %f, %f]\n", zIdx, key [0], key [1], key [2], - lastzIdx, lastKey [0], lastKey [1], lastKey [2]); -} + /* VCOPY(key, oct -> key(QueueHead(queue))); */ + if ((zIdx = OOC_KEY2MORTON(key, oct)) < lastzIdx) { + fprintf(stderr, + "OOC_BuildRecurse, make leaf: unsorted zIdx %020ld for " + "key [%f, %f, %f] (previous zIdx %020ld for " + "key [%f, %f, %f]\n", zIdx, key [0], key [1], key [2], + lastzIdx, lastKey [0], lastKey [1], lastKey [2] + ); + } -VCOPY(kidOrg, org); -OOC_OCTORIGIN(kidOrg, k, kidSize); -if (k < lastk || zIdx < lastzIdx) { - fprintf(stderr, - "OOC_BuildRecurse, make leaf: unsorted octant %d (last %d) with " - "bbox [%f-%f, %f-%f, %f-%f] for key [%f, %f, %f] with zIdx %020ld " - "(last [%f, %f, %f], zIdx %020ld)\n", - k, lastk, kidOrg [0], kidOrg [0] + kidSize, - kidOrg [1], kidOrg [1] + kidSize, kidOrg [2], kidOrg [2] + kidSize, - key [0], key [1], key [2], zIdx, - lastKey [0], lastKey [1], lastKey [2], lastzIdx); -} -lastk = k; -lastzIdx = zIdx; -VCOPY(lastKey, key); + VCOPY(kidOrg, org); + OOC_OCTORIGIN(kidOrg, k, kidSize); + if (k < lastk || zIdx < lastzIdx) { + fprintf(stderr, + "OOC_BuildRecurse, make leaf: unsorted octant %d (last %d) " + "with bbox [%f-%f, %f-%f, %f-%f] for key [%f, %f, %f] with " + "zIdx %020ld (last [%f, %f, %f], zIdx %020ld)\n", + k, lastk, kidOrg [0], kidOrg [0] + kidSize, + kidOrg [1], kidOrg [1] + kidSize, + kidOrg [2], kidOrg [2] + kidSize, + key [0], key [1], key [2], zIdx, + lastKey [0], lastKey [1], lastKey [2], lastzIdx + ); + } + lastk = k; + lastzIdx = zIdx; + VCOPY(lastKey, key); #endif /* Remove record from queue */ - QueuePop(queue, NULL); + QueuePop(queue, NULL); } /* Refill queue for next node(s) */ if (QueueFill(queue) < 0) { fputs("OOC_Build: failed input queue fill\n", stderr); return OOC_DATAIDX_ERR; } } return OOC_DATACNT(node); -} +} OOC_Octree *OOC_Build (OOC_Octree *oct, unsigned leafMax, unsigned maxDepth) /* Bottom-up construction of out-of-core octree in postorder traversal. The * octree oct is assumed to be initialised with its origin (oct -> org), * size (oct -> size), key callback (oct -> key), and its associated leaf * file (oct -> leafFile). * Records are read from the leafFile and assumed to be sorted in Z-order, * which defines an octree leaf ordering. Leaves (terminal nodes) are * constructed such that they contain <= leafMax records and have a maximum * depth of maxDepth. * Note that the following limits apply: * leafMax <= OOC_OCTCNT_MAX (see oococt.h) * maxDepth <= OOC_MORTON_BITS (see oocsort.h) * The maxDepth limit arises from the fact that the Z-ordering has a limited * resolution and will map node coordinates beyond a depth of * OOC_MORTON_BITS to the same Z-index, causing nodes to be potentially read * out of sequence and winding up in the wrong octree nodes. * On success, the octree pointer oct is returned, with the constructed * nodes in oct -> nodes, and the node count in oct -> numNodes. On * failure, NULL is returned. */ { OOC_BuildQueue queue; OOC_Node root; if (!oct || !oct -> size) { fputs("OOC_Build: octree not initialised", stderr); return NULL; } if (!oct -> leafFile) { fputs("OOC_Build: empty leaf file", stderr); return NULL; } oct -> leafMax = leafMax; oct -> maxDepth = maxDepth; queue.in = oct -> leafFile; rewind(queue.in); /* Init queue and fill from leaf file */ if (!QueueInit(&queue, oct -> recSize, leafMax + 1) || QueueFill(&queue) < 0) { fputs("OOC_Build: failed input queue init\n", stderr); return NULL; } /* Clear octree root and recurse */ OOC_CLEARNODE(&root); if (OOC_BuildRecurse(oct, &root, oct -> org, oct -> size, 0, &queue) == OOC_DATAIDX_ERR) return NULL; /* Finalise octree root */ if (!NewNode(oct)) { fputs("OOC_Build: failed to allocate octree root\n", stderr); return NULL; } OOC_SETROOT(oct, &root); /* Passing OOC_ROOT(oct) avoids annoying compiler warnings about (&root) * always evaluating to true when calling OOC_DATAIDX() */ oct -> numData = OOC_DATACNT(OOC_ROOT(oct)); return oct; } #endif /* NIX / PMAP_OOC */ + diff --git a/oocmorton.c b/oocmorton.c deleted file mode 100644 index 49be8e1..0000000 --- a/oocmorton.c +++ /dev/null @@ -1,44 +0,0 @@ -#ifndef lint -static const char RCSid[] = "$Id: oocmorton.c,v 2.2 2017/08/14 21:12:10 rschregle Exp $"; -#endif - - -/* - ========================================================================= - Routines to generate and compare Morton Codes, i.e. indices on space - filling Z-curves. - - 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: oocmorton.c,v 2.2 2017/08/14 21:12:10 rschregle Exp $ -*/ - - -#if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) -/* No Windoze support for now */ - -#include "oocmorton.h" - - - -OOC_MortonIdx OOC_Key2Morton (const FVECT key, const FVECT org, RREAL scale) -/* Compute Morton code (Z-curve index) of length 3 * OOC_MORTON_BITS bits - * for 3D key within bounding box defined by org and scaled to maximum index - * with scale */ -{ - unsigned i; - OOC_MortonIdx k [3]; - - /* Normalise key and map each dim to int of OOC_MORTON_BITS */ - for (i = 0; i < 3; i++) - k [i] = scale * (key [i] - org [i]); - - /* Interleave each dim with zeros and merge */ - return OOC_BitInterleave(k [0]) | OOC_BitInterleave(k [1]) << 1 | - OOC_BitInterleave(k [2]) << 2; -} - -#endif /* NIX / PMAP_OOC */ diff --git a/oocmorton.h b/oocmorton.h deleted file mode 100644 index bea3b16..0000000 --- a/oocmorton.h +++ /dev/null @@ -1,53 +0,0 @@ -/* - ====================================================================== - Routines to generate and compare Morton Codes, i.e. indices on space - filling Z-curves. - - 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: oocmorton.h,v 2.1 2016/05/17 17:39:47 rschregle Exp $ -*/ - -#ifndef OOC_MORTON_H - #define OOC_MORTON_H - - #include "fvect.h" - #include - - - /* Type to represent Morton codes (Z-curve indices) as sort keys */ - typedef uint64_t OOC_MortonIdx; - - - /* Number of bits per dimension for Morton code (Z-curve index) used to - * sort records. This corresponds to the maximum number of bounding box - * subdivisions which can be resolved per axis; further subdivisions will - * map to the same Morton code, thus invalidating the sort! */ - #define OOC_MORTON_BITS 21 - #define OOC_MORTON_MAX (((OOC_MortonIdx)1 << OOC_MORTON_BITS) - 1) - - - /* Interleave lower OOC_MORTON_BITS bits of k with 00, resulting in - * 3*OOC_MORTON_BITS bits. Optimised "bitmask hack" version. This code - * taken from: - * http://www.forceflow.be/2013/10/07/ - * morton-encodingdecoding-through-bit-interleaving-implementations/ */ - #define OOC_BitInterleave(k) \ - ((k) &= OOC_MORTON_MAX, \ - (k) = ((k) | (k) << 32) & 0x001f00000000ffff, \ - (k) = ((k) | (k) << 16) & 0x001f0000ff0000ff, \ - (k) = ((k) | (k) << 8) & 0x100f00f00f00f00f, \ - (k) = ((k) | (k) << 4) & 0x10c30c30c30c30c3, \ - (k) = ((k) | (k) << 2) & 0x1249249249249249) - - - /* Compute Morton code (Z-curve index) of length 3 * OOC_MORTON_BITS bits - * for 3D key within bounding box defined by org and scaled to maximum - * index with scale */ - OOC_MortonIdx OOC_Key2Morton (const FVECT key, const FVECT org, - RREAL scale); -#endif - \ No newline at end of file diff --git a/oococt.c b/oococt.c index 452f291..307b6f8 100644 --- a/oococt.c +++ b/oococt.c @@ -1,467 +1,482 @@ #ifndef lint static const char RCSid[] = "$Id: oococt.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $"; #endif /* ====================================================================== Out-of-core octree data structure Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF, #147053) ====================================================================== $Id: oococt.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $ */ #if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) /* No Windoze support for now */ #include "oococt.h" #include "rtio.h" #include #include void OOC_Null (OOC_Octree *oct) /* Init empty octree */ { oct -> maxNodes = oct -> numNodes = oct -> leafMax = oct -> maxDepth = oct -> numData = oct -> size = oct -> recSize = oct -> mortonScale = 0; oct -> org [0] = oct -> org [1] = oct -> org [2] = oct -> bound [0] = oct -> bound [1] = oct -> bound [2] = 0; oct -> key = NULL; oct -> nodes = NULL; oct -> leafFile = NULL; oct -> cache = NULL; } void OOC_Init (OOC_Octree *oct, unsigned recSize, FVECT org, RREAL size, - RREAL *(*key)(const void*), FILE *leafFile) + RREAL *(*key)(const void*), FILE *leafFile +) { oct -> maxNodes = oct -> numNodes = oct -> leafMax = oct -> maxDepth = oct -> numData = 0; VCOPY(oct -> org, org); oct -> size = oct -> bound[0] = oct -> bound[1] = oct -> bound[2] = size; VADD(oct -> bound, oct -> bound, org); - oct -> mortonScale = size > 0 ? OOC_MORTON_MAX / size : 0; + oct -> mortonScale = size > 0 ? MORTON3D_MAX / size : 0; oct -> recSize = recSize; oct -> key = key; oct -> nodes = NULL; oct -> leafFile = leafFile; oct -> cache = NULL; /* Cache currently initialised externally */ } OOC_Node *NewNode (OOC_Octree *oct) /* Allocate new octree node, enlarge array if necessary. Return pointer to new node or NULL if failed. */ { OOC_Node *n = NULL; if (oct -> numNodes >= OOC_NODEIDX_MAX) { /* Node index overflow */ fprintf(stderr, "OOC_NewNode: node index overflow (numNodes = %d)\n", - oct -> numNodes); + oct -> numNodes + ); return NULL; } if (++oct -> numNodes > oct -> maxNodes) { /* Need to allocate root or enlarge array */ oct -> maxNodes += OOC_BLKSIZE / sizeof(OOC_Node); oct -> nodes = realloc(oct -> nodes, - oct -> maxNodes * sizeof(OOC_Node)); + oct -> maxNodes * sizeof(OOC_Node) + ); if (!oct -> nodes) { perror("OOC_NewNode: couldn't realloc() nodes"); return NULL; } } n = oct -> nodes + oct -> numNodes - 1; n -> node.num = n -> node.kid = n -> node.branch = n -> node.type = 0; return n; } int OOC_GetData (OOC_Octree *oct, OOC_DataIdx idx, void *rec) /* Reads indexed data record from leaf file and copies it to rec, else * returns -1 on failure */ { if (!oct || !rec) return -1; if (idx >= oct -> numData) { fprintf(stderr, "OOC_GetData: invalid data record index %u\n", idx); return -1; } if (oct -> cache) { /* Retrieve record from leaf file via I/O cache */ void *cachedRec; if (!(cachedRec = OOC_CacheData(oct -> cache, oct -> leafFile, idx))) return -1; /* It's safer to copy the record to the caller's local buffer as a * returned pointer may be invalidated by subsequent calls if the * page it points to is swapped out */ memcpy(rec, cachedRec, oct -> recSize); } else { /* No I/O caching; do (SLOW!) per-record I/O from leaf file */ const unsigned long pos = (unsigned long)oct -> recSize * idx; if (pread(fileno(oct -> leafFile), rec, oct -> recSize, pos) != oct -> recSize) { perror("OOC_GetData: failed seek/read in leaf file"); return -1; } } return 0; } int OOC_InBBox (const OOC_Octree *oct, const FVECT org, RREAL size, - const FVECT key) + const FVECT key +) /* Check whether key lies within bounding box defined by org and size. * Compares Morton codes rather than the coordinates directly to account * for dicretisation error introduced by the former. */ { - const OOC_MortonIdx keyZ = OOC_KEY2MORTON(key, oct); - FVECT bbox = {size, size, size}; + const MortonIdx keyZ = OOC_KEY2MORTON(key, oct); + FVECT bbox = {size, size, size}; VADD(bbox, org, bbox); - return keyZ > OOC_KEY2MORTON(org, oct) && - keyZ < OOC_KEY2MORTON(bbox, oct); + return keyZ > OOC_KEY2MORTON(org, oct) && keyZ < OOC_KEY2MORTON(bbox, oct); } int OOC_Branch (const OOC_Octree *oct, const FVECT org, RREAL size, - const FVECT key, FVECT nuOrg) + const FVECT key, FVECT nuOrg +) /* Return index of octant containing key and corresponding origin if nuOrg * != NULL, or -1 if key lies outside all octants. Size is already assumed * to be halved, i.e. corresponding to the length of an octant axis. * Compares Morton codes rather than the coordinates directly to account * for dicretisation error introduced by the former. */ { int o; FVECT octOrg; for (o = 0; o < 8; o++) { VCOPY(octOrg, org); OOC_OCTORIGIN(octOrg, o, size); if (OOC_InBBox(oct, octOrg, size, key)) { if (nuOrg) VCOPY(nuOrg, octOrg); return o; } } return -1; } int OOC_BranchZ (const OOC_Octree *oct, const FVECT org, RREAL size, - OOC_MortonIdx keyZ, FVECT nuOrg) + MortonIdx keyZ, FVECT nuOrg +) /* Optimised version of OOC_Branch() with precomputed key Morton code */ { int o; const FVECT cent = {size, size, size}; FVECT octOrg, bbox; for (o = 0; o < 8; o++) { VCOPY(octOrg, org); OOC_OCTORIGIN(octOrg, o, size); VADD(bbox, octOrg, cent); if (keyZ > OOC_KEY2MORTON(octOrg, oct) && - keyZ < OOC_KEY2MORTON(bbox, oct)) { + keyZ < OOC_KEY2MORTON(bbox, oct) + ) { if (nuOrg) VCOPY(nuOrg, octOrg); return o; } } return -1; } OOC_DataIdx OOC_GetKid (const OOC_Octree *oct, OOC_Node **node, unsigned kid) /* For a given octree node, return the sum of the data counters for kids * [0..k-1]. On return, the node either points to the k'th kid, or * NULL if the kid is nonexistent or the node is a leaf */ { unsigned k; OOC_Node *kidNode = OOC_KID1(oct, *node); /* NULL if leaf */ OOC_DataIdx dataIdx = 0; for (k = 0; k < kid; k++) { if (OOC_ISLEAF(*node)) /* Sum data counters of leaf octants */ dataIdx += (*node) -> leaf.num [k]; else /* Sum data counter of inner node's nonempty kids and advance * pointer to sibling */ if (OOC_ISBRANCH(*node, k)) { dataIdx += OOC_DATACNT(kidNode); kidNode++; } } /* Update node pointer to selected kid (or NULL if nonexistent or node is * a leaf) */ *node = OOC_ISBRANCH(*node, kid) ? kidNode : NULL; return dataIdx; } #ifdef DEBUG_OOC -int OOC_Check (OOC_Octree *oct, const OOC_Node *node, - FVECT org, RREAL size, OOC_DataIdx dataIdx) -/* Traverse tree & check for valid nodes; oct -> leafFile must be open; - * return 0 if ok, otherwise -1 */ -{ - unsigned k; - FVECT kidOrg; - const RREAL kidSize = size * 0.5; - - if (!oct || !node) - return -1; - - if (OOC_ISLEAF(node)) { - /* Node is a leaf; check records in each octant */ - char rec [oct -> recSize]; /* Violates C89! */ - OOC_OctCnt d; - RREAL *key; - - for (k = 0; k < 8; k++) { - VCOPY(kidOrg, org); - OOC_OCTORIGIN(kidOrg, k, kidSize); + int OOC_Check (OOC_Octree *oct, const OOC_Node *node, + FVECT org, RREAL size, OOC_DataIdx dataIdx + ) + /* Traverse tree & check for valid nodes; oct -> leafFile must be open; + * return 0 if ok, otherwise -1 */ + { + unsigned k; + FVECT kidOrg; + const RREAL kidSize = size * 0.5; + + if (!oct || !node) + return -1; - for (d = node -> leaf.num [k]; d; d--, dataIdx++) { -#ifdef DEBUG_OOC_CHECK - fprintf(stderr, "OOC_Check: checking record %lu\n", - (unsigned long)dataIdx); -#endif - if (OOC_GetData(oct, dataIdx, rec)) { - fprintf(stderr, "OOC_Check: failed getting record at %lu\n", - (unsigned long)dataIdx); - return -1; - } + if (OOC_ISLEAF(node)) { + /* Node is a leaf; check records in each octant */ + char rec [oct -> recSize]; /* Violates C89! */ + OOC_OctCnt d; + RREAL *key; + + for (k = 0; k < 8; k++) { + VCOPY(kidOrg, org); + OOC_OCTORIGIN(kidOrg, k, kidSize); - key = oct -> key(rec); - if (!OOC_InBBox(oct, kidOrg, kidSize, key)) { - fprintf(stderr, "OOC_Check: key [%f, %f, %f] (zIdx %020lu) " - "outside bbox [[%f-%f], [%f-%f], [%f-%f]] " - "in octant %d (should be %d)\n", - key [0], key [1], key [2], - OOC_KEY2MORTON(key, oct), - kidOrg [0], kidOrg [0] + kidSize, - kidOrg [1], kidOrg [1] + kidSize, - kidOrg [2], kidOrg [2] + kidSize, - k, OOC_Branch(oct, org, kidSize, key, NULL)); - /* return -1; */ + for (d = node -> leaf.num [k]; d; d--, dataIdx++) { + #ifdef DEBUG_OOC_CHECK + fprintf(stderr, "OOC_Check: checking record %lu\n", + (unsigned long)dataIdx + ); + #endif + if (OOC_GetData(oct, dataIdx, rec)) { + fprintf(stderr, "OOC_Check: failed getting record at %lu\n", + (unsigned long)dataIdx + ); + return -1; + } + + key = oct -> key(rec); + if (!OOC_InBBox(oct, kidOrg, kidSize, key)) { + fprintf(stderr, "OOC_Check: key [%f, %f, %f] (zIdx %020lu) " + "outside bbox [[%f-%f], [%f-%f], [%f-%f]] " + "in octant %d (should be %d)\n", + key [0], key [1], key [2], + OOC_KEY2MORTON(key, oct), + kidOrg [0], kidOrg [0] + kidSize, + kidOrg [1], kidOrg [1] + kidSize, + kidOrg [2], kidOrg [2] + kidSize, + k, OOC_Branch(oct, org, kidSize, key, NULL) + ); + /* return -1; */ + } } } } - } - else { - /* Internal node; recurse on all kids */ - const OOC_Node *kid = OOC_KID1(oct, node); - OOC_DataIdx numData = 0; - - if (node -> node.kid >= oct -> numNodes) { - fprintf(stderr, "OOC_Check: invalid node index %u\n", - node -> node.kid); - return -1; - } + else { + /* Internal node; recurse on all kids */ + const OOC_Node *kid = OOC_KID1(oct, node); + OOC_DataIdx numData = 0; - if (!node -> node.num) { - fputs("OOC_Check: empty octree node\n", stderr); - return -1; - } - - for (k = 0; k < 8; k++) - if (OOC_ISBRANCH(node, k)) { - VCOPY(kidOrg, org); - OOC_OCTORIGIN(kidOrg, k, kidSize); - - if (OOC_Check(oct, kid, kidOrg, kidSize, dataIdx)) - return -1; + if (node -> node.kid >= oct -> numNodes) { + fprintf(stderr, "OOC_Check: invalid node index %u\n", + node -> node.kid + ); + return -1; + } - dataIdx += OOC_DATACNT(kid); - numData += OOC_DATACNT(kid); - kid++; + if (!node -> node.num) { + fputs("OOC_Check: empty octree node\n", stderr); + return -1; + } + + for (k = 0; k < 8; k++) + if (OOC_ISBRANCH(node, k)) { + VCOPY(kidOrg, org); + OOC_OCTORIGIN(kidOrg, k, kidSize); + + if (OOC_Check(oct, kid, kidOrg, kidSize, dataIdx)) + return -1; + + dataIdx += OOC_DATACNT(kid); + numData += OOC_DATACNT(kid); + kid++; + } + + if (node -> node.num != numData) { + fprintf(stderr, + "Parent/kid data count mismatch; expected %u, found %u\n", + node -> node.num, numData + ); + return -1; } - - if (node -> node.num != numData) { - fprintf(stderr, - "Parent/kid data count mismatch; expected %u, found %u\n", - node -> node.num, numData); - return -1; } + + return 0; } - - return 0; -} #endif - + int OOC_SaveOctree (const OOC_Octree *oct, FILE *out) /* Appends octree nodes to specified file along with metadata. Uses * RADIANCE's portable I/O routines. Returns 0 on success, else -1. Note * the internal representation of the nodes is platform specific as they * contain unions, therefore we use the portable I/O routines from the * RADIANCE library */ { int i; OOC_NodeIdx nc; OOC_Node *np = NULL; if (!oct || !oct->nodes || !oct->numData || !oct->numNodes || !out) { fputs("OOC_SaveOctree: empty octree\n", stderr); return -1; } /* Write octree origin, size, data record size, number of records, and * number of nodes */ for (i = 0; i < 3; i++) putflt(oct -> org [i], out); putflt(oct -> size, out); putint(oct -> recSize, sizeof(oct -> recSize), out); putint(oct -> numData, sizeof(oct -> numData), out); putint(oct -> numNodes, sizeof(oct -> numNodes), out); /* Write nodes to file */ for (nc = oct -> numNodes, np = oct -> nodes; nc; nc--, np++) { if (OOC_ISLEAF(np)) { /* Write leaf marker */ putint(-1, 1, out); /* Write leaf octant counters */ for (i = 0; i < 8; i++) putint(np -> leaf.num [i], sizeof(np -> leaf.num [i]), out); } else { /* Write node marker */ putint(0, 1, out); /* Write node, rounding field size up to nearest number of bytes */ putint(np -> node.kid, (OOC_NODEIDX_BITS + 7) >> 3, out); putint(np -> node.num, (OOC_DATAIDX_BITS + 7) >> 3, out); putint(np -> node.branch, 1, out); } if (ferror(out)) { fputs("OOC_SaveOctree: failed writing to node file", stderr); return -1; } } fflush(out); return 0; } int OOC_LoadOctree (OOC_Octree *oct, FILE *nodeFile, - RREAL *(*key)(const void*), FILE *leafFile) + RREAL *(*key)(const void*), FILE *leafFile +) /* Load octree nodes and metadata from nodeFile and associate with records * in leafFile. Uses RADIANCE's portable I/O routines. Returns 0 on * success, else -1. */ { int i; OOC_NodeIdx nc; OOC_Node *np = NULL; if (!oct || !nodeFile) return -1; OOC_Null(oct); /* Read octree origin, size, record size, #records and #nodes */ for (i = 0; i < 3; i++) oct -> org [i] = getflt(nodeFile); oct -> size = getflt(nodeFile); oct -> bound [0] = oct -> bound [1] = oct -> bound [2] = oct -> size; VADD(oct -> bound, oct -> bound, oct -> org); - oct -> mortonScale = OOC_MORTON_MAX / oct -> size; + oct -> mortonScale = MORTON3D_MAX / oct -> size; oct -> recSize = getint(sizeof(oct -> recSize), nodeFile); oct -> numData = getint(sizeof(oct -> numData), nodeFile); oct -> numNodes = getint(sizeof(oct -> numNodes), nodeFile); oct -> key = key; oct -> leafFile = leafFile; if (feof(nodeFile) || !oct -> numData || !oct -> numNodes) { fputs("OOC_LoadOctree: empty octree\n", stderr); return -1; } /* Allocate octree nodes */ if (!(oct -> nodes = calloc(oct -> numNodes, sizeof(OOC_Node)))) { fputs("OOC_LoadOctree: failed octree allocation\n", stderr); return -1; } /* Read nodes from file */ for (nc = oct -> numNodes, np = oct -> nodes; nc; nc--, np++) { OOC_CLEARNODE(np); /* Read node type marker */ if (getint(1, nodeFile)) { /* Read leaf octant counters and set node type */ for (i = 0; i < 8; i++) np -> leaf.num [i] = getint(sizeof(np -> leaf.num [i]), - nodeFile); - + nodeFile + ); + OOC_SETLEAF(np); } else { /* Read node, rounding field size up to nearest number of bytes */ np -> node.kid = getint((OOC_NODEIDX_BITS + 7) >> 3, nodeFile); np -> node.num = getint((OOC_DATAIDX_BITS + 7) >> 3, nodeFile); np -> node.branch = getint(1, nodeFile); } if (ferror(nodeFile)) { fputs("OOC_LoadOctree: failed reading from node file\n", stderr); return -1; } } return 0; } void OOC_Delete (OOC_Octree *oct) /* Self-destruct */ { free(oct -> nodes); fclose(oct -> leafFile); if (oct -> cache) OOC_DeleteCache(oct -> cache); } #endif /* NIX / PMAP_OOC */ + diff --git a/oococt.h b/oococt.h index b930328..360f22e 100644 --- a/oococt.h +++ b/oococt.h @@ -1,234 +1,256 @@ /* ====================================================================== Out-of-core octree data structure Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF, #147053) ====================================================================== $Id: oococt.h,v 2.3 2016/05/17 17:39:47 rschregle Exp $ */ -#ifndef OOC_OCT_H -#define OOC_OCT_H +#ifndef _OOC_OCT_H +#define _OOC_OCT_H - #include "oocmorton.h" + #include "morton.h" #include "ooccache.h" #include #include /* =================================================================== * CONSTS * =================================================================== */ /* Octree index & data counter types */ typedef uint32_t OOC_NodeIdx; typedef uint32_t OOC_DataIdx; typedef uint8_t OOC_OctCnt; /* Octree node field sizes for above */ #define OOC_NODEIDX_BITS 31 #define OOC_DATAIDX_BITS 32 /* Limits for above */ #define OOC_NODEIDX_MAX ((1L << OOC_NODEIDX_BITS) - 1) #define OOC_DATAIDX_MAX ((1L << OOC_DATAIDX_BITS) - 2) #define OOC_OCTCNT_MAX UINT8_MAX /* Reserved values for above */ #define OOC_DATAIDX_ERR (OOC_DATAIDX_MAX + 1) /* Block size for node allocation (8k) */ #define OOC_BLKSIZE (1L << 13) /* =================================================================== * UTILZ * =================================================================== */ /* Test/set leaf flag for node n */ #define OOC_ISLEAF(n) ((n) -> node.type) #define OOC_SETLEAF(n) ((n) -> node.type = 1) /* Test/set branch bit for octant k in node n */ - #define OOC_ISBRANCH(n,k) (!OOC_ISLEAF(n) && (n) -> node.branch & 1 << (k)) + #define OOC_ISBRANCH(n,k) ( \ + !OOC_ISLEAF(n) && (n) -> node.branch & 1 << (k) \ + ) #define OOC_SETBRANCH(n,k) ((n) -> node.branch |= 1 << (k)) /* Return index to node's 1st kid */ - #define OOC_KID1(o,n) (!OOC_ISLEAF(n) ? (o)->nodes + (n)->node.kid \ - : NULL) - + #define OOC_KID1(o,n) (!OOC_ISLEAF(n) \ + ? (o)->nodes + (n)->node.kid \ + : NULL \ + ) + /* Get index to last node in octree; -1 indicates empty octree */ #define OOC_ROOTIDX(o) ((o) -> nodes && (o) -> numNodes > 0 \ - ? (o) -> numNodes - 1 : -1) - + ? (o) -> numNodes - 1 \ + : -1 \ + ) + /* Get octree root; this is at the *end* of the node array! */ #define OOC_ROOT(o) (OOC_ROOTIDX(o) >= 0 \ - ? (o) -> nodes + OOC_ROOTIDX(o) : NULL) + ? (o) -> nodes + OOC_ROOTIDX(o) \ + : NULL \ + ) /* Copy node to octree root */ #define OOC_SETROOT(o,n) (memcpy(OOC_ROOT(o), (n), sizeof(OOC_Node))) /* Set all node fields to 0 */ #define OOC_CLEARNODE(n) (memset((n), 0, sizeof(OOC_Node))) /* Modify origin o for octant k of size s */ - #define OOC_OCTORIGIN(o,k,s) ((o) [0] += (s) * ((k) & 1), \ - (o) [1] += (s) * ((k) >> 1 & 1), \ - (o) [2] += (s) * ((k) >> 2 & 1)) + #define OOC_OCTORIGIN(o,k,s) ( \ + (o) [0] += (s) * ((k) & 1), \ + (o) [1] += (s) * ((k) >> 1 & 1), \ + (o) [2] += (s) * ((k) >> 2 & 1) \ + ) /* Get num records stored in node (also handles leaves and NULL pointer) */ - #define OOC_DATACNT(n) \ - ((n) ? OOC_ISLEAF(n) ? (n) -> leaf.num [0] + (n) -> leaf.num [1] + \ - (n) -> leaf.num [2] + (n) -> leaf.num [3] + \ - (n) -> leaf.num [4] + (n) -> leaf.num [5] + \ - (n) -> leaf.num [6] + (n) -> leaf.num [7] \ - : (n) -> node.num \ - : 0) + #define OOC_DATACNT(n) ((n) \ + ? OOC_ISLEAF(n) \ + ? (n) -> leaf.num [0] + (n) -> leaf.num [1] + \ + (n) -> leaf.num [2] + (n) -> leaf.num [3] + \ + (n) -> leaf.num [4] + (n) -> leaf.num [5] + \ + (n) -> leaf.num [6] + (n) -> leaf.num [7] \ + : (n) -> node.num \ + : 0 \ + ) /* Shortcuts for morton code from key and data record */ - #define OOC_KEY2MORTON(k, o) OOC_Key2Morton((k), (o) -> org, \ - (o) -> mortonScale) - #define OOC_REC2MORTON(r, o) OOC_Key2Morton((o) -> key(r), (o) -> org, \ - (o) -> mortonScale) + #define OOC_KEY2MORTON(k, o) Key2Morton3D( \ + (k), (o) -> org, (o) -> mortonScale \ + ) + #define OOC_REC2MORTON(r, o) Key2Morton3D( \ + (o) -> key(r), (o) -> org, (o) -> mortonScale \ + ) /* =================================================================== * DATA TYPEZ * =================================================================== */ /* Octree node */ typedef struct { union { struct { /* Interior node (node.type = 0) with: node.kid = index to 1st child node in octree array (a.k.a "Number One Son"), immediately followed by its nonzero siblings as indicated by branch node.num = num records stored in this subtree, node.branch = branch bitmask indicating nonzero kid nodes */ char type : 1; - OOC_NodeIdx kid : OOC_NODEIDX_BITS; + OOC_NodeIdx kid : OOC_NODEIDX_BITS; OOC_DataIdx num : OOC_DATAIDX_BITS; uint8_t branch; /* NOTE that we can't make the kid's index relative to parent * (which would be more compact), since we don't know the * parent's index a priori during the bottom-up construction! */ } node; struct { /* Leaf node (leaf.type = node.type = 1 with: leaf.num [k] = num records stored in octant k */ char type : 1; OOC_OctCnt num [8]; } leaf; }; } OOC_Node; /* Top level octree container */ typedef struct { FVECT org, bound; /* Cube origin (min. XYZ), size, and resulting bounds (max. XYZ) */ RREAL size, *(*key)(const void*); /* Callback for data rec coords */ RREAL mortonScale; /* Scale factor for generating Morton codes from coords */ OOC_Node *nodes; /* Pointer to node array */ /* ************************************* * **** NODES ARE STORED IN POSTFIX **** * **** ORDER, I.E. ROOT IS LAST! **** * *************************************/ OOC_DataIdx numData; /* Num records in leaves */ OOC_NodeIdx maxNodes, /* Num currently allocated nodes */ numNodes; /* Num used nodes (<= maxNodes) */ unsigned recSize, /* Size of data record in bytes */ leafMax, /* Max #records per leaf (for build) */ maxDepth; /* Max subdivision depth (for build) */ FILE *leafFile; /* File containing data in leaves */ OOC_Cache *cache; /* I/O cache for leafFile */ } OOC_Octree; /* DOAN' FORGET: NODES IN ***POSTFIX*** ORDAH!!! */ /* =================================================================== * FUNKZ * =================================================================== */ /* Init empty octree */ void OOC_Null (OOC_Octree *oct); /* Initialises octree for records of size recSize with 3D coordinates, * accessed via key() callback. The octree's bounding box is defined by * its origin org and size per axis, and must contain all keys. The * octree is associated with the specified leaf file containing the * records in Morton order */ void OOC_Init (OOC_Octree *oct, unsigned recSize, FVECT org, RREAL size, - RREAL *(*key)(const void*), FILE *leafFile); + RREAL *(*key)(const void*), FILE *leafFile + ); /* Allocate new octree node, enlarge array if necessary. Returns pointer to new node or NULL if failed. */ OOC_Node *NewNode (OOC_Octree *oct); /* Reads indexed data record from leaf file and copies it to rec, else * returns -1 on failure */ int OOC_GetData (OOC_Octree *oct, OOC_DataIdx idx, void *rec); - + /* Appends octree nodes to specified file along with metadata. Uses * RADIANCE's portable I/O routines. Returns 0 on success, else -1. */ int OOC_SaveOctree (const OOC_Octree *oct, FILE *nodeFile); /* Load octree nodes and metadata from nodeFile and associate with * records in leafFile. Uses RADIANCE's portable I/O routines. Returns * 0 on success, else -1. */ int OOC_LoadOctree (OOC_Octree *oct, FILE *nodeFile, - RREAL *(*key)(const void*), FILE *leafFile); + RREAL *(*key)(const void*), FILE *leafFile + ); #ifdef DEBUG_OOC /* Traverse tree & check for valid nodes; oct -> leafFile must be open; * return 0 if ok, otherwise -1 */ int OOC_Check (OOC_Octree *oct, const OOC_Node *node, - FVECT org, RREAL size, OOC_DataIdx dataIdx); + FVECT org, RREAL size, OOC_DataIdx dataIdx + ); #endif /* Check whether key lies within bounding box defined by org and size */ int OOC_InBBox (const OOC_Octree *oct, const FVECT org, RREAL size, - const FVECT key); + const FVECT key + ); /* Return index of octant containing key and corresponding origin if * nuOrg = NULL, or -1 if key lies outside all octants. Size is already * assumed to be halved, i.e. corresponding to the length of an octant * axis. */ int OOC_Branch (const OOC_Octree *oct, const FVECT org, RREAL size, - const FVECT key, FVECT nuOrg); + const FVECT key, FVECT nuOrg + ); /* Optimised version of OOC_Branch() with precomputed key Morton code */ int OOC_BranchZ (const OOC_Octree *oct, const FVECT org, RREAL size, - OOC_MortonIdx keyZ, FVECT nuOrg); + MortonIdx keyZ, FVECT nuOrg + ); /* For a given octree node, return the sum of the data counters for kids * [0..k-1]. On return, the node either points to the k'th kid, or * NULL if the kid is nonexistent or the node is a leaf */ OOC_DataIdx OOC_GetKid (const OOC_Octree *oct, OOC_Node **node, - unsigned k); + unsigned k) + ; /* Self-destruct */ void OOC_Delete (OOC_Octree *oct); /* DID WE MENTION NODES ARE IN POSTFIX ORDAH? */ -#endif +#endif /* _OOC_OCT_H */ + diff --git a/oocsort.c b/oocsort.c index 77a02bb..8f2ea8f 100644 --- a/oocsort.c +++ b/oocsort.c @@ -1,581 +1,597 @@ #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 */ #include "oocsort.h" -#include "oocmorton.h" +#include "morton.h" #include #include #include #include #include /* Priority queue node */ typedef struct { - OOC_MortonIdx pri; /* Record's priority (sort key) */ + MortonIdx pri; /* Record's priority (sort key) */ unsigned blk; /* Block containing record */ } OOC_SortQueueNode; /* Priority queue */ typedef struct { OOC_SortQueueNode *node; unsigned len, tail; } OOC_SortQueue; /* 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; static int OOC_KeyCompare (const void *rec1, const void *rec2) /* Comparison function for in-core quicksort */ { - OOC_MortonIdx pri1, pri2; + MortonIdx pri1, pri2; if (!rec1 || !rec2) return 0; - pri1 = OOC_Key2Morton(keyData.key(rec1), keyData.bbOrg, - keyData.mortonScale); - pri2 = OOC_Key2Morton(keyData.key(rec2), keyData.bbOrg, - keyData.mortonScale); + 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; 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; /* 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_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 -> 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); } } return 0; } #endif -static int OOC_SortPush (OOC_SortQueue *pq, OOC_MortonIdx pri, unsigned blk) +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; /* 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 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; - OOC_MortonIdx pri; + 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; root = kid; } 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; } #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) + 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 */ - OOC_MortonIdx pri; + 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); + 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; wait if necessary */ while (procCnt >= numProc && wait(&stat) >= 0) { if (!WIFEXITED(stat) || WEXITSTATUS(stat)) return -1; procCnt--; } /* Now fork kid process */ if (!(pid = fork())) { /* Recurse on subblocks with new input filehandle; */ if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk, - maxBlkSize, numProc, recSize, sortBuf, - pqueue, keyData)) + 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; } #ifdef DEBUG_OOC_FORK fprintf(stderr, "OOC_Sort: forking kid %d for block %d\n", - procCnt, b); -#endif + procCnt, b + ); +#endif /* Parent continues here */ procCnt++; } else { /* Recurse on subblock; without forking */ if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk, - maxBlkSize, numProc, recSize, sortBuf, - pqueue, keyData)) + 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 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); + "file %d\n", b + ); return -1; } fprintf(stderr, "OOC_Sort: tmp block file %d contains %ld rec\n", - b, ftell(blkFile [b]) / recSize); + b, ftell(blkFile [b]) / recSize + ); #endif rewind(blkFile [b]); 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 = OOC_Key2Morton(keyData -> key(rec), keyData -> bbOrg, - keyData -> mortonScale); + 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 the input block */ - fprintf(stderr, "OOC_Sort: unexpected end reading block file\n"); + fprintf(stderr, + "OOC_Sort: unexpected end reading block file\n" + ); return -1; } 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 insert in priority queue */ if (!OOC_SortPeek(blkFile [b], recSize, rec)) { /* Block not exhausted */ - pri = OOC_Key2Morton(keyData -> key(rec), keyData -> bbOrg, - keyData -> mortonScale); + 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); + 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)); } 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); + 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; - } + } /* Quicksort block in-core and write to output file */ qsort(sortBuf, blkLen, recSize, OOC_KeyCompare); 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 - } + 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*)) +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; /* 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.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 = OOC_MORTON_MAX / bbSize; + keyData.mortonScale = MORTON3D_MAX / bbSize; VCOPY(keyData.bbOrg, bbOrg); stat = OOC_SortRecurse(in, 0, numRec - 1, out, numBlk, blkSize, numProc, - recSize, sortBuf, &pqueue, &keyData); + recSize, sortBuf, &pqueue, &keyData + ); /* Cleanup */ free(pqueue.node); free(sortBuf); - return stat; + return stat; } #endif /* NIX / PMAP_OOC */ + diff --git a/pmapooc.c b/pmapooc.c index 7906c21..1d88b67 100644 --- a/pmapooc.c +++ b/pmapooc.c @@ -1,349 +1,349 @@ /* ====================================================================== 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; - OOC_MortonIdx zIdx, lastzIdx = 0; + 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) ) 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 -> lastContribOrg if precomputing contribution photon */ filtState.srcMod = isContribPmap(pmap) ? findmaterial(source [pmap -> lastContribOrg.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/pmapooc.h b/pmapooc.h index 328b3bd..f6c9bd2 100644 --- a/pmapooc.h +++ b/pmapooc.h @@ -1,90 +1,89 @@ /* ================================================================== 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") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ================================================================== $Id: pmapooc.h,v 1.1 2020/10/29 16:20:11 u-no-hoo Exp u-no-hoo $ */ #ifndef PMAPOOC_H #define PMAPOOC_H #include "oocnn.h" /* Suffixes for octree filenames */ /* #define PMAP_OOC_NODESUFFIX ".node" #define PMAP_OOC_HEAPSUFFIX ".heap" */ #define PMAP_OOC_LEAFSUFFIX ".leaf" /* Out-of-core octree constants */ #define PMAP_OOC_NUMBLK 32 /* Num blocks for external sort */ #define PMAP_OOC_BLKSIZE 1e8 /* Block size for external sort */ #define PMAP_OOC_LEAFMAX (OOC_OCTCNT_MAX) /* Max photons per leaf */ - #define PMAP_OOC_MAXDEPTH (OOC_MORTON_BITS) /* Max octree depth */ + #define PMAP_OOC_MAXDEPTH (MORTON3D_BITS) /* Max octree depth */ typedef OOC_SearchQueueNode PhotonSearchQueueNode; typedef OOC_SearchQueue PhotonSearchQueue; typedef OOC_Octree PhotonStorage; typedef unsigned PhotonIdx; /* Forward declarations to break dependency loop with pmapdata.h */ struct PhotonMap; void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc); /* Build out-of-core octree pmap -> store from photons in unsorted * heapfile pmap -> heap and generate nodes and leaf file with prefix * pmap -> fileName. Photon map construction may be parallelised if * numProc > 1, if supported. The heap is destroyed on return. */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out); /* Save photons in out-of-core octree to file. Return -1 on error, * else 0 */ int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *in); /* Load photons in out-of-core octree from file. Return -1 on error, * else 0 */ void OOC_InitFindPhotons (struct PhotonMap *pmap); /* Initialise NN search queue prior to calling kdT_FindPhotons() */ int OOC_FindPhotons ( struct PhotonMap* pmap, const FVECT pos, const FVECT norm ); /* Locate pmap -> squeue.len nearest photons to pos with similar normal * (NULL for volume photons) and return in search queue pmap -> squeue, * starting with the further photon at pmap -> squeue.node. Return -1 * if none found, else 0. */ - int OOC_Find1Photon ( - struct PhotonMap* pmap, const FVECT pos, const FVECT norm, - Photon *photon + int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, + const FVECT norm, Photon *photon ); /* Locate single nearest photon to pos with similar normal. Return -1 * if none found, else 0. */ int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon); /* Retrieve photon referenced by idx from leaf file and return -1 on - * error, else 0. */ + * error, else 0. */ - Photon *OOC_GetNearestPhoton ( - const PhotonSearchQueue *squeue, PhotonIdx idx + Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, + PhotonIdx idx ); /* Retrieve photon from NN search queue after OOC_FindPhotons() */ PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap); - /* Return index to first photon in octree */ + /* Return index to first photon in octree */ #endif