PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/core/mipmap.h
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