PtexExtractor/ptex/PtexMitchellFilter.cpp

PtexExtractor/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);
}