template <typename T>
T MIPMap<T>::Lookup(float s, float t, float ds0, float dt0,
                    float ds1, float dt1) const {
    if (doTrilinear) {
        PBRT_STARTED_TRILINEAR_TEXTURE_LOOKUP(s, t);
        T val = Lookup(s, t,
                       2.f * max(max(fabsf(ds0), fabsf(dt0)),
                                 max(fabsf(ds1), fabsf(dt1))));
        PBRT_FINISHED_TRILINEAR_TEXTURE_LOOKUP();
        return val;
    }
    PBRT_STARTED_EWA_TEXTURE_LOOKUP(s, t);
    // Compute ellipse minor and major axes
    if (ds0*ds0 + dt0*dt0 < ds1*ds1 + dt1*dt1) {
        swap(ds0, ds1);
        swap(dt0, dt1);
    }
    float majorLength = sqrtf(ds0*ds0 + dt0*dt0);
    float minorLength = sqrtf(ds1*ds1 + dt1*dt1);

    // Clamp ellipse eccentricity if too large
    if (minorLength * maxAnisotropy < majorLength && minorLength > 0.f) {
        float scale = majorLength / (minorLength * maxAnisotropy);
        ds1 *= scale;
        dt1 *= scale;
        minorLength *= scale;
    }
    if (minorLength == 0.f) {
        PBRT_FINISHED_EWA_TEXTURE_LOOKUP();
        PBRT_STARTED_TRILINEAR_TEXTURE_LOOKUP(s, t);
        T val = triangle(0, s, t);
        PBRT_FINISHED_TRILINEAR_TEXTURE_LOOKUP();
        return val;
    }

    // Choose level of detail for EWA lookup and perform EWA filtering
    float lod = max(0.f, nLevels - 1.f + Log2(minorLength));
    uint32_t ilod = Floor2Int(lod);
    PBRT_MIPMAP_EWA_FILTER(const_cast<MIPMap<T> *>(this), s, t, ds0, ds1, dt0, dt1, minorLength, majorLength, lod, nLevels);
    float d = lod - ilod;
    T val = (1.f - d) * EWA(ilod,   s, t, ds0, dt0, ds1, dt1) +
                   d  * EWA(ilod+1, s, t, ds0, dt0, ds1, dt1);
    PBRT_FINISHED_EWA_TEXTURE_LOOKUP();
    return val;
}
 

作用：

（EWA 的思想就是，一個點 的 橢圓區域 內的所有的點，都進行高斯過濾。）

The elliptically weighted average (EWA) algorithm fits an ellipse(橢圓) to the two axes in texture
space given by the texture coordinate differentials and then filters the texture with a
Gaussian filter function (Figure 10.15). It is widely regarded as one of the best texture
filtering algorithms in graphics and has been carefully derived from the basic principles
of sampling theory.

Unlike the triangle filter in the previous section, it can filter over
arbitrarily oriented regions of the texture, with some flexibility of having different filter
extents in different directions. This type of filter is known as anisotropic. This capability
greatly improves the quality of its results, since it can properly adapt to different sampling
rates along the two image axes.

Figure 10.15: The EWA filter applies a Gaussian filter to the texels in an elliptical area around the
evaluation point. The extent of the ellipse is such that its edge passes through the positions of the
adjacent texture samples as estimated by the texture coordinate partial derivatives.

 

細節

a.

    if (doTrilinear) {
        PBRT_STARTED_TRILINEAR_TEXTURE_LOOKUP(s, t);
        T val = Lookup(s, t,
                       2.f * max(max(fabsf(ds0), fabsf(dt0)),
                                 max(fabsf(ds1), fabsf(dt1))));
        PBRT_FINISHED_TRILINEAR_TEXTURE_LOOKUP();
        return val;
    }

作用：

先判斷是否設定了 doTrilinear，如果設定了的話，就直接使用 上一節的 《PBRT_V2 總結記錄 <51> MIPMap.Lookup (Triangle Filter)》，值得注意的話，上一節講到  利用 texel 的寬度來確定 MIPMap 的層數，那麼這裡就可以知道，這個 texel 的寬度 其實就是利用 引數 ds0，dt0，ds1,dt1 來確定的，這些引數的意義可以參考《PBRT_V2 總結記錄 <46> TextureMapping2D 和 UVMapping2D》，就是 dsdx，dtdx，dsdy, dtdy 偏導數，當 sample的位置，x,y 改變的時候，紋理座標u，v 怎麼改變。也就是 texel的 寬度用 紋理座標對螢幕空間pixel sample的偏導數決定。

 

