- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files
- Directories
- Examples

# PtexImporter/ptex/PtexMitchellFilter.cpp

    /*
    PTEX SOFTWARE
    Copyright 2009 Disney Enterprises, Inc.  All rights reserved

    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions are
    met:

      * Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

      * Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in
        the documentation and/or other materials provided with the
        distribution.

      * The names "Disney", "Walt Disney Pictures", "Walt Disney Animation
        Studios" or the names of its contributors may NOT be used to
        endorse or promote products derived from this software without
        specific prior written permission from Walt Disney Pictures.

    Disclaimer: THIS SOFTWARE IS PROVIDED BY WALT DISNEY PICTURES AND
    CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
    BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS
    FOR A PARTICULAR PURPOSE, NONINFRINGEMENT AND TITLE ARE DISCLAIMED.
    IN NO EVENT SHALL WALT DISNEY PICTURES, THE COPYRIGHT HOLDER OR
    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND BASED ON ANY
    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
    */

    #include "PtexPlatform.h"
    #include "PtexMitchellFilter.h"
    #include "PtexUtils.h"

    void PtexMitchellFilter::setSharpness(float sharpness)
    {
        _sharpness = sharpness;
        // compute Mitchell filter coefficients:
        // abs(x) < 1:
        //     1/6 * ((12 - 9*B - 6*C)*x^3 +
        //            (-18 + 12*B + 6*C)*x^2 +
        //        (6 - 2*B))
        //     == c[0]*x^3 + c[1]*x^2 + c[2]
        // abs(x) < 2:
        //       ((-B - 6*C)*x3 +
        //        (6*B + 30*C)*x2 +
        //        (-12*B - 48*C)*x +
        //        (8*B + 24*C)) :
        //     == c[3]*x^3 + c[4]*x^2 + c[5]*x + c[6]
        // else: 0

        double B = 1 - sharpness; // choose C = (1-B)/2
        _filter[0] = 1.5 - B;
        _filter[1] = 1.5 * B - 2.5;
        _filter[2] = 1 - (1./3) * B;
        _filter[3] = (1./3) * B - 0.5;
        _filter[4] = 2.5 - 1.5 * B;
        _filter[5] = 2 * B - 4;
        _filter[6] = 2 - (2./3) * B;
    }

    void PtexMitchellFilter::eval(float* result, int firstchan, int nchannels,
                      int faceid, float u, float v,
                      float uw1, float vw1, float uw2, float vw2,
                      float width, float blur)
    {
    #if 0
        // for debugging only!
        {
        // point sample highest res
        const FaceInfo& f = tx->getFaceInfo(faceid);
        int resu = f.res.u(), resv = f.res.v();
        int ui = PtexUtils::clamp(int(u * resu), 0, resu-1);
        int vi = PtexUtils::clamp(int(v * resv), 0, resv-1);
        tx->getPixel(faceid, ui, vi, result, firstchan, nchannels);
        return;
        }
    #endif

        /*
          Filter Scenarios:

          Regular - Local res available for all faces within kernel (and no e.p.)
          eval at local res only (local res = desired res clamped against local max)

          Blended - Desired res not available for all faces, or near e.p.
          blend between different resolution surfaces:
          a) local res: filter res clamped to res of local face
          b) uface res: local res clamped to res of u neighbor face
          c) vface res: local res clamped to res of v neighbor face
          d) cface res: local res clamped to res of u, v and corner faces
          Notes:
          - for all evals, adjacent faces are included if they have sufficent
          resolution; otherwise missing pixel values are taken from local face
          - reuse evals that have the same res

          u,v blend factors are 0..1 (smoothstepped) over blend region
          u blend region = 2 pixels given [vface ures at edge .. local ures mid-face]
          v blend region = 2 pixels given [uface vres at edge .. local vres mid-face]
          - smoothstep between edge/mid-face resolutions (mid-face = u or v == 0.5)
        */

        // find filter width as bounding box of vectors w1 and w2
        float uw = fabs(uw1) + fabs(uw2), vw = fabs(vw1) + fabs(vw2);

        if (!_ctx.prepare(result, firstchan, nchannels, _tx, faceid, u, v, uw, vw))
        return;

        double weight = OneValueInv(_ctx.dt);

        // get face
        const FaceInfo& f = _tx->getFaceInfo(faceid);
        _isConstant = f.isConstant();

        // clamp filter width to no larger than 0.25 (todo - handle larger filter widths)
        _ctx.uw = _ctx.uw * width + blur;
        _ctx.vw = _ctx.vw * width + blur;
        _ctx.uw = PtexUtils::min(_ctx.uw, 0.25f);
        _ctx.vw = PtexUtils::min(_ctx.vw, 0.25f);

        // clamp filter width to no smaller than a pixel
        _ctx.uw = PtexUtils::max(_ctx.uw, 1.0f/(f.res.u()));
        _ctx.vw = PtexUtils::max(_ctx.vw, 1.0f/(f.res.v()));

        // compute desired texture res based on filter width
        int ureslog2 = int(ceil(log2(1.0/_ctx.uw))),
        vreslog2 = int(ceil(log2(1.0/_ctx.vw)));
        _face.set(faceid, Res(ureslog2, vreslog2));

        // find neighboring faces and determine if neighborhood is constant
    #if 1
        getNeighborhood(f);
    #endif

        // if neighborhood is constant, just return constant value of face
        if (_isConstant) {
        PtexFaceData* data = _tx->getData(faceid, 0);
        if (data) {
            char* d = (char*) data->getData() + _ctx.firstchan*DataSize(_ctx.dt);
            Ptex::ConvertToFloat(_ctx.result, d, _ctx.dt, _ctx.nchannels);
            data->release();
        }
        return;
        }

        // if we don't have any neighbors, just eval face by itself
        if (!_uface && !_vface) {
        evalFaces(_face.res, weight, _ctx.uw, _ctx.vw);
            return;
        }

        // compute weights for the 4 possible surfaces
        double mweight = weight * (1 - _ublend) * (1 - _vblend); // main face
        double uweight = weight * _ublend * (1 - _vblend);       // u blend
        double vweight = weight * (1 - _ublend) * _vblend;       // v blend
        double cweight = weight * _ublend * _vblend;         // corner blend

        // eval surfaces using given weights
        // but first combine weights where surfaces are missing or no blending is required
        if (cweight) {
        if (_cface) {
            // combine w/ main if no blend needed
            if (!_cface.blend) mweight += cweight;
            // combine with an edge if res matches
            else if (_cface.res == _uface.res) uweight += cweight;
            else if (_cface.res == _vface.res) vweight += cweight;
            else evalFaces(_cface.res, cweight);
        }
        // else if (_cfaces.size() > 1) {
        //     // todo - near an e.p., blend using linear rbf
        // }
        else mweight += cweight;
        }
        if (uweight) {
        if (!_uface.blend) mweight += uweight;
        else if (_vface && (_uface.res == _vface.res)) vweight += uweight;
        else evalFaces(_uface.res, uweight);
        }
        if (vweight) {
        if (!_vface.blend) mweight += vweight;
        else evalFaces(_vface.res, vweight);
        }
        if (mweight) {
        evalFaces(_face.res, mweight, _ctx.uw, _ctx.vw);
        }
    }

    void PtexMitchellFilter::getNeighborhood(const FaceInfo& f)
    {
        _uface.clear(); _vface.clear(), _cface.clear(), _cfaces.clear();

        EdgeId ueid, veid;
        double udist, vdist;
        if (_ctx.u < .5) { ueid = e_left;   udist = _ctx.u; }
        else             { ueid = e_right;  udist = 1 - _ctx.u; }
        if (_ctx.v < .5) { veid = e_bottom; vdist = _ctx.v; }
        else             { veid = e_top;    vdist = 1 - _ctx.v; }

        // blend zone starts at 1.5 texels in (at the current res)
        // to keep the kernel from accessing samples from the low-res face
        static const double blendstart = 1.5; // dist from edge (in texels) to start blending
        static const double blendend = 2.5;   // dist from edge (in texels) to stop blending

        // get u and v neighbors and compute blend weights
        double ublendstart=0, ublendend=0, vblendstart=0, vblendend=0;
        int ufid = f.adjfaces[ueid], vfid = f.adjfaces[veid];
        const FaceInfo* uf = 0, * vf = 0;
        if (ufid != -1) {
        // compute blend distances in u dir
        double texel = 1.0/_face.res.u();
        ublendstart = PtexUtils::min(blendstart*texel, .375);
        ublendend = PtexUtils::min(blendend*texel, 0.5);
        // get u neighbor face
        uf = &_ctx.tx->getFaceInfo(ufid);
        _uface.set(ufid, uf->res, f.adjedge(ueid) - ueid + 2);
        // clamp res against desired res and check for blending
        _uface.clampres(_face.res);
        }
        if (vfid != -1)  {
        // compute blend distances in v dir
        double texel = 1.0/_face.res.v();
        vblendstart = PtexUtils::min(blendstart*texel, .375);
        vblendend = PtexUtils::min(blendend*texel, 0.5);
        // get v neighbor face
        vf = &_ctx.tx->getFaceInfo(vfid);
        _vface.set(vfid, vf->res, f.adjedge(veid) - veid + 2);
        // clamp res against desired res and check for blending
        _vface.clampres(_face.res);
        }

        // smoothstep blendwidths towards corner
        if (_uface && _vface) {
        // smoothstep ublendwidth towards adj ures
        if (_vface.res.ulog2 != _face.res.ulog2) {
            double texel = 1.0/_vface.res.u();
            double adjstart = PtexUtils::min(blendstart*texel, .375);
            double adjend = PtexUtils::min(blendend*texel, .5);
            double wblend = PtexUtils::smoothstep(vdist, vblendstart, vblendend);
            ublendstart = ublendstart * wblend + adjstart * (1-wblend);
            ublendend = ublendend * wblend + adjend * (1-wblend);
        }
        // smoothstep vblendwidth towards adj vres
        if (_uface.res.vlog2 != _face.res.vlog2) {
            double texel = 1.0/_uface.res.v();
            double adjstart = PtexUtils::min(blendstart*texel, .375);
            double adjend = PtexUtils::min(blendend*texel, .5);
            double wblend = PtexUtils::smoothstep(udist, ublendstart, ublendend);
            vblendstart = vblendstart * wblend + adjstart * (1-wblend);
            vblendend = vblendend * wblend + adjend * (1-wblend);
        }
        }

        // compute blend weights based on distances to edges
        bool nearu = _uface && (udist < ublendend);
        bool nearv = _vface && (vdist < vblendend);

        if (!nearu) {
        _ublend = 0;
        _uface.clear();
        }
        else {
        // in blend zone w/ u
        _ublend = 1 - PtexUtils::qsmoothstep(udist, ublendstart, ublendend);
        if (_isConstant && !uf->isConstant()) _isConstant = 0;
        }

        if (!nearv) {
        _vblend = 0;
        _vface.clear();
        }
        else {
        // in blend zone w/ v
        _vblend = 1 - PtexUtils::qsmoothstep(vdist, vblendstart, vblendend);
        if (_isConstant && !vf->isConstant()) _isConstant = 0;
        }

        // gather corner faces if needed
        _interior = false;
        if (nearu && nearv) {
        // gather faces around corner starting at uface
        _cfaces.reserve(8); // (could be any number, but typically just 1 or 2)

        int cfid = ufid;           // current face id (start at uface)
        const FaceInfo* cf = uf;       // current face info
        EdgeId ceid = f.adjedge(ueid);     // current edge id
        int rotate = _uface.rotate;    // cumulative rotation
        int dir = (ueid+1)%4==veid? 3 : 1; // ccw or cw
        int count = 0;             // runaway loop count

        while (count++ < 10) {
            // advance to next face
            int eid = EdgeId((ceid + dir) % 4);
            cfid = cf->adjfaces[eid];
            if (cfid == _vface.id || cfid == -1)
            // reached "vface" or boundary, stop
            break;
            ceid = cf->adjedge(eid);
            cf = &_ctx.tx->getFaceInfo(cfid);
            rotate += ceid - eid + 2;

            // record face (note: first face (uface) is skipped)
            _cfaces.push_back(Face());
            Face& face = _cfaces.back();
            face.set(cfid, cf->res, rotate);

            if (_isConstant && !cf->isConstant()) _isConstant = 0;
        }
        // if we reached the vface, corner is an interior point
        if (cfid == _vface.id) {
            _interior = true;
            // see if we're regular - i.e. we have a single corner face
            if (_cfaces.size() == 1) {
            _cface = _cfaces.front();

            // clamp res against u and v neighbors and check for blending
            _cface.clampres(_uface.res);
            _cface.clampres(_vface.res);

            // if either u or v needs blending, corner needs blending too
            if (_uface.blend || _vface.blend)
                _cface.blend = true;
            }
        }
        // otherwise corner is an e.p. on a mesh boundary
        else {
            // don't blend w/ corner faces for this case
            _cfaces.clear();
        }
        }

        // if all faces are constant, see if all have the same value
        if (_isConstant) {
        int pixelsize = DataSize(_ctx.dt) * _ctx.ntxchannels;
        PtexFaceData* data = _ctx.tx->getData(_face.id, 0);
        if (data) {
            void* constval = data->getData();
            if (_uface) {
            PtexFaceData* udata = _ctx.tx->getData(_uface.id, 0);
            if (udata) {
                if (0 != memcmp(constval, udata->getData(), pixelsize))
                _isConstant = 0;
                udata->release();
            }
            }
            if (_isConstant && _vface) {
            PtexFaceData* vdata = _ctx.tx->getData(_vface.id, 0);
            if (vdata) {
                if (0 != memcmp(constval, vdata->getData(), pixelsize))
                _isConstant = 0;
                vdata->release();
            }
            }
            if (_isConstant) {
            for (size_t i = 0, size = _cfaces.size(); i < size; i++) {
                PtexFaceData* cdata = _ctx.tx->getData(_cfaces[i].id, 0);
                if (cdata) {
                if (0 != memcmp(constval, cdata->getData(), pixelsize)) {
                    _isConstant = 0;
                    break;
                }
                cdata->release();
                }
            }
            }
            data->release();
        }
        }
    }

    void PtexMitchellFilter::evalFaces(Res res, double weight, float uw, float vw)
    {
        // initialize kernel
        // convert u,v and filter width to (fractional) pixels
        int ures = res.u(), vres = res.v();

        if (ures < 4 || vres < 4) {
        // can't use 4x4 mitchell, just do bilinear interp.
        // this should only happen for very small faces so smooth filtering shouldn't
        // be needed (fingers crossed)

        // todo - build 2x2 bilinear kernel
        PtexFilterKernel k;
        k.set(0, 0, 0, 1, 1, &weight, 0);
        k.apply(_face.id, 0, _ctx);
        return;
        }

        double u = _ctx.u * ures - 0.5, v = _ctx.v * vres - 0.5;
        uw *= ures; vw *= vres;

        // find integer pixel extent: [u,v] +/- [2*uw,2*vw]
        // (mitchell is 4 units wide for a 1 unit filter period)

        int u1 = int(ceil(u - 2*uw)), u2 = int(ceil(u + 2*uw));
        int v1 = int(ceil(v - 2*vw)), v2 = int(ceil(v + 2*vw));

        int kuw = u2-u1, kvw = v2-v1;
        if (kuw > 8 || kvw > 8) {
        // shouldn't happen - but just in case...
        assert(kuw <= 8 && kvw <= 8);
        return;
        }

        // compute kernel weights along u and v directions
        double* ukernel = (double*)alloca(kuw * sizeof(double));
        double* vkernel = (double*)alloca(kvw * sizeof(double));
        computeWeights(ukernel, (u1-u)/uw, 1.0/uw, kuw);
        computeWeights(vkernel, (v1-v)/vw, 1.0/vw, kvw);
        double scale = weight;

        // skip zero entries (will save a lot of work later)
        while (!ukernel[0])     { ukernel++; u1++; kuw--; }
        while (!ukernel[kuw-1]) { kuw--; }
        while (!vkernel[0])     { vkernel++; v1++; kvw--; }
        while (!vkernel[kvw-1]) { kvw--; }

        double sumu = 0; for (int i = 0; i < kuw; i++) sumu += ukernel[i];
        double sumv = 0; for (int i = 0; i < kvw; i++) sumv += vkernel[i];
        scale /= sumu * sumv;

        // compute tensor product to form rectangular kernel
        double* kbuffer = (double*) alloca(kuw*kvw*sizeof(double));
        double* kp = kbuffer;
        for (int i = 0; i < kvw; i++, kp += kuw) {
        double vk = vkernel[i] * scale;
        for (int j = 0; j < kuw; j++) kp[j] = ukernel[j]*vk;
        }
        PtexFilterKernel k; k.set(res, u1, v1, kuw, kvw, kbuffer, kuw);

        // split kernel across edges into face, u,v, and corner parts
        PtexFilterKernel ku, kv, kc;
        k.split(ku, kv, kc);

        if (ku || kv) {
        // merge kernel parts back in for missing/insufficient-res faces
        if (kc) {
            if (!_cface && _interior) {
            // valence-3 interior case
            // clear corner and renormalize kernel
            double amt = 1.0/(1 - kc.totalWeight()/weight);
            kc.clear();
            for (double *kp = kbuffer, *end = kbuffer + kuw*kvw; kp != end; kp++)
                *kp *= amt;
            }
            else if (!_cface || !(_cface.res >= res)) {
            // merge corner into u and/or v faces
            if (kv && _uface) {
                if (_vface) {
                // merge corner 50% into ku and kv faces
                ku.merge(kc, kv.eidval(), 0.5);
                kv.merge(kc, ku.eidval(), 0.5);
                }
                else {
                // merge corner into ku across v edge
                ku.merge(kc, kv.eidval());
                }
            } else {
                // merge corner into kv across u edge
                kv.merge(kc, ku.eidval());
            }
            }
        }
        // merge boundary edges into main kernel
        if (ku && (!_uface || !(_uface.res >= res)))
            k.merge(ku, ku.eidval());
        if (kv && (!_vface || !(_vface.res >= res)))
            k.merge(kv, kv.eidval());

        if (ku) ku.apply(_uface.id, _uface.rotate, _ctx);
        if (kv) kv.apply(_vface.id, _vface.rotate, _ctx);
        if (kc) kc.apply(_cface.id, _cface.rotate, _ctx);
        }

        // eval faces
        k.apply(_face.id, 0, _ctx);
    }