diff --git a/README-photonflow.txt b/README-photonflow.txt deleted file mode 100644 index e1e6ed4..0000000 --- a/README-photonflow.txt +++ /dev/null @@ -1,284 +0,0 @@ -This is the README file for the RADIANCE photon flow implementation. This -software was developed as part of the project "Three-Dimensional Light Flow", -hosted by Tokyo University of Science under the supervision of Prof. -Nozomu Yoshizawa, and supported by the JSPS Grants-in-Aid for Scientific -Research (KAKENHI) under the grant number JP19KK0115. - -$Revision: 1.3 $ $Date: 2020/11/16 22:41:09 $ - - - -INTRODUCTION - -These source files implement a modified volume photon density estimate in -RADIANCE to evaluate the illuminance of a physical lightfield using photon -flow. This code is an experimental variant of the RADIANCE photon map extension -intended for testing prior to committing it to the RADIANCE CVS as part of the -official distribution. Whether this (admittedly specialised) functionality -should be made available to the RADIANCE community as a default build option -is debatable and may be decided in due course. - -This code is essentially a C implementation of the prototype Python script -volPhotonIllum.py. As such, it is orders of magnitude faster and more -convenient, which comes to bear with large photon maps and many sensor -positions. One crucial difference between the RADIANCE and Python -implementations is that the former performs photon lookups in terms of a -fixed number of photons around the sensor position (nearest neighbour search), -while the latter does so in terms of a fixed radius (range search), and -therefore does not adapt to the local photon density, which is suboptimal. - -The illuminance evaluated from photon flow approximates Cuttle's cubic -illuminance for a virtual sensor position in a physical light field. -Photon flow constitutes volume photons deposited by mkpmap in a participating -medium (mist); due to the interaction of a photon with the medium, it will -be deposited at multiple locations along its path; the density of these -photons per path is governed by the participating medium's extinction. - -Evaluating cubic illuminance from photon flow entails performing volume -photon lookups for each of the 6 (virtual) cube faces corresponding to 3 -orthogonal axes in positive and negative directions. The resulting search -volume is a hemisphere, and photons incident from the respective cube -face (identified by it is normal or axis) are filtered on the fly based on -their stored incident directions. This results in six independent hemispherical -lookups with individual radii adapted to the density of the incident photons -for each axis. Although this incurs more overhead, it was found that this -approach yields consistently better results than a single lookup for all axes -with a posteriori filtering, particularly with asymmetric illuminance. - -It was also found that the illuminance scales with the extinction of the -participating medium, and is therefore corrected for the same. Generally, -smaller values yield more accurate results but incur longer photon distribution -times with mkpmap. - -The photon mapping code was modified to generate a variant of the existing -volume photon map, referred to as a "lightflow photon map", which corrects -the photon flux according to the extinction during it generation with mkpmap. - -When loaded with rtrace, a lightflow photon map triggers the cubic illuminance -evaluation via the function volumePhotonCubicIllum(). Unlike the existing -evaluation function for volume photons, volumePhotonDensity(), the lightflow -function does not apply the Henyey-Greenstein phase function to compute -outscattering (thus ignoring the participating medium), but instead performs -the per-axis lookups with a volume density estimate as described above. - -Because photon flow not only includes photons that have been reflected off one -or more surfaces (indirect contributions), but also photons that have just been -emitted from light sources (direct contributions), volumePhotonCubicIllum() is -is called by multiDirectPmap() in pmapsrc.c, which normally calls the direct -photon map density estimate (typically only used for debugging). To prevent -overcounting, the ambient calculation is therefore disabled in photon flow -mode (see ampPmap() in pmapamb.c). - - - -COMPILING - -The photon flow code is not compiled by default as part of the RADIANCE suite, -and must be enabled at compile time by defining PMAP_PHOTONFLOW. This may be -accelerated by only (re)compiling in ray/src/rt if a RADIANCE build already -exists: - - rmake OPT+=-DPMAP_PHOTONFLOW install - -Note that these options are independent of the underlying data structure used -to store the photon map. The default options build the legacy kd-tree. If a -larger, out-of-core photon map is desired (to accommodate >1 billion photons, -for example), this build option must be enabled by defining PMAP_OOC: - - rmake OPT+="-DPMAP_OOC -DPMAP_PHOTONFLOW" install - -For details, see the RADIANCE Out-Of-Core Photon Map Technical Report at -http://dx.doi.org/10.13140/RG.2.1.2158.9363 - -Similar to stanard (planar) photon mapping, photon flow supports kernel -density estimates by weighting individual photon contributions by their -inverse distance to the lookup point. This filter reduces bias, and can -be enabled with the following: - - rmake OPT+=-DPMAP_EPANECHNIKOV -DPMAP_OOC -DPMAP_PHOTONFLOW" install - -This command enables the 3D Epanechnikov kernel (1 - (d_i)^2), where d_i -is the normalised distance between the lookup point and photon (i). The bias -reduction is moderate, however, and ineffective with excessive lookup -bandwidths. - -Finally, the following checks verify that the photon flow code was -successfully built: - -1. 'mkpmap -help' lists the '-apV' option to generate lightflow photon maps, -2. 'rtrace -help' lists the '-aS' option to toggle between planar and - spherical illuminance evaluation from lightflow photons. - - - -USAGE - -A lightflow photon map representing the physical lightfield is generated with -mkpmap by specifying the extinction of a global participating medium either -on the command line, or in a RADIANCE scene description using a mist material -with boundary geometry (e.g. to limit photonflow to a specific region of -interest). The participating medium must not absorb nor scatter, so that -a photon's trajectory is not altered as it traverses the medium. Therefore as -a convenience, mkpmap implicitly sets the albedo and scattering eccentricity -to 1.0, which is equivalent to specifying the options -ma 1 1 1 (no absorption) -and -mg 1 (forward scattering): - - mkpmap -apV lightflow.vpm 1M [-apo port1 ... -apo portN] -apD 0.001 \ - -me .1 .1 .1 scene.oct - -The extinction -me is a free parameter here and determines the photon density -per path. Tests have shown that lower values improve accuracy but deposit -fewer photons per path, requiring mkpmap to emit more photon paths in order -to satisfy the given target photon count. - -There is a caveat with using higher values of -me, namely that the two-pass -photon distribution algorithm used by mkpmap to achieve the target number of -photons may fail. The fraction of photons emitted in the 1st pass (prepass) is -specified with -apD (predistribution factor). Its default value of 0.25 is -tuned for "well behaved" geometry and materials using global or caustic -photons. This value may be too high for volume photons in a dense participating -medium, and the photon distribution overflows in the prepass with a "photon -map overflow" error, meaning that the target photon count was already exceeded -(technically this does not yield an incorrect solution, but it's not what the -user wanted). When using photon flow, this value may need to be significantly -reduced from its default of the prepass overflows. - -Once the lightflow photon map has been generated, the RGB irradiance can -be evaluated using rtrace. This evaluation is enabled when loading the -lightflow photon map, which disables the normal ambient calculation. -Because the ambient calculation is entirely circumvented, the -ab parameter -is now irrelevant; the spherical volume photon density is evaluated immediately -at each sensor position. Each normal vector supplied with each sensor point -corresponds to the axis for each of the 6 cube faces, while the position is -(or should be!) identical for all 6 input records. Though redundant, this -parametrisation is consistent with rtrace convention. - - sensors.dat: - - <-u> - ... - - <-w> - - rtrace -h -I -ap lightflow.vpm 1000 scene.oct < sensors.dat - - Output: - - - ... - - - -In this example, the RGB irradiance E_i from photon flow is evaluated with -a bandwidth of 1000 photons for each axis i denoting one of the orthogonal -vectors {u, v, w} or its reverse. - -The cubic scalar illuminance Esr can be calculated from these values in -a separate script as a postprocess. The advantage of this discrete approach is -that it affords access to the RGB components of each axis irradiance E_i, which -would be lost if the scalar illuminance were calculated in rtrace directly. - -In addition, the mean spherical irradiance can be calculated with rtrace -to compare with the cubic scalar illuminance Esr. Rtrace accepts a boolean -option -aS to toggle between the planar irradiance evaluation for E_i as -above (corresponding to the default -aS-), and the mean spherical irradiance -evaluation using -aS+ (or just -aS): - - rtrace -h -I -ap lightflow.vpm 1000 -aS+ scene.oct < sph_sensor.dat - -Note that the normal vector specified for the sensor position is irrelevant -and ignored for this evaluation. - -As a further option, rtrace accepts an -aH option to toggle between spherical -and hemispherical lookups for lightflow photons.. In the latter case (-aH+), -photons located behind the plane are not considered, which can reduce bias -near "solid" boundaries (scene cube, geometry), behind which no photons are -deposited. However, this increases the search radius compared to a spherical -lookup, and can lead to proximity bias if the illumination gradient is high. -For lookup points primarily located in free space, it is recommended to use -spherical lookups (-aH-), which is the default. For example, the command - - rtrace -h -I -ap lightflow.vpm 1000 -aH+ scene.oct < wall_sensor.dat - -will locate lightflow photons on the near side of a solid wall. Note that -the -aH+ and -aS+ options are mutually exclusive! - -Finally, the lightflow photon distribution may be visualised by either dumping -it as a RADIANCE scene description, and rendering the photons as spheres, e.g. - - pmapdump -n 1m lightflow.vpm | objview - -or by dumping it as a point list for use in a point cloud viewer: - - pmapdump -a -f -n 1m lightflow.vpm > photonFlow.xyz - -In the latter example, pmapdump also outputs each photon's radiant flux -(in W) for further processing. See the pmapdump man page or the RADIANCE -photon map manual at http://dx.doi.org/10.13140/RG.2.1.4330.8405/1 -for details. - - - -TRANSIENT PHOTON FLOW - -An experimental variant of the light flow photon map that supports transient -(i.e. time-dependent) rendering is currently being tested. This option is -currently only available with the legacy kd-tree data structure (i.e. -PMAP_OOC is undefined at compile time), since it is straightforward to -extend to 4 dimensions. In addition, a surface-bound transient photon map -is available, which is essentially a 4D variant of the global photon map. - -Transient photons are tagged with their path length, which corresponds to -their time of flight. In order for this timestamp to be physically valid, -the user must specify the speed of light scaled to the scene geometry (i.e. -a multiple of approx. 3e8 m/s) when generating the photon map with mkpmap. -This constant is specified in addition to the photon map size with the -type-specific -apT (volume/lightflow) or -apt (surface-bound) options. -The commands - - mkpmap ... -apT lightflow.Tpm 1M 3e8 ... scene.oct - mkpmap ... -apt transient.tpm 1M 3e8 ... scene.oct - -generate ca. 1 million transient photons (lightflow and surface-bound) -with the standard speed of light in a vacuum for a scene calibrated in metres. -Transient lightflow can then be evaluated as a spatio-temporal function by -specifying a point in time in addition to a position: - - rtrace -h -I -ap lightflow.tpm 250 1.6e-9 scene.oct < sensor.dat - -In the above example, the irradiance/illuminance will be evaluated at -the given sensor position AND the given time: 1.6 nsec after photon emission -begins. - -This spatio-temporal distribution can also be visualised by passing -the parameters into rpict, which leads to surprising effects. To this end, -the Python script transpmap.py is being developed to render spatio-temporal -lightflow sequences and merge these into a video. For example, the commands - - transpmap.py -t 0 -T 8e-9 -d 1e-10 -p 10m -b 250 -n 8 cornell.vf cornell.oct - -will render transient lightflow from 10 million photons with a bandwidth of -250 photons on 8 cores within the time interval [0, 8 nsec] in increments -of 100 psec per frame. See 'transpmap.py -h' for the supported options. - - - -CONCLUSION - -Having reassessed the approach to evaluate cubic illuminance from volume -photons, revised the code accordingly, and tested it with asymmetric -illuminance distributions, I feel this release is easier to use, better -integrated into the existing RADIANCE framework, and more stable and mature. -However, the code, like the approach it implements, is still experimental and -needs some more testing in "real world" scenarious, so I welcome your feedback -in terms of its results, but also its usage. - -Once we are confident with the results, we should make the code public and -upload it to the RADIANCE CVS repository. Initially, though, it would be -prudent to still leave photon flow an option for advanced users only. - -Many thanks in advance for your cooperation and feedback. Happy photon flowing! - - ---Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) - diff --git a/Rmakefile b/Rmakefile index 9afc2cf..35846d2 100644 --- a/Rmakefile +++ b/Rmakefile @@ -1,559 +1,562 @@ # 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 pmcontribcache.o wavelet.o wavelet2.o wavelet3.o \ mrgbe.o PMRCSRC = pmapcontrib.c pmcontrib2.c pmcontrib3.c pmcontrib4.c \ pmcontribsort.c pmcontribcache.c wavelet.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 + set nonomatch; \ + rm -f $(PROGS) mrgbe-test wavelet*-test pmapcontrib-test *.o lint: $(RVSRC) $(LINT) $(LINTFLAGS) -DRVIEW $(RVSRC) $(LIBS) -# Preocmputed contrib pmap unit tests +# Precomputed contrib pmap unit tests wavelet-test: wavelet.c wavelet.h $(CC) $(CFLAGS) -DWAVELET_TEST_1D -DWAVELET_DBG \ -o wavelet-test -g wavelet.c $(MLIB) wavelet2-test: wavelet2.c wavelet2.h mrgbe.c mrgbe.h $(CC) $(CFLAGS) -DWAVELET_TEST_2D -DWAVELET_TEST_mRGBE -DWAVELET_DBG \ -o wavelet2-test -g wavelet2.c mrgbe.c $(MLIB) wavelet3-test: wavelet3.c wavelet2.c wavelet2.h mrgbe.c mrgbe.h $(CC) $(CFLAGS) -DWAVELET_TEST_2DPAD -DWAVELET_TEST_mRGBE -DWAVELET_DBG \ -o wavelet3-test -g wavelet3.c wavelet2.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 +wavelet.o: wavelet.c wavelet.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/Version.c b/Version.c index 2186455..7eded5e 100644 --- a/Version.c +++ b/Version.c @@ -1,5 +1,5 @@ /* * This file was created automatically during make. */ -char VersionID[]="RADIANCE 5.4a lastmod Sat 11 Dec 01:59:29 CET 2021 by u-no-hoo on ovo"; +char VersionID[]="RADIANCE 5.4a lastmod Thu 3 Feb 16:12:27 CET 2022 by u-no-hoo on ovo"; diff --git a/ambcomp.c b/ambcomp.c index b006118..5ed1f37 100644 --- a/ambcomp.c +++ b/ambcomp.c @@ -1,782 +1,780 @@ #ifndef lint -static const char RCSid[] = "$Id: ambcomp.c,v 2.86 2021/02/17 01:29:22 greg Exp $"; +static const char RCSid[] = "$Id: ambcomp.c,v 2.88 2021/12/15 01:38:50 greg Exp $"; #endif /* * Routines to compute "ambient" values using Monte Carlo * * Hessian calculations based on "Practical Hessian-Based Error Control * for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz * from ACM SIGGRAPH Asia 2012 conference proceedings. * * Added book-keeping optimization to avoid calculations that would * cancel due to traversal both directions on edges that are adjacent * to same-valued triangles. This cuts about half of Hessian math. * * Declarations of external symbols in ambient.h */ #include "copyright.h" #include "ray.h" #include "ambient.h" #include "random.h" #ifndef MINADIV #define MINADIV 7 /* minimum # divisions in each dimension */ #endif -extern void SDsquare2disk(double ds[2], double seedx, double seedy); - typedef struct { COLOR v; /* hemisphere sample value */ float d; /* reciprocal distance */ FVECT p; /* intersection point */ } AMBSAMP; /* sample value */ typedef struct { RAY *rp; /* originating ray sample */ int ns; /* number of samples per axis */ int sampOK; /* acquired full sample set? */ COLOR acoef; /* division contribution coefficient */ double acol[3]; /* accumulated color */ FVECT ux, uy; /* tangent axis unit vectors */ AMBSAMP sa[1]; /* sample array (extends struct) */ } AMBHEMI; /* ambient sample hemisphere */ #define AI(h,i,j) ((i)*(h)->ns + (j)) #define ambsam(h,i,j) (h)->sa[AI(h,i,j)] typedef struct { FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; double I1, I2; } FFTRI; /* vectors and coefficients for Hessian calculation */ static int ambcollision( /* proposed direciton collides? */ AMBHEMI *hp, int i, int j, FVECT dv ) { double cos_thresh; int ii, jj; /* min. spacing = 1/4th division */ cos_thresh = (PI/4.)/(double)hp->ns; cos_thresh = 1. - .5*cos_thresh*cos_thresh; /* check existing neighbors */ for (ii = i-1; ii <= i+1; ii++) { if (ii < 0) continue; if (ii >= hp->ns) break; for (jj = j-1; jj <= j+1; jj++) { AMBSAMP *ap; FVECT avec; double dprod; if (jj < 0) continue; if (jj >= hp->ns) break; if ((ii==i) & (jj==j)) continue; ap = &ambsam(hp,ii,jj); if (ap->d <= .5/FHUGE) continue; /* no one home */ VSUB(avec, ap->p, hp->rp->rop); dprod = DOT(avec, dv); if (dprod >= cos_thresh*VLEN(avec)) return(1); /* collision */ } } return(0); /* nothing to worry about */ } static int ambsample( /* initial ambient division sample */ AMBHEMI *hp, int i, int j, int n ) { AMBSAMP *ap = &ambsam(hp,i,j); RAY ar; int hlist[3], ii; - double spt[2], zd; + RREAL spt[2]; + double zd; /* generate hemispherical sample */ /* ambient coefficient for weight */ if (ambacc > FTINY) setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); else copycolor(ar.rcoef, hp->acoef); if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) return(0); if (ambacc > FTINY) { multcolor(ar.rcoef, hp->acoef); scalecolor(ar.rcoef, 1./AVGREFL); } hlist[0] = hp->rp->rno; hlist[1] = j; hlist[2] = i; multisamp(spt, 2, urand(ilhash(hlist,3)+n)); resample: - SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); + square2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); for (ii = 3; ii--; ) ar.rdir[ii] = spt[0]*hp->ux[ii] + spt[1]*hp->uy[ii] + zd*hp->rp->ron[ii]; checknorm(ar.rdir); /* avoid coincident samples */ if (!n && ambcollision(hp, i, j, ar.rdir)) { spt[0] = frandom(); spt[1] = frandom(); goto resample; /* reject this sample */ } dimlist[ndims++] = AI(hp,i,j) + 90171; rayvalue(&ar); /* evaluate ray */ ndims--; zd = raydistance(&ar); if (zd <= FTINY) return(0); /* should never happen */ multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ if (zd*ap->d < 1.0) /* new/closer distance? */ ap->d = 1.0/zd; if (!n) { /* record first vertex & value */ if (zd > 10.0*thescene.cusize + 1000.) zd = 10.0*thescene.cusize + 1000.; VSUM(ap->p, ar.rorg, ar.rdir, zd); copycolor(ap->v, ar.rcol); } else { /* else update recorded value */ hp->acol[RED] -= colval(ap->v,RED); hp->acol[GRN] -= colval(ap->v,GRN); hp->acol[BLU] -= colval(ap->v,BLU); zd = 1.0/(double)(n+1); scalecolor(ar.rcol, zd); zd *= (double)n; scalecolor(ap->v, zd); addcolor(ap->v, ar.rcol); } addcolor(hp->acol, ap->v); /* add to our sum */ return(1); } /* Estimate variance based on ambient division differences */ static float * getambdiffs(AMBHEMI *hp) { const double normf = 1./bright(hp->acoef); float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); float *ep; AMBSAMP *ap; double b, b1, d2; int i, j; if (earr == NULL) /* out of memory? */ return(NULL); /* sum squared neighbor diffs */ for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++) for (j = 0; j < hp->ns; j++, ap++, ep++) { b = bright(ap[0].v); if (i) { /* from above */ b1 = bright(ap[-hp->ns].v); d2 = b - b1; d2 *= d2*normf/(b + b1); ep[0] += d2; ep[-hp->ns] += d2; } if (!j) continue; /* from behind */ b1 = bright(ap[-1].v); d2 = b - b1; d2 *= d2*normf/(b + b1); ep[0] += d2; ep[-1] += d2; if (!i) continue; /* diagonal */ b1 = bright(ap[-hp->ns-1].v); d2 = b - b1; d2 *= d2*normf/(b + b1); ep[0] += d2; ep[-hp->ns-1] += d2; } /* correct for number of neighbors */ earr[0] *= 8./3.; earr[hp->ns-1] *= 8./3.; earr[(hp->ns-1)*hp->ns] *= 8./3.; earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.; for (i = 1; i < hp->ns-1; i++) { earr[i*hp->ns] *= 8./5.; earr[i*hp->ns + hp->ns-1] *= 8./5.; } for (j = 1; j < hp->ns-1; j++) { earr[j] *= 8./5.; earr[(hp->ns-1)*hp->ns + j] *= 8./5.; } return(earr); } /* Perform super-sampling on hemisphere (introduces bias) */ static void ambsupersamp(AMBHEMI *hp, int cnt) { float *earr = getambdiffs(hp); double e2rem = 0; float *ep; int i, j, n, nss; if (earr == NULL) /* just skip calc. if no memory */ return; /* accumulate estimated variances */ for (ep = earr + hp->ns*hp->ns; ep > earr; ) e2rem += *--ep; ep = earr; /* perform super-sampling */ for (i = 0; i < hp->ns; i++) for (j = 0; j < hp->ns; j++) { if (e2rem <= FTINY) goto done; /* nothing left to do */ nss = *ep/e2rem*cnt + frandom(); for (n = 1; n <= nss && ambsample(hp,i,j,n); n++) if (!--cnt) goto done; e2rem -= *ep++; /* update remainder */ } done: free(earr); } static AMBHEMI * samp_hemi( /* sample indirect hemisphere */ COLOR rcol, RAY *r, double wt ) { AMBHEMI *hp; double d; int n, i, j; /* insignificance check */ if (bright(rcol) <= FTINY) return(NULL); /* set number of divisions */ if (ambacc <= FTINY && wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight))) wt = d; /* avoid ray termination */ n = sqrt(ambdiv * wt) + 0.5; i = 1 + (MINADIV-1)*(ambacc > FTINY); if (n < i) /* use minimum number of samples? */ n = i; /* allocate sampling array */ hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); if (hp == NULL) error(SYSTEM, "out of memory in samp_hemi"); hp->rp = r; hp->ns = n; hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0; memset(hp->sa, 0, sizeof(AMBSAMP)*n*n); hp->sampOK = 0; /* assign coefficient */ copycolor(hp->acoef, rcol); d = 1.0/(n*n); scalecolor(hp->acoef, d); /* make tangent plane axes */ if (!getperpendicular(hp->ux, r->ron, 1)) error(CONSISTENCY, "bad ray direction in samp_hemi"); VCROSS(hp->uy, r->ron, hp->ux); /* sample divisions */ for (i = hp->ns; i--; ) for (j = hp->ns; j--; ) hp->sampOK += ambsample(hp, i, j, 0); copycolor(rcol, hp->acol); if (!hp->sampOK) { /* utter failure? */ free(hp); return(NULL); } if (hp->sampOK < hp->ns*hp->ns) { hp->sampOK *= -1; /* soft failure */ return(hp); } if (hp->sampOK <= MINADIV*MINADIV) return(hp); /* don't bother super-sampling */ n = ambssamp*wt + 0.5; if (n > 8) { /* perform super-sampling? */ ambsupersamp(hp, n); copycolor(rcol, hp->acol); } return(hp); /* all is well */ } /* Return brightness of farthest ambient sample */ static double back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3) { if (hp->sa[n1].d <= hp->sa[n2].d) { if (hp->sa[n1].d <= hp->sa[n3].d) return(colval(hp->sa[n1].v,CIEY)); return(colval(hp->sa[n3].v,CIEY)); } if (hp->sa[n2].d <= hp->sa[n3].d) return(colval(hp->sa[n2].v,CIEY)); return(colval(hp->sa[n3].v,CIEY)); } /* Compute vectors and coefficients for Hessian/gradient calcs */ static void comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1) { double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; int ii; VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop); VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop); VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p); VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); dot_e = DOT(ftp->e_i,ftp->e_i); dot_er = DOT(ftp->e_i, ftp->r_i); rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * sqrt( rdot_cp ); ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + dot_e*ftp->I1 )*0.5*rdot_cp; J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; for (ii = 3; ii--; ) ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii]; } /* Compose 3x3 matrix from two vectors */ static void compose_matrix(FVECT mat[3], FVECT va, FVECT vb) { mat[0][0] = 2.0*va[0]*vb[0]; mat[1][1] = 2.0*va[1]*vb[1]; mat[2][2] = 2.0*va[2]*vb[2]; mat[0][1] = mat[1][0] = va[0]*vb[1] + va[1]*vb[0]; mat[0][2] = mat[2][0] = va[0]*vb[2] + va[2]*vb[0]; mat[1][2] = mat[2][1] = va[1]*vb[2] + va[2]*vb[1]; } /* Compute partial 3x3 Hessian matrix for edge */ static void comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) { FVECT ncp; FVECT m1[3], m2[3], m3[3], m4[3]; double d1, d2, d3, d4; double I3, J3, K3; int i, j; /* compute intermediate coefficients */ d1 = 1.0/DOT(ftp->r_i,ftp->r_i); d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); d3 = 1.0/DOT(ftp->e_i,ftp->e_i); d4 = DOT(ftp->e_i, ftp->r_i); I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 ) / ( 4.0*DOT(ftp->rcp,ftp->rcp) ); J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); /* intermediate matrices */ VCROSS(ncp, nrm, ftp->e_i); compose_matrix(m1, ncp, ftp->rI2_eJ2); compose_matrix(m2, ftp->r_i, ftp->r_i); compose_matrix(m3, ftp->e_i, ftp->e_i); compose_matrix(m4, ftp->r_i, ftp->e_i); d1 = DOT(nrm, ftp->rcp); d2 = -d1*ftp->I2; d1 *= 2.0; for (i = 3; i--; ) /* final matrix sum */ for (j = 3; j--; ) { hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] + 2.0*J3*m4[i][j] ); hess[i][j] += d2*(i==j); hess[i][j] *= -1.0/PI; } } /* Reverse hessian calculation result for edge in other direction */ static void rev_hessian(FVECT hess[3]) { int i; for (i = 3; i--; ) { hess[i][0] = -hess[i][0]; hess[i][1] = -hess[i][1]; hess[i][2] = -hess[i][2]; } } /* Add to radiometric Hessian from the given triangle */ static void add2hessian(FVECT hess[3], FVECT ehess1[3], FVECT ehess2[3], FVECT ehess3[3], double v) { int i, j; for (i = 3; i--; ) for (j = 3; j--; ) hess[i][j] += v*( ehess1[i][j] + ehess2[i][j] + ehess3[i][j] ); } /* Compute partial displacement form factor gradient for edge */ static void comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm) { FVECT ncp; double f1; int i; f1 = 2.0*DOT(nrm, ftp->rcp); VCROSS(ncp, nrm, ftp->e_i); for (i = 3; i--; ) grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); } /* Reverse gradient calculation result for edge in other direction */ static void rev_gradient(FVECT grad) { grad[0] = -grad[0]; grad[1] = -grad[1]; grad[2] = -grad[2]; } /* Add to displacement gradient from the given triangle */ static void add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v) { int i; for (i = 3; i--; ) grad[i] += v*( egrad1[i] + egrad2[i] + egrad3[i] ); } /* Compute anisotropic radii and eigenvector directions */ static void eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3]) { double hess2[2][2]; FVECT a, b; double evalue[2], slope1, xmag1; int i; /* project Hessian to sample plane */ for (i = 3; i--; ) { a[i] = DOT(hessian[i], uv[0]); b[i] = DOT(hessian[i], uv[1]); } hess2[0][0] = DOT(uv[0], a); hess2[0][1] = DOT(uv[0], b); hess2[1][0] = DOT(uv[1], a); hess2[1][1] = DOT(uv[1], b); /* compute eigenvalue(s) */ i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]); if (i == 1) /* double-root (circle) */ evalue[1] = evalue[0]; if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) | ((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) { ra[0] = ra[1] = maxarad; return; } if (evalue[0] > evalue[1]) { ra[0] = sqrt(sqrt(4.0/evalue[0])); ra[1] = sqrt(sqrt(4.0/evalue[1])); slope1 = evalue[1]; } else { ra[0] = sqrt(sqrt(4.0/evalue[1])); ra[1] = sqrt(sqrt(4.0/evalue[0])); slope1 = evalue[0]; } /* compute unit eigenvectors */ if (fabs(hess2[0][1]) <= FTINY) return; /* uv OK as is */ slope1 = (slope1 - hess2[0][0]) / hess2[0][1]; xmag1 = sqrt(1.0/(1.0 + slope1*slope1)); for (i = 3; i--; ) { b[i] = xmag1*uv[0][i] + slope1*xmag1*uv[1][i]; a[i] = slope1*xmag1*uv[0][i] - xmag1*uv[1][i]; } VCOPY(uv[0], a); VCOPY(uv[1], b); } static void ambHessian( /* anisotropic radii & pos. gradient */ AMBHEMI *hp, FVECT uv[2], /* returned */ float ra[2], /* returned (optional) */ float pg[2] /* returned (optional) */ ) { static char memerrmsg[] = "out of memory in ambHessian()"; FVECT (*hessrow)[3] = NULL; FVECT *gradrow = NULL; FVECT hessian[3]; FVECT gradient; FFTRI fftr; int i, j; /* be sure to assign unit vectors */ VCOPY(uv[0], hp->ux); VCOPY(uv[1], hp->uy); /* clock-wise vertex traversal from sample POV */ if (ra != NULL) { /* initialize Hessian row buffer */ hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*(hp->ns-1)); if (hessrow == NULL) error(SYSTEM, memerrmsg); memset(hessian, 0, sizeof(hessian)); } else if (pg == NULL) /* bogus call? */ return; if (pg != NULL) { /* initialize form factor row buffer */ gradrow = (FVECT *)malloc(sizeof(FVECT)*(hp->ns-1)); if (gradrow == NULL) error(SYSTEM, memerrmsg); memset(gradient, 0, sizeof(gradient)); } /* compute first row of edges */ for (j = 0; j < hp->ns-1; j++) { comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1)); if (hessrow != NULL) comp_hessian(hessrow[j], &fftr, hp->rp->ron); if (gradrow != NULL) comp_gradient(gradrow[j], &fftr, hp->rp->ron); } /* sum each row of triangles */ for (i = 0; i < hp->ns-1; i++) { FVECT hesscol[3]; /* compute first vertical edge */ FVECT gradcol; comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0)); if (hessrow != NULL) comp_hessian(hesscol, &fftr, hp->rp->ron); if (gradrow != NULL) comp_gradient(gradcol, &fftr, hp->rp->ron); for (j = 0; j < hp->ns-1; j++) { FVECT hessdia[3]; /* compute triangle contributions */ FVECT graddia; double backg; backg = back_ambval(hp, AI(hp,i,j), AI(hp,i,j+1), AI(hp,i+1,j)); /* diagonal (inner) edge */ comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j)); if (hessrow != NULL) { comp_hessian(hessdia, &fftr, hp->rp->ron); rev_hessian(hesscol); add2hessian(hessian, hessrow[j], hessdia, hesscol, backg); } if (gradrow != NULL) { comp_gradient(graddia, &fftr, hp->rp->ron); rev_gradient(gradcol); add2gradient(gradient, gradrow[j], graddia, gradcol, backg); } /* initialize edge in next row */ comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j)); if (hessrow != NULL) comp_hessian(hessrow[j], &fftr, hp->rp->ron); if (gradrow != NULL) comp_gradient(gradrow[j], &fftr, hp->rp->ron); /* new column edge & paired triangle */ backg = back_ambval(hp, AI(hp,i+1,j+1), AI(hp,i+1,j), AI(hp,i,j+1)); comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1)); if (hessrow != NULL) { comp_hessian(hesscol, &fftr, hp->rp->ron); rev_hessian(hessdia); add2hessian(hessian, hessrow[j], hessdia, hesscol, backg); if (i < hp->ns-2) rev_hessian(hessrow[j]); } if (gradrow != NULL) { comp_gradient(gradcol, &fftr, hp->rp->ron); rev_gradient(graddia); add2gradient(gradient, gradrow[j], graddia, gradcol, backg); if (i < hp->ns-2) rev_gradient(gradrow[j]); } } } /* release row buffers */ if (hessrow != NULL) free(hessrow); if (gradrow != NULL) free(gradrow); if (ra != NULL) /* extract eigenvectors & radii */ eigenvectors(uv, ra, hessian); if (pg != NULL) { /* tangential position gradient */ pg[0] = DOT(gradient, uv[0]); pg[1] = DOT(gradient, uv[1]); } } /* Compute direction gradient from a hemispherical sampling */ static void ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) { AMBSAMP *ap; double dgsum[2]; int n; FVECT vd; double gfact; dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */ for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { /* use vector for azimuth + 90deg */ VSUB(vd, ap->p, hp->rp->rop); /* brightness over cosine factor */ gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); /* sine = proj_radius/vd_length */ dgsum[0] -= DOT(uv[1], vd) * gfact; dgsum[1] += DOT(uv[0], vd) * gfact; } dg[0] = dgsum[0] / (hp->ns*hp->ns); dg[1] = dgsum[1] / (hp->ns*hp->ns); } /* Compute potential light leak direction flags for cache value */ static uint32 ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1) { const double max_d = 1.0/(minarad*ambacc + 0.001); const double ang_res = 0.5*PI/hp->ns; const double ang_step = ang_res/((int)(16/PI*ang_res) + 1.01); double avg_d = 0; uint32 flgs = 0; FVECT vec; double u, v; double ang, a1; int i, j; /* don't bother for a few samples */ if (hp->ns < 8) return(0); /* check distances overhead */ for (i = hp->ns*3/4; i-- > hp->ns>>2; ) for (j = hp->ns*3/4; j-- > hp->ns>>2; ) avg_d += ambsam(hp,i,j).d; avg_d *= 4.0/(hp->ns*hp->ns); if (avg_d*r0 >= 1.0) /* ceiling too low for corral? */ return(0); if (avg_d >= max_d) /* insurance */ return(0); /* else circle around perimeter */ for (i = 0; i < hp->ns; i++) for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) { AMBSAMP *ap = &ambsam(hp,i,j); if ((ap->d <= FTINY) | (ap->d >= max_d)) continue; /* too far or too near */ VSUB(vec, ap->p, hp->rp->rop); u = DOT(vec, uv[0]); v = DOT(vec, uv[1]); if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v) continue; /* occluder outside ellipse */ ang = atan2a(v, u); /* else set direction flags */ for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step) flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0))); } return(flgs); } int doambient( /* compute ambient component */ COLOR rcol, /* input/output color */ RAY *r, double wt, FVECT uv[2], /* returned (optional) */ float ra[2], /* returned (optional) */ float pg[2], /* returned (optional) */ float dg[2], /* returned (optional) */ uint32 *crlp /* returned (optional) */ ) { AMBHEMI *hp = samp_hemi(rcol, r, wt); FVECT my_uv[2]; double d, K; AMBSAMP *ap; int i; /* clear return values */ if (uv != NULL) memset(uv, 0, sizeof(FVECT)*2); if (ra != NULL) ra[0] = ra[1] = 0.0; if (pg != NULL) pg[0] = pg[1] = 0.0; if (dg != NULL) dg[0] = dg[1] = 0.0; if (crlp != NULL) *crlp = 0; if (hp == NULL) /* sampling falure? */ return(0); if ((ra == NULL) & (pg == NULL) & (dg == NULL) || (hp->sampOK < 0) | (hp->ns < MINADIV)) { free(hp); /* Hessian not requested/possible */ return(-1); /* value-only return value */ } if ((d = bright(rcol)) > FTINY) { /* normalize Y values */ d = 0.99*(hp->ns*hp->ns)/d; K = 0.01; } else { /* or fall back on geometric Hessian */ K = 1.0; pg = NULL; dg = NULL; crlp = NULL; } ap = hp->sa; /* relative Y channel from here on... */ for (i = hp->ns*hp->ns; i--; ap++) colval(ap->v,CIEY) = bright(ap->v)*d + K; if (uv == NULL) /* make sure we have axis pointers */ uv = my_uv; /* compute radii & pos. gradient */ ambHessian(hp, uv, ra, pg); if (dg != NULL) /* compute direction gradient */ ambdirgrad(hp, uv, dg); if (ra != NULL) { /* scale/clamp radii */ if (pg != NULL) { if (ra[0]*(d = fabs(pg[0])) > 1.0) ra[0] = 1.0/d; if (ra[1]*(d = fabs(pg[1])) > 1.0) ra[1] = 1.0/d; if (ra[0] > ra[1]) ra[0] = ra[1]; } if (ra[0] < minarad) { ra[0] = minarad; if (ra[1] < minarad) ra[1] = minarad; } ra[0] *= d = 1.0/sqrt(wt); if ((ra[1] *= d) > 2.0*ra[0]) ra[1] = 2.0*ra[0]; if (ra[1] > maxarad) { ra[1] = maxarad; if (ra[0] > maxarad) ra[0] = maxarad; } /* flag encroached directions */ - if (crlp != NULL) + if (crlp != NULL) /* XXX doesn't update with changes to ambacc */ *crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc); if (pg != NULL) { /* cap gradient if necessary */ d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; if (d > 1.0) { d = 1.0/sqrt(d); pg[0] *= d; pg[1] *= d; } } } free(hp); /* clean up and return */ return(1); } - diff --git a/ambient.c b/ambient.c index a5f8e15..d091963 100644 --- a/ambient.c +++ b/ambient.c @@ -1,1050 +1,1054 @@ -static const char RCSid[] = "$Id: ambient.c,v 2.110 2021/02/19 22:05:46 greg Exp $"; +static const char RCSid[] = "$Id: ambient.c,v 2.111 2021/11/18 19:38:21 greg Exp $"; /* * ambient.c - routines dealing with ambient (inter-reflected) component. * * Declarations of external symbols in ambient.h */ #include "copyright.h" #include #include "platform.h" #include "ray.h" #include "otypes.h" #include "otspecial.h" #include "resolu.h" #include "ambient.h" #include "random.h" #include "pmapamb.h" #ifndef OCTSCALE #define OCTSCALE 1.0 /* ceil((valid rad.)/(cube size)) */ #endif extern char *shm_boundary; /* memory sharing boundary */ #ifndef MAXASET #define MAXASET 4095 /* maximum number of elements in ambient set */ #endif OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */ double maxarad; /* maximum ambient radius */ double minarad; /* minimum ambient radius */ static AMBTREE atrunk; /* our ambient trunk node */ static FILE *ambfp = NULL; /* ambient file pointer */ static int nunflshed = 0; /* number of unflushed ambient values */ #ifndef SORT_THRESH #ifdef SMLMEM #define SORT_THRESH ((16L<<20)/sizeof(AMBVAL)) #else #define SORT_THRESH ((64L<<20)/sizeof(AMBVAL)) #endif #endif #ifndef SORT_INTVL #define SORT_INTVL (SORT_THRESH<<1) #endif #ifndef MAX_SORT_INTVL #define MAX_SORT_INTVL (SORT_INTVL<<6) #endif static double avsum = 0.; /* computed ambient value sum (log) */ static unsigned int navsum = 0; /* number of values in avsum */ static unsigned int nambvals = 0; /* total number of indirect values */ static unsigned int nambshare = 0; /* number of values from file */ static unsigned long ambclock = 0; /* ambient access clock */ static unsigned long lastsort = 0; /* time of last value sort */ static long sortintvl = SORT_INTVL; /* time until next sort */ static FILE *ambinp = NULL; /* auxiliary file for input */ static long lastpos = -1; /* last flush position */ #define MAXACLOCK (1L<<30) /* clock turnover value */ /* * Track access times unless we are sharing ambient values * through memory on a multiprocessor, when we want to avoid * claiming our own memory (copy on write). Go ahead anyway * if more than two thirds of our values are unshared. * Compile with -Dtracktime=0 to turn this code off. */ #ifndef tracktime #define tracktime (shm_boundary == NULL || nambvals > 3*nambshare) #endif #define AMBFLUSH (BUFSIZ/AMBVALSIZ) #define newambval() (AMBVAL *)malloc(sizeof(AMBVAL)) #define tfunc(lwr, x, upr) (((x)-(lwr))/((upr)-(lwr))) static void initambfile(int creat); static void avsave(AMBVAL *av); static AMBVAL *avstore(AMBVAL *aval); static AMBTREE *newambtree(void); static void freeambtree(AMBTREE *atp); typedef void unloadtf_t(AMBVAL *); static unloadtf_t avinsert; static unloadtf_t av2list; static unloadtf_t avfree; static void unloadatree(AMBTREE *at, unloadtf_t *f); static int aposcmp(const void *avp1, const void *avp2); static int avlmemi(AMBVAL *avaddr); static void sortambvals(int always); static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang); static double sumambient(COLOR acol, RAY *r, FVECT rn, int al, AMBTREE *at, FVECT c0, double s); static int makeambient(COLOR acol, RAY *r, FVECT rn, int al); static int extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv, FVECT uvw[3]); #ifdef F_SETLKW static void aflock(int typ); #endif void setambres( /* set ambient resolution */ int ar ) { ambres = ar < 0 ? 0 : ar; /* may be done already */ /* set min & max radii */ if (ar <= 0) { minarad = 0; maxarad = thescene.cusize*0.2; } else { minarad = thescene.cusize / ar; maxarad = 64.0 * minarad; /* heuristic */ if (maxarad > thescene.cusize*0.2) maxarad = thescene.cusize*0.2; } if (minarad <= FTINY) minarad = 10.0*FTINY; if (maxarad <= minarad) maxarad = 64.0 * minarad; } void setambacc( /* set ambient accuracy */ double newa ) { static double olda; /* remember previous setting here */ newa *= (newa > 0); if (fabs(newa - olda) >= .05*(newa + olda)) { ambacc = newa; if (ambacc > FTINY && nambvals > 0) sortambvals(1); /* rebuild tree */ } } void setambient(void) /* initialize calculation */ { int readonly = 0; long flen; AMBVAL amb; /* make sure we're fresh */ ambdone(); /* init ambient limits */ setambres(ambres); setambacc(ambacc); if (ambfile == NULL || !ambfile[0]) return; if (ambacc <= FTINY) { sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened", ambfile); error(WARNING, errmsg); return; } /* open ambient file */ if ((ambfp = fopen(ambfile, "r+")) == NULL) readonly = (ambfp = fopen(ambfile, "r")) != NULL; if (ambfp != NULL) { initambfile(0); /* file exists */ lastpos = ftell(ambfp); while (readambval(&amb, ambfp)) avstore(&amb); nambshare = nambvals; /* share loaded values */ if (readonly) { sprintf(errmsg, "loaded %u values from read-only ambient file", nambvals); error(WARNING, errmsg); fclose(ambfp); /* close file so no writes */ ambfp = NULL; return; /* avoid ambsync() */ } /* align file pointer */ lastpos += (long)nambvals*AMBVALSIZ; flen = lseek(fileno(ambfp), (off_t)0, SEEK_END); if (flen != lastpos) { sprintf(errmsg, "ignoring last %ld values in ambient file (corrupted)", (flen - lastpos)/AMBVALSIZ); error(WARNING, errmsg); fseek(ambfp, lastpos, SEEK_SET); ftruncate(fileno(ambfp), (off_t)lastpos); } } else if ((ambfp = fopen(ambfile, "w+")) != NULL) { initambfile(1); /* else create new file */ fflush(ambfp); lastpos = ftell(ambfp); } else { sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile); error(SYSTEM, errmsg); } #ifdef F_SETLKW aflock(F_UNLCK); /* release file */ #endif } void ambdone(void) /* close ambient file and free memory */ { if (ambfp != NULL) { /* close ambient file */ ambsync(); fclose(ambfp); ambfp = NULL; if (ambinp != NULL) { fclose(ambinp); ambinp = NULL; } lastpos = -1; } /* free ambient tree */ unloadatree(&atrunk, avfree); /* reset state variables */ avsum = 0.; navsum = 0; nambvals = 0; nambshare = 0; ambclock = 0; lastsort = 0; sortintvl = SORT_INTVL; } void ambnotify( /* record new modifier */ OBJECT obj ) { static int hitlimit = 0; OBJREC *o; char **amblp; if (obj == OVOID) { /* starting over */ ambset[0] = 0; hitlimit = 0; return; } o = objptr(obj); if (hitlimit || !ismodifier(o->otype)) return; for (amblp = amblist; *amblp != NULL; amblp++) if (!strcmp(o->oname, *amblp)) { if (ambset[0] >= MAXASET) { error(WARNING, "too many modifiers in ambient list"); hitlimit++; return; /* should this be fatal? */ } insertelem(ambset, obj); return; } } void multambient( /* compute ambient component & multiply by coef. */ COLOR aval, RAY *r, FVECT nrm ) { + static double logAvgAbsorp = 1; static int rdepth = 0; /* ambient recursion */ COLOR acol, caustic; int i, ok; double d, l; /* PMAP: Factor in ambient from photon map, if enabled and ray is * ambient. Return as all ambient components accounted for, else * continue. */ if (ambPmap(aval, r, rdepth)) return; + if (logAvgAbsorp > 0) /* exclude in -aw to avoid growth */ + logAvgAbsorp = log(1.-AVGREFL); + /* PMAP: Factor in specular-diffuse ambient (caustics) from photon * map, if enabled and ray is primary, else caustic is zero. Continue * with RADIANCE ambient calculation */ copycolor(caustic, aval); ambPmapCaustic(caustic, r, rdepth); if (ambdiv <= 0) /* no ambient calculation */ goto dumbamb; /* check number of bounces */ if (rdepth >= ambounce) goto dumbamb; /* check ambient list */ if (ambincl != -1 && r->ro != NULL && ambincl != inset(ambset, r->ro->omod)) goto dumbamb; if (ambacc <= FTINY) { /* no ambient storage */ FVECT uvd[2]; float dgrad[2], *dgp = NULL; if (nrm != r->ron && DOT(nrm,r->ron) < 0.9999) dgp = dgrad; /* compute rotational grad. */ copycolor(acol, aval); rdepth++; ok = doambient(acol, r, r->rweight, uvd, NULL, NULL, dgp, NULL); rdepth--; if (!ok) goto dumbamb; if ((ok > 0) & (dgp != NULL)) { /* apply texture */ FVECT v1; VCROSS(v1, r->ron, nrm); d = 1.0; for (i = 3; i--; ) d += v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]); if (d >= 0.05) scalecolor(acol, d); } copycolor(aval, acol); /* PMAP: add in caustic */ addcolor(aval, caustic); return; } if (tracktime) /* sort to minimize thrashing */ sortambvals(0); /* interpolate ambient value */ setcolor(acol, 0.0, 0.0, 0.0); d = sumambient(acol, r, nrm, rdepth, &atrunk, thescene.cuorg, thescene.cusize); if (d > FTINY) { d = 1.0/d; scalecolor(acol, d); multcolor(aval, acol); /* PMAP: add in caustic */ addcolor(aval, caustic); return; } rdepth++; /* need to cache new value */ ok = makeambient(acol, r, nrm, rdepth-1); rdepth--; if (ok) { multcolor(aval, acol); /* computed new value */ /* PMAP: add in caustic */ addcolor(aval, caustic); return; } dumbamb: /* return global value */ if ((ambvwt <= 0) | (navsum == 0)) { multcolor(aval, ambval); /* PMAP: add in caustic */ addcolor(aval, caustic); return; } l = bright(ambval); /* average in computations */ if (l > FTINY) { - d = (log(l)*(double)ambvwt + avsum) / + d = (log(l)*(double)ambvwt + avsum + logAvgAbsorp*navsum) / (double)(ambvwt + navsum); d = exp(d) / l; scalecolor(aval, d); multcolor(aval, ambval); /* apply color of ambval */ } else { - d = exp( avsum / (double)navsum ); + d = exp( avsum/(double)navsum + logAvgAbsorp ); scalecolor(aval, d); /* neutral color */ } } /* Plug a potential leak where ambient cache value is occluded */ static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang) { const double cost70sq = 0.1169778; /* cos(70deg)^2 */ RAY rtst; FVECT vdif; double normdot, ndotd, nadotd; double a, b, c, t[2]; ang += 2.*PI*(ang < 0); /* check direction flags */ if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) ) return(0); /* * Generate test ray, targeting 20 degrees above sample point plane * along surface normal from cache position. This should be high * enough to miss local geometry we don't really care about. */ VSUB(vdif, ap->pos, r->rop); normdot = DOT(anorm, r->ron); ndotd = DOT(vdif, r->ron); nadotd = DOT(vdif, anorm); a = normdot*normdot - cost70sq; b = 2.0*(normdot*ndotd - nadotd*cost70sq); c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq; if (quadratic(t, a, b, c) != 2) return(1); /* should rarely happen */ if (t[1] <= FTINY) return(0); /* should fail behind test */ rayorigin(&rtst, SHADOW, r, NULL); VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */ rtst.rmax = normalize(rtst.rdir); /* short ray test */ while (localhit(&rtst, &thescene)) { /* check for occluder */ OBJREC *m = findmaterial(rtst.ro); if (m != NULL && !istransp(m->otype) && !isBSDFproxy(m) && (rtst.clipset == NULL || !inset(rtst.clipset, rtst.ro->omod))) return(1); /* plug light leak */ VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */ rtst.rmax -= rtst.rot; rayclear(&rtst); } return(0); /* seems we're OK */ } static double sumambient( /* get interpolated ambient value */ COLOR acol, RAY *r, FVECT rn, int al, AMBTREE *at, FVECT c0, double s ) { /* initial limit is 10 degrees plus ambacc radians */ const double minangle = 10.0 * PI/180.; double maxangle = minangle + ambacc; double wsum = 0.0; FVECT ck0; int i, j; AMBVAL *av; if (at->kid != NULL) { /* sum children first */ s *= 0.5; for (i = 0; i < 8; i++) { for (j = 0; j < 3; j++) { ck0[j] = c0[j]; if (1<rop[j] < ck0[j] - OCTSCALE*s) break; if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s) break; } if (j == 3) wsum += sumambient(acol, r, rn, al, at->kid+i, ck0, s); } /* good enough? */ if (wsum >= 0.05 && s > minarad*10.0) return(wsum); } /* adjust maximum angle */ if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6)) maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.; /* sum this node */ for (av = at->alist; av != NULL; av = av->next) { double u, v, d, delta_r2, delta_t2; COLOR ct; FVECT uvw[3]; /* record access */ if (tracktime) av->latick = ambclock; /* * Ambient level test */ if (av->lvl > al || /* list sorted, so this works */ (av->lvl == al) & (av->weight < 0.9*r->rweight)) break; /* * Direction test using unperturbed normal */ decodedir(uvw[2], av->ndir); d = DOT(uvw[2], r->ron); if (d <= 0.0) /* >= 90 degrees */ continue; delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */ if (delta_r2 >= maxangle*maxangle) continue; /* * Modified ray behind test */ VSUB(ck0, r->rop, av->pos); d = DOT(ck0, uvw[2]); if (d < -minarad*ambacc-.001) continue; d /= av->rad[0]; delta_t2 = d*d; if (delta_t2 >= ambacc*ambacc) continue; /* * Elliptical radii test based on Hessian */ decodedir(uvw[0], av->udir); VCROSS(uvw[1], uvw[2], uvw[0]); d = (u = DOT(ck0, uvw[0])) / av->rad[0]; delta_t2 += d*d; d = (v = DOT(ck0, uvw[1])) / av->rad[1]; delta_t2 += d*d; if (delta_t2 >= ambacc*ambacc) continue; /* * Test for potential light leak */ if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u))) continue; /* * Extrapolate value and compute final weight (hat function) */ if (!extambient(ct, av, r->rop, rn, uvw)) continue; d = tfunc(maxangle, sqrt(delta_r2), 0.0) * tfunc(ambacc, sqrt(delta_t2), 0.0); scalecolor(ct, d); addcolor(acol, ct); wsum += d; } return(wsum); } static int makeambient( /* make a new ambient value for storage */ COLOR acol, RAY *r, FVECT rn, int al ) { AMBVAL amb; FVECT uvw[3]; int i; amb.weight = 1.0; /* compute weight */ for (i = al; i-- > 0; ) amb.weight *= AVGREFL; if (r->rweight < 0.1*amb.weight) /* heuristic override */ amb.weight = 1.25*r->rweight; setcolor(acol, AVGREFL, AVGREFL, AVGREFL); /* compute ambient */ i = doambient(acol, r, amb.weight, uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral); scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */ if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */ return(i); /* store value */ VCOPY(amb.pos, r->rop); amb.ndir = encodedir(r->ron); amb.udir = encodedir(uvw[0]); amb.lvl = al; copycolor(amb.val, acol); /* insert into tree */ avsave(&amb); /* and save to file */ if (rn != r->ron) { /* texture */ VCOPY(uvw[2], r->ron); extambient(acol, &amb, r->rop, rn, uvw); } return(1); } static int extambient( /* extrapolate value at pv, nv */ COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv, FVECT uvw[3] ) { const double min_d = 0.05; const double max_d = 20.; static FVECT my_uvw[3]; FVECT v1; int i; double d = 1.0; /* zeroeth order */ if (uvw == NULL) { /* need local coordinates? */ decodedir(my_uvw[2], ap->ndir); decodedir(my_uvw[0], ap->udir); VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]); uvw = my_uvw; } for (i = 3; i--; ) /* gradient due to translation */ d += (pv[i] - ap->pos[i]) * (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]); VCROSS(v1, uvw[2], nv); /* gradient due to rotation */ for (i = 3; i--; ) d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]); if (d < min_d) /* clamp min/max scaling */ d = min_d; else if (d > max_d) d = max_d; copycolor(cr, ap->val); scalecolor(cr, d); return(d > min_d); } static void avinsert( /* insert ambient value in our tree */ AMBVAL *av ) { AMBTREE *at; AMBVAL *ap; AMBVAL avh; FVECT ck0; double s; int branch; int i; if (av->rad[0] <= FTINY) error(CONSISTENCY, "zero ambient radius in avinsert"); at = &atrunk; VCOPY(ck0, thescene.cuorg); s = thescene.cusize; while (s*(OCTSCALE/2) > av->rad[1]*ambacc) { if (at->kid == NULL) if ((at->kid = newambtree()) == NULL) error(SYSTEM, "out of memory in avinsert"); s *= 0.5; branch = 0; for (i = 0; i < 3; i++) if (av->pos[i] > ck0[i] + s) { ck0[i] += s; branch |= 1 << i; } at = at->kid + branch; } avh.next = at->alist; /* order by increasing level */ for (ap = &avh; ap->next != NULL; ap = ap->next) if ( ap->next->lvl > av->lvl || (ap->next->lvl == av->lvl) & (ap->next->weight <= av->weight) ) break; av->next = ap->next; ap->next = (AMBVAL*)av; at->alist = avh.next; } static void initambfile( /* initialize ambient file */ int cre8 ) { extern char *progname, *octname; static char *mybuf = NULL; #ifdef F_SETLKW aflock(cre8 ? F_WRLCK : F_RDLCK); #endif SET_FILE_BINARY(ambfp); if (mybuf == NULL) mybuf = (char *)bmalloc(BUFSIZ+8); setbuf(ambfp, mybuf); if (cre8) { /* new file */ newheader("RADIANCE", ambfp); fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ", progname, colval(ambval,RED), colval(ambval,GRN), colval(ambval,BLU), ambvwt, ambounce, ambacc); fprintf(ambfp, "-ad %d -as %d -ar %d ", ambdiv, ambssamp, ambres); if (octname != NULL) fputs(octname, ambfp); fputc('\n', ambfp); fprintf(ambfp, "SOFTWARE= %s\n", VersionID); fputnow(ambfp); fputformat(AMBFMT, ambfp); fputc('\n', ambfp); putambmagic(ambfp); } else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp)) error(USER, "bad ambient file"); } static void avsave( /* insert and save an ambient value */ AMBVAL *av ) { avstore(av); if (ambfp == NULL) return; if (writambval(av, ambfp) < 0) goto writerr; if (++nunflshed >= AMBFLUSH) if (ambsync() == EOF) goto writerr; return; writerr: error(SYSTEM, "error writing to ambient file"); } static AMBVAL * avstore( /* allocate memory and save aval */ AMBVAL *aval ) { AMBVAL *av; double d; if ((av = newambval()) == NULL) error(SYSTEM, "out of memory in avstore"); *av = *aval; av->latick = ambclock; av->next = NULL; nambvals++; d = bright(av->val); if (d > FTINY) { /* add to log sum for averaging */ avsum += log(d); navsum++; } avinsert(av); /* insert in our cache tree */ return(av); } #define ATALLOCSZ 512 /* #/8 trees to allocate at once */ static AMBTREE *atfreelist = NULL; /* free ambient tree structures */ static AMBTREE * newambtree(void) /* allocate 8 ambient tree structs */ { AMBTREE *atp, *upperlim; if (atfreelist == NULL) { /* get more nodes */ atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE)); if (atfreelist == NULL) return(NULL); /* link new free list */ upperlim = atfreelist + 8*(ATALLOCSZ-1); for (atp = atfreelist; atp < upperlim; atp += 8) atp->kid = atp + 8; atp->kid = NULL; } atp = atfreelist; atfreelist = atp->kid; memset(atp, 0, 8*sizeof(AMBTREE)); return(atp); } static void freeambtree( /* free 8 ambient tree structs */ AMBTREE *atp ) { atp->kid = atfreelist; atfreelist = atp; } static void unloadatree( /* unload an ambient value tree */ AMBTREE *at, unloadtf_t *f ) { AMBVAL *av; int i; /* transfer values at this node */ for (av = at->alist; av != NULL; av = at->alist) { at->alist = av->next; av->next = NULL; (*f)(av); } if (at->kid == NULL) return; for (i = 0; i < 8; i++) /* transfer and free children */ unloadatree(at->kid+i, f); freeambtree(at->kid); at->kid = NULL; } static struct avl { AMBVAL *p; unsigned long t; } *avlist1; /* ambient value list with ticks */ static AMBVAL **avlist2; /* memory positions for sorting */ static int i_avlist; /* index for lists */ static int alatcmp(const void *av1, const void *av2); static void avfree(AMBVAL *av) { free(av); } static void av2list( AMBVAL *av ) { #ifdef DEBUG if (i_avlist >= nambvals) error(CONSISTENCY, "too many ambient values in av2list1"); #endif avlist1[i_avlist].p = avlist2[i_avlist] = (AMBVAL*)av; avlist1[i_avlist++].t = av->latick; } static int alatcmp( /* compare ambient values for MRA */ const void *av1, const void *av2 ) { long lc = ((struct avl *)av2)->t - ((struct avl *)av1)->t; return(lc<0 ? -1 : lc>0 ? 1 : 0); } /* GW NOTE 2002/10/3: * I used to compare AMBVAL pointers, but found that this was the * cause of a serious consistency error with gcc, since the optimizer * uses some dangerous trick in pointer subtraction that * assumes pointers differ by exact struct size increments. */ static int aposcmp( /* compare ambient value positions */ const void *avp1, const void *avp2 ) { long diff = *(char * const *)avp1 - *(char * const *)avp2; if (diff < 0) return(-1); return(diff > 0); } static int avlmemi( /* find list position from address */ AMBVAL *avaddr ) { AMBVAL **avlpp; avlpp = (AMBVAL **)bsearch(&avaddr, avlist2, nambvals, sizeof(AMBVAL *), aposcmp); if (avlpp == NULL) error(CONSISTENCY, "address not found in avlmemi"); return(avlpp - avlist2); } static void sortambvals( /* resort ambient values */ int always ) { AMBTREE oldatrunk; AMBVAL tav, *tap, *pnext; int i, j; /* see if it's time yet */ if (!always && (ambclock++ < lastsort+sortintvl || nambvals < SORT_THRESH)) return; /* * The idea here is to minimize memory thrashing * in VM systems by improving reference locality. * We do this by periodically sorting our stored ambient * values in memory in order of most recently to least * recently accessed. This ordering was chosen so that new * ambient values (which tend to be less important) go into * higher memory with the infrequently accessed values. * Since we expect our values to need sorting less * frequently as the process continues, we double our * waiting interval after each call. * This routine is also called by setambacc() with * the "always" parameter set to 1 so that the ambient * tree will be rebuilt with the new accuracy parameter. */ if (tracktime) { /* allocate pointer arrays to sort */ avlist2 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *)); avlist1 = (struct avl *)malloc(nambvals*sizeof(struct avl)); } else { avlist2 = NULL; avlist1 = NULL; } if (avlist1 == NULL) { /* no time tracking -- rebuild tree? */ if (avlist2 != NULL) free(avlist2); if (always) { /* rebuild without sorting */ oldatrunk = atrunk; atrunk.alist = NULL; atrunk.kid = NULL; unloadatree(&oldatrunk, avinsert); } } else { /* sort memory by last access time */ /* * Sorting memory is tricky because it isn't contiguous. * We have to sort an array of pointers by MRA and also * by memory position. We then copy values in "loops" * to minimize memory hits. Nevertheless, we will visit * everyone at least twice, and this is an expensive process * when we're thrashing, which is when we need to do it. */ #ifdef DEBUG sprintf(errmsg, "sorting %u ambient values at ambclock=%lu...", nambvals, ambclock); eputs(errmsg); #endif i_avlist = 0; unloadatree(&atrunk, av2list); /* empty current tree */ #ifdef DEBUG if (i_avlist < nambvals) error(CONSISTENCY, "missing ambient values in sortambvals"); #endif qsort(avlist1, nambvals, sizeof(struct avl), alatcmp); qsort(avlist2, nambvals, sizeof(AMBVAL *), aposcmp); for (i = 0; i < nambvals; i++) { if (avlist1[i].p == NULL) continue; tap = avlist2[i]; tav = *tap; for (j = i; (pnext = avlist1[j].p) != tap; j = avlmemi(pnext)) { *(avlist2[j]) = *pnext; avinsert(avlist2[j]); avlist1[j].p = NULL; } *(avlist2[j]) = tav; avinsert(avlist2[j]); avlist1[j].p = NULL; } free(avlist1); free(avlist2); /* compute new sort interval */ sortintvl = ambclock - lastsort; if (sortintvl >= MAX_SORT_INTVL/2) sortintvl = MAX_SORT_INTVL; else sortintvl <<= 1; /* wait twice as long next */ #ifdef DEBUG eputs("done\n"); #endif } if (ambclock >= MAXACLOCK) ambclock = MAXACLOCK/2; lastsort = ambclock; } #ifdef F_SETLKW static void aflock( /* lock/unlock ambient file */ int typ ) { static struct flock fls; /* static so initialized to zeroes */ if (typ == fls.l_type) /* already called? */ return; fls.l_type = typ; do if (fcntl(fileno(ambfp), F_SETLKW, &fls) != -1) return; while (errno == EINTR); error(SYSTEM, "cannot (un)lock ambient file"); } int ambsync(void) /* synchronize ambient file */ { long flen; AMBVAL avs; int n; if (ambfp == NULL) /* no ambient file? */ return(0); /* gain appropriate access */ aflock(nunflshed ? F_WRLCK : F_RDLCK); /* see if file has grown */ if ((flen = lseek(fileno(ambfp), (off_t)0, SEEK_END)) < 0) goto seekerr; if ((n = flen - lastpos) > 0) { /* file has grown */ if (ambinp == NULL) { /* get new file pointer */ ambinp = fopen(ambfile, "rb"); if (ambinp == NULL) error(SYSTEM, "fopen failed in ambsync"); } if (fseek(ambinp, lastpos, SEEK_SET) < 0) goto seekerr; while (n >= AMBVALSIZ) { /* load contributed values */ if (!readambval(&avs, ambinp)) { sprintf(errmsg, "ambient file \"%s\" corrupted near character %ld", ambfile, flen - n); error(WARNING, errmsg); break; } avstore(&avs); n -= AMBVALSIZ; } lastpos = flen - n; /* check alignment */ if (n && lseek(fileno(ambfp), (off_t)lastpos, SEEK_SET) < 0) goto seekerr; } n = fflush(ambfp); /* calls write() at last */ lastpos += (long)nunflshed*AMBVALSIZ; aflock(F_UNLCK); /* release file */ nunflshed = 0; return(n); seekerr: error(SYSTEM, "seek failed in ambsync"); return(EOF); /* pro forma return */ } #else /* ! F_SETLKW */ int ambsync(void) /* flush ambient file */ { if (ambfp == NULL) return(0); nunflshed = 0; return(fflush(ambfp)); } #endif /* ! F_SETLKW */ diff --git a/m_bsdf.c b/m_bsdf.c index 943f639..5f1e56f 100644 --- a/m_bsdf.c +++ b/m_bsdf.c @@ -1,825 +1,803 @@ #ifndef lint -static const char RCSid[] = "$Id: m_bsdf.c,v 2.63 2021/03/27 20:08:58 greg Exp $"; +static const char RCSid[] = "$Id: m_bsdf.c,v 2.67 2021/12/07 23:49:50 greg Exp $"; #endif /* * Shading for materials with BSDFs taken from XML data files */ #include "copyright.h" #include "ray.h" #include "otypes.h" #include "ambient.h" #include "source.h" #include "func.h" #include "bsdf.h" #include "random.h" #include "pmapmat.h" /* - * Arguments to this material include optional diffuse colors. + * Arguments to this material include optional diffuse colors. * String arguments include the BSDF and function files. * For the MAT_BSDF type, a non-zero thickness causes the useful behavior * of translating transmitted rays this distance beneath the surface * (opposite the surface normal) to bypass any intervening geometry. * Translation only affects scattered, non-source-directed samples. * A non-zero thickness has the further side-effect that an unscattered * (view) ray will pass right through our material, making the BSDF * surface invisible and showing the proxied geometry instead. Thickness * has the further effect of turning off reflection on the reverse side so * rays heading in the opposite direction pass unimpeded through the BSDF * surface. A paired surface may be placed on the opposide side of * the detail geometry, less than this thickness away, if a two-way * proxy is desired. Note that the sign of the thickness is important. * A positive thickness hides geometry behind the BSDF surface and uses * front reflectance and transmission properties. A negative thickness * hides geometry in front of the surface when rays hit from behind, * and applies only the transmission and backside reflectance properties. * Reflection is ignored on the hidden side, as those rays pass through. * For the MAT_ABSDF type, we check for a strong "through" component. * Such a component will cause direct rays to pass through unscattered. * A separate test prevents over-counting by dropping samples that are * too close to this "through" direction. BSDFs with such a through direction * will also have a view component, meaning they are somewhat see-through. * A MAT_BSDF type with zero thickness behaves the same as a MAT_ABSDF * type with no strong through component. * The "up" vector for the BSDF is given by three variables, defined * (along with the thickness) by the named function file, or '.' if none. * Together with the surface normal, this defines the local coordinate * system for the BSDF. * We do not reorient the surface, so if the BSDF has no back-side * reflectance and none is given in the real arguments, a BSDF surface * with zero thickness will appear black when viewed from behind * unless backface visibility is on, when it becomes invisible. * The diffuse arguments are added to components in the BSDF file, * not multiplied. However, patterns affect this material as a multiplier * on everything except non-diffuse reflection. * * Arguments for MAT_ABSDF are: * 5+ BSDFfile ux uy uz funcfile transform * 0 * 0|3|6|9 rdf gdf bdf * rdb gdb bdb * rdt gdt bdt * * Arguments for MAT_BSDF are: * 6+ thick BSDFfile ux uy uz funcfile transform * 0 * 0|3|6|9 rdf gdf bdf * rdb gdb bdb * rdt gdt bdt */ /* * Note that our reverse ray-tracing process means that the positions * of incoming and outgoing vectors may be reversed in our calls * to the BSDF library. This is usually fine, since the bidirectional nature * of the BSDF (that's what the 'B' stands for) means it all works out. */ typedef struct { OBJREC *mp; /* material pointer */ RAY *pr; /* intersected ray */ FVECT pnorm; /* perturbed surface normal */ FVECT vray; /* local outgoing (return) vector */ double sr_vpsa[2]; /* sqrt of BSDF projected solid angle extrema */ RREAL toloc[3][3]; /* world to local BSDF coords */ RREAL fromloc[3][3]; /* local BSDF coords to world */ double thick; /* surface thickness */ COLOR cthru; /* "through" component for MAT_ABSDF */ COLOR cthru_surr; /* surround for "through" component */ SDData *sd; /* loaded BSDF data */ COLOR rdiff; /* diffuse reflection */ COLOR runsamp; /* BSDF hemispherical reflection */ COLOR tdiff; /* diffuse transmission */ COLOR tunsamp; /* BSDF hemispherical transmission */ } BSDFDAT; /* BSDF material data */ #define cvt_sdcolor(cv, svp) ccy2rgb(&(svp)->spec, (svp)->cieY, cv) typedef struct { double vy; /* brightness (for sorting) */ FVECT tdir; /* through sample direction (normalized) */ COLOR vcol; /* BTDF color */ } PEAKSAMP; /* BTDF peak sample */ /* Comparison function to put near-peak values in descending order */ static int cmp_psamp(const void *p1, const void *p2) { double diff = (*(const PEAKSAMP *)p1).vy - (*(const PEAKSAMP *)p2).vy; if (diff > 0) return(-1); if (diff < 0) return(1); return(0); } /* Compute "through" component color for MAT_ABSDF */ static void compute_through(BSDFDAT *ndp) { #define NDIR2CHECK 29 static const float dir2check[NDIR2CHECK][2] = { {0, 0}, {-0.6, 0}, {0, 0.6}, {0, -0.6}, {0.6, 0}, {-0.6, 0.6}, {-0.6, -0.6}, {0.6, 0.6}, {0.6, -0.6}, {-1.2, 0}, {0, 1.2}, {0, -1.2}, {1.2, 0}, {-1.2, 1.2}, {-1.2, -1.2}, {1.2, 1.2}, {1.2, -1.2}, {-1.8, 0}, {0, 1.8}, {0, -1.8}, {1.8, 0}, {-1.8, 1.8}, {-1.8, -1.8}, {1.8, 1.8}, {1.8, -1.8}, {-2.4, 0}, {0, 2.4}, {0, -2.4}, {2.4, 0}, }; - const double peak_over = 1.5; PEAKSAMP psamp[NDIR2CHECK]; SDSpectralDF *dfp; FVECT pdir; double tomega, srchrad; double tomsum, tomsurr; COLOR vpeak, vsurr; double vypeak; int i, ns; SDError ec; if (ndp->pr->rod > 0) dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb; else dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf; if (dfp == NULL) return; /* no specular transmission */ if (bright(ndp->pr->pcol) <= FTINY) return; /* pattern is black, here */ srchrad = sqrt(dfp->minProjSA); /* else evaluate peak */ for (i = 0; i < NDIR2CHECK; i++) { SDValue sv; psamp[i].tdir[0] = -ndp->vray[0] + dir2check[i][0]*srchrad; psamp[i].tdir[1] = -ndp->vray[1] + dir2check[i][1]*srchrad; psamp[i].tdir[2] = -ndp->vray[2]; normalize(psamp[i].tdir); - ec = SDevalBSDF(&sv, psamp[i].tdir, ndp->vray, ndp->sd); + ec = SDevalBSDF(&sv, ndp->vray, psamp[i].tdir, ndp->sd); if (ec) goto baderror; cvt_sdcolor(psamp[i].vcol, &sv); psamp[i].vy = sv.cieY; } qsort(psamp, NDIR2CHECK, sizeof(PEAKSAMP), cmp_psamp); if (psamp[0].vy <= FTINY) - return; /* zero area */ + return; /* zero BTDF here */ setcolor(vpeak, 0, 0, 0); setcolor(vsurr, 0, 0, 0); vypeak = tomsum = tomsurr = 0; /* combine top unique values */ ns = 0; for (i = 0; i < NDIR2CHECK; i++) { if (i && psamp[i].vy == psamp[i-1].vy) continue; /* assume duplicate sample */ - ec = SDsizeBSDF(&tomega, psamp[i].tdir, ndp->vray, + ec = SDsizeBSDF(&tomega, ndp->vray, psamp[i].tdir, SDqueryMin, ndp->sd); if (ec) goto baderror; - /* not really a peak? */ + + scalecolor(psamp[i].vcol, tomega); + /* not part of peak? */ if (tomega > 1.5*dfp->minProjSA || vypeak > 8.*psamp[i].vy*ns) { if (!i) return; /* abort */ - scalecolor(psamp[i].vcol, tomega); addcolor(vsurr, psamp[i].vcol); tomsurr += tomega; continue; } - scalecolor(psamp[i].vcol, tomega); addcolor(vpeak, psamp[i].vcol); tomsum += tomega; vypeak += psamp[i].vy; ++ns; } - if (vypeak*tomsurr < peak_over*bright(vsurr)*ns) - return; /* peak not peaky enough */ + if (tomsurr <= FTINY) /* no surround implies no peak */ + return; if ((vypeak/ns - (ndp->vray[2] > 0 ? ndp->sd->tLambFront.cieY - : ndp->sd->tLambBack.cieY)*(1./PI))*tomsum <= .001) - return; /* < 0.1% transmission */ + : ndp->sd->tLambBack.cieY)*(1./PI))*tomsum < .0005) + return; /* < 0.05% transmission */ copycolor(ndp->cthru, vpeak); /* already scaled by omega */ multcolor(ndp->cthru, ndp->pr->pcol); /* modify by pattern */ - if (tomsurr > FTINY) { /* surround contribution? */ - scalecolor(vsurr, 1./tomsurr); /* this one is avg. BTDF */ - copycolor(ndp->cthru_surr, vsurr); - multcolor(ndp->cthru_surr, ndp->pr->pcol); - } + scalecolor(vsurr, 1./tomsurr); /* surround is avg. BTDF */ + copycolor(ndp->cthru_surr, vsurr); + multcolor(ndp->cthru_surr, ndp->pr->pcol); return; baderror: objerror(ndp->mp, USER, transSDError(ec)); #undef NDIR2CHECK } /* Jitter ray sample according to projected solid angle and specjitter */ static void bsdf_jitter(FVECT vres, BSDFDAT *ndp, double sr_psa) { VCOPY(vres, ndp->vray); if (specjitter < 1.) sr_psa *= specjitter; if (sr_psa <= FTINY) return; vres[0] += sr_psa*(.5 - frandom()); vres[1] += sr_psa*(.5 - frandom()); normalize(vres); } /* Get BSDF specular for direct component, returning true if OK to proceed */ static int direct_specular_OK(COLOR cval, FVECT ldir, double omega, BSDFDAT *ndp) { - int nsamp; - double wtot = 0; - FVECT vsrc, vsmp, vjit; + int nsamp = 1; + int scnt = 0; + FVECT vsrc, vjit; double tomega, tomega2; double sf, tsr, sd[2]; COLOR csmp, cdiff; double diffY; SDValue sv; SDError ec; int i; /* in case we fail */ setcolor(cval, 0, 0, 0); /* transform source direction */ if (SDmapDir(vsrc, ndp->toloc, ldir) != SDEnone) return(0); /* will discount diffuse portion */ switch ((vsrc[2] > 0)<<1 | (ndp->vray[2] > 0)) { case 3: if (ndp->sd->rf == NULL) return(0); /* all diffuse */ sv = ndp->sd->rLambFront; break; case 0: if (ndp->sd->rb == NULL) return(0); /* all diffuse */ sv = ndp->sd->rLambBack; break; case 1: if ((ndp->sd->tf == NULL) & (ndp->sd->tb == NULL)) return(0); /* all diffuse */ sv = ndp->sd->tLambFront; break; case 2: if ((ndp->sd->tf == NULL) & (ndp->sd->tb == NULL)) return(0); /* all diffuse */ sv = ndp->sd->tLambBack; break; } if (sv.cieY > FTINY) { diffY = sv.cieY *= 1./PI; cvt_sdcolor(cdiff, &sv); } else { diffY = 0; setcolor(cdiff, 0, 0, 0); } - /* need projected solid angle */ - omega *= fabs(vsrc[2]); /* check indirect over-counting */ if ((vsrc[2] > 0) ^ (ndp->vray[2] > 0) && bright(ndp->cthru) > FTINY) { double dx = vsrc[0] + ndp->vray[0]; double dy = vsrc[1] + ndp->vray[1]; SDSpectralDF *dfp = (ndp->pr->rod > 0) ? ((ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb) : ((ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf) ; - if (dx*dx + dy*dy <= (2.5*4./PI)*(omega + dfp->minProjSA + - 2.*sqrt(omega*dfp->minProjSA))) { + tomega = omega*fabs(vsrc[2]); + if (dx*dx + dy*dy <= (2.5*4./PI)*(tomega + dfp->minProjSA + + 2.*sqrt(tomega*dfp->minProjSA))) { if (bright(ndp->cthru_surr) <= FTINY) return(0); copycolor(cval, ndp->cthru_surr); return(1); /* return non-zero surround BTDF */ } } ec = SDsizeBSDF(&tomega, ndp->vray, vsrc, SDqueryMin, ndp->sd); if (ec) goto baderror; - /* assign number of samples */ - sf = specjitter * ndp->pr->rweight; - if (tomega <= 0) - nsamp = 1; - else if (25.*tomega <= omega) - nsamp = 100.*sf + .5; - else - nsamp = 4.*sf*omega/tomega + .5; - nsamp += !nsamp; - sf = sqrt(omega); /* sample our source area */ - tsr = sqrt(tomega); + /* check if sampling BSDF */ + if ((tsr = sqrt(tomega)) > 0) { + nsamp = 4.*specjitter*ndp->pr->rweight + .5; + nsamp += !nsamp; + } + /* jitter to fuzz BSDF cells */ for (i = nsamp; i--; ) { - VCOPY(vsmp, vsrc); /* jitter query directions */ - if (nsamp > 1) { - multisamp(sd, 2, (i + frandom())/(double)nsamp); - vsmp[0] += (sd[0] - .5)*sf; - vsmp[1] += (sd[1] - .5)*sf; - normalize(vsmp); - } bsdf_jitter(vjit, ndp, tsr); /* compute BSDF */ - ec = SDevalBSDF(&sv, vjit, vsmp, ndp->sd); + ec = SDevalBSDF(&sv, vjit, vsrc, ndp->sd); if (ec) goto baderror; if (sv.cieY - diffY <= FTINY) continue; /* no specular part */ /* check for variable resolution */ - ec = SDsizeBSDF(&tomega2, vjit, vsmp, SDqueryMin, ndp->sd); + ec = SDsizeBSDF(&tomega2, vjit, vsrc, SDqueryMin, ndp->sd); if (ec) goto baderror; if (tomega2 < .12*tomega) continue; /* not safe to include */ cvt_sdcolor(csmp, &sv); -#if 0 - if (sf < 2.5*tsr) { /* weight by BSDF for small sources */ - scalecolor(csmp, sv.cieY); - wtot += sv.cieY; - } else -#endif - wtot += 1.; addcolor(cval, csmp); + ++scnt; } - if (wtot <= FTINY) /* no valid specular samples? */ + if (!scnt) /* no valid specular samples? */ return(0); - sf = 1./wtot; /* weighted average BSDF */ + sf = 1./scnt; /* weighted average BSDF */ scalecolor(cval, sf); /* subtract diffuse contribution */ for (i = 3*(diffY > FTINY); i--; ) if ((colval(cval,i) -= colval(cdiff,i)) < 0) colval(cval,i) = 0; return(1); baderror: objerror(ndp->mp, USER, transSDError(ec)); return(0); /* gratis return */ } /* Compute source contribution for BSDF (reflected & transmitted) */ static void dir_bsdf( COLOR cval, /* returned coefficient */ void *nnp, /* material data */ FVECT ldir, /* light source direction */ double omega /* light source size */ ) { BSDFDAT *np = (BSDFDAT *)nnp; double ldot; double dtmp; COLOR ctmp; setcolor(cval, 0, 0, 0); ldot = DOT(np->pnorm, ldir); if ((-FTINY <= ldot) & (ldot <= FTINY)) return; if (ldot > 0 && bright(np->rdiff) > FTINY) { /* * Compute diffuse reflected component */ copycolor(ctmp, np->rdiff); dtmp = ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ldot < 0 && bright(np->tdiff) > FTINY) { /* * Compute diffuse transmission */ copycolor(ctmp, np->tdiff); - dtmp = -ldot * omega * (1.0/PI); + dtmp = -ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ambRayInPmap(np->pr)) return; /* specular already in photon map */ /* * Compute specular scattering coefficient using BSDF */ if (!direct_specular_OK(ctmp, ldir, omega, np)) return; if (ldot < 0) { /* pattern for specular transmission */ multcolor(ctmp, np->pr->pcol); dtmp = -ldot * omega; } else dtmp = ldot * omega; scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } /* Compute source contribution for BSDF (reflected only) */ static void dir_brdf( COLOR cval, /* returned coefficient */ void *nnp, /* material data */ FVECT ldir, /* light source direction */ double omega /* light source size */ ) { BSDFDAT *np = (BSDFDAT *)nnp; double ldot; double dtmp; COLOR ctmp, ctmp1, ctmp2; setcolor(cval, 0, 0, 0); ldot = DOT(np->pnorm, ldir); if (ldot <= FTINY) return; if (bright(np->rdiff) > FTINY) { /* * Compute diffuse reflected component */ copycolor(ctmp, np->rdiff); dtmp = ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ambRayInPmap(np->pr)) return; /* specular already in photon map */ /* * Compute specular reflection coefficient using BSDF */ if (!direct_specular_OK(ctmp, ldir, omega, np)) return; dtmp = ldot * omega; scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } /* Compute source contribution for BSDF (transmitted only) */ static void dir_btdf( COLOR cval, /* returned coefficient */ void *nnp, /* material data */ FVECT ldir, /* light source direction */ double omega /* light source size */ ) { BSDFDAT *np = (BSDFDAT *)nnp; double ldot; double dtmp; COLOR ctmp; setcolor(cval, 0, 0, 0); ldot = DOT(np->pnorm, ldir); if (ldot >= -FTINY) return; if (bright(np->tdiff) > FTINY) { /* * Compute diffuse transmission */ copycolor(ctmp, np->tdiff); - dtmp = -ldot * omega * (1.0/PI); + dtmp = -ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ambRayInPmap(np->pr)) return; /* specular already in photon map */ /* * Compute specular scattering coefficient using BSDF */ if (!direct_specular_OK(ctmp, ldir, omega, np)) return; /* full pattern on transmission */ multcolor(ctmp, np->pr->pcol); dtmp = -ldot * omega; scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } /* Sample separate BSDF component */ static int sample_sdcomp(BSDFDAT *ndp, SDComponent *dcp, int xmit) { const int hasthru = (xmit && !(ndp->pr->crtype & (SPECULAR|AMBIENT)) && bright(ndp->cthru) > FTINY); int nstarget = 1; int nsent = 0; int n; SDError ec; SDValue bsv; double xrand; FVECT vsmp, vinc; RAY sr; /* multiple samples? */ if (specjitter > 1.5) { nstarget = specjitter*ndp->pr->rweight + .5; nstarget += !nstarget; } /* run through our samples */ for (n = 0; n < nstarget; n++) { if (nstarget == 1) { /* stratify random variable */ xrand = urand(ilhash(dimlist,ndims)+samplendx); if (specjitter < 1.) xrand = .5 + specjitter*(xrand-.5); } else { xrand = (n + frandom())/(double)nstarget; } SDerrorDetail[0] = '\0'; /* sample direction & coef. */ bsdf_jitter(vsmp, ndp, ndp->sr_vpsa[0]); VCOPY(vinc, vsmp); /* to compare after */ ec = SDsampComponent(&bsv, vsmp, xrand, dcp); if (ec) objerror(ndp->mp, USER, transSDError(ec)); if (bsv.cieY <= FTINY) /* zero component? */ break; if (hasthru) { /* check for view ray */ double dx = vinc[0] + vsmp[0]; double dy = vinc[1] + vsmp[1]; if (dx*dx + dy*dy <= ndp->sr_vpsa[0]*ndp->sr_vpsa[0]) continue; /* exclude view sample */ } /* map non-view sample->world */ if (SDmapDir(sr.rdir, ndp->fromloc, vsmp) != SDEnone) break; /* spawn a specular ray */ if (nstarget > 1) bsv.cieY /= (double)nstarget; cvt_sdcolor(sr.rcoef, &bsv); /* use sample color */ if (xmit) /* apply pattern on transmit */ multcolor(sr.rcoef, ndp->pr->pcol); if (rayorigin(&sr, SPECULAR, ndp->pr, sr.rcoef) < 0) { if (!n & (nstarget > 1)) { n = nstarget; /* avoid infinitue loop */ nstarget = nstarget*sr.rweight/minweight; if (n == nstarget) break; n = -1; /* moved target */ } continue; /* try again */ } if (xmit && ndp->thick != 0) /* need to offset origin? */ VSUM(sr.rorg, sr.rorg, ndp->pr->ron, -ndp->thick); rayvalue(&sr); /* send & evaluate sample */ multcolor(sr.rcol, sr.rcoef); addcolor(ndp->pr->rcol, sr.rcol); ++nsent; } return(nsent); } /* Sample non-diffuse components of BSDF */ static int sample_sdf(BSDFDAT *ndp, int sflags) { int hasthru = (sflags == SDsampSpT && !(ndp->pr->crtype & (SPECULAR|AMBIENT)) && bright(ndp->cthru) > FTINY); int n, ntotal = 0; double b = 0; SDSpectralDF *dfp; COLORV *unsc; if (sflags == SDsampSpT) { unsc = ndp->tunsamp; if (ndp->pr->rod > 0) dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb; else dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf; } else /* sflags == SDsampSpR */ { unsc = ndp->runsamp; if (ndp->pr->rod > 0) dfp = ndp->sd->rf; else dfp = ndp->sd->rb; } setcolor(unsc, 0, 0, 0); if (dfp == NULL) /* no specular component? */ return(0); if (hasthru) { /* separate view sample? */ RAY tr; if (rayorigin(&tr, TRANS, ndp->pr, ndp->cthru) == 0) { VCOPY(tr.rdir, ndp->pr->rdir); rayvalue(&tr); multcolor(tr.rcol, tr.rcoef); addcolor(ndp->pr->rcol, tr.rcol); ndp->pr->rxt = ndp->pr->rot + raydistance(&tr); ++ntotal; b = bright(ndp->cthru); } else hasthru = 0; } if (dfp->maxHemi - b <= FTINY) { /* have specular to sample? */ b = 0; } else { FVECT vjit; bsdf_jitter(vjit, ndp, ndp->sr_vpsa[1]); b = SDdirectHemi(vjit, sflags, ndp->sd) - b; if (b < 0) b = 0; } if (b <= specthresh+FTINY) { /* below sampling threshold? */ if (b > FTINY) { /* XXX no color from BSDF */ if (sflags == SDsampSpT) { copycolor(unsc, ndp->pr->pcol); scalecolor(unsc, b); } else /* no pattern on reflection */ setcolor(unsc, b, b, b); } return(ntotal); } dimlist[ndims] = (int)(size_t)ndp->mp; /* else sample specular */ ndims += 2; for (n = dfp->ncomp; n--; ) { /* loop over components */ dimlist[ndims-1] = n + 9438; ntotal += sample_sdcomp(ndp, &dfp->comp[n], sflags==SDsampSpT); } ndims -= 2; return(ntotal); } /* Color a ray that hit a BSDF material */ int m_bsdf(OBJREC *m, RAY *r) { int hasthick = (m->otype == MAT_BSDF); int hitfront; COLOR ctmp; SDError ec; FVECT upvec, vtmp; MFUNC *mf; BSDFDAT nd; /* check arguments */ if ((m->oargs.nsargs < hasthick+5) | (m->oargs.nfargs > 9) | (m->oargs.nfargs % 3)) objerror(m, USER, "bad # arguments"); /* record surface struck */ hitfront = (r->rod > 0); /* load cal file */ mf = hasthick ? getfunc(m, 5, 0x1d, 1) : getfunc(m, 4, 0xe, 1) ; setfunc(m, r); nd.thick = 0; /* set thickness */ if (hasthick) { nd.thick = evalue(mf->ep[0]); if ((-FTINY <= nd.thick) & (nd.thick <= FTINY)) nd.thick = 0; } /* check backface visibility */ if (!hitfront & !backvis) { raytrans(r); return(1); } /* check other rays to pass */ if (nd.thick != 0 && (r->crtype & SHADOW || !(r->crtype & (SPECULAR|AMBIENT)) || (nd.thick > 0) ^ hitfront)) { raytrans(r); /* hide our proxy */ return(1); } if (hasthick && r->crtype & SHADOW) /* early shadow check #1 */ return(1); nd.mp = m; nd.pr = r; /* get BSDF data */ nd.sd = loadBSDF(m->oargs.sarg[hasthick]); /* early shadow check #2 */ if (r->crtype & SHADOW && (nd.sd->tf == NULL) & (nd.sd->tb == NULL)) { SDfreeCache(nd.sd); return(1); } /* diffuse components */ if (hitfront) { cvt_sdcolor(nd.rdiff, &nd.sd->rLambFront); if (m->oargs.nfargs >= 3) { setcolor(ctmp, m->oargs.farg[0], m->oargs.farg[1], m->oargs.farg[2]); addcolor(nd.rdiff, ctmp); } cvt_sdcolor(nd.tdiff, &nd.sd->tLambFront); } else { cvt_sdcolor(nd.rdiff, &nd.sd->rLambBack); if (m->oargs.nfargs >= 6) { setcolor(ctmp, m->oargs.farg[3], m->oargs.farg[4], m->oargs.farg[5]); addcolor(nd.rdiff, ctmp); } cvt_sdcolor(nd.tdiff, &nd.sd->tLambBack); } if (m->oargs.nfargs >= 9) { /* add diffuse transmittance? */ setcolor(ctmp, m->oargs.farg[6], m->oargs.farg[7], m->oargs.farg[8]); addcolor(nd.tdiff, ctmp); } /* get modifiers */ raytexture(r, m->omod); /* modify diffuse values */ multcolor(nd.rdiff, r->pcol); multcolor(nd.tdiff, r->pcol); /* get up vector */ upvec[0] = evalue(mf->ep[hasthick+0]); upvec[1] = evalue(mf->ep[hasthick+1]); upvec[2] = evalue(mf->ep[hasthick+2]); /* return to world coords */ if (mf->fxp != &unitxf) { multv3(upvec, upvec, mf->fxp->xfm); nd.thick *= mf->fxp->sca; } if (r->rox != NULL) { multv3(upvec, upvec, r->rox->f.xfm); nd.thick *= r->rox->f.sca; } raynormal(nd.pnorm, r); /* compute local BSDF xform */ ec = SDcompXform(nd.toloc, nd.pnorm, upvec); if (!ec) { nd.vray[0] = -r->rdir[0]; nd.vray[1] = -r->rdir[1]; nd.vray[2] = -r->rdir[2]; ec = SDmapDir(nd.vray, nd.toloc, nd.vray); } if (ec) { objerror(m, WARNING, "Illegal orientation vector"); SDfreeCache(nd.sd); return(1); } setcolor(nd.cthru, 0, 0, 0); /* consider through component */ setcolor(nd.cthru_surr, 0, 0, 0); if (m->otype == MAT_ABSDF) { compute_through(&nd); if (r->crtype & SHADOW) { RAY tr; /* attempt to pass shadow ray */ SDfreeCache(nd.sd); if (rayorigin(&tr, TRANS, r, nd.cthru) < 0) return(1); /* no through component */ VCOPY(tr.rdir, r->rdir); rayvalue(&tr); /* transmit with scaling */ multcolor(tr.rcol, tr.rcoef); copycolor(r->rcol, tr.rcol); return(1); /* we're done */ } } ec = SDinvXform(nd.fromloc, nd.toloc); if (!ec) /* determine BSDF resolution */ ec = SDsizeBSDF(nd.sr_vpsa, nd.vray, NULL, SDqueryMin+SDqueryMax, nd.sd); if (ec) objerror(m, USER, transSDError(ec)); nd.sr_vpsa[0] = sqrt(nd.sr_vpsa[0]); nd.sr_vpsa[1] = sqrt(nd.sr_vpsa[1]); if (!hitfront) { /* perturb normal towards hit */ nd.pnorm[0] = -nd.pnorm[0]; nd.pnorm[1] = -nd.pnorm[1]; nd.pnorm[2] = -nd.pnorm[2]; } /* sample reflection */ sample_sdf(&nd, SDsampSpR); /* sample transmission */ sample_sdf(&nd, SDsampSpT); /* compute indirect diffuse */ copycolor(ctmp, nd.rdiff); addcolor(ctmp, nd.runsamp); if (bright(ctmp) > FTINY) { /* ambient from reflection */ if (!hitfront) flipsurface(r); multambient(ctmp, r, nd.pnorm); addcolor(r->rcol, ctmp); if (!hitfront) flipsurface(r); } copycolor(ctmp, nd.tdiff); addcolor(ctmp, nd.tunsamp); if (bright(ctmp) > FTINY) { /* ambient from other side */ FVECT bnorm; if (hitfront) flipsurface(r); bnorm[0] = -nd.pnorm[0]; bnorm[1] = -nd.pnorm[1]; bnorm[2] = -nd.pnorm[2]; if (nd.thick != 0) { /* proxy with offset? */ VCOPY(vtmp, r->rop); VSUM(r->rop, vtmp, r->ron, nd.thick); multambient(ctmp, r, bnorm); VCOPY(r->rop, vtmp); } else multambient(ctmp, r, bnorm); addcolor(r->rcol, ctmp); if (hitfront) flipsurface(r); } /* add direct component */ if ((bright(nd.tdiff) <= FTINY) & (nd.sd->tf == NULL) & (nd.sd->tb == NULL)) { direct(r, dir_brdf, &nd); /* reflection only */ } else if (nd.thick == 0) { direct(r, dir_bsdf, &nd); /* thin surface scattering */ } else { direct(r, dir_brdf, &nd); /* reflection first */ VCOPY(vtmp, r->rop); /* offset for transmitted */ VSUM(r->rop, vtmp, r->ron, -nd.thick); direct(r, dir_btdf, &nd); /* separate transmission */ VCOPY(r->rop, vtmp); } /* clean up */ SDfreeCache(nd.sd); return(1); } diff --git a/o_cone.c b/o_cone.c index 16d7aa1..4cf3345 100644 --- a/o_cone.c +++ b/o_cone.c @@ -1,143 +1,145 @@ #ifndef lint -static const char RCSid[] = "$Id: o_cone.c,v 2.10 2021/01/31 18:08:04 greg Exp $"; +static const char RCSid[] = "$Id: o_cone.c,v 2.11 2021/12/16 21:37:21 greg Exp $"; #endif /* * o_cone.c - routine to determine ray intersection with cones. */ #include "copyright.h" #include "ray.h" #include "otypes.h" #include "rtotypes.h" #include "cone.h" int o_cone( /* intersect ray with cone */ OBJREC *o, RAY *r ) { FVECT rox, rdx; double a, b, c; double root[2]; int nroots, rn; CONE *co; int i; /* get cone structure */ co = getcone(o, 1); if (co == NULL) objerror(o, INTERNAL, "unexpected illegal"); /* * To intersect a ray with a cone, we transform the * ray into the cone's normalized space. This greatly * simplifies the computation. * For a cone or cup, normalization results in the * equation: * * x*x + y*y - z*z == 0 * * For a cylinder or tube, the normalized equation is: * * x*x + y*y - r*r == 0 * * A normalized ring obeys the following set of equations: * * z == 0 && * x*x + y*y >= r0*r0 && * x*x + y*y <= r1*r1 */ /* transform ray */ multp3(rox, r->rorg, co->tm); multv3(rdx, r->rdir, co->tm); /* compute intersection */ if (o->otype == OBJ_CONE || o->otype == OBJ_CUP) { a = rdx[0]*rdx[0] + rdx[1]*rdx[1] - rdx[2]*rdx[2]; b = 2.0*(rdx[0]*rox[0] + rdx[1]*rox[1] - rdx[2]*rox[2]); c = rox[0]*rox[0] + rox[1]*rox[1] - rox[2]*rox[2]; } else if (o->otype == OBJ_CYLINDER || o->otype == OBJ_TUBE) { a = rdx[0]*rdx[0] + rdx[1]*rdx[1]; b = 2.0*(rdx[0]*rox[0] + rdx[1]*rox[1]); c = rox[0]*rox[0] + rox[1]*rox[1] - CO_R0(co)*CO_R0(co); } else { /* OBJ_RING */ if ((rdx[2] <= FTINY) & (rdx[2] >= -FTINY)) return(0); /* parallel */ root[0] = -rox[2]/rdx[2]; if (rayreject(o, r, root[0])) return(0); /* have better */ b = root[0]*rdx[0] + rox[0]; c = root[0]*rdx[1] + rox[1]; a = b*b + c*c; if (a < CO_R0(co)*CO_R0(co) || a > CO_R1(co)*CO_R1(co)) return(0); /* outside radii */ r->ro = o; r->rot = root[0]; VSUM(r->rop, r->rorg, r->rdir, r->rot); VCOPY(r->ron, co->ad); r->rod = -rdx[2]; + r->pert[0] = r->pert[1] = r->pert[2] = 0.0; + r->uv[0] = r->uv[1] = 0.0; r->rox = NULL; return(1); /* good */ } /* roots for cone, cup, cyl., tube */ nroots = quadratic(root, a, b, c); for (rn = 0; rn < nroots; rn++) { /* check real roots */ if (root[rn] <= FTINY) continue; /* too small */ if (root[rn] > r->rot + FTINY) break; /* too big */ /* check endpoints */ VSUM(rox, r->rorg, r->rdir, root[rn]); VSUB(rdx, rox, CO_P0(co)); b = DOT(rdx, co->ad); if (b < 0.0) continue; /* before p0 */ if (b > co->al) continue; /* after p1 */ if (rayreject(o, r, root[rn])) break; /* previous hit better */ r->ro = o; r->rot = root[rn]; VCOPY(r->rop, rox); /* get normal */ if (o->otype == OBJ_CYLINDER) a = CO_R0(co); else if (o->otype == OBJ_TUBE) a = -CO_R0(co); else { /* OBJ_CONE || OBJ_CUP */ c = CO_R1(co) - CO_R0(co); a = CO_R0(co) + b*c/co->al; if (o->otype == OBJ_CUP) { c = -c; a = -a; } } for (i = 0; i < 3; i++) r->ron[i] = (rdx[i] - b*co->ad[i])/a; if (o->otype == OBJ_CONE || o->otype == OBJ_CUP) for (i = 0; i < 3; i++) r->ron[i] = (co->al*r->ron[i] - c*co->ad[i]) / co->sl; a = DOT(r->ron, r->ron); if (a > 1.+FTINY || a < 1.-FTINY) { c = 1./(.5 + .5*a); /* avoid numerical error */ r->ron[0] *= c; r->ron[1] *= c; r->ron[2] *= c; } r->rod = -DOT(r->rdir, r->ron); r->pert[0] = r->pert[1] = r->pert[2] = 0.0; r->uv[0] = r->uv[1] = 0.0; r->rox = NULL; return(1); /* good */ } return(0); } diff --git a/pmapcontrib.c b/pmapcontrib.c index 58bf73f..7e15be3 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1535 +1,1533 @@ #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 $ */ #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; #ifdef PMAP_CONTRIB_DBG /* Sum/num of mRGBE encoding errors for wavelet coeffs */ static WAVELET_COEFF mrgbeDiffs = 0; static unsigned long mrgbeDiffCnt = 0; #endif -extern void SDdisk2square(double sq[2], double diskx, double disky); - static int ray2bin (const RAY *ray, unsigned scDim) /* Map ray dir (pointing away from origin) to its 1D bin in an (scDim x scDim) Shirley-Chiu square. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { 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; } #if 0 /* 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); + disk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return PMAP_CONTRIB_XY2LIN(scDim, scBin [0], scBin [1]); } #else /* 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 */ /* normalize(diskPlane); NEED TO NORMALISE? */ diskPlane [0] *= scRHS; diskx = DOT(ray -> rdir, diskPlane); disky = DOT(ray -> rdir, scUp); /* Apply Shirley-Chiu mapping of (diskx, disky) to square */ - SDdisk2square(scCoord, diskx, disky); + disk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return PMAP_CONTRIB_XY2LIN(scDim, scBin [0], scBin [1]); } #endif 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; unsigned numCoeffs; /* Warn if potential coefficient index overflow in mRGBE encoding. This requires getting the number of wavelet coefficients generated by the transform a priori. */ if (!(numCoeffs = padWaveletXform2(NULL, NULL, sqrt(binCnt), NULL))) { sprintf(errmsg, "can't determine number of wavelet coefficients " "for modifier %s", mod ); error(INTERNAL, errmsg); } if (numCoeffs * numCoeffs > PMAP_CONTRIB_MAXCOEFFS) { sprintf(errmsg, "mRGBE data field may overflow for modifier %s; reduce -bn " "and/or compression if contribution precomputation fails", mod ); error(WARNING, errmsg); } 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, sqrt(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 ) error(INTERNAL, "contribution source overflow"); } 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 -------------------- */ static void coeffSwap (PreComputedContribCoeff *c1, PreComputedContribCoeff *c2 ) { PreComputedContribCoeff tCoeff; tCoeff.coeff = c1 -> coeff; tCoeff.idx = c1 -> idx; c1 -> coeff = c2 -> coeff; c1 -> idx = c2 -> idx; c2 -> coeff = tCoeff.coeff; c2 -> idx = tCoeff.idx; } static int coeffPartition (PreComputedContribCoeff *coeffs, unsigned medianPos, unsigned left, unsigned right ) /* REVERSE partition coefficients by magnitude, such that coeffs [left..medianPos] >= coeffs [medianPos+1..right]. Returns median's position. */ { unsigned long l, r, m; WAVELET_COEFF lVal, rVal, mVal, tVal; #define COEFFVAL(c) (DOT((c) -> coeff, (c) -> coeff)) if (left < right) { /* Select pivot from median-of-three and move to photons [right] (a.k.a. Lomuto partitioning) */ l = left; r = right; m = l + ((r - l) >> 1); /* Avoids overflow vs. (l+r) >> 1 */ lVal = COEFFVAL(coeffs + l); rVal = COEFFVAL(coeffs + r); mVal = COEFFVAL(coeffs + m); if (mVal > lVal) { coeffSwap(coeffs + m, coeffs + l); tVal = mVal; mVal = lVal; lVal = tVal; } if (rVal > lVal) { coeffSwap(coeffs + r, coeffs + l); tVal = rVal; rVal = lVal; lVal = tVal; } if (mVal > rVal) { coeffSwap(coeffs + m, coeffs + r); tVal = mVal; mVal = rVal; rVal = tVal; } /* Pivot with key rVal is now in coeffs [right] */ /* l & r converge, swapping coefficients out of (reversed) order with respect to pivot. The convergence point l = r is the pivot's final position */ while (l < r) { while (l < r && COEFFVAL(coeffs + l) > rVal) l++; while (r > l && COEFFVAL(coeffs + r) <= rVal) r--; /* Coeffs out of order, must swap */ if (l < r) coeffSwap(coeffs + l, coeffs + r); }; /* Now l == r is pivot's final position; swap these coeffs */ coeffSwap(coeffs + l, coeffs + right); /* Recurse in partition containing median */ if (l > medianPos) return coeffPartition(coeffs, medianPos, left, l - 1); else if (l < medianPos) return coeffPartition(coeffs, medianPos, l + 1, right); else /* l == medianPos, partitioning done */ return l; } else return left; } static int coeffIdxCompare (const void *c1, const void *c2) /* Comparison function to sort thresholded coefficients by index */ { const PreComputedContribCoeff *tcoeff1 = c1, *tcoeff2 = c2; const unsigned v1 = tcoeff1 -> idx, v2 = tcoeff2 -> idx; if (v1 < v2) return -1; else if (v1 > v2) return 1; else return 0; } static int thresholdContribs (PreComputedContrib *preCompContrib) /* Threshold wavelet detail coefficients in preCompContrib -> waveletMatrix [approxDim..coeffDim-1] [approxDim..coeffDim-1] (where approxDim = WAVELET_PADD2_APPROXDIM) by keeping the (preCompContrib -> nCompressedCoeffs) largest of these and returning them in preCompContrib -> compressedCoeffs along with their original matrix indices. NOTE: The wavelet approximation coefficients preCompContrib -> waveletMatrix [0..approxDim-1] [0..approxDim-1] are excluded from thresholding to minimise compression artefacts. Returns 0 on success, else -1. */ { unsigned i, j, coeffDim, coeffIdx, nNzDetailCoeffs, nCompressedCoeffs, numThresh; WaveletMatrix2 waveletMatrix; PreComputedContribCoeff *threshCoeffs, *threshCoeffPtr; WAVELET_COEFF *coeffPtr; if (!preCompContrib || !(coeffDim = preCompContrib -> coeffDim) || !(threshCoeffs = preCompContrib -> compressedCoeffs) || !(waveletMatrix = preCompContrib -> waveletMatrix) ) /* Should be initialised by caller! */ return -1; /* Set up coefficient buffer for compression (thresholding), skipping * the approximation coefficients in the upper left of waveletMatrix, * which are not thresholded. Also skip zero coefficients (resulting * from padding), since these are already implicitly thresholded. The * original 2D matrix indices are linearised to 1D and saved with each * coefficient to restore the original sparse coefficient matrix. */ for (i = 0, threshCoeffPtr = threshCoeffs; i < coeffDim; i++) { /* Calc row pointer (=mult) in outer loop, inc in inner */ coeffIdx = PMAP_CONTRIB_XY2LIN(coeffDim, i, 0); for (j = 0; j < coeffDim; j++, coeffIdx++) { if ((i >= WAVELET_PADD2_APPROXDIM || j >= WAVELET_PADD2_APPROXDIM ) && !coeffIsZero(waveletMatrix [i] [j]) ) { /* Nonzero detail coefficient; set up pointer to coeff (sorts faster than 3 doubles) and save original (linearised) matrix index */ threshCoeffPtr -> idx = coeffIdx; threshCoeffPtr++ -> coeff = (WAVELET_COEFF*)&( waveletMatrix [i] [j] ); } } } /* Num of nonzero detail coeffs in buffer, number actually expected */ numThresh = threshCoeffPtr - threshCoeffs; nNzDetailCoeffs = WAVELET_PADD2_NUMDETAIL( preCompContrib -> nNonZeroCoeffs ); nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; /* If fewer nonzero detail coeff are in the threshold buffer than * anticipated, the loop below fills the remainder of the threshold * buffer with duplicates of a coefficient in the lower right of the * matrix, which is padding and guaranteed to be zero. This condition * can occur if the wavelet transform actually generated genuine zero * detail coefficients. Infact it's quite common if the wavelet * transformed contributions are quite sparse. */ i = j = coeffDim - 1; - coeffIdx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j), + coeffIdx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j); for (coeffPtr = (WAVELET_COEFF*)(waveletMatrix [i][j]); numThresh < nNzDetailCoeffs; numThresh++ ) { threshCoeffPtr -> idx = coeffIdx; threshCoeffPtr++ -> coeff = coeffPtr; } /* Partition coeffs in reverse order, such that all coeffs in threshCoeffs [0..nCompressedCoeffs-1] are larger than those in threshCoeffs [nCompressedCoeffs..nNzDetailCoeffs-1] */ coeffPartition(threshCoeffs, nCompressedCoeffs - 1, 0, nNzDetailCoeffs - 1 ); #ifdef PMAP_CONTRIB_DBG /* Check coefficient partitioning */ threshCoeffPtr = threshCoeffs + nCompressedCoeffs - 1; for (i = 0; i < nCompressedCoeffs; i++) if (COEFFVAL(threshCoeffs + i) < COEFFVAL(threshCoeffPtr)) error(CONSISTENCY, "invalid wavelet coefficient partitioning " "in thresholdContribs()" ); for (; i < nNzDetailCoeffs; i++) if (COEFFVAL(threshCoeffs + i) > COEFFVAL(threshCoeffPtr)) error(CONSISTENCY, "invalid wavelet coefficient partitioning " "in thresholdContribs()" ); #endif /* Sort partition containing nCompressedCoeffs coefficients by index * (ignoring the remaining coefficients since these are now dropped * due to tresholding) */ qsort(threshCoeffs, nCompressedCoeffs, sizeof(PreComputedContribCoeff), coeffIdxCompare ); return 0; } static int encodeContribs (PreComputedContrib *preCompContrib, float compressRatio ) /* Apply wavelet transform to input matrix preCompContrib -> waveletMatrix and compress according to compressRatio, storing thresholded and mRGBE-encoded coefficients in preCompContrib -> mrgbeCoeffs. Note that the approximation coefficients in the upper left of the matrix are not encoded, and returned in preCompContrib -> waveletMatrix [0..WAVELET_PADD2_APPROXDIM-1] [0..WAVELET_PADD2_APPROXDIM-1]. Returns 0 on success. */ { unsigned i, j, k, scDim, lastCoeffIdx; WaveletMatrix2 waveletMatrix, tWaveletMatrix; PreComputedContribCoeff *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; WAVELET_COEFF absCoeff; #ifdef PMAP_CONTRIB_DBG WaveletCoeff3 decCoeff; unsigned decIdx; #endif if (!preCompContrib || !preCompContrib -> mrgbeCoeffs || !preCompContrib -> compressedCoeffs || !(waveletMatrix = preCompContrib -> waveletMatrix) || !(tWaveletMatrix = preCompContrib -> tWaveletMatrix) || !(scDim = preCompContrib -> scDim) ) /* Should be initialised by the caller! */ return -1; #ifdef PMAP_CONTRIB_DBG for (i = 0; i < scDim; i++) { for (j = 0; j < scDim; j++) { for (k = 0; k < 3; k++) { #if 0 /* Set contributions to bins for debugging */ waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( PMAP_CONTRIB_XY2LIN(scDim, i, j) ); #elif 0 /* Replace contribs with "bump" function */ waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( (1. - fabs(i - scDim/2. + 0.5) / (scDim/2. - 0.5)) * (1. - fabs(j - scDim/2. + 0.5) / (scDim/2. - 0.5)) ); #elif 0 /* Inject reference contribs from rcontrib classic */ #include "rc-ref.c" if (preCompContrib -> nBins != PMAP_CONTRIB_REFSIZE) { sprintf(errmsg, "reference contribs require %d bins", PMAP_CONTRIB_REFSIZE ); error(USER, errmsg); } waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( refContrib [PMAP_CONTRIB_XY2LIN(scDim, i, j)] [k] ); #endif } #if 0 /* Dump contribs prior to encoding for debugging */ printf("% 7.3g\t", colorAvg(waveletMatrix [i] [j])); } putchar('\n'); } putchar('\n'); #else } } #endif #endif /* Do 2D wavelet transform on preCompContrib -> waveletMatrix */ if (padWaveletXform2(waveletMatrix, tWaveletMatrix, scDim, NULL) != preCompContrib -> coeffDim ) error(INTERNAL, "failed wavelet transform in encodeContribs()"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(preCompContrib) < 0) error(INTERNAL, "failed wavelet compression in encodeContribs()"); threshCoeffs = preCompContrib -> compressedCoeffs; /* 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 (i = 0; i < preCompContrib -> nCompressedCoeffs; i++) { for (k = 0; k < 3; k++) { #ifdef PMAP_CONTRIB_DBG #if 0 /* Replace wavelet coeff with its linear index for debugging */ threshCoeffs [i].coeff [k] = threshCoeffs [i].idx = i + 1; #endif #endif absCoeff = fabs(threshCoeffs [i].coeff [k]); if (absCoeff < mrgbeRange -> min [k]) mrgbeRange -> min [k] = absCoeff; if (absCoeff > mrgbeRange -> max [k]) mrgbeRange -> max [k] = absCoeff; } } if (preCompContrib -> nCompressedCoeffs == 1) /* Maximum compression with just 1 (!) compressed detail coeff (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 */ if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) return -1; mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; /* Encode wavelet detail coefficients to mRGBE */ for (i = lastCoeffIdx = 0; i < preCompContrib -> nCompressedCoeffs; lastCoeffIdx = threshCoeffs [i++].idx ) { /* HACK: To reduce the likelihood of overflowing the mRGBE data * field with the coefficient index, it is expressed incrementally * w.r.t. the previously encoded coefficient's index, instead of * absolutely. This implies threshCoeffs must be sorted by * coefficient index to ensure increments are positive, and to * minimise their magnitude. * Note that an overflow cannot be predicted beforehand, e.g. by * mkpmap when parsing the number of bins, as this depends on the * sparseness of the wavelet coefficients (which in turn depends on * the frequency distribution of the binned contributions), and the * fraction of those that are dropped (i.e. the compression ratio). * */ mrgbeCoeffs [i] = mRGBEencode(threshCoeffs [i].coeff, mrgbeRange, threshCoeffs [i].idx - lastCoeffIdx ); if (!mrgbeCoeffs [i].all) error(INTERNAL, "failed mRGBE encoding in encodeContribs()"); #ifdef PMAP_CONTRIB_DBG /* mRGBE encoding sanity check */ decIdx = mRGBEdecode(mrgbeCoeffs [i], mrgbeRange, decCoeff); #if 0 if (decIdx != threshCoeffs [i].idx - lastCoeffIdx || sqrt(dist2(decCoeff, threshCoeffs [i].coeff)) > 0.1 * colorAvg(mrgbeRange -> max) ) { sprintf(errmsg, "failed sanity check in encodeContribs()\n" "Encoded: [%.3g %.3g %.3g %d]\nDecoded: [%.3g %.3g %.3g %d]", threshCoeffs [i].coeff [0], threshCoeffs [i].coeff [1], threshCoeffs [i].coeff [2], threshCoeffs [i].idx - lastCoeffIdx, decCoeff [0], decCoeff [1], decCoeff [2], decIdx ); error(CONSISTENCY, errmsg); } #else mrgbeDiffs += sqrt(dist2(decCoeff, threshCoeffs [i].coeff)); mrgbeDiffCnt++; #endif for (k = 0; k < 3; k++) { absCoeff = fabs(threshCoeffs [i].coeff [k]) > 1e-6 ? decCoeff [k] / threshCoeffs [i].coeff [k] : 0; #if 0 if (absCoeff > 1e-6 && fabs(absCoeff - 1) > 1) { sprintf(errmsg, "mRGBE encoding out of tolerance in encodeContribs()\n" "Encoded: [%.3g %.3g %.3g]\nDecoded: [%.3g %.3g %.3g]", threshCoeffs [i].coeff [0], threshCoeffs [i].coeff [1], threshCoeffs [i].coeff [2], decCoeff [0], decCoeff [1], decCoeff [2] ); error(CONSISTENCY, errmsg); } #endif if (absCoeff < 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); */ } } #if 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). !!! HACK: THIS VERSION NORMALISES PER BIN !!! */ { 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, bin; float maxBinDist2 [contrib -> nbins]; /* 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; 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; #if 0 r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); #endif #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); norm = 2 / (PI * PI); #else /* Normalise accumulated flux by search area PI * r^2, including RADIANCE-specific lambertian factor PI */ // norm = 1 / (PI * PI * r2); norm = 1 / (PI * PI); #endif memset(maxBinDist2, 0, contrib -> nbins * sizeof(float)); for (i = 0, sqn = pmap -> squeue.node; i < pmap -> squeue.tail; i++, sqn++ ) { photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); bin = photonSrcBin(pmap, photon); if (sqn -> dist2 > maxBinDist2 [bin]) maxBinDist2 [bin] = sqn -> dist2; } /* Skip the extra photon */ for (i = 0, sqn = pmap -> squeue.node; i < pmap -> squeue.tail; i++, sqn++ ) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); bin = photonSrcBin(pmap, photon); r2 = maxBinDist2 [bin]; if (r2 > FTINY) { scalecolor(flux, norm / r2); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); addcolor(contrib -> cbin [bin], flux); if (!isnormal(contrib -> cbin [bin][0])) error(CONSISTENCY, "Inf/NaN contribs"); } } return contrib; } #else 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, bin; /* 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; 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)); #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 distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } #endif 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 */ freeWaveletMatrix2(preCompContrib -> waveletMatrix, preCompContrib -> coeffDim ); freeWaveletMatrix2(preCompContrib -> tWaveletMatrix, preCompContrib -> coeffDim ); /* Free thresholded coefficients */ free(preCompContrib -> compressedCoeffs); /* Free mRGBE encoded coefficients */ free(preCompContrib -> mrgbeCoeffs); free(preCompContrib -> waveletFname); if (preCompContrib -> cache) { /* Free contribution cache */ OOC_DeleteCache(preCompContrib -> cache); free(preCompContrib -> cache); } } /* 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 -> compressedCoeffs = 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 -> scDim = preCompContrib -> nBins = preCompContrib -> coeffDim = preCompContrib -> nCoeffs = preCompContrib -> nNonZeroCoeffs = preCompContrib -> nCompressedCoeffs = preCompContrib -> contribSize = 0; preCompContrib -> cache = NULL; } 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 p, numPreComp; unsigned i, j, k, coeffIdx, scDim, nBins, coeffDim, nCoeffs, nCompressedCoeffs, nNzDetailCoeffs; const double pInc = finalGather > FTINY ? 1 / finalGather : 0; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; COLR rgbe32 [WAVELET_PADD2_NUMAPPROX + 2]; LUENT *preCompContribNode; PhotonMap *preCompContribPmap, nuPmap; PreComputedContrib *preCompContrib; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); WaveletMatrix2 waveletMatrix; WaveletCoeff3 *coeffPtr; #ifdef PMAP_CONTRIB_DBG RAY dbgRay; #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; if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (p = 0; p < numPreComp; p++) { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering. * NOTE: multiplication with pInc expanded to prevent rounding * error which would occasionally cause the last photon to be * selected for precomputation again! */ pIdx = firstPhoton(pmap) + (unsigned long)( p * pInc + pmapRandom(pmap -> randState) * pInc ); 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 (i = 0; i < 3; i++) ray.ron [i] = photon.norm [i] / 127.0; /* Get contributions at photon position; retry with another photon if no contribs found */ binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); #ifdef PMAP_CONTRIB_DBG #if 0 /* Get all contribs from same photon for debugging */ /* Position must differ otherwise too many identical photon keys * will cause ooC octree to overflow! */ VCOPY(dbgRay.rop, photon.pos); getPhoton(pmap, 0, &photon); memcpy(&dbgRay, &ray, sizeof(RAY)); dbgRay.rsrc = photonSrcIdx(pmap, &photon); for (i = 0; i < 3; i++) dbgRay.ron [i] = photon.norm [i] / 127.0; binnedContribs = getPhotonContrib(pmap, &dbgRay, dbgRay.rcol); #endif #endif if (binnedContribs) { #if 0 if (!(p & PMAP_CNTUPDATE)) printf("Precomputing contribution photon %lu / %lu\n", p, numPreComp ); #endif /* 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; /* Set Shirley-Chiu square & wavelet matrix dimensions (number of bins resp. coefficients). */ preCompContrib -> scDim = scDim = sqrt( preCompContrib -> nBins = nBins = binnedContribs -> nbins ); if (scDim * scDim != nBins) /* Mkpmap shoulda took care of dis */ error(INTERNAL, "contribution bin count not " "integer square in preComputeContrib()" ); if (nBins > 1) { /* BINNING ENABLED; set up wavelet xform & compression. Get dimensions of wavelet coefficient matrix after boundary padding, and number of nonzero coefficients. The number of compressed (detail) coefficients is fixed and based on the number of NONZERO coefficients (minus approximations, see NOTE below) since zero coeffs are considered thresholded by default. !!! NOTE: THE APPROXIMATION COEFFICIENTS IN THE UPPER !!! LEFT OF THE MATRIX ARE _NEVER_ THRESHOLDED TO !!! MINIMISE COMPRESSION ARTEFACTS! THEY ARE ESSENTIAL !!! FOR RECONSTRUCTING THE ORIGINAL CONTRIBUTIONS. */ preCompContrib -> coeffDim = coeffDim = padWaveletXform2( NULL, NULL, scDim, &preCompContrib -> nNonZeroCoeffs ); preCompContrib -> nCoeffs = nCoeffs = coeffDim * coeffDim; nNzDetailCoeffs = WAVELET_PADD2_NUMDETAIL( preCompContrib -> nNonZeroCoeffs ); nCompressedCoeffs = (1 - compressRatio) * nNzDetailCoeffs; if (nCompressedCoeffs < 1) { error(WARNING, "maximum contribution compression, clamping number " "of wavelet coefficients in preComputeContrib()" ); nCompressedCoeffs = 1; } if (nCompressedCoeffs >= preCompContrib -> nCoeffs) { error(WARNING, "minimum contribution compression, clamping number " "of wavelet coefficients in preComputeContrib()" ); nCompressedCoeffs = nNzDetailCoeffs; } preCompContrib -> nCompressedCoeffs = nCompressedCoeffs; /* Lazily allocate primary and transposed wavelet * coefficient matrix */ preCompContrib -> waveletMatrix = waveletMatrix = allocWaveletMatrix2(coeffDim); preCompContrib -> tWaveletMatrix = allocWaveletMatrix2(coeffDim); if (!waveletMatrix || !preCompContrib -> tWaveletMatrix) error(SYSTEM, "out of memory allocating wavelet " "coefficient matrix in preComputeContrib()" ); /* Lazily allocate thresholded detail coefficients */ preCompContrib -> compressedCoeffs = calloc( nNzDetailCoeffs, sizeof(PreComputedContribCoeff) ); /* Lazily alloc mRGBE-encoded compressed wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc( nCompressedCoeffs, sizeof(mRGBE) ); if (!preCompContrib -> compressedCoeffs || !preCompContrib -> mrgbeCoeffs ) error(SYSTEM, "out of memory allocating compressed " "contribution buffer in preComputeContrib()" ); /* Set size of compressed contributions, in bytes */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( nCompressedCoeffs ); } } else { /* LUT node already exists for this modifier; use its data */ preCompContribPmap = (PhotonMap*)preCompContribNode -> data; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); scDim = preCompContrib -> scDim; nBins = preCompContrib -> nBins; nCoeffs = preCompContrib -> nCoeffs; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; waveletMatrix = preCompContrib -> waveletMatrix; } if (nBins > 1) { /* Binning enabled; copy binned contribs to wavelet xform * input matrix before compressing and encoding */ for (i = 0; i < scDim; i++) { /* Calc row pointer (=mult) in outer loop, inc in inner */ coeffPtr = binnedContribs -> cbin + PMAP_CONTRIB_XY2LIN(scDim, i, 0); #ifdef PMAP_CONTRIB_LOG /* Logarithmic contribs; copy & apply log per element */ for (j = 0; j < scDim; j++, coeffPtr++) for (k = 0; k < 3; k++) waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL((*coeffPtr) [k]); #else /* Linear contribs; copy matrix rows */ memcpy(waveletMatrix [i], *coeffPtr, scDim * WAVELET_COEFFSIZE ); #endif } /* Compress and encode contribs */ if (encodeContribs(preCompContrib, compressRatio)) error(INTERNAL, "failed contribution " "compression/encoding in preComputeContrib()" ); /* Encode wavelet approx coeffs in the upper left of the * wavelet matrix to 32-bit RGBE. These are not thresholded * to minimise compression artefacts. */ for (i = 0; i < WAVELET_PADD2_APPROXDIM; i++) { /* Calc row index (=mult) in outer loop, inc in inner */ coeffIdx = PMAP_CONTRIB_XY2LIN( WAVELET_PADD2_APPROXDIM, i, 0 ); for (j = 0; j < WAVELET_PADD2_APPROXDIM; j++, coeffIdx++) { setcolr(rgbe32 [coeffIdx], fabs(waveletMatrix [i] [j] [0]), fabs(waveletMatrix [i] [j] [1]), fabs(waveletMatrix [i] [j] [2]) ); /* HACK: depending on the boundary extension mode, some * approx coeff can be NEGATIVE (!), which 32-bit RGBE * doesn't accommodate. To get around this, we * sacrifice a bit in each mantissa for the sign. */ for (k = 0; k < 3; k++) rgbe32 [coeffIdx] [k] = PMAP_CONTRIB_SET_RGBE32_SGN( rgbe32 [coeffIdx] [k], waveletMatrix [i] [j] [k] ); } } /* Encode wavelet detail coeff range to 32-bit RGBE */ setcolr(rgbe32 [WAVELET_PADD2_NUMAPPROX], preCompContrib -> mrgbeRange.max [0], preCompContrib -> mrgbeRange.max [1], preCompContrib -> mrgbeRange.max [2] ); setcolr(rgbe32 [WAVELET_PADD2_NUMAPPROX + 1], preCompContrib -> mrgbeRange.min [0], preCompContrib -> mrgbeRange.min [1], preCompContrib -> mrgbeRange.min [2] ); /* Dump 32-bit RGBE approx coeffs followed by mRGBE range to * contribution heap file. NOTE: mRGBE range minimum and * maximum are reversed in the file to facilitate handling * the special (but pointless) case of a single compressed * coeff; if so, only the mRGBE maximum is dumped, since the * minimum is implicitly zero and can be omitted to save * space */ if (putbinary(rgbe32, sizeof(COLR), WAVELET_PADD2_NUMAPPROX + 1 + (nCompressedCoeffs > 1), preCompContribPmap -> contribHeap ) != WAVELET_PADD2_NUMAPPROX + 1 + (nCompressedCoeffs > 1) ) error(SYSTEM, "can't write wavelet approximation coefficients to " "contribution heap in preComputeContrib()" ); /* Now dump mRGBE encoded, compressed detail coefficients */ if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), nCompressedCoeffs, preCompContribPmap -> contribHeap ) != nCompressedCoeffs ) error(SYSTEM, "can't write wavelet detail coefficients to " "contribution heap in preComputeContrib()" ); } #if 0 /* HACK: signal newPhoton() to set precomputed photon's contribution source from ray -> rsrc */ /* XXX: DO WE STILL NEED THIS CRAP AFTER PRECOMP, SINCE CHILD PMAPS ARE PER-MODIFIER ANYWAY? */ preCompContribPmap -> lastContribSrc.srcIdx = -2; #endif /* 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 /* Bail out if no photons could be precomputed */ if (!pmap -> numPhotons) error(USER, "no contribution photons precomputed; try increasing -am" ); #ifdef PMAP_CONTRIB_DBG printf("*DBG* preComputeContrib: Avg mRGBE encoding diff = %g\n", mrgbeDiffs / mrgbeDiffCnt ); #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, scdim, 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 Shirley-Chiu dimension scDim>1, " "number of samples nsamp\n", stderr ); return -1; } if (!(scdim = atoi(argv [1]))) { fputs("Invalid Shirley-Chiu dimension\n", stderr); return -1; } 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, scdim); 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/portiotest.c b/portiotest.c deleted file mode 100644 index 05ea09e..0000000 --- a/portiotest.c +++ /dev/null @@ -1,17 +0,0 @@ -#include -#include -#include "rtio.h" - -int main () -{ - u_int32_t i = 0x88442211; - FILE *f = fopen("/tmp/portiotest.dat", "wb"); - - if (f) { - putint(i, 8, f); - fclose(f); - return 0; - } - else return -1; -} - diff --git a/raytrace.c b/raytrace.c index 4e91f70..e93a9fe 100644 --- a/raytrace.c +++ b/raytrace.c @@ -1,766 +1,771 @@ #ifndef lint -static const char RCSid[] = "$Id: raytrace.c,v 2.84 2021/01/31 20:55:04 greg Exp $"; +static const char RCSid[] = "$Id: raytrace.c,v 2.85 2021/06/09 18:21:10 greg Exp $"; #endif /* * raytrace.c - routines for tracing and shading rays. * * External symbols declared in ray.h */ #include "copyright.h" #include "ray.h" #include "source.h" #include "otypes.h" #include "otspecial.h" #include "random.h" #include "pmap.h" #define MAXCSET ((MAXSET+1)*2-1) /* maximum check set size */ RNUMBER raynum = 0; /* next unique ray number */ RNUMBER nrays = 0; /* number of calls to localhit */ static RREAL Lambfa[5] = {PI, PI, PI, 0.0, 0.0}; OBJREC Lamb = { OVOID, MAT_PLASTIC, "Lambertian", {NULL, Lambfa, 0, 5}, NULL }; /* a Lambertian surface */ OBJREC Aftplane; /* aft clipping plane object */ #define RAYHIT (-1) /* return value for intercepted ray */ static int raymove(FVECT pos, OBJECT *cxs, int dirf, RAY *r, CUBE *cu); static int checkhit(RAY *r, CUBE *cu, OBJECT *cxs); static void checkset(OBJECT *os, OBJECT *cs); int rayorigin( /* start new ray from old one */ RAY *r, int rt, const RAY *ro, const COLOR rc ) { double rw, re; /* assign coefficient/weight */ if (rc == NULL) { rw = 1.0; setcolor(r->rcoef, 1., 1., 1.); } else { rw = intens(rc); if (rw > 1.0) rw = 1.0; /* avoid calculation growth */ if (rc != r->rcoef) copycolor(r->rcoef, rc); } if ((r->parent = ro) == NULL) { /* primary ray */ r->rlvl = 0; r->rweight = rw; r->crtype = r->rtype = rt; r->rsrc = -1; r->clipset = NULL; r->revf = raytrace; copycolor(r->cext, cextinction); copycolor(r->albedo, salbedo); r->gecc = seccg; r->slights = NULL; } else { /* spawned ray */ if (ro->rot >= FHUGE*.99) { memset(r, 0, sizeof(RAY)); return(-1); /* illegal continuation */ } r->rlvl = ro->rlvl; if (rt & RAYREFL) { r->rlvl++; r->rsrc = -1; r->clipset = ro->clipset; r->rmax = 0.0; } else { r->rsrc = ro->rsrc; r->clipset = ro->newcset; r->rmax = ro->rmax <= FTINY ? 0.0 : ro->rmax - ro->rot; } r->revf = ro->revf; copycolor(r->cext, ro->cext); copycolor(r->albedo, ro->albedo); r->gecc = ro->gecc; r->slights = ro->slights; r->crtype = ro->crtype | (r->rtype = rt); VCOPY(r->rorg, ro->rop); r->rweight = ro->rweight * rw; /* estimate extinction */ re = colval(ro->cext,RED) < colval(ro->cext,GRN) ? colval(ro->cext,RED) : colval(ro->cext,GRN); if (colval(ro->cext,BLU) < re) re = colval(ro->cext,BLU); re *= ro->rot; if (re > 0.1) { if (re > 92.) { r->rweight = 0.0; } else { r->rweight *= exp(-re); } } } rayclear(r); if (r->rweight <= 0.0) /* check for expiration */ return(-1); if (r->crtype & SHADOW) /* shadow commitment */ return(0); /* ambient in photon map? */ if (ro != NULL && ro->crtype & AMBIENT) { if (causticPhotonMapping) return(-1); if (photonMapping && rt != TRANS) return(-1); } if ((maxdepth <= 0) & (rc != NULL)) { /* Russian roulette */ if (minweight <= 0.0) error(USER, "zero ray weight in Russian roulette"); if ((maxdepth < 0) & (r->rlvl > -maxdepth)) return(-1); /* upper reflection limit */ if (r->rweight >= minweight) return(0); if (frandom() > r->rweight/minweight) return(-1); rw = minweight/r->rweight; /* promote survivor */ scalecolor(r->rcoef, rw); r->rweight = minweight; return(0); } return((r->rweight >= minweight) & (r->rlvl <= abs(maxdepth)) ? 0 : -1); } void rayclear( /* clear a ray for (re)evaluation */ RAY *r ) { r->rno = raynum++; r->newcset = r->clipset; r->hitf = rayhit; r->robj = OVOID; r->ro = NULL; r->rox = NULL; r->rxt = r->rmt = r->rot = FHUGE; VCOPY(r->rop, r->rorg); r->ron[0] = -r->rdir[0]; r->ron[1] = -r->rdir[1]; r->ron[2] = -r->rdir[2]; r->rod = 1.0; r->pert[0] = r->pert[1] = r->pert[2] = 0.0; r->rflips = 0; r->uv[0] = r->uv[1] = 0.0; setcolor(r->pcol, 1.0, 1.0, 1.0); setcolor(r->mcol, 0.0, 0.0, 0.0); setcolor(r->rcol, 0.0, 0.0, 0.0); } void raytrace( /* trace a ray and compute its value */ RAY *r ) { if (localhit(r, &thescene)) raycont(r); /* hit local surface, evaluate */ else if (r->ro == &Aftplane) { r->ro = NULL; /* hit aft clipping plane */ r->rot = FHUGE; } else if (sourcehit(r)) rayshade(r, r->ro->omod); /* distant source */ if (trace != NULL) (*trace)(r); /* trace execution */ rayparticipate(r); /* for participating medium */ } void raycont( /* check for clipped object and continue */ RAY *r ) { if ((r->clipset != NULL && inset(r->clipset, r->ro->omod)) || !rayshade(r, r->ro->omod)) raytrans(r); } void raytrans( /* transmit ray as is */ RAY *r ) { RAY tr; rayorigin(&tr, TRANS, r, NULL); /* always continue */ VCOPY(tr.rdir, r->rdir); rayvalue(&tr); copycolor(r->mcol, tr.mcol); copycolor(r->rcol, tr.rcol); r->rmt = r->rot + tr.rmt; r->rxt = r->rot + tr.rxt; } int raytirrad( /* irradiance hack */ OBJREC *m, RAY *r ) { if (ofun[m->otype].flags & (T_M|T_X) && m->otype != MAT_CLIP) { if (istransp(m->otype) || isBSDFproxy(m)) { raytrans(r); return(1); } if (!islight(m->otype)) return((*ofun[Lamb.otype].funp)(&Lamb, r)); } return(0); /* not a qualifying surface */ } int rayshade( /* shade ray r with material mod */ RAY *r, int mod ) { int tst_irrad = do_irrad && !(r->crtype & ~(PRIMARY|TRANS)); OBJREC *m; r->rxt = r->rot; /* preset effective ray length */ for ( ; mod != OVOID; mod = m->omod) { m = objptr(mod); /****** unnecessary test since modifier() is always called if (!ismodifier(m->otype)) { sprintf(errmsg, "illegal modifier \"%s\"", m->oname); error(USER, errmsg); } ******/ /* hack for irradiance calculation */ if (tst_irrad && raytirrad(m, r)) return(1); if ((*ofun[m->otype].funp)(m, r)) return(1); /* materials call raytexture() */ } return(0); /* no material! */ } void rayparticipate( /* compute ray medium participation */ RAY *r ) { COLOR ce, ca; double re, ge, be; if (intens(r->cext) <= 1./FHUGE) return; /* no medium */ re = r->rot*colval(r->cext,RED); ge = r->rot*colval(r->cext,GRN); be = r->rot*colval(r->cext,BLU); if (r->crtype & SHADOW) { /* no scattering for sources */ re *= 1. - colval(r->albedo,RED); ge *= 1. - colval(r->albedo,GRN); be *= 1. - colval(r->albedo,BLU); } setcolor(ce, re<=FTINY ? 1. : re>92. ? 0. : exp(-re), ge<=FTINY ? 1. : ge>92. ? 0. : exp(-ge), be<=FTINY ? 1. : be>92. ? 0. : exp(-be)); multcolor(r->rcol, ce); /* path extinction */ if (r->crtype & SHADOW || intens(r->albedo) <= FTINY) return; /* no scattering */ /* PMAP: indirect inscattering accounted for by volume photons? */ if (!volumePhotonMapping) { setcolor(ca, colval(r->albedo,RED)*colval(ambval,RED)*(1.-colval(ce,RED)), colval(r->albedo,GRN)*colval(ambval,GRN)*(1.-colval(ce,GRN)), colval(r->albedo,BLU)*colval(ambval,BLU)*(1.-colval(ce,BLU))); addcolor(r->rcol, ca); /* ambient in scattering */ } srcscatter(r); /* source in scattering */ } void raytexture( /* get material modifiers */ RAY *r, OBJECT mod ) { OBJREC *m; /* execute textures and patterns */ for ( ; mod != OVOID; mod = m->omod) { m = objptr(mod); /****** unnecessary test since modifier() is always called if (!ismodifier(m->otype)) { sprintf(errmsg, "illegal modifier \"%s\"", m->oname); error(USER, errmsg); } ******/ if ((*ofun[m->otype].funp)(m, r)) { sprintf(errmsg, "conflicting material \"%s\"", m->oname); objerror(r->ro, USER, errmsg); } } } int raymixture( /* mix modifiers */ RAY *r, OBJECT fore, OBJECT back, double coef ) { RAY fr, br; double mfore, mback; int foremat, backmat; int i; /* bound coefficient */ if (coef > 1.0) coef = 1.0; else if (coef < 0.0) coef = 0.0; /* compute foreground and background */ foremat = backmat = 0; /* foreground */ fr = *r; if (coef > FTINY) { fr.rweight *= coef; scalecolor(fr.rcoef, coef); foremat = rayshade(&fr, fore); } /* background */ br = *r; if (coef < 1.0-FTINY) { br.rweight *= 1.0-coef; scalecolor(br.rcoef, 1.0-coef); backmat = rayshade(&br, back); } /* check for transparency */ if (backmat ^ foremat) { if (backmat && coef > FTINY) raytrans(&fr); else if (foremat && coef < 1.0-FTINY) raytrans(&br); } /* mix perturbations */ for (i = 0; i < 3; i++) r->pert[i] = coef*fr.pert[i] + (1.0-coef)*br.pert[i]; /* mix pattern colors */ scalecolor(fr.pcol, coef); scalecolor(br.pcol, 1.0-coef); copycolor(r->pcol, fr.pcol); addcolor(r->pcol, br.pcol); /* return value tells if material */ if (!foremat & !backmat) return(0); /* mix returned ray values */ scalecolor(fr.rcol, coef); scalecolor(br.rcol, 1.0-coef); copycolor(r->rcol, fr.rcol); addcolor(r->rcol, br.rcol); scalecolor(fr.mcol, coef); scalecolor(br.mcol, 1.0-coef); copycolor(r->mcol, fr.mcol); addcolor(r->mcol, br.mcol); mfore = bright(fr.mcol); mback = bright(br.mcol); r->rmt = mfore > mback ? fr.rmt : br.rmt; r->rxt = bright(fr.rcol)-mfore > bright(br.rcol)-mback ? fr.rxt : br.rxt; return(1); } double raydist( /* compute (cumulative) ray distance */ const RAY *r, int flags ) { double sum = 0.0; while (r != NULL && r->crtype&flags) { sum += r->rot; r = r->parent; } return(sum); } void raycontrib( /* compute (cumulative) ray contribution */ RREAL rc[3], const RAY *r, int flags ) { static int warnedPM = 0; rc[0] = rc[1] = rc[2] = 1.; while (r != NULL && r->crtype&flags) { int i = 3; while (i--) rc[i] *= colval(r->rcoef,i); /* check for participating medium */ if (!warnedPM && (bright(r->cext) > FTINY) | (bright(r->albedo) > FTINY)) { error(WARNING, "ray contribution calculation does not support participating media"); warnedPM++; } r = r->parent; } } double raynormal( /* compute perturbed normal for ray */ FVECT norm, RAY *r ) { double newdot; int i; /* The perturbation is added to the surface normal to obtain * the new normal. If the new normal would affect the surface * orientation wrt. the ray, a correction is made. The method is * still fraught with problems since reflected rays and similar * directions calculated from the surface normal may spawn rays behind * the surface. The only solution is to curb textures at high * incidence (namely, keep DOT(rdir,pert) < Rdot). */ for (i = 0; i < 3; i++) norm[i] = r->ron[i] + r->pert[i]; if (normalize(norm) == 0.0) { objerror(r->ro, WARNING, "illegal normal perturbation"); VCOPY(norm, r->ron); return(r->rod); } newdot = -DOT(norm, r->rdir); if ((newdot > 0.0) != (r->rod > 0.0)) { /* fix orientation */ for (i = 0; i < 3; i++) norm[i] += 2.0*newdot*r->rdir[i]; newdot = -newdot; } return(newdot); } void newrayxf( /* get new tranformation matrix for ray */ RAY *r ) { static struct xfn { struct xfn *next; FULLXF xf; } xfseed = { &xfseed }, *xflast = &xfseed; struct xfn *xp; const RAY *rp; /* * Search for transform in circular list that * has no associated ray in the tree. */ xp = xflast; for (rp = r->parent; rp != NULL; rp = rp->parent) if (rp->rox == &xp->xf) { /* xp in use */ xp = xp->next; /* move to next */ if (xp == xflast) { /* need new one */ xp = (struct xfn *)bmalloc(sizeof(struct xfn)); if (xp == NULL) error(SYSTEM, "out of memory in newrayxf"); /* insert in list */ xp->next = xflast->next; xflast->next = xp; break; /* we're done */ } rp = r; /* start check over */ } /* got it */ r->rox = &xp->xf; xflast = xp; } void flipsurface( /* reverse surface orientation */ RAY *r ) { r->rod = -r->rod; r->ron[0] = -r->ron[0]; r->ron[1] = -r->ron[1]; r->ron[2] = -r->ron[2]; r->pert[0] = -r->pert[0]; r->pert[1] = -r->pert[1]; r->pert[2] = -r->pert[2]; r->rflips++; } int rayreject( /* check if candidate hit is worse than current */ OBJREC *o, RAY *r, double t ) { OBJREC *mnew, *mray; if ((t <= FTINY) | (t > r->rot + FTINY)) return(1); if (t < r->rot - FTINY) /* is new hit significantly closer? */ return(0); /* coincident point, so decide... */ if (o == r->ro) return(1); /* shouldn't happen */ if (r->ro == NULL) return(0); /* ditto */ - if ((mnew = findmaterial(o)) == NULL) - return(1); /* new has no material */ - if ((mray = findmaterial(r->ro)) == NULL) + mnew = findmaterial(o); + mray = findmaterial(r->ro); /* check material transparencies */ + if (mnew == NULL) { + if (mray != NULL) + return(1); /* new has no material */ + } else if (mray == NULL) { return(0); /* old has no material(!) */ - if (istransp(mnew->otype)) - return(1); /* new is transparent */ - if (istransp(mray->otype)) + } else if (istransp(mnew->otype)) { + if (!istransp(mray->otype)) + return(1); /* new is transparent */ + } else if (istransp(mray->otype)) { return(0); /* old is transparent */ + } /* weakest priority to later modifier definition */ return (r->ro->omod >= o->omod); } void rayhit( /* standard ray hit test */ OBJECT *oset, RAY *r ) { OBJREC *o; int i; for (i = oset[0]; i > 0; i--) { o = objptr(oset[i]); if ((*ofun[o->otype].funp)(o, r)) r->robj = oset[i]; } } int localhit( /* check for hit in the octree */ RAY *r, CUBE *scene ) { OBJECT cxset[MAXCSET+1]; /* set of checked objects */ FVECT curpos; /* current cube position */ int sflags; /* sign flags */ double t, dt; int i; nrays++; /* increment trace counter */ sflags = 0; for (i = 0; i < 3; i++) { curpos[i] = r->rorg[i]; if (r->rdir[i] > 1e-7) sflags |= 1 << i; else if (r->rdir[i] < -1e-7) sflags |= 0x10 << i; } if (!sflags) { error(WARNING, "zero ray direction in localhit"); return(0); } /* start off assuming nothing hit */ if (r->rmax > FTINY) { /* except aft plane if one */ r->ro = &Aftplane; r->rot = r->rmax; VSUM(r->rop, r->rorg, r->rdir, r->rot); } /* find global cube entrance point */ t = 0.0; if (!incube(scene, curpos)) { /* find distance to entry */ for (i = 0; i < 3; i++) { /* plane in our direction */ if (sflags & 1<cuorg[i]; else if (sflags & 0x10<cuorg[i] + scene->cusize; else continue; /* distance to the plane */ dt = (dt - r->rorg[i])/r->rdir[i]; if (dt > t) t = dt; /* farthest face is the one */ } t += FTINY; /* fudge to get inside cube */ if (t >= r->rot) /* clipped already */ return(0); /* advance position */ VSUM(curpos, curpos, r->rdir, t); if (!incube(scene, curpos)) /* non-intersecting ray */ return(0); } cxset[0] = 0; raymove(curpos, cxset, sflags, r, scene); return((r->ro != NULL) & (r->ro != &Aftplane)); } static int raymove( /* check for hit as we move */ FVECT pos, /* current position, modified herein */ OBJECT *cxs, /* checked objects, modified by checkhit */ int dirf, /* direction indicators to speed tests */ RAY *r, CUBE *cu ) { int ax; double dt, t; if (istree(cu->cutree)) { /* recurse on subcubes */ CUBE cukid; int br, sgn; cukid.cusize = cu->cusize * 0.5; /* find subcube */ VCOPY(cukid.cuorg, cu->cuorg); br = 0; if (pos[0] >= cukid.cuorg[0]+cukid.cusize) { cukid.cuorg[0] += cukid.cusize; br |= 1; } if (pos[1] >= cukid.cuorg[1]+cukid.cusize) { cukid.cuorg[1] += cukid.cusize; br |= 2; } if (pos[2] >= cukid.cuorg[2]+cukid.cusize) { cukid.cuorg[2] += cukid.cusize; br |= 4; } for ( ; ; ) { cukid.cutree = octkid(cu->cutree, br); if ((ax = raymove(pos,cxs,dirf,r,&cukid)) == RAYHIT) return(RAYHIT); sgn = 1 << ax; if (sgn & dirf) /* positive axis? */ if (sgn & br) return(ax); /* overflow */ else { cukid.cuorg[ax] += cukid.cusize; br |= sgn; } else if (sgn & br) { cukid.cuorg[ax] -= cukid.cusize; br &= ~sgn; } else return(ax); /* underflow */ } /*NOTREACHED*/ } if (isfull(cu->cutree)) { if (checkhit(r, cu, cxs)) return(RAYHIT); } else if (r->ro == &Aftplane && incube(cu, r->rop)) return(RAYHIT); /* advance to next cube */ if (dirf&0x11) { dt = dirf&1 ? cu->cuorg[0] + cu->cusize : cu->cuorg[0]; t = (dt - pos[0])/r->rdir[0]; ax = 0; } else t = FHUGE; if (dirf&0x22) { dt = dirf&2 ? cu->cuorg[1] + cu->cusize : cu->cuorg[1]; dt = (dt - pos[1])/r->rdir[1]; if (dt < t) { t = dt; ax = 1; } } if (dirf&0x44) { dt = dirf&4 ? cu->cuorg[2] + cu->cusize : cu->cuorg[2]; dt = (dt - pos[2])/r->rdir[2]; if (dt < t) { t = dt; ax = 2; } } VSUM(pos, pos, r->rdir, t); return(ax); } static int checkhit( /* check for hit in full cube */ RAY *r, CUBE *cu, OBJECT *cxs ) { OBJECT oset[MAXSET+1]; objset(oset, cu->cutree); checkset(oset, cxs); /* avoid double-checking */ (*r->hitf)(oset, r); /* test for hit in set */ if (r->robj == OVOID) return(0); /* no scores yet */ return(incube(cu, r->rop)); /* hit OK if in current cube */ } static void checkset( /* modify checked set and set to check */ OBJECT *os, /* os' = os - cs */ OBJECT *cs /* cs' = cs + os */ ) { OBJECT cset[MAXCSET+MAXSET+1]; int i, j; int k; /* copy os in place, cset <- cs */ cset[0] = 0; k = 0; for (i = j = 1; i <= os[0]; i++) { while (j <= cs[0] && cs[j] < os[i]) cset[++cset[0]] = cs[j++]; if (j > cs[0] || os[i] != cs[j]) { /* object to check */ os[++k] = os[i]; cset[++cset[0]] = os[i]; } } if (!(os[0] = k)) /* new "to check" set size */ return; /* special case */ while (j <= cs[0]) /* get the rest of cs */ cset[++cset[0]] = cs[j++]; if (cset[0] > MAXCSET) /* truncate "checked" set if nec. */ cset[0] = MAXCSET; /* setcopy(cs, cset); */ /* copy cset back to cs */ os = cset; for (i = os[0]; i-- >= 0; ) *cs++ = *os++; } diff --git a/rc-ref.c b/rc-ref.c deleted file mode 100644 index c3e55f3..0000000 --- a/rc-ref.c +++ /dev/null @@ -1,67 +0,0 @@ -#define PMAP_CONTRIB_REFSIZE 64 -WaveletCoeff3 refContrib [PMAP_CONTRIB_REFSIZE] = { - { 4.409082e+02, 8.362935e+02, 8.362935e+02 }, - { 4.866695e+01, 9.851163e+01, 9.851163e+01 }, - { 2.049961e+01, 4.342454e+01, 4.342454e+01 }, - { 2.824923e+00, 6.434741e+00, 6.434741e+00 }, - { 1.830407e+00, 4.123382e+00, 4.123382e+00 }, - { 1.675290e+01, 3.534576e+01, 3.534576e+01 }, - { 4.701529e+01, 9.502763e+01, 9.502763e+01 }, - { 4.347457e+02, 8.244728e+02, 8.244728e+02 }, - { 7.044382e+02, 1.310281e+03, 1.310281e+03 }, - { 1.879216e+01, 3.796805e+01, 3.796805e+01 }, - { 7.240023e+00, 1.539449e+01, 1.539449e+01 }, - { 1.047041e+00, 2.365272e+00, 2.365272e+00 }, - { 6.742821e-01, 1.526210e+00, 1.526210e+00 }, - { 5.864430e+00, 1.256767e+01, 1.256767e+01 }, - { 1.728685e+01, 3.517820e+01, 3.517820e+01 }, - { 7.044838e+02, 1.310510e+03, 1.310510e+03 }, - { 7.703406e+02, 1.407557e+03, 1.407557e+03 }, - { 6.668084e+01, 1.278029e+02, 1.278029e+02 }, - { 6.110347e+00, 1.289413e+01, 1.289413e+01 }, - { 7.409838e-01, 1.674399e+00, 1.674399e+00 }, - { 5.349032e-01, 1.213084e+00, 1.213084e+00 }, - { 5.561961e+00, 1.182811e+01, 1.182811e+01 }, - { 7.973283e+01, 1.525534e+02, 1.525534e+02 }, - { 7.770507e+02, 1.419499e+03, 1.419499e+03 }, - { 7.860980e+02, 1.420810e+03, 1.420810e+03 }, - { 1.610307e+02, 3.036264e+02, 3.036264e+02 }, - { 7.655327e+00, 1.582977e+01, 1.582977e+01 }, - { 8.409774e-01, 1.863129e+00, 1.863129e+00 }, - { 6.460242e-01, 1.438277e+00, 1.438277e+00 }, - { 7.392244e+00, 1.545713e+01, 1.545713e+01 }, - { 1.139023e+02, 2.152189e+02, 2.152189e+02 }, - { 7.508683e+02, 1.358289e+03, 1.358289e+03 }, - { 8.130104e+02, 1.471817e+03, 1.471817e+03 }, - { 1.144158e+02, 2.158896e+02, 2.158896e+02 }, - { 7.170834e+00, 1.481093e+01, 1.481093e+01 }, - { 8.531226e-01, 1.887482e+00, 1.887482e+00 }, - { 6.437121e-01, 1.434679e+00, 1.434679e+00 }, - { 7.548666e+00, 1.577723e+01, 1.577723e+01 }, - { 1.140006e+02, 2.155748e+02, 2.155748e+02 }, - { 8.350603e+02, 1.511337e+03, 1.511337e+03 }, - { 7.734335e+02, 1.412687e+03, 1.412687e+03 }, - { 6.635901e+01, 1.270131e+02, 1.270131e+02 }, - { 5.797712e+00, 1.217059e+01, 1.217059e+01 }, - { 6.764597e-01, 1.516001e+00, 1.516001e+00 }, - { 5.276446e-01, 1.195057e+00, 1.195057e+00 }, - { 5.786208e+00, 1.228171e+01, 1.228171e+01 }, - { 6.685527e+01, 1.284225e+02, 1.284225e+02 }, - { 7.774189e+02, 1.420067e+03, 1.420067e+03 }, - { 7.031177e+02, 1.307338e+03, 1.307338e+03 }, - { 1.831548e+01, 3.670184e+01, 3.670184e+01 }, - { 6.543461e+00, 1.376826e+01, 1.376826e+01 }, - { 8.976740e-01, 2.003200e+00, 2.003200e+00 }, - { 7.976804e-01, 1.790019e+00, 1.790019e+00 }, - { 6.911056e+00, 1.472236e+01, 1.472236e+01 }, - { 1.879253e+01, 3.813792e+01, 3.813792e+01 }, - { 7.137211e+02, 1.328069e+03, 1.328069e+03 }, - { 4.371586e+02, 8.275908e+02, 8.275908e+02 }, - { 4.985108e+01, 9.986548e+01, 9.986548e+01 }, - { 1.919958e+01, 4.001629e+01, 4.001629e+01 }, - { 2.439399e+00, 5.400287e+00, 5.400287e+00 }, - { 2.640548e+00, 5.892877e+00, 5.892877e+00 }, - { 1.992110e+01, 4.182727e+01, 4.182727e+01 }, - { 5.089465e+01, 1.027822e+02, 1.027822e+02 }, - { 4.466842e+02, 8.461234e+02, 8.461234e+02 } -}; diff --git a/rgbecontrib.c b/rgbecontrib.c deleted file mode 100644 index 7c3767e..0000000 --- a/rgbecontrib.c +++ /dev/null @@ -1,282 +0,0 @@ -/* - ========================================================================= - Routines to encode/decode RGBE representation of normalised binned - contributions. See rgbecontrib.h for encoding details. - - Compile with -D_RGBE_CONTRIB_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 "rgbecontrib.h" -#include -#include -#include - - - -#define min(a,b) ((a) < (b) ? (a) : (b)) -#define max(a,b) ((a) > (b) ? (a) : (b)) - - - -RGBEContrib RGBEContribEnc (double rgb [3], double norm [3], unsigned bin) -/* Normalise signed float RGB contributions against abs(norm), and encode - * to RGBE with associated bin. Issues warning if norm < abs(rgb). */ -{ - #define sgn(a) ((a) >= 0 ? 1.0 : -1.0) - - double rgbNorm [3], mantNorm; - int i, exp; - RGBEContrib rgbeCont = {0, 0, 0, 0}; - - /* Check and set bin */ - if (bin > RGBE_CONTRIB_BINMAX) { - fprintf(stderr, "WARNING - RGBEContrib bin overflow\n"); - return rgbeCont; - } - else - rgbeCont.bin = bin; - - /* Normalise RGB. */ - for (i = 0; i < 3; i++) - if (norm [i] == 0) { - fprintf( - stderr, "WARNING - Zero RGBEContrib normalisation\n" - ); - return rgbeCont; - } - else { - rgbNorm [i] = rgb [i] / norm [i]; - - if (fabs(fabs(rgbNorm [i]) - 1) < RGBE_CONTRIB_MIN) - /* Normalisation to 1. This must be avoided by offsetting by - minimum value, otherwise frexp() returns a positive exponent. - Normalisation to >1 is caught by overflow check below */ - rgbNorm [i] = sgn(rgbNorm [i]) * (1 - RGBE_CONTRIB_MIN); - } - - /* Mantissa normalisation factor (= RGB maximum) */ - mantNorm = max( - max(fabs(rgbNorm [0]), fabs(rgbNorm [1])), fabs(rgbNorm [2]) - ); - - /* Check RGB range */ - if (mantNorm < RGBE_CONTRIB_MIN) { - /* Below encoding range, clamp to zero */ - fprintf(stderr, "WARNING - RGBEContrib underflow, clamping to 0\n"); - return rgbeCont; - } - - if (mantNorm > RGBE_CONTRIB_MAX) { - /* Exceeds encoding range, clamp to maximum */ - fprintf( - stderr, "WARNING - RGBEContrib overflow, clamping to +/-%.1f\n", - RGBE_CONTRIB_MAX - ); - - rgbeCont.red = (1 + sgn(rgbNorm [0])) * RGBE_CONTRIB_MANTMAX + 0.5; - rgbeCont.grn = (1 + sgn(rgbNorm [1])) * RGBE_CONTRIB_MANTMAX + 0.5; - rgbeCont.blu = (1 + sgn(rgbNorm [2])) * RGBE_CONTRIB_MANTMAX + 0.5; - - return rgbeCont; - } - - /* Get normalised mantissa and exponent to base 2 */ - mantNorm = frexp(mantNorm, &exp) * RGBE_CONTRIB_MANTMAX / mantNorm; - - /* Set RGB integer mantissas, adding excess (zero offset) and rounding - up or down, depending on sign. Clamp to RGBE_CONTRIB_MANTOVF to - prevent rounding beyond overflow into negative! */ - rgbeCont.red = min( - (int)(rgbNorm [0] * mantNorm + RGBE_CONTRIB_MANTMAX + 0.5), - RGBE_CONTRIB_MANTOVF - ); - rgbeCont.grn = min( - (int)(rgbNorm [1] * mantNorm + RGBE_CONTRIB_MANTMAX + 0.5), - RGBE_CONTRIB_MANTOVF - ); - rgbeCont.blu = min( - (int)(rgbNorm [2] * mantNorm + RGBE_CONTRIB_MANTMAX + 0.5), - RGBE_CONTRIB_MANTOVF - ); - - /* Drop the exponent sign, since implicitly negative */ - rgbeCont.exp = abs(exp); - - return rgbeCont; - - #undef sgn -} - - - -unsigned RGBEContribDec ( - RGBEContrib rgbeCont, double norm [3], double rgb [3] -) -/* Decode RGBE contribution, returning signed floating point contributions - * in array rgb after denormalisation against norm, and the associated - * bin as return value. */ -{ - #define sgn(a) ((a) >= RGBE_CONTRIB_MANTMAX ? 0.1 : -0.1) - - /* Mantissa normalisation factor; exponent implicitly negative */ - const double mantNorm = ldexp(1.0, -rgbeCont.exp) / RGBE_CONTRIB_MANTMAX; - - /* Decode and set RGB components, subtracting excess (zero offset) and - * adding some jitter to break up quantisation artefacts. - * Clamp to -RGBE_CONTRIB_MAX (-1) in case jitter underflows. - * Denormalise to original range with user-specified factor. */ - rgb [0] = norm [0] * max( - mantNorm * ( - rgbeCont.red - RGBE_CONTRIB_MANTMAX + drand48() * sgn(rgbeCont.red) - ), - -RGBE_CONTRIB_MAX - ); - rgb [1] = norm [1] * max( - mantNorm * ( - rgbeCont.grn - RGBE_CONTRIB_MANTMAX + drand48() * sgn(rgbeCont.grn) - ), - -RGBE_CONTRIB_MAX - ); - rgb [2] = norm [2] * max( - mantNorm * ( - rgbeCont.blu - RGBE_CONTRIB_MANTMAX + drand48() * sgn(rgbeCont.blu) - ), - -RGBE_CONTRIB_MAX - ); - - /* Return bin */ - return rgbeCont.bin; - - #undef sgn -} - - - -#ifdef _RGBE_CONTRIB_TEST - /* Build standalone to run unit test */ - - #define RGBMAX 1 - - int main (int argc, char **argv) - { - double rgb [3] = {0, 0, 0}, rgbDec [3], rgbNorm [3] = {0, 0, 0}, - relDiff, rmse, rmseSum; - const short rndState0 [3] = {11, 23, 31}; - short encFlag = 0, rndState [3]; - unsigned numTests, i, j, k, bin, binDec; - RGBEContrib rgbeCont; - - /* Test range checks */ - puts("Check bin overflow"); - RGBEContribEnc(rgb, rgbNorm, RGBE_CONTRIB_BINMAX << 1); - puts("\nCheck zero normalisation"); - RGBEContribEnc(rgb, rgbNorm, 0); - puts("\nCheck underflow"); - rgbNorm [0] = rgbNorm [1] = rgbNorm [2] = 1; - rgb [0] = RGBE_CONTRIB_MIN / 2; - RGBEContribEnc(rgb, rgbNorm, 0); - puts("\nCheck overflow (positive)"); - rgb [0] = 2; - RGBEContribEnc(rgb, rgbNorm, 0); - puts("\nCheck overflow (negative)"); - rgb [0] = -2; - RGBEContribEnc(rgb, rgbNorm, 0); - putchar('\n'); - - if (argc > 1 && (numTests = atoi(argv [1])) > 0) { - rgbNorm [0] = rgbNorm [1] = rgbNorm [2] = 0; - /* Iterate twice with same random RGB tuples; gather normalisation - factors in 1st pass, then do actual encode/decode tests using - found normalisation 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++) { - bin = (unsigned)(RGBE_CONTRIB_BINMAX * erand48(rndState)); - #if 1 - for (i = 0; i < 3; i++) { - if (i > 0) - /* Correlate RGB */ - rgb [i] = rgb [0] * (0.5 + 0.5 * erand48(rndState)); - else - rgb [i] = RGBMAX * (2 * erand48(rndState) - 1); - } - #else - #if 1 - rgb [0] = rgb [1] = rgb [2] = RGBMAX * ( - 2 * erand48(rndState) - 1 - ); - #else - rgb [0] = 0.9772; - rgb [1] = 0.4369; - rgb [2] = 0.1869; - #endif - #endif - - if (encFlag) { - /* Encode */ - /* - rgb [0] = rgbNorm [0]; - rgb [1] = 0.0077; - rgb [2] = 0.0058; - */ - rgbeCont = RGBEContribEnc(rgb, rgbNorm, bin); - 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], bin, - rgbeCont.red, rgbeCont.grn, rgbeCont.blu, rgbeCont.exp - ); - - /* Decode */ - binDec = RGBEContribDec(rgbeCont, rgbNorm, 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], binDec, 100 * rmse - ); - } - else { - /* Update RGB normalisation factor */ - for (i = 0; i < 3; i++) - if (fabs(rgb [i]) > rgbNorm [i]) - rgbNorm [i] = fabs(rgb [i]); - } - } - - /* Dump results */ - if (!encFlag) - printf( - "RGB Normalisation = (%.3f, %.3f, %.3f)\n\n", - rgbNorm [0], rgbNorm [1], rgbNorm [2] - ); - else - printf("\nAvg RMSE = %.1f%%\n", 100 * rmseSum / numTests); - } while (encFlag = !encFlag); - } - - return 0; - } -#endif diff --git a/rgbecontrib.h b/rgbecontrib.h deleted file mode 100644 index 7be7693..0000000 --- a/rgbecontrib.h +++ /dev/null @@ -1,74 +0,0 @@ -/* - ========================================================================= - Reduced RGBE representation of normalised coefficients. Per default, this - consists of three signed 5-bit mantissas for each RGB colour channel and - a common 5-bit exponent (= 20 bits). This allocation can be modified - within a 32-bit envelope if more precision is required. - - The floating point RGB contributions are assumed to lie in the range - [-1, 1], with values in the interval [-RGBE_CONTRIB_MIN, - RGBE_CONTRIB_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 _RGBE_CONTRIB_H - #define _RGBE_CONTRIB_H - - #include - - - /* RGBE Mantissa / exponent / bin bits and their encoding ranges - (NOTE: The binary shifts can overflow if the number of bits is - increased; in this case the much slower pow(2, ..) can be used. */ - #define RGBE_CONTRIB_MANTBITS 5 - #define RGBE_CONTRIB_EXPBITS 5 - #define RGBE_CONTRIB_BINBITS 12 - - /* Signed maximum for mantissa */ - #define RGBE_CONTRIB_MANTMAX (1 << (RGBE_CONTRIB_MANTBITS - 1)) - /* Absolute overflow limit for mantissa (sign flips if exceeded) */ - #define RGBE_CONTRIB_MANTOVF ((RGBE_CONTRIB_MANTMAX << 1) - 1) - #define RGBE_CONTRIB_MIN (1.0 / (1L << (1 << RGBE_CONTRIB_EXPBITS))) - #define RGBE_CONTRIB_MAX 1.0 - #define RGBE_CONTRIB_BINMAX ((1U << RGBE_CONTRIB_BINBITS) - 1) - - - typedef union { - struct { - /* Compressed RGBE representation of binned contributions, - * consisting of a 4-bit mantissa + sign per RGB, a common 5-bit - * exponent, and a 12-bit bin number. */ - unsigned red: RGBE_CONTRIB_MANTBITS, - grn: RGBE_CONTRIB_MANTBITS, - blu: RGBE_CONTRIB_MANTBITS, - exp: RGBE_CONTRIB_EXPBITS, - bin: RGBE_CONTRIB_BINBITS; - }; - - uint32_t rgbeBin; - } RGBEContrib; - - - RGBEContrib RGBEContribEnc ( - double rgb [3], double norm [3], unsigned bin - ); - /* Normalise signed float RGB contributions against abs(norm), and encode - * to RGBE with associated bin. Issues warning if abs(norm) < abs(rgb). */ - - unsigned RGBEContribDec ( - RGBEContrib rgbeCont, double norm [3], double rgb [3] - ); - /* Decode RGBE contribution, returning signed floating point contributions - * in array rgb after denormalisation against abs(norm), and the associated - * bin as return value. */ -#endif diff --git a/rpict.c b/rpict.c index 9d5b5da..05f5be9 100644 --- a/rpict.c +++ b/rpict.c @@ -1,818 +1,781 @@ #ifndef lint -static const char RCSid[] = "$Id: rpict.c,v 2.94 2020/07/20 15:52:30 greg Exp $"; +static const char RCSid[] = "$Id: rpict.c,v 2.97 2021/12/03 18:10:48 greg Exp $"; #endif /* * rpict.c - routines and variables for picture generation. */ #include "copyright.h" #include #include "platform.h" #ifdef NON_POSIX #ifdef MINGW #include #endif #else #ifdef BSD #include #include #else #include #include #endif #endif #include #include #include "ray.h" #include "paths.h" #include "ambient.h" #include "view.h" #include "random.h" #include "paths.h" #include "hilbert.h" #include "pmapbias.h" #include "pmapdiag.h" #define RFTEMPLATE "rfXXXXXX" #ifndef SIGCONT #ifdef SIGIO /* XXX can we live without this? */ #define SIGCONT SIGIO #endif #endif CUBE thescene; /* our scene */ OBJECT nsceneobjs; /* number of objects in our scene */ int dimlist[MAXDIM]; /* sampling dimensions */ int ndims = 0; /* number of sampling dimensions */ int samplendx; /* sample index number */ void (*addobjnotify[])() = {ambnotify, NULL}; VIEW ourview = STDVIEW; /* view parameters */ int hresolu = 512; /* horizontal resolution */ int vresolu = 512; /* vertical resolution */ double pixaspect = 1.0; /* pixel aspect ratio */ int psample = 4; /* pixel sample size */ double maxdiff = .05; /* max. difference for interpolation */ double dstrpix = 0.67; /* square pixel distribution */ double mblur = 0.; /* motion blur parameter */ double dblur = 0.; /* depth-of-field blur parameter */ void (*trace)() = NULL; /* trace call */ int do_irrad = 0; /* compute irradiance? */ int rand_samp = 0; /* pure Monte Carlo sampling? */ double dstrsrc = 0.0; /* square source distribution */ double shadthresh = .05; /* shadow threshold */ double shadcert = .5; /* shadow certainty */ int directrelay = 1; /* number of source relays */ int vspretest = 512; /* virtual source pretest density */ int directvis = 1; /* sources visible? */ double srcsizerat = .25; /* maximum ratio source size/dist. */ COLOR cextinction = BLKCOLOR; /* global extinction coefficient */ COLOR salbedo = BLKCOLOR; /* global scattering albedo */ double seccg = 0.; /* global scattering eccentricity */ double ssampdist = 0.; /* scatter sampling distance */ double specthresh = .15; /* specular sampling threshold */ double specjitter = 1.; /* specular sampling jitter */ int backvis = 1; /* back face visibility */ int maxdepth = 7; /* maximum recursion depth */ double minweight = 1e-3; /* minimum ray weight */ char *ambfile = NULL; /* ambient file name */ COLOR ambval = BLKCOLOR; /* ambient value */ int ambvwt = 0; /* initial weight for ambient value */ double ambacc = 0.2; /* ambient accuracy */ int ambres = 64; /* ambient resolution */ int ambdiv = 512; /* ambient divisions */ int ambssamp = 128; /* ambient super-samples */ int ambounce = 0; /* ambient bounces */ char *amblist[AMBLLEN]; /* ambient include/exclude list */ int ambincl = -1; /* include == 1, exclude == 0 */ int ralrm = 0; /* seconds between reports */ double pctdone = 0.0; /* percentage done */ time_t tlastrept = 0L; /* time at last report */ time_t tstart; /* starting time */ #define MAXDIV 16 /* maximum sample size */ #define pixjitter() (.5+dstrpix*(.5-frandom())) int hres, vres; /* resolution for this frame */ static VIEW lastview; /* the previous view input */ static void report(int); static int nextview(FILE *fp); static void render(char *zfile, char *oldfile); static void fillscanline(COLOR *scanline, float *zline, char *sd, int xres, int y, int xstep); static void fillscanbar(COLOR *scanbar[], float *zbar[], int xres, int y, int ysize); static int fillsample(COLOR *colline, float *zline, int x, int y, int xlen, int ylen, int b); static double pixvalue(COLOR col, int x, int y); static int salvage(char *oldfile); static int pixnumber(int x, int y, int xres, int yres); void quit(int code) /* quit program */ { if (code) /* report status */ report(0); #ifndef NON_POSIX headclean(); /* delete header file */ pfclean(); /* clean up persist files */ #endif exit(code); } #ifndef NON_POSIX static void report(int dummy) /* report progress */ { char bcStat [128]; double u, s; #ifdef BSD struct rusage rubuf; #else double period = 1.0 / 60.0; struct tms tbuf; #endif tlastrept = time((time_t *)NULL); #ifdef BSD getrusage(RUSAGE_SELF, &rubuf); u = rubuf.ru_utime.tv_sec + rubuf.ru_utime.tv_usec*1e-6; s = rubuf.ru_stime.tv_sec + rubuf.ru_stime.tv_usec*1e-6; getrusage(RUSAGE_CHILDREN, &rubuf); u += rubuf.ru_utime.tv_sec + rubuf.ru_utime.tv_usec*1e-6; s += rubuf.ru_stime.tv_sec + rubuf.ru_stime.tv_usec*1e-6; #else times(&tbuf); #ifdef _SC_CLK_TCK period = 1.0 / sysconf(_SC_CLK_TCK); #endif u = ( tbuf.tms_utime + tbuf.tms_cutime ) * period; s = ( tbuf.tms_stime + tbuf.tms_cstime ) * period; #endif /* PMAP: Get photon map bias compensation statistics */ pmapBiasCompReport(bcStat); sprintf(errmsg, "%lu rays, %s%4.2f%% after %.3fu %.3fs %.3fr hours on %s (PID %d)\n", nrays, bcStat, pctdone, u*(1./3600.), s*(1./3600.), (tlastrept-tstart)*(1./3600.), myhostname(), getpid()); eputs(errmsg); #ifdef SIGCONT signal(SIGCONT, report); #endif } #else static void report(int dummy) /* report progress */ { char bcStat [128]; tlastrept = time((time_t *)NULL); /* PMAP: Get photon map bias compensation statistics */ pmapBiasCompReport(bcStat); sprintf(errmsg, "%lu rays, %s%4.2f%% after %5.4f hours\n", nrays, bcStat, pctdone, (tlastrept-tstart)/3600.0); eputs(errmsg); } #endif void rpict( /* generate image(s) */ int seq, char *pout, char *zout, char *prvr ) /* * If seq is greater than zero, then we will render a sequence of * images based on view parameter strings read from the standard input. * If pout is NULL, then all images will be sent to the standard ouput. * If seq is greater than zero and prvr is an integer, then it is the * frame number at which rendering should begin. Preceeding view parameter * strings will be skipped in the input. * If pout and prvr are the same, prvr is renamed to avoid overwriting. * Note that pout and zout should contain %d format specifications for * sequenced file naming. */ { char fbuf[128], fbuf2[128]; int npicts; char *cp; RESOLU rs; double pa; /* check sampling */ if (psample < 1) psample = 1; else if (psample > MAXDIV) { sprintf(errmsg, "pixel sampling reduced from %d to %d", psample, MAXDIV); error(WARNING, errmsg); psample = MAXDIV; } /* get starting frame */ if (seq <= 0) seq = 0; else if (prvr != NULL && isint(prvr)) { int rn; /* skip to specified view */ if ((rn = atoi(prvr)) < seq) error(USER, "recover frame less than start frame"); if (pout == NULL) error(USER, "missing output file specification"); for ( ; seq < rn; seq++) if (nextview(stdin) == EOF) error(USER, "unexpected EOF on view input"); setview(&ourview); prvr = fbuf; /* mark for renaming */ } if ((pout != NULL) & (prvr != NULL)) { sprintf(fbuf, pout, seq); if (!strcmp(prvr, fbuf)) { /* rename */ strcpy(fbuf2, fbuf); for (cp = fbuf2; *cp; cp++) ; while (cp > fbuf2 && !ISDIRSEP(cp[-1])) cp--; strcpy(cp, RFTEMPLATE); prvr = mktemp(fbuf2); if (rename(fbuf, prvr) < 0) { if (errno == ENOENT) { /* ghost file */ sprintf(errmsg, "new output file \"%s\"", fbuf); error(WARNING, errmsg); prvr = NULL; } else { /* serious error */ sprintf(errmsg, "cannot rename \"%s\" to \"%s\"", fbuf, prvr); error(SYSTEM, errmsg); } } } } npicts = 0; /* render sequence */ do { if (seq && nextview(stdin) == EOF) break; pctdone = 0.0; if (pout != NULL) { int myfd; close(1); /* reassign stdout */ sprintf(fbuf, pout, seq); tryagain: errno = 0; /* exclusive open */ if ((myfd = open(fbuf, O_WRONLY|O_CREAT|O_EXCL, 0666)) < 0) { if ((errno != EEXIST) | (prvr != NULL) || !strcmp(fbuf, pout)) { sprintf(errmsg, "cannot open output file \"%s\"", fbuf); error(SYSTEM, errmsg); } setview(&ourview); continue; /* don't clobber */ } if (myfd != 1) { unlink(fbuf); goto tryagain; /* leave it open */ } SET_FD_BINARY(1); dupheader(); } hres = hresolu; vres = vresolu; pa = pixaspect; if (prvr != NULL) { if (viewfile(prvr, &ourview, &rs) <= 0) { sprintf(errmsg, "cannot recover view parameters from \"%s\"", prvr); error(WARNING, errmsg); } else { pa = 0.0; hres = scanlen(&rs); vres = numscans(&rs); } } if ((cp = setview(&ourview)) != NULL) error(USER, cp); normaspect(viewaspect(&ourview), &pa, &hres, &vres); if (seq) { if (ralrm > 0) { fprintf(stderr, "FRAME %d:", seq); fprintview(&ourview, stderr); putc('\n', stderr); fflush(stderr); } printf("FRAME=%d\n", seq); } fputs(VIEWSTR, stdout); fprintview(&ourview, stdout); putchar('\n'); if ((pa < .99) | (pa > 1.01)) fputaspect(pa, stdout); fputnow(stdout); + fputprims(stdprims, stdout); fputformat(COLRFMT, stdout); putchar('\n'); if (zout != NULL) sprintf(cp=fbuf, zout, seq); else cp = NULL; render(cp, prvr); prvr = NULL; npicts++; } while (seq++); /* check that we did something */ if (npicts == 0) error(WARNING, "no output produced"); } static int nextview( /* get next view from fp */ FILE *fp ) { char linebuf[256]; lastview = ourview; while (fgets(linebuf, sizeof(linebuf), fp) != NULL) if (isview(linebuf) && sscanview(&ourview, linebuf) > 0) return(0); return(EOF); } static void render( /* render the scene */ char *zfile, char *oldfile ) { COLOR *scanbar[MAXDIV+1]; /* scanline arrays of pixel values */ float *zbar[MAXDIV+1]; /* z values */ char *sampdens; /* previous sample density */ int ypos; /* current scanline */ int ystep; /* current y step size */ int hstep; /* h step size */ int zfd; COLOR *colptr; float *zptr; int i; /* check for empty image */ if ((hres <= 0) | (vres <= 0)) { error(WARNING, "empty output picture"); fprtresolu(0, 0, stdout); return; } /* allocate scanlines */ for (i = 0; i <= psample; i++) { scanbar[i] = (COLOR *)malloc(hres*sizeof(COLOR)); if (scanbar[i] == NULL) goto memerr; } hstep = (psample*140+49)/99; /* quincunx sampling */ ystep = (psample*99+70)/140; if (hstep > 2) { i = hres/hstep + 2; if ((sampdens = (char *)malloc(i)) == NULL) goto memerr; while (i--) sampdens[i] = hstep; } else sampdens = NULL; /* open z-file */ if (zfile != NULL) { if ((zfd = open(zfile, O_WRONLY|O_CREAT, 0666)) == -1) { sprintf(errmsg, "cannot open z-file \"%s\"", zfile); error(SYSTEM, errmsg); } SET_FD_BINARY(zfd); for (i = 0; i <= psample; i++) { zbar[i] = (float *)malloc(hres*sizeof(float)); if (zbar[i] == NULL) goto memerr; } } else { zfd = -1; for (i = 0; i <= psample; i++) zbar[i] = NULL; } /* write out boundaries */ fprtresolu(hres, vres, stdout); /* recover file and compute first */ i = salvage(oldfile); if (i >= vres) goto alldone; if ((zfd != -1) & (i > 0) && lseek(zfd, (off_t)i*hres*sizeof(float), SEEK_SET) < 0) error(SYSTEM, "z-file seek error in render"); pctdone = 100.0*i/vres; if (ralrm > 0) /* report init stats */ report(0); #ifdef SIGCONT else signal(SIGCONT, report); #endif ypos = vres-1 - i; /* initialize sampling */ if (directvis) init_drawsources(psample); fillscanline(scanbar[0], zbar[0], sampdens, hres, ypos, hstep); /* compute scanlines */ for (ypos -= ystep; ypos > -ystep; ypos -= ystep) { /* bottom adjust? */ if (ypos < 0) { ystep += ypos; ypos = 0; } colptr = scanbar[ystep]; /* move base to top */ scanbar[ystep] = scanbar[0]; scanbar[0] = colptr; zptr = zbar[ystep]; zbar[ystep] = zbar[0]; zbar[0] = zptr; /* fill base line */ fillscanline(scanbar[0], zbar[0], sampdens, hres, ypos, hstep); /* fill bar */ fillscanbar(scanbar, zbar, hres, ypos, ystep); if (directvis) /* add bitty sources */ drawsources(scanbar, zbar, 0, hres, ypos, ystep); /* write it out */ #ifdef SIGCONT signal(SIGCONT, SIG_IGN); /* don't interrupt writes */ #endif for (i = ystep; i > 0; i--) { if (zfd != -1 && write(zfd, (char *)zbar[i], hres*sizeof(float)) < hres*sizeof(float)) goto writerr; if (fwritescan(scanbar[i], hres, stdout) < 0) goto writerr; } if (fflush(stdout) == EOF) goto writerr; /* record progress */ pctdone = 100.0*(vres-1-ypos)/vres; if (ralrm > 0 && time((time_t *)NULL) >= tlastrept+ralrm) report(0); #ifdef SIGCONT else signal(SIGCONT, report); #endif } /* clean up */ #ifdef SIGCONT signal(SIGCONT, SIG_IGN); #endif if (zfd != -1 && write(zfd, (char *)zbar[0], hres*sizeof(float)) < hres*sizeof(float)) goto writerr; fwritescan(scanbar[0], hres, stdout); if (fflush(stdout) == EOF) goto writerr; alldone: if (zfd != -1) { if (close(zfd) == -1) goto writerr; for (i = 0; i <= psample; i++) free((void *)zbar[i]); } for (i = 0; i <= psample; i++) free((void *)scanbar[i]); if (sampdens != NULL) free(sampdens); pctdone = 100.0; if (ralrm > 0) report(0); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif return; writerr: error(SYSTEM, "write error in render"); memerr: error(SYSTEM, "out of memory in render"); } static void fillscanline( /* fill scan at y */ COLOR *scanline, float *zline, char *sd, int xres, int y, int xstep ) { static int nc = 0; /* number of calls */ int bl = xstep, b = xstep; double z; int i; z = pixvalue(scanline[0], 0, y); if (zline) zline[0] = z; /* zig-zag start for quincunx pattern */ for (i = ++nc & 1 ? xstep : xstep/2; i < xres-1+xstep; i += xstep) { if (i >= xres) { xstep += xres-1-i; i = xres-1; } z = pixvalue(scanline[i], i, y); if (zline) zline[i] = z; if (sd) b = sd[0] > sd[1] ? sd[0] : sd[1]; if (i <= xstep) b = fillsample(scanline, zline, 0, y, i, 0, b/2); else b = fillsample(scanline+i-xstep, zline ? zline+i-xstep : (float *)NULL, i-xstep, y, xstep, 0, b/2); if (sd) *sd++ = nc & 1 ? bl : b; bl = b; } if (sd && nc & 1) *sd = bl; } static void fillscanbar( /* fill interior */ COLOR *scanbar[], float *zbar[], int xres, int y, int ysize ) { COLOR vline[MAXDIV+1]; float zline[MAXDIV+1]; int b = ysize; int i, j; for (i = 0; i < xres; i++) { copycolor(vline[0], scanbar[0][i]); copycolor(vline[ysize], scanbar[ysize][i]); if (zbar[0]) { zline[0] = zbar[0][i]; zline[ysize] = zbar[ysize][i]; } b = fillsample(vline, zbar[0] ? zline : (float *)NULL, i, y, 0, ysize, b/2); for (j = 1; j < ysize; j++) copycolor(scanbar[j][i], vline[j]); if (zbar[0]) for (j = 1; j < ysize; j++) zbar[j][i] = zline[j]; } } static int fillsample( /* fill interior points */ COLOR *colline, float *zline, int x, int y, int xlen, int ylen, int b ) { double ratio; double z; COLOR ctmp; int ncut; int len; if (xlen > 0) /* x or y length is zero */ len = xlen; else len = ylen; if (len <= 1) /* limit recursion */ return(0); if (b > 0 || (zline && 2.*fabs(zline[0]-zline[len]) > maxdiff*(zline[0]+zline[len])) || bigdiff(colline[0], colline[len], maxdiff)) { z = pixvalue(colline[len>>1], x + (xlen>>1), y + (ylen>>1)); if (zline) zline[len>>1] = z; ncut = 1; } else { /* interpolate */ copycolor(colline[len>>1], colline[len]); ratio = (double)(len>>1) / len; scalecolor(colline[len>>1], ratio); if (zline) zline[len>>1] = zline[len] * ratio; ratio = 1.0 - ratio; copycolor(ctmp, colline[0]); scalecolor(ctmp, ratio); addcolor(colline[len>>1], ctmp); if (zline) zline[len>>1] += zline[0] * ratio; ncut = 0; } /* recurse */ ncut += fillsample(colline, zline, x, y, xlen>>1, ylen>>1, (b-1)/2); ncut += fillsample(colline+(len>>1), zline ? zline+(len>>1) : (float *)NULL, x+(xlen>>1), y+(ylen>>1), xlen-(xlen>>1), ylen-(ylen>>1), b/2); return(ncut); } static double pixvalue( /* compute pixel value */ COLOR col, /* returned color */ int x, /* pixel position */ int y ) { - extern void SDsquare2disk(double ds[2], double seedx, double seedy); RAY thisray; FVECT lorg, ldir; - double hpos, vpos, vdist, lmax; + double hpos, vpos, lmax; int i; /* compute view ray */ setcolor(col, 0.0, 0.0, 0.0); hpos = (x+pixjitter())/hres; vpos = (y+pixjitter())/vres; if ((thisray.rmax = viewray(thisray.rorg, thisray.rdir, &ourview, hpos, vpos)) < -FTINY) return(0.0); - - vdist = ourview.vdist; /* set pixel index */ samplendx = pixnumber(x,y,hres,vres); /* optional motion blur */ if (lastview.type && mblur > FTINY && (lmax = viewray(lorg, ldir, &lastview, hpos, vpos)) >= -FTINY) { double d = mblur*(.5-urand(281+samplendx)); thisray.rmax = (1.-d)*thisray.rmax + d*lmax; for (i = 3; i--; ) { thisray.rorg[i] = (1.-d)*thisray.rorg[i] + d*lorg[i]; thisray.rdir[i] = (1.-d)*thisray.rdir[i] + d*ldir[i]; } if (normalize(thisray.rdir) == 0.0) return(0.0); - vdist = (1.-d)*vdist + d*lastview.vdist; - } - /* optional depth-of-field */ - if (dblur > FTINY) { - double vc, df[2]; - /* random point on disk */ - SDsquare2disk(df, frandom(), frandom()); - df[0] *= .5*dblur; - df[1] *= .5*dblur; - if ((ourview.type == VT_PER) | (ourview.type == VT_PAR)) { - double adj = 1.0; - if (ourview.type == VT_PER) - adj /= DOT(thisray.rdir, ourview.vdir); - df[0] /= sqrt(ourview.hn2); - df[1] /= sqrt(ourview.vn2); - for (i = 3; i--; ) { - vc = ourview.vp[i] + adj*vdist*thisray.rdir[i]; - thisray.rorg[i] += df[0]*ourview.hvec[i] + - df[1]*ourview.vvec[i] ; - thisray.rdir[i] = vc - thisray.rorg[i]; - } - } else { /* non-standard view case */ - double dfd = PI/4.*dblur*(.5 - frandom()); - if ((ourview.type != VT_ANG) & (ourview.type != VT_PLS)) { - if (ourview.type != VT_CYL) - df[0] /= sqrt(ourview.hn2); - df[1] /= sqrt(ourview.vn2); - } - for (i = 3; i--; ) { - vc = ourview.vp[i] + vdist*thisray.rdir[i]; - thisray.rorg[i] += df[0]*ourview.hvec[i] + - df[1]*ourview.vvec[i] + - dfd*ourview.vdir[i] ; - thisray.rdir[i] = vc - thisray.rorg[i]; - } - } - if (normalize(thisray.rdir) == 0.0) - return(0.0); } + /* depth-of-field */ + if (!jitteraperture(thisray.rorg, thisray.rdir, &ourview, dblur)) + return(0.0); rayorigin(&thisray, PRIMARY, NULL, NULL); rayvalue(&thisray); /* trace ray */ copycolor(col, thisray.rcol); /* return color */ return(raydistance(&thisray)); /* return distance */ } static int salvage( /* salvage scanlines from killed program */ char *oldfile ) { COLR *scanline; FILE *fp; int x, y; if (oldfile == NULL) goto gotzip; if ((fp = fopen(oldfile, "r")) == NULL) { sprintf(errmsg, "cannot open recover file \"%s\"", oldfile); error(WARNING, errmsg); goto gotzip; } SET_FILE_BINARY(fp); /* discard header */ getheader(fp, NULL, NULL); /* get picture size */ if (!fscnresolu(&x, &y, fp)) { sprintf(errmsg, "bad recover file \"%s\" - not removed", oldfile); error(WARNING, errmsg); fclose(fp); goto gotzip; } if ((x != hres) | (y != vres)) { sprintf(errmsg, "resolution mismatch in recover file \"%s\"", oldfile); error(USER, errmsg); } scanline = (COLR *)malloc(hres*sizeof(COLR)); if (scanline == NULL) error(SYSTEM, "out of memory in salvage"); for (y = 0; y < vres; y++) { if (freadcolrs(scanline, hres, fp) < 0) break; if (fwritecolrs(scanline, hres, stdout) < 0) goto writerr; } if (fflush(stdout) == EOF) goto writerr; free((void *)scanline); fclose(fp); unlink(oldfile); return(y); gotzip: if (fflush(stdout) == EOF) error(SYSTEM, "error writing picture header"); return(0); writerr: sprintf(errmsg, "write error during recovery of \"%s\"", oldfile); error(SYSTEM, errmsg); return -1; /* pro forma return */ } static int pixnumber( /* compute pixel index (screen door) */ int x, int y, int xres, int yres ) { unsigned nbits = 0; bitmask_t coord[2]; if (xres < yres) xres = yres; while (xres > 0) { xres >>= 1; ++nbits; } coord[0] = x; coord[1] = y; return ((int)hilbert_c2i(2, nbits, coord)); } diff --git a/rv2.c b/rv2.c index b265fe4..9807e62 100644 --- a/rv2.c +++ b/rv2.c @@ -1,892 +1,893 @@ #ifndef lint -static const char RCSid[] = "$Id: rv2.c,v 2.72 2021/02/12 00:41:19 greg Exp $"; +static const char RCSid[] = "$Id: rv2.c,v 2.73 2021/11/14 17:30:02 greg Exp $"; #endif /* * rv2.c - command routines used in tracing a view. * * External symbols declared in rpaint.h */ #include "copyright.h" #include #include #include "platform.h" #include "rtprocess.h" /* win_popen() */ #include "paths.h" #include "ray.h" #include "ambient.h" #include "otypes.h" #include "otspecial.h" #include "rpaint.h" extern int psample; /* pixel sample size */ extern double maxdiff; /* max. sample difference */ #define CTRL(c) ((c)-'@') #define sscanvec(s,v) (sscanf(s,FVFORMAT,v,v+1,v+2)==3) extern char rifname[128]; /* rad input file name */ extern char *progname; extern char *octname; void getframe( /* get a new frame */ char *s ) { if (getrect(s, &pframe) < 0) return; pdepth = 0; } void getrepaint( /* get area and repaint */ char *s ) { RECT box; if (getrect(s, &box) < 0) return; paintrect(&ptrunk, &box); } void getview( /* get/show/save view parameters */ char *s ) { FILE *fp; char buf[128]; char *fname; int change = 0; VIEW nv = ourview; while (isspace(*s)) s++; if (*s == '-') { /* command line parameters */ if (sscanview(&nv, s)) newview(&nv); else error(COMMAND, "bad view option(s)"); return; } if (nextword(buf, sizeof(buf), s) != NULL) { /* write to a file */ if ((fname = getpath(buf, NULL, 0)) == NULL || (fp = fopen(fname, "a")) == NULL) { sprintf(errmsg, "cannot open \"%s\"", buf); error(COMMAND, errmsg); return; } fputs(progname, fp); fprintview(&ourview, fp); fputs(sskip(s), fp); putc('\n', fp); fclose(fp); return; } sprintf(buf, "view type (%c): ", ourview.type); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (buf[0] && buf[0] != ourview.type) { nv.type = buf[0]; change++; } sprintf(buf, "view point (%.6g %.6g %.6g): ", ourview.vp[0], ourview.vp[1], ourview.vp[2]); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (sscanvec(buf, nv.vp)) change++; sprintf(buf, "view direction (%.6g %.6g %.6g): ", ourview.vdir[0]*ourview.vdist, ourview.vdir[1]*ourview.vdist, ourview.vdir[2]*ourview.vdist); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (sscanvec(buf, nv.vdir)) { nv.vdist = 1.; change++; } sprintf(buf, "view up (%.6g %.6g %.6g): ", ourview.vup[0], ourview.vup[1], ourview.vup[2]); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (sscanvec(buf, nv.vup)) change++; sprintf(buf, "view horiz and vert size (%.6g %.6g): ", ourview.horiz, ourview.vert); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (sscanf(buf, "%lf %lf", &nv.horiz, &nv.vert) == 2) change++; sprintf(buf, "fore and aft clipping plane (%.6g %.6g): ", ourview.vfore, ourview.vaft); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (sscanf(buf, "%lf %lf", &nv.vfore, &nv.vaft) == 2) change++; sprintf(buf, "view shift and lift (%.6g %.6g): ", ourview.hoff, ourview.voff); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == CTRL('C')) return; if (sscanf(buf, "%lf %lf", &nv.hoff, &nv.voff) == 2) change++; if (change) newview(&nv); } void lastview( /* return to a previous view */ char *s ) { char buf[128]; char *fname; int success; VIEW nv; /* get parameters from a file */ if (nextword(buf, sizeof(buf), s) != NULL) { nv = stdview; if ((fname = getpath(buf, "", R_OK)) == NULL || (success = viewfile(fname, &nv, NULL)) == -1) { sprintf(errmsg, "cannot open \"%s\"", buf); error(COMMAND, errmsg); return; } if (!success) error(COMMAND, "wrong file format"); else newview(&nv); return; } if (oldview.type == 0) { /* no old view! */ error(COMMAND, "no previous view"); return; } nv = ourview; ourview = oldview; oldview = nv; newimage(NULL); } void saveview( /* save view to rad file */ char *s ) { char view[64]; char *fname; FILE *fp; if (*atos(view, sizeof(view), s)) { if (isint(view)) { error(COMMAND, "cannot write view by number"); return; } s = sskip(s); } if (nextword(rifname, sizeof(rifname), s) == NULL && !rifname[0]) { error(COMMAND, "no previous rad file"); return; } if ((fname = getpath(rifname, NULL, 0)) == NULL || (fp = fopen(fname, "a")) == NULL) { sprintf(errmsg, "cannot open \"%s\"", rifname); error(COMMAND, errmsg); return; } fputs("view= ", fp); fputs(view, fp); fprintview(&ourview, fp); putc('\n', fp); fclose(fp); } void loadview( /* load view from rad file */ char *s ) { char buf[512]; char *fname; FILE *fp; VIEW nv; strcpy(buf, "rad -n -s -V -v "); if (sscanf(s, "%s", buf+strlen(buf)) == 1) s = sskip(s); else strcat(buf, "1"); if (nextword(rifname, sizeof(rifname), s) == NULL && !rifname[0]) { error(COMMAND, "no previous rad file"); return; } if ((fname = getpath(rifname, "", R_OK)) == NULL) { sprintf(errmsg, "cannot access \"%s\"", rifname); error(COMMAND, errmsg); return; } sprintf(buf+strlen(buf), " %s", fname); if ((fp = popen(buf, "r")) == NULL) { error(COMMAND, "cannot run rad"); return; } buf[0] = '\0'; fgets(buf, sizeof(buf), fp); pclose(fp); nv = stdview; if (!sscanview(&nv, buf)) { error(COMMAND, "rad error -- no such view?"); return; } newview(&nv); } void getaim( /* aim camera */ char *s ) { VIEW nv = ourview; double zfact; if (getinterest(s, 1, nv.vdir, &zfact) < 0) return; zoomview(&nv, zfact); newview(&nv); } void getfocus( /* set focus distance */ char *s ) { char buf[64]; double dist; if (sscanf(s, "%lf", &dist) < 1) { int x, y; RAY thisray; if (dev->getcur == NULL) return; (*dev->comout)("Pick focus point\n"); if ((*dev->getcur)(&x, &y) == ABORT) return; if ((thisray.rmax = viewray(thisray.rorg, thisray.rdir, &ourview, (x+.5)/hresolu, (y+.5)/vresolu)) < -FTINY) { error(COMMAND, "not on image"); return; } rayorigin(&thisray, PRIMARY, NULL, NULL); if (!localhit(&thisray, &thescene)) { error(COMMAND, "not a local object"); return; } dist = thisray.rot; } else if (dist <= .0) { error(COMMAND, "focus distance must be positive"); return; } ourview.vdist = dist; sprintf(buf, "Focus distance set to %f\n", dist); (*dev->comout)(buf); } void getmove( /* move camera */ char *s ) { FVECT vc; double mag; if (getinterest(s, 0, vc, &mag) < 0) return; moveview(0.0, 0.0, mag, vc); } void getrotate( /* rotate camera */ char *s ) { VIEW nv = ourview; double angle, elev, zfact; elev = 0.0; zfact = 1.0; if (sscanf(s, "%lf %lf %lf", &angle, &elev, &zfact) < 1) { error(COMMAND, "missing angle"); return; } spinvector(nv.vdir, ourview.vdir, ourview.vup, angle*(PI/180.)); if (elev != 0.0) geodesic(nv.vdir, nv.vdir, nv.vup, elev*(PI/180.), GEOD_RAD); zoomview(&nv, zfact); newview(&nv); } void getpivot( /* pivot viewpoint */ char *s ) { FVECT vc; double angle, elev, mag; elev = 0.0; if (sscanf(s, "%lf %lf", &angle, &elev) < 1) { error(COMMAND, "missing angle"); return; } if (getinterest(sskip2(s,2), 0, vc, &mag) < 0) return; moveview(angle, elev, mag, vc); } void getorigin( /* origin viewpoint */ char *s ) { VIEW nv = ourview; double d; /* get new view origin */ if (sscanf(s, "%lf %lf", &d, &d) == 1) { /* just moving some distance */ VSUM(nv.vp, nv.vp, nv.vdir, d); } else if (!sscanvec(s, nv.vp)) { int x, y; /* need to pick origin */ RAY thisray; if (dev->getcur == NULL) return; (*dev->comout)("Pick point on surface for new origin\n"); if ((*dev->getcur)(&x, &y) == ABORT) return; if ((thisray.rmax = viewray(thisray.rorg, thisray.rdir, &ourview, (x+.5)/hresolu, (y+.5)/vresolu)) < -FTINY) { error(COMMAND, "not on image"); return; } rayorigin(&thisray, PRIMARY, NULL, NULL); if (!localhit(&thisray, &thescene)) { error(COMMAND, "not a local object"); return; } if (thisray.rod < 0.0) /* don't look through other side */ flipsurface(&thisray); VSUM(nv.vp, thisray.rop, thisray.ron, 20.0*FTINY); VCOPY(nv.vdir, thisray.ron); } else if (!sscanvec(sskip2(s,3), nv.vdir) || normalize(nv.vdir) == 0.0) VCOPY(nv.vdir, ourview.vdir); d = DOT(nv.vdir, nv.vup); /* need different up vector? */ if (d*d >= 1.-2.*FTINY) { int i; nv.vup[0] = nv.vup[1] = nv.vup[2] = 0.0; for (i = 3; i--; ) if (nv.vdir[i]*nv.vdir[i] < 0.34) break; nv.vup[i] = 1.; } newview(&nv); } void getexposure( /* get new exposure */ char *s ) { char buf[128]; char *cp; int x, y; PNODE *p = &ptrunk; int adapt = 0; double e = 1.0; for (cp = s; isspace(*cp); cp++) ; if (*cp == '@') { adapt++; while (isspace(*++cp)) ; } if (*cp == '\0') { /* normalize to point */ if (dev->getcur == NULL) return; (*dev->comout)("Pick point for exposure\n"); if ((*dev->getcur)(&x, &y) == ABORT) return; p = findrect(x, y, &ptrunk, -1); } else { if (*cp == '=') { /* absolute setting */ p = NULL; e = 1.0/exposure; for (cp++; isspace(*cp); cp++) ; if (*cp == '\0') { /* interactive */ sprintf(buf, "exposure (%f): ", exposure); (*dev->comout)(buf); (*dev->comin)(buf, NULL); for (cp = buf; isspace(*cp); cp++) ; if (*cp == '\0') return; } } if (*cp == '+' || *cp == '-') /* f-stops */ e *= pow(2.0, atof(cp)); else /* multiplier */ e *= atof(cp); } if (p != NULL) { /* relative setting */ compavg(p); if (bright(p->v) < 1e-15) { error(COMMAND, "cannot normalize to zero"); return; } if (adapt) e *= 106./pow(1.219+pow(luminance(p->v)/exposure,.4),2.5)/exposure; else e *= 0.5 / bright(p->v); } if (e <= FTINY || fabs(1.0 - e) <= FTINY) return; scalepict(&ptrunk, e); exposure *= e; redraw(); } typedef union {int i; double d; COLOR C;} *MyUptr; int getparam( /* get variable from user */ char *str, char *dsc, int typ, void *p ) { MyUptr ptr = (MyUptr)p; int i0; double d0, d1, d2; char buf[48]; switch (typ) { case 'i': /* integer */ if (sscanf(str, "%d", &i0) != 1) { (*dev->comout)(dsc); sprintf(buf, " (%d): ", ptr->i); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (sscanf(buf, "%d", &i0) != 1) return(0); } if (ptr->i == i0) return(0); ptr->i = i0; break; case 'r': /* real */ if (sscanf(str, "%lf", &d0) != 1) { (*dev->comout)(dsc); sprintf(buf, " (%.6g): ", ptr->d); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (sscanf(buf, "%lf", &d0) != 1) return(0); } if (FABSEQ(ptr->d, d0)) return(0); ptr->d = d0; break; case 'b': /* boolean */ if (sscanf(str, "%1s", buf) != 1) { (*dev->comout)(dsc); sprintf(buf, "? (%c): ", ptr->i ? 'y' : 'n'); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (buf[0] == '\0') return(0); } if (strchr("yY+1tTnN-0fF", buf[0]) == NULL) return(0); i0 = strchr("yY+1tT", buf[0]) != NULL; if (ptr->i == i0) return(0); ptr->i = i0; break; case 'C': /* color */ if (sscanf(str, "%lf %lf %lf", &d0, &d1, &d2) != 3) { (*dev->comout)(dsc); sprintf(buf, " (%.6g %.6g %.6g): ", colval(ptr->C,RED), colval(ptr->C,GRN), colval(ptr->C,BLU)); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (sscanf(buf, "%lf %lf %lf", &d0, &d1, &d2) != 3) return(0); } if (FABSEQ(colval(ptr->C,RED), d0) && FABSEQ(colval(ptr->C,GRN), d1) && FABSEQ(colval(ptr->C,BLU), d2)) return(0); setcolor(ptr->C, d0, d1, d2); break; default: return(0); /* shouldn't happen */ } newparam++; return(1); } void setparam( /* get/set program parameter */ char *s ) { int prev_newp = newparam; char buf[128]; if (s[0] == '\0') { (*dev->comout)( "aa ab ad ar as av aw b bv dc dv dj ds dt i lr lw me ma mg ms ps pt ss st u: "); (*dev->comin)(buf, NULL); s = buf; } switch (s[0]) { case 'u': /* uncorrelated sampling */ getparam(s+1, "uncorrelated sampling", 'b', (void *)&rand_samp); break; case 'l': /* limit */ switch (s[1]) { case 'w': /* weight */ getparam(s+2, "limit weight", 'r', (void *)&minweight); break; case 'r': /* reflection */ getparam(s+2, "limit reflection", 'i', (void *)&maxdepth); break; default: goto badparam; } break; case 'd': /* direct */ switch (s[1]) { case 'j': /* jitter */ getparam(s+2, "direct jitter", 'r', (void *)&dstrsrc); break; case 'c': /* certainty */ getparam(s+2, "direct certainty", 'r', (void *)&shadcert); break; case 't': /* threshold */ getparam(s+2, "direct threshold", 'r', (void *)&shadthresh); break; case 'v': /* visibility */ getparam(s+2, "direct visibility", 'b', (void *)&directvis); break; case 's': /* sampling */ getparam(s+2, "direct sampling", 'r', (void *)&srcsizerat); break; default: goto badparam; } break; case 'b': /* back faces or black and white */ switch (s[1]) { case 'v': /* back face visibility */ getparam(s+2, "back face visibility", 'b', (void *)&backvis); break; case '\0': /* black and white */ case ' ': case 'y': case 'Y': case 't': case 'T': case '1': case '+': case 'n': case 'N': case 'f': case 'F': case '0': case '-': getparam(s+1, "black and white", 'b', (void *)&greyscale); newparam = prev_newp; break; default: goto badparam; } break; case 'i': /* irradiance */ getparam(s+1, "irradiance", 'b', (void *)&do_irrad); break; case 'a': /* ambient */ switch (s[1]) { case 'v': /* value */ getparam(s+2, "ambient value", 'C', (void *)ambval); break; case 'w': /* weight */ getparam(s+2, "ambient value weight", 'i', (void *)&ambvwt); break; case 'a': /* accuracy */ if (getparam(s+2, "ambient accuracy", 'r', (void *)&ambacc)) setambacc(ambacc); break; case 'd': /* divisions */ getparam(s+2, "ambient divisions", 'i', (void *)&ambdiv); break; case 's': /* samples */ getparam(s+2, "ambient super-samples", 'i', (void *)&ambssamp); break; case 'b': /* bounces */ getparam(s+2, "ambient bounces", 'i', (void *)&ambounce); break; case 'r': if (getparam(s+2, "ambient resolution", 'i', (void *)&ambres)) setambres(ambres); break; default: goto badparam; } break; case 'm': /* medium */ switch (s[1]) { case 'e': /* extinction coefficient */ getparam(s+2, "extinction coefficient", 'C', (void *)cextinction); break; case 'a': /* scattering albedo */ getparam(s+2, "scattering albedo", 'C', (void *)salbedo); break; case 'g': /* scattering eccentricity */ getparam(s+2, "scattering eccentricity", 'r', (void *)&seccg); break; case 's': /* sampling distance */ getparam(s+2, "mist sampling distance", 'r', (void *)&ssampdist); break; default: goto badparam; } break; case 'p': /* pixel */ switch (s[1]) { case 's': /* sample */ if (getparam(s+2, "pixel sample", 'i', (void *)&psample)) pdepth = 0; break; case 't': /* threshold */ if (getparam(s+2, "pixel threshold", 'r', (void *)&maxdiff)) pdepth = 0; break; default: goto badparam; } newparam = prev_newp; break; case 's': /* specular */ switch (s[1]) { case 's': /* sampling */ getparam(s+2, "specular sampling", 'r', (void *)&specjitter); break; case 't': /* threshold */ getparam(s+2, "specular threshold", 'r', (void *)&specthresh); break; default: goto badparam; } break; case '\0': /* nothing */ break; default:; badparam: *sskip(s) = '\0'; sprintf(errmsg, "%s: unknown variable", s); error(COMMAND, errmsg); break; } } void traceray( /* trace a single ray */ char *s ) { RAY thisray; char buf[512]; thisray.rmax = 0.0; if (!sscanvec(s, thisray.rorg) || !sscanvec(sskip2(s,3), thisray.rdir)) { int x, y; if (dev->getcur == NULL) return; (*dev->comout)("Pick ray\n"); if ((*dev->getcur)(&x, &y) == ABORT) return; if ((thisray.rmax = viewray(thisray.rorg, thisray.rdir, &ourview, (x+.5)/hresolu, (y+.5)/vresolu)) < -FTINY) { error(COMMAND, "not on image"); return; } } else if (normalize(thisray.rdir) == 0.0) { error(COMMAND, "zero ray direction"); return; } ray_trace(&thisray); if (thisray.ro == NULL) (*dev->comout)("ray hit nothing"); else { OBJREC *mat = NULL; OBJREC *mod = NULL; char matspec[256]; OBJREC *ino; matspec[0] = '\0'; if (thisray.ro->omod != OVOID) { mod = objptr(thisray.ro->omod); mat = findmaterial(thisray.ro); } if (thisray.rod < 0.0) strcpy(matspec, "back of "); if (mod != NULL) { strcat(matspec, mod->oname); if (mat != mod && mat != NULL) sprintf(matspec+strlen(matspec), " (%s)", mat->oname); } else strcat(matspec, VOIDID); sprintf(buf, "ray hit %s %s \"%s\"", matspec, ofun[thisray.ro->otype].funame, thisray.ro->oname); if ((ino = objptr(thisray.robj)) != thisray.ro) sprintf(buf+strlen(buf), " in %s \"%s\"", ofun[ino->otype].funame, ino->oname); (*dev->comout)(buf); (*dev->comin)(buf, NULL); if (thisray.rot >= FHUGE*.99) (*dev->comout)("at infinity"); else { sprintf(buf, "at (%.6g %.6g %.6g) (%.6g)", thisray.rop[0], thisray.rop[1], thisray.rop[2], raydistance(&thisray)); (*dev->comout)(buf); } (*dev->comin)(buf, NULL); sprintf(buf, "value (%.5g %.5g %.5g) (%.3gL)", colval(thisray.rcol,RED), colval(thisray.rcol,GRN), colval(thisray.rcol,BLU), luminance(thisray.rcol)); (*dev->comout)(buf); } (*dev->comin)(buf, NULL); } void writepict( /* write the picture to a file */ char *s ) { static char buf[128]; char *fname; FILE *fp; COLR *scanline; int y; /* XXX relies on words.c 2.11 behavior */ if (nextword(buf, sizeof(buf), s) == NULL && !buf[0]) { error(COMMAND, "no file"); return; } if ((fname = getpath(buf, NULL, 0)) == NULL || (fp = fopen(fname, "w")) == NULL) { sprintf(errmsg, "cannot open \"%s\"", buf); error(COMMAND, errmsg); return; } SET_FILE_BINARY(fp); (*dev->comout)("writing \""); (*dev->comout)(fname); (*dev->comout)("\"...\n"); /* write header */ newheader("RADIANCE", fp); fputs(progname, fp); fprintview(&ourview, fp); if (octname != NULL) fprintf(fp, " %s\n", octname); else putc('\n', fp); fprintf(fp, "SOFTWARE= %s\n", VersionID); fputnow(fp); if (exposure != 1.0) fputexpos(exposure, fp); if (dev->pixaspect != 1.0) fputaspect(dev->pixaspect, fp); + fputprims(stdprims, fp); fputformat(COLRFMT, fp); putc('\n', fp); fprtresolu(hresolu, vresolu, fp); scanline = (COLR *)malloc(hresolu*sizeof(COLR)); if (scanline == NULL) { error(COMMAND, "not enough memory!"); fclose(fp); unlink(fname); return; } for (y = vresolu-1; y >= 0; y--) { getpictcolrs(y, scanline, &ptrunk, hresolu, vresolu); if (fwritecolrs(scanline, hresolu, fp) < 0) break; } free((void *)scanline); if (fclose(fp) < 0) error(COMMAND, "write error"); } diff --git a/source.c b/source.c index f5b8df6..7103c63 100644 --- a/source.c +++ b/source.c @@ -1,784 +1,784 @@ #ifndef lint -static const char RCSid[] = "$Id: source.c,v 2.76 2021/02/01 16:19:49 greg Exp $"; +static const char RCSid[] = "$Id: source.c,v 2.77 2021/11/24 19:08:51 greg Exp $"; #endif /* * source.c - routines dealing with illumination sources. * * External symbols declared in source.h */ #include "ray.h" #include "otypes.h" #include "otspecial.h" #include "rtotypes.h" #include "source.h" #include "random.h" #include "pmapsrc.h" #include "pmapmat.h" #ifndef MAXSSAMP #define MAXSSAMP 16 /* maximum samples per ray */ #endif /* * Structures used by direct() */ typedef struct { int sno; /* source number */ FVECT dir; /* source direction */ COLOR coef; /* material coefficient */ COLOR val; /* contribution */ } CONTRIB; /* direct contribution */ typedef struct { int sndx; /* source index (to CONTRIB array) */ float brt; /* brightness (for comparison) */ } CNTPTR; /* contribution pointer */ static CONTRIB *srccnt; /* source contributions in direct() */ static CNTPTR *cntord; /* source ordering in direct() */ static int maxcntr = 0; /* size of contribution arrays */ static int cntcmp(const void *p1, const void *p2); void marksources(void) /* find and mark source objects */ { int foundsource = 0; int i; OBJREC *o, *m; int ns; /* initialize dispatch table */ initstypes(); /* find direct sources */ for (i = 0; i < nsceneobjs; i++) { o = objptr(i); if (!issurface(o->otype) || o->omod == OVOID) continue; /* find material */ m = findmaterial(o); if (m == NULL) continue; if (m->otype == MAT_CLIP) { markclip(m); /* special case for antimatter */ continue; } if (!islight(m->otype)) continue; /* not source modifier */ if (m->oargs.nfargs != (m->otype == MAT_GLOW ? 4 : m->otype == MAT_SPOT ? 7 : 3)) objerror(m, USER, "bad # arguments"); if (m->oargs.farg[0] <= FTINY && (m->oargs.farg[1] <= FTINY) & (m->oargs.farg[2] <= FTINY)) continue; /* don't bother */ if (m->otype == MAT_GLOW && o->otype != OBJ_SOURCE && m->oargs.farg[3] <= FTINY) { foundsource += (ambounce > 0); continue; /* don't track these */ } if (sfun[o->otype].of == NULL || sfun[o->otype].of->setsrc == NULL) objerror(o, USER, "illegal material"); if ((ns = newsource()) < 0) goto memerr; setsource(&source[ns], o); if (m->otype == MAT_GLOW) { source[ns].sflags |= SPROX; source[ns].sl.prox = m->oargs.farg[3]; if (source[ns].sflags & SDISTANT) { source[ns].sflags |= SSKIP; foundsource += (ambounce > 0); } } else if (m->otype == MAT_SPOT) { if (source[ns].sflags & SDISTANT) objerror(o, WARNING, "distant source is a spotlight"); source[ns].sflags |= SSPOT; if ((source[ns].sl.s = makespot(m)) == NULL) goto memerr; if (source[ns].sflags & SFLAT && !checkspot(source[ns].sl.s,source[ns].snorm)) { objerror(o, WARNING, "invalid spotlight direction"); source[ns].sflags |= SSKIP; } } foundsource += !(source[ns].sflags & SSKIP); } if (!foundsource) { error(WARNING, "no light sources found"); return; } #if SHADCACHE for (ns = 0; ns < nsources; ns++) /* initialize obstructor cache */ initobscache(ns); #endif /* PMAP: disable virtual sources */ if (!photonMapping) markvirtuals(); /* find and add virtual sources */ /* allocate our contribution arrays */ maxcntr = nsources + MAXSPART; /* start with this many */ srccnt = (CONTRIB *)malloc(maxcntr*sizeof(CONTRIB)); cntord = (CNTPTR *)malloc(maxcntr*sizeof(CNTPTR)); if ((srccnt != NULL) & (cntord != NULL)) return; memerr: error(SYSTEM, "out of memory in marksources"); } void distantsources(void) /* only mark distant sources */ { int i; OBJREC *o, *m; int ns; /* initialize dispatch table */ initstypes(); /* sources needed for sourcehit() */ for (i = 0; i < nsceneobjs; i++) { o = objptr(i); if ((o->otype != OBJ_SOURCE) | (o->omod == OVOID)) continue; /* find material */ m = findmaterial(o); if (m == NULL) continue; if (!islight(m->otype)) continue; /* not source modifier */ if (m->oargs.nfargs != (m->otype == MAT_GLOW ? 4 : m->otype == MAT_SPOT ? 7 : 3)) objerror(m, USER, "bad # arguments"); if (m->oargs.farg[0] <= FTINY && (m->oargs.farg[1] <= FTINY) & (m->oargs.farg[2] <= FTINY)) continue; /* don't bother */ if (sfun[o->otype].of == NULL || sfun[o->otype].of->setsrc == NULL) objerror(o, USER, "illegal material"); if ((ns = newsource()) < 0) error(SYSTEM, "out of memory in distantsources"); setsource(&source[ns], o); if (m->otype == MAT_GLOW) { source[ns].sflags |= SPROX|SSKIP; source[ns].sl.prox = m->oargs.farg[3]; } else if (m->otype == MAT_SPOT) objerror(o, WARNING, "distant source is a spotlight"); } } void freesources(void) /* free all source structures */ { if (nsources > 0) { #if SHADCACHE while (nsources--) freeobscache(&source[nsources]); #endif free((void *)source); source = NULL; nsources = 0; } markclip(NULL); if (maxcntr <= 0) return; free((void *)srccnt); srccnt = NULL; free((void *)cntord); cntord = NULL; maxcntr = 0; } int srcray( /* send a ray to a source, return domega */ RAY *sr, /* returned source ray */ RAY *r, /* ray which hit object */ SRCINDEX *si /* source sample index */ ) { double d; /* distance to source */ SRCREC *srcp; rayorigin(sr, SHADOW, r, NULL); /* ignore limits */ if (r == NULL) sr->rmax = 0.0; while ((d = nextssamp(sr, si)) != 0.0) { sr->rsrc = si->sn; /* remember source */ srcp = source + si->sn; if (srcp->sflags & SDISTANT) { if (srcp->sflags & SSPOT && spotout(sr, srcp->sl.s)) continue; return(1); /* sample OK */ } /* local source */ /* check proximity */ if (srcp->sflags & SPROX && d > srcp->sl.prox) continue; /* check angle */ if (srcp->sflags & SSPOT) { if (spotout(sr, srcp->sl.s)) continue; /* adjust solid angle */ si->dom *= d*d; d += srcp->sl.s->flen; si->dom /= d*d; } return(1); /* sample OK */ } return(0); /* no more samples */ } void srcvalue( /* punch ray to source and compute value */ RAY *r ) { SRCREC *sp; sp = &source[r->rsrc]; if (sp->sflags & SVIRTUAL) { /* virtual source */ /* check intersection */ if (!(*ofun[sp->so->otype].funp)(sp->so, r)) return; if (!rayshade(r, r->ro->omod)) /* compute contribution */ goto nomat; rayparticipate(r); return; } /* compute intersection */ if (sp->sflags & SDISTANT ? sourcehit(r) : (*ofun[sp->so->otype].funp)(sp->so, r)) { if (sp->sa.success >= 0) sp->sa.success++; if (!rayshade(r, r->ro->omod)) /* compute contribution */ goto nomat; rayparticipate(r); return; } /* we missed our mark! */ if (sp->sa.success < 0) return; /* bitched already */ sp->sa.success -= AIMREQT; if (sp->sa.success >= 0) return; /* leniency */ sprintf(errmsg, "aiming failure for light source \"%s\"", sp->so->oname); error(WARNING, errmsg); /* issue warning */ return; nomat: objerror(r->ro, USER, "material not found"); } static int transillum( /* check if material is transparent illum */ OBJREC *m ) { m = findmaterial(m); if (m == NULL) return(1); if (m->otype != MAT_ILLUM) return(0); return(!m->oargs.nsargs || !strcmp(m->oargs.sarg[0], VOIDID)); } int sourcehit( /* check to see if ray hit distant source */ RAY *r ) { int glowsrc = -1; int transrc = -1; int first, last; int i; if (r->rsrc >= 0) { /* check only one if aimed */ first = last = r->rsrc; } else { /* otherwise check all */ first = 0; last = nsources-1; } for (i = first; i <= last; i++) { if ((source[i].sflags & (SDISTANT|SVIRTUAL)) != SDISTANT) continue; /* * Check to see if ray is within * solid angle of source. */ if (2.*PI*(1. - DOT(source[i].sloc,r->rdir)) > source[i].ss2) continue; - /* is it the only possibility? */ - if (first == last) { + /* is it what we aimed for? */ + if (i == r->rsrc) { r->ro = source[i].so; break; } /* * If it's a glow or transparent illum, just remember it. */ if (source[i].sflags & SSKIP) { if (glowsrc < 0) glowsrc = i; continue; } if (transillum(source[i].so)) { if (transrc < 0) transrc = i; continue; } r->ro = source[i].so; /* otherwise, use first hit */ break; } /* * Do we need fallback? */ if (r->ro == NULL) { if (transrc >= 0 && r->crtype & (AMBIENT|SPECULAR)) return(0); /* avoid overcounting */ if (glowsrc >= 0) r->ro = source[glowsrc].so; else return(0); /* nothing usable */ } /* * Assign object index */ r->robj = objndx(r->ro); return(1); } static int cntcmp( /* contribution compare (descending) */ const void *p1, const void *p2 ) { const CNTPTR *sc1 = (const CNTPTR *)p1; const CNTPTR *sc2 = (const CNTPTR *)p2; if (sc1->brt > sc2->brt) return(-1); if (sc1->brt < sc2->brt) return(1); return(0); } void direct( /* add direct component */ RAY *r, /* ray that hit surface */ srcdirf_t *f, /* direct component coefficient function */ void *p /* data for f */ ) { int sn; CONTRIB *scp; SRCINDEX si; int nshadcheck, ncnts; int nhits; double prob, ourthresh, hwt; RAY sr; /* PMAP: Factor in direct photons (primarily for debugging/validation) */ /* PMAP: Also primary hook in lightflow mode, since this accounts for direct _and_ indirect component */ if (directPhotonMapping #ifdef PMAP_PHOTONFLOW || lightFlowPhotonMapping #endif ) { (*f)(r -> rcol, p, r -> ron, PI); multDirectPmap(r); return; } /* NOTE: srccnt and cntord global so no recursion */ if (nsources <= 0) return; /* no sources?! */ /* potential contributions */ initsrcindex(&si); for (sn = 0; srcray(&sr, r, &si); sn++) { if (sn >= maxcntr) { maxcntr = sn + MAXSPART; srccnt = (CONTRIB *)realloc((void *)srccnt, maxcntr*sizeof(CONTRIB)); cntord = (CNTPTR *)realloc((void *)cntord, maxcntr*sizeof(CNTPTR)); if ((srccnt == NULL) | (cntord == NULL)) error(SYSTEM, "out of memory in direct"); } cntord[sn].sndx = sn; scp = srccnt + sn; scp->sno = sr.rsrc; #if SHADCACHE /* check shadow cache */ if (si.np == 1 && srcblocked(&sr)) { cntord[sn].brt = 0.0; continue; } #endif /* compute coefficient */ (*f)(scp->coef, p, sr.rdir, si.dom); cntord[sn].brt = intens(scp->coef); if (cntord[sn].brt <= 0.0) continue; VCOPY(scp->dir, sr.rdir); copycolor(sr.rcoef, scp->coef); /* compute potential */ sr.revf = srcvalue; rayvalue(&sr); multcolor(sr.rcol, sr.rcoef); copycolor(scp->val, sr.rcol); cntord[sn].brt = bright(sr.rcol); } /* sort contributions */ qsort(cntord, sn, sizeof(CNTPTR), cntcmp); { /* find last */ int l, m; ncnts = l = sn; sn = 0; while ((m = (sn + ncnts) >> 1) != l) { if (cntord[m].brt > 0.0) sn = m; else ncnts = m; l = m; } } if (ncnts == 0) return; /* no contributions! */ /* accumulate tail */ for (sn = ncnts-1; sn > 0; sn--) cntord[sn-1].brt += cntord[sn].brt; /* compute number to check */ nshadcheck = pow((double)ncnts, shadcert) + .5; /* modify threshold */ ourthresh = shadthresh / r->rweight; /* test for shadows */ for (nhits = 0, hwt = 0.0, sn = 0; sn < ncnts; hwt += (double)source[scp->sno].nhits / (double)source[scp->sno].ntests, sn++) { /* check threshold */ if (sn >= MINSHADCNT && (sn+nshadcheck>=ncnts ? cntord[sn].brt : cntord[sn].brt-cntord[sn+nshadcheck].brt) < ourthresh*bright(r->rcol)) break; scp = srccnt + cntord[sn].sndx; /* test for hit */ rayorigin(&sr, SHADOW, r, NULL); copycolor(sr.rcoef, scp->coef); VCOPY(sr.rdir, scp->dir); sr.rsrc = scp->sno; /* keep statistics */ if (source[scp->sno].ntests++ > 0xfffffff0) { source[scp->sno].ntests >>= 1; source[scp->sno].nhits >>= 1; } if (localhit(&sr, &thescene) && ( sr.ro != source[scp->sno].so || source[scp->sno].sflags & SFOLLOW )) { /* follow entire path */ raycont(&sr); if (trace != NULL) (*trace)(&sr); /* trace execution */ if (bright(sr.rcol) <= FTINY) { #if SHADCACHE if ((scp <= srccnt || scp[-1].sno != scp->sno) && (scp >= srccnt+ncnts-1 || scp[1].sno != scp->sno)) srcblocker(&sr); #endif continue; /* missed! */ } rayparticipate(&sr); multcolor(sr.rcol, sr.rcoef); copycolor(scp->val, sr.rcol); } else if (trace != NULL && (source[scp->sno].sflags & (SDISTANT|SVIRTUAL|SFOLLOW)) == (SDISTANT|SFOLLOW) && sourcehit(&sr) && rayshade(&sr, sr.ro->omod)) { (*trace)(&sr); /* trace execution */ /* skip call to rayparticipate() & scp->val update */ } /* add contribution if hit */ addcolor(r->rcol, scp->val); nhits++; source[scp->sno].nhits++; } /* source hit rate */ if (hwt > FTINY) hwt = (double)nhits / hwt; else hwt = 0.5; #ifdef DEBUG sprintf(errmsg, "%d tested, %d untested, %f conditional hit rate\n", sn, ncnts-sn, hwt); eputs(errmsg); #endif /* add in untested sources */ for ( ; sn < ncnts; sn++) { scp = srccnt + cntord[sn].sndx; prob = hwt * (double)source[scp->sno].nhits / (double)source[scp->sno].ntests; if (prob < 1.0) scalecolor(scp->val, prob); addcolor(r->rcol, scp->val); } } void srcscatter( /* compute source scattering into ray */ RAY *r ) { int oldsampndx; int nsamps; RAY sr; SRCINDEX si; double t, d; double re, ge, be; COLOR cvext; int i, j; if (r->rot >= FHUGE*.99 || r->gecc >= 1.-FTINY) return; /* this can never work */ /* PMAP: do unconditional inscattering for volume photons */ if (!volumePhotonMapping && (r->slights == NULL || r->slights[0] == 0)) return; if (ssampdist <= FTINY || (nsamps = r->rot/ssampdist + .5) < 1) nsamps = 1; #if MAXSSAMP else if (nsamps > MAXSSAMP) nsamps = MAXSSAMP; #endif oldsampndx = samplendx; samplendx = random()&0x7fff; /* randomize */ for (i = volumePhotonMapping ? 1 : r->slights[0]; i > 0; i--) { /* for each source OR once if volume photon map enabled */ for (j = 0; j < nsamps; j++) { /* for each sample position */ samplendx++; t = r->rot * (j+frandom())/nsamps; /* extinction */ re = t*colval(r->cext,RED); ge = t*colval(r->cext,GRN); be = t*colval(r->cext,BLU); setcolor(cvext, re > 92. ? 0. : exp(-re), ge > 92. ? 0. : exp(-ge), be > 92. ? 0. : exp(-be)); if (intens(cvext) <= FTINY) break; /* too far away */ sr.rorg[0] = r->rorg[0] + r->rdir[0]*t; sr.rorg[1] = r->rorg[1] + r->rdir[1]*t; sr.rorg[2] = r->rorg[2] + r->rdir[2]*t; if (!volumePhotonMapping) { initsrcindex(&si); /* sample ray to this source */ si.sn = r->slights[i]; nopart(&si, &sr); if (!srcray(&sr, NULL, &si) || sr.rsrc != r->slights[i]) continue; /* no path */ #if SHADCACHE if (srcblocked(&sr)) /* check shadow cache */ continue; #endif copycolor(sr.cext, r->cext); copycolor(sr.albedo, r->albedo); sr.gecc = r->gecc; sr.slights = r->slights; rayvalue(&sr); /* eval. source ray */ if (bright(sr.rcol) <= FTINY) { #if SHADCACHE srcblocker(&sr); /* add blocker to cache */ #endif continue; } if (r->gecc <= FTINY) /* compute P(theta) */ d = 1.; else { d = DOT(r->rdir, sr.rdir); d = 1. + r->gecc*r->gecc - 2.*r->gecc*d; d = (1. - r->gecc*r->gecc) / (d*sqrt(d)); } /* other factors */ d *= si.dom * r->rot / (4.*PI*nsamps); scalecolor(sr.rcol, d); } else { /* PMAP: Add ambient inscattering from * volume photons; note we reverse the * incident ray direction since we're * now in *backward* raytracing mode! */ sr.rdir [0] = -r -> rdir [0]; sr.rdir [1] = -r -> rdir [1]; sr.rdir [2] = -r -> rdir [2]; sr.gecc = r -> gecc; inscatterVolumePmap(&sr, sr.rcol); scalecolor(sr.rcol, r -> rot / nsamps); } multcolor(sr.rcol, r->cext); multcolor(sr.rcol, r->albedo); multcolor(sr.rcol, cvext); addcolor(r->rcol, sr.rcol); /* add it in */ } } samplendx = oldsampndx; } /**************************************************************** * The following macros were separated from the m_light() routine * because they are very nasty and difficult to understand. */ /* illumblock * * * We cannot allow an illum to pass to another illum, because that * would almost certainly constitute overcounting. * However, we do allow an illum to pass to another illum * that is actually going to relay to a virtual light source. * We also prevent an illum from passing to a glow; this provides a * convenient mechanism for defining detailed light source * geometry behind (or inside) an effective radiator. */ static int weaksrcmat(OBJREC *m) /* identify material */ { m = findmaterial(m); if (m == NULL) return(0); return((m->otype==MAT_ILLUM) | (m->otype==MAT_GLOW)); } #define illumblock(m, r) (!(source[r->rsrc].sflags&SVIRTUAL) && \ r->rod > 0.0 && \ weaksrcmat(source[r->rsrc].so)) /* wrongsource * * * This source is the wrong source (ie. overcounted) if we are * aimed to a different source than the one we hit and the one * we hit is not an illum that should be passed. */ #define wrongsource(m, r) (r->rsrc>=0 && source[r->rsrc].so!=r->ro && \ (m->otype!=MAT_ILLUM || illumblock(m,r))) /* distglow * * * A distant glow is an object that sometimes acts as a light source, * but is too far away from the test point to be one in this case. * (Glows with negative radii should NEVER participate in illumination.) */ #define distglow(m, r, d) (m->otype==MAT_GLOW && \ m->oargs.farg[3] >= -FTINY && \ d > m->oargs.farg[3]) /* badcomponent * * * We must avoid counting light sources in the ambient calculation, * since the direct component is handled separately. Therefore, any * ambient ray which hits an active light source must be discarded. * The same is true for stray specular samples, since the specular * contribution from light sources is calculated separately. */ /* PMAP: Also avoid counting sources via transferred ambient rays (e.g. * through glass) when photon mapping is enabled, as these indirect * components are already accounted for. */ #define badcomponent(m, r) (srcRayInPmap(r) || \ (r->crtype&(AMBIENT|SPECULAR) && \ !(r->crtype&SHADOW || r->rod < 0.0 || \ /* not 100% correct */ distglow(m, r, r->rot)))) /* passillum * * * An illum passes to another material type when we didn't hit it * on purpose (as part of a direct calculation), or it is relaying * a virtual light source. */ #define passillum(m, r) (m->otype==MAT_ILLUM && \ (r->rsrc<0 || source[r->rsrc].so!=r->ro || \ source[r->rsrc].sflags&SVIRTUAL)) /* srcignore * * * The -dv flag is normally on for sources to be visible. */ #define srcignore(m, r) !(directvis || r->crtype&SHADOW || \ distglow(m, r, raydist(r,PRIMARY))) int m_light( /* ray hit a light source */ OBJREC *m, RAY *r ) { /* check for over-counting */ if (badcomponent(m, r)) { setcolor(r->rcoef, 0.0, 0.0, 0.0); return(1); } if (wrongsource(m, r)) { setcolor(r->rcoef, 0.0, 0.0, 0.0); return(1); } /* check for passed illum */ if (passillum(m, r)) { if (m->oargs.nsargs && strcmp(m->oargs.sarg[0], VOIDID)) return(rayshade(r,lastmod(objndx(m),m->oargs.sarg[0]))); raytrans(r); return(1); } /* check for invisibility */ if (srcignore(m, r)) { setcolor(r->rcoef, 0.0, 0.0, 0.0); return(1); } /* otherwise treat as source */ /* check for behind */ if (r->rod < 0.0) return(1); /* check for outside spot */ if (m->otype==MAT_SPOT && spotout(r, makespot(m))) return(1); /* get distribution pattern */ raytexture(r, m->omod); /* get source color */ setcolor(r->rcol, m->oargs.farg[0], m->oargs.farg[1], m->oargs.farg[2]); /* modify value */ multcolor(r->rcol, r->pcol); return(1); } diff --git a/text.c b/text.c index 75da56d..af935b8 100644 --- a/text.c +++ b/text.c @@ -1,352 +1,354 @@ #ifndef lint -static const char RCSid[] = "$Id: text.c,v 2.27 2016/03/22 03:56:17 greg Exp $"; +static const char RCSid[] = "$Id: text.c,v 2.28 2021/11/19 22:51:31 greg Exp $"; #endif /* * text.c - functions for text patterns and mixtures. */ #include "copyright.h" #include "ray.h" #include "paths.h" #include "otypes.h" #include "rtotypes.h" #include "font.h" /* * A text pattern is specified as the text (a file or line), * the upper left anchor point, the right motion vector, the down * motion vector, and the foreground and background brightness. * For a file, the description is as follows: * * modifier brighttext id * 2 fontfile textfile * 0 * 11+ * Ax Ay Az * Rx Ry Rz * Dx Dy Dz * foreground background * [spacing] * * For a single line, we use: * * modifier brighttext id * N+2 fontfile . This is a line with N words... * 0 * 11+ * Ax Ay Az * Rx Ry Rz * Dx Dy Dz * foreground background * [spacing] * * Colortext is identical, except colors are given rather than * brightnesses. * * Mixtext has foreground and background modifiers: * * modifier mixtext id * 4+ foremod backmod fontfile text.. * 0 * 9+ * Ax Ay Az * Rx Ry Rz * Dx Dy Dz * [spacing] */ #define fndx(m) ((m)->otype==MIX_TEXT ? 2 : 0) #define tndx(m) ((m)->otype==MIX_TEXT ? 3 : 1) #define sndx(m) ((m)->otype==PAT_BTEXT ? 11 : \ (m)->otype==PAT_CTEXT ? 15 : 9) typedef struct tline { struct tline *next; /* pointer to next line */ short *spc; /* character spacing */ int width; /* total line width */ /* followed by the string */ } TLINE; #define TLSTR(l) ((char *)((l)+1)) typedef struct { FVECT right, down; /* right and down unit vectors */ FONT *f; /* our font */ TLINE tl; /* line list */ } TEXT; static TLINE * tlalloc(char *s); static TEXT * gettext(OBJREC *tm); static int intext(FVECT p, OBJREC *m); static int inglyph(double x, double y, GLYPH *gl); int do_text( OBJREC *m, RAY *r ) { FVECT v; int foreground; /* get transformed position */ if (r->rox != NULL) multp3(v, r->rop, r->rox->b.xfm); else VCOPY(v, r->rop); /* check if we are within a text glyph */ foreground = intext(v, m); /* modify */ if (m->otype == MIX_TEXT) { OBJECT omod; char *modname = m->oargs.sarg[foreground ? 0 : 1]; if (!strcmp(modname, VOIDID)) omod = OVOID; else if ((omod = lastmod(objndx(m), modname)) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", modname); objerror(m, USER, errmsg); } if (rayshade(r, omod)) { if (m->omod != OVOID) objerror(m, USER, "inappropriate modifier"); return(1); } } else if (m->otype == PAT_BTEXT) { if (foreground) scalecolor(r->pcol, m->oargs.farg[9]); else scalecolor(r->pcol, m->oargs.farg[10]); } else { /* PAT_CTEXT */ COLOR cval; if (foreground) setcolor(cval, m->oargs.farg[9], m->oargs.farg[10], m->oargs.farg[11]); else setcolor(cval, m->oargs.farg[12], m->oargs.farg[13], m->oargs.farg[14]); multcolor(r->pcol, cval); } return(0); } static TLINE * tlalloc( /* allocate and assign text line */ char *s ) { int siz; TLINE *tl; siz = strlen(s) + 1; if ((tl=(TLINE *)malloc(sizeof(TLINE)+siz)) == NULL || (tl->spc=(short *)malloc(siz*sizeof(short))) == NULL) error(SYSTEM, "out of memory in tlalloc"); tl->next = NULL; strcpy(TLSTR(tl), s); return(tl); } static TEXT * gettext( /* get text structure for material */ OBJREC *tm ) { #define R (tm->oargs.farg+3) #define D (tm->oargs.farg+6) FVECT DxR; double d; FILE *fp; char linbuf[512]; TEXT *t; int i; TLINE *tlp; char *s; if ((t = (TEXT *)tm->os) != NULL) return(t); /* check arguments */ if (tm->oargs.nsargs - tndx(tm) < 1 || tm->oargs.nfargs < sndx(tm)) objerror(tm, USER, "bad # arguments"); if ((t = (TEXT *)malloc(sizeof(TEXT))) == NULL) error(SYSTEM, "out of memory in gettext"); /* compute vectors */ fcross(DxR, D, R); fcross(t->right, DxR, D); d = DOT(t->right,t->right); if (d <= FTINY*FTINY*FTINY*FTINY) objerror(tm, USER, "illegal motion vector"); d = DOT(D,D)/d; for (i = 0; i < 3; i++) t->right[i] *= d; fcross(t->down, R, DxR); d = DOT(R,R)/DOT(t->down,t->down); for (i = 0; i < 3; i++) t->down[i] *= d; /* get text */ tlp = &t->tl; if (tm->oargs.nsargs - tndx(tm) > 1) { /* single line */ s = linbuf; for (i = tndx(tm)+1; i < tm->oargs.nsargs; i++) { strcpy(s, tm->oargs.sarg[i]); s += strlen(s); *s++ = ' '; } *--s = '\0'; tlp->next = tlalloc(linbuf); tlp = tlp->next; } else { /* text file */ if ((s = getpath(tm->oargs.sarg[tndx(tm)], getrlibpath(), R_OK)) == NULL) { sprintf(errmsg, "cannot find text file \"%s\"", tm->oargs.sarg[tndx(tm)]); error(SYSTEM, errmsg); } if ((fp = fopen(s, "r")) == NULL) { sprintf(errmsg, "cannot open text file \"%s\"", s); error(SYSTEM, errmsg); } while (fgets(linbuf, sizeof(linbuf), fp) != NULL) { s = linbuf + strlen(linbuf) - 1; if (*s == '\n') *s = '\0'; tlp->next = tlalloc(linbuf); tlp = tlp->next; } fclose(fp); } tlp->next = NULL; /* get the font */ t->f = getfont(tm->oargs.sarg[fndx(tm)]); + if (!t->f) + objerror(tm, USER, "font load error"); /* compute character spacing */ i = sndx(tm); d = i < tm->oargs.nfargs ? tm->oargs.farg[i] : 0.0; i = d * 255.0; t->tl.width = 0; for (tlp = t->tl.next; tlp != NULL; tlp = tlp->next) { if (i < 0) tlp->width = squeeztext(tlp->spc, TLSTR(tlp), t->f, -i); else if (i > 0) tlp->width = proptext(tlp->spc, TLSTR(tlp), t->f, i, 3); else tlp->width = uniftext(tlp->spc, TLSTR(tlp), t->f); if (tlp->width > t->tl.width) t->tl.width = tlp->width; } /* we're done */ tm->os = (char *)t; return(t); #undef R #undef D } void freetext( /* free text structures associated with m */ OBJREC *m ) { TEXT *tp; TLINE *tlp; tp = (TEXT *)m->os; if (tp == NULL) return; while ((tlp = tp->tl.next) != NULL) { tp->tl.next = tlp->next; free((void *)tlp->spc); free((void *)tlp); } freefont(tp->f); /* release font reference */ free((void *)tp); m->os = NULL; } static int intext( /* check to see if p is in text glyph */ FVECT p, OBJREC *m ) { TEXT *tp; TLINE *tlp; FVECT v; double y, x; int i, h; /* first, compute position in text */ tp = gettext(m); v[0] = p[0] - m->oargs.farg[0]; v[1] = p[1] - m->oargs.farg[1]; v[2] = p[2] - m->oargs.farg[2]; x = DOT(v, tp->right); i = sndx(m); if (i < m->oargs.nfargs) x *= tp->f->mwidth + 255.*fabs(m->oargs.farg[i]); else x *= 255.; h = x; i = y = DOT(v, tp->down); if (x < 0.0 || y < 0.0) return(0); x -= (double)h; y = ((i+1) - y)*255.; /* find the line position */ for (tlp = tp->tl.next; tlp != NULL; tlp = tlp->next) if (--i < 0) break; if (tlp == NULL || h >= tlp->width) return(0); for (i = 0; (h -= tlp->spc[i]) >= 0; i++) if (h < 255 && inglyph(h+x, y, tp->f->fg[TLSTR(tlp)[i]&0xff])) return(1); return(0); } static int inglyph( /* (x,y) within font glyph gl? */ double x, /* real coordinates in range [0,255) */ double y, GLYPH *gl ) { int n, ncross; int xlb, ylb; int tv; GORD *p0, *p1; if (gl == NULL) return(0); xlb = x; ylb = y; if (gl->left > xlb || gl->right <= xlb || /* check extent */ gl->bottom > ylb || gl->top <= ylb) return(0); xlb = xlb<<1 | 1; /* add 1/2 to test points... */ ylb = ylb<<1 | 1; /* ...so no equal comparisons */ n = gl->nverts; /* get # of vertices */ p0 = gvlist(gl) + 2*(n-1); /* connect last to first */ p1 = gvlist(gl); ncross = 0; /* positive x axis cross test */ while (n--) { if ((p0[1]<<1 > ylb) ^ (p1[1]<<1 > ylb)) { tv = (p0[0]<<1 > xlb) | ((p1[0]<<1 > xlb) << 1); if (tv == 03) ncross++; else if (tv) ncross += (p1[1] > p0[1]) ^ ((p0[1]-y)*(p1[0]-x) > (p0[0]-x)*(p1[1]-y)); } p0 = p1; p1 += 2; } return(ncross & 01); }