b. 

    // Compute ellipse minor and major axes
    if (ds0*ds0 + dt0*dt0 < ds1*ds1 + dt1*dt1) {
        swap(ds0, ds1);
        swap(dt0, dt1);
    }
    float majorLength = sqrtf(ds0*ds0 + dt0*dt0);
    float minorLength = sqrtf(ds1*ds1 + dt1*dt1);

作用：

（其實上面所說的，dsdx，dtdx，dsdy, dtdy 偏導數 來確定 哪一個是 elipse的 major 軸，哪一個是 minor 軸）

The screen space partial derivatives(偏導數) of the texture coordinates define the axes of the
ellipse. The lookup method starts out by determining which of the two axes is the major
axis (the longer of the two) and which is the minor, swapping them if needed so that
(ds0,dt0) is the major axis. The length of the minor axis will be used shortly to select a
MIP map level.

 

c.

    // Clamp ellipse eccentricity if too large
    if (minorLength * maxAnisotropy < majorLength && minorLength > 0.f) {
        float scale = majorLength / (minorLength * maxAnisotropy);
        ds1 *= scale;
        dt1 *= scale;
        minorLength *= scale;
    }
    if (minorLength == 0.f) {
        PBRT_FINISHED_EWA_TEXTURE_LOOKUP();
        PBRT_STARTED_TRILINEAR_TEXTURE_LOOKUP(s, t);
        T val = triangle(0, s, t);
        PBRT_FINISHED_TRILINEAR_TEXTURE_LOOKUP();
        return val;
    }

作用：

（離心率(major / minor ) 高的話，會表明一個又長又瘦的橢圓，而且意味著 這個橢圓會包含很多 texel，所以，通過 增加 minor 軸的長度來降低離心率,減少橢圓包含texel的數目）

Next the eccentricity(離心率) of the ellipse(橢圓) is computed—the ratio of the length of the major axis
to the length of the minor axis. A large eccentricity(離心率) indicates(表明) a very long and skinny
ellipse. Because this method filters texels from a MIP map level chosen based on the
length of the minor axis, highly eccentric ellipses mean that a large number of texels need
to be filtered. To avoid this expense (and to ensure that any EWA lookup takes a bounded
amount of time), the length of the minor axis may be increased to limit the eccentricity.
The result may be an increase in blurring, although this effect usually isn’t noticeable in
practice.

 

d.

    // Choose level of detail for EWA lookup and perform EWA filtering
    float lod = max(0.f, nLevels - 1.f + Log2(minorLength));
    uint32_t ilod = Floor2Int(lod);
    PBRT_MIPMAP_EWA_FILTER(const_cast<MIPMap<T> *>(this), s, t, ds0, ds1, dt0, dt1, minorLength, majorLength, lod, nLevels);
    float d = lod - ilod;
    T val = (1.f - d) * EWA(ilod,   s, t, ds0, dt0, ds1, dt1) +
                   d  * EWA(ilod+1, s, t, ds0, dt0, ds1, dt1);
    PBRT_FINISHED_EWA_TEXTURE_LOOKUP();
    return val;

作用：

（EWA 利用 minor 軸的長度 作為 texel 的 寬度，來確定 在哪一層進行 過濾，計算方法與 上一節的 《PBRT_V2 總結記錄 <51> MIPMap.Lookup (Triangle Filter)》基本一致，找到 level 層之後，思路其實也與 上一節的 基本一致，只不過是 上一節的 triangle 換了 EWA ，也就是，不使用 雙線性插值 混合 第 level 層 （s,t）座標附近4個texel 的值，而是使用 EWA的方法）

Like the triangle filter, the EWA filter uses the image pyramid to reduce the number of
texels to be filtered for a particular texture lookup, choosing a MIP map level based on the
length of the minor axis. Given the limited eccentricity(離心率) of the ellipse due to the clamping
above, the total number of texels used is thus bounded. Given the length of the minor
axis, the computation to find the appropriate pyramid level is the same as was used for the
triangle filter. Similarly, the implementation here blends between the filtered results at the
two levels around the computed level of detail, again to reduce artifacts from transitions(過度)
from one level to another.

 

f.

MIPMap::EWA()

