/disk1/scratch/aselle/projects/seexpr/SeExpr/src/SeExpr/SeCurve.cpp

Go to the documentation of this file.

/*
* (c) Disney Enterprises, Inc.  All rights reserved.
*
* This file is licensed under the terms of the Microsoft Public License (MS-PL)
* as defined at: http://opensource.org/licenses/MS-PL.
*
* A complete copy of this license is included in this distribution as the file
* LICENSE.
*/
#include "SeExpression.h"
#include "SeExprBuiltins.h"
#include <cfloat>
#include <cassert>
#include <algorithm>

#include "SeCurve.h"

namespace SeExpr{

template<> double SeCurve<double>::comp(const double& val,const int)
{return val;}

template<> double SeCurve<SeVec3d>::comp(const SeVec3d& val,const int i)
{return val[i];}
    
template<class T> bool SeCurve<T>::
cvLessThan(const CV &cv1, const CV &cv2)
{
    return cv1._pos < cv2._pos;
}

template<class T> SeCurve<T>::
SeCurve()
    :cacheCV(0),prepared(false)
{
    _cvData.push_back(CV(-FLT_MAX,T(),kNone));
    _cvData.push_back(CV(FLT_MAX,T(),kNone));
}

template<class T> void SeCurve<T>::
addPoint(double position,const T& val,InterpType type)
{
    prepared=false;
    _cvData.push_back(CV(position,val,type));
}

template<class T> void SeCurve<T>::
preparePoints()
{
    prepared=true;
    cacheCV=0;
    // sort
    std::sort(_cvData.begin(), _cvData.end(), cvLessThan);
    
    // Setup boundary conditions on sentinel values
    CV& end=*(_cvData.end()-1);
    CV& begin=*(_cvData.begin());
    int realCVs=_cvData.size()-2;
    assert(realCVs>=0);
    if(realCVs>0){
        begin._val=_cvData[1]._val;
        begin._deriv=T();
        begin._interp=kNone;
        int lastIndex=_cvData.size()-1;
        end._val=_cvData[lastIndex-1]._val;
        end._deriv=T();
        end._interp=kNone;
    }else{
        begin._pos=end._pos=0;
        begin._val=end._val=T();
        begin._interp=kNone;
        begin._deriv=end._deriv=T();
    }
            
    // Initialize "Catmull-Rom" derivatives (centered differences)
    for(unsigned int i=1;i<_cvData.size()-1;i++){
        _cvData[i]._deriv=(_cvData[i+1]._val-_cvData[i-1]._val)
            /(_cvData[i+1]._pos-_cvData[i-1]._pos);
    }

    // Fix extrema by going through all intervals
    for(unsigned int i=0;i<_cvData.size()-1;i++){
        if(_cvData[i]._interp==kMonotoneSpline){
            double h=_cvData[i+1]._pos-_cvData[i]._pos;
            if(h==0)
                _cvData[i]._deriv=_cvData[i+1]._deriv=T();
            else{
                T delta=(_cvData[i+1]._val-_cvData[i]._val)/h;
                clampCurveSegment(delta,_cvData[i]._deriv,
                    _cvData[i+1]._deriv);
            }
        }
    }
}

// TODO: this function and the next could be merged with template magic
//       but it might be simpler to just have two copies!
template<class T> T SeCurve<T>::
getValue(const double param) const
{
    assert(prepared);
    // find the cv data point index just greater than the desired param
    const int numPoints = _cvData.size();
    const CV *cvDataBegin = &_cvData[0];
    int index = std::upper_bound(cvDataBegin, cvDataBegin + numPoints, 
        CV(param, T(),kLinear),cvLessThan) - cvDataBegin;
    index=std::max(1,std::min(index,numPoints-1));

    const float t0 = _cvData[index - 1]._pos;
    const T k0 = _cvData[index - 1]._val;
    const InterpType interp = _cvData[index - 1]._interp;
    const float t1 = _cvData[index]._pos;
    const T k1 = _cvData[index]._val;
    switch (interp) {
        case kNone:
            return k0;
            break;
        case kLinear:
            {
                double u = (param-t0)/(t1-t0);
                return k0 + u*(k1-k0);
            }
            break;
        case kSmooth:
            {
                double u = (param-t0)/(t1-t0);
                return k0*(u-1)*(u-1)*(2*u+1) + k1*u*u*(3 - 2*u);
            }
            break;
        case kSpline:
        case kMonotoneSpline:
            {
                double x=param-_cvData[index-1]._pos; // xstart
                double h=_cvData[index]._pos-_cvData[index-1]._pos; // xend-xstart
                T y=_cvData[index-1]._val; // f(xstart)
                T delta=_cvData[index]._val-_cvData[index-1]._val; // f(xend)-f(xstart)
                T d1=_cvData[index-1]._deriv; // f'(xstart)
                T d2=_cvData[index]._deriv; // f'(xend)
                return (x*(delta*(3*h - 2*x)*x + h*(-h + x)*(-(d1*h) + (d1 + d2)*x)))
                    /(h*h*h) + y;
            }
            break;
        default:
            assert(false);
            return T();
            break;
    }    
}

// TODO: this function and the previous could be merged with template magic
//       but it might be simpler to just have two copies!
template<class T> double SeCurve<T>::
getChannelValue(const double param,int channel) const{
    assert(prepared);
    // find the cv data point index just greater than the desired param
    const int numPoints = _cvData.size();
    const CV *cvDataBegin = &_cvData[0];
    int index = std::upper_bound(cvDataBegin, cvDataBegin + numPoints, 
        CV(param, T(),kLinear),cvLessThan) - cvDataBegin;
    index=std::max(1,std::min(index,numPoints-1));
    
    const float t0 = _cvData[index - 1]._pos;
    const double k0 = comp(_cvData[index - 1]._val,channel);
    const InterpType interp = _cvData[index - 1]._interp;
    const float t1 = _cvData[index]._pos;
    const double k1 = comp(_cvData[index]._val,channel);
    switch (interp) {
        case kNone:
            return k0;
            break;
        case kLinear:
            {
                double u = (param-t0)/(t1-t0);
                return k0 + u*(k1-k0);
            }
            break;
        case kSmooth:
            // standard cubic interpolation
            {
                double u = (param-t0)/(t1-t0);
                return k0*(u-1)*(u-1)*(2*u+1) + k1*u*u*(3 - 2*u);
            }
            break;
        case kSpline:
        case kMonotoneSpline:
            {
                double x=param-_cvData[index-1]._pos; // xstart
                double h=_cvData[index]._pos-_cvData[index-1]._pos; // xend-xstart
                double y=comp(_cvData[index-1]._val,channel); // f(xtart)
                double delta=comp(_cvData[index]._val,channel)
                    -comp(_cvData[index-1]._val,channel); // f(xend)-f(xtart)
                double d1=comp(_cvData[index-1]._deriv,channel); // f'(xtart)
                double d2=comp(_cvData[index]._deriv,channel); // f'(xend)
                
                return (x*(delta*(3*h - 2*x)*x + h*(-h + x)*(-(d1*h) + (d1 + d2)*x)))
                    /(h*h*h) + y;
            }
            break;
        default:
            assert(false);
            return 0;
            break;
    }
}

template<class T> typename SeCurve<T>::CV SeCurve<T>::
getLowerBoundCV(const double param) const
{
    assert(prepared);
    const CV *cvDataBegin = &_cvData[0];
    int numPoints=_cvData.size();
    int index = std::upper_bound(cvDataBegin, cvDataBegin + numPoints, 
    CV(param, T(),kLinear),cvLessThan) - cvDataBegin;
    index=std::max(1,std::min(index,numPoints-1));
    if(index-1 > 0) return _cvData[index-1];
    return _cvData[index];
}

template<class T> bool SeCurve<T>::
interpTypeValid(InterpType interp)
{
      return interp==kNone || interp==kLinear || interp==kSmooth 
          || interp==kSpline || interp==kMonotoneSpline;
}

template<>
inline void SeCurve<double>::clampCurveSegment(const double& delta,double& d1,double& d2)
{
    if(delta==0) d1=d2=0;
    else{
        d1=SeExpr::clamp(d1/delta,0,3)*delta;
        d2=SeExpr::clamp(d2/delta,0,3)*delta;
    }
}

template<>
void SeCurve<SeVec3d>::clampCurveSegment(const SeVec3d& delta,SeVec3d& d1,SeVec3d& d2)
{
    for(int i=0;i<3;i++){
        if(delta[i]==0) d1[i]=d2[i]=0;
        else{
            d1[i]=SeExpr::clamp(d1[i]/delta[i],0,3)*delta[i];
            d2[i]=SeExpr::clamp(d2[i]/delta[i],0,3)*delta[i];
        }
    }
}

template class SeCurve<SeVec3d>;
template class SeCurve<double>;

}

---