00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "volume.h"
00026
00027 COREDLL
00028 float PhaseIsotropic(const Vector &, const Vector &) {
00029 return 1.f / (4.f * M_PI);
00030 }
00031 COREDLL float PhaseRayleigh(const Vector &w, const Vector &wp) {
00032 float costheta = Dot(w, wp);
00033 return 3.f/(16.f*M_PI) * (1 + costheta * costheta);
00034 }
00035 COREDLL float PhaseMieHazy(const Vector &w, const Vector &wp) {
00036 float costheta = Dot(w, wp);
00037 return (0.5f + 4.5f * powf(0.5 * (1.f + costheta), 8.f)) / (4.f*M_PI);
00038 }
00039 COREDLL float PhaseMieMurky(const Vector &w, const Vector &wp) {
00040 float costheta = Dot(w, wp);
00041 return (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f)) / (4.f*M_PI);
00042 }
00043 COREDLL
00044 float PhaseHG(const Vector &w, const Vector &wp, float g) {
00045 float costheta = Dot(w, wp);
00046 return 1.f / (4.f * M_PI) * (1.f - g*g) /
00047 powf(1.f + g*g - 2.f * g * costheta, 1.5f);
00048 }
00049 COREDLL
00050 float PhaseSchlick(const Vector &w,
00051 const Vector &wp, float g) {
00052 float k = 1.55f * g - .55f * g * g * g;
00053 float kcostheta = k * Dot(w, wp);
00054 return 1.f / (4.f * M_PI) * (1.f - k*k) /
00055 ((1.f - kcostheta) * (1.f - kcostheta));
00056 }
00057 Spectrum VolumeRegion::sigma_t(const Point &p,
00058 const Vector &w) const {
00059 return sigma_a(p, w) + sigma_s(p, w);
00060 }
00061 DensityRegion::DensityRegion(const Spectrum &sa,
00062 const Spectrum &ss, float gg,
00063 const Spectrum &emit,
00064 const Transform &VolumeToWorld)
00065 : sig_a(sa), sig_s(ss), le(emit), g(gg) {
00066 WorldToVolume = VolumeToWorld.GetInverse();
00067 }
00068 AggregateVolume::
00069 AggregateVolume(const vector<VolumeRegion *> &r) {
00070 regions = r;
00071 for (u_int i = 0; i < regions.size(); ++i)
00072 bound = Union(bound, regions[i]->WorldBound());
00073 }
00074 Spectrum AggregateVolume::sigma_a(const Point &p,
00075 const Vector &w) const {
00076 Spectrum s(0.);
00077 for (u_int i = 0; i < regions.size(); ++i)
00078 s += regions[i]->sigma_a(p, w);
00079 return s;
00080 }
00081 Spectrum AggregateVolume::sigma_s(const Point &p, const Vector &w) const {
00082 Spectrum s(0.);
00083 for (u_int i = 0; i < regions.size(); ++i)
00084 s += regions[i]->sigma_s(p, w);
00085 return s;
00086 }
00087 Spectrum AggregateVolume::Lve(const Point &p, const Vector &w) const {
00088 Spectrum L(0.);
00089 for (u_int i = 0; i < regions.size(); ++i)
00090 L += regions[i]->Lve(p, w);
00091 return L;
00092 }
00093 float AggregateVolume::p(const Point &p, const Vector &w, const Vector &wp) const {
00094 float ph = 0, sumWt = 0;
00095 for (u_int i = 0; i < regions.size(); ++i) {
00096 float sigt = regions[i]->sigma_t(p, w).y();
00097 if (sigt != 0.f) {
00098 float wt = regions[i]->sigma_s(p, w).y() / sigt;
00099 sumWt += wt;
00100 ph += wt * regions[i]->p(p, w, wp);
00101 }
00102 }
00103 return ph / sumWt;
00104 }
00105 Spectrum AggregateVolume::sigma_t(const Point &p, const Vector &w) const {
00106 Spectrum s(0.);
00107 for (u_int i = 0; i < regions.size(); ++i)
00108 s += regions[i]->sigma_t(p, w);
00109 return s;
00110 }
00111 Spectrum AggregateVolume::Tau(const Ray &ray, float step, float offset) const {
00112 Spectrum t(0.);
00113 for (u_int i = 0; i < regions.size(); ++i)
00114 t += regions[i]->Tau(ray, step, offset);
00115 return t;
00116 }
00117 bool AggregateVolume::IntersectP(const Ray &ray,
00118 float *t0, float *t1) const {
00119 *t0 = INFINITY;
00120 *t1 = -INFINITY;
00121 for (u_int i = 0; i < regions.size(); ++i) {
00122 float tr0, tr1;
00123 if (regions[i]->IntersectP(ray, &tr0, &tr1)) {
00124 *t0 = min(*t0, tr0);
00125 *t1 = max(*t1, tr1);
00126 }
00127 }
00128 return (*t0 < *t1);
00129 }
00130 AggregateVolume::~AggregateVolume() {
00131 for (u_int i = 0; i < regions.size(); ++i)
00132 delete regions[i];
00133 }
00134 BBox AggregateVolume::WorldBound() const {
00135 return bound;
00136 }
00137 Spectrum DensityRegion::Tau(const Ray &r,
00138 float stepSize, float u) const {
00139 float t0, t1;
00140 float length = r.d.Length();
00141 if (length == 0.f) return 0.f;
00142 Ray rn(r.o, r.d / length,
00143 r.mint * length,
00144 r.maxt * length);
00145 if (!IntersectP(rn, &t0, &t1)) return 0.;
00146 Spectrum tau(0.);
00147 t0 += u * stepSize;
00148 while (t0 < t1) {
00149 tau += sigma_t(rn(t0), -rn.d);
00150 t0 += stepSize;
00151 }
00152 return tau * stepSize;
00153 }