PBRT
|
00001 00002 /* 00003 pbrt source code Copyright(c) 1998-2012 Matt Pharr and Greg Humphreys. 00004 00005 This file is part of pbrt. 00006 00007 Redistribution and use in source and binary forms, with or without 00008 modification, are permitted provided that the following conditions are 00009 met: 00010 00011 - Redistributions of source code must retain the above copyright 00012 notice, this list of conditions and the following disclaimer. 00013 00014 - Redistributions in binary form must reproduce the above copyright 00015 notice, this list of conditions and the following disclaimer in the 00016 documentation and/or other materials provided with the distribution. 00017 00018 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 00019 IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 00020 TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 00021 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 00022 HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00023 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 00024 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00025 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00026 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00027 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 00028 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00029 00030 */ 00031 00032 #if defined(_MSC_VER) 00033 #pragma once 00034 #endif 00035 00036 #ifndef PBRT_CORE_MIPMAP_H 00037 #define PBRT_CORE_MIPMAP_H 00038 00039 // core/mipmap.h* 00040 #include "pbrt.h" 00041 #include "spectrum.h" 00042 #include "texture.h" 00043 00044 // MIPMap Declarations 00045 typedef enum { 00046 TEXTURE_REPEAT, 00047 TEXTURE_BLACK, 00048 TEXTURE_CLAMP 00049 } ImageWrap; 00050 template <typename T> class MIPMap { 00051 public: 00052 // MIPMap Public Methods 00053 MIPMap() { pyramid = NULL; width = height = nLevels = 0; } 00054 MIPMap(uint32_t xres, uint32_t yres, const T *data, bool doTri = false, 00055 float maxAniso = 8.f, ImageWrap wrapMode = TEXTURE_REPEAT); 00056 ~MIPMap(); 00057 uint32_t Width() const { return width; } 00058 uint32_t Height() const { return height; } 00059 uint32_t Levels() const { return nLevels; } 00060 const T &Texel(uint32_t level, int s, int t) const; 00061 T Lookup(float s, float t, float width = 0.f) const; 00062 T Lookup(float s, float t, float ds0, float dt0, 00063 float ds1, float dt1) const; 00064 private: 00065 // MIPMap Private Methods 00066 struct ResampleWeight; 00067 ResampleWeight *resampleWeights(uint32_t oldres, uint32_t newres) { 00068 Assert(newres >= oldres); 00069 ResampleWeight *wt = new ResampleWeight[newres]; 00070 float filterwidth = 2.f; 00071 for (uint32_t i = 0; i < newres; ++i) { 00072 // Compute image resampling weights for _i_th texel 00073 float center = (i + .5f) * oldres / newres; 00074 wt[i].firstTexel = Floor2Int((center - filterwidth) + 0.5f); 00075 for (int j = 0; j < 4; ++j) { 00076 float pos = wt[i].firstTexel + j + .5f; 00077 wt[i].weight[j] = Lanczos((pos - center) / filterwidth); 00078 } 00079 00080 // Normalize filter weights for texel resampling 00081 float invSumWts = 1.f / (wt[i].weight[0] + wt[i].weight[1] + 00082 wt[i].weight[2] + wt[i].weight[3]); 00083 for (uint32_t j = 0; j < 4; ++j) 00084 wt[i].weight[j] *= invSumWts; 00085 } 00086 return wt; 00087 } 00088 float clamp(float v) { return Clamp(v, 0.f, INFINITY); } 00089 RGBSpectrum clamp(const RGBSpectrum &v) { return v.Clamp(0.f, INFINITY); } 00090 SampledSpectrum clamp(const SampledSpectrum &v) { return v.Clamp(0.f, INFINITY); } 00091 T triangle(uint32_t level, float s, float t) const; 00092 T EWA(uint32_t level, float s, float t, float ds0, float dt0, float ds1, float dt1) const; 00093 00094 // MIPMap Private Data 00095 bool doTrilinear; 00096 float maxAnisotropy; 00097 ImageWrap wrapMode; 00098 struct ResampleWeight { 00099 int firstTexel; 00100 float weight[4]; 00101 }; 00102 BlockedArray<T> **pyramid; 00103 uint32_t width, height, nLevels; 00104 #define WEIGHT_LUT_SIZE 128 00105 static float *weightLut; 00106 }; 00107 00108 00109 00110 // MIPMap Method Definitions 00111 template <typename T> 00112 MIPMap<T>::MIPMap(uint32_t sres, uint32_t tres, const T *img, bool doTri, 00113 float maxAniso, ImageWrap wm) { 00114 doTrilinear = doTri; 00115 maxAnisotropy = maxAniso; 00116 wrapMode = wm; 00117 T *resampledImage = NULL; 00118 if (!IsPowerOf2(sres) || !IsPowerOf2(tres)) { 00119 // Resample image to power-of-two resolution 00120 uint32_t sPow2 = RoundUpPow2(sres), tPow2 = RoundUpPow2(tres); 00121 00122 // Resample image in $s$ direction 00123 ResampleWeight *sWeights = resampleWeights(sres, sPow2); 00124 resampledImage = new T[sPow2 * tPow2]; 00125 00126 // Apply _sWeights_ to zoom in $s$ direction 00127 for (uint32_t t = 0; t < tres; ++t) { 00128 for (uint32_t s = 0; s < sPow2; ++s) { 00129 // Compute texel $(s,t)$ in $s$-zoomed image 00130 resampledImage[t*sPow2+s] = 0.; 00131 for (int j = 0; j < 4; ++j) { 00132 int origS = sWeights[s].firstTexel + j; 00133 if (wrapMode == TEXTURE_REPEAT) 00134 origS = Mod(origS, sres); 00135 else if (wrapMode == TEXTURE_CLAMP) 00136 origS = Clamp(origS, 0, sres-1); 00137 if (origS >= 0 && origS < (int)sres) 00138 resampledImage[t*sPow2+s] += sWeights[s].weight[j] * 00139 img[t*sres + origS]; 00140 } 00141 } 00142 } 00143 delete[] sWeights; 00144 00145 // Resample image in $t$ direction 00146 ResampleWeight *tWeights = resampleWeights(tres, tPow2); 00147 T *workData = new T[tPow2]; 00148 for (uint32_t s = 0; s < sPow2; ++s) { 00149 for (uint32_t t = 0; t < tPow2; ++t) { 00150 workData[t] = 0.; 00151 for (uint32_t j = 0; j < 4; ++j) { 00152 int offset = tWeights[t].firstTexel + j; 00153 if (wrapMode == TEXTURE_REPEAT) offset = Mod(offset, tres); 00154 else if (wrapMode == TEXTURE_CLAMP) offset = Clamp(offset, 0, tres-1); 00155 if (offset >= 0 && offset < (int)tres) 00156 workData[t] += tWeights[t].weight[j] * 00157 resampledImage[offset*sPow2 + s]; 00158 } 00159 } 00160 for (uint32_t t = 0; t < tPow2; ++t) 00161 resampledImage[t*sPow2 + s] = clamp(workData[t]); 00162 } 00163 delete[] workData; 00164 delete[] tWeights; 00165 img = resampledImage; 00166 sres = sPow2; 00167 tres = tPow2; 00168 } 00169 width = sres; 00170 height = tres; 00171 // Initialize levels of MIPMap from image 00172 nLevels = 1 + Log2Int(float(max(sres, tres))); 00173 pyramid = new BlockedArray<T> *[nLevels]; 00174 00175 // Initialize most detailed level of MIPMap 00176 pyramid[0] = new BlockedArray<T>(sres, tres, img); 00177 for (uint32_t i = 1; i < nLevels; ++i) { 00178 // Initialize $i$th MIPMap level from $i-1$st level 00179 uint32_t sRes = max(1u, pyramid[i-1]->uSize()/2); 00180 uint32_t tRes = max(1u, pyramid[i-1]->vSize()/2); 00181 pyramid[i] = new BlockedArray<T>(sRes, tRes); 00182 00183 // Filter four texels from finer level of pyramid 00184 for (uint32_t t = 0; t < tRes; ++t) 00185 for (uint32_t s = 0; s < sRes; ++s) 00186 (*pyramid[i])(s, t) = .25f * 00187 (Texel(i-1, 2*s, 2*t) + Texel(i-1, 2*s+1, 2*t) + 00188 Texel(i-1, 2*s, 2*t+1) + Texel(i-1, 2*s+1, 2*t+1)); 00189 } 00190 if (resampledImage) delete[] resampledImage; 00191 // Initialize EWA filter weights if needed 00192 if (!weightLut) { 00193 weightLut = AllocAligned<float>(WEIGHT_LUT_SIZE); 00194 for (int i = 0; i < WEIGHT_LUT_SIZE; ++i) { 00195 float alpha = 2; 00196 float r2 = float(i) / float(WEIGHT_LUT_SIZE - 1); 00197 weightLut[i] = expf(-alpha * r2) - expf(-alpha); 00198 } 00199 } 00200 } 00201 00202 00203 template <typename T> 00204 const T &MIPMap<T>::Texel(uint32_t level, int s, int t) const { 00205 Assert(level < nLevels); 00206 const BlockedArray<T> &l = *pyramid[level]; 00207 // Compute texel $(s,t)$ accounting for boundary conditions 00208 switch (wrapMode) { 00209 case TEXTURE_REPEAT: 00210 s = Mod(s, l.uSize()); 00211 t = Mod(t, l.vSize()); 00212 break; 00213 case TEXTURE_CLAMP: 00214 s = Clamp(s, 0, l.uSize() - 1); 00215 t = Clamp(t, 0, l.vSize() - 1); 00216 break; 00217 case TEXTURE_BLACK: { 00218 static const T black = 0.f; 00219 if (s < 0 || s >= (int)l.uSize() || 00220 t < 0 || t >= (int)l.vSize()) 00221 return black; 00222 break; 00223 } 00224 } 00225 PBRT_ACCESSED_TEXEL(const_cast<MIPMap<T> *>(this), level, s, t); 00226 return l(s, t); 00227 } 00228 00229 00230 template <typename T> 00231 MIPMap<T>::~MIPMap() { 00232 for (uint32_t i = 0; i < nLevels; ++i) 00233 delete pyramid[i]; 00234 delete[] pyramid; 00235 } 00236 00237 00238 template <typename T> 00239 T MIPMap<T>::Lookup(float s, float t, float width) const { 00240 // Compute MIPMap level for trilinear filtering 00241 float level = nLevels - 1 + Log2(max(width, 1e-8f)); 00242 00243 // Perform trilinear interpolation at appropriate MIPMap level 00244 PBRT_MIPMAP_TRILINEAR_FILTER(const_cast<MIPMap<T> *>(this), s, t, width, level, nLevels); 00245 if (level < 0) 00246 return triangle(0, s, t); 00247 else if (level >= nLevels - 1) 00248 return Texel(nLevels-1, 0, 0); 00249 else { 00250 uint32_t iLevel = Floor2Int(level); 00251 float delta = level - iLevel; 00252 return (1.f-delta) * triangle(iLevel, s, t) + 00253 delta * triangle(iLevel+1, s, t); 00254 } 00255 } 00256 00257 00258 template <typename T> 00259 T MIPMap<T>::triangle(uint32_t level, float s, float t) const { 00260 level = Clamp(level, 0, nLevels-1); 00261 s = s * pyramid[level]->uSize() - 0.5f; 00262 t = t * pyramid[level]->vSize() - 0.5f; 00263 int s0 = Floor2Int(s), t0 = Floor2Int(t); 00264 float ds = s - s0, dt = t - t0; 00265 return (1.f-ds) * (1.f-dt) * Texel(level, s0, t0) + 00266 (1.f-ds) * dt * Texel(level, s0, t0+1) + 00267 ds * (1.f-dt) * Texel(level, s0+1, t0) + 00268 ds * dt * Texel(level, s0+1, t0+1); 00269 } 00270 00271 00272 template <typename T> 00273 T MIPMap<T>::Lookup(float s, float t, float ds0, float dt0, 00274 float ds1, float dt1) const { 00275 if (doTrilinear) { 00276 PBRT_STARTED_TRILINEAR_TEXTURE_LOOKUP(s, t); 00277 T val = Lookup(s, t, 00278 2.f * max(max(fabsf(ds0), fabsf(dt0)), 00279 max(fabsf(ds1), fabsf(dt1)))); 00280 PBRT_FINISHED_TRILINEAR_TEXTURE_LOOKUP(); 00281 return val; 00282 } 00283 PBRT_STARTED_EWA_TEXTURE_LOOKUP(s, t); 00284 // Compute ellipse minor and major axes 00285 if (ds0*ds0 + dt0*dt0 < ds1*ds1 + dt1*dt1) { 00286 swap(ds0, ds1); 00287 swap(dt0, dt1); 00288 } 00289 float majorLength = sqrtf(ds0*ds0 + dt0*dt0); 00290 float minorLength = sqrtf(ds1*ds1 + dt1*dt1); 00291 00292 // Clamp ellipse eccentricity if too large 00293 if (minorLength * maxAnisotropy < majorLength && minorLength > 0.f) { 00294 float scale = majorLength / (minorLength * maxAnisotropy); 00295 ds1 *= scale; 00296 dt1 *= scale; 00297 minorLength *= scale; 00298 } 00299 if (minorLength == 0.f) { 00300 PBRT_FINISHED_EWA_TEXTURE_LOOKUP(); 00301 PBRT_STARTED_TRILINEAR_TEXTURE_LOOKUP(s, t); 00302 T val = triangle(0, s, t); 00303 PBRT_FINISHED_TRILINEAR_TEXTURE_LOOKUP(); 00304 return val; 00305 } 00306 00307 // Choose level of detail for EWA lookup and perform EWA filtering 00308 float lod = max(0.f, nLevels - 1.f + Log2(minorLength)); 00309 uint32_t ilod = Floor2Int(lod); 00310 PBRT_MIPMAP_EWA_FILTER(const_cast<MIPMap<T> *>(this), s, t, ds0, ds1, dt0, dt1, minorLength, majorLength, lod, nLevels); 00311 float d = lod - ilod; 00312 T val = (1.f - d) * EWA(ilod, s, t, ds0, dt0, ds1, dt1) + 00313 d * EWA(ilod+1, s, t, ds0, dt0, ds1, dt1); 00314 PBRT_FINISHED_EWA_TEXTURE_LOOKUP(); 00315 return val; 00316 } 00317 00318 00319 template <typename T> 00320 T MIPMap<T>::EWA(uint32_t level, float s, float t, float ds0, float dt0, 00321 float ds1, float dt1) const { 00322 if (level >= nLevels) return Texel(nLevels-1, 0, 0); 00323 // Convert EWA coordinates to appropriate scale for level 00324 s = s * pyramid[level]->uSize() - 0.5f; 00325 t = t * pyramid[level]->vSize() - 0.5f; 00326 ds0 *= pyramid[level]->uSize(); 00327 dt0 *= pyramid[level]->vSize(); 00328 ds1 *= pyramid[level]->uSize(); 00329 dt1 *= pyramid[level]->vSize(); 00330 00331 // Compute ellipse coefficients to bound EWA filter region 00332 float A = dt0*dt0 + dt1*dt1 + 1; 00333 float B = -2.f * (ds0*dt0 + ds1*dt1); 00334 float C = ds0*ds0 + ds1*ds1 + 1; 00335 float invF = 1.f / (A*C - B*B*0.25f); 00336 A *= invF; 00337 B *= invF; 00338 C *= invF; 00339 00340 // Compute the ellipse's $(s,t)$ bounding box in texture space 00341 float det = -B*B + 4.f*A*C; 00342 float invDet = 1.f / det; 00343 float uSqrt = sqrtf(det * C), vSqrt = sqrtf(A * det); 00344 int s0 = Ceil2Int (s - 2.f * invDet * uSqrt); 00345 int s1 = Floor2Int(s + 2.f * invDet * uSqrt); 00346 int t0 = Ceil2Int (t - 2.f * invDet * vSqrt); 00347 int t1 = Floor2Int(t + 2.f * invDet * vSqrt); 00348 00349 // Scan over ellipse bound and compute quadratic equation 00350 T sum(0.); 00351 float sumWts = 0.f; 00352 for (int it = t0; it <= t1; ++it) { 00353 float tt = it - t; 00354 for (int is = s0; is <= s1; ++is) { 00355 float ss = is - s; 00356 // Compute squared radius and filter texel if inside ellipse 00357 float r2 = A*ss*ss + B*ss*tt + C*tt*tt; 00358 if (r2 < 1.) { 00359 float weight = weightLut[min(Float2Int(r2 * WEIGHT_LUT_SIZE), 00360 WEIGHT_LUT_SIZE-1)]; 00361 sum += Texel(level, is, it) * weight; 00362 sumWts += weight; 00363 } 00364 } 00365 } 00366 return sum / sumWts; 00367 } 00368 00369 00370 template <typename T> float *MIPMap<T>::weightLut = NULL; 00371 00372 #endif // PBRT_CORE_MIPMAP_H