template <typename T>
T MIPMap<T>::EWA(uint32_t level, float s, float t, float ds0, float dt0,
                 float ds1, float dt1) const {
    if (level >= nLevels) return Texel(nLevels-1, 0, 0);
    // Convert EWA coordinates to appropriate scale for level
    s = s * pyramid[level]->uSize() - 0.5f;
    t = t * pyramid[level]->vSize() - 0.5f;
    ds0 *= pyramid[level]->uSize();
    dt0 *= pyramid[level]->vSize();
    ds1 *= pyramid[level]->uSize();
    dt1 *= pyramid[level]->vSize();

    // Compute ellipse coefficients to bound EWA filter region
    float A = dt0*dt0 + dt1*dt1 + 1;
    float B = -2.f * (ds0*dt0 + ds1*dt1);
    float C = ds0*ds0 + ds1*ds1 + 1;
    float invF = 1.f / (A*C - B*B*0.25f);
    A *= invF;
    B *= invF;
    C *= invF;

    // Compute the ellipse's $(s,t)$ bounding box in texture space
    float det = -B*B + 4.f*A*C;
    float invDet = 1.f / det;
    float uSqrt = sqrtf(det * C), vSqrt = sqrtf(A * det);
    int s0 = Ceil2Int (s - 2.f * invDet * uSqrt);
    int s1 = Floor2Int(s + 2.f * invDet * uSqrt);
    int t0 = Ceil2Int (t - 2.f * invDet * vSqrt);
    int t1 = Floor2Int(t + 2.f * invDet * vSqrt);

    // Scan over ellipse bound and compute quadratic equation
    T sum(0.);
    float sumWts = 0.f;
    for (int it = t0; it <= t1; ++it) {
        float tt = it - t;
        for (int is = s0; is <= s1; ++is) {
            float ss = is - s;
            // Compute squared radius and filter texel if inside ellipse
            float r2 = A*ss*ss + B*ss*tt + C*tt*tt;
            if (r2 < 1.) {
                float weight = weightLut[min(Float2Int(r2 * WEIGHT_LUT_SIZE),
                                             WEIGHT_LUT_SIZE-1)];
                sum += Texel(level, is, it) * weight;
                sumWts += weight;
            }
        }
    }
    return sum / sumWts;
}

作用：

（從一個橢圓中計算一個包圍盒，判斷包圍盒裡面的texel是否在橢圓中，利用 高斯分佈，得到每一個texel 在 橢圓中的 權值，再進行 加權平均）

The MIPMap::EWA() method actually applies the filter at a particular level.

 

細節

a.

    s = s * pyramid[level]->uSize() - 0.5f;
    t = t * pyramid[level]->vSize() - 0.5f;
    ds0 *= pyramid[level]->uSize();
    dt0 *= pyramid[level]->vSize();
    ds1 *= pyramid[level]->uSize();
    dt1 *= pyramid[level]->vSize();

作用：

（s，t,和 解析度範圍，從 [0,1] -> [0, imageResolution],）

This method first converts from texture coordinates in [0, 1] to coordinates and differentials
in terms of the resolution of the chosen MIP map level. It also subtracts 0.5 from
the continuous position coordinate to align the sample point with the discrete texel coordinates,
as was done in MIPMap::triangle().

 

b.

    // Compute ellipse coefficients to bound EWA filter region
    float A = dt0*dt0 + dt1*dt1 + 1;
    float B = -2.f * (ds0*dt0 + ds1*dt1);
    float C = ds0*ds0 + ds1*ds1 + 1;
    float invF = 1.f / (A*C - B*B*0.25f);
    A *= invF;
    B *= invF;
    C *= invF;

作用：

（計算 橢圓方程的 A， B， C，1/F 引數，這個橢圓的中心預設在原點，兩個軸為(ds0, dt0), (ds1, dt1)， 看上面的程式碼可以發現 (ds0,dt0) 是 major 軸，(ds1,dt1) 是 minor 軸， 滿足下面方程的 點都是在橢圓裡面）

It next computes the coefficients of the implicit equation for the ellipse(橢圓) with axes
(ds0,dt0) and (ds1,dt1) and centered at the origin. Placing the ellipse at the origin

rather than at (s , t) simplifies the implicit equation and the computation of its coefficients,
and can be easily corrected for when the equation is evaluated later. The general
form of the implicit equation for all points (s , t) inside such an ellipse is

although it is more computationally efficient to divide through by F and express this as

 

