diff --git a/Rmakefile b/Rmakefile index 6509647..81f363d 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,547 +1,553 @@ # 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 + pmcontribsort.o wavelet2.o wavelet3.o mrgbe.o PMRCSRC = pmapcontrib.c pmcontrib2.c pmcontrib3.c pmcontrib4.c \ - pmcontribsort.c wavelet2.c mrgbe.c + pmcontribsort.c wavelet2.c wavelet3.c mrgbe.c HEADERS = ambient.h ray.h data.h otspecial.h source.h # # What this makefile produces: # PROGS = $(DESTDIR)/rtrace $(DESTDIR)/rpict $(DESTDIR)/rvu $(DESTDIR)/rcontrib \ $(DESTDIR)/lookamb $(DESTDIR)/mkpmap $(DESTDIR)/pmapdump all: $(PROGS) $(RCLIB) $(SPECIAL) install: all rayinit.cal $(INSTALL) $(PROGS) $(INSTDIR) rm -f $(LIBDIR)/rayinit.cal cp rayinit.cal $(LIBDIR) ogl: clean: set nonomatch; rm -f $(PROGS) *.o lint: $(RVSRC) $(LINT) $(LINTFLAGS) -DRVIEW $(RVSRC) $(LIBS) # Preocmputed contrib pmap unit tests -wavelet-test: wavelet2.c wavelet2.h mrgbe.c mrgbe.h +wavelet2-test: wavelet2.c wavelet2.h mrgbe.c mrgbe.h $(CC) $(CFLAGS) -DWAVELET_TEST -DWAVELET_TEST_mRGBE -DWAVELET_DBG \ - -o wavelet-test -g wavelet2.c mrgbe.c $(MLIB) + -o wavelet2-test -g wavelet2.c mrgbe.c $(MLIB) + +wavelet3-test: wavelet3.c wavelet2.h mrgbe.c mrgbe.h + $(CC) $(CFLAGS) -DWAVELET_TEST -DWAVELET_TEST_mRGBE -DWAVELET_DBG \ + -o wavelet3-test -g wavelet3.c mrgbe.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 +wavelet3.o: wavelet3.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/wavelet.c b/wavelet.c index f0cf392..690c058 100644 --- a/wavelet.c +++ b/wavelet.c @@ -1,595 +1,595 @@ /* ========================================================================= - Daubechies D4 and Haar wavelet transforms on arbitrary sized arrays (i.e. - not restricted to powers of 2) of 3-tuples. Array boundaries are - extended cyclically. + Daubechies D4 and Haar wavelet transforms on arbitrary sized arrays + of 3-tuples. Array boundaries are extended as constants, which introduces + padding coefficients. The implementation here is mainly distilled from the following sources: -- Filip Wasilewski's pywavelets: https://github.com/PyWavelets/pywt/tree/master/pywt/_extensions/c [D4 with boundary condition handling via padding, arbitrary sizes] -- Ingrid Daubechies and Wim Sweldens, "Factoring Wavelet Tranforms into Lifting Steps", section 7.5: https://services.math.duke.edu/~ingrid/publications/J_Four_Anal_Appl_4_p247.pdf [Lifting scheme on D4] -- Ian Kaplan's webpages on wavelets and signal processing:: http://bearcave.com/misl/misl_tech/wavelets/daubechies/index.html http://bearcave.com/misl/misl_tech/wavelets/lifting/basiclift.html [Lifting scheme on D4] -- Stefan Moebius' TurboWavelets.Net: https://github.com/codeprof/TurboWavelets.Net/tree/master/TurboWavelets [Boundary condition handling with arbitrary sizes, backreferencing] -- "Numerical Recipes in C", chapter 13.10: http://cacs.usc.edu/education/cs653/NRC-c13.10-Wavelets.pdf [Linear algebra (matrix mult.) implementation of D4] 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 "wavelet.h" #include #include #include #include #define COEFSIZE (3 * sizeof(Coeff)) #define min(a, b) ((a) < (b) ? (a) : (b)) /* Daubechies D4 wavelet coeffs */ static const Coeff norm = 0.25 / sqrt(2), sqrt3 = sqrt(3), h [4] = { norm * (1 + sqrt3), norm * (3 + sqrt3), norm * (3 - sqrt3), norm * (1 - sqrt3) }, g [4] = { h [3], -h [2], h [1], -h [0] }, /* Summed coeffs for padding at left boundary */ hBound = h [0] + h [1] + h [2], gBound = g [0] + g [1] + g [2]; /* Haar wavelet coeffs */ static const Coeff c0 = 1 / sqrt(2); static int tupleCmp (const void *t0, const void *t1) /* Comparison callback for sorting 3-tuples; reduce by summation, assuming unsigned. */ { int i; Coeff sum0 = 0, sum1 = 0, diff; for (i = 0; i < 3; i++) { sum0 += (*(Coeff **)t0) [i]; sum1 += (*(Coeff **)t1) [i]; } diff = sum0 - sum1; return diff > 0 ? 1 : diff < 0 ? -1 : 0; } void waveletPresort (Coeff (*x) [3], int xIdx [], int xLen) /* Presort x to obtain a symmetric ("gaussian") distribution. This minimises discontinuities in the signal, thus minimising the resulting detail coefficients over the range of the signal. This essentially optimises the compression. Returns resorted signal in x and the permuted indices in xIdx to subsequently reconstruct the original order via waveletPostsort().*/ { int i, j; Coeff *xp [xLen], *yp [xLen], y [xLen][3]; /* Initialise array of pointers to sort x via indirection, thus permuting the pointers rather than the signal itself. */ for (i = 0; i < xLen; i++) xp [i] = (Coeff *)&x [i]; /* Sort pointers in xp by the values they refer to (see deref in tupleCmp() above) */ qsort(xp, xLen, sizeof(Coeff *), tupleCmp); /* Redistribute xp symmetrically from both boundaries towards centre. Chuck into yp.*/ for (i = 0; i < xLen - 1; i += 2) { j = i >> 1; yp [j] = xp [i]; yp [xLen - j - 1] = xp [i + 1]; } if (xLen & 1) /* Last partial iteration if odd size */ yp [(i >> 1) + 1] = xp [i]; /* Transfer values referenced by yp into a temp array (to prevent clobbering x and SEGFAULTing), and transfer their corresponding indices into xIdx. */ for (i = 0; i < xLen; i++) { memcpy(y [i], yp [i], COEFSIZE); /* XXX: Assume tuples in x are allocated consecutively! */ xIdx [i] = (yp [i] - x [0]) / 3; } /* Transfer temp to x */ memcpy(x, y, xLen * COEFSIZE); } static void haarStep (Coeff (*x) [3], int len, Coeff (*y) [3]) /* Single step (at fixed resolution) of a forward Haar wavelet transform on array x of arbitrary size len containing 3-tuples. Note this is a 1D transform, and the triplets per array element are _not_ decorrelated, but transformed independently. The wavelet coefficients are returned in the preallocated array y. */ { const unsigned hl = -(-len >> 1); int i, hi, j; if (len >= 2) { /* First and intermediate iterations */ for (i = 0; i <= len - 2; i += 2) { hi = i >> 1; for (j = 0; j < 3; j++) { /* Smooth */ y [hi][j] = c0 * (x [i][j] + x [i + 1][j]); /* Detail */ y [hl + hi][j] = c0 * (x [i][j] - x [i + 1][j]); } } if (len & 1) { /* Odd length; pad with last value as additional approx term */ for (j = 0; j < 3; j++) y [hl - 1][j] = x [len - 1][j]; } #ifdef _WAVELET_DBG printf("FWD HAAR XFORM (len = %d)\n", len); for (i = 0; i < len; i++) printf("%.3f\t-->\t%.3f\n", x [i][0], y [i][0]); putchar('\n'); #endif } } static void haarInvStep (Coeff (*y) [3], int len, Coeff (*x) [3]) /* Single step (at fixed resolution) of an inverse Haar wavelet transform on coefficient array y of size len containing 3-tuples. This is essentially reverses the forward transform above; see additional comments there for implementation details. The reconstructed signal is return in the preallocated array x. */ { const unsigned hl = -(-len >> 1); int i, hi, j; if (len >= 2) { for (i = 0; i <= len - 2; i += 2) { hi = i >> 1; for (j = 0; j < 3; j++) { x [i][j] = c0 * (y [hi][j] + y [hl + hi][j]); x [i + 1][j] = c0 * (y [hi][j] - y [hl + hi][j]); } } if (len & 1) { /* Odd length; pad additional value from last approx term */ for (j = 0; j < 3; j++) x [len - 1][j] = y [hl - 1][j]; } #ifdef _WAVELET_DBG printf("INV HAAR XFORM (len = %d)\n", len); for (i = 0; i < len; i++) printf("%.3f\t-->\t%.3f\n", y [i][0], x [i][0]); putchar('\n'); #endif } } static unsigned d4Step (Coeff (*x) [3], int xLen, Coeff (*y) [3]) /* Single step (at fixed resolution) of forward Daubechies D4 wavelet transform on array x of arbitrary size xLen containing 3-tuples. Note this is a 1D transform, and the triplets per array element are _not_ decorrelated, but transformed independently. The approximation and detail coefficients are returned in the preallocated output array y. This array is larger than xLen as the transform introduces padding coefficients to handle boundary effects. This size is obtained as return value, and can be interrogated beforehand without performing a transform by passing either x or y as NULL pointer. */ { unsigned i, j, k, hi; const unsigned yLen = xLen + 3 >> 1, /* Output length + boundary pad */ hy = yLen; /* Offset for detail coeffs */ if (x && y) { if (xLen >= 2) { #ifdef _WAVELET_DBG /* Zero output array to facilitate debugging */ memset(y, 0, 2 * yLen * COEFSIZE); #endif /* First iteration (i = -2); pad coeffs at left edge with constant * boundary extension (note repeated x [0]) */ for (j = 0; j < 3; j++) { /* Approximation */ y [0][j] = hBound * x [0][j] + h [3] * x [1][j]; /* Detail */ y [hy][j] = gBound * x [0][j] + g [3] * x [1][j]; } /* Intermediate and last iterations (i >= 0) */ for (i = 0; i < xLen; i += 2) { hi = (i >> 1) + 1; for (j = 0; j < 3; j++) { y [hi][j] = y [hy + hi][j] = 0; /* Convolve up to wavelet or input boundary, whichever comes * first */ for (k = i; k < min(i + 4, xLen); k++) { /* Approximation */ y [hi][j] += h [k - i] * x [k][j]; /* Detail */ y [hy + hi][j] += g [k - i] * x [k][j]; } /* If necessary, convolve beyond input boundary with constant * extension (note repeated x [xLen - 1]). This is skipped if * the wavelet boundary was already reached */ for (; k < i + 4; k++) { /* Approximation */ y [hi][j] += h [k - i] * x [xLen - 1][j]; /* Detail */ y [hy + hi][j] += g [k - i] * x [xLen - 1][j]; } } } #ifdef _WAVELET_DBG printf( "D4 XFORM + boundary padding (len = %d -> %d)\n", xLen, yLen << 1 ); for (i = 0; i < yLen << 1; i++) { if (i < xLen) printf("%.3f\t-->\t", x [i][0]); else printf("\t\t"); printf("%.3f\n", y [i][0]); } putchar('\n'); #endif } } return yLen; } static unsigned d4InvStep (Coeff (*y) [3], int yLen, Coeff (*x) [3]) /* Single step (at fixed resolution) of inverse Daubechies D4 wavelet transform on coefficient array y of size yLen (i.e. yLen/2 approx and detail coeffs each) containing 3-tuples. This essentially reverses the forward transform above; see additional comments there for implementation details. The reconstructed signal is returned in the preallocated array x. This array is smaller than yLen as boundary padding coefficients are collapsed in the inverse transform. This size is obtained as returned value, and can be interrogated beforehand without performing a transform by passing either x or y as NULL pointer. */ { unsigned i, j, k, hi, xLen = yLen - 2; /* Max size of output array x */ const unsigned hy = yLen >> 1; /* Offset for detail coeffs */ if (x && y) { if (yLen >= 2) { /* Reconstruct from approximation and detail coefficients */ for (i = 0; i < xLen - 1; i += 2) { hi = i >> 1; for (j = 0; j < 3; j++) { x [i][j] = h [2] * y [hi][j] + g [2] * y [hy + hi][j] + h [0] * y [hi + 1][j] + g [0] * y [hy + hi + 1][j]; x [i + 1][j] = h [3] * y [hi][j] + g [3] * y [hy + hi][j] + h [1] * y [hi + 1][j] + g [1] * y [hy + hi + 1][j]; } } xLen = i; #ifdef _WAVELET_DBG printf("D4 INV XFORM (len = %d --> %d)\n", yLen, xLen); for (i = 0; i < yLen; i++) { printf("%.3f\t", y [i][0]); if (i < xLen) printf("-->\t%.3f\n", x [i][0]); else putchar('\n'); } putchar('\n'); #endif } } return xLen; } unsigned waveletXform (Coeff (*x) [3], unsigned xLen, Coeff (*y) [3]) /* Perform full multiresolution forward wavelet transform on array x of length xLen containing 3-tuples. Note no intra-tuple transform occurs (1-D only). The wavelet coefficients are returned in preallocated array y, containing the coarsest approximation in y [0] followed by details in order of increasing resolution/frequency. Note y must accommodate more than xLen coefficients to accommodate padding for additional boundary coefficients. The number of coefficients can be interrogated beforehand by calling the function with either x or y set to NULL. Note also that x is used as a buffer for intermediate results, and is therefore modified. The returned value is the number of output coefficients in y, or 0 if the transform failed, e.g. if the input is too short. */ { /* Output size (number of coeffs, including boundary padding), output * array offset and pointers */ unsigned numCoeffs = 0, yOff, yLen, l; Coeff (*yOffPtr) [3] = NULL; /* Skip transform if input too short */ if (xLen > 3) { /* Figure out total number of coefficients (= minimum size of output * array including boundary padding). */ for ( numCoeffs = 0, l = xLen; l > 3; numCoeffs += l = d4Step(NULL, l, NULL) ); numCoeffs += l; if (x && y) { /* Apply Daubechies D4 transform, accounting for boundary padding */ for (l = xLen, yOff = numCoeffs; l > 3; l = yLen) { /* Get num coeffs and calculate offset into output array * (subtracting num coeffs once more to account for approx * coeffs, which will be replaced in next iteration) */ yOff -= yLen = d4Step(NULL, l, NULL); yOffPtr = y + yOff - yLen; /* Resulting wavelet coefficients in y are prepended to those * from previous iteration; consequently, returned num coeffs is * input size for next iteration */ d4Step(x, l, yOffPtr); /* Copy approx coeffs from y into x as input for next iteration */ memcpy(x, yOffPtr, yLen * COEFSIZE); } /* Apply Haar transform in penultimate and final (since len is odd) * iterations to decompose to a single approximation and detail * coefficient at coarsest resolution at y [0] and y [1], * respectively. Obviously no boundary padding required here. */ for (; l >= 2; l--) { haarStep(x, l, y); /* Copy coeffs from y into x as input for next iteration */ if (l > 2) memcpy(x, y, l * COEFSIZE); } } } /* All coeffs now in y; return their number */ return numCoeffs; } unsigned waveletInvXform ( Coeff (*y) [3], unsigned yLen, Coeff (*x) [3], unsigned xLen0 ) /* Perform full multiresolution inverse wavelet transform on array y of * length yLen containing wavelet coefficients as 3-tuples. Note no * intra-tuple transform occurs (1-D only). The original signal length xLen0 * is required for the reconstruction. * * The reconstructed signal is returned in the preallocated array x. This * is shorter than yLen since the additional boundary coefficients in y are * collapsed. Note y is used as buffer for intermediate results, and is * therefore modified. * * The returned value is the reconstructed signal length, which should match * xLen, or 0 if the transform failed. */ { unsigned xLen = 0, yOff, i, j, l, nextl, nIter; /* The temporary output buffer xTmp is passed to d4InvStep() for * intermediate xform steps. The additional space is required if xLen0 * is odd, as an additional coefficient is reconstructed (which is * subsequently dropped). If x were to be passed directly to * d4InvStep() under these circumstances, it would overflow. */ Coeff (*yOffPtr) [3] = NULL, xTmp [xLen0 + 1][3]; /* Figure out number of iterations based on original signal length */ for (nIter = 0, l = xLen0; l > 3; nIter++, l = d4Step(NULL, l, NULL)); /* Invert Haar transform at two coarsest levels */ for (l = 2; l <= 3; l++) { haarInvStep(y, l, x); /* Copy reconstructed signal into y for next iteration */ memcpy(y, x, l * COEFSIZE); } /* Invert Daubechies D4 transform */ for (i = nIter; i; i--) { /* Figure out length of current detail coeffs based on number of * iterations */ for ( j = 0, yOff = yLen, l = nextl = xLen0; j < i; j++, nextl = l, yOff -= l = d4Step(NULL, l, NULL) ); #ifdef _WAVELET_DBG /* Zero output buffer to facilitate debugging */ memset(xTmp, 0, sizeof(xTmp)); #endif yOffPtr = y + yOff; /* Drop excess reconstructed approximation coeffs to match number of * detail coeffs in next iteration */ xLen = min(d4InvStep(yOffPtr - l, l << 1, xTmp), nextl); /* Copy reconstructed coeffs/signal from output buffer xTmp into y for * next iteration, or into x after last iteration */ memcpy(i > 1 ? yOffPtr + l - xLen : x, xTmp, xLen * COEFSIZE); } return xLen; } #ifdef _WAVELET_TEST #include #define ERRSTRLEN 1024 int main (int argc, char *argv []) { int i, j, xLen0, yLen, xLen, *xIdx; Coeff (*xOrg) [3], (*x) [3], (*y) [3] = NULL; char errstr [ERRSTRLEN]; if (argc < 2) { fputs("Missing array length\n", stderr); return -1; } if (!(xLen0 = atoi(argv [1]))) { fputs("Invalid array length\n", stderr); return -1; } /* Allocate arrays for original and reconstruction (num coeffs is * obtained from waveletXform() before allocating output array) */ if ( !(xOrg = calloc(xLen0, COEFSIZE)) || !(xIdx = calloc(xLen0, sizeof(int))) || !(x = calloc(xLen0, COEFSIZE)) || !(y = calloc(waveletXform(NULL, xLen0, NULL), COEFSIZE)) ) { fputs("Failed array allocation\n", stderr); return -1; } /* Init with random data */ srand48(111); for (i = 0; i < xLen0; i++) { #if 1 for (j = 0; j < 3; j++) xOrg [i][j] = drand48(); #else xOrg [i][0] = xOrg [i][1] = xOrg [i][2] = drand48(); #endif } /* Optimise signal distribution */ waveletPresort(xOrg, xIdx, xLen0); #ifdef _WAVELET_DBG /* Dump presorted signal */ puts("-----------------------------------------------------------"); puts("Presorted Original & Indices"); for (i = 0; i < xLen0; i++) printf( "%.3f\t%d\n", xOrg [i][0] + xOrg [i][1] + xOrg [i][2], xIdx [i] ); puts("-----------------------------------------------------------\n"); #endif /* Transform and reconstruct */ memcpy(x, xOrg, xLen0 * COEFSIZE); if (!(yLen = waveletXform(x, xLen0, y))) { puts("Forward xform failed\n"); return -1; } #ifdef _WAVELET_DBG /* Dump coefficients */ puts("-----------------------------------------------------------\n"); printf("%d coefficients:\n\n", yLen); for (i = 0; i < yLen; i++) printf("%.3f\t%.3f\t%.3f\n", y [i][0], y [i][1], y [i][2]); puts("\n-----------------------------------------------------------\n"); #endif if (!(xLen = waveletInvXform(y, yLen, x, xLen0))) { puts("Inverse xform failed\n"); return -1; } if (xLen != xLen0) { printf( "Original length (%d) != reconstructed length (%d)\n", xLen0, xLen ); return -1; } #ifdef _WAVELET_DBG puts("-----------------------------------------------------------\n"); #endif puts("ORIG\t vs.\tINV XFORM"); for (i = 0; i < xLen; i++) #if 1 printf( "%.3f\t%.3f\t%.3f\t\t%.3f\t%.3f\t%.3f\n", xOrg [i][0], xOrg [i][1], xOrg [i][2], x [i][0], x [i][1], x [i][2] ); #else printf("%.3f\t\t%.3f\n", xOrg [i][0], x [i][0]); #endif free(xOrg); free(xIdx); free(x); free(y); return 0; } #endif diff --git a/wavelet2.c b/wavelet2.c index ce56a84..d7d99b9 100644 --- a/wavelet2.c +++ b/wavelet2.c @@ -1,832 +1,846 @@ /* ========================================================================= - 2-dimensional wavelet transform on (2^l) x (2^l) sized arrays of 3-tuples, - where l > 1. + 2-dimensional Daubechies D4 and Haar wavelet transforms on (2^l) x (2^l) + sized arrays of 3-tuples, where l > 1. Array boundaries are extended as + periodic, which does not require padding coefficients. 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$ */ #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 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] }; */ +/* Summed coeffs for padding at left boundary */ +static const WAVELET_COEFF h4Bound = h4 [0] + h4 [1] + h4 [2], + g4Bound = g4 [0] + g4 [1] + g4 [2]; -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. -*/ + +WaveletMatrix allocWaveletMatrix (unsigned len) +/* Allocate and return a 2D coefficient array of size (len x len). + Returns NULL if allocation failed. */ { - const unsigned len = 1 << l; unsigned i; WaveletMatrix y = NULL; - if (l >= 1) { + if (len >= 2) { 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) +void freeWaveletMatrix (WaveletMatrix y, unsigned len) /* - Free previously allocated 2D coefficient array y of size (2^l) x (2^l) + Free previously allocated 2D coefficient array y of size (len x len) */ { 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) + static void zeroCoeffs (WaveletMatrix y, unsigned len) /* 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, float thresh + unsigned yLen, unsigned ytLen, float thresh ) - /* Dump arrays y and yt side-by-side to stdout (skip yt if NULL) */ + /* Dump arrays y and yt of dims yLen resp. ytLen 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( + for (i = 0; i < max(yLen, ytLen); i++) { + if (y && i < yLen) { + for (j = 0; j < yLen; j++) + printf( #ifdef WAVELET_FINALHAAR - (i | j) && coeffThresh(y [i][j], thresh) + (i | j) && coeffThresh(y [i][j], thresh) #else - (i > 1 || j > 1) && coeffThresh(y [i][j], thresh) + (i > 1 || j > 1) && coeffThresh(y [i][j], thresh) #endif - ? "[% 5.2f]\t" : " % 5.2f\t", - coeffAvg(y [i][j]) - ); + ? "[% 5.2f]\t" : " % 5.2f\t", + coeffAvg(y [i][j]) + ); + } - if (yt) { + if (yt && i < ytLen) { printf(" -->\t"); - for (j = 0; j < len; j++) + for (j = 0; j < ytLen; j++) printf("% 5.2f\t", coeffAvg(yt [i][j])); } putchar('\n'); } } static float rmseCoeffs (const WaveletMatrix y1, const WaveletMatrix y2, - unsigned l + unsigned len ) /* 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++) { #if 0 d = (coeffAvg(y1 [i][j]) - coeffAvg(y2 [i][j])) / coeffAvg(y1 [i][j]); #else d = coeffAvg(y1 [i][j]) - coeffAvg(y2 [i][j]); #endif 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). + 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; +#ifdef WAVELET_DBG + static unsigned axis = 0; +#endif if (len < 2 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs(yt, l); + zeroCoeffs(yt, len); #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, 0); + dumpCoeffs(y, yt, len, 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; +#ifdef WAVELET_DBG + static unsigned axis = 1; +#endif if (len < 2 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs(yt, l); + zeroCoeffs(yt, len); #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, 0); + dumpCoeffs(y, yt, len, 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; +#ifdef WAVELET_DBG + static unsigned axis = 0; +#endif if (len < 4 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs(yt, l); + zeroCoeffs(yt, len); #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, 0); + dumpCoeffs(y, yt, len, 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; +#ifdef WAVELET_DBG + static unsigned axis = 1; +#endif if (len < 4 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs(yt, l); + zeroCoeffs(yt, len); #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, 0); + dumpCoeffs(y, yt, len, 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..1][0..1] (or just y [0][0] with WAVELET_FINALHAAR defined) 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; } #ifdef WAVELET_FINALHAAR /* 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; #endif /* 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; } #ifdef WAVELET_FINALHAAR /* 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; #endif 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 #ifdef WAVELET_TEST_mRGBE #include "mrgbe.h" #endif + #define WAVELET_TEST_INIT 0 + int main (int argc, char *argv []) { int i, j, k, l; unsigned len, numThresh = 0; WaveletMatrix y0 = NULL, y = NULL; FILE *dataFile = NULL; WAVELET_COEFF inData, thresh = 0; #ifdef WAVELET_TEST_mRGBE #define HUGE 1e10 mRGBE mrgbeCoeff; mRGBERange mrgbeRange; WaveletMatrix ymrgbe = NULL; #endif if (argc < 2) { 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)) #ifdef WAVELET_TEST_mRGBE || !(ymrgbe = allocWaveletMatrix(l)) #endif ) { fprintf(stderr, "Failed allocating %dx%d array\n", len, len); return -1; } if (argc > 3) { /* Load data from file; length must not exceed allocated */ 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, " %lf", &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] = inData; } else { fprintf(stderr, "Error reading from data file %s\n", argv [2] ); fclose(dataFile); return -1; } } } fclose(dataFile); } else { /* Init input */ srand48(111); - for (i = 0; i < len; i++) - for (j = 0; j < len; j++) - #if 0 - /* Random data, channel-independent */ - for (k = 0; k < 3; k++) - y0 [i][j][k] = y [i][j][k] = drand48(); - #endif - #if 0 - /* Random data, indentical for all channels */ - 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 - #if 0 /* Monotonically increasing along axis 0 */ - y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = - y [i][j][0] = y [i][j][1] = y [i][j][2] = i; - #endif - #if 0 /* Monotonically increasing along axis 1 */ - y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = - y [i][j][0] = y [i][j][1] = y [i][j][2] = j; - #endif - #if 0 /* Monotonically increasing along both axes */ - y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = - y [i][j][0] = y [i][j][1] = y [i][j][2] = i * j; - #endif - #if 1 /* Monotonically increasing by linear index */ - y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = - y [i][j][0] = y [i][j][1] = y [i][j][2] = i * len + j; + for (i = 0; i < len; i++) { + for (j = 0; j < len; j++) { + for (k = 0; k < 3; k++) { + y0 [i][j][k] = y [i][j][k] = + #if WAVELET_TEST_INIT == 0 + /* Random data, channel-independent */ + drand48(); + #elif WAVELET_TEST_INIT == 1 + /* Random data, indentical for all channels */ + k ? y [i][j][k - 1] : drand48(); + #elif WAVELET_TEST_INIT == 2 + /* Monotonically increasing along axis 0 */ + i; + #elif WAVELET_TEST_INIT == 3 + /* Monotonically increasing along axis 1 */ + j; + #elif WAVELET_TEST_INIT == 4 + /* Monotonically increasing along both axes */ + i * j; + #elif WAVELET_TEST_INIT == 5 + /* Monotonically increasing by linear index */ + i * len + j; #endif + } + } + } } #ifdef WAVELET_DBG puts("----------------------- FWD XFORM -------------------------\n"); #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..1][0..1] (or just y [0][0] if WAVELET_FINALHAAR is * defned) are omitted as they are the coarsest approximation * coefficients. */ for (i = 0; i < len; i++) for (j = 0; j < len; j++) { #ifdef WAVELET_FINALHAAR if ((i | j) && coeffThresh(y [i][j], thresh)) { #else if ((i > 1 || j > 1) && coeffThresh(y [i][j], thresh)) { #endif 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("----------------------- COEFFICIENTS -------------------------\n"); - dumpCoeffs(y, NULL, l, thresh); + dumpCoeffs(y, NULL, len, thresh); if (numThresh) printf("\n%d/%d coefficients thresholded = %.1f%% compression", numThresh, len * len, 100. * numThresh / (len * len) ); #endif #ifdef WAVELET_TEST_mRGBE /* Test mRGBE coefficient encoding */ mrgbeRange.min [0] = mrgbeRange.min [1] = mrgbeRange.min [2] = HUGE; mrgbeRange.max [0] = mrgbeRange.max [1] = mrgbeRange.max [2] = 0; /* Find min/max coeff and init mRGBE range */ for (i = 0; i < len; i++) for (j = 0; j < len; j++) for (k = 0; k < 3; k++) { inData = fabs(y [i][j][k]); if (inData < mrgbeRange.min [k]) mrgbeRange.min [k] = inData; if (inData > mrgbeRange.max [k]) mrgbeRange.max [k] = inData; }; mRGBEinit(&mrgbeRange, mrgbeRange.min, mrgbeRange.max); /* Encode coeffs to mRGBE and back, but preserve approximation * coefficient at y [0][0] */ for (i = 0; i < len; i++) for (j = 0; j < len; j++) if (i | j) { mrgbeCoeff = mRGBEencode(y [i][j], &mrgbeRange, 0); mRGBEdecode(mrgbeCoeff, &mrgbeRange, ymrgbe [i][j]); } else for (k = 0; k < 3; k++) ymrgbe [i][j][k] = y [i][j][k]; #ifdef WAVELET_DBG /* Dump mRGBE-decoded coefficients */ puts("\n\n-------------------- mRGBE COEFFICIENTS ----------------------\n"); - dumpCoeffs(y, ymrgbe, l, 0); + dumpCoeffs(y, ymrgbe, len, 0); #endif #endif #ifdef WAVELET_DBG puts("\n----------------------- INV XFORM -------------------------\n"); #endif /* Inverse xform (also using mRGBE coeffs if enabled) */ if (waveletInvXform2(y, NULL, l) #ifdef WAVELET_TEST_mRGBE || waveletInvXform2(ymrgbe, NULL, l) #endif ) { fputs("\nInverse xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG puts("\n--------------------- ORIG vs. INV XFORM ------------------------\n"); - dumpCoeffs(y0, y, l, 0); + dumpCoeffs(y0, y, len, 0); #endif - printf("\nAvg RMSE = %.2f\n", rmseCoeffs(y0, y, l)); + printf("\nAvg RMSE = %.2f\n", rmseCoeffs(y0, y, len)); #ifdef WAVELET_TEST_mRGBE #ifdef WAVELET_DBG puts("\n------------------ ORIG vs. INV XFORM + mRGBE -------------------\n"); - dumpCoeffs(y0, ymrgbe, l, 0); + dumpCoeffs(y0, ymrgbe, len, 0); #endif - printf("\nAvg RMSE with mRGBE enc = %.2f\n", rmseCoeffs(y0, ymrgbe, l)); + printf("\nAvg RMSE with mRGBE enc = %.2f\n", + rmseCoeffs(y0, ymrgbe, len) + ); #endif freeWaveletMatrix(y0, l); freeWaveletMatrix(y, l); #ifdef WAVELET_TEST_mRGBE freeWaveletMatrix(ymrgbe, l); #endif return 0; } #endif diff --git a/wavelet2.h b/wavelet2.h index aa2b274..ef4d974 100644 --- a/wavelet2.h +++ b/wavelet2.h @@ -1,97 +1,96 @@ /* ========================================================================= - 2-dimensional wavelet transform on (2^l) x (2^l) sized arrays of 3-tuples, - where l >= 1. + Definitions for wavelet transforms on arrays of 3-tuples. 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 _WAVELET2_H #define _WAVELET2_H /* Run a final Haar wavelet transform on the DB4 coefficients at the * coarsest resolution to obtain a single approximation coefficient (as * opposed to 2 x 2 with just DB4. */ #define WAVELET_FINALHAAR /* Wavelet coefficient type; defaults to double if not already defined. (may be overridden in compiler command line, for example). NOTE: single precision float may improve cache performance during 2D transform with large arrays. */ #ifndef WAVELET_COEFF #define WAVELET_COEFF double #endif /* Wavelet matrix defs; these are stored as arrays of pointers (a.k.a. * Iliffe vectors), as this proved to be more performant than a flattened * representation. Encapsulating the coefficient's triplets in a struct * prevents the compiler from introducing another layer of indirection. * */ typedef WAVELET_COEFF WaveletCoeff3 [3]; #define WAVELET_COEFFSIZE (sizeof(WaveletCoeff3)) typedef WaveletCoeff3 **WaveletMatrix; 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. */ void freeWaveletMatrix (WaveletMatrix y, unsigned l); /* Free previously allocated 2D coefficient array y of size (2^l) x (2^l) */ 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. Ideally, yt should be preallocated to the maximum expected size. The dimensions of a preallocated yt are not checked; this is the caller's responsibility. Returns 0 on success, else -1. */ 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. Ideally, yt should be preallocated to the maximum expected size. The dimensions of a preallocated yt are not checked; this is the caller's responsibility. The reconstructed signal is returned in array y. Returns 0 on success, else -1. */ #endif diff --git a/wavelet3.c b/wavelet3.c new file mode 100644 index 0000000..11937b8 --- /dev/null +++ b/wavelet3.c @@ -0,0 +1,663 @@ +/* + ========================================================================= + Daubechies D4 and Haar wavelet transforms on arbitrary sized arrays + of 3-tuples. Array boundaries are extended as constants, which introduces + padding coefficients. + + The implementation here is mainly distilled from the following sources: + -- Filip Wasilewski's pywavelets: + https://github.com/PyWavelets/pywt/tree/master/pywt/_extensions/c + [D4 with boundary condition handling via padding, arbitrary sizes] + -- Ingrid Daubechies and Wim Sweldens, "Factoring Wavelet Tranforms into + Lifting Steps", section 7.5: + https://services.math.duke.edu/~ingrid/publications/J_Four_Anal_Appl_4_p247.pdf + [Lifting scheme on D4] + -- Ian Kaplan's webpages on wavelets and signal processing:: + http://bearcave.com/misl/misl_tech/wavelets/daubechies/index.html + http://bearcave.com/misl/misl_tech/wavelets/lifting/basiclift.html + [Lifting scheme on D4] + -- Stefan Moebius' TurboWavelets.Net: + https://github.com/codeprof/TurboWavelets.Net/tree/master/TurboWavelets + [Boundary condition handling with arbitrary sizes, backreferencing] + -- "Numerical Recipes in C", chapter 13.10: + http://cacs.usc.edu/education/cs653/NRC-c13.10-Wavelets.pdf + [Linear algebra (matrix mult.) implementation of D4] + + 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 "wavelet2.h" +#include +#include +#include +#include + + +#define min(a, b) ((a) < (b) ? (a) : (b)) + + +/* 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] }; */ + +/* Summed coeffs for padding at left boundary */ +static const WAVELET_COEFF h4Bound = h4 [0] + h4 [1] + h4 [2], + g4Bound = g4 [0] + g4 [1] + g4 [2]; + + + +static unsigned padHaarStep2 (WaveletMatrix y, WaveletMatrix yt, + unsigned len +) +/* + Single step of forward 2D Haar wavelet transform with padding on array y + of arbitrary size (len x len) containing 3-tuples. This variant uses + constant boundary extension. 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 array dimension len on success, else 0. +*/ +{ + const unsigned hLen = -(-len >> 1); + int h, i, j, k; +#ifdef WAVELET_DBG + static unsigned axis = 0; +#endif + + if (len < 2 || !y || !yt) + /* Input shorter than wavelet support, or no input/output */ + return 0; + +#ifdef WAVELET_DBG + zeroCoeffs(yt, len); +#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 - 2; 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] + ); + } + } + + if (len & 1) { + /* Odd length; pad with last value as additional approx term */ + for (k = 0; k < 3; k++) + /* Smooth/approx/avg/lowpass */ + yt [hLen - 1] [i] [k] = y [i] [len - 1] [k]; + } + } + +#ifdef WAVELET_DBG + printf("%s FWD HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); + dumpCoeffs(y, yt, len, len, 0); + putchar('\n'); + axis ^= 1; +#endif + + return len; +} + + + +static unsigned padHaarInvStep2 (WaveletMatrix y, WaveletMatrix yt, + unsigned len +) +/* + Single step of inverse 2D Haar wavelet transform with padding on array y + of arbitrary size (len x len) containing 3-tuples. 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 array dimension len on success, else 0. +*/ +{ + const unsigned hLen = -(-len >> 1); + int h, i, j, k; +#ifdef WAVELET_DBG + static unsigned axis = 1; +#endif + + if (len < 2 || !y || !yt) + /* Too few coeffs for reconstruction, or no input/output */ + return 0; + +#ifdef WAVELET_DBG + zeroCoeffs(yt, len); +#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 - 2; 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] /* Detail */ + ); + + yt [i] [j + 1] [k] = h2 * ( + y [h ] [i] [k] - /* Avg */ + y [hLen + h] [i] [k] /* Detail */ + ); + } + } + + if (len & 1) { + /* Odd length; pad additional value from last avg term */ + for (k = 0; k < 3; k++) + yt [i] [len - 1] [k] = y [hLen - 1] [i] [k]; + } + } + +#ifdef WAVELET_DBG + printf("%s INV HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); + dumpCoeffs(y, yt, len, len, 0); + putchar('\n'); + axis ^= 1; +#endif + + return len; +} + + + +static unsigned padD4Step2 (WaveletMatrix y, WaveletMatrix yt, + unsigned yLen +) +/* + Single step of forward 2D Daubechies D4 wavelet transform with boundary + padding on array y containing (yLen x yLen) samples of 3-tuples. This + variant uses constant boundary extension, which results in padding + coefficients. 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. 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. + + Both the input array y and output array yt must be large enough to also + accommodate the additional padding coefficients introduced at the + boundaries. The number of output coefficient pairs ytLen returned in yt + is obtained as return value, and can be interrogated beforehand without + performing a transform by passing either y or yt as NULL. +*/ +{ + unsigned h, i, j, k, l; + /* Number of output coefficient pairs ytLen incl. padding at boundaries, + offset hyt for detail coeffs */ + const unsigned ytLen = (yLen + 3) >> 1, hyt = ytLen; +#ifdef WAVELET_DBG + static unsigned axis = 0; +#endif + + if (y && yt && yLen >= 2) { +#ifdef WAVELET_DBG + /* Zero output array to facilitate debugging */ + zeroCoeffs(yt, ytLen << 1); +#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 < yLen; i++) { + /* First iteration (j = -2); pad coeffs at left edge with constant + * boundary extension (note multiple factoring of y [i] [0] via + * aggregated coefficients h4Bound and g4Bound). */ + for (k = 0; k < 3; k++) { + /* Smooth/approx/avg/lowpass */ + yt [0 ] [i] [k] = + h4Bound * y [i] [0] [k] + + h4 [3] * y [i] [1] [k]; + + /* Detail/diff/highpass */ + yt [hyt] [i] [k] = + g4Bound * y [i] [0] [k] + + g4 [3] * y [i] [1] [k]; + } + + /* Intermediate and last iterations (j >= 0) */ + for (j = 0; j < yLen; j += 2) { + h = (j >> 1) + 1; + + for (k = 0; k < 3; k++) { + yt [h] [i] [k] = yt [hyt + h] [i] [k] = 0; + + /* Convolve up to boundary of wavelet or input, whichever + * comes first */ + for (l = j; l < min(j + 4, yLen); l++) { + /* Smooth/approx/avg/lowpass */ + y [h ] [i] [k] += h4 [l - j] * y [i] [l] [k]; + /* Detail/diff/highpass */ + y [hyt + h] [i] [k] += g4 [l - j] * y [i] [l] [k]; + } + + /* If necessary, convolve beyond input boundary with constant + * extension (note multiple factoring of y [i] [0] via + * aggregated coefficients h4Bound and g4Bound). This is + * skipped if the wavelet boundary was already reached */ + for (; l < j + 4; l++) { + /* Smooth/approx/avg/lowpass */ + yt [h ] [i] [k] += h4 [l - j] * y [i] [yLen - 1] [k]; + /* Detail/diff/highpass */ + yt [hyt + h] [i] [k] += g4 [l - j] * y [i] [yLen - 1] [k]; + } + } + } + } + +#ifdef WAVELET_DBG + printf("%s FWD D4 (%d x %d) -> (%d x %d)\n", axis ? "VERT" : "HORIZ", + yLen, yLen, ytLen << 1, ytLen << 1 + ); + dumpCoeffs(y, yt, yLen << 1, ytLen << 1, 0); + putchar('\n'); + axis ^= 1; +#endif + } + + return ytLen; +} + + + +static unsigned padD4InvStep2 (WaveletMatrix y, WaveletMatrix yt, + unsigned yLen +) +/* + Single step of inverse 2D Daubechies D4 wavelet transform with boundary + padding on array y containng (yLen x yLen) avg and detail coefficients in + 3-tuples. 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. 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. + + The reconstructed signal is returned in the preallocated array yt. The + number of output coefficients is smaller than the input as padding + coefficients are collapsed at the boundaries during the inverse + transform. The number of output coefficients is obtained as + returned value, and can be interrogated beforehand without performing a + transform by passing either y or yt as NULL. +*/ +{ + unsigned h, i, j, k, l, + ytLen = yLen - 2; /* Max size of output array yt */ + const unsigned hy = yLen >> 1; /* Offset for detail coeffs */ +#ifdef WAVELET_DBG + static unsigned axis = 0; +#endif + + if (y && yt && yLen >= 2) { + /* Reconstruct from avg and detail coefficients */ + for (i = 0; i < yLen; i++); + for (j = 0; j < yLen - 1; j += 2) { + h = j >> 1; + + for (k = 0; k < 3; k++) { + yt [i] [j ] [k] = + h4 [2] * y [h ] [i] [k] + /* Avg */ + g4 [2] * y [hy + h ] [i] [k] + /* Detail */ + h4 [0] * y [h + 1 ] [i] [k] + /* Avg */ + g4 [0] * y [hy + h + 1] [i] [k]; /* Detail */ + + yt [i] [j + 1] [k] = + h4 [3] * y [h ] [i] [k] + /* Avg */ + g4 [3] * y [hy + h ] [i] [k] + /* Detail */ + h4 [1] * y [h + 1 ] [i] [k] + /* Avg */ + g4 [1] * y [hy + h + 1] [i] [k]; /* Detail */ + } + } + + ytLen = j; + +#ifdef WAVELET_DBG + printf("%s INV D4 (%d x %d) -> (%d x %d)\n", axis ? "VERT" : "HORIZ", + yLen, yLen, ytLen, ytLen + ); + dumpCoeffs(y, yt, yLen, ytLen, 0); + putchar('\n'); + axis ^= 1; +#endif + } + + return ytLen; +} + + + +unsigned padWaveletXform2 (WaveletMatrix y, WaveletMatrix yt, + unsigned yLen +) +/* Perform full 2D multiresolution forward wavelet transform with padded + boundary coefficients on array y containing (yLen x yLen) samples of + 3-tuples. Note no intra-tuple transform occurs. An additional array yt + of identical dimension to y is required as buffer for intermediate + results. + + The wavelet coefficients are returned in array y, containing the coarsest + average coefficients in y [0][0] followed by horizontal/vertical detail + coefficients in order of increasing resolution/frequency. + + NOTE: Both y and yt must be large enough to accommodate the additional + padding coefficients introduced at the array boundaries during the + transform. It is the caller's responsibility to preallocate and + dimension these arrays accordingly with allocWaveletMatrix(). The + maximum expected number of coefficients (which generally exceeds the + input dimension yLen!) must be interrogated beforehand by calling the + function with y or yt set to NULL. + + The returned value is the maximum number of output coefficients, or 0 if + the transform failed. +*/ +{ + /* Output size (number of coeffs, including boundary padding coeffs), + output array offset and pointers */ + unsigned numCoeffs, totalCoeffs = 0, ytOff, ytLen; + + /* Skip transform if input too short */ + if (yLen > 3) { + /* Figure out total number of coefficients (= maximum size of output + array including boundary padding). */ + for (numCoeffs = yLen, totalCoeffs = 0; + numCoeffs > 3; + totalCoeffs += numCoeffs = padD4Step2(NULL, NULL, numCoeffs) + ); + totalCoeffs += numCoeffs; + + if (y && !yt) { + /* Apply Daubechies D4 transform, accounting for boundary padding */ + +#if 0 + 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; + } +#endif + + for (numCoeffs = yLen, ytOff = totalCoeffs; + numCoeffs > 3; + numCoeffs = ytLen + ) { + /* Get num coeffs and calculate offset into output array + * (subtracting num coeffs once more to account for avg + * coeffs, which will be replaced in next iteration) */ + ytOff -= ytLen = padD4Step2(NULL, NULL, numCoeffs); + yOffPtr = y + ytOff - ytLen; + /* Resulting wavelet coefficients in y are prepended to those + * from previous iteration; consequently, returned num coeffs is + * input size for next iteration */ + d4Step(y, yt, numCoeffs, ytOffset); + /* Copy approx coeffs from y into x as input for next iteration */ + memcpy(x, yOffPtr, yLen * WAVELET_COEFFSIZE); + } + + /* Apply Haar transform in penultimate and final (since len is odd) + * iterations to decompose to a single approximation and detail + * coefficient at coarsest resolution at y [0] and y [1], + * respectively. Obviously no boundary padding required here. */ + for (; l >= 2; l--) { + haarStep(x, l, y); + /* Copy coeffs from y into x as input for next iteration */ + if (l > 2) + memcpy(x, y, l * COEFSIZE); + } + } + } + + /* All coeffs now in y; return their number */ + return numCoeffs; +} + + + +unsigned waveletInvXform ( + Coeff (*y) [3], unsigned yLen, Coeff (*x) [3], unsigned xLen0 +) +/* Perform full 2D multiresolution inverse wavelet transform on array y of + * axis dimension yLen containing wavelet coefficients as 3-tuples. Note no + * intra-tuple transform occurs (1-D only). The original signal length + * xLen0 is required for the reconstruction. + * + * The reconstructed signal is returned in the preallocated array x. This + * is shorter than yLen since the additional boundary coefficients in y are + * collapsed. Note y is used as buffer for intermediate results, and is + * therefore modified. + * + * The returned value is the reconstructed signal length, which should match + * xLen, or 0 if the transform failed. + */ +{ + unsigned xLen = 0, yOff, i, j, l, nextl, nIter; + /* The temporary output buffer xTmp is passed to d4InvStep() for + * intermediate xform steps. The additional space is required if xLen0 + * is odd, as an additional coefficient is reconstructed (which is + * subsequently dropped). If x were to be passed directly to + * d4InvStep() under these circumstances, it would overflow. */ + Coeff (*yOffPtr) [3] = NULL, xTmp [xLen0 + 1][3]; + + /* Figure out number of iterations based on original signal length */ + for (nIter = 0, l = xLen0; l > 3; nIter++, l = d4Step(NULL, l, NULL)); + + /* Invert Haar transform at two coarsest levels */ + for (l = 2; l <= 3; l++) { + haarInvStep(y, l, x); + /* Copy reconstructed signal into y for next iteration */ + memcpy(y, x, l * COEFSIZE); + } + + /* Invert Daubechies D4 transform */ + for (i = nIter; i; i--) { + /* Figure out length of current detail coeffs based on number of + * iterations */ + for ( + j = 0, yOff = yLen, l = nextl = xLen0; + j < i; + j++, nextl = l, yOff -= l = d4Step(NULL, l, NULL) + ); + +#ifdef _WAVELET_DBG + /* Zero output buffer to facilitate debugging */ + memset(xTmp, 0, sizeof(xTmp)); +#endif + + yOffPtr = y + yOff; + /* Drop excess reconstructed approximation coeffs to match number of + * detail coeffs in next iteration */ + xLen = min(d4InvStep(yOffPtr - l, l << 1, xTmp), nextl); + + /* Copy reconstructed coeffs/signal from output buffer xTmp into y for + * next iteration, or into x after last iteration */ + memcpy(i > 1 ? yOffPtr + l - xLen : x, xTmp, xLen * COEFSIZE); + } + + return xLen; +} + + + +#ifdef _WAVELET_TEST + #include + + #define ERRSTRLEN 1024 + + int main (int argc, char *argv []) + { + int i, j, xLen0, yLen, xLen, *xIdx; + Coeff (*xOrg) [3], (*x) [3], (*y) [3] = NULL; + char errstr [ERRSTRLEN]; + + if (argc < 2) { + fputs("Missing array length\n", stderr); + return -1; + } + if (!(xLen0 = atoi(argv [1]))) { + fputs("Invalid array length\n", stderr); + return -1; + } + + /* Allocate arrays for original and reconstruction (num coeffs is + * obtained from waveletXform() before allocating output array) */ + if ( + !(xOrg = calloc(xLen0, COEFSIZE)) || + !(xIdx = calloc(xLen0, sizeof(int))) || + !(x = calloc(xLen0, COEFSIZE)) || + !(y = calloc(waveletXform(NULL, xLen0, NULL), COEFSIZE)) + ) { + fputs("Failed array allocation\n", stderr); + return -1; + } + + /* Init with random data */ + srand48(111); + for (i = 0; i < xLen0; i++) { + #if 1 + for (j = 0; j < 3; j++) + xOrg [i][j] = drand48(); + #else + xOrg [i][0] = xOrg [i][1] = xOrg [i][2] = drand48(); + #endif + } + + /* Optimise signal distribution */ + waveletPresort(xOrg, xIdx, xLen0); + + #ifdef _WAVELET_DBG + /* Dump presorted signal */ + puts("-----------------------------------------------------------"); + puts("Presorted Original & Indices"); + + for (i = 0; i < xLen0; i++) + printf( + "%.3f\t%d\n", + xOrg [i][0] + xOrg [i][1] + xOrg [i][2], + xIdx [i] + ); + + puts("-----------------------------------------------------------\n"); + #endif + + /* Transform and reconstruct */ + memcpy(x, xOrg, xLen0 * COEFSIZE); + if (!(yLen = waveletXform(x, xLen0, y))) { + puts("Forward xform failed\n"); + return -1; + } + + #ifdef _WAVELET_DBG + /* Dump coefficients */ + puts("-----------------------------------------------------------\n"); + printf("%d coefficients:\n\n", yLen); + for (i = 0; i < yLen; i++) + printf("%.3f\t%.3f\t%.3f\n", y [i][0], y [i][1], y [i][2]); + puts("\n-----------------------------------------------------------\n"); + #endif + + if (!(xLen = waveletInvXform(y, yLen, x, xLen0))) { + puts("Inverse xform failed\n"); + return -1; + } + + if (xLen != xLen0) { + printf( + "Original length (%d) != reconstructed length (%d)\n", + xLen0, xLen + ); + return -1; + } + + #ifdef _WAVELET_DBG + puts("-----------------------------------------------------------\n"); + #endif + + puts("ORIG\t vs.\tINV XFORM"); + for (i = 0; i < xLen; i++) + #if 1 + printf( + "%.3f\t%.3f\t%.3f\t\t%.3f\t%.3f\t%.3f\n", + xOrg [i][0], xOrg [i][1], xOrg [i][2], + x [i][0], x [i][1], x [i][2] + ); + #else + printf("%.3f\t\t%.3f\n", xOrg [i][0], x [i][0]); + #endif + + free(xOrg); + free(xIdx); + free(x); + free(y); + return 0; + } +#endif