diff --git a/Rmakefile b/Rmakefile index f921052..6bb0917 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,548 +1,548 @@ # 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 \ $(PMRCOBJS) 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 $(PMRCSRC) PMRCOBJS = pmapcontrib.o pmcontrib2.o pmcontrib3.o pmcontrib4.o \ pmcontribsort.o wavelet2.o mrgbe.o PMRCSRC = pmapcontrib.c pmcontrib2.c pmcontrib3.c pmcontrib4.c \ pmcontribsort.c wavelet2.c mrgbe.c HEADERS = ambient.h ray.h data.h otspecial.h source.h # # What this makefile produces: # PROGS = $(DESTDIR)/rtrace $(DESTDIR)/rpict $(DESTDIR)/rvu $(DESTDIR)/rcontrib \ $(DESTDIR)/lookamb $(DESTDIR)/mkpmap $(DESTDIR)/pmapdump \ 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 + -o wavelet-test -g wavelet2.c $(MLIB) 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 -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 $(RLIB) $(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 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/mrgbe.c b/mrgbe.c index d675d07..3356877 100644 --- a/mrgbe.c +++ b/mrgbe.c @@ -1,394 +1,396 @@ /* ========================================================================= mRGBE (miniRGBE): reduced RGBE encoding of normalised RGB floats plus associated integer payload data. See mrgbe.h for encoding details. Compile with -DmRGBE_TEST to enable unit test (requires libm) 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$ */ #include "mrgbe.h" #include #include #include #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define ABS(a) ((a) >= 0 ? (a) : -(a)) mRGBERange *mRGBEinit (mRGBERange *range, double rgbMin [3], double rgbMax [3] ) /* Initialise range state for mRGBEencode()/mRGBEdecode() with the specified * float input ranges [rgbMin, rgbMax] in _absolute_ values per colour * channel. This sets the offset and normalisation for the encoding. */ { unsigned i; for (i = 0; i < 3; i++) { if (rgbMin [i] >= rgbMax [i]) { fprintf(stderr, "ERROR - invalid mRGBE range\n"); return NULL; } range -> min [i] = fabs(rgbMin [i]); range -> max [i] = fabs(rgbMax [i]); #if 1 - /* Reduce minimum so sign is preserved after offset during encoding */ + /* Reduce minimum within mRGBE resolution to avoid erroneous sign + reversals after offset during encoding */ range -> min [i] = MAX(0, range -> min [i] - mRGBE_MIN); #endif range -> norm [i] = mRGBE_RANGE / (range -> max [i] - range -> min [i]); if (range -> norm [i] < mRGBE_MIN) { fprintf(stderr, "ERROR - Zero mRGBE normalisation\n"); return NULL; } } return range; } mRGBE mRGBEencode (double rgb [3], const mRGBERange *range, unsigned data) /* Encode signed float RGB tuple to mRGBE along with payload data. The * supplied range state must be previously initialised with mRGBEinit(). * A warning is issued if rgb lies outside the initialised range. */ { #define SGN(a, b) ((a) >= 0 ? (b) : -(b)) mRGBE mrgbe = {0, 0, 0, 0}; double rgbRange [3], rgbNorm [3], mantNorm; int i, exp; /* Check payload data */ if (data > mRGBE_DATAMAX) { fprintf(stderr, "ERROR - mRGBE data overflow"); return mrgbe; } if (!range) { fprintf(stderr, "ERROR - mRGBE range not initialised\n"); return mrgbe; } /* Offset and normalise RGB according to initialised range; note rgbNorm * is in ABSOLUTE VALUES */ for (i = 0; i < 3; i++) { rgbNorm [i] = fabs(rgb [i]); /* Check range */ if (rgbNorm [i] < range -> min [i] || rgbNorm [i] > range -> max [i]) { fprintf(stderr, "ERROR - mRGBE input outside range\n"); return mrgbe; } rgbNorm [i] = (rgbNorm [i] - range -> min [i]) * range -> norm [i]; if (fabs(rgbNorm [i] - mRGBE_MAX) < mRGBE_MIN) /* Normalisation to maximum. This must be avoided by offsetting by minimum, otherwise frexp() returns a positive exponent. Normalisation to >mRGBE_MAX is caught by overflow check below */ rgbNorm [i] = mRGBE_RANGE; } /* Mantissa normalisation factor (= RGB maximum) */ mantNorm = MAX(MAX(rgbNorm [0], rgbNorm [1]), rgbNorm [2]); /* Get normalised mantissa and exponent to base 2 */ mantNorm = mantNorm < mRGBE_MIN ? 0 : frexp(mantNorm, &exp) * mRGBE_MANTMAX / mantNorm; /* Set RGB integer mantissas, adding excess (zero offset) and rounding up or down, depending on sign. Set sign according to original value. Clamp to mRGBE_MANTOVF to prevent rounding beyond overflow into negative! */ mrgbe.red = MIN((int)( SGN(rgb [0], rgbNorm [0] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX ), mRGBE_MANTOVF ); mrgbe.grn = MIN((int)( SGN(rgb [1], rgbNorm [1] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX ), mRGBE_MANTOVF ); mrgbe.blu = MIN((int)( SGN(rgb [2], rgbNorm [2] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX ), mRGBE_MANTOVF ); /* Drop the exponent sign, since implicitly negative */ mrgbe.exp = abs(exp); /* Set payload data */ mrgbe.dat = data; return mrgbe; #undef SGN } unsigned mRGBEdecode (mRGBE mrgbe, const mRGBERange *range, double rgb [3]) /* Decode mRGBE, returning signed float RGB tuple in array rgb. The * associated payload data is decoded as return value. The supplied range * state must be previously initialised with mRGBEinit(). */ { #define SGN(a, b) ((a) >= mRGBE_MANTMAX ? (b) : -(b)) #define JITTER (mRGBE_JITTER * drand48()) double mantNorm; if (!range) { fprintf(stderr, "ERROR - mRGBE range not initialised\n"); return 0; } /* Mantissa normalisation factor; exponent is implicitly negative */ mantNorm = ldexp(1.0, -mrgbe.exp) / mRGBE_MANTMAX; /* Decode and set RGB components, subtracting excess (zero offset) and * adding jitter to break up quantisation artefacts. Clamp to -mRGBE_MAX * in case jitter underflows. (Note jitter is added/subtracted for * negative/positive values to compensate rounding down/up during * encoding). Finally denormalise and offset to original range with * saved state. */ /* rgb [0] = MAX( (mrgbe.red - mRGBE_MANTMAX - SGN(mrgbe.red, JITTER)) * mantNorm / range -> norm [0] + SGN(mrgbe.red, range -> min [0]), -mRGBE_MAX / range -> norm [0] ); */ /* rgb [0] = (mrgbe.red - mRGBE_MANTMAX - SGN(mrgbe.red, JITTER)) * mantNorm / range -> norm [0] + SGN(mrgbe.red, range -> min [0]); rgb [1] = (mrgbe.grn - mRGBE_MANTMAX - SGN(mrgbe.grn, JITTER)) * mantNorm / range -> norm [1] + SGN(mrgbe.grn, range -> min [1]); rgb [2] = (mrgbe.blu - mRGBE_MANTMAX - SGN(mrgbe.blu, JITTER)) * mantNorm / range -> norm [2] + SGN(mrgbe.blu, range -> min [2]); */ rgb [0] = SGN(mrgbe.red, 1) * ( ABS((signed)mrgbe.red - mRGBE_MANTMAX + JITTER) * mantNorm / range -> norm [0] + range -> min [0] ); rgb [1] = SGN(mrgbe.grn, 1) * ( ABS((signed)mrgbe.grn - mRGBE_MANTMAX + JITTER) * mantNorm / range -> norm [1] + range -> min [1] ); rgb [2] = SGN(mrgbe.blu, 1) * ( ABS((signed)mrgbe.blu - mRGBE_MANTMAX + JITTER) * mantNorm / range -> norm [2] + range -> min [2] ); /* Return payload data */ return mrgbe.dat; #undef SGN #undef JITTER } #ifdef mRGBE_TEST /* Build standalone to run unit test */ #define RGBMIN 1 #define RGBMAX 15 int main (int argc, char **argv) { double rgb [3] = {0, 0, 0}, rgbDec [3], rgbMin [3] = {0, 0, 0}, rgbMax [3] = {1, 1, 1}, rgbAbs, relDiff, rmse, rmseSum; unsigned numTests, i, j, k, data, dataDec; const short rndState0 [3] = {11, 23, 31}; short encFlag = 0, rndState [3]; mRGBERange mrgbeRange; mRGBE mrgbe; #if 0 rgb [0] = -1.481; rgb [1] = -0.437; rgb [2] = -1.392; rgbMin [0] = 1.003; rgbMin [1] = 0.436; rgbMin [2] = 0.299; rgbMax [0] = 14.863; rgbMax [1] = 13.505; rgbMax [2] = 13.359; mRGBEinit(&mrgbeRange, rgbMin, rgbMax); mrgbe = mRGBEencode(rgb, &mrgbeRange, 0); printf("% 7.4f\t% 7.4f\t% 7.4f\t%4d\t-->\t" "%02d:%02d:%02d:%02d\t-->\t", rgb [0], rgb [1], rgb [2], 0, mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp ); /* Decode */ dataDec = mRGBEdecode(mrgbe, &mrgbeRange, rgbDec); for (j = rmse = 0; j < 3; j++) { relDiff = 1 - rgbDec [j] / rgb [j]; rmse += relDiff * relDiff; } rmseSum += rmse = sqrt(rmse / 3); printf("% 7.4f\t% 7.4f\t% 7.4f\t%4d\tRMSE = %.1f%%\n", rgbDec [0], rgbDec [1], rgbDec [2], dataDec, 100 * rmse ); return 0; #endif if (argc > 1 && (numTests = atoi(argv [1])) > 0) { /* Test range checks */ mRGBEinit(&mrgbeRange, rgbMin, rgbMax); puts("Check zero"); rgb [0] = rgb [1] = rgb [2] = 0; mrgbe = mRGBEencode(rgb, &mrgbeRange, 0); mRGBEdecode(mrgbe, &mrgbeRange, rgbDec); printf("% 7.3f\t% 7.3f\t% 7.3f\t-->\t%02d:%02d:%02d\t-->" "\t% 7.3f\t% 7.3f\t%7.3f\n", rgb [0], rgb [1], rgb [2], mrgbe.red, mrgbe.grn, mrgbe.blu, rgbDec [0], rgbDec [1], rgbDec [2] ); puts("Check data overflow"); rgb [0] = rgb [1] = rgb [2] = mRGBE_MAX / 2; mRGBEencode(rgb, &mrgbeRange, mRGBE_DATAMAX << 1); puts("\nCheck underflow"); rgb [0] = rgb [1] = rgb [2] = mRGBE_MIN / 2; mRGBEencode(rgb, &mrgbeRange, 0); puts("\nCheck overflow (positive)"); rgb [0] = rgb [1] = rgb [2] = 2 * mRGBE_MAX; mRGBEencode(rgb, &mrgbeRange, 0); puts("\nCheck overflow (negative)"); rgb [0] = rgb [1] = rgb [2] = -2 * mRGBE_MAX; mRGBEencode(rgb, &mrgbeRange, 0); puts("\nCheck empty range"); rgbMax [0] = rgbMax [1] = rgbMax [2] = 0; mRGBEinit(&mrgbeRange, rgbMin, rgbMax); putchar('\n'); rgbMin [0] = rgbMin [1] = rgbMin [2] = RGBMAX; rgbMax [0] = rgbMax [1] = rgbMax [2] = RGBMIN; /* Iterate twice with same random RGB tuples; gather minimum and * maximum in 1st pass, then initialise and run encode/decode tests * in 2nd pass. */ do { /* Initialise RNG state so this can be restored to obtain the same RGB tuple sequence in both passes. This state needs to be kept separate from the RNG used for jittering in the encode/decode routines. */ for (i = 0; i < 3; i++) rndState [i] = rndState0 [i]; /* Random encode vs. decode */ for (k = 0, rmseSum = 0; k < numTests; k++) { data = (unsigned)(mRGBE_DATAMAX * erand48(rndState)); for (i = 0; i < 3; i++) { rgb [i] = (erand48(rndState) > 0.5 ? 1 : -1) * ((RGBMAX - RGBMIN) * erand48(rndState) + RGBMIN); #if 1 /* Correlate RGB */ if (i > 0) rgb [i] = rgb [0] * (0.1 + 0.9 * erand48(rndState)); #endif } if (encFlag) { printf("% 4d:\t", k); /* Encode */ mrgbe = mRGBEencode(rgb, &mrgbeRange, data); printf("% 7.3f\t% 7.3f\t% 7.3f\t%4d\t-->\t" "%02d:%02d:%02d:%02d\t-->\t", rgb [0], rgb [1], rgb [2], data, mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp ); /* Decode */ dataDec = mRGBEdecode(mrgbe, &mrgbeRange, rgbDec); for (j = rmse = 0; j < 3; j++) { relDiff = 1 - rgbDec [j] / rgb [j]; rmse += relDiff * relDiff; } rmseSum += rmse = sqrt(rmse / 3); printf("% 7.3f\t% 7.3f\t% 7.3f\t%4d\tRMSE = %.1f%% %s\n", rgbDec [0], rgbDec [1], rgbDec [2], dataDec, 100 * rmse, rmse > 0.1 ? "!" : "" ); #if 1 /* Stop at sign reversal; this error isn't acceptable * under any circumstances! */ for (j = 0; j < 3; j++) if (rgb [j] / rgbDec [j] < 0) { puts("**** SIGN REVERSAL!!! ****\n"); return -1; } #endif } else { /* Update RGB minimum and maximum */ for (i = 0; i < 3; i++) { rgbAbs = fabs(rgb [i]); if (rgbAbs > rgbMax [i]) rgbMax [i] = rgbAbs; if (rgbAbs < rgbMin [i]) rgbMin [i] = rgbAbs; #if 0 + /* NOTE: this increases RMSE if RGBMIN > 0 */ rgbMin [i] = 0; #endif } } } /* Dump results */ if (!encFlag) { /* Initialise mRGBE encoding state with min/max */ if (!mRGBEinit(&mrgbeRange, rgbMin, rgbMax)) return -1; printf( "RGB Range = [%.3f, %.3f, %.3f] - [%.3f, %.3f, %.3f]\n\n", rgbMin [0], rgbMin [1], rgbMin [2], rgbMax [0], rgbMax [1], rgbMax [2] ); } else printf("\nAvg RMSE = %.2f%%\n", 100 * rmseSum / numTests); } while (encFlag = !encFlag); return 0; } else { printf("Usage: %s \n", argv [0]); return -1; } } #endif diff --git a/mrgbe.h b/mrgbe.h index 35c2881..1c6ba38 100644 --- a/mrgbe.h +++ b/mrgbe.h @@ -1,98 +1,98 @@ /* ========================================================================= mRGBE (miniRGBE): reduced RGBE encoding of normalised RGB floats plus associated integer payload data. By default, this encoding consists of: - three signed 5-bit mantissas for each RGB colour channel, - a common 5-bit exponent, - and an associated 12-bit integer datum. However, this allocation can be modified within a 32-bit envelope if more precision is required by modifying mRGBE_MANTBITS, mRGBE_EXPBITS and mRGBE_DATABITS. The encoded float RGB tuples are assumed to lie in the range [-1, 1], with values in the interval [-mRGBE_MIN, mRGBE_MIN] clamped to 0, since these cannot be resolved. The encoded values are therefore normalised against a given normalisation factor (e.g. absolute maximum) in order to maximise the useful range. 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 _mRGBE_H #define _mRGBE_H #include /* Mantissa / exponent / payload data bits and their encoding ranges. NOTE: binary shifts can overflow if the number of bits is increased; in this case shifts must either be promoted to 64 bits (possibly at the expense of portability), or replaced by much slower pow(2, ..) */ #define mRGBE_MANTBITS 5 #define mRGBE_EXPBITS 5 #define mRGBE_DATABITS 12 /* Rounding amplitude for encoding NOTE: unit test indicates this actually _increases_ RMSE! */ #define mRGBE_ROUND 0 /* Jitter amplitude for decoding */ - #define mRGBE_JITTER 0.75 + #define mRGBE_JITTER 0.5 /* Signed maximum for mantissa */ #define mRGBE_MANTMAX (1L << (mRGBE_MANTBITS - 1)) /* Absolute overflow limit for mantissa (sign flips if exceeded) */ #define mRGBE_MANTOVF ((mRGBE_MANTMAX << 1) - 1) #define mRGBE_MIN (1.0 / (1L << (1 << mRGBE_EXPBITS))) #define mRGBE_MAX 1.0 #define mRGBE_RANGE (mRGBE_MAX - mRGBE_MIN) #define mRGBE_DATAMAX ((1U << mRGBE_DATABITS) - 1) typedef union { struct { /* (mini)RGBE representation consisting of a 4-bit mantissa + sign * per RGB, a common 5-bit exponent, and a 12-bit payload data */ unsigned red: mRGBE_MANTBITS, grn: mRGBE_MANTBITS, blu: mRGBE_MANTBITS, exp: mRGBE_EXPBITS, dat: mRGBE_DATABITS; }; uint32_t all; } mRGBE; /* State for mRGBE encoding to offset and normalise floats in order to * optimally utilise the encoding range */ typedef struct { double min [3], max [3], norm [3]; } mRGBERange; mRGBERange *mRGBEinit (mRGBERange *range, double rgbMin [3], double rgbMax [3] ); /* Initialise range state for mRGBEencode()/mRGBEdecode() with * the specified float input ranges [rgbMin, rgbMax] in _absolute_ values * per colour channel. This sets the offset and normalisation for the * encoding. */ mRGBE mRGBEencode (double rgb [3], const mRGBERange *range, unsigned data); /* Encode signed float RGB tuple to mRGBE along with payload data. The * supplied range state must be previously initialised with mRGBEinit(). * A warning is issued if rgb lies outside the initialised range. */ unsigned mRGBEdecode (mRGBE mrgbe, const mRGBERange *range, double rgb [3]); /* Decode mRGBE, returning signed float RGB tuple in array rgb. The * associated payload data is decoded as return value. The supplied * range state must be previously initialised with mRGBEinit(). */ #endif diff --git a/pmapcontrib.c b/pmapcontrib.c index 6a7490c..3104a0b 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1104 +1,1115 @@ #ifndef lint static const char RCSid[] = "$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines handle contribution binning, compression and encoding, 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: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $ */ /* TODO: Isolate sanity checks with PMAP_CONTRIB_DBG once stable! */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmutil.h" #include "otspecial.h" #include "otypes.h" #include "lookup.h" #ifdef PMAP_CONTRIB /* The following are convenient placeholders interfacing to mkpmap and rcontrib. They can be externally set via initPmapContribTab() and referenced within the contrib pmap modules. These variables can then be safely ignored by rtrace/rpict/rvu, without annoying linking errors. */ /* Global pointer to rcontrib's contribution binning LUT */ LUTAB *pmapContribTab = NULL; /* Contribution/coefficient mode flag */ int *pmapContribMode; extern void SDdisk2square(double sq[2], double diskx, double disky); static int logDim (unsigned size) /* Return log2(sqrt(size)) = l, where size = (2^l)(2^l). Is there a faster way (e.g. binary search)? */ { unsigned i, sz, dim; if (!size) return 0; for (i = 0, dim = 0, sz = size; sz >>= 2; dim++); return dim; } static int xy2bin (unsigned l, int x, int y) /* Serialise 2D coordinates in range (2^l) x (2^l) to 1D bin. Returns -1 if coordinates invalid */ { return x < 0 || y < 0 ? -1 : (x << l) + y; } static void bin2xy (unsigned l, int bin, int *x, int *y) /* Deserialise 1D bin to 2D coordinates in range (2^l) x (2^l). Returns -1 if bin invalid */ { /* Obviously this is faster than integer division/modulo */ *x = bin < 0 ? -1 : bin >> l; *y = bin < 0 ? -1 : bin & PMAP_CONTRIB_SCDIM(l) - 1; } static int ray2bin (const RAY *ray, unsigned nbins) /* Map ray dir (pointing away from origin) to its 1D Shirley-Chiu bin, where nbins = (2^l) x (2^l) for l > 1. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { const unsigned l = logDim(nbins), scDim = PMAP_CONTRIB_SCDIM(l); static int scRHS, varInit = 0; static FVECT scNorm, scUp; unsigned scBin [2]; FVECT diskPlane; RREAL dz, diskx, disky, rad, diskd2, scCoord [2]; if (!varInit) { /* Lazy init shirley-Chiu mapping orientation from function variables */ scRHS = varvalue(PMAP_CONTRIB_SCRHS); scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX); scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY); scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ); scUp [0] = varvalue(PMAP_CONTRIB_SCUPX); scUp [1] = varvalue(PMAP_CONTRIB_SCUPY); scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ); varInit ^= 1; } /* Map incident direction to disk */ dz = DOT(ray -> rdir, scNorm); /* normal proj */ if (dz > 0) { fcross(diskPlane, scUp, scNorm); /* y-axis in plane, perp to up */ diskPlane [0] *= scRHS; diskx = DOT(ray -> rdir, diskPlane); disky = DOT(ray -> rdir, scUp) - dz * DOT(ray -> rdir, scUp); diskd2 = diskx * diskx + disky * disky; /* in-plane length^2 */ rad = dz>FTINY ? sqrt((1 - dz*dz) / diskd2) : 0; /* radial factor */ disky *= rad; /* diskx, disky now in range [-1, 1] */ /* Apply Shirley-Chiu mapping of (diskx, disky) to square */ SDdisk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return xy2bin(l, scBin [0], scBin [1]); } else return -1; } /* ------------------ CONTRIBSRC STUFF --------------------- */ #ifndef PMAP_CONTRIB_TEST MODCONT *addContribModifier (LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, int binCnt ) /* Add light source modifier mod to contribution lookup table contribTab, and update numContribs. Return initialised contribution data for this modifier. This code adapted from rcontrib.c:addmodifier(). */ { LUENT *lutEntry = lu_find(contribTab, mod); MODCONT *contrib; EPNODE *eBinVal; if (lutEntry -> data) { /* Reject duplicate modifiers */ sprintf(errmsg, "duplicate light source modifier %s", mod); error(USER, errmsg); } if (*numContribs >= MAXMODLIST) { sprintf(errmsg, "too many modifiers (max. %d)", MAXMODLIST); error(INTERNAL, errmsg); } lutEntry -> key = mod; if (binCnt <= 0) { sprintf(errmsg, "undefined/invalid bin count for modifier %s", mod); error(USER, errmsg); } /* Allocate and init contributions */ contrib = (MODCONT *)malloc( sizeof(MODCONT) + sizeof(DCOLOR) * (binCnt - 1) ); if (!contrib) error(SYSTEM, "out of memory in addContribModifier()"); contrib -> outspec = NULL; contrib -> modname = mod; contrib -> params = binParm; contrib -> nbins = binCnt; contrib -> binv = eBinVal; contrib -> bin0 = 0; memset(contrib -> cbin, 0, sizeof(DCOLOR) * binCnt); lutEntry -> data = lutEntry -> data = (char *)contrib; ++(*numContribs); return(contrib); } void addContribModfile (LUTAB *contribTab, unsigned *numContribs, char *modFile, char *binParm, int binCnt ) /* Add light source modifiers from file modFile to contribution lookup table * contribTab, and update numContribs. * NOTE: This code is adapted from rcontrib.c */ { char *mods [MAXMODLIST]; int i; /* Find file and store strings */ i = wordfile(mods, MAXMODLIST, getpath(modFile, getrlibpath(), R_OK)); if (i < 0) { sprintf(errmsg, "can't open modifier file %s", modFile); error(SYSTEM, errmsg); } if (*numContribs + i >= MAXMODLIST - 1) { sprintf(errmsg, "too many modifiers (max. %d) in file %s", MAXMODLIST - 1, modFile ); error(INTERNAL, errmsg); } for (i = 0; mods [i]; i++) /* Add each modifier */ addContribModifier(contribTab, numContribs, mods [i], binParm, binCnt ); } static int contribSourceBin (LUTAB *contribs, const RAY *ray) /* Map contribution source ray to its bin for light source ray -> rsrc, using Shirley-Chiu disk-to-square mapping. Return invalid bin -1 if mapping failed. */ { const SRCREC *src; const OBJREC *srcMod; const MODCONT *srcCont; RAY srcRay; int bin, i; /* Check we have a valid ray and contribution LUT */ if (!ray || !contribs) return -1; src = &source [ray -> rsrc]; srcMod = findmaterial(src -> so); srcCont = (MODCONT*)lu_find(contribs, srcMod -> oname) -> data; if (!srcCont) /* Not interested in this source (modifier not in contrib LUT) */ return -1; /* Set up shadow ray pointing to source for disk2square mapping */ rayorigin(&srcRay, SHADOW, NULL, NULL); srcRay.rsrc = ray -> rsrc; VCOPY(srcRay.rorg, ray -> rop); for (i = 0; i < 3; i++) srcRay.rdir [i] = -ray -> rdir [i]; if (!(src -> sflags & SDISTANT ? sourcehit(&srcRay) : (*ofun[srcMod -> otype].funp)(srcMod, &srcRay) )) /* (Redundant?) sanity check for valid source ray? */ return -1; #if 0 worldfunc(RCCONTEXT, &srcRay); #endif set_eparams((char*)srcCont -> params); if ((bin = ray2bin(&srcRay, srcCont -> nbins)) < 0) error(WARNING, "Ignoring invalid bin in contribSourceBin()"); return bin; } PhotonContribSourceIdx newPhotonContribSource (PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap ) /* Add contribution source ray for emitted contribution photon and save * light source index and binned direction. The current contribution * source is stored in pmap -> lastContribSrc. If the previous contrib * source spawned photons (i.e. has srcIdx >= 0), it's appended to * contribSrcHeap. If contribSrcRay == NULL, the current contribution * source is still flushed, but no new source is set. Returns updated * contribution source counter pmap -> numContribSrc. */ { if (!pmap || !contribSrcHeap) return 0; /* Check if last contribution source has spawned photons (srcIdx >= 0, * see newPhoton()), in which case we save it to the heap file before * clobbering it. (Note this is short-term storage, so we doan' need * da portable I/O stuff here). */ if (pmap -> lastContribSrc.srcIdx >= 0) { if (!fwrite(&pmap -> lastContribSrc, sizeof(PhotonContribSource), 1, contribSrcHeap )) error(SYSTEM, "failed writing photon contrib source in " "newPhotonContribSource()" ); pmap -> numContribSrc++; if (!pmap -> numContribSrc || pmap -> numContribSrc > PMAP_MAXCONTRIBSRC ); } if (contribSrcRay) { /* Mark this contribution source unused with a negative source index until its path spawns a photon (see newPhoton()) */ pmap -> lastContribSrc.srcIdx = -1; /* Map ray to bin in anticipation that this contrib source will be used, since the ray will be lost once a photon is spawned */ pmap -> lastContribSrc.srcBin = contribSourceBin( pmapContribTab, contribSrcRay ); if (pmap -> lastContribSrc.srcBin < 0) { /* Warn if invalid bin, but trace photon nonetheless. It will count as emitted to prevent bias, but will not be stored in newPhoton(), as it contributes zero flux */ sprintf(errmsg, "invalid bin for light source %s, " "contribution photons will be discarded", source [contribSrcRay -> rsrc].so -> oname ); error(WARNING, errmsg); } } return pmap -> numContribSrc; } PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ) /* Consolidate per-subprocess contribution source heaps into array * pmap -> contribSrc. Returns offset for contribution source index * linearisation in pmap -> numContribSrc. The heap files in * contribSrcHeap are closed on return. */ { PhotonContribSourceIdx heapLen; unsigned heap; if (!pmap || !contribSrcHeap || !contribSrcOfs || !numHeaps) return 0; pmap -> numContribSrc = 0; for (heap = 0; heap < numHeaps; heap++) { contribSrcOfs [heap] = pmap -> numContribSrc; if (fseek(contribSrcHeap [heap], 0, SEEK_END) < 0) error(SYSTEM, "failed photon contrib source seek " "in buildContribSources()" ); pmap -> numContribSrc += heapLen = ftell(contribSrcHeap [heap]) / sizeof(PhotonContribSource); if (!(pmap -> contribSrc = realloc(pmap -> contribSrc, pmap -> numContribSrc * sizeof(PhotonContribSource) ))) error(SYSTEM, "failed photon contrib source alloc " "in buildContribSources()" ); rewind(contribSrcHeap [heap]); if (fread(pmap -> contribSrc + contribSrcOfs [heap], sizeof(PhotonContribSource), heapLen, contribSrcHeap [heap] ) != heapLen) error(SYSTEM, "failed reading photon contrib source " "in buildContribSources()" ); fclose(contribSrcHeap [heap]); unlink(contribSrcHeapFname [heap]); } return pmap -> numContribSrc; } /* ----------------- CONTRIB PMAP ENCODING STUFF -------------------- */ typedef struct { WAVELET_COEFF *coeff; unsigned bin; } EncodedPreComputedContrib; static int coeffCompare (const void *c1, const void *c2) /* Comparison function to REVERSE sort thresholded coefficients */ { const EncodedPreComputedContrib *tcoeff1 = c1, *tcoeff2 = c2; /* Use dot product as magnitude to compare _absolute_ values */ const WAVELET_COEFF v1 = DOT(tcoeff1 -> coeff, tcoeff1 -> coeff), v2 = DOT(tcoeff2 -> coeff, tcoeff2 -> coeff); if (v1 < v2) return 1; else if (v1 > v2) return -1; else return 0; } static int thresholdContribs (MODCONT *binnedContribs, PreComputedContrib *preCompContrib ) /* Threshold binned detail coefficients in binnedContribs -> cbin [1..] by keeping the (preCompContrib -> nCompressedBins) largest and returning these in preContrib -> compressedBins along with their original bin order. NOTE: binnedContribs -> cbin [0] is the average wavelet coefficient and excluded from thresholding. Returns 0 on success. */ { unsigned b; EncodedPreComputedContrib *threshCoeffs; if (!binnedContribs || !preCompContrib || !preCompContrib -> compressedBins ) /* Should be initialised by caller! */ return -1; /* Set up coefficient buffer for compression (thresholding), skipping * first coeff since it's the average and therefore not thresholded */ threshCoeffs = (EncodedPreComputedContrib*)( preCompContrib -> compressedBins ); for (b = 0; b < binnedContribs -> nbins - 1; b++) { /* Set up pointers to coeffs (sorted faster than 3 floats/doubles) and remember original bin position prior to sorting. NOTE: The encoded bin position preserves the offset for the 1st coefficient i.e. the average! */ threshCoeffs [b].coeff = (WAVELET_COEFF*)( &(binnedContribs -> cbin [b + 1]) ); threshCoeffs [b].bin = b + 1; } /* REVERSE sort coeffs by magnitude; the non-thresholded coeffs will then start at threshCoeffs [0] */ qsort(threshCoeffs, binnedContribs -> nbins - 1, sizeof(EncodedPreComputedContrib), coeffCompare ); return 0; } static int encodeContribs (MODCONT *binnedContribs, PreComputedContrib *preCompContrib, float compressRatio ) /* Apply wavelet transform to binnedContribs -> cbin and compress according to compressRatio, storing thresholded and mRGBE-encoded coefficients in preCompContrib -> mrgbeCoeffs. Note that the average coefficient is not encoded, and returned in binnedContribs -> cbin [0]. Returns 0 on success. */ { unsigned b, i; EncodedPreComputedContrib *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; WAVELET_COEFF absCoeff; #ifdef PMAP_CONTRIB_DBG WaveletCoeff3 decCoeff; unsigned decBin; #endif if (!binnedContribs || !preCompContrib || !preCompContrib -> compressedBins || !preCompContrib -> mrgbeCoeffs ) /* Should be initialised by the caller! */ return -1; /* Do 2D wavelet transform on preCompContrib -> waveletMatrix (which actually maps to binnedContribs -> cbin) */ if (waveletXform2(preCompContrib -> waveletMatrix, preCompContrib -> tWaveletMatrix, preCompContrib -> l ) < 0 ) error(INTERNAL, "failed wavelet transform in encodeContribs()"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(binnedContribs, preCompContrib) < 0) error(INTERNAL, "failed wavelet compression in encodeContribs()"); threshCoeffs = (EncodedPreComputedContrib*)( preCompContrib -> compressedBins ); /* Init per-channel coefficient range for mRGBE encoding */ mrgbeRange = &preCompContrib -> mrgbeRange; setcolor(mrgbeRange -> min, FHUGE, FHUGE, FHUGE); setcolor(mrgbeRange -> max, 0, 0, 0); /* Update per-channel coefficient range */ for (b = 0; b < preCompContrib -> nCompressedBins; b++) { for (i = 0; i < 3; i++) { #ifdef PMAP_CONTRIB_DBG + #if 0 /* Encode bin as coeff instead of actual coeff itself */ threshCoeffs [b].coeff [i] = b + 1; + #endif #endif absCoeff = fabs(threshCoeffs [b].coeff [i]); if (absCoeff < mrgbeRange -> min [i]) mrgbeRange -> min [i] = absCoeff; if (absCoeff > mrgbeRange -> max [i]) mrgbeRange -> max [i] = absCoeff; } } if (preCompContrib -> nCompressedBins == 1) /* Maximum compression with just 1 (!) compressed bin (is this even * useful?), so mRGBE range is undefined since min & max coincide; * set minimum to 0, maximum to the single remaining coefficient */ setcolor(mrgbeRange -> min, 0, 0, 0); /* Init mRGBE coefficient normalisation from range */ mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max); mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; /* Encode wavelet detail coefficients to mRGBE */ for (b = 0; b < preCompContrib -> nCompressedBins; b++) { mrgbeCoeffs [b] = mRGBEencode(threshCoeffs [b].coeff, mrgbeRange, threshCoeffs [b].bin ); #ifdef PMAP_CONTRIB_DBG /* Encoding sanity check */ decBin = mRGBEdecode(mrgbeCoeffs [b], mrgbeRange, decCoeff); + #if 0 if (decBin != threshCoeffs [b].bin || - dist2(decCoeff, threshCoeffs [b].coeff) > 1e-3 + sqrt(dist2(decCoeff, threshCoeffs [b].coeff)) > 1.7 ) error(CONSISTENCY, "failed sanity check in encodeContribs()"); - + #endif for (i = 0; i < 3; i++) if (decCoeff [i] / threshCoeffs [b].coeff [i] < 0) error(CONSISTENCY, "coefficient sign reversal in encodeContribs()" ); #endif } return 0; } static void initContribHeap (PhotonMap *pmap) /* Initialise precomputed contribution heap */ { int fdFlags; if (!pmap -> contribHeap) { /* Open heap file for precomputed contributions */ mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); if (!(pmap -> contribHeap = fopen(pmap -> contribHeapFname, "w+b"))) error(SYSTEM, "failed opening precomputed contribution " "heap file in initContribHeap()" ); #ifdef F_SETFL /* XXX is there an alternative needed for Wind0z? */ fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL); fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> heap), 0); */ } } static MODCONT *getPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Locate photons near ray -> rop which originated from the light source modifier indexed by ray -> rsrc, and accumulate their contributions in pmap -> srcContrib. Return total contributions in irrad and pointer to binned contributions (or NULL if none were found). */ { const OBJREC *srcMod = findmaterial(source [ray->rsrc].so); MODCONT *contrib = (MODCONT*)lu_find( pmapContribTab, srcMod -> oname ) -> data; Photon *photon; COLOR flux; PhotonSearchQueueNode *sqn; float r2, norm; unsigned i; /* Zero bins for this source modifier */ memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins); setcolor(irrad, 0, 0, 0); if (!contrib || !pmap -> maxGather) /* Modifier not in LUT or zero bandwidth */ return NULL; /* Lookup photons */ pmap -> squeue.tail = 0; /* Pass light source index to filter in findPhotons() via pmap -> lastContribSrc */ pmap -> lastContribSrc.srcIdx = ray -> rsrc; findPhotons(pmap, ray); /* Need at least 2 photons */ if (pmap -> squeue.tail < 2) { #ifdef PMAP_NONEFOUND sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", ray -> ro ? ray -> ro -> oname : "", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); error(WARNING, errmsg); #endif return NULL; } /* Avg radius^2 between furthest two photons to improve accuracy */ sqn = pmap -> squeue.node + 1; r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); /* XXX: OMIT extra normalisation factor 1/PI for ambient calc? */ #ifdef PMAP_EPANECHNIKOV /* Normalise accumulated flux by Epanechnikov kernel integral in 2D (see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986, and Wann Jensen, "Realistic Image Synthesis using Photon Mapping"), include RADIANCE-specific lambertian factor PI */ norm = 2 / (PI * PI * r2); #else /* Normalise accumulated flux by search area PI * r^2, including RADIANCE-specific lambertian factor PI */ norm = 1 / (PI * PI * r2); #endif /* Skip the extra photon */ for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); scalecolor(flux, norm); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on photon distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); + #ifdef PMAP_CONTRIB_DBG + /* Set photon flux to bin for debugging */ + setcolor(contrib -> cbin [photonSrcBin(pmap, photon)], + photonSrcBin(pmap, photon), + photonSrcBin(pmap, photon), + photonSrcBin(pmap, photon) + ); + #endif addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } void freePreCompContribNode (void *p) /* Clean up precomputed contribution LUT entry */ { PhotonMap *preCompContribPmap = (PhotonMap*)p; PreComputedContrib *preCompContrib; if (preCompContribPmap) { preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); if (preCompContrib) { /* Free primary and transposed wavelet matrices */ free(preCompContrib -> waveletMatrix); freeWaveletMatrix(preCompContrib -> tWaveletMatrix, preCompContrib -> l ); /* Free thresholded coefficients */ free(preCompContrib -> compressedBins); /* Free mRGBE encoded coefficients */ free(preCompContrib -> mrgbeCoeffs); free(preCompContrib -> waveletFname); } /* Clean up precomputed contrib photon map */ deletePhotons(preCompContribPmap); free(preCompContribPmap); } } void initPreComputedContrib (PreComputedContrib *preCompContrib) /* Initialise precomputed contribution container in photon map */ { preCompContrib -> waveletFname = NULL; preCompContrib -> waveletFile = NULL; preCompContrib -> waveletMatrix = preCompContrib -> tWaveletMatrix = NULL; preCompContrib -> compressedBins = NULL; setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0); setcolor(preCompContrib -> mrgbeRange.max, 0, 0, 0); setcolor(preCompContrib -> mrgbeRange.norm, 0, 0, 0); preCompContrib -> mrgbeCoeffs = NULL; preCompContrib -> nBins = preCompContrib -> nCompressedBins = preCompContrib -> l = preCompContrib -> scDim = preCompContrib -> contribSize = 0; } void preComputeContrib (PhotonMap *pmap) /* Precompute contributions for a random subset of (finalGather * pmap -> numPhotons) photons, and return the per-modifier precomputed contribution photon maps in LUT preCompContribTab, discarding the original photons. */ { unsigned long i, numPreComp; unsigned b, j; long nCompressedBins; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; COLR mrgbeRange32 [2]; LUENT *preCompContribNode; PhotonMap *preCompContribPmap, nuPmap; PreComputedContrib *preCompContrib; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); #ifdef PMAP_CONTRIB_DBG RAY dbgRay; #endif if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Init new parent photon map and set output filename */ initPhotonMap(&nuPmap, pmap -> type); nuPmap.fileName = pmap -> fileName; /* Set contrib/coeff mode in new parent photon map */ nuPmap.contribMode = pmap -> contribMode; /* Allocate and init LUT containing per-modifier child photon maps */ if (!(nuPmap.preCompContribTab = malloc(sizeof(LUTAB)))) error(SYSTEM, "out of memory allocating LUT in preComputeContrib()" ); memcpy(nuPmap.preCompContribTab, &lutInit, sizeof(LUTAB)); /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repComplete = numPreComp = finalGather * pmap -> numPhotons; repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)( (i + pmapRandom(pmap -> randState)) / finalGather ); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection and normal at photon position. Set emitting light source index from origin. */ VCOPY(ray.rop, photon.pos); ray.rsrc = photonSrcIdx(pmap, &photon); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; #ifndef PMAP_CONTRIB_DBG /* Get contributions at photon position; retry with another photon if no contribs found */ binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); #else /* Get all contribs from same photon for debugging */ getPhoton(pmap, 0, &photon); memcpy(&dbgRay, &ray, sizeof(RAY)); VCOPY(dbgRay.rop, photon.pos); dbgRay.rsrc = photonSrcIdx(pmap, &photon); for (j = 0; j < 3; j++) dbgRay.ron [j] = photon.norm [j] / 127.0; binnedContribs = getPhotonContrib(pmap, &dbgRay, dbgRay.rcol); #endif if (binnedContribs) { /* Found contributions at photon position, so generate * precomputed photon */ if (!(preCompContribNode = lu_find(nuPmap.preCompContribTab, binnedContribs -> modname) ) ) error(SYSTEM, "out of memory allocating LUT entry in " "preComputeContrib()" ); if (!preCompContribNode -> key) { /* New LUT node for precomputed contribs for this modifier */ preCompContribNode -> key = (char*)binnedContribs -> modname; preCompContribNode -> data = (char*)( preCompContribPmap = malloc(sizeof(PhotonMap)) ); if (preCompContribPmap) { /* Init new child photon map and its contributions */ initPhotonMap(preCompContribPmap, PMAP_TYPE_CONTRIB_CHILD ); initPhotonHeap(preCompContribPmap); initContribHeap(preCompContribPmap); preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib = malloc(sizeof(PreComputedContrib)) ); initPreComputedContrib(preCompContrib); } if (!preCompContribPmap || !preCompContrib) error(SYSTEM, "out of memory allocating new photon map " "in preComputeContrib()" ); /* Set output filename from parent photon map */ preCompContribPmap -> fileName = nuPmap.fileName; /* Set contrib/coeff mode from parent photon map */ preCompContribPmap -> contribMode = nuPmap.contribMode; /* Get Shirley-Chiu square & wavelet matrix resolution, bins per dimension, and number of bins */ preCompContrib -> l = logDim(binnedContribs -> nbins); preCompContrib -> scDim = PMAP_CONTRIB_SCDIM( preCompContrib -> l ); preCompContrib -> nBins = preCompContrib -> scDim * preCompContrib -> scDim; if (preCompContrib -> nBins != binnedContribs -> nbins) /* Shouldn't happen */ error(INTERNAL, "bin count mismatch in preComputeContrib()" ); if (preCompContrib -> nBins > 1) { /* BINNING ENABLED; set up wavelet xform & compression. Number of compressed coeffs/bins is fixed (subtract one as the average coefficient is not thresholded) */ nCompressedBins = (preCompContrib -> nBins - 1) * (1 - compressRatio); if (nCompressedBins < 1) { error(WARNING, "maximum contribution compression, clamping number " "of bins in preComputeContrib()" ); nCompressedBins = 1; } if (nCompressedBins >= preCompContrib -> nBins) { error(WARNING, "minimum contribution compression, clamping number " "of bins in preComputeContrib()" ); nCompressedBins = preCompContrib -> nBins - 1; } preCompContrib -> nCompressedBins = nCompressedBins; /* Lazily allocate 2D primary wavelet matrix and map its rows to linear array binnedContribs -> cbin */ if (!(preCompContrib -> waveletMatrix = calloc( preCompContrib -> scDim, sizeof(WaveletCoeff3*)) ) ) error(SYSTEM, "out of memory allocating primary " "wavelet matrix in preComputeContrib()" ); for (b = 0; b < preCompContrib -> scDim; b++) /* Point to each row in existing 1D contrib array */ preCompContrib -> waveletMatrix [b] = &binnedContribs -> cbin [b * preCompContrib -> scDim]; /* Lazily allocate transposed wavelet matrix */ preCompContrib -> tWaveletMatrix = allocWaveletMatrix(preCompContrib -> l); if (!preCompContrib -> tWaveletMatrix) error(SYSTEM, "out of memory allocating transposed wavelet " "matrix in preComputeContrib()" ); /* Lazily allocate thresholded detail coefficients (minus the average coeff) */ preCompContrib -> compressedBins = calloc( preCompContrib -> nBins - 1, sizeof(EncodedPreComputedContrib) ); /* Lazily alloc mRGBE-encoded compressed wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc( preCompContrib -> nCompressedBins, sizeof(mRGBE) ); if (!preCompContrib -> compressedBins || !preCompContrib -> mrgbeCoeffs ) error(SYSTEM, "out of memory allocating compressed " "coefficients in preComputeContrib()" ); /* Set size of compressed contributions, in bytes */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( preCompContrib -> nCompressedBins ); } } else { /* LUT node already exists for this modifier; use its data */ preCompContribPmap = (PhotonMap*)preCompContribNode -> data; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); } if (binnedContribs -> nbins > 1) { /* Binning enabled; compress & encode binned contribs */ if (encodeContribs(binnedContribs, preCompContrib, compressRatio ) ) error(INTERNAL, "failed contribution " "compression/encoding in preComputeContrib()" ); /* Dump encoded bins to contrib heap file, prepended by mRGBE range in 32-bit RGBE. NOTE: minimum and maximum are reversed in the file to facilitate handling of single compressed bins (see below) */ setcolr(mrgbeRange32 [0], preCompContrib -> mrgbeRange.max [0], preCompContrib -> mrgbeRange.max [1], preCompContrib -> mrgbeRange.max [2] ); setcolr(mrgbeRange32 [1], preCompContrib -> mrgbeRange.min [0], preCompContrib -> mrgbeRange.min [1], preCompContrib -> mrgbeRange.min [2] ); /* Dump only mRGBE maximum if 1 compressed bin (maximum * compression), since minimum is implicitly zero and can * be omitted to save space */ putbinary(mrgbeRange32, sizeof(COLR), 1 + (preCompContrib -> nCompressedBins > 1), preCompContribPmap -> contribHeap ); if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), preCompContrib -> nCompressedBins, preCompContribPmap -> contribHeap ) != preCompContrib -> nCompressedBins ) error(SYSTEM, "failed writing to coefficients to " "contribution heap in preComputeContrib()" ); } /* Set photon flux to coarsest average wavelet coefficient NOTE: binnedContrib -> cbin [0] == preCompContrib -> waveletMatrix [0][0] == preCompContrib -> tWaveletMatrix [0][0]. IF BINNING IS DISABLED, this is the total contribution from all directions */ copycolor(ray.rcol, binnedContribs -> cbin [0]); /* HACK: signal newPhoton() to set precomputed photon's contribution source from ray -> rsrc */ /* XXX: DO WE STILL NEED THIS SHIT AFTER PRECOMP, SINCE CHILD PMAPS ARE PER-MODIFIER ANYWAY? */ preCompContribPmap -> lastContribSrc.srcIdx = -2; /* Append photon to new heap from ray and increment total count for all modifiers in parent photon map */ newPhoton(preCompContribPmap, &ray); nuPmap.numPhotons++; } /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Trash original pmap and its leaf file, and replace with new parent, which now acts as container for the per-modifier child pmaps */ unlink(pmap -> store.leafFname); deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Build per-modifier precomputed photon maps from their contribution heaps */ lu_doall(pmap -> preCompContribTab, buildPreCompContribPmap, NULL); } #else /* ------------------------------------------------------------------- U N I T T E S T S ------------------------------------------------------------------- */ #include #include #include "func.h" int main (int argc, char *argv []) { unsigned i, l, nbins, nsamp, numTheta = 0, numPhi = 0; double t, p; RAY ray; int bin; if (argc < 3) { fprintf(stderr, "%s [=; ..]\n", argv [0]); fprintf(stderr, "Default variable defs: %s\n", PMAP_CONTRIB_SCDEFAULTS ); fputs("\nMissing resolution l>1, number of samples nsamp\n", stderr ); return -1; } if (!(l = atoi(argv [1]))) { fputs("Invalid resolution\n", stderr); return -1; } else nbins = PMAP_CONTRIB_SCBINS(l); if (!(nsamp = atoi(argv [2]))) { fputs("Invalid num samples\n", stderr); return -1; } else { numTheta = (int)(sqrt(nsamp) / 2); numPhi = 4* numTheta; } /* (Doan' need to) Init cal func routines for binning */ #if 0 initfunc(); setcontext(RCCONTEXT); #endif /* Compile default orientation variables for contribution binning */ scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); /* Compile custom orientation variabls from command line */ for (i = 3; i < argc; i++) scompile(argv [i], NULL, 0); for (i = 0; i < nsamp; i++) { #if 0 /* Random */ t = 0.5 * PI * drand48(); p = 2 * PI * drand48(); #else /* Stratified */ t = 0.5 * PI * ((i % numTheta) + drand48()) / numTheta; p = 2.0 * PI * ((i / numTheta) + drand48()) / numPhi; #endif ray.rdir [0] = sin(t) * cos(p); ray.rdir [1] = sin(t) * sin(p); ray.rdir [2] = cos(t); bin = ray2bin(&ray, nbins); printf("%.3f\t%.3f\t%.3f\t-->\t%d\n", ray.rdir [0], ray.rdir [1], ray.rdir [2], bin ); } return 0; } #endif #endif /* PMAP_CONTRIB */ diff --git a/wavelet2.c b/wavelet2.c index 8830359..5726691 100644 --- a/wavelet2.c +++ b/wavelet2.c @@ -1,668 +1,727 @@ /* ========================================================================= 2-dimensional wavelet transform on (2^l) x (2^l) sized arrays of 3-tuples, where l > 1. Compile with -DWAVELET_TEST to build standalone unit tests. 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$ */ -/* TODO: threshold in unit test */ - - - #include "wavelet2.h" #include #include #include #include -#define min(a, b) ((a) < (b) ? (a) : (b)) -#define coeffAvg(a) (((a) [0] + (a) [1] + (a) [2]) / 3) +#define min(a, b) ((a) < (b) ? (a) : (b)) +#define coeffAvg(a) (((a) [0] + (a) [1] + (a) [2]) / 3) +#define coeffThresh(c, t) ( \ + fabs((c) [0]) < (t) || fabs((c) [1]) < (t) || fabs((c) [2]) < (t) \ +) + /* The following defs are const, but strict compilers pretend to be dumber than they are, and refuse to init from a func (in this case sqrt(2) and sqrt(3)), or other consts */ #define SQRT2 1.41421356 #define SQRT3 1.73205081 #define H4NORM (0.25 / SQRT2) /* Haar wavelet coeffs */ static const WAVELET_COEFF h2 = 1 / SQRT2; /* Daubechies D4 wavelet coeffs */ static const WAVELET_COEFF h4 [4] = { H4NORM * (1 + SQRT3), H4NORM * (3 + SQRT3), H4NORM * (3 - SQRT3), H4NORM * (1 - SQRT3) }; static const WAVELET_COEFF g4 [4] = { H4NORM * (1 - SQRT3), -H4NORM * (3 - SQRT3), H4NORM * (3 + SQRT3), -H4NORM * (1 + SQRT3) }; /* g4 [4] = { h4 [3], -h4 [2], h4 [1], -h4 [0] }; */ WaveletMatrix allocWaveletMatrix (unsigned l) /* Allocate and return a 2D coefficient array of size (2^l) x (2^l), where l >= 1. Returns NULL if allocation failed. */ { const unsigned len = 1 << l; unsigned i; WaveletMatrix y = NULL; if (l >= 1) { if (!(y = calloc(len, sizeof(WaveletCoeff3*)))) return NULL; for (i = 0; i < len; i++) if (!(y [i] = calloc(len, WAVELET_COEFFSIZE))) return NULL; } return y; } void freeWaveletMatrix (WaveletMatrix y, unsigned l) /* Free previously allocated 2D coefficient array y of size (2^l) x (2^l) */ { unsigned i, j; const unsigned len = 1 << l; if (y) { for (i = 0; i < len; i++) free(y [i]); free(y); } } #ifdef WAVELET_DBG static void zeroCoeffs (WaveletMatrix y, unsigned l) /* Zero output array to facilitate debugging */ { const unsigned len = 1 << l; unsigned i; for (i = 0; i < len; i++) memset(y [i], 0, len * WAVELET_COEFFSIZE); } static void dumpCoeffs (const WaveletMatrix y, const WaveletMatrix yt, - unsigned l + unsigned l, float thresh ) /* Dump arrays y and yt side-by-side to stdout (skip yt if NULL) */ { const unsigned len = 1 << l; unsigned i, j, k; for (i = 0; i < len; i++) { for (j = 0; j < len; j++) - printf("%4.3f\t", coeffAvg(y [i][j])); + printf( + coeffThresh(y [i][j], thresh) && (i | j) + ? "[% 3.2f]\t" : " % 3.2f \t", + coeffAvg(y [i][j]) + ); if (yt) { - printf("-->\t"); + printf(" -->\t"); for (j = 0; j < len; j++) - printf("%4.3f\t", coeffAvg(yt [i][j])); + printf("% 4.3f\t", coeffAvg(yt [i][j])); } putchar('\n'); } } + + + + static float rmseCoeffs (const WaveletMatrix y1, const WaveletMatrix y2, + unsigned l + ) + /* Calculate RMSE between matrices y and y0 */ + { + const unsigned len = 1 << l; + unsigned i, j; + float d, rmse = 0; + + for (i = 0; i < len; i++) + for (j = 0; j < len; j++) { + d = (coeffAvg(y1 [i][j]) - coeffAvg(y2 [i][j])) / + coeffAvg(y1 [i][j]); + rmse += d * d; + } + + return sqrt(rmse / ((float)len * len)); + } #endif static int haarStep (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of forward 2D Haar wavelet transform on array y of size (2^l) x (2^l) containing 3-tuples, where l >= 1. The transform is performed over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). Note the triplets per array element are _not_ decorrelated, but transformed independently. The wavelet coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the generated coefficients vertically, and prepares the array for another transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been decorrelated. Returns 0 on success, else -1. */ { static unsigned axis = 0; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 2 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: yt is transposed on the fly such that the next function call * transforms over the alternate axis. This is done by simply swapping * the indices during assignment */ for (i = 0; i < len; i++) { for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] = h2 * ( y [i] [j ] [k] + y [i] [j + 1] [k] ); /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = h2 * ( y [i] [j ] [k] - y [i] [j + 1] [k] ); } } } #ifdef WAVELET_DBG printf("%s FWD HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); - dumpCoeffs(y, yt, l); + dumpCoeffs(y, yt, l, 0); putchar('\n'); axis ^= 1; #endif return 0; } static int haarInvStep (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of inverse 2D Haar wavelet transform on coefficient array y of size (2^l) x (2^l) containing 3-tuples, where l >= 1. This reverses the forward transform above. The transform is inverted over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). The inverted coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the inverted coefficients vertically, and prepares the array for another inverse transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent inverse transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been inversely transformed. Returns 0 on success, else -1. */ { static unsigned axis = 1; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 2 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: i, j are swapped relative to the forward transform, as axis * order is now reversed. */ /* NOTE: yt is transposed on the fly such that the next function call * inverts over the alternate axis. This is done by simply swapping * the indices during assignment */ for (i = 0; i < len; i++) { for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { yt [i] [j ] [k] = h2 * ( y [h ] [i] [k] + /* Avg */ y [hlen + h] [i] [k] /* Diff */ ); yt [i] [j + 1] [k] = h2 * ( y [h ] [i] [k] - /* Avg */ y [hlen + h] [i] [k] /* Diff */ ); } } } #ifdef WAVELET_DBG printf("%s INV HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); - dumpCoeffs(y, yt, l); + dumpCoeffs(y, yt, l, 0); putchar('\n'); axis ^= 1; #endif return 0; } static int d4Step (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of forward 2D Daubechies D4 wavelet transform on array y of size (2^l) x (2^l) containing 3-tuples, where l >= 2. The transform is performed over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). Note the triplets per array element are _not_ decorrelated, but transformed independently. The wavelet coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the generated coefficients vertically, and prepares the array for another transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been decorrelated. Returns 0 on success, else -1. */ { static unsigned axis = 0; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 4 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: yt is transposed on the fly such that the next function call * transforms over the alternate axis. This is done by simply swapping * the indices during assignment */ for (i = 0; i < len; i++) { /* Transform until upper boundary */ for (j = 0; j < len - 2; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] = h4 [0] * y [i] [j ] [k] + h4 [1] * y [i] [j + 1] [k] + h4 [2] * y [i] [j + 2] [k] + h4 [3] * y [i] [j + 3] [k]; /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = g4 [0] * y [i] [j ] [k] + g4 [1] * y [i] [j + 1] [k] + g4 [2] * y [i] [j + 2] [k] + g4 [3] * y [i] [j + 3] [k]; } } /* Transform at upper boundary with wraparound. Note j is set to last index from previous loop */ h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] = h4 [0] * y [i] [j ] [k] + h4 [1] * y [i] [j + 1] [k] + h4 [2] * y [i] [0 ] [k] + h4 [3] * y [i] [1 ] [k]; /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = g4 [0] * y [i] [j ] [k] + g4 [1] * y [i] [j + 1] [k] + g4 [2] * y [i] [0 ] [k] + g4 [3] * y [i] [1 ] [k]; } } #ifdef WAVELET_DBG printf("%s FWD D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); - dumpCoeffs(y, yt, l); + dumpCoeffs(y, yt, l, 0); putchar('\n'); axis ^= 1; #endif return 0; } static int d4InvStep (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of inverse 2D Daubechies D4 wavelet transform on coefficient array y of size (2^l) x (2^l) containing 3-tuples, where l >= 2. This reverses the forward transform above. The transform is inverted over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). The inverted coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the inverted coefficients vertically, and prepares the array for another inverse transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent inverse transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been inversely transformed. Returns 0 on success, else -1. */ { static unsigned axis = 1; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 4 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: i, j are swapped relative to the forward transform, as axis * order is now reversed. */ /* NOTE: yt is transposed on the fly such that the next function call * inverts over the alternate axis. This is done by simply swapping * the indices during assignment */ for (i = 0; i < len; i++) { /* Invert at lower boundary with wraparound */ for (k = 0; k < 3; k++) { yt [i] [0] [k] = h4 [2] * y [hlen - 1] [i] [k] + /* Last avg */ g4 [2] * y [len - 1] [i] [k] + /* Last diff */ h4 [0] * y [0 ] [i] [k] + /* First avg */ g4 [0] * y [hlen ] [i] [k]; /* First diff */ yt [i] [1] [k] = h4 [3] * y [hlen - 1] [i] [k] + g4 [3] * y [len - 1] [i] [k] + h4 [1] * y [0 ] [i] [k] + g4 [1] * y [hlen ] [i] [k]; } /* Invert until upper boundary */ for (j = 2; j < len; j += 2) { h = (j >> 1) - 1; for (k = 0; k < 3; k++) { yt [i] [j ] [k] = h4 [2] * y [h ] [i] [k] + /* Avg */ g4 [2] * y [hlen + h ] [i] [k] + /* Diff */ h4 [0] * y [h + 1 ] [i] [k] + /* Next avg */ g4 [0] * y [hlen + h + 1] [i] [k]; /* Next diff */ yt [i] [j + 1] [k] = h4 [3] * y [h ] [i] [k] + g4 [3] * y [hlen + h ] [i] [k] + h4 [1] * y [h + 1 ] [i] [k] + g4 [1] * y [hlen + h + 1] [i] [k]; } } } #ifdef WAVELET_DBG printf("%s INV D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); - dumpCoeffs(y, yt, l); + dumpCoeffs(y, yt, l, 0); putchar('\n'); axis ^= 1; #endif return 0; } int waveletXform2 (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Perform full 2D multiresolution forward wavelet transform on array y of size (2^l) x (2^l) containing original signal as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. The wavelet coefficients are returned in array y, containing the coarsest approximation in y [0][0] followed by horizontal/vertical details in order of increasing resolution/frequency. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. It is recommended to preallocate yt to the maximum expected size. The dimensions of yt are not checked; this is the caller's responsibility. Returns 0 on success, else -1. */ { const unsigned len = 1 << l; unsigned li; WaveletMatrix ytloc = NULL; /* Skip transform if input too short or missing */ if (l < 1 || !y) return -1; if (!yt) /* No buffer supplied; allocate one on demand */ if (!(yt = ytloc = allocWaveletMatrix(l))) { fprintf(stderr, "ERROR - Failed allocating %dx%d buffer array" " in WaveletXform2()", len, len); return -1; } for (li = l; li > 1; li--) { /* Apply horizontal & vertical Daubechies D4 transform, swapping input and transposed output array */ if (d4Step(y, yt, li) || d4Step(yt, y, li)) return -1; } /* Apply horizontal & vertical Haar transform at coarsest resolution (li==1) to obtain single approximation coefficient at y [0][0]; all other coeffs are details. */ if (haarStep(y, yt, li) || haarStep(yt, y, li)) return -1; /* NOTE: All coefficients now in y */ if (ytloc) /* Free yt if allocated on demand */ freeWaveletMatrix(ytloc, l); return 0; } int waveletInvXform2 (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Perform full 2D multiresolution inverse wavelet transform on array y of size (2^l) x (2^l) containing wavelet coefficients as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. It is recommended to preallocate yt to the maximum expected size. The dimensions of yt are not checked; this is the caller's responsibility. The reconstructed signal is returned in array y. Returns 0 on success, else -1. */ { const unsigned len = 1 << l; unsigned li; WaveletMatrix ytloc = NULL; /* Skip inverse transform if input too short or missing */ if (l < 1 || !y) return -1; if (!yt) /* No buffer supplied; allocate one on demand */ if (!(yt = ytloc = allocWaveletMatrix(l))) { fprintf(stderr, "ERROR - Failed allocating %dx%d buffer array" " in WaveletInvXform2()", len, len); return -1; } /* Invert horizontal & vertical Haar transform at coarsest level (li==1), swapping input and transposed output array */ if (haarInvStep(y, yt, 1) || haarInvStep(yt, y, 1)) return -1; for (li = 2; li <= l; li++) { /* Invert horizontal & vertical Daubechies D4 transform, swapping input and transposed output arrays */ if (d4InvStep(y, yt, li) || d4InvStep(yt, y, li)) return -1; } /* NOTE: Reconstructed signal now in y */ if (ytloc) /* Free yt if allocated on demand */ freeWaveletMatrix(ytloc, l); return 0; } #ifdef WAVELET_TEST #include + int main (int argc, char *argv []) { int i, j, k, l; - unsigned len; + unsigned len, numThresh = 0; WaveletMatrix y0 = NULL, y = NULL; FILE *dataFile = NULL; - float inData; + float inData, thresh = 0; if (argc < 2) { - fprintf(stderr, "%s [dataFile]\n", argv [0]); - fputs("Missing array resolution l > 1\n", stderr); + fprintf(stderr, "%s [threshold] [dataFile]\n", argv [0]); + fputs("Missing array resolution l > 1, " + "compression threshold >= 0\n", stderr + ); return -1; } if (!(l = atoi(argv [1])) || l < 1) { fputs("Invalid array resolution l\n", stderr); return -1; } else len = 1 << l; + if (argc > 2 && (thresh = atof(argv [2])) < 0) { + fprintf(stderr, "Invalid threshold %.3f\n", thresh); + return -1; + } /* Allocate arrays for original and reconstruction */ if (!(y0 = allocWaveletMatrix(l)) || !(y = allocWaveletMatrix(l))) { fprintf(stderr, "Failed allocating %dx%d array\n", len, len); return -1; } - if (argc > 2) { + if (argc > 3) { /* Load data from file; length must not exceed allocated */ - if (!(dataFile = fopen(argv [2], "r"))) { - fprintf(stderr, "Failed opening data file %s\n", argv [2]); + if (!(dataFile = fopen(argv [3], "r"))) { + fprintf(stderr, "Failed opening data file %s\n", argv [3]); return -1; } for (i = 0; i < len; i++) { for (j = 0; j < len; j++) { if (feof(dataFile)) { fprintf(stderr, "Premature end of file reading data from %s\n", argv [2] ); fclose(dataFile); return -1; } /* Read next float, skipping any leading whitespace */ if (fscanf(dataFile, " %f", &inData)) { y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = y [i][j][0] = y [i][j][1] = y [i][j][2] = (WAVELET_COEFF)inData; } else { fprintf(stderr, "Error reading from data file %s\n", argv [2] ); fclose(dataFile); return -1; } } } fclose(dataFile); } else { /* Init with random data */ srand48(111); for (i = 0; i < len; i++) for (j = 0; j < len; j++) - #if 1 + #if 0 for (k = 0; k < 3; k++) y0 [i][j][k] = y [i][j][k] = drand48(); #else y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = y [i][j][0] = y [i][j][1] = y [i][j][2] = drand48(); #endif } /* Forward xform */ if (waveletXform2(y, NULL, l)) { fputs("Forward xform failed\n", stderr); return -1; } + /* Threshold coefficients; we use hard thresholding as it's easier + * to implement than soft thresholding, which requires sorting the + * coefficients */ + /* NOTE: y [0][0] is omitted as it's the coarsest + * smooth/avg/approx/lowpass coefficient! */ + for (i = 0; i < len; i++) + for (j = 0; j < len; j++) { + if (coeffThresh(y [i][j], thresh) && (i | j)) { + y [i][j][0] = y [i][j][1] = y [i][j][2] = 0; + numThresh++; + #if 0 + /* Replace thresholded values with random noise in range + [-threshold, threshold] */ + y [i][j][0] = y [i][j][1] = y [i][j][2] = thresh * ( + 2 * drand48() - 1 + ); + #endif + } + } + #ifdef WAVELET_DBG /* Dump coefficients */ puts("-----------------------------------------------------------\n"); - dumpCoeffs(y, NULL, l); + if (numThresh) + printf("%d/%d coefficients thresholded\n", + numThresh, len * len + ); + + dumpCoeffs(y, NULL, l, thresh); puts("\n-----------------------------------------------------------\n"); #endif /* Inverse xform */ if (waveletInvXform2(y, NULL, l)) { fputs("Inverse xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG puts("-----------------------------------------------------------\n"); puts("ORIG vs. INV XFORM"); - dumpCoeffs(y0, y, l); + dumpCoeffs(y0, y, l, 0); #endif - + + printf("\nAvg RMSE = %.2f\n", rmseCoeffs(y0, y, l)); + freeWaveletMatrix(y0, l); freeWaveletMatrix(y, l); return 0; } #endif