c.

    // Compute the ellipse's $(s,t)$ bounding box in texture space
    float det = -B*B + 4.f*A*C;
    float invDet = 1.f / det;
    float uSqrt = sqrtf(det * C), vSqrt = sqrtf(A * det);
    int s0 = Ceil2Int (s - 2.f * invDet * uSqrt);
    int s1 = Floor2Int(s + 2.f * invDet * uSqrt);
    int t0 = Ceil2Int (t - 2.f * invDet * vSqrt);
    int t1 = Floor2Int(t + 2.f * invDet * vSqrt);

作用：

（這裡就是計算 橢圓 裡面的 包圍盒，min，max 的兩個點 (s0, t0) (s1, t1)）

The next step is to find the axis-aligned bounding box in discrete integer texel coordinates
of the texels that are potentially inside the ellipse. The EWA algorithm loops over all
of these candidate texels, filtering the contributions of those that are in fact inside the
ellipse. The bounding box is found by determining the minimum and maximum values
that the ellipse takes in the s and t directions. These extrema can be calculated by finding
the partial derivatives ∂e/∂s and ∂e/∂t , finding their solutions for s = 0 and t = 0, and
adding the offset to the ellipse center. For brevity, we will not include the derivation for
these expressions here.

 

d.

    T sum(0.);
    float sumWts = 0.f;
    for (int it = t0; it <= t1; ++it) {
        float tt = it - t;
        for (int is = s0; is <= s1; ++is) {
            float ss = is - s;
            // Compute squared radius and filter texel if inside ellipse
            float r2 = A*ss*ss + B*ss*tt + C*tt*tt;
            if (r2 < 1.) {
                float weight = weightLut[min(Float2Int(r2 * WEIGHT_LUT_SIZE),
                                             WEIGHT_LUT_SIZE-1)];
                sum += Texel(level, is, it) * weight;
                sumWts += weight;
            }
        }
    }
    return sum / sumWts;

作用：

（ 橢圓裡面的 包圍盒 確定下來了，那麼就可以開始遍歷這個包圍盒 所有的 texel，因為橢圓預設中心點在原點，所有會涉及到平移的操作，如何判斷哪一些texel是在橢圓裡面？ 因為 橢圓方程式 中的 F，表示的意思就是， 1個texel離橢圓中心的距離 與 橢圓中心沿著這個texel的方向到 橢圓邊的距離 的比值平方，其實就是 簡單理解就類似，圓內一個點 與 半徑的比值的感覺，如果這個比值是  小於 1 的話，就表示 這個點在 橢圓內， 如果這個texel是在橢圓內的話，就可以計算出這個texel相對於 橢圓重點的 權重值weight）

Now that the bounding box is known, the EWA algorithm loops over the texels, transforming
each one to the coordinate system where the texture lookup point (s , t) is at
the origin with a translation(平移). It then evaluates the ellipse equation to see if the texel is
inside the ellipse (Figure 10.16).The value of the implicit ellipse equation e(s , t) gives
the squared ratio of the distance from the texel to the ellipse center to the distance from

the ellipse edge to the center along the line through the texel. If it is inside, the weight
of the texel is computed with a Gaussian centered at the middle of the ellipse. The final
filtered value returned is a weighted sum over texels (s‘, t’) inside the ellipse, where f is
the Gaussian filter function:

Figure 10.16: Finding the r2 Ellipse Value for the EWA Filter Table Lookup.

 

A nice feature of the implicit equation e(s , t) is that its value at a particular texel is the
squared ratio of the distance from the center of the ellipse to the texel to the distance
from the center of the ellipse to the ellipse boundary along the line through that texel
(Figure 10.16). This value is used to index into a precomputed lookup table of Gaussian
filter function values.

 

e.

    #define WEIGHT_LUT_SIZE 128
    static float *weightLut;

    if (!weightLut) {
        weightLut = AllocAligned<float>(WEIGHT_LUT_SIZE);
        for (int i = 0; i < WEIGHT_LUT_SIZE; ++i) {
            float alpha = 2;
            float r2 = float(i) / float(WEIGHT_LUT_SIZE - 1);
            weightLut[i] = expf(-alpha * r2) - expf(-alpha);
        }
    }

作用：

(weightLut 在建構函式中進行初始化，高斯過濾)

The lookup table is initialized the first time a MIPMap is constructed. Because it will be
indexed with squared distances from the filter center r^2, each entry stores a value e^(−αr),
rather than e^(−αr^2)
.So that the filter function goes to zero at the end of its extent rather than having an abrupt(生硬)
step, expf(-alpha) is subtracted from the filter values here.